1

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

3

4 u'''Ellipsoidal and spherical earth models.

5

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

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

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

9

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

11

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

13

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

62'''

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

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

65

66 from pygeodesy.basics import copysign0, 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, IntersectionError, _ValueError, _xkwds_not

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

73 hypot1, hypot2, sqrt3

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

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

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

77 _Intl1924_, _incompatible_, _invalid_, _Krassovski1940_, \

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

79 _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.05.23'

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): # see .albers._atanhee

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/C++/doc/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/C++/doc/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 hartzell4(self, pov, los=None):

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

1025 from a Point-Of-View in space.

1026

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

1028 or L{Vector3d}).

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

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

1031

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

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

1034 I{distance} C{h} from B{C{pov}} to C{(x, y, z)} along B{C{los}},

1035 all in C{meter}, conventionally.

1036

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

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

1039 outside this ellipsoid or points in an opposite

1040 direction.

1041

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

1043

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

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

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

1047 '''

1048 try:

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

1050 except Exception as x:

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

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

1053

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

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

1056 this ellipsoid's surface.

1057

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

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

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

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

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

1063

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

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

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

1067

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

1069

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

1071

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

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

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

1075 '''

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

1077 r = v.length

1078

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

1080 if r < EPS0: # EPS

1081 v = v.times(_0_0)

1082 h = -a

1083

1084 elif self.isSpherical:

1085 v = v.times(a / r)

1086 h = r - a

1087

1088 elif normal: # perpendicular to ellipsoid

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

1090 if x < EPS0: # PYCHOK no cover

1091 z = copysign0(b, v.z)

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

1093 h = y - b # polar

1094 elif y < EPS0: # PYCHOK no cover

1095 t = a / r

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

1097 h = x - a # equatorial

1098 else: # normal in 1st quadrant

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

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

1101 h = t.minus(v).length

1102

1103 else: # radial to ellipsoid's center

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

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

1106 v = v.times(t)

1107 h = r * (_1_0 - t)

1108

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

1110 name=self.height4.__name__)

1111

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

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

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

1115 '''

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

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

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

1119

1120 @Property_RO

1121 def isEllipsoidal(self):

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

1123 '''

1124 return self.f != 0

1125

1126 @Property_RO

1127 def isOblate(self):

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

1129 spherical otherwise.

1130 '''

1131 return self.f > 0

1132

1133 @Property_RO

1134 def isProlate(self):

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

1136 spherical otherwise.

1137 '''

1138 return self.f < 0

1139

1140 @Property_RO

1141 def isSpherical(self):

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

1143 '''

1144 return self.f == 0

1145

1146 def _Kseries(self, *AB8Ks):

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

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

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

1150

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

1152 coefficient tuples.

1153

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

1155

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

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

1158 '''

1159 k = self.KsOrder

1160 if self.n:

1161 ns = fpowers(self.n, k)

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

1163 else:

1164 ks = _0_0s(k)

1165 return ks

1166

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

1168 def KsOrder(self):

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

1170 '''

1171 return self._KsOrder

1172

1173 @KsOrder.setter # PYCHOK setter!

1174 def KsOrder(self, order):

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

1176

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

1178 '''

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

1180 raise _ValueError(order=order)

1181 if order != self._KsOrder:

1182 Ellipsoid.AlphaKs._update(self)

1183 Ellipsoid.BetaKs._update(self)

1184 self._KsOrder = order

1185

1186 @Property_RO

1187 def L(self):

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

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

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

1191 '''

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

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

1194

1195 @Property_RO

1196 def _L_90(self):

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

1198 '''

1199 return self.L / _90_0

1200

1201 def Llat(self, lat):

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

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

1204

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

1206

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

1208 hemisphere (C{meter}).

1209 '''

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

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

1212

1213 Lmeridian = Llat # meridional distance

1214

1215 @deprecated_Property_RO

1216 def majoradius(self): # PYCHOK no cover

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

1218 return self.a

1219

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

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

1222 along a parallel of (geodetic) latitude.

1223

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

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

1226

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

1228

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

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

1231

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

1233 '''

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

1235

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

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

1238 along a parallel of (geodetic) latitude.

1239

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

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

1242

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

1244

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

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

1247

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

1249 '''

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

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

1252

1253 @deprecated_Property_RO

1254 def minoradius(self): # PYCHOK no cover

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

1256 return self.b

1257

1258 @Property_RO

1259 def n(self):

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

1261 '''

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

