1

2# -*- coding: utf-8 -*-

3

4 u'''Ellipsoidal and spherical earth models.

5

6 Classes L{a_f2Tuple}, L{Ellipsoid} and L{Ellipsoid2}, an L{Ellipsoids} registry and

7 2 dozen functions to convert I{equatorial} radius, I{polar} radius, I{eccentricities},

8 I{flattenings} and I{inverse flattening}.

9

10 See module L{datums} for L{Datum} and L{Transform} information and other details.

11

12 Following is the list of predefined L{Ellipsoid}s, all instantiated lazily.

13

14@var Ellipsoids.Airy1830: Ellipsoid(name='Airy1830', a=6377563.396, b=6356256.90923729, f_=299.3249646, f=0.00334085, f2=0.00335205, n=0.00167322, e=0.08167337, e2=0.00667054, e21=0.99332946, e22=0.00671533, e32=0.00334643, A=6366914.60892522, L=10001126.0807165, R1=6370461.23374576, R2=6370459.65470808, R3=6370453.30994572, Rbiaxial=6366919.065224, Rtriaxial=6372243.45317691)

15@var Ellipsoids.AiryModified: Ellipsoid(name='AiryModified', a=6377340.189, b=6356034.44793853, f_=299.3249646, f=0.00334085, f2=0.00335205, n=0.00167322, e=0.08167337, e2=0.00667054, e21=0.99332946, e22=0.00671533, e32=0.00334643, A=6366691.77461988, L=10000776.05340819, R1=6370238.27531284, R2=6370236.69633043, R3=6370230.35179013, Rbiaxial=6366696.2307627, Rtriaxial=6372020.43236847)

16@var Ellipsoids.ATS1977: Ellipsoid(name='ATS1977', a=6378135, b=6356750.30492159, f_=298.257, f=0.00335281, f2=0.00336409, n=0.00167922, e=0.08181922, e2=0.00669438, e21=0.99330562, e22=0.0067395, e32=0.00335843, A=6367447.14116695, L=10001962.58040571, R1=6371006.7683072, R2=6371005.17780873, R3=6370998.78689182, Rbiaxial=6367451.62986519, Rtriaxial=6372795.55363648)

17@var Ellipsoids.Australia1966: Ellipsoid(name='Australia1966', a=6378160, b=6356774.71919531, f_=298.25, f=0.00335289, f2=0.00336417, n=0.00167926, e=0.08182018, e2=0.00669454, e21=0.99330546, e22=0.00673966, e32=0.00335851, A=6367471.84853228, L=10002001.39064442, R1=6371031.5730651, R2=6371029.9824858, R3=6371023.59124344, Rbiaxial=6367476.337459, Rtriaxial=6372820.40754721)

18@var Ellipsoids.Bessel1841: Ellipsoid(name='Bessel1841', a=6377397.155, b=6356078.962818, f_=299.1528128, f=0.00334277, f2=0.00335398, n=0.00167418, e=0.08169683, e2=0.00667437, e21=0.99332563, e22=0.00671922, e32=0.00334836, A=6366742.52023395, L=10000855.76443237, R1=6370291.09093933, R2=6370289.51012659, R3=6370283.15821523, Rbiaxial=6366746.98155108, Rtriaxial=6372074.29334012)

19@var Ellipsoids.BesselModified: Ellipsoid(name='BesselModified', a=6377492.018, b=6356173.5087127, f_=299.1528128, f=0.00334277, f2=0.00335398, n=0.00167418, e=0.08169683, e2=0.00667437, e21=0.99332563, e22=0.00671922, e32=0.00334836, A=6366837.22474766, L=10001004.52593463, R1=6370385.84823756, R2=6370384.26740131, R3=6370377.91539546, Rbiaxial=6366841.68613115, Rtriaxial=6372169.07716325)

20@var Ellipsoids.CGCS2000: Ellipsoid(name='CGCS2000', a=6378137, b=6356752.31414036, f_=298.2572221, f=0.00335281, f2=0.00336409, n=0.00167922, e=0.08181919, e2=0.00669438, e21=0.99330562, e22=0.0067395, e32=0.00335843, A=6367449.14577105, L=10001965.72923046, R1=6371008.77138012, R2=6371007.18088352, R3=6371000.78997414, Rbiaxial=6367453.63446401, Rtriaxial=6372797.55593326)

21@var Ellipsoids.Clarke1866: Ellipsoid(name='Clarke1866', a=6378206.4, b=6356583.8, f_=294.97869821, f=0.00339008, f2=0.00340161, n=0.00169792, e=0.08227185, e2=0.00676866, e21=0.99323134, e22=0.00681478, e32=0.00339582, A=6367399.68916978, L=10001888.04298286, R1=6370998.86666667, R2=6370997.240633, R3=6370990.70659881, Rbiaxial=6367404.2783313, Rtriaxial=6372807.62791066)

22@var Ellipsoids.Clarke1880: Ellipsoid(name='Clarke1880', a=6378249.145, b=6356514.86954978, f_=293.465, f=0.00340756, f2=0.00341921, n=0.00170669, e=0.0824834, e2=0.00680351, e21=0.99319649, e22=0.00685012, e32=0.00341337, A=6367386.64398051, L=10001867.55164747, R1=6371004.38651659, R2=6371002.74366963, R3=6370996.1419165, Rbiaxial=6367391.2806777, Rtriaxial=6372822.52526083)

23@var Ellipsoids.Clarke1880IGN: Ellipsoid(name='Clarke1880IGN', a=6378249.2, b=6356515, f_=293.46602129, f=0.00340755, f2=0.0034192, n=0.00170668, e=0.08248326, e2=0.00680349, e21=0.99319651, e22=0.00685009, e32=0.00341336, A=6367386.73667336, L=10001867.69724907, R1=6371004.46666667, R2=6371002.82383112, R3=6370996.22212395, Rbiaxial=6367391.37333829, Rtriaxial=6372822.59907505)

24@var Ellipsoids.Clarke1880Mod: Ellipsoid(name='Clarke1880Mod', a=6378249.145, b=6356514.96639549, f_=293.46630766, f=0.00340755, f2=0.0034192, n=0.00170668, e=0.08248322, e2=0.00680348, e21=0.99319652, e22=0.00685009, e32=0.00341335, A=6367386.69236201, L=10001867.62764495, R1=6371004.4187985, R2=6371002.77596616, R3=6370996.17427195, Rbiaxial=6367391.32901784, Rtriaxial=6372822.5494103)

25@var Ellipsoids.CPM1799: Ellipsoid(name='CPM1799', a=6375738.7, b=6356671.92557493, f_=334.39, f=0.00299052, f2=0.00299949, n=0.0014975, e=0.07727934, e2=0.0059721, e21=0.9940279, e22=0.00600798, e32=0.00299499, A=6366208.88184784, L=10000017.52721564, R1=6369383.10852498, R2=6369381.8434158, R3=6369376.76247022, Rbiaxial=6366212.45090321, Rtriaxial=6370977.3559758)

26@var Ellipsoids.Delambre1810: Ellipsoid(name='Delambre1810', a=6376428, b=6355957.92616372, f_=311.5, f=0.00321027, f2=0.00322061, n=0.00160772, e=0.08006397, e2=0.00641024, e21=0.99358976, e22=0.0064516, e32=0.00321543, A=6366197.07684334, L=9999998.98395793, R1=6369604.64205457, R2=6369603.18419749, R3=6369597.32739068, Rbiaxial=6366201.19059818, Rtriaxial=6371316.64722284)

27@var Ellipsoids.Engelis1985: Ellipsoid(name='Engelis1985', a=6378136.05, b=6356751.32272154, f_=298.2566, f=0.00335282, f2=0.0033641, n=0.00167922, e=0.08181928, e2=0.00669439, e21=0.99330561, e22=0.00673951, e32=0.00335844, A=6367448.17507971, L=10001964.20447208, R1=6371007.80757385, R2=6371006.21707085, R3=6370999.82613573, Rbiaxial=6367452.66379074, Rtriaxial=6372796.59560563)

28@var Ellipsoids.Everest1969: Ellipsoid(name='Everest1969', a=6377295.664, b=6356094.667915, f_=300.8017, f=0.00332445, f2=0.00333554, n=0.00166499, e=0.08147298, e2=0.00663785, e21=0.99336215, e22=0.0066822, e32=0.00332998, A=6366699.57839501, L=10000788.3115495, R1=6370228.665305, R2=6370227.10178537, R3=6370220.81951618, Rbiaxial=6366703.99082487, Rtriaxial=6372002.02812501)

29@var Ellipsoids.Everest1975: Ellipsoid(name='Everest1975', a=6377299.151, b=6356098.14512013, f_=300.8017255, f=0.00332445, f2=0.00333554, n=0.00166499, e=0.08147298, e2=0.00663785, e21=0.99336215, e22=0.0066822, e32=0.00332997, A=6366703.06049924, L=10000793.78122603, R1=6370232.14904004, R2=6370230.58551983, R3=6370224.30324826, Rbiaxial=6366707.47293076, Rtriaxial=6372005.51267879)

30@var Ellipsoids.Fisher1968: Ellipsoid(name='Fisher1968', a=6378150, b=6356768.33724438, f_=298.3, f=0.00335233, f2=0.00336361, n=0.00167898, e=0.08181333, e2=0.00669342, e21=0.99330658, e22=0.00673853, e32=0.00335795, A=6367463.65604381, L=10001988.52191361, R1=6371022.77908146, R2=6371021.18903735, R3=6371014.79995035, Rbiaxial=6367468.14345752, Rtriaxial=6372811.30979281)

31@var Ellipsoids.GEM10C: Ellipsoid(name='GEM10C', a=6378137, b=6356752.31424783, f_=298.2572236, f=0.00335281, f2=0.00336409, n=0.00167922, e=0.08181919, e2=0.00669438, e21=0.99330562, e22=0.0067395, e32=0.00335843, A=6367449.14582474, L=10001965.7293148, R1=6371008.77141594, R2=6371007.18091936, R3=6371000.79001005, Rbiaxial=6367453.63451765, Rtriaxial=6372797.55596006)

32@var Ellipsoids.GPES: Ellipsoid(name='GPES', a=6378135, b=6378135, f_=0, f=0, f2=0, n=0, e=0, e2=0, e21=1, e22=0, e32=0, A=6378135, L=10018751.02980197, R1=6378135, R2=6378135, R3=6378135, Rbiaxial=6378135, Rtriaxial=6378135)

33@var Ellipsoids.GRS67: Ellipsoid(name='GRS67', a=6378160, b=6356774.51609071, f_=298.24716743, f=0.00335292, f2=0.0033642, n=0.00167928, e=0.08182057, e2=0.00669461, e21=0.99330539, e22=0.00673973, e32=0.00335854, A=6367471.74706533, L=10002001.2312605, R1=6371031.50536357, R2=6371029.91475409, R3=6371023.52339015, Rbiaxial=6367476.23607738, Rtriaxial=6372820.3568989)

34@var Ellipsoids.GRS80: Ellipsoid(name='GRS80', a=6378137, b=6356752.31414035, f_=298.2572221, f=0.00335281, f2=0.00336409, n=0.00167922, e=0.08181919, e2=0.00669438, e21=0.99330562, e22=0.0067395, e32=0.00335843, A=6367449.14577104, L=10001965.72923046, R1=6371008.77138012, R2=6371007.18088351, R3=6371000.78997414, Rbiaxial=6367453.634464, Rtriaxial=6372797.55593326)

35@var Ellipsoids.Helmert1906: Ellipsoid(name='Helmert1906', a=6378200, b=6356818.16962789, f_=298.3, f=0.00335233, f2=0.00336361, n=0.00167898, e=0.08181333, e2=0.00669342, e21=0.99330658, e22=0.00673853, e32=0.00335795, A=6367513.57227074, L=10002066.93013953, R1=6371072.7232093, R2=6371071.13315272, R3=6371064.74401563, Rbiaxial=6367518.05971963, Rtriaxial=6372861.26794141)

36@var Ellipsoids.IAU76: Ellipsoid(name='IAU76', a=6378140, b=6356755.28815753, f_=298.257, f=0.00335281, f2=0.00336409, n=0.00167922, e=0.08181922, e2=0.00669438, e21=0.99330562, e22=0.0067395, e32=0.00335843, A=6367452.13278844, L=10001970.42122641, R1=6371011.76271918, R2=6371010.17221946, R3=6371003.78129754, Rbiaxial=6367456.6214902, Rtriaxial=6372800.54945074)

37@var Ellipsoids.IERS1989: Ellipsoid(name='IERS1989', a=6378136, b=6356751.30156878, f_=298.257, f=0.00335281, f2=0.00336409, n=0.00167922, e=0.08181922, e2=0.00669438, e21=0.99330562, e22=0.0067395, e32=0.00335843, A=6367448.13949125, L=10001964.14856985, R1=6371007.76718959, R2=6371006.17669088, R3=6370999.78577297, Rbiaxial=6367452.62819019, Rtriaxial=6372796.55279934)

38@var Ellipsoids.IERS1992TOPEX: Ellipsoid(name='IERS1992TOPEX', a=6378136.3, b=6356751.61659215, f_=298.25722356, f=0.00335281, f2=0.00336409, n=0.00167922, e=0.08181919, e2=0.00669438, e21=0.99330562, e22=0.0067395, e32=0.00335843, A=6367448.44699641, L=10001964.63159782, R1=6371008.07219738, R2=6371006.48170097, R3=6371000.09079236, Rbiaxial=6367452.93568883, Rtriaxial=6372796.85654541)

39@var Ellipsoids.IERS2003: Ellipsoid(name='IERS2003', a=6378136.6, b=6356751.85797165, f_=298.25642, f=0.00335282, f2=0.0033641, n=0.00167922, e=0.0818193, e2=0.0066944, e21=0.9933056, e22=0.00673951, e32=0.00335844, A=6367448.71771058, L=10001965.05683465, R1=6371008.35265722, R2=6371006.76215217, R3=6371000.37120877, Rbiaxial=6367453.20642742, Rtriaxial=6372797.14192686)

40@var Ellipsoids.Intl1924: Ellipsoid(name='Intl1924', a=6378388, b=6356911.94612795, f_=297, f=0.003367, f2=0.00337838, n=0.00168634, e=0.08199189, e2=0.00672267, e21=0.99327733, e22=0.00676817, e32=0.00337267, A=6367654.50005758, L=10002288.29898944, R1=6371229.31537598, R2=6371227.71133444, R3=6371221.26587487, Rbiaxial=6367659.02704315, Rtriaxial=6373025.77129687)

41@var Ellipsoids.Intl1967: Ellipsoid(name='Intl1967', a=6378157.5, b=6356772.2, f_=298.24961539, f=0.0033529, f2=0.00336418, n=0.00167926, e=0.08182023, e2=0.00669455, e21=0.99330545, e22=0.00673967, e32=0.00335852, A=6367469.33894446, L=10001997.44859308, R1=6371029.06666667, R2=6371027.47608389, R3=6371021.08482752, Rbiaxial=6367473.827881, Rtriaxial=6372817.9027631)

42@var Ellipsoids.Krassovski1940: Ellipsoid(name='Krassovski1940', a=6378245, b=6356863.01877305, f_=298.3, f=0.00335233, f2=0.00336361, n=0.00167898, e=0.08181333, e2=0.00669342, e21=0.99330658, e22=0.00673853, e32=0.00335795, A=6367558.49687498, L=10002137.49754285, R1=6371117.67292435, R2=6371116.08285656, R3=6371109.69367439, Rbiaxial=6367562.98435553, Rtriaxial=6372906.23027515)

43@var Ellipsoids.Krassowsky1940: Ellipsoid(name='Krassowsky1940', a=6378245, b=6356863.01877305, f_=298.3, f=0.00335233, f2=0.00336361, n=0.00167898, e=0.08181333, e2=0.00669342, e21=0.99330658, e22=0.00673853, e32=0.00335795, A=6367558.49687498, L=10002137.49754285, R1=6371117.67292435, R2=6371116.08285656, R3=6371109.69367439, Rbiaxial=6367562.98435553, Rtriaxial=6372906.23027515)

44@var Ellipsoids.Maupertuis1738: Ellipsoid(name='Maupertuis1738', a=6397300, b=6363806.28272251, f_=191, f=0.0052356, f2=0.00526316, n=0.00262467, e=0.10219488, e2=0.01044379, e21=0.98955621, e22=0.01055402, e32=0.00524931, A=6380564.13011837, L=10022566.69846922, R1=6386135.42757417, R2=6386131.54144847, R3=6386115.8862823, Rbiaxial=6380575.11882818, Rtriaxial=6388943.03218495)

45@var Ellipsoids.Mercury1960: Ellipsoid(name='Mercury1960', a=6378166, b=6356784.28360711, f_=298.3, f=0.00335233, f2=0.00336361, n=0.00167898, e=0.08181333, e2=0.00669342, e21=0.99330658, e22=0.00673853, e32=0.00335795, A=6367479.62923643, L=10002013.61254591, R1=6371038.76120237, R2=6371037.17115427, R3=6371030.78205124, Rbiaxial=6367484.1166614, Rtriaxial=6372827.29640037)

46@var Ellipsoids.Mercury1968Mod: Ellipsoid(name='Mercury1968Mod', a=6378150, b=6356768.33724438, f_=298.3, f=0.00335233, f2=0.00336361, n=0.00167898, e=0.08181333, e2=0.00669342, e21=0.99330658, e22=0.00673853, e32=0.00335795, A=6367463.65604381, L=10001988.52191361, R1=6371022.77908146, R2=6371021.18903735, R3=6371014.79995035, Rbiaxial=6367468.14345752, Rtriaxial=6372811.30979281)

47@var Ellipsoids.NWL1965: Ellipsoid(name='NWL1965', a=6378145, b=6356759.76948868, f_=298.25, f=0.00335289, f2=0.00336417, n=0.00167926, e=0.08182018, e2=0.00669454, e21=0.99330546, e22=0.00673966, e32=0.00335851, A=6367456.87366841, L=10001977.86818326, R1=6371016.58982956, R2=6371014.999254, R3=6371008.60802667, Rbiaxial=6367461.36258457, Rtriaxial=6372805.42010473)

48@var Ellipsoids.OSU86F: Ellipsoid(name='OSU86F', a=6378136.2, b=6356751.51693008, f_=298.2572236, f=0.00335281, f2=0.00336409, n=0.00167922, e=0.08181919, e2=0.00669438, e21=0.99330562, e22=0.0067395, e32=0.00335843, A=6367448.3471653, L=10001964.47478349, R1=6371007.97231003, R2=6371006.38181364, R3=6370999.99090513, Rbiaxial=6367452.83585765, Rtriaxial=6372796.75662978)

49@var Ellipsoids.OSU91A: Ellipsoid(name='OSU91A', a=6378136.3, b=6356751.6165948, f_=298.2572236, f=0.00335281, f2=0.00336409, n=0.00167922, e=0.08181919, e2=0.00669438, e21=0.99330562, e22=0.0067395, e32=0.00335843, A=6367448.44699773, L=10001964.63159991, R1=6371008.07219827, R2=6371006.48170186, R3=6371000.09079324, Rbiaxial=6367452.93569015, Rtriaxial=6372796.85654607)

50@var Ellipsoids.Plessis1817: Ellipsoid(name='Plessis1817', a=6376523, b=6355862.93325557, f_=308.64, f=0.00324002, f2=0.00325055, n=0.00162264, e=0.08043347, e2=0.00646954, e21=0.99353046, e22=0.00651167, e32=0.00324527, A=6366197.15710739, L=9999999.11003639, R1=6369636.31108519, R2=6369634.82608583, R3=6369628.85999668, Rbiaxial=6366201.34758009, Rtriaxial=6371364.26393357)

51@var Ellipsoids.PZ90: Ellipsoid(name='PZ90', a=6378136, b=6356751.36174571, f_=298.2578393, f=0.0033528, f2=0.00336408, n=0.00167922, e=0.08181911, e2=0.00669437, e21=0.99330563, e22=0.00673948, e32=0.00335842, A=6367448.16955443, L=10001964.19579298, R1=6371007.78724857, R2=6371006.1967588, R3=6370999.80587691, Rbiaxial=6367452.65822809, Rtriaxial=6372796.56780569)

52@var Ellipsoids.SGS85: Ellipsoid(name='SGS85', a=6378136, b=6356751.30156878, f_=298.257, f=0.00335281, f2=0.00336409, n=0.00167922, e=0.08181922, e2=0.00669438, e21=0.99330562, e22=0.0067395, e32=0.00335843, A=6367448.13949125, L=10001964.14856985, R1=6371007.76718959, R2=6371006.17669087, R3=6370999.78577297, Rbiaxial=6367452.62819019, Rtriaxial=6372796.55279934)

53@var Ellipsoids.SoAmerican1969: Ellipsoid(name='SoAmerican1969', a=6378160, b=6356774.71919531, f_=298.25, f=0.00335289, f2=0.00336417, n=0.00167926, e=0.08182018, e2=0.00669454, e21=0.99330546, e22=0.00673966, e32=0.00335851, A=6367471.84853228, L=10002001.39064442, R1=6371031.5730651, R2=6371029.98248581, R3=6371023.59124344, Rbiaxial=6367476.337459, Rtriaxial=6372820.40754721)

54@var Ellipsoids.Sphere: Ellipsoid(name='Sphere', a=6371008.771415, b=6371008.771415, f_=0, f=0, f2=0, n=0, e=0, e2=0, e21=1, e22=0, e32=0, A=6371008.771415, L=10007557.17611675, R1=6371008.771415, R2=6371008.771415, R3=6371008.771415, Rbiaxial=6371008.771415, Rtriaxial=6371008.771415)

55@var Ellipsoids.SphereAuthalic: Ellipsoid(name='SphereAuthalic', a=6371000, b=6371000, f_=0, f=0, f2=0, n=0, e=0, e2=0, e21=1, e22=0, e32=0, A=6371000, L=10007543.39801029, R1=6371000, R2=6371000, R3=6371000, Rbiaxial=6371000, Rtriaxial=6371000)

56@var Ellipsoids.SpherePopular: Ellipsoid(name='SpherePopular', a=6378137, b=6378137, f_=0, f=0, f2=0, n=0, e=0, e2=0, e21=1, e22=0, e32=0, A=6378137, L=10018754.17139462, R1=6378137, R2=6378137, R3=6378137, Rbiaxial=6378137, Rtriaxial=6378137)

57@var Ellipsoids.Struve1860: Ellipsoid(name='Struve1860', a=6378298.3, b=6356657.14266956, f_=294.73, f=0.00339294, f2=0.00340449, n=0.00169935, e=0.0823065, e2=0.00677436, e21=0.99322564, e22=0.00682056, e32=0.00339869, A=6367482.31832549, L=10002017.83655714, R1=6371084.58088985, R2=6371082.95208988, R3=6371076.40691418, Rbiaxial=6367486.91530791, Rtriaxial=6372894.90029454)

58@var Ellipsoids.WGS60: Ellipsoid(name='WGS60', a=6378165, b=6356783.28695944, f_=298.3, f=0.00335233, f2=0.00336361, n=0.00167898, e=0.08181333, e2=0.00669342, e21=0.99330658, e22=0.00673853, e32=0.00335795, A=6367478.63091189, L=10002012.0443814, R1=6371037.76231981, R2=6371036.17227197, R3=6371029.78316994, Rbiaxial=6367483.11833616, Rtriaxial=6372826.29723739)

59@var Ellipsoids.WGS66: Ellipsoid(name='WGS66', a=6378145, b=6356759.76948868, f_=298.25, f=0.00335289, f2=0.00336417, n=0.00167926, e=0.08182018, e2=0.00669454, e21=0.99330546, e22=0.00673966, e32=0.00335851, A=6367456.87366841, L=10001977.86818326, R1=6371016.58982956, R2=6371014.999254, R3=6371008.60802667, Rbiaxial=6367461.36258457, Rtriaxial=6372805.42010473)

60@var Ellipsoids.WGS72: Ellipsoid(name='WGS72', a=6378135, b=6356750.52001609, f_=298.26, f=0.00335278, f2=0.00336406, n=0.0016792, e=0.08181881, e2=0.00669432, e21=0.99330568, e22=0.00673943, e32=0.0033584, A=6367447.24862383, L=10001962.74919858, R1=6371006.84000536, R2=6371005.24953886, R3=6370998.8587507, Rbiaxial=6367451.7372317, Rtriaxial=6372795.60727472)

61@var Ellipsoids.WGS84: Ellipsoid(name='WGS84', a=6378137, b=6356752.31424518, f_=298.25722356, f=0.00335281, f2=0.00336409, n=0.00167922, e=0.08181919, e2=0.00669438, e21=0.99330562, e22=0.0067395, e32=0.00335843, A=6367449.14582341, L=10001965.72931272, R1=6371008.77141506, R2=6371007.18091847, R3=6371000.79000916, Rbiaxial=6367453.63451633, Rtriaxial=6372797.5559594)

62'''

