Coverage for pygeodesy/ellipsoids.py: 96%

744 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-08-20 14:00 -0400

1 

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

3 

4u'''Ellipsoidal and spherical earth models. 

5 

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

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

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

9 

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

11 

12Following 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 

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

65 

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

67from 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! 

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

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

73 hypot1, hypot2, sqrt3 

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

75from 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_ 

81from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS 

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

83 _NamedEnumItem, _NamedTuple, _Pass 

84from pygeodesy.namedTuples import Distance2Tuple, Vector3Tuple, Vector4Tuple 

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

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

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

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

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

90 

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

92 

93__all__ = _ALL_LAZY.ellipsoids 

94__version__ = '23.08.20' 

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 

102def _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 

108def _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 

123class 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 

173class 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 

188class 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 

202class 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(*self._Geodesics, exact=g) 

1025 if g.ellipsoid != self: 

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

1027 return g 

1028 

1029 @property_RO 

1030 def _Geodesics(self): 

1031 '''(INTERNAL) Get all C{Geodesic...} classes, I{once}. 

1032 ''' 

1033 Ellipsoid._Geodesics = t = (_MODS.geodesicw.Geodesic, # overwrite property_RO 

1034 _MODS.geodesicx.GeodesicExact, 

1035 _MODS.geodsolve.GeodesicSolve) 

1036 return t 

1037 

1038 @Property_RO 

1039 def geodesicw(self): 

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

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

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

1043 package is installed. 

1044 ''' 

1045 # if not self.isEllipsoidal: 

1046 # raise _IsnotError(_ellipsoidal_, ellipsoid=self) 

1047 return _MODS.geodesicw.Geodesic(self) 

1048 

1049 @Property_RO 

1050 def geodesicx(self): 

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

1052 ''' 

1053 # if not self.isEllipsoidal: 

1054 # raise _IsnotError(_ellipsoidal_, ellipsoid=self) 

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

1056 

1057 @property 

1058 def geodsolve(self): 

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

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

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

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

1063 ''' 

1064 # if not self.isEllipsoidal: 

1065 # raise _IsnotError(_ellipsoidal_, ellipsoid=self) 

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

1067 

1068 @geodsolve.setter # PYCHOK setter! 

1069 def geodsolve(self, path): 

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

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

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

1073 ''' 

1074 self._geodsolve = path 

1075 

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

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

1078 from a Point-Of-View in space. 

1079 

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

1081 or L{Vector3d}). 

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

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

1084 

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

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

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

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

1089 

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

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

1092 outside this ellipsoid or points in an opposite 

1093 direction. 

1094 

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

1096 

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

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

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

1100 ''' 

1101 try: 

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

1103 except Exception as x: 

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

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

1106 

1107 @Property_RO 

1108 def _hash(self): 

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

1110 

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

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

1113 this ellipsoid's surface. 

1114 

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

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

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

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

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

1120 

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

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

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

1124 

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

1126 

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

1128 

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

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

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

1132 ''' 

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

1134 r = v.length 

1135 

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

1137 if r < EPS0: # EPS 

1138 v = v.times(_0_0) 

1139 h = -a 

1140 

1141 elif self.isSpherical: 

1142 v = v.times(a / r) 

1143 h = r - a 

1144 

1145 elif normal: # perpendicular to ellipsoid 

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

1147 if x < EPS0: # PYCHOK no cover 

1148 z = copysign0(b, v.z) 

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

1150 h = y - b # polar 

1151 elif y < EPS0: # PYCHOK no cover 

1152 t = a / r 

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

1154 h = x - a # equatorial 

1155 else: # normal in 1st quadrant 

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

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

1158 h = t.minus(v).length 

1159 

1160 else: # radial to ellipsoid's center 

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

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

1163 v = v.times(t) 

1164 h = r * (_1_0 - t) 

1165 

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

1167 name=self.height4.__name__) 

1168 

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

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

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

1172 ''' 

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

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

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

1176 

1177 @Property_RO 

1178 def isEllipsoidal(self): 

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

1180 ''' 

1181 return self.f != 0 

1182 

1183 @Property_RO 

1184 def isOblate(self): 

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

1186 spherical otherwise. 

1187 ''' 

1188 return self.f > 0 

1189 

1190 @Property_RO 

1191 def isProlate(self): 

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

1193 spherical otherwise. 

1194 ''' 

1195 return self.f < 0 

1196 

1197 @Property_RO 

1198 def isSpherical(self): 

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

1200 ''' 

1201 return self.f == 0 

1202 

1203 def _Kseries(self, *AB8Ks): 

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

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

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

1207 

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

1209 coefficient tuples. 

1210 

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

1212 

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

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

1215 ''' 

1216 k = self.KsOrder 

1217 if self.n: 

1218 ns = fpowers(self.n, k) 

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

1220 else: 

1221 ks = _0_0s(k) 

1222 return ks 

1223 

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

1225 def KsOrder(self): 

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

1227 ''' 

1228 return self._KsOrder 

1229 

1230 @KsOrder.setter # PYCHOK setter! 

1231 def KsOrder(self, order): 

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

1233 

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

1235 ''' 

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

1237 raise _ValueError(order=order) 

1238 if order != self._KsOrder: 

1239 Ellipsoid.AlphaKs._update(self) 

1240 Ellipsoid.BetaKs._update(self) 

1241 self._KsOrder = order 

1242 

1243 @Property_RO 

1244 def L(self): 

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

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

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

1248 ''' 

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

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

1251 

1252 @Property_RO 

1253 def _L_90(self): 

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

1255 ''' 

1256 return self.L / _90_0 

1257 

1258 def Llat(self, lat): 

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

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

1261 

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

1263 

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

1265 hemisphere (C{meter}). 

1266 ''' 

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

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

1269 

1270 Lmeridian = Llat # meridional distance 

1271 

1272 @deprecated_Property_RO 

1273 def majoradius(self): # PYCHOK no cover 

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

1275 return self.a 

1276 

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

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

1279 along a parallel of (geodetic) latitude. 

1280 

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

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

1283 

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

1285 

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

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

1288 

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

1290 ''' 

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

1292 

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

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

1295 along a parallel of (geodetic) latitude. 

1296 

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

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

1299 

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

1301 

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

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

1304 

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

1306 ''' 

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

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

1309 

1310 @deprecated_Property_RO 

1311 def minoradius(self): # PYCHOK no cover 

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

