1

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

3

4 u'''Ellipsoidal and spherical earth models.

5

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

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

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

9

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

11

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

13

14@var Ellipsoids.Airy1830: Ellipsoid(name='Airy1830', a=6377563.396, b=6356256.90923729, f_=299.3249646, f=0.00334085, f2=0.00335205, n=0.00167322, e=0.08167337, e2=0.00667054, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, e22=0.0067395, e32=0.00335843, A=6367449.14582341, L=10001965.72931272, R1=6371008.77141506, R2=6371007.18091847, R3=6371000.79000916, Rbiaxial=6367453.63451633, Rtriaxial=6372797.5559594)

62'''

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

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

65

66 from pygeodesy.basics import copysign0, isint

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

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

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

70 _0_25, _3_0 # PYCHOK used!

71 from pygeodesy.errors import _AssertionError, _ValueError, _xkwds_not

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

73 hypot1, hypot2, sqrt3

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

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

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

77 _Intl1924_, _incompatible_, _invalid_, _Krassovski1940_, \

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

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

80 from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS

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

82 _NamedEnumItem, _NamedTuple, _Pass

83 from pygeodesy.namedTuples import Distance2Tuple, Vector3Tuple, Vector4Tuple

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

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

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

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

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

89

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

91

92__all__ = _ALL_LAZY.ellipsoids

93__version__ = '23.04.11'

94

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

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

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

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

99

100

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

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

103 '''

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

105

106

107 def _s2_c2(phi):

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

109 '''

110 if phi:

111 s2 = sin(phi)**2

112 if s2 > EPS:

113 c2 = _1_0 - s2

114 if c2 > EPS:

115 if c2 < EPS1:

116 return s2, c2

117 else:

118 return _1_0, _0_0 # phi == PI_2

119 return _0_0, _1_0 # phi == 0

120

121

122 class a_f2Tuple(_NamedTuple):

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

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

125

126 @see: Class L{Ellipsoid2}.

127 '''

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

129 _Units_ = (_Pass, _Pass)

130

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

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

133

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

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

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

137

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

139

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

141

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

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

144 '''

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

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

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

148 f = _f_0_0

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

150

151 @Property_RO

152 def b(self):

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

154 '''

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

156

157 def ellipsoid(self, name=NN):

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

159

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

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

162 '''

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

164

165 @Property_RO

166 def f_(self):

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

168 '''

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

170

171

172 class Circle4Tuple(_NamedTuple):

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

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

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

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

177

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

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

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

181 equals C{lat}.

182 '''

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

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

185

186

187 class Curvature2Tuple(_NamedTuple):

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

189 C{meter}, conventionally.

190 '''

191 _Names_ = (_meridional_, _prime_vertical_)

192 _Units_ = ( Meter, Meter)

193

194 @property_RO

195 def transverse(self):

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

197 '''

198 return self.prime_vertical

199

200

201 class Ellipsoid(_NamedEnumItem):

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

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

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

205 '''

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

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

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

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

210

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

212

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

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

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

216

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

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

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

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

221 spherical).

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

223

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

225

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

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

228

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

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

231 '''

232 ff_ = f, f_ # assertion below

233 try:

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

235 if not _isfinite(a):

236 raise ValueError(_SPACE_(_a_, _not_finite_))

237

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

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

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

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

242 elif f is not None:

243 f = Float(f=f)

244 b = a_f2b(a, f)

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

246 elif f_:

247 f_ = Float(f_=f_)

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

249 f = f_2f(f_)

250 else: # only a, spherical

251 f_ = f = 0

252 b = a # superfluous

253

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

255 raise ValueError(_SPACE_(_f_, _invalid_))

256 if not _isfinite(b):

257 raise ValueError(_SPACE_(_b_, _not_finite_))

258

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

260 b = a

261 f = _f_0_0

262 f_ = _f__0_0

263

264 except (TypeError, ValueError) as x:

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

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

267 raise _ValueError(t, cause=x)

268

269 self._a = a

270 self._b = b

271 self._f = f

272 self._f_ = f_

273

274 self._register(Ellipsoids, name)

275

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

277 d = dict(eps=_TOL)

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

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

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

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

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

283

284 def __eq__(self, other):

285 '''Compare this and an other ellipsoid.

286

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

288

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

290 '''

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

292 self.a == other.a and

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

294

295 @Property_RO

296 def a(self):

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

298 '''

299 return self._a

300

301 equatoradius = a # = Requatorial

302

303 @Property_RO

304 def a2(self):

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

306 '''

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

308

309 @Property_RO

310 def a2_(self):

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

312 '''

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

314

315 @Property_RO

316 def a_b(self):

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

318 '''

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

320

321 @Property_RO

322 def a2_b(self):

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

324

325 @see: U{Radii of Curvature

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

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

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

329

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

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

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

333 '''

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

335

336 @Property_RO

337 def a2_b2(self):

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

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

340 '''

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

342

343 @Property_RO

344 def a_f(self):

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

346 '''

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

348

349 @Property_RO

350 def A(self):

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

352

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

354 '''

355 A, n = self.a, self.n

356 if n:

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

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

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

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

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

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

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

364 return A

365

366 @Property_RO

367 def _albersCyl(self):

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

369 '''

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

371

372 @Property_RO

373 def AlphaKs(self):

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

375 '''

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

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

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

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

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

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

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

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

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

385 _T(1424729850961/743921418240)) # PYCHOK unaligned

386

387 @Property_RO

388 def area(self):

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

390

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

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

393 '''

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

395

396 @Property_RO

397 def areax(self):

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

399 accurate for very I{oblate} ellipsoids.

400

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

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

403 '''

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

405

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

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

408 '''

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

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

411 t = (v, _vs_, val)

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

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

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

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

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

417 eps=eps, f0=f0, **name_value))

418

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

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

421

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

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

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

425

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

427

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

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

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

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

432 '''

433 if self.f:

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

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

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

437

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

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

440

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

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

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

444

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

446

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

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

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

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

451 '''

452 if self.f:

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

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

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

456

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

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

459

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

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

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

463

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

465

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

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

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

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

470 '''

471 if self.f:

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

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

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

475

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

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

478

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

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

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

482

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

484

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

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

487 thereof is the original geodetic latitude.

488

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

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

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

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

493 '''

494 if self.f:

495 r = Phi_(lat, clip=0)

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

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

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

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

500

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

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

503

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

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

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

507

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

509

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

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

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

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

514 '''

515 if self.f:

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

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

518

519 auxReduced = auxParametric # synonymous

520

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

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

523

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

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

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

527

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

529

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

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

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

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