63# make sure int/int division yields float quotient, see .basics

64 from __future__ import division as _; del _ # PYCHOK semicolon

65

66 from pygeodesy.basics import copysign0, isbool, isint, _xinstanceof

67 from pygeodesy.constants import EPS, EPS0, EPS02, EPS1, INF, NINF, PI4, PI_2, PI_3, R_M, R_MA, R_FM, \

68 _EPSqrt, _EPStol as _TOL, _floatuple as _T, _isfinite, _SQRT2_2, \

69 _0_0s, _0_0, _0_5, _1_0, _1_EPS, _2_0, _4_0, _90_0, \

70 _0_25, _3_0 # PYCHOK used!

71 from pygeodesy.errors import _AssertionError, IntersectionError, _ValueError, _xkwds_not

72 from pygeodesy.fmath import cbrt, cbrt2, fdot, Fhorner, fpowers, Fsum, hypot, hypot_, \

73 hypot1, hypot2, sqrt3

74# from pygeodesy.fsums import Fsum # from .fmath

75 from pygeodesy.interns import NN, _a_, _Airy1830_, _AiryModified_, _b_, _Bessel1841_, _beta_, \

76 _Clarke1866_, _Clarke1880IGN_, _DOT_, _f_, _GRS80_, _height_, \

77 _Intl1924_, _incompatible_, _invalid_, _Krassovski1940_, \

78 _Krassowsky1940_, _meridional_, _lat_, _negative_, _not_finite_, \

79 _prime_vertical_, _radius_, _Sphere_, _SPACE_, _vs_, _WGS72_, \

80 _WGS84_

81 from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS

82 from pygeodesy.named import _lazyNamedEnumItem as _lazy, _NamedEnum, \

83 _NamedEnumItem, _NamedTuple, _Pass

84 from pygeodesy.namedTuples import Distance2Tuple, Vector3Tuple, Vector4Tuple

85 from pygeodesy.props import deprecated_Property_RO, Property_RO, property_doc_, property_RO

86 from pygeodesy.streprs import Fmt, fstr, instr, strs, unstr

87 from pygeodesy.units import Bearing_, Distance, Float, Float_, Height, Lam_, Lat, Meter, \

88 Meter2, Meter3, Phi, Phi_, Radius, Radius_, Scalar

89 from pygeodesy.utily import atand, atan2b, atan2d, degrees90, m2radians, radians2m, sincos2d

90

91 from math import asinh, atan, atanh, cos, degrees, exp, fabs, radians, sin, sinh, sqrt, tan

92

93__all__ = _ALL_LAZY.ellipsoids

94__version__ = '23.08.04'

95

96_f_0_0 = Float(f =_0_0) # zero flattening

97_f__0_0 = Float(f_=_0_0) # zero inverse flattening

98# see U{WGS84_f<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Constants.html>}

99_f__WGS84 = Float(f_=_1_0 / (1000000000 / 298257223563)) # 298.25722356299997 vs 298.257223563

100

101

102 def _aux(lat, inverse, auxLat, clip=90):

103 '''Return a named auxiliary latitude in C{degrees}.

104 '''

105 return Lat(lat, clip=clip, name=_lat_ if inverse else auxLat.__name__)

106

107

108 def _s2_c2(phi):

109 '''(INTERNAL) Return 2-tuple C{(sin(B{phi})**2, cos(B{phi})**2)}.

110 '''

111 if phi:

112 s2 = sin(phi)**2

113 if s2 > EPS:

114 c2 = _1_0 - s2

115 if c2 > EPS:

116 if c2 < EPS1:

117 return s2, c2

118 else:

119 return _1_0, _0_0 # phi == PI_2

120 return _0_0, _1_0 # phi == 0

121

122

123 class a_f2Tuple(_NamedTuple):

124 '''2-Tuple C{(a, f)} specifying an ellipsoid by I{equatorial}

125 radius C{a} in C{meter} and scalar I{flattening} C{f}.

126

127 @see: Class L{Ellipsoid2}.

128 '''

129 _Names_ = (_a_, _f_) # name 'f' not 'f_'

130 _Units_ = (_Pass, _Pass)

131

132 def __new__(cls, a, f, **name):

133 '''New L{a_f2Tuple} ellipsoid specification.

134

135 @arg a: Equatorial radius (C{scalar} > 0).

136 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

137 @kwarg name: Optional name (C{str}).

138

139 @return: An L{a_f2Tuple}C{(a, f)} instance.

140

141 @raise UnitError: Invalid B{C{a}} or B{C{f}}.

142

143 @note: C{abs(B{f}) < EPS} is forced to C{B{f}=0}, I{spherical}.

144 Negative C{B{f}} produces a I{prolate} ellipsoid.

145 '''

146 a = Radius_(a=a) # low=EPS, high=None

147 f = Float_( f=f, low=None, high=EPS1)

148 if fabs(f) < EPS: # force spherical

149 f = _f_0_0

150 return _NamedTuple.__new__(cls, a, f, **name)

151

152 @Property_RO

153 def b(self):

154 '''Get the I{polar} radius (C{meter}), M{a * (1 - f)}.

155 '''

156 return a_f2b(self.a, self.f) # PYCHOK .a and .f

157

158 def ellipsoid(self, name=NN):

159 '''Return an L{Ellipsoid} for this 2-tuple C{(a, f)}.

160

161 @raise NameError: A registered C{ellipsoid} with the

162 same B{C{name}} already exists.

163 '''

164 return Ellipsoid(self.a, f=self.f, name=name or self.name) # PYCHOK .a and .f

165

166 @Property_RO

167 def f_(self):

168 '''Get the I{inverse} flattening (C{float}), M{1 / f} == M{a / (a - b)}.

169 '''

170 return f2f_(self.f) # PYCHOK .f

171

172

173 class Circle4Tuple(_NamedTuple):

174 '''4-Tuple C{(radius, height, lat, beta)} of the C{radius} and C{height},

175 both conventionally in C{meter} of a parallel I{circle of latitude} at

176 (geodetic) latitude C{lat} and the I{parametric (or reduced) auxiliary

177 latitude} C{beta}, both in C{degrees90}.

178

179 The C{height} is the (signed) distance along the z-axis between the

180 parallel and the equator. At near-polar C{lat}s, the C{radius} is C{0},

181 the C{height} is the ellipsoid's (signed) polar radius and C{beta}

182 equals C{lat}.

183 '''

184 _Names_ = (_radius_, _height_, _lat_, _beta_)

185 _Units_ = ( Radius, Height, Lat, Lat)

186

187

188 class Curvature2Tuple(_NamedTuple):

189 '''2-Tuple C{(meridional, prime_vertical)} of radii of curvature, both in

190 C{meter}, conventionally.

191 '''

192 _Names_ = (_meridional_, _prime_vertical_)

193 _Units_ = ( Meter, Meter)

194

195 @property_RO

196 def transverse(self):

197 '''Get this I{prime_vertical}, aka I{transverse} radius of curvature.

198 '''

199 return self.prime_vertical

200

201

202 class Ellipsoid(_NamedEnumItem):

203 '''Ellipsoid with I{equatorial} and I{polar} radii, I{flattening}, I{inverse

204 flattening} and other, often used, I{cached} attributes, supporting

205 I{oblate} and I{prolate} ellipsoidal and I{spherical} earth models.

206 '''

207 _a = 0 # equatorial radius, semi-axis (C{meter})

208 _b = 0 # polar radius, semi-axis (C{meter}): a * (f - 1) / f

209 _f = 0 # (1st) flattening: (a - b) / a

210 _f_ = 0 # inverse flattening: 1 / f = a / (a - b)

211

212 _geodsolve = NN # means, use PYGEODESY_GEODSOLVE

213 _KsOrder = 8 # Krüger series order (4, 6 or 8)

214 _rhumbsolve = NN # means, use PYGEODESY_RHUMBSOLVE

215

216 def __init__(self, a, b=None, f_=None, f=None, name=NN):

217 '''New L{Ellipsoid} from the I{equatorial} radius and either the

218 I{polar} radius or I{inverse flattening} or I{flattening}.

219

220 @arg a: Equatorial radius, semi-axis (C{meter}).

221 @arg b: Optional polar radius, semi-axis (C{meter}).

222 @arg f_: Inverse flattening: M{a / (a - b)} (C{float} >>> 1.0).

223 @arg f: Flattening: M{(a - b) / a} (C{float}, near zero for

224 spherical).

225 @kwarg name: Optional, unique name (C{str}).

226

227 @raise NameError: Ellipsoid with the same B{C{name}} already exists.

228

229 @raise ValueError: Invalid B{C{a}}, B{C{b}}, B{C{f_}} or B{C{f}} or

230 B{C{f_}} and B{C{f}} are incompatible.

231

232 @note: M{abs(f_) > 1 / EPS} or M{abs(1 / f_) < EPS} is forced

233 to M{1 / f_ = 0}, spherical.

234 '''

235 ff_ = f, f_ # assertion below

236 try:

237 a = Radius_(a=a) # low=EPS

238 if not _isfinite(a):

239 raise ValueError(_SPACE_(_a_, _not_finite_))

240

241 if b: # not in (_0_0, None)

242 b = Radius_(b=b) # low=EPS

243 f = a_b2f(a, b) if f is None else Float(f=f)

244 f_ = f2f_(f) if f_ is None else Float(f_=f_)

245 elif f is not None:

246 f = Float(f=f)

247 b = a_f2b(a, f)

248 f_ = f2f_(f) if f_ is None else Float(f_=f_)

249 elif f_:

250 f_ = Float(f_=f_)

251 b = a_f_2b(a, f_) # a * (f_ - 1) / f_

252 f = f_2f(f_)

253 else: # only a, spherical

254 f_ = f = 0

255 b = a # superfluous

256

257 if not f < _1_0: # sanity check, see .ecef.Ecef.__init__

258 raise ValueError(_SPACE_(_f_, _invalid_))

259 if not _isfinite(b):

260 raise ValueError(_SPACE_(_b_, _not_finite_))

261

262 if fabs(f) < EPS or a == b or not f_: # spherical

263 b = a

264 f = _f_0_0

265 f_ = _f__0_0

266

267 except (TypeError, ValueError) as x:

268 d = _xkwds_not(None, b=b, f_=f_, f=f)

269 t = instr(self, a=a, name=name, **d)

270 raise _ValueError(t, cause=x)

271

272 self._a = a

273 self._b = b

274 self._f = f

275 self._f_ = f_

276

277 self._register(Ellipsoids, name)

278

279 if f and f_: # see .test/testEllipsoidal.py

280 d = dict(eps=_TOL)

281 if None in ff_: # both f_ and f given

282 d.update(Error=_ValueError, txt=_incompatible_)

283 self._assert(_1_0 / f, f_=f_, **d)

284 self._assert(_1_0 / f_, f =f, **d)

285 self._assert(self.b2_a2, e21=self.e21, eps=EPS)

286

287 def __eq__(self, other):