1313 return self.b 

1314 

1315 @Property_RO 

1316 def n(self): 

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

1318 ''' 

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

1320 

1321 flattening = f 

1322 flattening1st = f 

1323 flattening2nd = f2 

1324 flattening3rd = n 

1325 

1326 polaradius = b # Rpolar 

1327 

1328# @Property_RO 

1329# def Q(self): 

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

1331# 

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

1333# 

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

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

1336# ''' 

1337# n = self.n 

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

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

1340 

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

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

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

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

1345 

1346 @deprecated_Property_RO 

1347 def quarteradius(self): # PYCHOK no cover 

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

1349 return self.L 

1350 

1351 @Property_RO 

1352 def R1(self): 

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

1354 

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

1356 and method C{Rgeometric}. 

1357 ''' 

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

1359 return Radius(R1=r) 

1360 

1361 Rmean = R1 

1362 

1363 @Property_RO 

1364 def R2(self): 

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

1366 

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

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

1369 ''' 

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

1371 

1372 Rauthalic = R2 

1373 

1374# @Property_RO 

1375# def R2(self): 

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

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

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

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

1380 

1381 @Property_RO 

1382 def R2x(self): 

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

1384 

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

1386 ''' 

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

1388 

1389 Rauthalicx = R2x 

1390 

1391 @Property_RO 

1392 def R3(self): 

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

1394 

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

1396 ''' 

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

1398 return Radius(R3=r) 

1399 

1400 Rvolumetric = R3 

1401 

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

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

1404 along a parallel of (geodetic) latitude. 

1405 

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

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

1408 

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

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

1411 

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

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

1414 

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

1416 ''' 

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

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

1419 

1420 @Property_RO 

1421 def Rbiaxial(self): 

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

1423 

1424 @see: C{Rtriaxial} 

1425 ''' 

1426 a, b = self.a, self.b 

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

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

1429 return Radius(Rbiaxial=q) 

1430 

1431 Requatorial = a # for consistent naming 

1432 

1433 def Rgeocentric(self, lat): 

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

1435 

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

1437 

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

1439 

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

1441 