1263

1264 flattening = f

1265 flattening1st = f

1266 flattening2nd = f2

1267 flattening3rd = n

1268

1269 polaradius = b # Rpolar

1270

1271# @Property_RO

1272# def Q(self):

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

1274#

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

1276#

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

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

1279# '''

1280# n = self.n

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

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

1283

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

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

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

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

1288

1289 @deprecated_Property_RO

1290 def quarteradius(self): # PYCHOK no cover

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

1292 return self.L

1293

1294 @Property_RO

1295 def R1(self):

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

1297

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

1299 and method C{Rgeometric}.

1300 '''

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

1302 return Radius(R1=r)

1303

1304 Rmean = R1

1305

1306 @Property_RO

1307 def R2(self):

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

1309

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

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

1312 '''

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

1314

1315 Rauthalic = R2

1316

1317# @Property_RO

1318# def R2(self):

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

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

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

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

1323

1324 @Property_RO

1325 def R2x(self):

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

1327

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

1329 '''

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

1331

1332 Rauthalicx = R2x

1333

1334 @Property_RO

1335 def R3(self):

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

1337

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

1339 '''

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

1341 return Radius(R3=r)

1342

1343 Rvolumetric = R3

1344

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

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

1347 along a parallel of (geodetic) latitude.

1348

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

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

1351

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

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

1354

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

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

1357

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

1359 '''

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

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

1362

1363 @Property_RO

1364 def Rbiaxial(self):

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

1366

1367 @see: C{Rtriaxial}

1368 '''

1369 a, b = self.a, self.b

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

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

1372 return Radius(Rbiaxial=q)

1373

1374 Requatorial = a # for consistent naming

1375

1376 def Rgeocentric(self, lat):

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

1378

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

1380

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

1382

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

1384