288 '''Compare this and an other ellipsoid.

289

290 @arg other: The other ellipsoid (L{Ellipsoid} or L{Ellipsoid2}).

291

292 @return: C{True} if equal, C{False} otherwise.

293 '''

294 return self is other or (isinstance(other, Ellipsoid) and

295 self.a == other.a and

296 (self.f == other.f or self.b == other.b))

297

298 def __hash__(self):

299 return self._hash # memoized

300

301 @Property_RO

302 def a(self):

303 '''Get the I{equatorial} radius, semi-axis (C{meter}).

304 '''

305 return self._a

306

307 equatoradius = a # = Requatorial

308

309 @Property_RO

310 def a2(self):

311 '''Get the I{equatorial} radius I{squared} (C{meter} I{squared}), M{a**2}.

312 '''

313 return Meter2(a2=self.a**2)

314

315 @Property_RO

316 def a2_(self):

317 '''Get the inverse of the I{equatorial} radius I{squared} (C{meter} I{squared}), M{1 / a**2}.

318 '''

319 return Float(a2_=_1_0 / self.a2)

320

321 @Property_RO

322 def a_b(self):

323 '''Get the ratio I{equatorial} over I{polar} radius (C{float}), M{a / b} == M{1 / (1 - f)}.

324 '''

325 return Float(a_b=self.a / self.b if self.f else _1_0)

326

327 @Property_RO

328 def a2_b(self):

329 '''Get the I{polar} meridional (or polar) radius of curvature (C{meter}), M{a**2 / b}.

330

331 @see: U{Radii of Curvature

332 <https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>}

333 and U{Moritz, H. (1980), Geodetic Reference System 1980

334 <https://WikiPedia.org/wiki/Earth_radius#cite_note-Moritz-2>}.

335

336 @note: Symbol C{c} is used by IUGG and IERS for the U{polar radius of curvature

337 <https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>}, see L{c2}

338 and L{R2} or L{Rauthalic}.

339 '''

340 return Radius(a2_b=self.a2 / self.b if self.f else self.a) # = rocPolar

341

342 @Property_RO

343 def a2_b2(self):

344 '''Get the ratio I{equatorial} over I{polar} radius I{squared} (C{float}),

345 M{(a / b)**2} == M{1 / (1 - e**2)} == M{1 / (1 - e2)} == M{1 / e21}.

346 '''

347 return Float(a2_b2=self.a_b**2 if self.f else _1_0)

348

349 @Property_RO

350 def a_f(self):

351 '''Get the I{equatorial} radius and I{flattening} (L{a_f2Tuple}), see method C{toEllipsoid2}.

352 '''

353 return a_f2Tuple(self.a, self.f, name=self.name)

354

355 @Property_RO

356 def A(self):

357 '''Get the UTM I{meridional (or rectifying)} radius (C{meter}).

358

359 @see: I{Meridian arc unit} U{Q<https://StudyLib.net/doc/7443565/>}.

360 '''

361 A, n = self.a, self.n

362 if n:

363 d = (n + _1_0) * 1048576 / A

364 if d: # use 6 n**2 terms, half-way between the _KsOrder's 4, 6, 8

365 # <https://GeographicLib.SourceForge.io/C++/doc/tmseries30.html>

366 # <https://GeographicLib.SourceForge.io/C++/doc/transversemercator.html> and

367 # <https://www.MyGeodesy.id.AU/documents/Karney-Krueger%20equations.pdf> (3)

368 # A *= fhorner(n**2, 1048576, 262144, 16384, 4096, 1600, 784, 441) / 1048576) / (1 + n)

369 A = Radius(A=Fhorner(n**2, 1048576, 262144, 16384, 4096, 1600, 784, 441).fover(d))

370 return A

371

372 @Property_RO

373 def _albersCyl(self):

374 '''(INTERNAL) Helper for C{auxAuthalic}.

375 '''

376 return _MODS.albers.AlbersEqualAreaCylindrical(datum=self, name=self.name)

377

378 @Property_RO

379 def AlphaKs(self):

380 '''Get the I{Krüger} U{Alpha series coefficients<https://GeographicLib.SourceForge.io/C++/doc/tmseries30.html>} (C{KsOrder}C{-tuple}).

381 '''

382 return self._Kseries( # XXX int/int quotients may require from __future__ import division as _; del _ # PYCHOK semicolon

383 # n n**2 n**3 n**4 n**5 n**6 n**7 n**8

384 _T(1/2, -2/3, 5/16, 41/180, -127/288, 7891/37800, 72161/387072, -18975107/50803200),

385 _T(13/48, -3/5, 557/1440, 281/630, -1983433/1935360, 13769/28800, 148003883/174182400), # PYCHOK unaligned

386 _T(61/240, -103/140, 15061/26880, 167603/181440, -67102379/29030400, 79682431/79833600), # PYCHOK unaligned

387 _T(49561/161280, -179/168, 6601661/7257600, 97445/49896, -40176129013/7664025600), # PYCHOK unaligned

388 _T(34729/80640, -3418889/1995840, 14644087/9123840, 2605413599/622702080), # PYCHOK unaligned

389 _T(212378941/319334400, -30705481/10378368, 175214326799/58118860800), # PYCHOK unaligned

390 _T(1522256789/1383782400, -16759934899/3113510400), # PYCHOK unaligned

391 _T(1424729850961/743921418240)) # PYCHOK unaligned

392

393 @Property_RO

394 def area(self):

395 '''Get the ellipsoid's surface area (C{meter} I{squared}), M{4 * PI * c2}.

396

397 @see: Properties L{areax}, L{c2} and L{R2} and functions

398 L{ellipsoidalExact.areaOf} and L{ellipsoidalKarney.areaOf}.

399 '''

400 return Meter2(area=self.c2 * PI4)

401

402 @Property_RO

403 def areax(self):

404 '''Get the ellipsoid's surface area (C{meter} I{squared}), M{4 * PI * c2x}, more

405 accurate for very I{oblate} ellipsoids.

406

407 @see: Properties L{area}, L{c2x} and L{R2x}, class L{GeodesicExact} and

408 functions L{ellipsoidalExact.areaOf} and L{ellipsoidalKarney.areaOf}.

409 '''

410 return Meter2(areax=self.c2x * PI4)

411

412 def _assert(self, val, eps=_TOL, f0=_0_0, Error=_AssertionError, txt=NN, **name_value):

413 '''(INTERNAL) Assert a C{name=value} vs C{val}.

414 '''

415 for n, v in name_value.items():

416 if fabs(v - val) > eps: # PYCHOK no cover

417 t = (v, _vs_, val)

418 t = _SPACE_.join(strs(t, prec=12, fmt=Fmt.g))

419 t = Fmt.EQUAL(self._DOT_(n), t)

420 raise Error(t, txt=txt or Fmt.exceeds_eps(eps))

421 return Float(v if self.f else f0, name=n)

422 raise Error(unstr(self._DOT_(self._assert.__name__), val,

423 eps=eps, f0=f0, **name_value))

424

425 def auxAuthalic(self, lat, inverse=False):

426 '''Compute the I{authalic} auxiliary latitude or the I{inverse} thereof.

427

428 @arg lat: The geodetic (or I{authalic}) latitude (C{degrees90}).

429 @kwarg inverse: If C{True}, B{C{lat}} is the I{authalic} and

430 return the geodetic latitude (C{bool}).

431

432 @return: The I{authalic} (or geodetic) latitude in C{degrees90}.

433

434 @see: U{Inverse-/AuthalicLatitude<https://GeographicLib.SourceForge.io/

435 html/classGeographicLib_1_1Ellipsoid.html>}, U{Authalic latitude

436 <https://WikiPedia.org/wiki/Latitude#Authalic_latitude>}, and

437 U{Snyder<https://Pubs.USGS.gov/pp/1395/report.pdf>}, p 16.

438 '''

439 if self.f:

440 f = self._albersCyl._tanf if inverse else self._albersCyl._txif # PYCHOK attr

441 lat = atand(f(tan(Phi_(lat)))) # PYCHOK attr

442 return _aux(lat, inverse, Ellipsoid.auxAuthalic)

443

444 def auxConformal(self, lat, inverse=False):

445 '''Compute the I{conformal} auxiliary latitude or the I{inverse} thereof.

446

447 @arg lat: The geodetic (or I{conformal}) latitude (C{degrees90}).

448 @kwarg inverse: If C{True}, B{C{lat}} is the I{conformal} and

449 return the geodetic latitude (C{bool}).

450

451 @return: The I{conformal} (or geodetic) latitude in C{degrees90}.

452

453 @see: U{Inverse-/ConformalLatitude<https://GeographicLib.SourceForge.io/

454 html/classGeographicLib_1_1Ellipsoid.html>}, U{Conformal latitude

455 <https://WikiPedia.org/wiki/Latitude#Conformal_latitude>}, and

456 U{Snyder<https://Pubs.USGS.gov/pp/1395/report.pdf>}, pp 15-16.

457 '''

458 if self.f:

459 f = self.es_tauf if inverse else self.es_taupf # PYCHOK attr

460 lat = atand(f(tan(Phi_(lat)))) # PYCHOK attr

461 return _aux(lat, inverse, Ellipsoid.auxConformal)

462

463 def auxGeocentric(self, lat, inverse=False):

464 '''Compute the I{geocentric} auxiliary latitude or the I{inverse} thereof.

465

466 @arg lat: The geodetic (or I{geocentric}) latitude (C{degrees90}).

467 @kwarg inverse: If C{True}, B{C{lat}} is the geocentric and

468 return the I{geocentric} latitude (C{bool}).

469

470 @return: The I{geocentric} (or geodetic) latitude in C{degrees90}.

471

472 @see: U{Inverse-/GeocentricLatitude<https://GeographicLib.SourceForge.io/

473 html/classGeographicLib_1_1Ellipsoid.html>}, U{Geocentric latitude

474 <https://WikiPedia.org/wiki/Latitude#Geocentric_latitude>}, and

475 U{Snyder<https://Pubs.USGS.gov/pp/1395/report.pdf>}, pp 17-18.

476 '''

477 if self.f:

478 f = self.a2_b2 if inverse else self.b2_a2

479 lat = atand(f * tan(Phi_(lat)))

480 return _aux(lat, inverse, Ellipsoid.auxGeocentric)

481

482 def auxIsometric(self, lat, inverse=False):

483 '''Compute the I{isometric} auxiliary latitude or the I{inverse} thereof.

484

485 @arg lat: The geodetic (or I{isometric}) latitude (C{degrees}).

486 @kwarg inverse: If C{True}, B{C{lat}} is the I{isometric} and

487 return the geodetic latitude (C{bool}).

488

489 @return: The I{isometric} (or geodetic) latitude in C{degrees}.

490

491 @note: The I{isometric} latitude for geodetic C{+/-90} is far

492 outside the C{[-90..+90]} range but the inverse

493 thereof is the original geodetic latitude.

494

495 @see: U{Inverse-/IsometricLatitude<https://GeographicLib.SourceForge.io/

496 html/classGeographicLib_1_1Ellipsoid.html>}, U{Isometric latitude

497 <https://WikiPedia.org/wiki/Latitude#Isometric_latitude>}, and

498 U{Snyder<https://Pubs.USGS.gov/pp/1395/report.pdf>}, pp 15-16.

499 '''

500 if self.f:

501 r = Phi_(lat, clip=0)

502 lat = degrees(atan(self.es_tauf(sinh(r))) if inverse else

503 asinh(self.es_taupf(tan(r))))

504 # clip=0, since auxIsometric(+/-90) is far outside [-90..+90]

505 return _aux(lat, inverse, Ellipsoid.auxIsometric, clip=0)

506

507 def auxParametric(self, lat, inverse=False):

508 '''Compute the I{parametric} auxiliary latitude or the I{inverse} thereof.

509

510 @arg lat: The geodetic (or I{parametric}) latitude (C{degrees90}).

511 @kwarg inverse: If C{True}, B{C{lat}} is the I{parametric} and

512 return the geodetic latitude (C{bool}).

513

514 @return: The I{parametric} (or geodetic) latitude in C{degrees90}.

515

516 @see: U{Inverse-/ParametricLatitude<https://GeographicLib.SourceForge.io/

517 html/classGeographicLib_1_1Ellipsoid.html>}, U{Parametric latitude

518 <https://WikiPedia.org/wiki/Latitude#Parametric_(or_reduced)_latitude>},

519 and U{Snyder<https://Pubs.USGS.gov/pp/1395/report.pdf>}, p 18.

520 '''

521 if self.f:

522 lat = self._beta(Lat(lat), inverse=inverse)

523 return _aux(lat, inverse, Ellipsoid.auxParametric)

524

525 auxReduced = auxParametric # synonymous

526

527 def auxRectifying(self, lat, inverse=False):

528 '''Compute the I{rectifying} auxiliary latitude or the I{inverse} thereof.

529

530 @arg lat: The geodetic (or I{rectifying}) latitude (C{degrees90}).

531 @kwarg inverse: If C{True}, B{C{lat}} is the I{rectifying} and

532 return the geodetic latitude (C{bool}).

533

534 @return: The I{rectifying} (or geodetic) latitude in C{degrees90}.

535

536 @see: U{Inverse-/RectifyingLatitude<https://GeographicLib.SourceForge.io/

537 html/classGeographicLib_1_1Ellipsoid.html>}, U{Rectifying latitude

538 <https://WikiPedia.org/wiki/Latitude#Rectifying_latitude>}, and

539 U{Snyder<https://Pubs.USGS.gov/pp/1395/report.pdf>}, pp 16-17.

540 '''

541 if self.f:

542 lat = Lat(lat)

543 if 0 < fabs(lat) < _90_0:

544 if inverse:

545 e = self._elliptic_e22

546 d = degrees90(e.fEinv(e.cE * lat / _90_0))

547 lat = self.auxParametric(d, inverse=True)

548 else:

549 lat = _90_0 * self.Llat(lat) / self.L

550 return _aux(lat, inverse, Ellipsoid.auxRectifying)

551

552 @Property_RO

553 def b(self):

554 '''Get the I{polar} radius, semi-axis (C{meter}).

555 '''

556 return self._b

557

558 polaradius = b # = Rpolar

559

560 @Property_RO

561 def b_a(self):

562 '''Get the ratio I{polar} over I{equatorial} radius (C{float}), M{b / a == f1 == 1 - f}.

563

564 @see: Property L{f1}.

565 '''

566 return self._assert(self.b / self.a, b_a=self.f1, f0=_1_0)

567

568 @Property_RO

569 def b2(self):

570 '''Get the I{polar} radius I{squared} (C{float}), M{b**2}.

571 '''

572 return Meter2(b2=self.b**2)

573

574 @Property_RO

575 def b2_a(self):

576 '''Get the I{equatorial} meridional radius of curvature (C{meter}), M{b**2 / a}, see C{rocMeridional}C{(0)}.

577

578 @see: U{Radii of Curvature<https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>}.

579 '''

580 return Radius(b2_a=self.b2 / self.a if self.f else self.b)

581

582 @Property_RO

583 def b2_a2(self):

584 '''Get the ratio I{polar} over I{equatorial} radius I{squared} (C{float}), M{(b / a)**2}

585 == M{(1 - f)**2} == M{1 - e**2} == C{e21}.

586 '''

587 return Float(b2_a2=self.b_a**2 if self.f else _1_0)

588

589 def _beta(self, lat, inverse=False):

590 '''(INTERNAL) Get the I{parametric (or reduced) auxiliary latitude} or inverse thereof.

591 '''

592 s, c = sincos2d(lat) # like Karney's tand(lat)

593 s *= self.a_b if inverse else self.b_a

594 return atan2d(s, c) # == atand(s / c) if c else copysign0(_90_0, lat)

595

596 @Property_RO

597 def BetaKs(self):

598 '''Get the I{Krüger} U{Beta series coefficients<https://GeographicLib.SourceForge.io/C++/doc/tmseries30.html>} (C{KsOrder}C{-tuple}).

599 '''

600 return self._Kseries( # XXX int/int quotients may require from __future__ import division as _; del _ # PYCHOK semicolon

601 # n n**2 n**3 n**4 n**5 n**6 n**7 n**8

602 _T(1/2, -2/3, 37/96, -1/360, -81/512, 96199/604800, -5406467/38707200, 7944359/67737600),

603 _T(1/48, 1/15, -437/1440, 46/105, -1118711/3870720, 51841/1209600, 24749483/348364800), # PYCHOK unaligned

604 _T(17/480, -37/840, -209/4480, 5569/90720, 9261899/58060800, -6457463/17740800), # PYCHOK unaligned

605 _T(4397/161280, -11/504, -830251/7257600, 466511/2494800, 324154477/7664025600), # PYCHOK unaligned

606 _T(4583/161280, -108847/3991680, -8005831/63866880, 22894433/124540416), # PYCHOK unaligned

607 _T(20648693/638668800, -16363163/518918400, -2204645983/12915302400), # PYCHOK unaligne

608 _T(219941297/5535129600, -497323811/12454041600), # PYCHOK unaligned

609 _T(191773887257/3719607091200)) # PYCHOK unaligned

610

611 @deprecated_Property_RO

612 def c(self): # PYCHOK no cover

613 '''DEPRECATED, use property C{R2} or C{Rauthalic}.'''

614 return self.R2

615

616 @Property_RO

617 def c2(self):

618 '''Get the I{authalic} earth radius I{squared} (C{meter} I{squared}).

619

620 @see: Properties L{c2x}, L{area}, L{R2}, L{Rauthalic}, I{Karney's} U{equation (60)

621 <https://Link.Springer.com/article/10.1007%2Fs00190-012-0578-z>} and C++ U{Ellipsoid.Area

622 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Ellipsoid.html>},

623 U{Authalic radius<https://WikiPedia.org/wiki/Earth_radius#Authalic_radius>}, U{Surface area

624 <https://WikiPedia.org/wiki/Ellipsoid>} and U{surface area

625 <https://www.Numericana.com/answer/geometry.htm#oblate>}.

626 '''

627 return self._c2f(False)