1442 @see: U{Geocentric Radius 

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

1444 ''' 

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

1446 if a and self.f: 

1447 if fabs(a) < PI_2: 

1448 s2, c2 = _s2_c2(a) 

1449 b2_a2_s2 = self.b2_a2 * s2 

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

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

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

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

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

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

1456 else: 

1457 r = self.b 

1458 return Radius(Rgeocentric=r) 

1459 

1460 @Property_RO 

1461 def Rgeometric(self): 

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

1463 

1464 @see: C{R1}. 

1465 ''' 

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

1467 return Radius(Rgeometric=g) 

1468 

1469 def rhumb_(self, exact=True): 

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

1471 

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

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

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

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

1476 

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

1478 

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

1480 

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

1482 ''' 

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

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

1485 else: 

1486 r = exact 

1487 _xinstanceof(*self._Rhumbs, exact=r) 

1488 if r.ellipsoid != self: 

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

1490 return r 

1491 

1492 @Property_RO 

1493 def rhumbaux(self): 

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

1495 ''' 

1496 # if not self.isEllipsoidal: 

1497 # raise _IsnotError(_ellipsoidal_, ellipsoid=self) 

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

1499 

1500 @property_RO 

1501 def _Rhumbs(self): 

1502 '''(INTERNAL) Get all C{Rhumb...} classes, I{once}. 

1503 ''' 

1504 Ellipsoid._Rhumbs = t = (_MODS.rhumbaux.RhumbAux, # overwrite property_RO 

1505 _MODS.rhumbsolve.RhumbSolve, 

1506 _MODS.rhumbx.Rhumb) 

1507 return t 

1508 

1509 @property 

1510 def rhumbsolve(self): 

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

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

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

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

1515 ''' 

1516 # if not self.isEllipsoidal: 

1517 # raise _IsnotError(_ellipsoidal_, ellipsoid=self) 

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

1519 

1520 @rhumbsolve.setter # PYCHOK setter! 

1521 def rhumbsolve(self, path): 

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

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

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

1525 ''' 

1526 self._rhumbsolve = path 

1527 

1528 @Property_RO 

1529 def rhumbx(self): 

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

1531 ''' 

1532 # if not self.isEllipsoidal: 

1533 # raise _IsnotError(_ellipsoidal_, ellipsoid=self) 

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

1535 

1536 def Rlat(self, lat): 

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

1538 

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

1540 

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

1542 

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

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

1545 

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

1547 

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

1549 

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

1551 

1552 @see: Method C{circle4}. 

1553 ''' 

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

1555 r = self.a 

1556 if self.f and lat: # .isEllipsoidal 

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

1558 r = Radius(Rlat=r) 

1559 return r 

1560 

1561 Rpolar = b # for consistent naming 

1562 

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

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

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

1566 

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

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

1569 to use an alternate formula. 

1570 

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

1572 

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

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

1575 

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

1577 ''' 

1578 if not self.f: # .isSpherical 

1579 n = self.a 

1580 elif ca is None: 

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

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

1583 elif ca: # derived from EcefYou.forward 

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

1585 n = self.a / h 

1586 elif sa: 

1587 n = self.a2_b / fabs(sa) 

1588 else: 

1589 n = self.a 

1590 return n 

1591 

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

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

1594 radii of curvature of (geodetic) latitude. 

1595 

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

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

1598 

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

1600 the radii of curvature. 

1601 

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

1603 

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

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

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

1607 Earth_radius#Radii_of_curvature>}. 

1608 ''' 

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

1610 

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

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

1613 curvature of (geodetic) latitude. 

1614 

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

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

1617 

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

1619 radii of curvature. 

1620 

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

1622 

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

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

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

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

1627 ''' 

1628 a = fabs(Phi(phi)) 

1629 if self.f: 

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

1631 if r > EPS02: 

1632 n = self.a / sqrt(r) 

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

1634 else: 

1635 m = n = _0_0 # PYCHOK attr 

1636 else: 

1637 m = n = self.a 

1638 if scaled and a: 

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

1640 return Curvature2Tuple(Radius(rocMeridional=m), 

1641 Radius(rocPrimeVertical=n)) 

1642 

1643 def rocBearing(self, lat, bearing): 

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

1645 latitude and compass direction. 

1646 

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

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

1649 

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

1651 

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

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

1654 

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

1656 

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

1658 ''' 

1659 if self.f: 

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

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

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

1663 c2 *= n / m 

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

1665 s2 *= m / n 

1666 n = m 

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

1668 else: 

1669 b = self.b # == self.a 

1670 return Radius(rocBearing=b) 

1671 

1672 @Property_RO 

1673 def rocEquatorial2(self): 

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

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

1676 

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

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

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

1680 ''' 

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

1682 Radius(rocPrimeVertical0=self.a)) 

1683 

1684 def rocGauss(self, lat): 

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

1686 

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

1688 

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

1690 

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

1692 

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

1694 Earth_radius#Radii_of_curvature>} 

1695 ''' 

1696 # using ... 

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

1698 # return sqrt(m * n) 

1699 # ... requires 1 or 2 sqrt 

1700 g = self.b 

1701 if self.f: 

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

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

1704 return Radius(rocGauss=g) 

1705 

1706 def rocMean(self, lat): 

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

1708 

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

1710 

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

1712 

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

1714 

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

1716 Earth_radius#Radii_of_curvature>} 

1717 ''' 

1718 if self.f: 

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

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

1721 else: 

1722 m = self.a 

1723 return Radius(rocMean=m) 

1724 

1725 def rocMeridional(self, lat): 

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

1727 

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

1729 

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

1731 

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

1733 

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

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

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

1737 ''' 

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

1739 self.rocEquatorial2.meridional 

1740 

1741 rocPolar = a2_b # synonymous 

1742 

1743 def rocPrimeVertical(self, lat): 

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

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

1746 

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

1748 

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

1750 

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

1752 

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

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

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

1756 Earth_radius#Radii_of_curvature>}. 

1757 ''' 

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

1759 self.rocEquatorial2.prime_vertical 

1760 

1761 rocTransverse = rocPrimeVertical # synonymous 

1762 

1763 @deprecated_Property_RO 

1764 def Rquadratic(self): # PYCHOK no cover 

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

1766 return self.Rbiaxial 

1767 

1768 @deprecated_Property_RO 

1769 def Rr(self): # PYCHOK no cover 

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

1771 return self.Rrectifying 

1772 

1773 @Property_RO 

1774 def Rrectifying(self): 

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

1776 

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

1778 ''' 

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

1780 return Radius(Rrectifying=r) 

1781 

1782 @deprecated_Property_RO 

1783 def Rs(self): # PYCHOK no cover 

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

1785 return self.Rgeometric 

1786 

1787 @Property_RO 

1788 def Rtriaxial(self): 

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

1790 

1791 @see: C{Rbiaxial} 

1792 ''' 

1793 a, b = self.a, self.b 

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

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

1796 return Radius(Rtriaxial=q) 

1797 

1798 def toEllipsoid2(self, name=NN): 

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

1800 

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

1802 

1803 @see: Property C{a_f}. 

1804 ''' 

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

1806 

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

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

1809 

1810 @kwarg prec: Number of decimal digits, unstripped (C{int}). 

1811 @kwarg name: Override name (C{str}) or C{None} to exclude 

1812 this ellipsoid's name. 

1813 

1814 @return: This C{Ellipsoid}'s attributes (C{str}). 

1815 ''' 

1816 E = Ellipsoid 

1817 t = E.a, E.b, E.f_, E.f, E.f2, E.n, E.e, E.e2, E.e21, E.e22, E.e32, \ 

1818 E.A, E.L, E.R1, E.R2, E.R3, E.Rbiaxial, E.Rtriaxial 

1819 return self._instr(name, prec, props=t) 

1820 

1821 def toTriaxial(self, name=NN): 

1822 '''Convert this ellipsoid to a L{Triaxial_}. 

1823 

1824 @return: A L{Triaxial_} or L{Triaxial} with the C{X} axis 

1825 pointing east and C{Z} pointing north. 

1826 

1827 @see: Method L{Triaxial_.toEllipsoid}. 

1828 ''' 

1829 T = self._triaxial 

1830 return T.copy(name=name) if name else T 

1831 

1832 @Property_RO 

1833 def _triaxial(self): 

1834 '''(INTERNAL) Get this ellipsoid's un-/ordered C{Triaxial/_}. 

1835 ''' 

1836 a, b, m = self.a, self.b, _MODS.triaxials 

1837 T = m.Triaxial if a > b else m.Triaxial_ 

1838 return T(a, a, b, name=self.name) 

1839 

1840 @Property_RO 

1841 def volume(self): 

1842 '''Get the ellipsoid's I{volume} (C{meter**3}), M{4 / 3 * PI * R3**3}. 

1843 

1844 @see: C{R3}. 

1845 ''' 

1846 return Meter3(volume=self.a2 * self.b * PI_3 * _4_0) 

1847 

1848 

1849class Ellipsoid2(Ellipsoid): 

1850 '''An L{Ellipsoid} specified by I{equatorial} radius and I{flattening}. 

1851 ''' 

1852 def __init__(self, a, f, name=NN): 

1853 '''New L{Ellipsoid2}. 

1854 

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

1856 @arg f: Flattening: (C{float} < 1.0, negative for I{prolate}). 

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

1858 

1859 @raise NameError: Ellipsoid with that B{C{name}} already exists. 

1860 

1861 @raise ValueError: Invalid B{C{a}} or B{C{f}}. 

1862 

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

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

1865 ''' 

1866 if f is None and isinstance(a, Ellipsoid): 

1867 Ellipsoid.__init__(self, a.a, f =a.f, 

1868 b=a.b, f_=a.f_, name=name) 

1869 else: 

1870 Ellipsoid.__init__(self, a, f=f, name=name) 

1871 

1872 

1873def _spherical_a_b(a, b): 

1874 '''(INTERNAL) C{True} for spherical or invalid C{a} or C{b}. 

1875 ''' 

1876 return a < EPS0 or b < EPS0 or fabs(a - b) < EPS0 

1877 

1878 

1879def _spherical_f(f): 

1880 '''(INTERNAL) C{True} for spherical or invalid C{f}. 

1881 ''' 

1882 return fabs(f) < EPS or f > EPS1 

1883 

1884 

1885def _spherical_f_(f_): 

1886 '''(INTERNAL) C{True} for spherical or invalid C{f_}. 

1887 ''' 

1888 return fabs(f_) < EPS or fabs(f_) > _1_EPS 

1889 

1890 

1891def a_b2e(a, b): 

1892 '''Return C{e}, the I{1st eccentricity} for a given I{equatorial} and I{polar} radius. 

1893 

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

1895 @arg b: Polar radius (C{scalar} > 0). 

1896 

1897 @return: The I{unsigned}, (1st) eccentricity (C{float} or C{0}), 

1898 M{sqrt(1 - (b / a)**2)}. 

1899 

1900 @note: The result is always I{non-negative} and C{0} for I{near-spherical} ellipsoids. 

1901 ''' 

1902 return Float(e=sqrt(fabs(a_b2e2(a, b)))) # == sqrt(fabs(a - b) * (a + b)) / a) 

1903 

1904 

1905def a_b2e2(a, b): 

1906 '''Return C{e2}, the I{1st 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}, (1st) eccentricity I{squared} (C{float} or C{0}), 

1912 M{1 - (b / a)**2}. 

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(e2=_0_0 if _spherical_a_b(a, b) else ((a - b) * (a + b) / a**2)) 

1918 

1919 

1920def a_b2e22(a, b): 

1921 '''Return C{e22}, the I{2nd 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}, 2nd eccentricity I{squared} (C{float} or C{0}), 

1927 M{(a / b)**2 - 1}. 

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 return Float(e22=_0_0 if _spherical_a_b(a, b) else ((a - b) * (a + b) / b**2)) 

1933 

1934 

1935def a_b2e32(a, b): 

1936 '''Return C{e32}, the I{3rd eccentricity squared} for a given I{equatorial} and I{polar} radius. 

1937 

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

1939 @arg b: Polar radius (C{scalar} > 0). 

1940 

1941 @return: The I{signed}, 3rd eccentricity I{squared} (C{float} or C{0}), 

1942 M{(a**2 - b**2) / (a**2 + b**2)}. 

1943 

1944 @note: The result is positive for I{oblate}, negative for I{prolate} 

1945 or C{0} for I{near-spherical} ellipsoids. 

1946 ''' 

1947 a2, b2 = a**2, b**2 

1948 return Float(e32=_0_0 if _spherical_a_b(a2, b2) else ((a2 - b2) / (a2 + b2))) 

1949 

1950 

1951def a_b2f(a, b): 

1952 '''Return C{f}, the I{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 flattening (C{float} or C{0}), M{(a - b) / a}. 

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 - b) / a) 