1385 @see: U{Geocentric Radius

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

1387 '''

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

1389 if a and self.f:

1390 if fabs(a) < PI_2:

1391 s2, c2 = _s2_c2(a)

1392 b2_a2_s2 = self.b2_a2 * s2

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

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

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

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

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

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

1399 else:

1400 r = self.b

1401 return Radius(Rgeocentric=r)

1402

1403 @Property_RO

1404 def Rgeometric(self):

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

1406

1407 @see: C{R1}.

1408 '''

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

1410 return Radius(Rgeometric=g)

1411

1412 @Property_RO

1413 def rhumbsolve(self):

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

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

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

1417 variable C{PYGEODESY_RHUMBSOLVE}.

1418 '''

1419 # if not self.isEllipsoidal:

1420 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1422

1423 @Property_RO

1424 def rhumbx(self):

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

1426 '''

1427 # if not self.isEllipsoidal:

1428 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1430

1431 def Rlat(self, lat):

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

1433

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

1435

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

1437

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

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

1440

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

1442

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

1444

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

1446

1447 @see: Method C{circle4}.

1448 '''

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

1450 r = self.a

1451 if self.f and lat: # .isEllipsoidal

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

1453 r = Radius(Rlat=r)

1454 return r

1455

1456 Rpolar = b # for consistent naming

1457

1458 @deprecated_Property_RO

1459 def Rquadratic(self): # PYCHOK no cover

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

1461 return self.Rbiaxial

1462

1463 @deprecated_Property_RO

1464 def Rr(self): # PYCHOK no cover

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

1466 return self.Rrectifying

1467

1468 @Property_RO

1469 def Rrectifying(self):

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

1471

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

1473 '''

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

1475 return Radius(Rrectifying=r)

1476

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

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

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

1480

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

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

1483 to use an alternate formula.

1484

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

1486

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

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

1489

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

1491 '''

1492 if not self.f: # .isSpherical

1493 n = self.a

1494 elif ca is None:

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

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

1497 elif ca: # derived from EcefYou.forward

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

1499 n = self.a / h

1500 elif sa:

1501 n = self.a2_b / fabs(sa)

1502 else:

1503 n = self.a

1504 return n

1505

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

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

1508 radii of curvature of (geodetic) latitude.

1509

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

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

1512

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

1514 the radii of curvature.

1515

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

1517

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

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

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

1521 Earth_radius#Radii_of_curvature>}.

1522 '''

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

1524

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

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

1527 curvature of (geodetic) latitude.

1528

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

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

1531

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

1533 radii of curvature.

1534

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

1536

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

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

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

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

1541 '''

1542 a = fabs(Phi(phi))

1543 if self.f:

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

1545 if r > EPS02:

1546 n = self.a / sqrt(r)

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

1548 else:

1549 m = n = _0_0 # PYCHOK attr

1550 else:

1551 m = n = self.a

1552 if scaled and a:

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

1554 return Curvature2Tuple(Radius(rocMeridional=m),

1555 Radius(rocPrimeVertical=n))

1556

1557 def rocBearing(self, lat, bearing):

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

1559 latitude and compass direction.

1560

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

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

1563

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

1565

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

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

1568

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

1570

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

1572 '''

1573 if self.f:

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

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

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

1577 c2 *= n / m

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

1579 s2 *= m / n

1580 n = m

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

1582 else:

1583 b = self.b # == self.a

1584 return Radius(rocBearing=b)

1585

1586 @Property_RO

1587 def rocEquatorial2(self):

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

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

1590

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

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

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

1594 '''

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

1596 Radius(rocPrimeVertical0=self.a))

1597

1598 def rocGauss(self, lat):

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

1600

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

1602

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

1604

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

1606

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

1608 Earth_radius#Radii_of_curvature>}

1609 '''

1610 # using ...

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

1612 # return sqrt(m * n)

1613 # ... requires 1 or 2 sqrt

1614 g = self.b

1615 if self.f:

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

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

1618 return Radius(rocGauss=g)

1619

1620 def rocMean(self, lat):

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

1622

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

1624

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

1626

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

1628

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

1630 Earth_radius#Radii_of_curvature>}

1631 '''

1632 if self.f:

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

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

1635 else:

1636 m = self.a

1637 return Radius(rocMean=m)

1638

1639 def rocMeridional(self, lat):

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

1641

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

1643

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

1645

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

1647

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

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

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

1651 '''

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

1653 self.rocEquatorial2.meridional

1654

1655 rocPolar = a2_b # synonymous

1656

1657 def rocPrimeVertical(self, lat):

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

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

1660

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

1662

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

1664

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

1666

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

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

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

1670 Earth_radius#Radii_of_curvature>}.

1671 '''

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

1673 self.rocEquatorial2.prime_vertical

1674

1675 rocTransverse = rocPrimeVertical # synonymous

1676

1677 @deprecated_Property_RO

1678 def Rs(self): # PYCHOK no cover

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

1680 return self.Rgeometric

1681

1682 @Property_RO

1683 def Rtriaxial(self):

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

1685

1686 @see: C{Rbiaxial}

1687 '''

1688 a, b = self.a, self.b

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

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

1691 return Radius(Rtriaxial=q)

1692

1693 def toEllipsoid2(self, name=NN):

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

1695

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

1697

1698 @see: Property C{a_f}.

1699 '''

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

1701

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

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

1704

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

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

1707 this ellipsoid's name.

1708

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

1710 '''

1711 E = Ellipsoid

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

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

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

1715

1716 def toTriaxial(self, name=NN):

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

1718

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

1720 pointing east and C{Z} pointing north.

1721

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

1723 '''

1724 T = self._triaxial

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

1726

1727 @Property_RO

1728 def _triaxial(self):

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

1730 '''

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

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

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

1734

1735 @Property_RO

1736 def volume(self):

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

1738

1739 @see: C{R3}.

1740 '''

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

1742

1743

1744 class Ellipsoid2(Ellipsoid):

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

1746 '''

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

1748 '''New L{Ellipsoid2}.

1749

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

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

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

1753

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

1755

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

1757

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

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

1760 '''

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

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

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

1764 else:

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

1766

1767

1768 def _spherical_a_b(a, b):

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

1770 '''

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

1772

1773

1774 def _spherical_f(f):

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

1776 '''

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

1778

1779

1780 def _spherical_f_(f_):

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

1782 '''

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

1784

1785

1786 def a_b2e(a, b):

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

1788

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

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

1791

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

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

1794

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

1796 '''

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

1798

1799

1800 def a_b2e2(a, b):

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

1802

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

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

1805

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

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

1808

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

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

1811 '''

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

1813

1814

1815 def a_b2e22(a, b):

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

1817

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

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

1820

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

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

1823

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

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

1826 '''

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

1828

1829

1830 def a_b2e32(a, b):

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

1832

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

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

1835

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

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

1838

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

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

1841 '''

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

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

1844

1845

1846 def a_b2f(a, b):

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

1848

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

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

1851

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

1853

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

1855 for I{near-spherical} ellipsoids.

1856 '''

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

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

1859

1860

1861 def a_b2f_(a, b):

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

1863

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

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

1866

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

1868

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

1870 for I{near-spherical} ellipsoids.

1871 '''

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

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

1874

1875

1876 def a_b2f2(a, b):

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

1878

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

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

1881

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

1883

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

1885 for I{near-spherical} ellipsoids.

1886 '''

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

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

1889

1890

1891 def a_b2n(a, b):

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

1893

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

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

1896

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

1898

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

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

1901 '''

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

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

1904

1905

1906 def a_f2b(a, f):

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

1908

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

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

1911

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

1913 '''

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

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

1916

1917

1918 def a_f_2b(a, f_):

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

1920

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

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

1923

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

1925 '''

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

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

1928

1929

1930 def b_f2a(b, f):

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

1932

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

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

1935

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

1937 '''

1938 t = _1_0 - f

1939 a = b if fabs(t < EPS) else (b / t)

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

1941

1942

1943 def b_f_2a(b, f_):

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

1945

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

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

1948

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

1950 '''

1951 t = f_ - _1_0

1952 a = b if _spherical_f_(f_) or fabs(t - f_) < EPS \

1953 or fabs(t) < EPS else (b * f_ / t)

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

1955

1956

1957 def e2f(e):

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

1959

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

1961

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

1963

1964 @see: Function L{e22f}.

1965 '''

1966 return e22f(e**2)

1967

1968

1969 def e22f(e2):

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

1971

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

1973

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

1975 '''

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

1977

1978

1979 def f2e2(f):

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

1981

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

1983

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

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

1986

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

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

1989

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

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

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

1993 '''

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

1995

1996

1997 def f2e22(f):

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

1999

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

2001

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

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

2004

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

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

2007

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

2009 html/classGeographicLib_1_1Ellipsoid.html>}.