534 '''

535 if self.f:

536 lat = Lat(lat)

537 if 0 < fabs(lat) < _90_0:

538 if inverse:

539 e = self._elliptic_e22

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

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

542 else:

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

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

545

546 @Property_RO

547 def b(self):

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

549 '''

550 return self._b

551

552 polaradius = b # = Rpolar

553

554 @Property_RO

555 def b_a(self):

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

557

558 @see: Property L{f1}.

559 '''

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

561

562 @Property_RO

563 def b2(self):

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

565 '''

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

567

568 @Property_RO

569 def b2_a(self):

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

571

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

573 '''

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

575

576 @Property_RO

577 def b2_a2(self):

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

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

580 '''

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

582

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

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

585 '''

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

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

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

589

590 @Property_RO

591 def BetaKs(self):

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

593 '''

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

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

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

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

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

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

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

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

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

603 _T(191773887257/3719607091200)) # PYCHOK unaligned

604

605 @deprecated_Property_RO

606 def c(self): # PYCHOK no cover

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

608 return self.R2

609

610 @Property_RO

611 def c2(self):

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

613

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

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

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

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

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

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

620 '''

621 return self._c2f(False)

622

623 @Property_RO

624 def c2x(self):

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

626 ellipsoids.

627

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

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

630 '''

631 return self._c2f(True)

632

633 def _c2f(self, c2x):

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

635 '''

636 f, c2 = self.f, self.b2

637 if f:

638 e = self.e

639 if e > EPS0:

640 if f > 0: # .isOblate

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

642 elif f < 0: # .isProlate

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

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

645 return c2

646

647 def circle4(self, lat):

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

649

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

651

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

653 instance.

654

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

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

657

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

659

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

661

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

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

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

665 and method C{Rlat}.

666 '''

667 lat = Lat(lat)

668 if lat:

669 b = lat

670 if fabs(lat) < _90_0:

671 if self.f:

672 b = self._beta(lat)

673 z, r = sincos2d(b)

674 r *= self.a

675 z *= self.b

676 else: # near-polar

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

678 else: # equator

679 r = self.a

680 z = lat = b = _0_0

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

682

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

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

685 along a parallel of (geodetic) latitude.

686

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

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

689

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

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

692

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

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

695

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

697 '''

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

699

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

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

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

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

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

705

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

707 and only between points not near-polar}.

708

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

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

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

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

713

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

715 in same units as this ellipsoid's axes.

716

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

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

719

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

721 '''

722 phi0 = Phi_(lat0=lat0)

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

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

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

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

727

728 @Property_RO

729 def e(self):

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

731

732 @see: Property L{es}.

733 '''

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

735

736 @deprecated_Property_RO

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

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

739 return self.e21

740

741# @Property_RO

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

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

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

745

746 @Property_RO

747 def e2(self):

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

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

750 '''

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

752

753 @Property_RO

754 def e2abs(self):

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

756 '''

757 return fabs(self.e2)

758

759 @Property_RO

760 def e21(self):

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

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

763 '''

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

765

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

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

768

769 @Property_RO

770 def e22(self):

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

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

773 '''

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

775

776 @Property_RO

777 def e22abs(self):

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

779 '''

780 return fabs(self.e22)

781

782 @Property_RO

783 def e32(self):

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

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

786 '''

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

788

789 @Property_RO

790 def e32abs(self):

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

792 '''

793 return fabs(self.e32)

794

795 @Property_RO

796 def e4(self):

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

798 '''

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

800

801 eccentricity = e # eccentricity

802# eccentricity2 = e2 # eccentricity squared

803 eccentricity1st2 = e2 # first eccentricity squared

804 eccentricity2nd2 = e22 # second eccentricity squared

805 eccentricity3rd2 = e32 # third eccentricity squared

806

807 def ecef(self, Ecef=None):

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

809

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

811

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

813

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

815

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

817 '''

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

819

820 @Property_RO

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

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

823 '''

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

825 return _MODS.elliptic.Elliptic(-e12)

826

827 @Property_RO

828 def _elliptic_e22(self): # aka ._elliptic_ep2

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

830 '''

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

832

833 equatoradius = a # Requatorial

834

835 def e2s(self, s):

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

837

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

839

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

841

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

843 '''

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

845

846 def e2s2(self, s):

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

848

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

850

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

852

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

854 '''

855 r = _1_0

856 if self.e2:

857 try:

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

859 if r < 0:

860 raise ValueError(_negative_)

861 except (TypeError, ValueError) as x:

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

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

864 return r

865

866 @Property_RO

867 def es(self):

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

869

870 @see: Property L{e}.

871 '''

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

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

874

875 def es_atanh(self, x):

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

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

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

879

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

881

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

883 html/classGeographicLib_1_1Math.html>}.

884 '''

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

886

887 def _es_atanh(self, x):

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

889 '''

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

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

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

893 _0_0)) # .isSpherical

894

895 @Property_RO

896 def es_c(self):

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

898 '''

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

900

901 def es_tauf(self, taup):

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

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

904

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

906 SourceForge.io/html/classGeographicLib_1_1Math.html>} and

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

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

909 '''

910 t = Scalar(taup=taup)

911 if self.f: # .isEllipsoidal

912 a = fabs(t)

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

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

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

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

917 if fabs(d) < s:

918 break

919 t = Scalar(tauf=T)

920 return t

921

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

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

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

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

926 '''

927 e = self._1_e21

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

929 _tf2 = self._es_taupf2

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

931 a, h = _tf2(T)

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

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

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

935 yield T, i, d

936

937 def es_taupf(self, tau):

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

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

940

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

942 SourceForge.io/html/classGeographicLib_1_1Math.html>}.

943 '''

944 t = Scalar(tau=tau)

945 if self.f: # .isEllipsoidal

946 t, _ = self._es_taupf2(t)

947 t = Scalar(taupf=t)

948 return t

949

950 def _es_taupf2(self, tau):

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

952 '''

953 if _isfinite(tau):

954 h = hypot1(tau)

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

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

957 else:

958 a, h = tau, INF

959 return a, h

960

961 @Property_RO

962 def _exp_es_atanh_1(self):

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

964 '''

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

966

967 @Property_RO

968 def f(self):

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

970 '''

971 return self._f

972

973 @Property_RO

974 def f_(self):

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

976 '''

977 return self._f_

978

979 @Property_RO

980 def f1(self):

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

982

983 @see: Property L{b_a}.

984 '''

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

986

987 @Property_RO

988 def f2(self):

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

990 '''

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

992

993 @Property_RO

994 def geodesic(self):

995 '''Get this ellipsoid's I{wrapped} U{geodesic.Geodesic

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

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