628

629 @Property_RO

630 def c2x(self):

631 '''Get the I{authalic} earth radius I{squared} (C{meter} I{squared}), more accurate for very I{oblate}

632 ellipsoids.

633

634 @see: Properties L{c2}, L{areax}, L{R2x}, L{Rauthalicx}, class L{GeodesicExact} and I{Karney}'s comments at C++

635 attribute U{GeodesicExact._c2<https://GeographicLib.SourceForge.io/C++/doc/GeodesicExact_8cpp_source.html>}.

636 '''

637 return self._c2f(True)

638

639 def _c2f(self, c2x):

640 '''(INTERNAL) Helper for C{.c2} and C{.c2x}.

641 '''

642 f, c2 = self.f, self.b2

643 if f:

644 e = self.e

645 if e > EPS0:

646 if f > 0: # .isOblate

647 c2 *= (asinh(sqrt(self.e22abs)) if c2x else atanh(e)) / e

648 elif f < 0: # .isProlate

649 c2 *= atan(e) / e # XXX asin?

650 c2 = Meter2(c2=(self.a2 + c2) * _0_5)

651 return c2

652

653 def circle4(self, lat):

654 '''Get the equatorial or a parallel I{circle of latitude}.

655

656 @arg lat: Geodetic latitude (C{degrees90}, C{str}).

657

658 @return: A L{Circle4Tuple}C{(radius, height, lat, beta)}

659 instance.

660

661 @raise RangeError: Latitude B{C{lat}} outside valid range and

662 L{pygeodesy.rangerrors} set to C{True}.

663

664 @raise TypeError: Invalid B{C{lat}}.

665

666 @raise ValueError: Invalid B{C{lat}}.

667

668 @see: Definition of U{I{p} and I{z} under B{Parametric (or reduced) latitude}

669 <https://WikiPedia.org/wiki/Latitude>}, I{Karney's} C++ U{CircleRadius and CircleHeight

670 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Ellipsoid.html>}

671 and method C{Rlat}.

672 '''

673 lat = Lat(lat)

674 if lat:

675 b = lat

676 if fabs(lat) < _90_0:

677 if self.f:

678 b = self._beta(lat)

679 z, r = sincos2d(b)

680 r *= self.a

681 z *= self.b

682 else: # near-polar

683 r, z = _0_0, copysign0(self.b, lat)

684 else: # equator

685 r = self.a

686 z = lat = b = _0_0

687 return Circle4Tuple(r, z, lat, b)

688

689 def degrees2m(self, deg, lat=0):

690 '''Convert an angle to the distance along the equator or

691 along a parallel of (geodetic) latitude.

692

693 @arg deg: The angle (C{degrees}).

694 @kwarg lat: Parallel latitude (C{degrees90}, C{str}).

695

696 @return: Distance (C{meter}, same units as the equatorial

697 and polar radii) or C{0} for near-polar B{C{lat}}.

698

699 @raise RangeError: Latitude B{C{lat}} outside valid range and

700 L{pygeodesy.rangerrors} set to C{True}.

701

702 @raise ValueError: Invalid B{C{deg}} or B{C{lat}}.

703 '''

704 return self.radians2m(radians(deg), lat=lat)

705

706 def distance2(self, lat0, lon0, lat1, lon1):

707 '''I{Approximate} the distance and (initial) bearing between

708 two points based on the U{local, flat earth approximation

709 <https://www.EdWilliams.org/avform.htm#flat>} aka U{Hubeny

710 <https://www.OVG.AT/de/vgi/files/pdf/3781/>} formula.

711

712 I{Suitable only for distances of several hundred Km or Miles

713 and only between points not near-polar}.

714

715 @arg lat0: From latitude (C{degrees}).

716 @arg lon0: From longitude (C{degrees}).

717 @arg lat1: To latitude (C{degrees}).

718 @arg lon1: To longitude (C{degrees}).

719

720 @return: A L{Distance2Tuple}C{(distance, initial)} with C{distance}

721 in same units as this ellipsoid's axes.

722

723 @note: The meridional and prime_vertical radii of curvature are

724 taken and scaled I{at the initial latitude}, see C{roc2}.

725

726 @see: Function L{pygeodesy.flatLocal}/L{pygeodesy.hubeny}.

727 '''

728 phi0 = Phi_(lat0=lat0)

729 m, n = self.roc2_(phi0, scaled=True)

730 m *= Phi_(lat1=lat1) - phi0

731 n *= Lam_(lon1=lon1) - Lam_(lon0=lon0)

732 return Distance2Tuple(hypot(m, n), atan2b(n, m))

733

734 @Property_RO

735 def e(self):

736 '''Get the I{unsigned, (1st) eccentricity} (C{float}), M{sqrt(1 - (b / a)**2))}, see C{a_b2e}.

737

738 @see: Property L{es}.

739 '''

740 return Float(e=sqrt(self.e2abs) if self.e2 else _0_0)

741

742 @deprecated_Property_RO

743 def e12(self): # see property ._e12

744 '''DEPRECATED, use property C{e21}.'''

745 return self.e21

746

747# @Property_RO

748# def _e12(self): # see property ._elliptic_e12

749# # (INTERNAL) until e12 above can be replaced with e21.

750# return self.e2 / (_1_0 - self.e2) # see I{Karney}'s Ellipsoid._e12 = e2 / (1 - e2)

751

752 @Property_RO

753 def e2(self):

754 '''Get the I{signed, (1st) eccentricity squared} (C{float}), M{f * (2 - f)

755 == 1 - (b / a)**2}, see C{a_b2e2}.

756 '''

757 return self._assert(a_b2e2(self.a, self.b), e2=f2e2(self.f))

758

759 @Property_RO

760 def e2abs(self):

761 '''Get the I{unsigned, (1st) eccentricity squared} (C{float}).

762 '''

763 return fabs(self.e2)

764

765 @Property_RO

766 def e21(self):

767 '''Get 1 less I{1st eccentricity squared} (C{float}), M{1 - e**2}

768 == M{1 - e2} == M{(1 - f)**2} == M{b**2 / a**2}, see C{b2_a2}.

769 '''

770 return self._assert((_1_0 - self.f)**2, e21=_1_0 - self.e2, f0=_1_0)

771

772# _e2m = e21 # see I{Karney}'s Ellipsoid._e2m = 1 - _e2

773 _1_e21 = a2_b2 # == M{1 / e21} == M{1 / (1 - e**2)}

774

775 @Property_RO

776 def e22(self):

777 '''Get the I{signed, 2nd eccentricity squared} (C{float}), M{e2 / (1 - e2)

778 == e2 / (1 - f)**2 == (a / b)**2 - 1}, see C{a_b2e22}.

779 '''

780 return self._assert(a_b2e22(self.a, self.b), e22=f2e22(self.f))

781

782 @Property_RO

783 def e22abs(self):

784 '''Get the I{unsigned, 2nd eccentricity squared} (C{float}).

785 '''

786 return fabs(self.e22)

787

788 @Property_RO

789 def e32(self):

790 '''Get the I{signed, 3rd eccentricity squared} (C{float}), M{e2 / (2 - e2)

791 == (a**2 - b**2) / (a**2 + b**2)}, see C{a_b2e32}.

792 '''

793 return self._assert(a_b2e32(self.a, self.b), e32=f2e32(self.f))

794

795 @Property_RO

796 def e32abs(self):

797 '''Get the I{unsigned, 3rd eccentricity squared} (C{float}).

798 '''

799 return fabs(self.e32)

800

801 @Property_RO

802 def e4(self):

803 '''Get the I{unsignd, (1st) eccentricity} to 4th power (C{float}), M{e**4 == e2**2}.

804 '''

805 return Float(e4=self.e2**2 if self.e2 else _0_0)

806

807 eccentricity = e # eccentricity

808# eccentricity2 = e2 # eccentricity squared

809 eccentricity1st2 = e2 # first eccentricity squared

810 eccentricity2nd2 = e22 # second eccentricity squared

811 eccentricity3rd2 = e32 # third eccentricity squared

812

813 def ecef(self, Ecef=None):

814 '''Return U{ECEF<https://WikiPedia.org/wiki/ECEF>} converter.

815

816 @kwarg Ecef: ECEF class to use, default L{EcefKarney}.

817

818 @return: An ECEF converter for this C{ellipsoid}.

819

820 @raise TypeError: Invalid B{C{Ecef}}.

821

822 @see: Module L{pygeodesy.ecef}.

823 '''

824 return _MODS.ecef._4Ecef(self, Ecef)

825

826 @Property_RO

827 def _elliptic_e12(self): # see I{Karney}'s Ellipsoid._e12

828 '''(INTERNAL) Elliptic helper for C{rhumbx.Rhumb}.

829 '''

830 e12 = self.e2 / (_1_0 - self.e2) # NOT DEPRECATED .e12!

831 return _MODS.elliptic.Elliptic(-e12)

832

833 @Property_RO

834 def _elliptic_e22(self): # aka ._elliptic_ep2

835 '''(INTERNAL) Elliptic helper for C{auxRectifying}, C{L}, C{Llat}.

836 '''

837 return _MODS.elliptic.Elliptic(-self.e22abs) # complex

838

839 equatoradius = a # Requatorial

840

841 def e2s(self, s):

842 '''Compute norm M{sqrt(1 - e2 * s**2)}.

843

844 @arg s: Sine value (C{scalar}).

845

846 @return: Norm (C{float}).

847

848 @raise ValueError: Invalid B{C{s}}.

849 '''

850 return sqrt(self.e2s2(s)) if self.e2 else _1_0

851

852 def e2s2(self, s):

853 '''Compute M{1 - e2 * s**2}.

854

855 @arg s: Sine value (C{scalar}).

856

857 @return: Result (C{float}).

858

859 @raise ValueError: Invalid B{C{s}}.

860 '''

861 r = _1_0

862 if self.e2:

863 try:

864 r -= self.e2 * Scalar(s=s)**2

865 if r < 0:

866 raise ValueError(_negative_)

867 except (TypeError, ValueError) as x:

868 t = self._DOT_(Ellipsoid.e2s2.__name__)

869 raise _ValueError(t, s, cause=x)

870 return r

871

872 @Property_RO

873 def es(self):

874 '''Get the I{signed (1st) eccentricity} (C{float}).

875

876 @see: Property L{e}.

877 '''

878 # note, self.e is always non-negative

879 return Float(es=copysign0(self.e, self.f)) # see .ups

880

881 def es_atanh(self, x):

882 '''Compute M{es * atanh(es * x)} or M{-es * atan(es * x)}

883 for I{oblate} respectively I{prolate} ellipsoids where

884 I{es} is the I{signed} (1st) eccentricity.

885

886 @raise ValueError: Invalid B{C{x}}.

887

888 @see: Function U{Math::eatanhe<https://GeographicLib.SourceForge.io/

889 html/classGeographicLib_1_1Math.html>}.

890 '''

891 return self._es_atanh(Scalar(x=x)) if self.f else _0_0

892

893 def _es_atanh(self, x): # see .albers._atanhee, .AuxLat._atanhee

894 '''(INTERNAL) Helper for .es_atanh, ._es_taupf2 and ._exp_es_atanh.

895 '''

896 es = self.es # signOf(es) == signOf(f)

897 return es * (atanh(es * x) if es > 0 else # .isOblate

898 (-atan(es * x) if es < 0 else # .isProlate

899 _0_0)) # .isSpherical

900

901 @Property_RO

902 def es_c(self):

903 '''Get M{(1 - f) * exp(es_atanh(1))} (C{float}), M{b_a * exp(es_atanh(1))}.

904 '''

905 return Float(es_c=(self._exp_es_atanh_1 * self.b_a) if self.f else _1_0)

906

907 def es_tauf(self, taup):

908 '''Compute I{Karney}'s U{equations (19), (20) and (21)

909 <https://ArXiv.org/abs/1002.1417>}.

910

911 @see: I{Karney}'s C++ method U{Math::tauf<https://GeographicLib.

912 SourceForge.io/C++/doc/classGeographicLib_1_1Math.html>} and

913 and I{Veness}' JavaScript method U{toLatLon<https://www.

914 Movable-Type.co.UK/scripts/latlong-utm-mgrs.html>}.

915 '''

916 t = Scalar(taup=taup)

917 if self.f: # .isEllipsoidal

918 a = fabs(t)

919 T = t * (self._exp_es_atanh_1 if a > 70 else self._1_e21)

920 if fabs(T * _EPSqrt) < _2_0: # handles +/- INF and NAN

921 s = (a * _TOL) if a > _1_0 else _TOL

922 for T, _, d in self._es_tauf3(t, T): # max 2

923 if fabs(d) < s:

924 break

925 t = Scalar(tauf=T)

926 return t

927

928 def _es_tauf3(self, taup, T, N=9): # in .utm.Utm._toLLEB

929 '''(INTERNAL) Yield a 3-tuple C{(τi, iteration, delta)} for at most

930 B{C{N}} Newton iterations, converging rapidly except when C{delta}

931 toggles on +/-1.12e-16 or +/-4.47e-16, see C{.utm.Utm._toLLEB}.

932 '''

933 e = self._1_e21

934 _F2_ = Fsum(T).fsum2_ # τ0

935 _tf2 = self._es_taupf2

936 for i in range(1, N + 1):

937 a, h = _tf2(T)

938 d = (taup - a) * (e + T**2) / (hypot1(a) * h)

939 # = (taup - a) / hypot1(a) / ((e + T**2) / h)

940 T, d = _F2_(d) # τi, (τi - τi-1)

941 yield T, i, d

942

943 def es_taupf(self, tau):

944 '''Compute I{Karney}'s U{equations (7), (8) and (9)

945 <https://ArXiv.org/abs/1002.1417>}.

946

947 @see: I{Karney}'s C++ method U{Math::taupf<https://GeographicLib.

948 SourceForge.io/C++/doc/classGeographicLib_1_1Math.html>}.

949 '''

950 t = Scalar(tau=tau)

951 if self.f: # .isEllipsoidal

952 t, _ = self._es_taupf2(t)

953 t = Scalar(taupf=t)

954 return t

955

956 def _es_taupf2(self, tau):

957 '''(INTERNAL) Return 2-tuple C{(es_taupf(tau), hypot1(tau))}.

958 '''

959 if _isfinite(tau):

960 h = hypot1(tau)

961 s = sinh(self._es_atanh(tau / h))

962 a = hypot1(s) * tau - h * s

963 else:

964 a, h = tau, INF

965 return a, h

966

967 @Property_RO

968 def _exp_es_atanh_1(self):

969 '''(INTERNAL) Helper for .es_c and .es_tauf.

970 '''

971 return exp(self._es_atanh(_1_0)) if self.es else _1_0

972

973 @Property_RO

974 def f(self):

975 '''Get the I{flattening} (C{float}), M{(a - b) / a}, C{0} for spherical, negative for prolate.

976 '''

977 return self._f

978

979 @Property_RO

980 def f_(self):

981 '''Get the I{inverse flattening} (C{float}), M{1 / f} == M{a / (a - b)}, C{0} for spherical, see C{a_b2f_}.

982 '''

983 return self._f_

984

985 @Property_RO

986 def f1(self):

987 '''Get the I{1 - flattening} (C{float}), M{f1 == 1 - f == b / a}.

988

989 @see: Property L{b_a}.

990 '''

991 return Float(f1=_1_0 - self.f)

992

993 @Property_RO

994 def f2(self):

995 '''Get the I{2nd flattening} (C{float}), M{(a - b) / b == f / (1 - f)}, C{0} for spherical, see C{a_b2f2}.

996 '''

997 return self._assert(self.a_b - _1_0, f2=f2f2(self.f))

998

999 @deprecated_Property_RO

1000 def geodesic(self):

1001 '''DEPRECATED, use property C{geodesicw}.'''

1002 return self.geodesicw

1003

1004 def geodesic_(self, exact=True):

1005 '''Get the an I{exact} C{Geodesic...} instance for this ellipsoid.

1006

1007 @kwarg exact: If C{bool} return L{GeodesicExact}C{(exact=B{exact},

1008 ...)}, otherwise a C{geodesic} for this ellipsoid

1009 (C{.geodesicw}, C{.geodesicx}, C{.geodsolve}) or an

1010 other L{Geodesic}, L{GeodesicExact} or L{GeodesicSolve}

1011 instance for I{this} ellipsoid.

1012

1013 @return: The C{exact} geodesic (C{Geodesic...}).

1014

1015 @raise TypeError: Invalid B{C{exact}}.

1016

1017 @raise ValueError: Incompatible B{C{exact}} ellipsoid.

1018 '''

1019 if isbool(exact): # for consistenccy with C{.rhumb_}

1020 g = _MODS.geodesicx.GeodesicExact(self, C4order=30 if exact else 24,

1021 name=self.name)

1022 else:

1023 g = exact

1024 _xinstanceof(_MODS.geodesicw.Geodesic,

1025 _MODS.geodesicx.GeodesicExact,

1026 _MODS.geodsolve.GeodesicSolve, exact=g)

1027 if g.ellipsoid != self:

1028 raise _ValueError(exact=g, ellipsosid=g.ellipsoid)

1029 return g

1030

1031 @Property_RO

1032 def geodesicw(self):

1033 '''Get this ellipsoid's I{wrapped} U{geodesicw.Geodesic

1034 <https://GeographicLib.SourceForge.io/Python/doc/code.html>}, provided

1035 I{Karney}'s U{geographiclib<https://PyPI.org/project/geographiclib>}

1036 package is installed.

1037 '''

1038 # if not self.isEllipsoidal:

1039 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

1040 return _MODS.geodesicw.Geodesic(self)

1041

1042 @Property_RO

1043 def geodesicx(self):

1044 '''Get this ellipsoid's I{exact} L{GeodesicExact}.

1045 '''

1046 # if not self.isEllipsoidal:

1047 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

1048 return _MODS.geodesicx.GeodesicExact(self, name=self.name)

1049

1050 @property

1051 def geodsolve(self):