2010 '''

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

2012 t = (_1_0 - f)**2

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

2014

2015

2016 def f2e32(f):

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

2018

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

2020

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

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

2023

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

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

2026

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

2028 html/classGeographicLib_1_1Ellipsoid.html>}.

2029 '''

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

2031 e2 = f2e2(f)

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

2033

2034

2035 def f_2f(f_):

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

2037

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

2039

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

2041

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

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

2044 '''

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

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

2047

2048

2049 def f2f_(f):

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

2051

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

2053

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

2055

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

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

2058 '''

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

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

2061

2062

2063 def f2f2(f):

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

2065

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

2067

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

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{Eccentricity conversions<https://GeographicLib.SourceForge.io/

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

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

2076 '''

2077 t = _1_0 - f

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

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

2080

2081

2082 def f2n(f):

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

2084

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

2086

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

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

2089

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

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

2092

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

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

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

2096 '''

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

2098

2099

2100 def n2e2(n):

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

2102

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

2104

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

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

2107

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

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

2110

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

2112 '''

2113 t = (n + _1_0)**2

2114 return Float(e2=_0_0 if fabs(n) < EPS else

2115 (NINF if t < EPS else (_4_0 * n / t)))

2116

2117

2118 def n2f(n):

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

2120

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

2122

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

2124

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

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

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

2128 '''

2129 t = n + _1_0

2130 f = 0 if fabs(n) < EPS else (NINF if t < EPS else (_2_0 * n / t))

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

2132

2133

2134 def n2f_(n):

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

2136

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

2138

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

2140

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

2142 '''

2143 return f2f_(n2f(n))

2144

2145

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

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

2148

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

2150 '''

2151 a, b = E.a, E.b

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

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

2154

2155 a2 = a - b * E.b_a

2156 b2 = b - a * E.a_b

2157 tx = ty = _SQRT2_2

2158 for i in range(10): # max 5

2159 ex = a2 * tx**3

2160 ey = b2 * ty**3

2161

2162 qx = px - ex

2163 qy = py - ey

2164 q = hypot(qx, qy)

2165 if q < EPS0: # PYCHOK no cover

2166 break

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

2168

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

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

2171 t = hypot(ty, tx)

2172 if t < EPS0: # PYCHOK no cover

2173 break

2174 tx = tx / t # /= chokes PyChecker

2175 ty = ty / t

2176 if max(fabs(sx - tx), fabs(sy - ty)) < EPS:

2177 break

2178

2179 tx *= a / px

2180 ty *= b / py

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

2182

2183

2184 class Ellipsoids(_NamedEnum):

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

2186 to accommodate the L{_LazyNamedEnumItem} properties.

2187 '''

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

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

2190 '''

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

2192

2193 Ellipsoids = Ellipsoids(Ellipsoid) # PYCHOK singleton

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2232 GRS80 = _lazy(_GRS80_, *_T(R_MA, 6356752.314140347, 298.257222101)), # IUGG, ITRS, ETRS89

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2272)