998 package is installed.

999 '''

1000 # if not self.isEllipsoidal:

1001 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

1002 return _MODS.geodesicw.Geodesic(self)

1003

1004 @Property_RO

1005 def geodesicx(self):

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

1007 '''

1008 # if not self.isEllipsoidal:

1009 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1011

1012 @Property_RO

1013 def geodsolve(self):

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

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

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

1017 variable C{PYGEODESY_GEODSOLVE}.

1018 '''

1019 # if not self.isEllipsoidal:

1020 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

1021 return _MODS.geodsolve.GeodesicSolve(self, name=self.name)

1022

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

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

1025 this ellipsoid's surface.

1026

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

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

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

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

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

1032

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

1034 C{x}, C{y} and C{z} of the projection on or the intersection

1035 with and with the height C{h} above the ellipsoid's surface

1036 in C{meter}, conventionally.

1037

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

1039

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

1041

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

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

1044 method L{Triaxial.height4}.

1045 '''

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

1047 r = v.length

1048

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

1050 if r < EPS0: # EPS

1051 v = v.times(_0_0)

1052 h = -a

1053

1054 elif self.isSpherical:

1055 v = v.times(a / r)

1056 h = r - a

1057

1058 elif normal: # perpendicular to ellipsoid

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

1060 if x < EPS0: # PYCHOK no cover

1061 z = copysign0(b, v.z)

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

1063 h = y - b # polar

1064 elif y < EPS0: # PYCHOK no cover

1065 t = a / r

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

1067 h = x - a # equatorial

1068 else: # normal in 1st quadrant

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

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

1071 h = t.minus(v).length

1072

1073 else: # radial to ellipsoid's center

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

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

1076 v = v.times(t)

1077 h = r * (_1_0 - t)

1078

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

1080 name=self.height4.__name__)

1081

1082 def _hubeny_2(self, phi2, phi1, lam21):

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

1084 but returning the I{angular} distance in C{radians squared}.

1085 '''

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

1087 return hypot2(m * (phi2 - phi1), n * lam21) * self.a2_

1088

1089 @Property_RO

1090 def isEllipsoidal(self):

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

1092 '''

1093 return self.f != 0

1094

1095 @Property_RO

1096 def isOblate(self):

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

1098 spherical otherwise.

1099 '''

1100 return self.f > 0

1101

1102 @Property_RO

1103 def isProlate(self):

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

1105 spherical otherwise.

1106 '''

1107 return self.f < 0

1108

1109 @Property_RO

1110 def isSpherical(self):

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

1112 '''

1113 return self.f == 0

1114

1115 def _Kseries(self, *AB8Ks):

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

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

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

1119

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

1121 coefficient tuples.

1122

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

1124

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

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

1127 '''

1128 k = self.KsOrder

1129 if self.n:

1130 ns = fpowers(self.n, k)

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

1132 else:

1133 ks = _0_0s(k)

1134 return ks

1135

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

1137 def KsOrder(self):

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

1139 '''

1140 return self._KsOrder

1141

1142 @KsOrder.setter # PYCHOK setter!

1143 def KsOrder(self, order):

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

1145

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

1147 '''

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

1149 raise _ValueError(order=order)

1150 if order != self._KsOrder:

1151 Ellipsoid.AlphaKs._update(self)

1152 Ellipsoid.BetaKs._update(self)

1153 self._KsOrder = order

1154

1155 @Property_RO

1156 def L(self):

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

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

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

1160 '''

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

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

1163

1164 @Property_RO

1165 def _L_90(self):

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

1167 '''

1168 return self.L / _90_0

1169

1170 def Llat(self, lat):

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

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

1173

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

1175

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

1177 hemisphere (C{meter}).

1178 '''

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

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

1181

1182 Lmeridian = Llat # meridional distance

1183

1184 @deprecated_Property_RO

1185 def majoradius(self): # PYCHOK no cover

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

1187 return self.a

1188

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

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

1191 along a parallel of (geodetic) latitude.

1192

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

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

1195

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

1197

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

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

1200

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

1202 '''

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

1204

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

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

1207 along a parallel of (geodetic) latitude.

1208

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

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

1211

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

1213

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

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

1216

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

1218 '''

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

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

1221

1222 @deprecated_Property_RO

1223 def minoradius(self): # PYCHOK no cover

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

1225 return self.b

1226

1227 @Property_RO

1228 def n(self):

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

1230 '''

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

1232

1233 flattening = f

1234 flattening1st = f

1235 flattening2nd = f2

1236 flattening3rd = n

1237

1238 polaradius = b # Rpolar

1239

1240# @Property_RO

1241# def Q(self):

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

1243#

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

1245#

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

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

1248# '''

1249# n = self.n

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

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

1252

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

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

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

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

1257

1258 @deprecated_Property_RO

1259 def quarteradius(self): # PYCHOK no cover

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

1261 return self.L

1262

1263 @Property_RO

1264 def R1(self):

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

1266

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

1268 and method C{Rgeometric}.

1269 '''

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

1271 return Radius(R1=r)

1272

1273 Rmean = R1

1274

1275 @Property_RO

1276 def R2(self):

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

1278

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

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

1281 '''

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

1283

1284 Rauthalic = R2

1285

1286# @Property_RO

1287# def R2(self):

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

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

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

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

1292

1293 @Property_RO

1294 def R2x(self):

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

1296

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

1298 '''

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

1300

1301 Rauthalicx = R2x

1302

1303 @Property_RO

1304 def R3(self):

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

1306

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

1308 '''

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

1310 return Radius(R3=r)

1311

1312 Rvolumetric = R3

1313

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

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

1316 along a parallel of (geodetic) latitude.

1317

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

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

1320

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

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

1323

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

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

1326

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

1328 '''

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

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

1331

1332 @Property_RO

1333 def Rbiaxial(self):

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

1335

1336 @see: C{Rtriaxial}

1337 '''

1338 a, b = self.a, self.b

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

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

1341 return Radius(Rbiaxial=q)

1342

1343 Requatorial = a # for consistent naming

1344

1345 def Rgeocentric(self, lat):

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

1347

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

1349

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

1351

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

1353