1052 '''Get this ellipsoid's L{GeodesicSolve}, the I{wrapper} around utility

1053 U{GeodSolve<https://GeographicLib.SourceForge.io/C++/doc/GeodSolve.1.html>},

1054 provided the path to the C{GeodSolve} executable is specified with env

1055 variable C{PYGEODESY_GEODSOLVE} or re-/set with this property..

1056 '''

1057 # if not self.isEllipsoidal:

1058 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

1059 return _MODS.geodsolve.GeodesicSolve(self, path=self._geodsolve, name=self.name)

1060

1061 @geodsolve.setter # PYCHOK setter!

1062 def geodsolve(self, path):

1063 '''Re-/set the (fully qualified) path to the U{GeodSolve

1064 <https://GeographicLib.SourceForge.io/C++/doc/GeodSolve.1.html>} executable,

1065 overriding env variable C{PYGEODESY_GEODSOLVE} (C{str}).

1066 '''

1067 self._geodsolve = path

1068

1069 def hartzell4(self, pov, los=None):

1070 '''Compute the intersection of this ellipsoid's surface and a Line-Of-Sight

1071 from a Point-Of-View in space.

1072

1073 @arg pov: Point-Of-View outside this ellipsoid (C{Cartesian}, L{Ecef9Tuple}

1074 or L{Vector3d}).

1075 @kwarg los: Line-Of-Sight, I{direction} to this ellipsoid (L{Vector3d}) or

1076 C{None} to point to this ellipsoid's center.

1077

1078 @return: L{Vector4Tuple}C{(x, y, z, h)} with the cartesian coordinates C{x},

1079 C{y} and C{z} of the projection on or the intersection with this

1080 ellipsoid and the I{distance} C{h} from B{C{pov}} to C{(x, y, z)}

1081 along B{C{los}}, all in C{meter}, conventionally.

1082

1083 @raise IntersectionError: Null B{C{pov}} or B{C{los}} vector, or B{C{pov}}

1084 is inside this ellipsoid or B{C{los}} points

1085 outside this ellipsoid or points in an opposite

1086 direction.

1087

1088 @raise TypeError: Invalid B{C{pov}} or B{C{los}}.

1089

1090 @see: U{I{Satellite Line-of-Sight Intersection with Earth}<https://StephenHartzell.

1091 Medium.com/satellite-line-of-sight-intersection-with-earth-d786b4a6a9b6>} and

1092 methods L{Ellipsoid.height4} and L{Triaxial.hartzell4}.

1093 '''

1094 try:

1095 v, d = _MODS.triaxials._hartzell3d2(pov, los, self._triaxial)

1096 except Exception as x:

1097 raise IntersectionError(pov=pov, los=los, cause=x)

1098 return Vector4Tuple(v.x, v.y, v.z, d, name=self.hartzell4.__name__)

1099

1100 @Property_RO

1101 def _hash(self):

1102 return hash((self.a, self.f))

1103

1104 def height4(self, xyz, normal=True):

1105 '''Compute the projection on and the height of a cartesian above or below

1106 this ellipsoid's surface.

1107

1108 @arg xyz: The cartesian (C{Cartesian}, L{Ecef9Tuple}, L{Vector3d},

1109 L{Vector3Tuple} or L{Vector4Tuple}).

1110 @kwarg normal: If C{True}, the projection is perpendicular to (the nearest

1111 point on) this ellipsoid's surface, otherwise the C{radial}

1112 line to this ellipsoid's center (C{bool}).

1113

1114 @return: L{Vector4Tuple}C{(x, y, z, h)} with the cartesian coordinates C{x},

1115 C{y} and C{z} of the projection on and with the height C{h} above

1116 the ellipsoid's surface, all in C{meter}, conventionally.

1117

1118 @raise ValueError: Null B{C{xyz}}.

1119

1120 @raise TypeError: Non-cartesian B{C{xyz}}.

1121

1122 @see: U{Distance to<https://StackOverflow.com/questions/22959698/distance-from-given-point-to-given-ellipse>}

1123 and U{intersection with<https://MathWorld.wolfram.com/Ellipse-LineIntersection.html>} an ellipse and

1124 methods L{Ellipsoid.hartzell4} and L{Triaxial.height4}.

1125 '''

1126 v = _MODS.vector3d._otherV3d(xyz=xyz)

1127 r = v.length

1128

1129 a, b, i = self.a, self.b, None

1130 if r < EPS0: # EPS

1131 v = v.times(_0_0)

1132 h = -a

1133

1134 elif self.isSpherical:

1135 v = v.times(a / r)

1136 h = r - a

1137

1138 elif normal: # perpendicular to ellipsoid

1139 x, y = hypot(v.x, v.y), fabs(v.z)

1140 if x < EPS0: # PYCHOK no cover

1141 z = copysign0(b, v.z)

1142 v = Vector3Tuple(v.x, v.y, z)

1143 h = y - b # polar

1144 elif y < EPS0: # PYCHOK no cover

1145 t = a / r

1146 v = v.times_(t, t, 0) # force z=0.0

1147 h = x - a # equatorial

1148 else: # normal in 1st quadrant

1149 x, y, i = _normalTo3(x, y, self)

1150 t, v = v, v.times_(x, x, y)

1151 h = t.minus(v).length

1152

1153 else: # radial to ellipsoid's center

1154 h = hypot_(a * v.z, b * v.x, b * v.y)

1155 t = (a * b / h) if h > EPS0 else _0_0 # EPS

1156 v = v.times(t)

1157 h = r * (_1_0 - t)

1158

1159 return Vector4Tuple(v.x, v.y, v.z, h, iteration=i,

1160 name=self.height4.__name__)

1161

1162 def _hubeny_2(self, phi2, phi1, lam21, scaled=True, squared=True):

1163 '''(INTERNAL) like function C{pygeodesy.flatLocal_}/C{pygeodesy.hubeny_},

1164 returning the I{angular} distance in C{radians squared} or C{radians}

1165 '''

1166 m, n = self.roc2_((phi2 + phi1) * _0_5, scaled=scaled)

1167 h, r = (hypot2, self.a2_) if squared else (hypot, _1_0 / self.a)

1168 return h(m * (phi2 - phi1), n * lam21) * r

1169

1170 @Property_RO

1171 def isEllipsoidal(self):

1172 '''Is this model I{ellipsoidal} (C{bool})?

1173 '''

1174 return self.f != 0

1175

1176 @Property_RO

1177 def isOblate(self):

1178 '''Is this ellipsoid I{oblate} (C{bool})? I{Prolate} or

1179 spherical otherwise.

1180 '''

1181 return self.f > 0

1182

1183 @Property_RO

1184 def isProlate(self):

1185 '''Is this ellipsoid I{prolate} (C{bool})? I{Oblate} or

1186 spherical otherwise.

1187 '''

1188 return self.f < 0

1189

1190 @Property_RO

1191 def isSpherical(self):

1192 '''Is this ellipsoid I{spherical} (C{bool})?

1193 '''

1194 return self.f == 0

1195

1196 def _Kseries(self, *AB8Ks):

1197 '''(INTERNAL) Compute the 4-, 6- or 8-th order I{Krüger} Alpha

1198 or Beta series coefficients per I{Karney}'s U{equations (35)

1199 and (36)<https://ArXiv.org/pdf/1002.1417v3.pdf>}.

1200

1201 @arg AB8Ks: 8-Tuple of 8-th order I{Krüger} Alpha or Beta series

1202 coefficient tuples.

1203

1204 @return: I{Krüger} series coefficients (L{KsOrder}C{-tuple}).

1205

1206 @see: I{Karney}'s 30-th order U{TMseries30

1207 <https://GeographicLib.SourceForge.io/C++/doc/tmseries30.html>}.

1208 '''

1209 k = self.KsOrder

1210 if self.n:

1211 ns = fpowers(self.n, k)

1212 ks = tuple(fdot(AB8Ks[i][:k-i], *ns[i:]) for i in range(k))

1213 else:

1214 ks = _0_0s(k)

1215 return ks

1216

1217 @property_doc_(''' the I{Krüger} series' order (C{int}), see properties C{AlphaKs}, C{BetaKs}.''')

1218 def KsOrder(self):

1219 '''Get the I{Krüger} series' order (C{int} 4, 6 or 8).

1220 '''

1221 return self._KsOrder

1222

1223 @KsOrder.setter # PYCHOK setter!

1224 def KsOrder(self, order):

1225 '''Set the I{Krüger} series' order (C{int} 4, 6 or 8).

1226

1227 @raise ValueError: Invalid B{C{order}}.

1228 '''

1229 if not (isint(order) and order in (4, 6, 8)):

1230 raise _ValueError(order=order)

1231 if order != self._KsOrder:

1232 Ellipsoid.AlphaKs._update(self)

1233 Ellipsoid.BetaKs._update(self)

1234 self._KsOrder = order

1235

1236 @Property_RO

1237 def L(self):

1238 '''Get the I{quarter meridian} C{L}, aka the C{polar distance}

1239 along a meridian between the equator and a pole (C{meter}),

1240 M{b * Elliptic(-e2 / (1 - e2)).E} or M{b * PI / 2}.

1241 '''

1242 r = self._elliptic_e22.cE if self.f else PI_2

1243 return Distance(L=self.b * r)

1244

1245 @Property_RO

1246 def _L_90(self):

1247 '''Get the I{quarter meridian} per degree (C{meter}), M{self.L / 90}.

1248 '''

1249 return self.L / _90_0

1250

1251 def Llat(self, lat):

1252 '''Return the I{meridional length}, the distance along a meridian

1253 between the equator and a (geodetic) latitude, see C{L}.

1254

1255 @arg lat: Geodetic latitude (C{degrees90}).

1256

1257 @return: The meridional length at B{C{lat}}, negative on southern

1258 hemisphere (C{meter}).

1259 '''

1260 r = self._elliptic_e22.fEd(self.auxParametric(lat)) if self.f else Phi_(lat)

1261 return Distance(Llat=self.b * r)

1262

1263 Lmeridian = Llat # meridional distance

1264

1265 @deprecated_Property_RO

1266 def majoradius(self): # PYCHOK no cover

1267 '''DEPRECATED, use property C{a} or C{Requatorial}.'''

1268 return self.a

1269

1270 def m2degrees(self, distance, lat=0):

1271 '''Convert a distance to an angle along the equator or

1272 along a parallel of (geodetic) latitude.

1273

1274 @arg distance: Distance (C{meter}).

1275 @kwarg lat: Parallel latitude (C{degrees90}, C{str}).

1276

1277 @return: Angle (C{degrees}) or C{INF} for near-polar B{C{lat}}.

1278

1279 @raise RangeError: Latitude B{C{lat}} outside valid range and

1280 L{pygeodesy.rangerrors} set to C{True}.

1281

1282 @raise ValueError: Invalid B{C{distance}} or B{C{lat}}.

1283 '''

1284 return degrees(self.m2radians(distance, lat=lat))

1285

1286 def m2radians(self, distance, lat=0):

1287 '''Convert a distance to an angle along the equator or

1288 along a parallel of (geodetic) latitude.

1289

1290 @arg distance: Distance (C{meter}).

1291 @kwarg lat: Parallel latitude (C{degrees90}, C{str}).

1292

1293 @return: Angle (C{radians}) or C{INF} for near-polar B{C{lat}}.

1294

1295 @raise RangeError: Latitude B{C{lat}} outside valid range and

1296 L{pygeodesy.rangerrors} set to C{True}.

1297

1298 @raise ValueError: Invalid B{C{distance}} or B{C{lat}}.

1299 '''

1300 r = self.circle4(lat).radius if lat else self.a

1301 return m2radians(distance, radius=r, lat=0)

1302

1303 @deprecated_Property_RO

1304 def minoradius(self): # PYCHOK no cover

1305 '''DEPRECATED, use property C{b}, C{polaradius} or C{Rpolar}.'''

1306 return self.b

1307

1308 @Property_RO

1309 def n(self):

1310 '''Get the I{3rd flattening} (C{float}), M{f / (2 - f) == (a - b) / (a + b)}, see C{a_b2n}.

1311 '''

1312 return self._assert(a_b2n(self.a, self.b), n=f2n(self.f))

1313

1314 flattening = f

1315 flattening1st = f

1316 flattening2nd = f2

1317 flattening3rd = n

1318

1319 polaradius = b # Rpolar

1320

1321# @Property_RO

1322# def Q(self):

1323# '''Get the I{meridian arc unit} C{Q}, the mean, meridional length I{per radian} C({float}).

1324#

1325# @note: C{Q * PI / 2} ≈ C{L}, the I{quarter meridian}.

1326#

1327# @see: Property C{A} and U{Engsager, K., Poder, K.<https://StudyLib.net/doc/7443565/

1328# a-highly-accurate-world-wide-algorithm-for-the-transverse...>}.

1329# '''

1330# n = self.n

1331# d = (n + _1_0) / self.a

1332# return Float(Q=Fhorner(n**2, _1_0, _0_25, _1_16th, _0_25).fover(d) if d else self.b)

1333

1334# # Moritz, H. <https://Geodesy.Geology.Ohio-State.EDU/course/refpapers/00740128.pdf>

1335# # Q = (1 - 3/4 * e'2 + 45/64 * e'4 - 175/256 * e'6 + 11025/16384 * e'8) * rocPolar

1336# # = (4 + e'2 * (-3 + e'2 * (45/16 + e'2 * (-175/64 + e'2 * 11025/4096)))) * rocPolar / 4

1337# return Fhorner(self.e22, 4, -3, 45 / 16, -175 / 64, 11025 / 4096).fover(4 / self.rocPolar)

1338

1339 @deprecated_Property_RO

1340 def quarteradius(self): # PYCHOK no cover

1341 '''DEPRECATED, use property C{L} or method C{Llat}.'''

1342 return self.L

1343

1344 @Property_RO

1345 def R1(self):

1346 '''Get the I{mean} earth radius per I{IUGG} (C{meter}), M{(2 * a + b) / 3 == a * (1 - f / 3)}.

1347

1348 @see: U{Earth radius<https://WikiPedia.org/wiki/Earth_radius>}

1349 and method C{Rgeometric}.

1350 '''

1351 r = Fsum(self.a, self.a, self.b).fover(_3_0) if self.f else self.a

1352 return Radius(R1=r)

1353

1354 Rmean = R1

1355

1356 @Property_RO

1357 def R2(self):

1358 '''Get the I{authalic} earth radius (C{meter}), M{sqrt(c2)}.

1359

1360 @see: C{R2x}, C{c2}, C{area} and U{Earth radius

1361 <https://WikiPedia.org/wiki/Earth_radius>}.

1362 '''

1363 return Radius(R2=sqrt(self.c2) if self.f else self.a)

1364

1365 Rauthalic = R2

1366

1367# @Property_RO

1368# def R2(self):

1369# # Moritz, H. <https://Geodesy.Geology.Ohio-State.EDU/course/refpapers/00740128.pdf>

1370# # R2 = (1 - 2/3 * e'2 + 26/45 * e'4 - 100/189 * e'6 + 7034/14175 * e'8) * rocPolar

1371# # = (3 + e'2 * (-2 + e'2 * (26/15 + e'2 * (-100/63 + e'2 * 7034/4725)))) * rocPolar / 3

1372# return Fhorner(self.e22, 3, -2, 26 / 15, -100 / 63, 7034 / 4725).fover(3 / self.rocPolar)

1373

1374 @Property_RO

1375 def R2x(self):

1376 '''Get the I{authalic} earth radius (C{meter}), M{sqrt(c2x)}.

1377

1378 @see: C{R2}, C{c2x} and C{areax}.

1379 '''

1380 return Radius(R2x=sqrt(self.c2x) if self.f else self.a)

1381

1382 Rauthalicx = R2x

1383

1384 @Property_RO

1385 def R3(self):

1386 '''Get the I{volumetric} earth radius (C{meter}), M{(a * a * b)**(1/3)}.

1387

1388 @see: U{Earth radius<https://WikiPedia.org/wiki/Earth_radius>} and C{volume}.

1389 '''

1390 r = (cbrt(self.b_a) * self.a) if self.f else self.a

1391 return Radius(R3=r)

1392

1393 Rvolumetric = R3

1394

1395 def radians2m(self, rad, lat=0):

1396 '''Convert an angle to the distance along the equator or

1397 along a parallel of (geodetic) latitude.

1398

1399 @arg rad: The angle (C{radians}).

1400 @kwarg lat: Parallel latitude (C{degrees90}, C{str}).

1401

1402 @return: Distance (C{meter}, same units as the equatorial

1403 and polar radii) or C{0} for near-polar B{C{lat}}.

1404

1405 @raise RangeError: Latitude B{C{lat}} outside valid range and

1406 L{pygeodesy.rangerrors} set to C{True}.

1407

1408 @raise ValueError: Invalid B{C{rad}} or B{C{lat}}.

1409 '''

1410 r = self.circle4(lat).radius if lat else self.a

1411 return radians2m(rad, radius=r, lat=0)

1412

1413 @Property_RO

1414 def Rbiaxial(self):

1415 '''Get the I{biaxial, quadratic} mean earth radius (C{meter}), M{sqrt((a**2 + b**2) / 2)}.

1416

1417 @see: C{Rtriaxial}

1418 '''

1419 a, b = self.a, self.b

1420 q = (sqrt(_0_5 + self.b2_a2 * _0_5) * a) if a > b else (

1421 (sqrt(_0_5 + self.a2_b2 * _0_5) * b) if a < b else a)

1422 return Radius(Rbiaxial=q)

1423

1424 Requatorial = a # for consistent naming

1425

1426 def Rgeocentric(self, lat):

1427 '''Compute the I{geocentric} earth radius of (geodetic) latitude.

1428

1429 @arg lat: Latitude (C{degrees90}).

1430

1431 @return: Geocentric earth radius (C{meter}).

1432

1433 @raise ValueError: Invalid B{C{lat}}.

1434

1435 @see: U{Geocentric Radius

1436 <https://WikiPedia.org/wiki/Earth_radius#Geocentric_radius>}

1437 '''