1963 return _f_0_0 if _spherical_f(f) else Float(f=f) 

1964 

1965 

1966def a_b2f_(a, b): 

1967 '''Return C{f_}, the I{inverse 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 inverse flattening (C{float} or C{0}), M{a / (a - 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 f_ = 0 if _spherical_a_b(a, b) else (a / float(a - b)) 

1978 return _f__0_0 if _spherical_f_(f_) else Float(f_=f_) 

1979 

1980 

1981def a_b2f2(a, b): 

1982 '''Return C{f2}, the I{2nd 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}, 2nd flattening (C{float} or C{0}), M{(a - b) / b}. 

1988 

1989 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0} 

1990 for I{near-spherical} ellipsoids. 

1991 ''' 

1992 t = 0 if _spherical_a_b(a, b) else float(a - b) 

1993 return Float(f2=_0_0 if fabs(t) < EPS0 else (t / b)) 

1994 

1995 

1996def a_b2n(a, b): 

1997 '''Return C{n}, the I{3rd flattening} for a given I{equatorial} and I{polar} radius. 

1998 

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

2000 @arg b: Polar radius (C{scalar} > 0). 

2001 

2002 @return: The I{signed}, 3rd flattening (C{float} or C{0}), M{(a - b) / (a + b)}. 

2003 

2004 @note: The result is positive for I{oblate}, negative for I{prolate} 

2005 or C{0} for I{near-spherical} ellipsoids. 

2006 ''' 

2007 t = 0 if _spherical_a_b(a, b) else float(a - b) 

2008 return Float(n=_0_0 if fabs(t) < EPS0 else (t / (a + b))) 

2009 

2010 

2011def a_f2b(a, f): 

2012 '''Return C{b}, the I{polar} radius for a given I{equatorial} radius and I{flattening}. 

2013 

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

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

2016 

2017 @return: The polar radius (C{float}), M{a * (1 - f)}. 

2018 ''' 

2019 b = a if _spherical_f(f) else (a * (_1_0 - f)) 

2020 return Radius_(b=a if _spherical_a_b(a, b) else b) 

2021 

2022 

2023def a_f_2b(a, f_): 

2024 '''Return C{b}, the I{polar} radius for a given I{equatorial} radius and I{inverse flattening}. 

2025 

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

2027 @arg f_: Inverse flattening (C{scalar} >>> 1). 

2028 

2029 @return: The polar radius (C{float}), M{a * (f_ - 1) / f_}. 

2030 ''' 

2031 b = a if _spherical_f_(f_) else (a * (f_ - _1_0) / f_) 

2032 return Radius_(b=a if _spherical_a_b(a, b) else b) 

2033 

2034 

2035def b_f2a(b, f): 

2036 '''Return C{a}, the I{equatorial} radius for a given I{polar} radius and I{flattening}. 

2037 

2038 @arg b: Polar radius (C{scalar} > 0). 

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

2040 

2041 @return: The equatorial radius (C{float}), M{b / (1 - f)}. 

2042 ''' 

2043 t = _1_0 - f 

2044 a = b if fabs(t) < EPS0 else (b / t) 

2045 return Radius_(a=b if _spherical_a_b(a, b) else a) 

2046 

2047 

2048def b_f_2a(b, f_): 

2049 '''Return C{a}, the I{equatorial} radius for a given I{polar} radius and I{inverse flattening}. 

2050 

2051 @arg b: Polar radius (C{scalar} > 0). 

2052 @arg f_: Inverse flattening (C{scalar} >>> 1). 

2053 

2054 @return: The equatorial radius (C{float}), M{b * f_ / (f_ - 1)}. 

2055 ''' 

2056 t = f_ - _1_0 

2057 a = b if _spherical_f_(f_) or fabs(t - f_) < EPS0 \ 

2058 or fabs(t) < EPS0 else (b * f_ / t) 

2059 return Radius_(a=b if _spherical_a_b(a, b) else a) 

2060 

2061 

2062def e2f(e): 

2063 '''Return C{f}, the I{flattening} for a given I{1st eccentricity}. 

2064 

2065 @arg e: The (1st) eccentricity (0 <= C{float} < 1) 