1354 @see: U{Geocentric Radius

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

1356 '''

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

1358 if a and self.f:

1359 if fabs(a) < PI_2:

1360 s2, c2 = _s2_c2(a)

1361 b2_a2_s2 = self.b2_a2 * s2

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

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

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

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

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

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

1368 else:

1369 r = self.b

1370 return Radius(Rgeocentric=r)

1371

1372 @Property_RO

1373 def Rgeometric(self):

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

1375

1376 @see: C{R1}.

1377 '''

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

1379 return Radius(Rgeometric=g)

1380

1381 @Property_RO

1382 def rhumbsolve(self):

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

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

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

1386 variable C{PYGEODESY_RHUMBSOLVE}.

1387 '''

1388 # if not self.isEllipsoidal:

1389 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

1390 return _MODS.rhumbsolve.RhumbSolve(self, name=self.name)

1391

1392 @Property_RO

1393 def rhumbx(self):

1394 '''Get this ellipsoid's L{Rhumb}.

1395 '''

1396 # if not self.isEllipsoidal:

1397 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1399

1400 def Rlat(self, lat):

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

1402

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

1404

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

1406

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

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

1409

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

1411

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

1413

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

1415

1416 @see: Method C{circle4}.

1417 '''

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

1419 r = self.a

1420 if self.f and lat: # .isEllipsoidal

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

1422 r = Radius(Rlat=r)

1423 return r

1424

1425 Rpolar = b # for consistent naming

1426

1427 @deprecated_Property_RO

1428 def Rquadratic(self): # PYCHOK no cover

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

1430 return self.Rbiaxial

1431

1432 @deprecated_Property_RO

1433 def Rr(self): # PYCHOK no cover

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

1435 return self.Rrectifying

1436

1437 @Property_RO

1438 def Rrectifying(self):

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

1440

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

1442 '''

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

1444 return Radius(Rrectifying=r)

1445

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

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

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

1449

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

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

1452 to use an alternate formula.

1453

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

1455

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

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

1458

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

1460 '''

1461 if not self.f: # .isSpherical

1462 n = self.a

1463 elif ca is None:

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

1465 n = (self.a / sqrt(r)) if r > EPS02 else _0_0

1466 elif ca: # derived from EcefYou.forward

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

1468 n = self.a / h

1469 elif sa:

1470 n = self.a2_b / fabs(sa)

1471 else:

1472 n = self.a

1473 return n

1474

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

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

1477 radii of curvature of (geodetic) latitude.

1478

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

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

1481

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

1483 the radii of curvature.

1484

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

1486

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

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

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

1490 Earth_radius#Radii_of_curvature>}.

1491 '''

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

1493

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

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

1496 curvature of (geodetic) latitude.

1497

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

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

1500

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

1502 radii of curvature.

1503

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

1505

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

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

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

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

1510 '''

1511 a = fabs(Phi(phi))

1512 if self.f:

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

1514 if r > EPS02:

1515 n = self.a / sqrt(r)

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

1517 else:

1518 m = n = _0_0 # PYCHOK attr

1519 else:

1520 m = n = self.a

1521 if scaled and a:

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

1523 return Curvature2Tuple(Radius(rocMeridional=m),

1524 Radius(rocPrimeVertical=n))

1525

1526 def rocBearing(self, lat, bearing):

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

1528 latitude and compass direction.

1529

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

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

1532

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

1534

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

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

1537

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

1539

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

1541 '''

1542 if self.f:

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

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

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

1546 c2 *= n / m

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

1548 s2 *= m / n

1549 n = m

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

1551 else:

1552 b = self.b # == self.a

1553 return Radius(rocBearing=b)

1554

1555 @Property_RO

1556 def rocEquatorial2(self):

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

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

1559

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

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

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

1563 '''

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

1565 Radius(rocPrimeVertical0=self.a))

1566

1567 def rocGauss(self, lat):

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

1569

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

1571

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

1573

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

1575

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

1577 Earth_radius#Radii_of_curvature>}

1578 '''

1579 # using ...

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

1581 # return sqrt(m * n)

1582 # ... requires 1 or 2 sqrt

1583 g = self.b

1584 if self.f:

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

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

1587 return Radius(rocGauss=g)

1588

1589 def rocMean(self, lat):

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

1591

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

1593

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

1595

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

1597

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

1599 Earth_radius#Radii_of_curvature>}

1600 '''

1601 if self.f:

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

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

1604 else:

1605 m = self.a

1606 return Radius(rocMean=m)

1607

1608 def rocMeridional(self, lat):

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

1610

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

1612

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

1614

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

1616

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

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

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

1620 '''

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

1622 self.rocEquatorial2.meridional

1623

1624 rocPolar = a2_b # synonymous

1625

1626 def rocPrimeVertical(self, lat):

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

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

1629

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

1631

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

1633

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

1635

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

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

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

1639 Earth_radius#Radii_of_curvature>}.

1640 '''

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

1642 self.rocEquatorial2.prime_vertical

1643

1644 rocTransverse = rocPrimeVertical # synonymous

1645

1646 @deprecated_Property_RO

1647 def Rs(self): # PYCHOK no cover

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

1649 return self.Rgeometric

1650

1651 @Property_RO

1652 def Rtriaxial(self):

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

1654

1655 @see: C{Rbiaxial}

1656 '''

1657 a, b = self.a, self.b

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

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

1660 return Radius(Rtriaxial=q)

1661

1662 def toEllipsoid2(self, name=NN):

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

1664

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

1666

1667 @see: Property C{a_f}.

1668 '''

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

1670

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

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

1673

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

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

1676 this ellipsoid's name.

1677

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

1679 '''

1680 E = Ellipsoid

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

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

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

1684

1685 def toTriaxial(self, name=NN):

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

1687

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

1689 pointing east and C{Z} pointing north.

1690

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

1692 '''

1693 a, b, n = self.a, self.b, (name or self.name)

1694 return _MODS.triaxials.Triaxial_(a, a, b, name=n) if a < b else \

1695 _MODS.triaxials.Triaxial( a, a, b, name=n)

1696

1697 @Property_RO

1698 def volume(self):

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

1700

1701 @see: C{R3}.

1702 '''

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

1704

1705

1706 class Ellipsoid2(Ellipsoid):

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

1708 '''

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

1710 '''New L{Ellipsoid2}.

1711

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

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

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

1715

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

1717

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

1719

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

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

1722 '''

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

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

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

1726 else:

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

1728

1729

1730 def _spherical_a_b(a, b):

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

1732 '''

1733 return a < EPS or b < EPS or fabs(a - b) < EPS

1734

1735

1736 def _spherical_f(f):

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

1738 '''

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

1740

1741

1742 def _spherical_f_(f_):

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

1744 '''

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

1746

1747

1748 def a_b2e(a, b):

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

1750

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

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

1753

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

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

1756

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

1758 '''

1759 return Float(e=sqrt(fabs(a_b2e2(a, b))))

1760

1761

1762 def a_b2e2(a, b):