1438 r, a = self.a, Phi_(lat)

1439 if a and self.f:

1440 if fabs(a) < PI_2:

1441 s2, c2 = _s2_c2(a)

1442 b2_a2_s2 = self.b2_a2 * s2

1443 # R == sqrt((a2**2 * c2 + b2**2 * s2) / (a2 * c2 + b2 * s2))

1444 # == sqrt(a2**2 * (c2 + (b2 / a2)**2 * s2) / (a2 * (c2 + b2 / a2 * s2)))

1445 # == sqrt(a2 * (c2 + (b2 / a2)**2 * s2) / (c2 + (b2 / a2) * s2))

1446 # == a * sqrt((c2 + b2_a2 * b2_a2 * s2) / (c2 + b2_a2 * s2))

1447 # == a * sqrt((c2 + b2_a2 * b2_a2_s2) / (c2 + b2_a2_s2))

1448 r *= sqrt((c2 + b2_a2_s2 * self.b2_a2) / (c2 + b2_a2_s2))

1449 else:

1450 r = self.b

1451 return Radius(Rgeocentric=r)

1452

1453 @Property_RO

1454 def Rgeometric(self):

1455 '''Get the I{geometric} mean earth radius (C{meter}), M{sqrt(a * b)}.

1456

1457 @see: C{R1}.

1458 '''

1459 g = sqrt(self.a * self.b) if self.f else self.a

1460 return Radius(Rgeometric=g)

1461

1462 def rhumb_(self, exact=True):

1463 '''Get the an I{exact} C{Rhumb...} instance for this ellipsoid.

1464

1465 @kwarg exact: If C{bool} or C{None} return L{Rhumb}C{(exact=B{exact},

1466 ...)}, otherwise a C{rhumb} for this ellipsoid (C{.rhumbaux},

1467 C{.rhumbsolve}, C{.rhumbx}) or an other L{Rhumb}, L{RhumbAux}

1468 or L{RhumbSolve} instance for I{this} ellipsoid.

1469

1470 @return: The C{exact} rhumb (C{Rhumb...}).

1471

1472 @raise TypeError: Invalid B{C{exact}}.

1473

1474 @raise ValueError: Incompatible B{C{exact}} ellipsoid.

1475 '''

1476 if isbool(exact): # use Rhumb for backward compatibility

1477 r = _MODS.rhumbx.Rhumb(self, exact=exact, name=self.name)

1478 else:

1479 _xinstanceof(_MODS.rhumbaux.RhumbAux,

1480 _MODS.rhumbsolve.RhumbSolve,

1481 _MODS.rhumbx.Rhumb, exact=r)

1482 if r.ellipsoid != self:

1483 raise _ValueError(exact=r, ellipsosid=r.ellipsoid)

1484 return r

1485

1486 @Property_RO

1487 def rhumbaux(self):

1488 '''Get this ellipsoid's I{Auxiliary} C{rhumbaux.RhumbAux}.

1489 '''

1490 # if not self.isEllipsoidal:

1491 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

1492 return _MODS.rhumbaux.RhumbAux(self, name=self.name)

1493

1494 @property

1495 def rhumbsolve(self):

1496 '''Get this ellipsoid's L{RhumbSolve}, the I{wrapper} around utility

1497 U{RhumbSolve<https://GeographicLib.SourceForge.io/C++/doc/GeodSolve.1.html>},

1498 provided the path to the C{RhumbSolve} executable is specified with env

1499 variable C{PYGEODESY_RHUMBSOLVE} or re-/set with this property.

1500 '''

1501 # if not self.isEllipsoidal:

1502 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

1503 return _MODS.rhumbsolve.RhumbSolve(self, path=self._rhumbsolve, name=self.name)

1504

1505 @rhumbsolve.setter # PYCHOK setter!

1506 def rhumbsolve(self, path):

1507 '''Re-/set the (fully qualified) path to the U{RhumbSolve

1508 <https://GeographicLib.SourceForge.io/C++/doc/GeodSolve.1.html>} executable,

1509 overriding env variable C{PYGEODESY_RHUMBSOLVE} (C{str}).

1510 '''

1511 self._rhumbsolve = path

1512

1513 @Property_RO

1514 def rhumbx(self):

1515 '''Get this ellipsoid's I{Elliptic} C{rhumbx.Rhumb}.

1516 '''

1517 # if not self.isEllipsoidal:

1518 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

1519 return _MODS.rhumbx.Rhumb(self, name=self.name)

1520

1521 def Rlat(self, lat):

1522 '''I{Approximate} the earth radius of (geodetic) latitude.

1523

1524 @arg lat: Latitude (C{degrees90}).

1525

1526 @return: Approximate earth radius (C{meter}).

1527

1528 @raise RangeError: Latitude B{C{lat}} outside valid range and

1529 L{pygeodesy.rangerrors} set to C{True}.

1530

1531 @raise TypeError: Invalid B{C{lat}}.

1532

1533 @raise ValueError: Invalid B{C{lat}}.

1534

1535 @note: C{Rlat(B{90})} equals C{Rpolar}.

1536

1537 @see: Method C{circle4}.

1538 '''

1539 # r = a - (a - b) * |lat| / 90

1540 r = self.a

1541 if self.f and lat: # .isEllipsoidal

1542 r -= (r - self.b) * fabs(Lat(lat)) / _90_0

1543 r = Radius(Rlat=r)

1544 return r

1545

1546 Rpolar = b # for consistent naming

1547

1548 def roc1_(self, sa, ca=None):

1549 '''Compute the I{prime-vertical}, I{normal} radius of curvature

1550 of (geodetic) latitude, I{unscaled}.

1551

1552 @arg sa: Sine of the latitude (C{float}, [-1.0..+1.0]).

1553 @kwarg ca: Optional cosine of the latitude (C{float}, [-1.0..+1.0])

1554 to use an alternate formula.

1555

1556 @return: The prime-vertical radius of curvature (C{float}).

1557

1558 @note: The delta between both formulae with C{Ellipsoids.WGS84}

1559 is less than 2 nanometer over the entire latitude range.

1560

1561 @see: Method L{roc2_} and class L{EcefYou}.

1562 '''

1563 if not self.f: # .isSpherical

1564 n = self.a

1565 elif ca is None:

1566 r = self.e2s2(sa) # see .roc2_ and _EcefBase._forward

1567 n = sqrt(self.a2 / r) if r > EPS02 else _0_0

1568 elif ca: # derived from EcefYou.forward

1569 h = hypot(ca, self.b_a * sa) if sa else fabs(ca)

1570 n = self.a / h

1571 elif sa:

1572 n = self.a2_b / fabs(sa)

1573 else:

1574 n = self.a

1575 return n

1576

1577 def roc2(self, lat, scaled=False):

1578 '''Compute the I{meridional} and I{prime-vertical}, I{normal}

1579 radii of curvature of (geodetic) latitude.

1580

1581 @arg lat: Latitude (C{degrees90}).

1582 @kwarg scaled: Scale prime_vertical by C{cos(radians(B{lat}))} (C{bool}).

1583

1584 @return: An L{Curvature2Tuple}C{(meridional, prime_vertical)} with

1585 the radii of curvature.

1586

1587 @raise ValueError: Invalid B{C{lat}}.

1588

1589 @see: Methods L{roc2_} and L{roc1_}, U{Local, flat earth approximation

1590 <https://www.EdWilliams.org/avform.htm#flat>} and meridional and

1591 prime vertical U{Radii of Curvature<https://WikiPedia.org/wiki/

1592 Earth_radius#Radii_of_curvature>}.

1593 '''

1594 return self.roc2_(Phi_(lat), scaled=scaled)

1595

1596 def roc2_(self, phi, scaled=False):

1597 '''Compute the I{meridional} and I{prime-vertical}, I{normal} radii of

1598 curvature of (geodetic) latitude.

1599

1600 @arg phi: Latitude (C{radians}).

1601 @kwarg scaled: Scale prime_vertical by C{cos(B{phi})} (C{bool}).

1602

1603 @return: An L{Curvature2Tuple}C{(meridional, prime_vertical)} with the

1604 radii of curvature.

1605

1606 @raise ValueError: Invalid B{C{phi}}.

1607

1608 @see: Methods L{roc2} and L{roc1_}, property L{rocEquatorial2}, U{Local,

1609 flat earth approximation<https://www.EdWilliams.org/avform.htm#flat>}

1610 and the meridional and prime vertical U{Radii of Curvature

1611 <https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>}.

1612 '''

1613 a = fabs(Phi(phi))

1614 if self.f:

1615 r = self.e2s2(sin(a))

1616 if r > EPS02:

1617 n = self.a / sqrt(r)

1618 m = n * self.e21 / r # PYCHOK attr

1619 else:

1620 m = n = _0_0 # PYCHOK attr

1621 else:

1622 m = n = self.a

1623 if scaled and a:

1624 n *= cos(a) if a < PI_2 else _0_0

1625 return Curvature2Tuple(Radius(rocMeridional=m),

1626 Radius(rocPrimeVertical=n))

1627

1628 def rocBearing(self, lat, bearing):

1629 '''Compute the I{directional} radius of curvature of (geodetic)

1630 latitude and compass direction.

1631

1632 @arg lat: Latitude (C{degrees90}).

1633 @arg bearing: Direction (compass C{degrees360}).

1634

1635 @return: Directional radius of curvature (C{meter}).

1636

1637 @raise RangeError: Latitude B{C{lat}} outside valid range and

1638 L{pygeodesy.rangerrors} set to C{True}.

1639

1640 @raise ValueError: Invalid B{C{lat}} or B{C{bearing}}.

1641

1642 @see: U{Radii of Curvature<https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>}

1643 '''

1644 if self.f:

1645 s2, c2 = _s2_c2(Bearing_(bearing))

1646 m, n = self.roc2_(Phi_(lat))

1647 if n < m: # == n / (c2 * n / m + s2)

1648 c2 *= n / m

1649 elif m < n: # == m / (c2 + s2 * m / n)

1650 s2 *= m / n

1651 n = m

1652 b = n / (c2 + s2) # == 1 / (c2 / m + s2 / n)

1653 else:

1654 b = self.b # == self.a

1655 return Radius(rocBearing=b)

1656

1657 @Property_RO

1658 def rocEquatorial2(self):

1659 '''Get the I{meridional} and I{prime-vertical}, I{normal} radii of curvature

1660 at the equator as L{Curvature2Tuple}C{(meridional, prime_vertical)}.

1661

1662 @see: Methods L{rocMeridional} and L{rocPrimeVertical}, properties L{b2_a},

1663 L{a2_b}, C{rocPolar} and polar and equatorial U{Radii of Curvature

1664 <https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>}.

1665 '''

1666 return Curvature2Tuple(Radius(rocMeridional0=self.b2_a if self.f else self.a),

1667 Radius(rocPrimeVertical0=self.a))

1668

1669 def rocGauss(self, lat):

1670 '''Compute the I{Gaussian} radius of curvature of (geodetic) latitude.

1671

1672 @arg lat: Latitude (C{degrees90}).

1673

1674 @return: Gaussian radius of curvature (C{meter}).

1675

1676 @raise ValueError: Invalid B{C{lat}}.

1677

1678 @see: Non-directional U{Radii of Curvature<https://WikiPedia.org/wiki/

1679 Earth_radius#Radii_of_curvature>}

1680 '''

1681 # using ...

1682 # m, n = self.roc2_(Phi_(lat))

1683 # return sqrt(m * n)

1684 # ... requires 1 or 2 sqrt

1685 g = self.b

1686 if self.f:

1687 s2, c2 = _s2_c2(Phi_(lat))

1688 g = g / (c2 + self.b2_a2 * s2)

1689 return Radius(rocGauss=g)

1690

1691 def rocMean(self, lat):

1692 '''Compute the I{mean} radius of curvature of (geodetic) latitude.

1693

1694 @arg lat: Latitude (C{degrees90}).

1695

1696 @return: Mean radius of curvature (C{meter}).

1697

1698 @raise ValueError: Invalid B{C{lat}}.

1699

1700 @see: Non-directional U{Radii of Curvature<https://WikiPedia.org/wiki/

1701 Earth_radius#Radii_of_curvature>}

1702 '''

1703 if self.f:

1704 m, n = self.roc2_(Phi_(lat))

1705 m *= n * _2_0 / (m + n) # == 2 / (1 / m + 1 / n)

1706 else:

1707 m = self.a

1708 return Radius(rocMean=m)

1709

1710 def rocMeridional(self, lat):

1711 '''Compute the I{meridional} radius of curvature of (geodetic) latitude.

1712

1713 @arg lat: Latitude (C{degrees90}).

1714

1715 @return: Meridional radius of curvature (C{meter}).

1716

1717 @raise ValueError: Invalid B{C{lat}}.

1718

1719 @see: Methods L{roc2} and L{roc2_}, U{Local, flat earth approximation

1720 <https://www.EdWilliams.org/avform.htm#flat>} and U{Radii of

1721 Curvature<https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>}.

1722 '''

1723 return self.roc2_(Phi_(lat)).meridional if lat else \

1724 self.rocEquatorial2.meridional

1725

1726 rocPolar = a2_b # synonymous

1727

1728 def rocPrimeVertical(self, lat):

1729 '''Compute the I{prime-vertical}, I{normal} radius of curvature of

1730 (geodetic) latitude, aka the I{transverse} radius of curvature.

1731

1732 @arg lat: Latitude (C{degrees90}).

1733

1734 @return: Prime-vertical radius of curvature (C{meter}).

1735

1736 @raise ValueError: Invalid B{C{lat}}.

1737

1738 @see: Methods L{roc2}, L{roc2_} and L{roc1_}, U{Local, flat earth

1739 approximation<https://www.EdWilliams.org/avform.htm#flat>} and

1740 U{Radii of Curvature<https://WikiPedia.org/wiki/

1741 Earth_radius#Radii_of_curvature>}.

1742 '''

1743 return self.roc2_(Phi_(lat)).prime_vertical if lat else \

1744 self.rocEquatorial2.prime_vertical

1745

1746 rocTransverse = rocPrimeVertical # synonymous

1747

1748 @deprecated_Property_RO

1749 def Rquadratic(self): # PYCHOK no cover

1750 '''DEPRECATED, use property C{Rbiaxial} or C{Rtriaxial}.'''

1751 return self.Rbiaxial

1752

1753 @deprecated_Property_RO

1754 def Rr(self): # PYCHOK no cover

1755 '''DEPRECATED, use property C{Rrectifying}.'''

1756 return self.Rrectifying

1757

1758 @Property_RO

1759 def Rrectifying(self):

1760 '''Get the I{rectifying} earth radius (C{meter}), M{((a**(3/2) + b**(3/2)) / 2)**(2/3)}.

1761

1762 @see: U{Earth radius<https://WikiPedia.org/wiki/Earth_radius>}.

1763 '''

1764 r = (cbrt2((_1_0 + sqrt3(self.b_a)) * _0_5) * self.a) if self.f else self.a

1765 return Radius(Rrectifying=r)

1766

1767 @deprecated_Property_RO

1768 def Rs(self): # PYCHOK no cover

1769 '''DEPRECATED, use property C{Rgeometric}.'''

1770 return self.Rgeometric

1771

1772 @Property_RO

1773 def Rtriaxial(self):

1774 '''Get the I{triaxial, quadratic} mean earth radius (C{meter}), M{sqrt((3 * a**2 + b**2) / 4)}.

1775

1776 @see: C{Rbiaxial}

1777 '''

1778 a, b = self.a, self.b

1779 q = (sqrt((_3_0 + self.b2_a2) * _0_25) * a) if a > b else (

1780 (sqrt((_3_0 * self.a2_b2 + _1_0) * _0_25) * b) if a < b else a)

1781 return Radius(Rtriaxial=q)

1782

1783 def toEllipsoid2(self, name=NN):

1784 '''Get a copy of this ellipsoid as an L{Ellipsoid2}.

1785

1786 @kwarg name: Optional, unique name (C{str}).

1787

1788 @see: Property C{a_f}.

1789 '''

1790 return Ellipsoid2(self, None, name=name)

1791

1792 def toStr(self, prec=8, name=NN, **unused): # PYCHOK expected

1793 '''Return this ellipsoid as a text string.

1794

1795 @kwarg prec: Number of decimal digits, unstripped (C{int}).

1796 @kwarg name: Override name (C{str}) or C{None} to exclude

1797 this ellipsoid's name.

1798

1799 @return: This C{Ellipsoid}'s attributes (C{str}).

1800 '''

1801 E = Ellipsoid

1802 t = E.a, E.b, E.f_, E.f, E.f2, E.n, E.e, E.e2, E.e21, E.e22, E.e32, \

1803 E.A, E.L, E.R1, E.R2, E.R3, E.Rbiaxial, E.Rtriaxial

1804 return self._instr(name, prec, props=t)

1805

1806 def toTriaxial(self, name=NN):

1807 '''Convert this ellipsoid to a L{Triaxial_}.

1808

1809 @return: A L{Triaxial_} or L{Triaxial} with the C{X} axis

1810 pointing east and C{Z} pointing north.

1811

1812 @see: Method L{Triaxial_.toEllipsoid}.

1813 '''

1814 T = self._triaxial

1815 return T.copy(name=name) if name else T

1816

1817 @Property_RO

1818 def _triaxial(self):

1819 '''(INTERNAL) Get this ellipsoid's un-/ordered C{Triaxial/_}.

1820 '''

1821 a, b, m = self.a, self.b, _MODS.triaxials

1822 T = m.Triaxial if a > b else m.Triaxial_

1823 return T(a, a, b, name=self.name)

1824

1825 @Property_RO

1826 def volume(self):

1827 '''Get the ellipsoid's I{volume} (C{meter**3}), M{4 / 3 * PI * R3**3}.

1828

1829 @see: C{R3}.

1830 '''

1831 return Meter3(volume=self.a2 * self.b * PI_3 * _4_0)

1832

1833

1834 class Ellipsoid2(Ellipsoid):