2066 

2067 @return: The flattening (C{float} or C{0}). 

2068 

2069 @see: Function L{e22f}. 

2070 ''' 

2071 return e22f(e**2) 

2072 

2073 

2074def e22f(e2): 

2075 '''Return C{f}, the I{flattening} for a given I{1st eccentricity squared}. 

2076 

2077 @arg e2: The (1st) eccentricity I{squared}, I{signed} (L{NINF} < C{float} < 1) 

2078 

2079 @return: The flattening (C{float} or C{0}), M{e2 / (sqrt(e2 - 1) + 1)}. 

2080 ''' 

2081 return Float(f=e2 / (sqrt(_1_0 - e2) + _1_0)) if e2 else _f_0_0 

2082 

2083 

2084def f2e2(f): 

2085 '''Return C{e2}, the I{1st eccentricity squared} for a given I{flattening}. 

2086 

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

2088 

2089 @return: The I{signed}, (1st) eccentricity I{squared} (C{float} < 1), 

2090 M{f * (2 - f)}. 

2091 

2092 @note: The result is positive for I{oblate}, negative for I{prolate} 

2093 or C{0} for I{near-spherical} ellipsoids. 

2094 

2095 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/ 

2096 html/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening 

2097 <https://WikiPedia.org/wiki/Flattening>}. 

2098 ''' 

2099 return Float(e2=_0_0 if _spherical_f(f) else (f * (_2_0 - f))) 

2100 

2101 

2102def f2e22(f): 

2103 '''Return C{e22}, the I{2nd eccentricity squared} for a given I{flattening}. 

2104 

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

2106 

2107 @return: The I{signed}, 2nd eccentricity I{squared} (C{float} > -1 or 

2108 C{INF}), M{f * (2 - f) / (1 - f)**2}. 

2109 

2110 @note: The result is positive for I{oblate}, negative for I{prolate} 

2111 or C{0} for near-spherical ellipsoids. 

2112 

2113 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/ 

2114 html/classGeographicLib_1_1Ellipsoid.html>}. 

2115 ''' 

2116 # e2 / (1 - e2) == f * (2 - f) / (1 - f)**2 

2117 t = (_1_0 - f)**2 

2118 return Float(e22=INF if t < EPS0 else (f2e2(f) / t)) # PYCHOK type 

2119 

2120 

2121def f2e32(f): 

2122 '''Return C{e32}, the I{3rd eccentricity squared} for a given I{flattening}. 

2123 

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

2125 

2126 @return: The I{signed}, 3rd eccentricity I{squared} (C{float}), 

2127 M{f * (2 - f) / (1 + (1 - f)**2)}. 

2128 

2129 @note: The result is positive for I{oblate}, negative for I{prolate} 

2130 or C{0} for I{near-spherical} ellipsoids. 

2131 

2132 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/ 

2133 html/classGeographicLib_1_1Ellipsoid.html>}. 

2134 ''' 

2135 # e2 / (2 - e2) == f * (2 - f) / (1 + (1 - f)**2) 

2136 e2 = f2e2(f) 

2137 return Float(e32=e2 / (_2_0 - e2)) 

2138 

2139 

2140def f_2f(f_): 

2141 '''Return C{f}, the I{flattening} for a given I{inverse flattening}. 

2142 

2143 @arg f_: Inverse flattening (C{scalar} >>> 1). 

2144 

2145 @return: The flattening (C{float} or C{0}), M{1 / f_}. 

2146 

2147 @note: The result is positive for I{oblate}, negative for I{prolate} 

2148 or C{0} for I{near-spherical} ellipsoids. 

2149 ''' 

2150 f = 0 if _spherical_f_(f_) else _1_0 / f_ 

2151 return _f_0_0 if _spherical_f(f) else Float(f=f) # PYCHOK type 

2152 

2153 

2154def f2f_(f): 

2155 '''Return C{f_}, the I{inverse flattening} for a given I{flattening}. 

2156 

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

2158 

2159 @return: The inverse flattening (C{float} or C{0}), M{1 / f}. 

2160 

2161 @note: The result is positive for I{oblate}, negative for I{prolate} 

2162 or C{0} for I{near-spherical} ellipsoids. 

2163 ''' 

2164 f_ = 0 if _spherical_f(f) else _1_0 / f 

2165 return _f__0_0 if _spherical_f_(f_) else Float(f_=f_) # PYCHOK type 

2166 

2167 

2168def f2f2(f): 

2169 '''Return C{f2}, the I{2nd flattening} for a given I{flattening}. 

2170 

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

2172 

2173 @return: The I{signed}, 2nd flattening (C{float} or C{INF}), M{f / (1 - f)}. 

2174 

2175 @note: The result is positive for I{oblate}, negative for I{prolate} 

2176 or C{0} for I{near-spherical} ellipsoids. 

2177 

2178 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/ 

2179 html/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening 

2180 <https://WikiPedia.org/wiki/Flattening>}. 

2181 ''' 

2182 t = _1_0 - f 

2183 return Float(f2=_0_0 if _spherical_f(f) else (INF if fabs(t) < EPS 

2184 else (f / t))) # PYCHOK type 

2185 

2186 

2187def f2n(f): 

2188 '''Return C{n}, the I{3rd flattening} for a given I{flattening}. 

2189 

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

2191 

2192 @return: The I{signed}, 3rd flattening (-1 <= C{float} < 1), 

2193 M{f / (2 - f)}. 

2194 

2195 @note: The result is positive for I{oblate}, negative for I{prolate} 

2196 or C{0} for I{near-spherical} ellipsoids. 

2197 

2198 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/ 

2199 html/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening 

2200 <https://WikiPedia.org/wiki/Flattening>}. 

2201 ''' 

2202 return Float(n=_0_0 if _spherical_f(f) else (f / float(_2_0 - f))) 

2203 

2204 

2205def n2e2(n): 

2206 '''Return C{e2}, the I{1st eccentricity squared} for a given I{3rd flattening}. 

2207 

2208 @arg n: The 3rd flattening (-1 <= C{scalar} < 1). 

2209 

2210 @return: The I{signed}, (1st) eccentricity I{squared} (C{float} or NINF), 

2211 M{4 * n / (1 + n)**2}. 

2212 

2213 @note: The result is positive for I{oblate}, negative for I{prolate} 

2214 or C{0} for I{near-spherical} ellipsoids. 