2273

2274

2275 if __name__ == '__main__':

2276

2277 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_

2278 from pygeodesy.lazily import printf

2279 from pygeodesy.named import nameof

2280

2281 for E in (Ellipsoids.WGS84, Ellipsoids.GRS80, # NAD83,

2282 Ellipsoids.Sphere, Ellipsoids.SpherePopular,

2283 Ellipsoid(Ellipsoids.WGS84.b, Ellipsoids.WGS84.a, name='_Prolate')):

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

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

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

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

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

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

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

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

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

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

2294

2295 # __doc__ of this file, force all into registry

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

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

2298

2299# **) MIT License

2300#

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

2302#

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

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

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

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

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

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

2309#

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

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

2312#

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

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

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

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

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

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

2319# OTHER DEALINGS IN THE SOFTWARE.

2320

2321# % python3 -m pygeodesy.ellipsoids

2322

2323# Ellipsoids.WGS84: name='WGS84', a=6378137, b=6356752.3142451793, f_=298.257223563, f=0.0033528107, f2=0.0033640898, n=0.0016792204, e=0.0818191908, e2=0.00669438, e21=0.99330562, e22=0.0067394967, e32=0.0033584313, A=6367449.1458234144, L=10001965.7293127216, R1=6371008.7714150595, R2=6371007.1809184738, R3=6371000.7900091587, Rbiaxial=6367453.6345163295, Rtriaxial=6372797.5559594007

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

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

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

2327# KsOrder 8

2328

2329# Ellipsoids.GRS80: name='GRS80', a=6378137, b=6356752.3141403468, f_=298.257222101, f=0.0033528107, f2=0.0033640898, n=0.0016792204, e=0.081819191, e2=0.00669438, e21=0.99330562, e22=0.0067394968, e32=0.0033584313, A=6367449.1457710434, L=10001965.7292304579, R1=6371008.7713801153, R2=6371007.1808835147, R3=6371000.7899741363, Rbiaxial=6367453.6344640013, Rtriaxial=6372797.5559332585

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

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

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

2333# KsOrder 8

2334

2335# Ellipsoids.Sphere: name='Sphere', a=6371008.7714149999, b=6371008.7714149999, f_=0, f=0, f2=0, n=0, e=0, e2=0, e21=1, e22=0, e32=0, A=6371008.7714149999, L=10007557.1761167478, R1=6371008.7714149999, R2=6371008.7714149999, R3=6371008.7714149999, Rbiaxial=6371008.7714149999, Rtriaxial=6371008.7714149999

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

2337# AlphaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2338# BetaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2339# KsOrder 8

2340

2341# Ellipsoids.SpherePopular: name='SpherePopular', a=6378137, b=6378137, f_=0, f=0, f2=0, n=0, e=0, e2=0, e21=1, e22=0, e32=0, A=6378137, L=10018754.171394622, R1=6378137, R2=6378137, R3=6378137, Rbiaxial=6378137, Rtriaxial=6378137

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

2343# AlphaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2344# BetaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2345# KsOrder 8

2346

2347# Ellipsoids._Prolate: name='_Prolate', a=6356752.3142451793, b=6378137, f_=-297.257223563, f=-0.0033640898, f2=-0.0033528107, n=-0.0016792204, e=0.0820944379, e2=-0.0067394967, e21=1.0067394967, e22=-0.00669438, e32=-0.0033584313, A=6367449.1458234144, L=10035500.5204500332, R1=6363880.5428301189, R2=6363878.9413582645, R3=6363872.5644020075, Rbiaxial=6367453.6345163295, Rtriaxial=6362105.2243882557

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

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

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

2351# KsOrder 8