1835 '''An L{Ellipsoid} specified by I{equatorial} radius and I{flattening}.

1836 '''

1837 def __init__(self, a, f, name=NN):

1838 '''New L{Ellipsoid2}.

1839

1840 @arg a: Equatorial radius, semi-axis (C{meter}).

1841 @arg f: Flattening: (C{float} < 1.0, negative for I{prolate}).

1842 @kwarg name: Optional, unique name (C{str}).

1843

1844 @raise NameError: Ellipsoid with that B{C{name}} already exists.

1845

1846 @raise ValueError: Invalid B{C{a}} or B{C{f}}.

1847

1848 @note: C{abs(B{f}) < EPS} is forced to C{B{f}=0}, I{spherical}.

1849 Negative C{B{f}} produces a I{prolate} ellipsoid.

1850 '''

1851 if f is None and isinstance(a, Ellipsoid):

1852 Ellipsoid.__init__(self, a.a, f =a.f,

1853 b=a.b, f_=a.f_, name=name)

1854 else:

1855 Ellipsoid.__init__(self, a, f=f, name=name)

1856

1857

1858 def _spherical_a_b(a, b):

1859 '''(INTERNAL) C{True} for spherical or invalid C{a} or C{b}.

1860 '''

1861 return a < EPS0 or b < EPS0 or fabs(a - b) < EPS0

1862

1863

1864 def _spherical_f(f):

1865 '''(INTERNAL) C{True} for spherical or invalid C{f}.

1866 '''

1867 return fabs(f) < EPS or f > EPS1

1868

1869

1870 def _spherical_f_(f_):

1871 '''(INTERNAL) C{True} for spherical or invalid C{f_}.

1872 '''

1873 return fabs(f_) < EPS or fabs(f_) > _1_EPS

1874

1875

1876 def a_b2e(a, b):

1877 '''Return C{e}, the I{1st eccentricity} for a given I{equatorial} and I{polar} radius.

1878

1879 @arg a: Equatorial radius (C{scalar} > 0).

1880 @arg b: Polar radius (C{scalar} > 0).

1881

1882 @return: The I{unsigned}, (1st) eccentricity (C{float} or C{0}),

1883 M{sqrt(1 - (b / a)**2)}.

1884

1885 @note: The result is always I{non-negative} and C{0} for I{near-spherical} ellipsoids.

1886 '''

1887 return Float(e=sqrt(fabs(a_b2e2(a, b)))) # == sqrt(fabs(a - b) * (a + b)) / a)

1888

1889

1890 def a_b2e2(a, b):

1891 '''Return C{e2}, the I{1st eccentricity squared} for a given I{equatorial} and I{polar} radius.

1892

1893 @arg a: Equatorial radius (C{scalar} > 0).

1894 @arg b: Polar radius (C{scalar} > 0).

1895

1896 @return: The I{signed}, (1st) eccentricity I{squared} (C{float} or C{0}),

1897 M{1 - (b / a)**2}.

1898

1899 @note: The result is positive for I{oblate}, negative for I{prolate}

1900 or C{0} for I{near-spherical} ellipsoids.

1901 '''

1902 return Float(e2=_0_0 if _spherical_a_b(a, b) else ((a - b) * (a + b) / a**2))

1903

1904

1905 def a_b2e22(a, b):

1906 '''Return C{e22}, the I{2nd eccentricity squared} for a given I{equatorial} and I{polar} radius.

1907

1908 @arg a: Equatorial radius (C{scalar} > 0).

1909 @arg b: Polar radius (C{scalar} > 0).

1910

1911 @return: The I{signed}, 2nd eccentricity I{squared} (C{float} or C{0}),

1912 M{(a / b)**2 - 1}.

1913

1914 @note: The result is positive for I{oblate}, negative for I{prolate}

1915 or C{0} for I{near-spherical} ellipsoids.

1916 '''

1917 return Float(e22=_0_0 if _spherical_a_b(a, b) else ((a - b) * (a + b) / b**2))

1918

1919

1920 def a_b2e32(a, b):

1921 '''Return C{e32}, the I{3rd eccentricity squared} for a given I{equatorial} and I{polar} radius.

1922

1923 @arg a: Equatorial radius (C{scalar} > 0).

1924 @arg b: Polar radius (C{scalar} > 0).

1925

1926 @return: The I{signed}, 3rd eccentricity I{squared} (C{float} or C{0}),

1927 M{(a**2 - b**2) / (a**2 + b**2)}.

1928

1929 @note: The result is positive for I{oblate}, negative for I{prolate}

1930 or C{0} for I{near-spherical} ellipsoids.

1931 '''

1932 a2, b2 = a**2, b**2

1933 return Float(e32=_0_0 if _spherical_a_b(a2, b2) else ((a2 - b2) / (a2 + b2)))

1934

1935

1936 def a_b2f(a, b):

1937 '''Return C{f}, the I{flattening} for a given I{equatorial} and I{polar} radius.

1938

1939 @arg a: Equatorial radius (C{scalar} > 0).

1940 @arg b: Polar radius (C{scalar} > 0).

1941

1942 @return: The flattening (C{float} or C{0}), M{(a - b) / a}.

1943

1944 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0}

1945 for I{near-spherical} ellipsoids.

1946 '''

1947 f = 0 if _spherical_a_b(a, b) else ((a - b) / a)

1948 return _f_0_0 if _spherical_f(f) else Float(f=f)

1949

1950

1951 def a_b2f_(a, b):

1952 '''Return C{f_}, the I{inverse flattening} for a given I{equatorial} and I{polar} radius.

1953

1954 @arg a: Equatorial radius (C{scalar} > 0).

1955 @arg b: Polar radius (C{scalar} > 0).

1956

1957 @return: The inverse flattening (C{float} or C{0}), M{a / (a - b)}.

1958

1959 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0}

1960 for I{near-spherical} ellipsoids.

1961 '''

1962 f_ = 0 if _spherical_a_b(a, b) else (a / float(a - b))

1963 return _f__0_0 if _spherical_f_(f_) else Float(f_=f_)

1964

1965

1966 def a_b2f2(a, b):

1967 '''Return C{f2}, the I{2nd flattening} for a given I{equatorial} and I{polar} radius.

1968

1969 @arg a: Equatorial radius (C{scalar} > 0).

1970 @arg b: Polar radius (C{scalar} > 0).

1971

1972 @return: The I{signed}, 2nd flattening (C{float} or C{0}), M{(a - b) / b}.

1973

1974 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0}

1975 for I{near-spherical} ellipsoids.

1976 '''

1977 t = 0 if _spherical_a_b(a, b) else float(a - b)

1978 return Float(f2=_0_0 if fabs(t) < EPS0 else (t / b))

1979

1980

1981 def a_b2n(a, b):

1982 '''Return C{n}, the I{3rd flattening} for a given I{equatorial} and I{polar} radius.

1983

1984 @arg a: Equatorial radius (C{scalar} > 0).

1985 @arg b: Polar radius (C{scalar} > 0).

1986

1987 @return: The I{signed}, 3rd flattening (C{float} or C{0}), M{(a - b) / (a + b)}.

1988

1989 @note: The result is positive for I{oblate}, negative for I{prolate}

1990 or C{0} for I{near-spherical} ellipsoids.

1991 '''

1992 t = 0 if _spherical_a_b(a, b) else float(a - b)

1993 return Float(n=_0_0 if fabs(t) < EPS0 else (t / (a + b)))

1994

1995

1996 def a_f2b(a, f):

1997 '''Return C{b}, the I{polar} radius for a given I{equatorial} radius and I{flattening}.

1998

1999 @arg a: Equatorial radius (C{scalar} > 0).

2000 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2001

2002 @return: The polar radius (C{float}), M{a * (1 - f)}.

2003 '''

2004 b = a if _spherical_f(f) else (a * (_1_0 - f))

2005 return Radius_(b=a if _spherical_a_b(a, b) else b)

2006

2007

2008 def a_f_2b(a, f_):

2009 '''Return C{b}, the I{polar} radius for a given I{equatorial} radius and I{inverse flattening}.

2010

2011 @arg a: Equatorial radius (C{scalar} > 0).

2012 @arg f_: Inverse flattening (C{scalar} >>> 1).

2013

2014 @return: The polar radius (C{float}), M{a * (f_ - 1) / f_}.

2015 '''

2016 b = a if _spherical_f_(f_) else (a * (f_ - _1_0) / f_)

2017 return Radius_(b=a if _spherical_a_b(a, b) else b)

2018

2019

2020 def b_f2a(b, f):

2021 '''Return C{a}, the I{equatorial} radius for a given I{polar} radius and I{flattening}.

2022

2023 @arg b: Polar radius (C{scalar} > 0).

2024 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2025

2026 @return: The equatorial radius (C{float}), M{b / (1 - f)}.

2027 '''

2028 t = _1_0 - f

2029 a = b if fabs(t) < EPS0 else (b / t)

2030 return Radius_(a=b if _spherical_a_b(a, b) else a)

2031

2032

2033 def b_f_2a(b, f_):

2034 '''Return C{a}, the I{equatorial} radius for a given I{polar} radius and I{inverse flattening}.

2035

2036 @arg b: Polar radius (C{scalar} > 0).

2037 @arg f_: Inverse flattening (C{scalar} >>> 1).

2038

2039 @return: The equatorial radius (C{float}), M{b * f_ / (f_ - 1)}.

2040 '''

2041 t = f_ - _1_0

2042 a = b if _spherical_f_(f_) or fabs(t - f_) < EPS0 \

2043 or fabs(t) < EPS0 else (b * f_ / t)

2044 return Radius_(a=b if _spherical_a_b(a, b) else a)

2045

2046

2047 def e2f(e):

2048 '''Return C{f}, the I{flattening} for a given I{1st eccentricity}.

2049

2050 @arg e: The (1st) eccentricity (0 <= C{float} < 1)

2051

2052 @return: The flattening (C{float} or C{0}).

2053

2054 @see: Function L{e22f}.

2055 '''

2056 return e22f(e**2)

2057

2058

2059 def e22f(e2):

2060 '''Return C{f}, the I{flattening} for a given I{1st eccentricity squared}.

2061

2062 @arg e2: The (1st) eccentricity I{squared}, I{signed} (L{NINF} < C{float} < 1)

2063

2064 @return: The flattening (C{float} or C{0}), M{e2 / (sqrt(e2 - 1) + 1)}.

2065 '''

2066 return Float(f=e2 / (sqrt(_1_0 - e2) + _1_0)) if e2 else _f_0_0

2067

2068

2069 def f2e2(f):

2070 '''Return C{e2}, the I{1st eccentricity squared} for a given I{flattening}.

2071

2072 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2073

2074 @return: The I{signed}, (1st) eccentricity I{squared} (C{float} < 1),

2075 M{f * (2 - f)}.

2076

2077 @note: The result is positive for I{oblate}, negative for I{prolate}

2078 or C{0} for I{near-spherical} ellipsoids.

2079

2080 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2081 html/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2082 <https://WikiPedia.org/wiki/Flattening>}.

2083 '''

2084 return Float(e2=_0_0 if _spherical_f(f) else (f * (_2_0 - f)))

2085

2086

2087 def f2e22(f):

2088 '''Return C{e22}, the I{2nd eccentricity squared} for a given I{flattening}.

2089

2090 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2091

2092 @return: The I{signed}, 2nd eccentricity I{squared} (C{float} > -1 or

2093 C{INF}), M{f * (2 - f) / (1 - f)**2}.

2094

2095 @note: The result is positive for I{oblate}, negative for I{prolate}

2096 or C{0} for near-spherical ellipsoids.

2097

2098 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2099 html/classGeographicLib_1_1Ellipsoid.html>}.

2100 '''

2101 # e2 / (1 - e2) == f * (2 - f) / (1 - f)**2

2102 t = (_1_0 - f)**2

2103 return Float(e22=INF if t < EPS0 else (f2e2(f) / t)) # PYCHOK type

2104

2105

2106 def f2e32(f):

2107 '''Return C{e32}, the I{3rd eccentricity squared} for a given I{flattening}.

2108

2109 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2110

2111 @return: The I{signed}, 3rd eccentricity I{squared} (C{float}),

2112 M{f * (2 - f) / (1 + (1 - f)**2)}.

2113

2114 @note: The result is positive for I{oblate}, negative for I{prolate}

2115 or C{0} for I{near-spherical} ellipsoids.

2116

2117 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2118 html/classGeographicLib_1_1Ellipsoid.html>}.

2119 '''

2120 # e2 / (2 - e2) == f * (2 - f) / (1 + (1 - f)**2)

2121 e2 = f2e2(f)

2122 return Float(e32=e2 / (_2_0 - e2))

2123

2124

2125 def f_2f(f_):

2126 '''Return C{f}, the I{flattening} for a given I{inverse flattening}.

2127

2128 @arg f_: Inverse flattening (C{scalar} >>> 1).

2129

2130 @return: The flattening (C{float} or C{0}), M{1 / f_}.

2131

2132 @note: The result is positive for I{oblate}, negative for I{prolate}

2133 or C{0} for I{near-spherical} ellipsoids.

2134 '''

2135 f = 0 if _spherical_f_(f_) else _1_0 / f_

2136 return _f_0_0 if _spherical_f(f) else Float(f=f) # PYCHOK type

2137

2138

2139 def f2f_(f):

2140 '''Return C{f_}, the I{inverse flattening} for a given I{flattening}.

2141

2142 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2143

2144 @return: The inverse flattening (C{float} or C{0}), M{1 / f}.

2145

2146 @note: The result is positive for I{oblate}, negative for I{prolate}

2147 or C{0} for I{near-spherical} ellipsoids.

2148 '''

2149 f_ = 0 if _spherical_f(f) else _1_0 / f

2150 return _f__0_0 if _spherical_f_(f_) else Float(f_=f_) # PYCHOK type

2151

2152

2153 def f2f2(f):

2154 '''Return C{f2}, the I{2nd flattening} for a given I{flattening}.

2155

2156 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2157

2158 @return: The I{signed}, 2nd flattening (C{float} or C{INF}), M{f / (1 - f)}.

2159

2160 @note: The result is positive for I{oblate}, negative for I{prolate}

2161 or C{0} for I{near-spherical} ellipsoids.

2162

2163 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2164 html/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2165 <https://WikiPedia.org/wiki/Flattening>}.

2166 '''

2167 t = _1_0 - f

2168 return Float(f2=_0_0 if _spherical_f(f) else (INF if fabs(t) < EPS

2169 else (f / t))) # PYCHOK type

2170

2171

2172 def f2n(f):

2173 '''Return C{n}, the I{3rd flattening} for a given I{flattening}.

2174

2175 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2176

2177 @return: The I{signed}, 3rd flattening (-1 <= C{float} < 1),

2178 M{f / (2 - f)}.

2179

2180 @note: The result is positive for I{oblate}, negative for I{prolate}

2181 or C{0} for I{near-spherical} ellipsoids.

2182

2183 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2184 html/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2185 <https://WikiPedia.org/wiki/Flattening>}.

2186 '''

2187 return Float(n=_0_0 if _spherical_f(f) else (f / float(_2_0 - f)))

2188

2189

2190 def n2e2(n):

2191 '''Return C{e2}, the I{1st eccentricity squared} for a given I{3rd flattening}.

2192

2193 @arg n: The 3rd flattening (-1 <= C{scalar} < 1).

2194

2195 @return: The I{signed}, (1st) eccentricity I{squared} (C{float} or NINF),

2196 M{4 * n / (1 + n)**2}.

2197

2198 @note: The result is positive for I{oblate}, negative for I{prolate}

2199 or C{0} for I{near-spherical} ellipsoids.

2200

2201 @see: U{Flattening<https://WikiPedia.org/wiki/Flattening>}.

2202 '''

2203 t = (n + _1_0)**2

2204 return Float(e2=_0_0 if fabs(n) < EPS0 else

2205 (NINF if t < EPS0 else (_4_0 * n / t)))

2206

2207

2208 def n2f(n):

2209 '''Return C{f}, the I{flattening} for a given I{3rd flattening}.

2210

2211 @arg n: The 3rd flattening (-1 <= C{scalar} < 1).

2212

2213 @return: The flattening (C{float} or NINF), M{2 * n / (1 + n)}.

2214

2215 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2216 html/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2217 <https://WikiPedia.org/wiki/Flattening>}.

2218 '''

2219 t = n + _1_0

2220 f = 0 if fabs(n) < EPS0 else (NINF if t < EPS0 else (_2_0 * n / t))

2221 return _f_0_0 if _spherical_f(f) else Float(f=f)

2222

2223

2224 def n2f_(n):

2225 '''Return C{f_}, the I{inverse flattening} for a given I{3rd flattening}.

2226

2227 @arg n: The 3rd flattening (-1 <= C{scalar} < 1).

2228

2229 @return: The inverse flattening (C{float} or C{0}), M{1 / f}.

2230

2231 @see: L{n2f} and L{f2f_}.

2232 '''

2233 return f2f_(n2f(n))

2234

2235

2236 def _normalTo3(px, py, E): # in .height4 above

2237 '''(INTERNAL) Nearest point on a 2-D ellipse in 1st quadrant.

2238

2239 @see: Functions C{.triaxial._normalTo4} and C{-To5}.

2240 '''

2241 a, b = E.a, E.b

2242 if min(px, py, a, b) < EPS0:

2243 raise _AssertionError(px=px, py=py, a=a, b=b, E=E)

2244

2245 a2 = a - b * E.b_a

2246 b2 = b - a * E.a_b

2247 tx = ty = _SQRT2_2

2248 for i in range(16): # max 5

2249 ex = a2 * tx**3

2250 ey = b2 * ty**3

2251

2252 qx = px - ex

2253 qy = py - ey

2254 q = hypot(qx, qy)

2255 if q < EPS0: # PYCHOK no cover

2256 break

2257 r = hypot(ex - tx * a, ey - ty * b) / q

2258

2259 sx, tx = tx, min(_1_0, max(0, (ex + qx * r) / a))