2215 

2216 @see: U{Flattening<https://WikiPedia.org/wiki/Flattening>}. 

2217 ''' 

2218 t = (n + _1_0)**2 

2219 return Float(e2=_0_0 if fabs(n) < EPS0 else 

2220 (NINF if t < EPS0 else (_4_0 * n / t))) 

2221 

2222 

2223def n2f(n): 

2224 '''Return C{f}, the I{flattening} for a given I{3rd flattening}. 

2225 

2226 @arg n: The 3rd flattening (-1 <= C{scalar} < 1). 

2227 

2228 @return: The flattening (C{float} or NINF), M{2 * n / (1 + n)}. 

2229 

2230 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/ 

2231 html/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening 

2232 <https://WikiPedia.org/wiki/Flattening>}. 

2233 ''' 

2234 t = n + _1_0 

2235 f = 0 if fabs(n) < EPS0 else (NINF if t < EPS0 else (_2_0 * n / t)) 

2236 return _f_0_0 if _spherical_f(f) else Float(f=f) 

2237 

2238 

2239def n2f_(n): 

2240 '''Return C{f_}, the I{inverse flattening} for a given I{3rd flattening}. 

2241 

2242 @arg n: The 3rd flattening (-1 <= C{scalar} < 1). 

2243 

2244 @return: The inverse flattening (C{float} or C{0}), M{1 / f}. 

2245 

2246 @see: L{n2f} and L{f2f_}. 

2247 ''' 

2248 return f2f_(n2f(n)) 

2249 

2250 

2251def _normalTo3(px, py, E): # in .height4 above 

2252 '''(INTERNAL) Nearest point on a 2-D ellipse in 1st quadrant. 

2253 

2254 @see: Functions C{.triaxial._normalTo4} and C{-To5}. 

2255 ''' 

2256 a, b = E.a, E.b 

2257 if min(px, py, a, b) < EPS0: 

2258 raise _AssertionError(px=px, py=py, a=a, b=b, E=E) 

2259 

2260 a2 = a - b * E.b_a 

2261 b2 = b - a * E.a_b 

2262 tx = ty = _SQRT2_2 

2263 for i in range(16): # max 5 

2264 ex = a2 * tx**3 

2265 ey = b2 * ty**3 

2266 

2267 qx = px - ex 

2268 qy = py - ey 

2269 q = hypot(qx, qy) 

2270 if q < EPS0: # PYCHOK no cover 

2271 break 

2272 r = hypot(ex - tx * a, ey - ty * b) / q 

2273 

2274 sx, tx = tx, min(_1_0, max(0, (ex + qx * r) / a)) 

2275 sy, ty = ty, min(_1_0, max(0, (ey + qy * r) / b)) 

2276 t = hypot(ty, tx) 

2277 if t < EPS0: # PYCHOK no cover 

2278 break 

2279 tx = tx / t # /= chokes PyChecker 

2280 ty = ty / t 

2281 if fabs(sx - tx) < EPS and \ 

2282 fabs(sy - ty) < EPS: 

2283 break 

2284 

2285 tx *= a / px 

2286 ty *= b / py 

2287 return tx, ty, i # x and y as fractions 

2288 

2289 

2290class Ellipsoids(_NamedEnum): 

2291 '''(INTERNAL) L{Ellipsoid} registry, I{must} be a sub-class 

2292 to accommodate the L{_LazyNamedEnumItem} properties. 

2293 ''' 

2294 def _Lazy(self, a, b, f_, **kwds): 

2295 '''(INTERNAL) Instantiate the L{Ellipsoid}. 

2296 ''' 

2297 return Ellipsoid(a, b=b, f_=f_, **kwds) 

2298 

2299Ellipsoids = Ellipsoids(Ellipsoid) # PYCHOK singleton 

2300'''Some pre-defined L{Ellipsoid}s, all I{lazily} instantiated.''' 

2301# <https://www.GNU.org/software/gama/manual/html_node/Supported-ellipsoids.html> 

2302# <https://GSSC.ESA.int/navipedia/index.php/Reference_Frames_in_GNSS> 

2303# <https://kb.OSU.edu/dspace/handle/1811/77986> 

2304# <https://www.IBM.com/docs/en/db2/11.5?topic=systems-supported-spheroids> 

2305# <https://w3.Energistics.org/archive/Epicentre/Epicentre_v3.0/DataModel/LogicalDictionary/StandardValues/ellipsoid.html> 

2306# <https://GitHub.com/locationtech/proj4j/blob/master/src/main/java/org/locationtech/proj4j/datum/Ellipsoid.java> 

2307Ellipsoids._assert( # <https://WikiPedia.org/wiki/Earth_ellipsoid> 

2308 Airy1830 = _lazy(_Airy1830_, *_T(6377563.396, _0_0, 299.3249646)), # b=6356256.909 

2309 AiryModified = _lazy(_AiryModified_, *_T(6377340.189, _0_0, 299.3249646)), # b=6356034.448 

2310# APL4_9 = _lazy('APL4_9', *_T(6378137.0, _0_0, 298.24985392)), # Appl. Phys. Lab. 1965 

2311# ANS = _lazy('ANS', *_T(6378160.0, _0_0, 298.25)), # Australian Nat. Spheroid 

2312# AN_SA96 = _lazy('AN_SA96', *_T(6378160.0, _0_0, 298.24985392)), # Australian Nat. South America 

2313 Australia1966 = _lazy('Australia1966', *_T(6378160.0, _0_0, 298.25)), # b=6356774.7192 

2314 ATS1977 = _lazy('ATS1977', *_T(6378135.0, _0_0, 298.257)), # "Average Terrestrial System" 

2315 Bessel1841 = _lazy(_Bessel1841_, *_T(6377397.155, 6356078.962818, 299.152812797)), 

2316 BesselModified = _lazy('BesselModified', *_T(6377492.018, _0_0, 299.1528128)), 

2317# BesselNamibia = _lazy('BesselNamibia', *_T(6377483.865, _0_0, 299.1528128)), 

2318 CGCS2000 = _lazy('CGCS2000', *_T(R_MA, _0_0, 298.257222101)), # BeiDou Coord System (BDC) 

2319# Clarke1858 = _lazy('Clarke1858', *_T(6378293.639, _0_0, 294.260676369)), 

2320 Clarke1866 = _lazy(_Clarke1866_, *_T(6378206.4, 6356583.8, 294.978698214)), 

2321 Clarke1880 = _lazy('Clarke1880', *_T(6378249.145, 6356514.86954978, 293.465)), 

2322 Clarke1880IGN = _lazy(_Clarke1880IGN_, *_T(6378249.2, 6356515.0, 293.466021294)), 

2323 Clarke1880Mod = _lazy('Clarke1880Mod', *_T(6378249.145, 6356514.96639549, 293.466307656)), # aka Clarke1880Arc 

2324 CPM1799 = _lazy('CPM1799', *_T(6375738.7, 6356671.92557493, 334.39)), # Comm. des Poids et Mesures 

2325 Delambre1810 = _lazy('Delambre1810', *_T(6376428.0, 6355957.92616372, 311.5)), # Belgium 

2326 Engelis1985 = _lazy('Engelis1985', *_T(6378136.05, 6356751.32272154, 298.2566)), 

2327# Everest1830 = _lazy('Everest1830', *_T(6377276.345, _0_0, 300.801699997)), 

2328# Everest1948 = _lazy('Everest1948', *_T(6377304.063, _0_0, 300.801699997)), 

2329# Everest1956 = _lazy('Everest1956', *_T(6377301.243, _0_0, 300.801699997)), 

2330 Everest1969 = _lazy('Everest1969', *_T(6377295.664, 6356094.667915, 300.801699997)), 

2331 Everest1975 = _lazy('Everest1975', *_T(6377299.151, 6356098.14512013, 300.8017255)), 

2332 Fisher1968 = _lazy('Fisher1968', *_T(6378150.0, 6356768.33724438, 298.3)), 

2333# Fisher1968Mod = _lazy('Fisher1968Mod', *_T(6378155.0, _0_0, 298.3)), 

2334 GEM10C = _lazy('GEM10C', *_T(R_MA, 6356752.31424783, 298.2572236)), 

2335 GPES = _lazy('GPES', *_T(6378135.0, 6356750.0, _0_0)), # "Gen. Purpose Earth Spheroid" 

2336 GRS67 = _lazy('GRS67', *_T(6378160.0, _0_0, 298.247167427)), # Lucerne b=6356774.516 

2337# GRS67Truncated = _lazy('GRS67Truncated', *_T(6378160.0, _0_0, 298.25)), 

2338 GRS80 = _lazy(_GRS80_, *_T(R_MA, 6356752.314140347, 298.25722210088)), # IUGG, ITRS, ETRS89 

2339# Hayford1924 = _lazy('Hayford1924', *_T(6378388.0, 6356911.94612795, None)), # aka Intl1924 f_=297 

2340 Helmert1906 = _lazy('Helmert1906', *_T(6378200.0, 6356818.16962789, 298.3)), 

2341# Hough1960 = _lazy('Hough1960', *_T(6378270.0, _0_0, 297.0)), 

2342 IAU76 = _lazy('IAU76', *_T(6378140.0, _0_0, 298.257)), # Int'l Astronomical Union 

2343 IERS1989 = _lazy('IERS1989', *_T(6378136.0, _0_0, 298.257)), # b=6356751.302 

2344 IERS1992TOPEX = _lazy('IERS1992TOPEX', *_T(6378136.3, 6356751.61659215, 298.257223563)), # IERS/TOPEX/Poseidon/McCarthy 

2345 IERS2003 = _lazy('IERS2003', *_T(6378136.6, 6356751.85797165, 298.25642)), 

2346 Intl1924 = _lazy(_Intl1924_, *_T(6378388.0, _0_0, 297.0)), # aka Hayford b=6356911.9462795 

2347 Intl1967 = _lazy('Intl1967', *_T(6378157.5, 6356772.2, 298.24961539)), # New Int'l 

2348 Krassovski1940 = _lazy(_Krassovski1940_, *_T(6378245.0, 6356863.01877305, 298.3)), # spelling 

2349 Krassowsky1940 = _lazy(_Krassowsky1940_, *_T(6378245.0, 6356863.01877305, 298.3)), # spelling 

2350# Kaula = _lazy('Kaula', *_T(6378163.0, _0_0, 298.24)), # Kaula 1961 

2351# Lerch = _lazy('Lerch', *_T(6378139.0, _0_0, 298.257)), # Lerch 1979 

2352 Maupertuis1738 = _lazy('Maupertuis1738', *_T(6397300.0, 6363806.28272251, 191.0)), # France 

2353 Mercury1960 = _lazy('Mercury1960', *_T(6378166.0, 6356784.28360711, 298.3)), 

2354 Mercury1968Mod = _lazy('Mercury1968Mod', *_T(6378150.0, 6356768.33724438, 298.3)), 

2355# MERIT = _lazy('MERIT', *_T(6378137.0, _0_0, 298.257)), # MERIT 1983 

2356# NWL10D = _lazy('NWL10D', *_T(6378135.0, _0_0, 298.26)), # Naval Weapons Lab. 

2357 NWL1965 = _lazy('NWL1965', *_T(6378145.0, 6356759.76948868, 298.25)), # Naval Weapons Lab. 

2358# NWL9D = _lazy('NWL9D', *_T(6378145.0, 6356759.76948868, 298.25)), # NWL1965 

2359 OSU86F = _lazy('OSU86F', *_T(6378136.2, 6356751.51693008, 298.2572236)), 

2360 OSU91A = _lazy('OSU91A', *_T(6378136.3, 6356751.6165948, 298.2572236)), 

2361# Plessis1817 = _lazy('Plessis1817', *_T(6397523.0, 6355863.0, 153.56512242)), # XXX incorrect? 

2362 Plessis1817 = _lazy('Plessis1817', *_T(6376523.0, 6355862.93325557, 308.64)), # XXX IGN France 1972 

2363 PZ90 = _lazy('PZ90', *_T(6378136.0, _0_0, 298.257839303)), # GLOSNASS PZ-90 and PZ-90.11 

2364# SEAsia = _lazy('SEAsia', *_T(6378155.0, _0_0, 298.3)), # SouthEast Asia 

2365 SGS85 = _lazy('SGS85', *_T(6378136.0, 6356751.30156878, 298.257)), # Soviet Geodetic System 

2366 SoAmerican1969 = _lazy('SoAmerican1969', *_T(6378160.0, 6356774.71919531, 298.25)), # South American 

2367 Struve1860 = _lazy('Struve1860', *_T(6378298.3, 6356657.14266956, 294.73)), 

2368# Walbeck = _lazy('Walbeck', *_T(6376896.0, _0_0, 302.78)), 

2369# WarOffice = _lazy('WarOffice', *_T(6378300.0, _0_0, 296.0)), 

2370 WGS60 = _lazy('WGS60', *_T(6378165.0, 6356783.28695944, 298.3)), 

2371 WGS66 = _lazy('WGS66', *_T(6378145.0, 6356759.76948868, 298.25)), 

2372 WGS72 = _lazy(_WGS72_, *_T(6378135.0, _0_0, 298.26)), # b=6356750.52 

2373 WGS84 = _lazy(_WGS84_, *_T(R_MA, _0_0, _f__WGS84)), # GPS b=6356752.3142451793 

2374# Prolate = _lazy('Prolate', *_T(6356752.3, R_MA, _0_0)), 

2375 Sphere = _lazy(_Sphere_, *_T(R_M, R_M, _0_0)), # pseudo 

2376 SphereAuthalic = _lazy('SphereAuthalic', *_T(R_FM, R_FM, _0_0)), # pseudo 

2377 SpherePopular = _lazy('SpherePopular', *_T(R_MA, R_MA, _0_0)) # EPSG:3857 Spheroid 

2378) 

2379 

2380_EWGS84 = Ellipsoids.WGS84 # (INTERNAL) shared 

2381 

2382if __name__ == '__main__': 

2383 

2384 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_ 

2385 from pygeodesy import nameof, printf 

2386 

2387 for E in (_EWGS84, Ellipsoids.GRS80, # NAD83, 

2388 Ellipsoids.Sphere, Ellipsoids.SpherePopular, 

2389 Ellipsoid(_EWGS84.b, _EWGS84.a, name='_Prolate')): 

2390 e = f2n(E.f) - E.n 

2391 printf('# %s: %s', _DOT_('Ellipsoids', E.name), E.toStr(prec=10), nl=1) 

2392 printf('# e=%s, f_=%s, f=%s, n=%s (%s)', fstr(E.e, prec=13, fmt=Fmt.e), 

2393 fstr(E.f_, prec=13, fmt=Fmt.e), 

2394 fstr(E.f, prec=13, fmt=Fmt.e), 

2395 fstr(E.n, prec=13, fmt=Fmt.e), 

2396 fstr(e, prec=9, fmt=Fmt.e)) 

2397 printf('# %s %s', Ellipsoid.AlphaKs.name, fstr(E.AlphaKs, prec=20)) 

2398 printf('# %s %s', Ellipsoid.BetaKs.name, fstr(E.BetaKs, prec=20)) 

2399 printf('# %s %s', nameof(Ellipsoid.KsOrder), E.KsOrder) # property 

2400 

2401 # __doc__ of this file, force all into registry 

2402 t = [NN] + Ellipsoids.toRepr(all=True, asorted=True).split(_NL_) 

2403 printf(_NLATvar_.join(i.strip(_COMMA_) for i in t)) 

2404 

2405# % python3 -m pygeodesy.ellipsoids 

2406 

2407# 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 

2408# e=8.1819190842622e-02, f_=2.98257223563e+02, f=3.3528106647475e-03, n=1.6792203863837e-03 (0.0e+00) 

2409# AlphaKs 0.00083773182062446994, 0.00000076085277735725, 0.00000000119764550324, 0.00000000000242917068, 0.00000000000000571182, 0.0000000000000000148, 0.00000000000000000004, 0.0 

2410# BetaKs 0.00083773216405794875, 0.0000000590587015222, 0.00000000016734826653, 0.00000000000021647981, 0.00000000000000037879, 0.00000000000000000072, 0.0, 0.0 

2411# KsOrder 8 

2412 

2413# Ellipsoids.GRS80: name='GRS80', a=6378137, b=6356752.3141403468, f_=298.2572221009, 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 

2414# e=8.1819191042833e-02, f_=2.9825722210088e+02, f=3.3528106811837e-03, n=1.6792203946295e-03 (0.0e+00) 

2415# AlphaKs 0.00083773182472890429, 0.00000076085278481561, 0.00000000119764552086, 0.00000000000242917073, 0.00000000000000571182, 0.0000000000000000148, 0.00000000000000000004, 0.0 

2416# BetaKs 0.0008377321681623882, 0.00000005905870210374, 0.000000000167348269, 0.00000000000021647982, 0.00000000000000037879, 0.00000000000000000072, 0.0, 0.0 

2417# KsOrder 8 

2418 

2419# 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 

2420# e=0.0e+00, f_=0.0e+00, f=0.0e+00, n=0.0e+00 (0.0e+00) 

2421# AlphaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 

2422# BetaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 

2423# KsOrder 8 

2424 

2425# 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 

2426# e=0.0e+00, f_=0.0e+00, f=0.0e+00, n=0.0e+00 (0.0e+00) 

2427# AlphaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 

2428# BetaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 

2429# KsOrder 8 

2430 

2431# 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 

2432# e=8.2094437949696e-02, f_=-2.97257223563e+02, f=-3.3640898209765e-03, n=-1.6792203863837e-03 (0.0e+00) 

2433# AlphaKs -0.00084149152514366627, 0.00000076653480614871, -0.00000000120934503389, 0.0000000000024576225, -0.00000000000000578863, 0.00000000000000001502, -0.00000000000000000004, 0.0 

2434# BetaKs -0.00084149187224351817, 0.00000005842735196773, -0.0000000001680487236, 0.00000000000021706261, -0.00000000000000038002, 0.00000000000000000073, -0.0, 0.0 

2435# KsOrder 8 

2436 

2437# **) MIT License 

2438# 

2439# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved. 

2440# 

2441# Permission is hereby granted, free of charge, to any person obtaining a 

2442# copy of this software and associated documentation files (the "Software"), 

2443# to deal in the Software without restriction, including without limitation 

2444# the rights to use, copy, modify, merge, publish, distribute, sublicense, 

2445# and/or sell copies of the Software, and to permit persons to whom the 

2446# Software is furnished to do so, subject to the following conditions: 

2447# 

2448# The above copyright notice and this permission notice shall be included 

2449# in all copies or substantial portions of the Software. 

2450# 

2451# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 

2452# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 

2453# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 

2454# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR 

2455# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, 

2456# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR 

2457# OTHER DEALINGS IN THE SOFTWARE.