1763 '''Return C{e2}, the I{1st eccentricity squared} for a given I{equatorial} and I{polar} radius.

1764

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

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

1767

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

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

1770

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

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

1773 '''

1774 return Float(e2=_0_0 if _spherical_a_b(a, b) else (_1_0 - (b / a)**2))

1775

1776

1777 def a_b2e22(a, b):

1778 '''Return C{e22}, the I{2nd eccentricity squared} for a given I{equatorial} and I{polar} radius.

1779

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

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

1782

1783 @return: The I{signed}, 2nd eccentricity I{squared} (C{float} or C{0}),

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

1785

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

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

1788 '''

1789 return Float(e22=_0_0 if _spherical_a_b(a, b) else ((a / b)**2 - _1_0))

1790

1791

1792 def a_b2e32(a, b):

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

1794

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

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

1797

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

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

1800

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

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

1803 '''

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

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

1806

1807

1808 def a_b2f(a, b):

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

1810

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

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

1813

1814 @return: The flattening (C{float} or C{0}), M{(a - b) / a}.

1815

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

1817 for I{near-spherical} ellipsoids.

1818 '''

1819 f = 0 if _spherical_a_b(a, b) else ((a - b) / a)

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

1821

1822

1823 def a_b2f_(a, b):

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

1825

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

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

1828

1829 @return: The inverse flattening (C{float} or C{0}), M{a / (a - b)}.

1830

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

1832 for I{near-spherical} ellipsoids.

1833 '''

1834 f_ = 0 if _spherical_a_b(a, b) else (a / float(a - b))

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

1836

1837

1838 def a_b2f2(a, b):

1839 '''Return C{f2}, the I{2nd flattening} for a given I{equatorial} and I{polar} radius.

1840

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

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

1843

1844 @return: The I{signed}, 2nd flattening (C{float} or C{0}), M{(a - b) / b}.

1845

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

1847 for I{near-spherical} ellipsoids.

1848 '''

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

1850 return Float(f2=_0_0 if fabs(t) < EPS else (t / b))

1851

1852

1853 def a_b2n(a, b):

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

1855

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

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

1858

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

1860

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

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

1863 '''

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

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

1866

1867

1868 def a_f2b(a, f):

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

1870

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

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

1873

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

1875 '''

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

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

1878

1879

1880 def a_f_2b(a, f_):

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

1882

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

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

1885

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

1887 '''

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

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

1890

1891

1892 def b_f2a(b, f):

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

1894

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

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

1897

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

1899 '''

1900 t = _1_0 - f

1901 a = b if fabs(t < EPS) else (b / t)

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

1903

1904

1905 def b_f_2a(b, f_):

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

1907

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

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

1910

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

1912 '''

1913 t = f_ - _1_0

1914 a = b if _spherical_f_(f_) or fabs(t - f_) < EPS \

1915 or fabs(t) < EPS else (b * f_ / t)

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

1917

1918

1919 def e2f(e):

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

1921

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

1923

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

1925

1926 @see: Function L{e22f}.

1927 '''

1928 return e22f(e**2)

1929

1930

1931 def e22f(e2):

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

1933

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

1935

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

1937 '''

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

1939

1940

1941 def f2e2(f):

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

1943

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

1945

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

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

1948

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

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

1951

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

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

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

1955 '''

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

1957

1958

1959 def f2e22(f):

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

1961

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

1963

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

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

1966

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

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

1969

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

1971 html/classGeographicLib_1_1Ellipsoid.html>}.

1972 '''

1973 # e2 / (1 - e2) == f * (2 - f) / (1 - f)**2

1974 t = (_1_0 - f)**2

1975 return Float(e22=INF if t < EPS else (f2e2(f) / t)) # PYCHOK type

1976

1977

1978 def f2e32(f):

1979 '''Return C{e32}, the I{3rd eccentricity squared} for a given I{flattening}.

1980

1981 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

1982

1983 @return: The I{signed}, 3rd eccentricity I{squared} (C{float}),

1984 M{f * (2 - f) / (1 + (1 - f)**2)}.

1985

1986 @note: The result is positive for I{oblate}, negative for I{prolate}

1987 or C{0} for I{near-spherical} ellipsoids.

1988

1989 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

1990 html/classGeographicLib_1_1Ellipsoid.html>}.

1991 '''

1992 # e2 / (2 - e2) == f * (2 - f) / (1 + (1 - f)**2)

1993 e2 = f2e2(f)

1994 return Float(e32=e2 / (_2_0 - e2))

1995

1996

1997 def f_2f(f_):

1998 '''Return C{f}, the I{flattening} for a given I{inverse flattening}.

1999

2000 @arg f_: Inverse flattening (C{scalar} >>> 1).

2001

2002 @return: The flattening (C{float} or C{0}), M{1 / f_}.

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 f = 0 if _spherical_f_(f_) else _1_0 / f_

2008 return _f_0_0 if _spherical_f(f) else Float(f=f) # PYCHOK type

2009

2010

2011 def f2f_(f):

2012 '''Return C{f_}, the I{inverse flattening} for a given I{flattening}.

2013

2014 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2015

2016 @return: The inverse flattening (C{float} or C{0}), M{1 / f}.

2017

2018 @note: The result is positive for I{oblate}, negative for I{prolate}

2019 or C{0} for I{near-spherical} ellipsoids.

2020 '''

2021 f_ = 0 if _spherical_f(f) else _1_0 / f

2022 return _f__0_0 if _spherical_f_(f_) else Float(f_=f_) # PYCHOK type

2023

2024

2025 def f2f2(f):

2026 '''Return C{f2}, the I{2nd flattening} for a given I{flattening}.

2027

2028 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2029

2030 @return: The I{signed}, 2nd flattening (C{float} or C{INF}), M{f / (1 - f)}.

2031

2032 @note: The result is positive for I{oblate}, negative for I{prolate}

2033 or C{0} for I{near-spherical} ellipsoids.

2034

2035 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2036 html/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2037 <https://WikiPedia.org/wiki/Flattening>}.

2038 '''

2039 t = _1_0 - f

2040 return Float(f2=_0_0 if _spherical_f(f) else (INF if fabs(t) < EPS

2041 else (f / t))) # PYCHOK type

2042

2043

2044 def f2n(f):

2045 '''Return C{n}, the I{3rd flattening} for a given I{flattening}.

2046

2047 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2048

2049 @return: The I{signed}, 3rd flattening (-1 <= C{float} < 1),

2050 M{f / (2 - f)}.

2051

2052 @note: The result is positive for I{oblate}, negative for I{prolate}

2053 or C{0} for I{near-spherical} ellipsoids.

2054