2260 sy, ty = ty, min(_1_0, max(0, (ey + qy * r) / b))

2261 t = hypot(ty, tx)

2262 if t < EPS0: # PYCHOK no cover

2263 break

2264 tx = tx / t # /= chokes PyChecker

2265 ty = ty / t

2266 if fabs(sx - tx) < EPS and \

2267 fabs(sy - ty) < EPS:

2268 break

2269

2270 tx *= a / px

2271 ty *= b / py

2272 return tx, ty, i # x and y as fractions

2273

2274

2275 class Ellipsoids(_NamedEnum):

2276 '''(INTERNAL) L{Ellipsoid} registry, I{must} be a sub-class

2277 to accommodate the L{_LazyNamedEnumItem} properties.

2278 '''

2279 def _Lazy(self, a, b, f_, **kwds):

2280 '''(INTERNAL) Instantiate the L{Ellipsoid}.

2281 '''

2282 return Ellipsoid(a, b=b, f_=f_, **kwds)

2283

2284 Ellipsoids = Ellipsoids(Ellipsoid) # PYCHOK singleton

2285'''Some pre-defined L{Ellipsoid}s, all I{lazily} instantiated.'''

2286# <https://www.GNU.org/software/gama/manual/html_node/Supported-ellipsoids.html>

2287# <https://GSSC.ESA.int/navipedia/index.php/Reference_Frames_in_GNSS>

2288# <https://kb.OSU.edu/dspace/handle/1811/77986>

2289# <https://www.IBM.com/docs/en/db2/11.5?topic=systems-supported-spheroids>

2290# <https://w3.Energistics.org/archive/Epicentre/Epicentre_v3.0/DataModel/LogicalDictionary/StandardValues/ellipsoid.html>

2291# <https://GitHub.com/locationtech/proj4j/blob/master/src/main/java/org/locationtech/proj4j/datum/Ellipsoid.java>

2292 Ellipsoids._assert( # <https://WikiPedia.org/wiki/Earth_ellipsoid>

2293 Airy1830 = _lazy(_Airy1830_, *_T(6377563.396, _0_0, 299.3249646)), # b=6356256.909

2294 AiryModified = _lazy(_AiryModified_, *_T(6377340.189, _0_0, 299.3249646)), # b=6356034.448

2295# APL4_9 = _lazy('APL4_9', *_T(6378137.0, _0_0, 298.24985392)), # Appl. Phys. Lab. 1965

2296# ANS = _lazy('ANS', *_T(6378160.0, _0_0, 298.25)), # Australian Nat. Spheroid

2297# AN_SA96 = _lazy('AN_SA96', *_T(6378160.0, _0_0, 298.24985392)), # Australian Nat. South America

2298 Australia1966 = _lazy('Australia1966', *_T(6378160.0, _0_0, 298.25)), # b=6356774.7192

2299 ATS1977 = _lazy('ATS1977', *_T(6378135.0, _0_0, 298.257)), # "Average Terrestrial System"

2300 Bessel1841 = _lazy(_Bessel1841_, *_T(6377397.155, 6356078.962818, 299.152812797)),

2301 BesselModified = _lazy('BesselModified', *_T(6377492.018, _0_0, 299.1528128)),

2302# BesselNamibia = _lazy('BesselNamibia', *_T(6377483.865, _0_0, 299.1528128)),

2303 CGCS2000 = _lazy('CGCS2000', *_T(R_MA, _0_0, 298.257222101)), # BeiDou Coord System (BDC)

2304# Clarke1858 = _lazy('Clarke1858', *_T(6378293.639, _0_0, 294.260676369)),

2305 Clarke1866 = _lazy(_Clarke1866_, *_T(6378206.4, 6356583.8, 294.978698214)),

2306 Clarke1880 = _lazy('Clarke1880', *_T(6378249.145, 6356514.86954978, 293.465)),

2307 Clarke1880IGN = _lazy(_Clarke1880IGN_, *_T(6378249.2, 6356515.0, 293.466021294)),

2308 Clarke1880Mod = _lazy('Clarke1880Mod', *_T(6378249.145, 6356514.96639549, 293.466307656)), # aka Clarke1880Arc

2309 CPM1799 = _lazy('CPM1799', *_T(6375738.7, 6356671.92557493, 334.39)), # Comm. des Poids et Mesures

2310 Delambre1810 = _lazy('Delambre1810', *_T(6376428.0, 6355957.92616372, 311.5)), # Belgium

2311 Engelis1985 = _lazy('Engelis1985', *_T(6378136.05, 6356751.32272154, 298.2566)),

2312# Everest1830 = _lazy('Everest1830', *_T(6377276.345, _0_0, 300.801699997)),

2313# Everest1948 = _lazy('Everest1948', *_T(6377304.063, _0_0, 300.801699997)),

2314# Everest1956 = _lazy('Everest1956', *_T(6377301.243, _0_0, 300.801699997)),

2315 Everest1969 = _lazy('Everest1969', *_T(6377295.664, 6356094.667915, 300.801699997)),

2316 Everest1975 = _lazy('Everest1975', *_T(6377299.151, 6356098.14512013, 300.8017255)),

2317 Fisher1968 = _lazy('Fisher1968', *_T(6378150.0, 6356768.33724438, 298.3)),

2318# Fisher1968Mod = _lazy('Fisher1968Mod', *_T(6378155.0, _0_0, 298.3)),

2319 GEM10C = _lazy('GEM10C', *_T(R_MA, 6356752.31424783, 298.2572236)),

2320 GPES = _lazy('GPES', *_T(6378135.0, 6356750.0, _0_0)), # "Gen. Purpose Earth Spheroid"

2321 GRS67 = _lazy('GRS67', *_T(6378160.0, _0_0, 298.247167427)), # Lucerne b=6356774.516

2322# GRS67Truncated = _lazy('GRS67Truncated', *_T(6378160.0, _0_0, 298.25)),

2323 GRS80 = _lazy(_GRS80_, *_T(R_MA, 6356752.314140347, 298.257222101)), # IUGG, ITRS, ETRS89

2324# Hayford1924 = _lazy('Hayford1924', *_T(6378388.0, 6356911.94612795, None)), # aka Intl1924 f_=297

2325 Helmert1906 = _lazy('Helmert1906', *_T(6378200.0, 6356818.16962789, 298.3)),

2326# Hough1960 = _lazy('Hough1960', *_T(6378270.0, _0_0, 297.0)),

2327 IAU76 = _lazy('IAU76', *_T(6378140.0, _0_0, 298.257)), # Int'l Astronomical Union

2328 IERS1989 = _lazy('IERS1989', *_T(6378136.0, _0_0, 298.257)), # b=6356751.302

2329 IERS1992TOPEX = _lazy('IERS1992TOPEX', *_T(6378136.3, 6356751.61659215, 298.257223563)), # IERS/TOPEX/Poseidon/McCarthy

2330 IERS2003 = _lazy('IERS2003', *_T(6378136.6, 6356751.85797165, 298.25642)),

2331 Intl1924 = _lazy(_Intl1924_, *_T(6378388.0, _0_0, 297.0)), # aka Hayford b=6356911.9462795

2332 Intl1967 = _lazy('Intl1967', *_T(6378157.5, 6356772.2, 298.24961539)), # New Int'l

2333 Krassovski1940 = _lazy(_Krassovski1940_, *_T(6378245.0, 6356863.01877305, 298.3)), # spelling

2334 Krassowsky1940 = _lazy(_Krassowsky1940_, *_T(6378245.0, 6356863.01877305, 298.3)), # spelling

2335# Kaula = _lazy('Kaula', *_T(6378163.0, _0_0, 298.24)), # Kaula 1961

2336# Lerch = _lazy('Lerch', *_T(6378139.0, _0_0, 298.257)), # Lerch 1979

2337 Maupertuis1738 = _lazy('Maupertuis1738', *_T(6397300.0, 6363806.28272251, 191.0)), # France

2338 Mercury1960 = _lazy('Mercury1960', *_T(6378166.0, 6356784.28360711, 298.3)),

2339 Mercury1968Mod = _lazy('Mercury1968Mod', *_T(6378150.0, 6356768.33724438, 298.3)),

2340# MERIT = _lazy('MERIT', *_T(6378137.0, _0_0, 298.257)), # MERIT 1983

2341# NWL10D = _lazy('NWL10D', *_T(6378135.0, _0_0, 298.26)), # Naval Weapons Lab.

2342 NWL1965 = _lazy('NWL1965', *_T(6378145.0, 6356759.76948868, 298.25)), # Naval Weapons Lab.

2343# NWL9D = _lazy('NWL9D', *_T(6378145.0, 6356759.76948868, 298.25)), # NWL1965

2344 OSU86F = _lazy('OSU86F', *_T(6378136.2, 6356751.51693008, 298.2572236)),

2345 OSU91A = _lazy('OSU91A', *_T(6378136.3, 6356751.6165948, 298.2572236)),

2346# Plessis1817 = _lazy('Plessis1817', *_T(6397523.0, 6355863.0, 153.56512242)), # XXX incorrect?

2347 Plessis1817 = _lazy('Plessis1817', *_T(6376523.0, 6355862.93325557, 308.64)), # XXX IGN France 1972

2348 PZ90 = _lazy('PZ90', *_T(6378136.0, _0_0, 298.257839303)), # GLOSNASS PZ-90 and PZ-90.11

2349# SEAsia = _lazy('SEAsia', *_T(6378155.0, _0_0, 298.3)), # SouthEast Asia

2350 SGS85 = _lazy('SGS85', *_T(6378136.0, 6356751.30156878, 298.257)), # Soviet Geodetic System

2351 SoAmerican1969 = _lazy('SoAmerican1969', *_T(6378160.0, 6356774.71919531, 298.25)), # South American

2352 Struve1860 = _lazy('Struve1860', *_T(6378298.3, 6356657.14266956, 294.73)),

2353# Walbeck = _lazy('Walbeck', *_T(6376896.0, _0_0, 302.78)),

2354# WarOffice = _lazy('WarOffice', *_T(6378300.0, _0_0, 296.0)),

2355 WGS60 = _lazy('WGS60', *_T(6378165.0, 6356783.28695944, 298.3)),

2356 WGS66 = _lazy('WGS66', *_T(6378145.0, 6356759.76948868, 298.25)),

2357 WGS72 = _lazy(_WGS72_, *_T(6378135.0, _0_0, 298.26)), # b=6356750.52

2358 WGS84 = _lazy(_WGS84_, *_T(R_MA, _0_0, _f__WGS84)), # GPS b=6356752.3142451793

2359# Prolate = _lazy('Prolate', *_T(6356752.3, R_MA, _0_0)),

2360 Sphere = _lazy(_Sphere_, *_T(R_M, R_M, _0_0)), # pseudo

2361 SphereAuthalic = _lazy('SphereAuthalic', *_T(R_FM, R_FM, _0_0)), # pseudo

2362 SpherePopular = _lazy('SpherePopular', *_T(R_MA, R_MA, _0_0)) # EPSG:3857 Spheroid

2363)

2364

2365

2366 if __name__ == '__main__':

2367

2368 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_

2369 from pygeodesy.lazily import printf

2370 from pygeodesy.named import nameof

2371

2372 for E in (Ellipsoids.WGS84, Ellipsoids.GRS80, # NAD83,

2373 Ellipsoids.Sphere, Ellipsoids.SpherePopular,

2374 Ellipsoid(Ellipsoids.WGS84.b, Ellipsoids.WGS84.a, name='_Prolate')):

2375 e = f2n(E.f) - E.n

2376 printf('# %s: %s', _DOT_('Ellipsoids', E.name), E.toStr(prec=10), nl=1)

2377 printf('# e=%s, f_=%s, f=%s, n=%s (%s)', fstr(E.e, prec=13, fmt=Fmt.e),

2378 fstr(E.f_, prec=13, fmt=Fmt.e),

2379 fstr(E.f, prec=13, fmt=Fmt.e),

2380 fstr(E.n, prec=13, fmt=Fmt.e),

2381 fstr(e, prec=9, fmt=Fmt.e))

2382 printf('# %s %s', Ellipsoid.AlphaKs.name, fstr(E.AlphaKs, prec=20))

2383 printf('# %s %s', Ellipsoid.BetaKs.name, fstr(E.BetaKs, prec=20))

2384 printf('# %s %s', nameof(Ellipsoid.KsOrder), E.KsOrder) # property

2385

2386 # __doc__ of this file, force all into registry

2387 t = [NN] + Ellipsoids.toRepr(all=True, asorted=True).split(_NL_)

2388 printf(_NLATvar_.join(i.strip(_COMMA_) for i in t))

2389

2390# % python3 -m pygeodesy.ellipsoids

2391

2392# Ellipsoids.WGS84: name='WGS84', a=6378137, b=6356752.3142451793, f_=298.257223563, f=0.0033528107, f2=0.0033640898, n=0.0016792204, e=0.0818191908, e2=0.00669438, e21=0.99330562, e22=0.0067394967, e32=0.0033584313, A=6367449.1458234144, L=10001965.7293127216, R1=6371008.7714150595, R2=6371007.1809184738, R3=6371000.7900091587, Rbiaxial=6367453.6345163295, Rtriaxial=6372797.5559594007

2393# e=8.1819190842622e-02, f_=2.98257223563e+02, f=3.3528106647475e-03, n=1.6792203863837e-03 (0.0e+00)

2394# AlphaKs 0.00083773182062446994, 0.00000076085277735725, 0.00000000119764550324, 0.00000000000242917068, 0.00000000000000571182, 0.0000000000000000148, 0.00000000000000000004, 0.0

2395# BetaKs 0.00083773216405794875, 0.0000000590587015222, 0.00000000016734826653, 0.00000000000021647981, 0.00000000000000037879, 0.00000000000000000072, 0.0, 0.0

2396# KsOrder 8

2397

2398# Ellipsoids.GRS80: name='GRS80', a=6378137, b=6356752.3141403468, f_=298.257222101, f=0.0033528107, f2=0.0033640898, n=0.0016792204, e=0.081819191, e2=0.00669438, e21=0.99330562, e22=0.0067394968, e32=0.0033584313, A=6367449.1457710434, L=10001965.7292304579, R1=6371008.7713801153, R2=6371007.1808835147, R3=6371000.7899741363, Rbiaxial=6367453.6344640013, Rtriaxial=6372797.5559332585

2399# e=8.1819191042833e-02, f_=2.98257222101e+02, f=3.3528106811837e-03, n=1.6792203946295e-03 (0.0e+00)

2400# AlphaKs 0.00083773182472890429, 0.00000076085278481561, 0.00000000119764552086, 0.00000000000242917073, 0.00000000000000571182, 0.0000000000000000148, 0.00000000000000000004, 0.0

2401# BetaKs 0.0008377321681623882, 0.00000005905870210374, 0.000000000167348269, 0.00000000000021647982, 0.00000000000000037879, 0.00000000000000000072, 0.0, 0.0

2402# KsOrder 8

2403

2404# Ellipsoids.Sphere: name='Sphere', a=6371008.7714149999, b=6371008.7714149999, f_=0, f=0, f2=0, n=0, e=0, e2=0, e21=1, e22=0, e32=0, A=6371008.7714149999, L=10007557.1761167478, R1=6371008.7714149999, R2=6371008.7714149999, R3=6371008.7714149999, Rbiaxial=6371008.7714149999, Rtriaxial=6371008.7714149999

2405# e=0.0e+00, f_=0.0e+00, f=0.0e+00, n=0.0e+00 (0.0e+00)

2406# AlphaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2407# BetaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2408# KsOrder 8

2409

2410# Ellipsoids.SpherePopular: name='SpherePopular', a=6378137, b=6378137, f_=0, f=0, f2=0, n=0, e=0, e2=0, e21=1, e22=0, e32=0, A=6378137, L=10018754.171394622, R1=6378137, R2=6378137, R3=6378137, Rbiaxial=6378137, Rtriaxial=6378137

2411# e=0.0e+00, f_=0.0e+00, f=0.0e+00, n=0.0e+00 (0.0e+00)

2412# AlphaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2413# BetaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2414# KsOrder 8

2415

2416# Ellipsoids._Prolate: name='_Prolate', a=6356752.3142451793, b=6378137, f_=-297.257223563, f=-0.0033640898, f2=-0.0033528107, n=-0.0016792204, e=0.0820944379, e2=-0.0067394967, e21=1.0067394967, e22=-0.00669438, e32=-0.0033584313, A=6367449.1458234144, L=10035500.5204500332, R1=6363880.5428301189, R2=6363878.9413582645, R3=6363872.5644020075, Rbiaxial=6367453.6345163295, Rtriaxial=6362105.2243882557

2417# e=8.2094437949696e-02, f_=-2.97257223563e+02, f=-3.3640898209765e-03, n=-1.6792203863837e-03 (0.0e+00)

2418# AlphaKs -0.00084149152514366627, 0.00000076653480614871, -0.00000000120934503389, 0.0000000000024576225, -0.00000000000000578863, 0.00000000000000001502, -0.00000000000000000004, 0.0

2419# BetaKs -0.00084149187224351817, 0.00000005842735196773, -0.0000000001680487236, 0.00000000000021706261, -0.00000000000000038002, 0.00000000000000000073, -0.0, 0.0

2420# KsOrder 8

2421

2422# **) MIT License

2423#

2424# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.

2425#

2426# Permission is hereby granted, free of charge, to any person obtaining a

2427# copy of this software and associated documentation files (the "Software"),

2428# to deal in the Software without restriction, including without limitation

2429# the rights to use, copy, modify, merge, publish, distribute, sublicense,

2430# and/or sell copies of the Software, and to permit persons to whom the

2431# Software is furnished to do so, subject to the following conditions:

2432#

2433# The above copyright notice and this permission notice shall be included

2434# in all copies or substantial portions of the Software.

2435#

2436# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS

2437# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

2438# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL

2439# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR

2440# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,

2441# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR

2442# OTHER DEALINGS IN THE SOFTWARE.