2055 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2056 html/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2057 <https://WikiPedia.org/wiki/Flattening>}.

2058 '''

2059 return Float(n=_0_0 if _spherical_f(f) else (f / float(_2_0 - f)))

2060

2061

2062 def n2e2(n):

2063 '''Return C{e2}, the I{1st eccentricity squared} for a given I{3rd flattening}.

2064

2065 @arg n: The 3rd flattening (-1 <= C{scalar} < 1).

2066

2067 @return: The I{signed}, (1st) eccentricity I{squared} (C{float} or NINF),

2068 M{4 * n / (1 + n)**2}.

2069

2070 @note: The result is positive for I{oblate}, negative for I{prolate}

2071 or C{0} for I{near-spherical} ellipsoids.

2072

2073 @see: U{Flattening<https://WikiPedia.org/wiki/Flattening>}.

2074 '''

2075 t = (n + _1_0)**2

2076 return Float(e2=_0_0 if fabs(n) < EPS else

2077 (NINF if t < EPS else (_4_0 * n / t)))

2078

2079

2080 def n2f(n):

2081 '''Return C{f}, the I{flattening} for a given I{3rd flattening}.

2082

2083 @arg n: The 3rd flattening (-1 <= C{scalar} < 1).

2084

2085 @return: The flattening (C{float} or NINF), M{2 * n / (1 + n)}.

2086

2087 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2088 html/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2089 <https://WikiPedia.org/wiki/Flattening>}.

2090 '''

2091 t = n + _1_0

2092 f = 0 if fabs(n) < EPS else (NINF if t < EPS else (_2_0 * n / t))

2093 return _f_0_0 if _spherical_f(f) else Float(f=f)

2094

2095

2096 def n2f_(n):

2097 '''Return C{f_}, the I{inverse flattening} for a given I{3rd flattening}.

2098

2099 @arg n: The 3rd flattening (-1 <= C{scalar} < 1).

2100

2101 @return: The inverse flattening (C{float} or C{0}), M{1 / f}.

2102

2103 @see: L{n2f} and L{f2f_}.

2104 '''

2105 return f2f_(n2f(n))

2106

2107

2108 def _normalTo3(px, py, E): # in .height4 above

2109 '''(INTERNAL) Nearest point on a 2-D ellipse in 1st quadrant.

2110

2111 @see: Functions C{.triaxial._normalTo4} and C{-To5}.

2112 '''

2113 a, b = E.a, E.b

2114 if min(px, py, a, b) < EPS0:

2115 raise _AssertionError(px=px, py=py, a=a, b=b, E=E)

2116

2117 a2 = a - b * E.b_a

2118 b2 = b - a * E.a_b

2119 tx = ty = _SQRT2_2

2120 for i in range(10): # max 5

2121 ex = a2 * tx**3

2122 ey = b2 * ty**3

2123

2124 qx = px - ex

2125 qy = py - ey

2126 q = hypot(qx, qy)

2127 if q < EPS0: # PYCHOK no cover

2128 break

2129 r = hypot(ex - tx * a, ey - ty * b) / q

2130

2131 sx, tx = tx, min(_1_0, max(0, (ex + qx * r) / a))

2132 sy, ty = ty, min(_1_0, max(0, (ey + qy * r) / b))

2133 t = hypot(ty, tx)

2134 if t < EPS0: # PYCHOK no cover

2135 break

2136 tx = tx / t # /= chokes PyChecker

2137 ty = ty / t

2138 if max(fabs(sx - tx), fabs(sy - ty)) < EPS:

2139 break

2140

2141 tx *= a / px

2142 ty *= b / py

2143 return tx, ty, i # x and y as fractions

2144

2145

2146 class Ellipsoids(_NamedEnum):

2147 '''(INTERNAL) L{Ellipsoid} registry, I{must} be a sub-class

2148 to accommodate the L{_LazyNamedEnumItem} properties.

2149 '''

2150 def _Lazy(self, a, b, f_, **kwds):

2151 '''(INTERNAL) Instantiate the L{Ellipsoid}.

2152 '''

2153 return Ellipsoid(a, b=b, f_=f_, **kwds)

2154

2155 Ellipsoids = Ellipsoids(Ellipsoid) # PYCHOK singleton

2156'''Some pre-defined L{Ellipsoid}s, all I{lazily} instantiated.'''

2157# <https://www.GNU.org/software/gama/manual/html_node/Supported-ellipsoids.html>

2158# <https://GSSC.ESA.int/navipedia/index.php/Reference_Frames_in_GNSS>

2159# <https://kb.OSU.edu/dspace/handle/1811/77986>

2160# <https://www.IBM.com/docs/en/db2/11.5?topic=systems-supported-spheroids>

2161# <https://w3.Energistics.org/archive/Epicentre/Epicentre_v3.0/DataModel/LogicalDictionary/StandardValues/ellipsoid.html>

2162# <https://GitHub.com/locationtech/proj4j/blob/master/src/main/java/org/locationtech/proj4j/datum/Ellipsoid.java>

2163 Ellipsoids._assert( # <https://WikiPedia.org/wiki/Earth_ellipsoid>

2164 Airy1830 = _lazy(_Airy1830_, *_T(6377563.396, _0_0, 299.3249646)), # b=6356256.909

2165 AiryModified = _lazy(_AiryModified_, *_T(6377340.189, _0_0, 299.3249646)), # b=6356034.448

2166# APL4_9 = _lazy('APL4_9', *_T(6378137.0, _0_0, 298.24985392)), # Appl. Phys. Lab. 1965

2167# ANS = _lazy('ANS', *_T(6378160.0, _0_0, 298.25)), # Australian Nat. Spheroid

2168# AN_SA96 = _lazy('AN_SA96', *_T(6378160.0, _0_0, 298.24985392)), # Australian Nat. South America

2169 Australia1966 = _lazy('Australia1966', *_T(6378160.0, _0_0, 298.25)), # b=6356774.7192

2170 ATS1977 = _lazy('ATS1977', *_T(6378135.0, _0_0, 298.257)), # "Average Terrestrial System"

2171 Bessel1841 = _lazy(_Bessel1841_, *_T(6377397.155, 6356078.962818, 299.152812797)),

2172 BesselModified = _lazy('BesselModified', *_T(6377492.018, _0_0, 299.1528128)),

2173# BesselNamibia = _lazy('BesselNamibia', *_T(6377483.865, _0_0, 299.1528128)),

2174 CGCS2000 = _lazy('CGCS2000', *_T(R_MA, _0_0, 298.257222101)), # BeiDou Coord System (BDC)

2175# Clarke1858 = _lazy('Clarke1858', *_T(6378293.639, _0_0, 294.260676369)),

2176 Clarke1866 = _lazy(_Clarke1866_, *_T(6378206.4, 6356583.8, 294.978698214)),

2177 Clarke1880 = _lazy('Clarke1880', *_T(6378249.145, 6356514.86954978, 293.465)),

2178 Clarke1880IGN = _lazy(_Clarke1880IGN_, *_T(6378249.2, 6356515.0, 293.466021294)),

2179 Clarke1880Mod = _lazy('Clarke1880Mod', *_T(6378249.145, 6356514.96639549, 293.466307656)), # aka Clarke1880Arc

2180 CPM1799 = _lazy('CPM1799', *_T(6375738.7, 6356671.92557493, 334.39)), # Comm. des Poids et Mesures

2181 Delambre1810 = _lazy('Delambre1810', *_T(6376428.0, 6355957.92616372, 311.5)), # Belgium

2182 Engelis1985 = _lazy('Engelis1985', *_T(6378136.05, 6356751.32272154, 298.2566)),

2183# Everest1830 = _lazy('Everest1830', *_T(6377276.345, _0_0, 300.801699997)),

2184# Everest1948 = _lazy('Everest1948', *_T(6377304.063, _0_0, 300.801699997)),

2185# Everest1956 = _lazy('Everest1956', *_T(6377301.243, _0_0, 300.801699997)),

2186 Everest1969 = _lazy('Everest1969', *_T(6377295.664, 6356094.667915, 300.801699997)),

2187 Everest1975 = _lazy('Everest1975', *_T(6377299.151, 6356098.14512013, 300.8017255)),

2188 Fisher1968 = _lazy('Fisher1968', *_T(6378150.0, 6356768.33724438, 298.3)),

2189# Fisher1968Mod = _lazy('Fisher1968Mod', *_T(6378155.0, _0_0, 298.3)),

2190 GEM10C = _lazy('GEM10C', *_T(R_MA, 6356752.31424783, 298.2572236)),

2191 GPES = _lazy('GPES', *_T(6378135.0, 6356750.0, _0_0)), # "Gen. Purpose Earth Spheroid"

2192 GRS67 = _lazy('GRS67', *_T(6378160.0, _0_0, 298.247167427)), # Lucerne b=6356774.516

2193# GRS67Truncated = _lazy('GRS67Truncated', *_T(6378160.0, _0_0, 298.25)),

2194 GRS80 = _lazy(_GRS80_, *_T(R_MA, 6356752.314140347, 298.257222101)), # IUGG, ITRS, ETRS89

2195# Hayford1924 = _lazy('Hayford1924', *_T(6378388.0, 6356911.94612795, None)), # aka Intl1924 f_=297

2196 Helmert1906 = _lazy('Helmert1906', *_T(6378200.0, 6356818.16962789, 298.3)),

2197# Hough1960 = _lazy('Hough1960', *_T(6378270.0, _0_0, 297.0)),

2198 IAU76 = _lazy('IAU76', *_T(6378140.0, _0_0, 298.257)), # Int'l Astronomical Union

2199 IERS1989 = _lazy('IERS1989', *_T(6378136.0, _0_0, 298.257)), # b=6356751.302

2200 IERS1992TOPEX = _lazy('IERS1992TOPEX', *_T(6378136.3, 6356751.61659215, 298.257223563)), # IERS/TOPEX/Poseidon/McCarthy

2201 IERS2003 = _lazy('IERS2003', *_T(6378136.6, 6356751.85797165, 298.25642)),

2202 Intl1924 = _lazy(_Intl1924_, *_T(6378388.0, _0_0, 297.0)), # aka Hayford b=6356911.9462795

2203 Intl1967 = _lazy('Intl1967', *_T(6378157.5, 6356772.2, 298.24961539)), # New Int'l

2204 Krassovski1940 = _lazy(_Krassovski1940_, *_T(6378245.0, 6356863.01877305, 298.3)), # spelling

2205 Krassowsky1940 = _lazy(_Krassowsky1940_, *_T(6378245.0, 6356863.01877305, 298.3)), # spelling

2206# Kaula = _lazy('Kaula', *_T(6378163.0, _0_0, 298.24)), # Kaula 1961

2207# Lerch = _lazy('Lerch', *_T(6378139.0, _0_0, 298.257)), # Lerch 1979

2208 Maupertuis1738 = _lazy('Maupertuis1738', *_T(6397300.0, 6363806.28272251, 191.0)), # France

2209 Mercury1960 = _lazy('Mercury1960', *_T(6378166.0, 6356784.28360711, 298.3)),

2210 Mercury1968Mod = _lazy('Mercury1968Mod', *_T(6378150.0, 6356768.33724438, 298.3)),

2211# MERIT = _lazy('MERIT', *_T(6378137.0, _0_0, 298.257)), # MERIT 1983

2212# NWL10D = _lazy('NWL10D', *_T(6378135.0, _0_0, 298.26)), # Naval Weapons Lab.

2213 NWL1965 = _lazy('NWL1965', *_T(6378145.0, 6356759.76948868, 298.25)), # Naval Weapons Lab.

2214# NWL9D = _lazy('NWL9D', *_T(6378145.0, 6356759.76948868, 298.25)), # NWL1965

2215 OSU86F = _lazy('OSU86F', *_T(6378136.2, 6356751.51693008, 298.2572236)),

2216 OSU91A = _lazy('OSU91A', *_T(6378136.3, 6356751.6165948, 298.2572236)),

2217# Plessis1817 = _lazy('Plessis1817', *_T(6397523.0, 6355863.0, 153.56512242)), # XXX incorrect?

2218 Plessis1817 = _lazy('Plessis1817', *_T(6376523.0, 6355862.93325557, 308.64)), # XXX IGN France 1972

2219 PZ90 = _lazy('PZ90', *_T(6378136.0, _0_0, 298.257839303)), # GLOSNASS PZ-90 and PZ-90.11

2220# SEAsia = _lazy('SEAsia', *_T(6378155.0, _0_0, 298.3)), # SouthEast Asia

2221 SGS85 = _lazy('SGS85', *_T(6378136.0, 6356751.30156878, 298.257)), # Soviet Geodetic System

2222 SoAmerican1969 = _lazy('SoAmerican1969', *_T(6378160.0, 6356774.71919531, 298.25)), # South American

2223 Struve1860 = _lazy('Struve1860', *_T(6378298.3, 6356657.14266956, 294.73)),

2224# Walbeck = _lazy('Walbeck', *_T(6376896.0, _0_0, 302.78)),

2225# WarOffice = _lazy('WarOffice', *_T(6378300.0, _0_0, 296.0)),

2226 WGS60 = _lazy('WGS60', *_T(6378165.0, 6356783.28695944, 298.3)),

2227 WGS66 = _lazy('WGS66', *_T(6378145.0, 6356759.76948868, 298.25)),

2228 WGS72 = _lazy(_WGS72_, *_T(6378135.0, _0_0, 298.26)), # b=6356750.52

2229 WGS84 = _lazy(_WGS84_, *_T(R_MA, _0_0, _f__WGS84)), # GPS b=6356752.3142451793

2230# Prolate = _lazy('Prolate', *_T(6356752.3, R_MA, _0_0)),

2231 Sphere = _lazy(_Sphere_, *_T(R_M, R_M, _0_0)), # pseudo

2232 SphereAuthalic = _lazy('SphereAuthalic', *_T(R_FM, R_FM, _0_0)), # pseudo

2233 SpherePopular = _lazy('SpherePopular', *_T(R_MA, R_MA, _0_0)) # EPSG:3857 Spheroid

2234)

2235

2236

2237 if __name__ == '__main__':

2238

2239 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_

2240 from pygeodesy.lazily import printf

2241 from pygeodesy.named import nameof

2242

2243 for E in (Ellipsoids.WGS84, Ellipsoids.GRS80, # NAD83,

2244 Ellipsoids.Sphere, Ellipsoids.SpherePopular,

2245 Ellipsoid(Ellipsoids.WGS84.b, Ellipsoids.WGS84.a, name='_Prolate')):

2246 e = f2n(E.f) - E.n

2247 printf('# %s: %s', _DOT_('Ellipsoids', E.name), E.toStr(prec=10), nl=1)

2248 printf('# e=%s, f_=%s, f=%s, n=%s (%s)', fstr(E.e, prec=13, fmt=Fmt.e),

2249 fstr(E.f_, prec=13, fmt=Fmt.e),

2250 fstr(E.f, prec=13, fmt=Fmt.e),

2251 fstr(E.n, prec=13, fmt=Fmt.e),

2252 fstr(e, prec=9, fmt=Fmt.e))

2253 printf('# %s %s', Ellipsoid.AlphaKs.name, fstr(E.AlphaKs, prec=20))

2254 printf('# %s %s', Ellipsoid.BetaKs.name, fstr(E.BetaKs, prec=20))

2255 printf('# %s %s', nameof(Ellipsoid.KsOrder), E.KsOrder) # property

2256

2257 # __doc__ of this file, force all into registry

2258 t = [NN] + Ellipsoids.toRepr(all=True, asorted=True).split(_NL_)

2259 printf(_NLATvar_.join(i.strip(_COMMA_) for i in t))

2260

2261# **) MIT License

2262#

2263# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.

2264#

2265# Permission is hereby granted, free of charge, to any person obtaining a

2266# copy of this software and associated documentation files (the "Software"),

2267# to deal in the Software without restriction, including without limitation

2268# the rights to use, copy, modify, merge, publish, distribute, sublicense,

2269# and/or sell copies of the Software, and to permit persons to whom the

2270# Software is furnished to do so, subject to the following conditions:

2271#

2272# The above copyright notice and this permission notice shall be included

2273# in all copies or substantial portions of the Software.

2274#

2275# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS

2276# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

2277# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL

2278# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR

2279# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,

2280# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR

2281# OTHER DEALINGS IN THE SOFTWARE.

2282

2283# % python3 -m pygeodesy.ellipsoids

2284

2285# 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, 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

2286# e=8.1819190842622e-02, f_=2.98257223563e+02, f=3.3528106647475e-03, n=1.6792203863837e-03 (0.0e+00)

2287# AlphaKs 0.00083773182062446994, 0.00000076085277735725, 0.00000000119764550324, 0.00000000000242917068, 0.00000000000000571182, 0.0000000000000000148, 0.00000000000000000004, 0.0

2288# BetaKs 0.00083773216405794875, 0.0000000590587015222, 0.00000000016734826653, 0.00000000000021647981, 0.00000000000000037879, 0.00000000000000000072, 0.0, 0.0

2289# KsOrder 8

2290

2291# Ellipsoids.GRS80: name='GRS80', a=6378137, b=6356752.3141403468, f_=298.257222101, f=0.0033528107, f2=0.0033640898, n=0.0016792204, e=0.081819191, e2=0.00669438, 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

2292# e=8.1819191042833e-02, f_=2.98257222101e+02, f=3.3528106811837e-03, n=1.6792203946295e-03 (0.0e+00)

2293# AlphaKs 0.00083773182472890429, 0.00000076085278481561, 0.00000000119764552086, 0.00000000000242917073, 0.00000000000000571182, 0.0000000000000000148, 0.00000000000000000004, 0.0

2294# BetaKs 0.0008377321681623882, 0.00000005905870210374, 0.000000000167348269, 0.00000000000021647982, 0.00000000000000037879, 0.00000000000000000072, 0.0, 0.0

2295# KsOrder 8

2296

2297# Ellipsoids.Sphere: name='Sphere', a=6371008.7714149999, b=6371008.7714149999, f_=0, f=0, f2=0, n=0, e=0, e2=0, 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

2298# e=0.0e+00, f_=0.0e+00, f=0.0e+00, n=0.0e+00 (0.0e+00)

2299# AlphaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2300# BetaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2301# KsOrder 8

2302

2303# Ellipsoids.SpherePopular: name='SpherePopular', a=6378137, b=6378137, f_=0, f=0, f2=0, n=0, e=0, e2=0, e22=0, e32=0, A=6378137, L=10018754.171394622, R1=6378137, R2=6378137, R3=6378137, Rbiaxial=6378137, Rtriaxial=6378137

2304# e=0.0e+00, f_=0.0e+00, f=0.0e+00, n=0.0e+00 (0.0e+00)

2305# AlphaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2306# BetaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2307# KsOrder 8

2308

2309# 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, 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

2310# e=8.2094437949696e-02, f_=-2.97257223563e+02, f=-3.3640898209765e-03, n=-1.6792203863837e-03 (0.0e+00)

2311# AlphaKs -0.00084149152514366627, 0.00000076653480614871, -0.00000000120934503389, 0.0000000000024576225, -0.00000000000000578863, 0.00000000000000001502, -0.00000000000000000004, 0.0

2312# BetaKs -0.00084149187224351817, 0.00000005842735196773, -0.0000000001680487236, 0.00000000000021706261, -0.00000000000000038002, 0.00000000000000000073, -0.0, 0.0

2313# KsOrder 8