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.62764496, 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.4212264, 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.63159783, 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.04438139, 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@var Ellipsoids.WGS84_NGS: Ellipsoid(name='WGS84_NGS', 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)

63'''

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

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

66

67 from pygeodesy.basics import copysign0, isbool, isint

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

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

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

71 _0_25, _3_0 # PYCHOK used!

72 from pygeodesy.errors import _AssertionError, IntersectionError, _ValueError, _xattr, _xkwds_not

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

74 hypot1, hypot2, sqrt3

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

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

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

78 _Intl1924_, _incompatible_, _invalid_, _Krassovski1940_, \

79 _Krassowsky1940_, _meridional_, _lat_, _negative_, _not_, \

80 _not_finite_, _prime_vertical_, _radius_, _Sphere_, _SPACE_, \

81 _vs_, _WGS72_, _WGS84_

82# from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS # from .named

83 from pygeodesy.named import _lazyNamedEnumItem as _lazy, _NamedEnum, _NamedEnumItem, \

84 _NamedTuple, _Pass, _ALL_LAZY, _MODS

85 from pygeodesy.namedTuples import Distance2Tuple, Vector3Tuple, Vector4Tuple

86 from pygeodesy.props import deprecated_Property_RO, Property_RO, property_doc_, \

87 deprecated_property_RO, property_RO

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

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

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

91 from pygeodesy.utily import atan1, atan1d, atan2b, degrees90, m2radians, radians2m, sincos2d

92

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

94

95__all__ = _ALL_LAZY.ellipsoids

96__version__ = '24.03.22'

97

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

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

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

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

102

103

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

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

106 '''

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

108

109

110 def _s2_c2(phi):

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

112 '''

113 if phi:

114 s2 = sin(phi)**2

115 if s2 > EPS:

116 c2 = _1_0 - s2

117 if c2 > EPS:

118 if c2 < EPS1:

119 return s2, c2

120 else:

121 return _1_0, _0_0 # phi == PI_2

122 return _0_0, _1_0 # phi == 0

123

124

125 class a_f2Tuple(_NamedTuple):

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

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

128

129 @see: Class L{Ellipsoid2}.

130 '''

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

132 _Units_ = (_Pass, _Pass)

133

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

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

136

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

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

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

140

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

142

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

144

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

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

147 '''

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

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

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

151 f = _f_0_0

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

153

154 @Property_RO

155 def b(self):

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

157 '''

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

159

160 def ellipsoid(self, name=NN):

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

162

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

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

165 '''

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

167

168 @Property_RO

169 def f_(self):

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

171 '''

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

173

174

175 class Circle4Tuple(_NamedTuple):

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

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

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

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

180

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

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

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

184 equals C{lat}.

185 '''

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

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

188

189

190 class Curvature2Tuple(_NamedTuple):

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

192 C{meter}, conventionally.

193 '''

194 _Names_ = (_meridional_, _prime_vertical_)

195 _Units_ = ( Meter, Meter)

196

197 @property_RO

198 def transverse(self):

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

200 '''

201 return self.prime_vertical

202

203

204 class Ellipsoid(_NamedEnumItem):

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

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

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

208 '''

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

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

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

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

213

214 _geodsolve = NN # means, use PYGEODESY_GEODSOLVE

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

216 _rhumbsolve = NN # means, use PYGEODESY_RHUMBSOLVE

217

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

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

220 the I{polar} radius or I{inverse flattening} or I{flattening}.

221

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

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

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

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

226 spherical).

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

228

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

230

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

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

233

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

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

236 '''

237 ff_ = f, f_ # assertion below

238 try:

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

240 if not _isfinite(a):

241 raise ValueError(_SPACE_(_a_, _not_finite_))

242

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

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

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

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

247 elif f is not None:

248 f = Float(f=f)

249 b = a_f2b(a, f)

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

251 elif f_:

252 f_ = Float(f_=f_)

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

254 f = f_2f(f_)

255 else: # only a, spherical

256 f_ = f = 0

257 b = a # superfluous

258

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

260 raise ValueError(_SPACE_(_f_, _invalid_))

261 if not _isfinite(b):

262 raise ValueError(_SPACE_(_b_, _not_finite_))

263

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

265 b = a

266 f = _f_0_0

267 f_ = _f__0_0

268

269 except (TypeError, ValueError) as x:

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

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

272 raise _ValueError(t, cause=x)

273

274 self._a = a

275 self._b = b

276 self._f = f

277 self._f_ = f_

278

279 self._register(Ellipsoids, name)

280

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

282 d = dict(eps=_TOL)

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

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

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

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

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

288

289 def __eq__(self, other):

290 '''Compare this and an other ellipsoid.

291

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

293

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

295 '''

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

297 self.a == other.a and

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

299

300 def __hash__(self):

301 return self._hash # memoized

302

303 @Property_RO

304 def a(self):

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

306 '''

307 return self._a

308

309 equatoradius = a # = Requatorial

310

311 @Property_RO

312 def a2(self):

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

314 '''

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

316

317 @Property_RO

318 def a2_(self):

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

320 '''

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

322

323 @Property_RO

324 def a_b(self):

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

326 '''

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

328

329 @Property_RO

330 def a2_b(self):

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

332

333 @see: U{Radii of Curvature

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

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

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

337

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

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

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

341 '''

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

343

344 @Property_RO

345 def a2_b2(self):

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

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

348 '''

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

350

351 @Property_RO

352 def a_f(self):

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

354 '''

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

356

357 @Property_RO

358 def A(self):

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

360

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

362 '''

363 A, n = self.a, self.n

364 if n:

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

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

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

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

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

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

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

372 return A

373

374 @Property_RO

375 def _albersCyl(self):

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

377 '''

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

379

380 @Property_RO

381 def AlphaKs(self):

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

383 '''

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

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

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

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

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

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

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

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

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

393 _T(1424729850961/743921418240)) # PYCHOK unaligned

394

395 @Property_RO

396 def area(self):

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

398

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

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

401 '''

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

403

404 @Property_RO

405 def areax(self):

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

407 accurate for very I{oblate} ellipsoids.

408

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

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

411 '''

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

413

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

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

416 '''

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

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

419 t = (v, _vs_, val)

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

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

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

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

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

425 eps=eps, f0=f0, **name_value))

426

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

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

429

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

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

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

433

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

435

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

437 C++/doc/classGeographicLib_1_1Ellipsoid.html>}, U{Authalic latitude

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

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

440 '''

441 if self.f:

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

443 lat = atan1d(f(tan(Phi_(lat)))) # PYCHOK attr

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

445

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

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

448

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

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

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

452

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

454

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

456 C++/doc/classGeographicLib_1_1Ellipsoid.html>}, U{Conformal latitude

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

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

459 '''

460 if self.f:

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

462 lat = atan1d(f(tan(Phi_(lat)))) # PYCHOK attr

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

464

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

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

467

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

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

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

471

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

473

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

475 C++/doc/classGeographicLib_1_1Ellipsoid.html>}, U{Geocentric latitude

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

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

478 '''

479 if self.f:

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

481 lat = atan1d(f * tan(Phi_(lat)))

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

483

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

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

486

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

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

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

490

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

492

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

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

495 thereof is the original geodetic latitude.

496

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

498 C++/doc/classGeographicLib_1_1Ellipsoid.html>}, U{Isometric latitude

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

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

501 '''

502 if self.f:

503 r = Phi_(lat, clip=0)

504 lat = degrees(atan1(self.es_tauf(sinh(r))) if inverse else

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

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

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

508

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

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

511

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

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

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

515

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

517

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

519 C++/doc/classGeographicLib_1_1Ellipsoid.html>}, U{Parametric latitude

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

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

522 '''

523 if self.f:

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

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

526

527 auxReduced = auxParametric # synonymous

528

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

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

531

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

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

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

535

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

537

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

539 C++/doc/classGeographicLib_1_1Ellipsoid.html>}, U{Rectifying latitude

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

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

542 '''

543 if self.f:

544 lat = Lat(lat)

545 if 0 < fabs(lat) < _90_0:

546 if inverse:

547 e = self._elliptic_e22

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

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

550 else:

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

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

553

554 @Property_RO

555 def b(self):

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

557 '''

558 return self._b

559

560 polaradius = b # = Rpolar

561

562 @Property_RO

563 def b_a(self):

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

565

566 @see: Property L{f1}.

567 '''

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

569

570 @Property_RO

571 def b2(self):

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

573 '''

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

575

576 @Property_RO

577 def b2_a(self):

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

579

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

581 '''

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

583

584 @Property_RO

585 def b2_a2(self):

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

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

588 '''

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

590

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

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

593 '''

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

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

596 return atan1d(s, c)

597

598 @Property_RO

599 def BetaKs(self):

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

601 '''

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

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

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

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

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

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

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

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

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

611 _T(191773887257/3719607091200)) # PYCHOK unaligned

612

613 @deprecated_Property_RO

614 def c(self): # PYCHOK no cover

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

616 return self.R2

617

618 @Property_RO

619 def c2(self):

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

621

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

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

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

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

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

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

628 '''

629 return self._c2f(False)

630

631 @Property_RO

632 def c2x(self):

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

634 ellipsoids.

635

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

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

638 '''

639 return self._c2f(True)

640

641 def _c2f(self, c2x):

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

643 '''

644 f, c2 = self.f, self.b2

645 if f:

646 e = self.e

647 if e > EPS0:

648 if f > 0: # .isOblate

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

650 elif f < 0: # .isProlate

651 c2 *= atan1(e) / e # XXX asin?

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

653 return c2

654

655 def circle4(self, lat):

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

657

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

659

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

661 instance.

662

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

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

665

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

667

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

669

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

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

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

673 and method C{Rlat}.

674 '''

675 lat = Lat(lat)

676 if lat:

677 b = lat

678 if fabs(lat) < _90_0:

679 if self.f:

680 b = self._beta(lat)

681 z, r = sincos2d(b)

682 r *= self.a

683 z *= self.b

684 else: # near-polar

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

686 else: # equator

687 r = self.a

688 z = lat = b = _0_0

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

690

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

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

693 along a parallel of (geodetic) latitude.

694

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

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

697

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

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

700

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

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

703

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

705 '''

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

707

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

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

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

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

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

713

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

715 and only between points not near-polar}.

716

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

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

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

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

721

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

723 in same units as this ellipsoid's axes.

724

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

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

727

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

729 '''

730 phi0 = Phi_(lat0=lat0)

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

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

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

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

735

736 @Property_RO

737 def e(self):

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

739

740 @see: Property L{es}.

741 '''

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

743

744 @deprecated_Property_RO

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

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

747 return self.e21

748

749# @Property_RO

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

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

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

753

754 @Property_RO

755 def e2(self):

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

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

758 '''

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

760

761 @Property_RO

762 def e2abs(self):

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

764 '''

765 return fabs(self.e2)

766

767 @Property_RO

768 def e21(self):

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

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

771 '''

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

773

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

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

776

777 @Property_RO

778 def e22(self):

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

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

781 '''

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

783

784 @Property_RO

785 def e22abs(self):

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

787 '''

788 return fabs(self.e22)

789

790 @Property_RO

791 def e32(self):

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

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

794 '''

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

796

797 @Property_RO

798 def e32abs(self):

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

800 '''

801 return fabs(self.e32)

802

803 @Property_RO

804 def e4(self):

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

806 '''

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

808

809 eccentricity = e # eccentricity

810# eccentricity2 = e2 # eccentricity squared

811 eccentricity1st2 = e2 # first eccentricity squared

812 eccentricity2nd2 = e22 # second eccentricity squared

813 eccentricity3rd2 = e32 # third eccentricity squared

814

815 def ecef(self, Ecef=None):

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

817

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

819

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

821

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

823

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

825 '''

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

827

828 @Property_RO

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

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

831 '''

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

833 return _MODS.elliptic.Elliptic(e12)

834

835 @Property_RO

836 def _elliptic_e22(self): # aka ._elliptic_ep2

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

838 '''

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

840

841 equatoradius = a # Requatorial

842

843 def e2s(self, s):

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

845

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

847

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

849

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

851 '''

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

853

854 def e2s2(self, s):

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

856

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

858

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

860

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

862 '''

863 r = _1_0

864 if self.e2:

865 try:

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

867 if r < 0:

868 raise ValueError(_negative_)

869 except (TypeError, ValueError) as x:

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

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

872 return r

873

874 @Property_RO

875 def es(self):

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

877

878 @see: Property L{e}.

879 '''

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

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

882

883 def es_atanh(self, x):

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

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

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

887

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

889

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

891 C++/doc/classGeographicLib_1_1Math.html>}.

892 '''

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

894

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

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

897 '''

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

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

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

901 _0_0)) # .isSpherical

902

903 @Property_RO

904 def es_c(self):

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

906 '''

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

908

909 def es_tauf(self, taup):

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

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

912

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

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

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

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

917 '''

918 t = Scalar(taup=taup)

919 if self.f: # .isEllipsoidal

920 a = fabs(t)

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

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

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

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

925 if fabs(d) < s:

926 break

927 t = Scalar(tauf=T)

928 return t

929

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

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

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

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

934 '''

935 e = self._1_e21

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

937 _tf2 = self._es_taupf2

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

939 a, h = _tf2(T)

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

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

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

943 yield T, i, d

944

945 def es_taupf(self, tau):

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

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

948

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

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

951 '''

952 t = Scalar(tau=tau)

953 if self.f: # .isEllipsoidal

954 t, _ = self._es_taupf2(t)

955 t = Scalar(taupf=t)

956 return t

957

958 def _es_taupf2(self, tau):

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

960 '''

961 if _isfinite(tau):

962 h = hypot1(tau)

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

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

965 else:

966 a, h = tau, INF

967 return a, h

968

969 @Property_RO

970 def _exp_es_atanh_1(self):

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

972 '''

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

974

975 @Property_RO

976 def f(self):

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

978 '''

979 return self._f

980

981 @Property_RO

982 def f_(self):

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

984 '''

985 return self._f_

986

987 @Property_RO

988 def f1(self):

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

990

991 @see: Property L{b_a}.

992 '''

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

994

995 @Property_RO

996 def f2(self):

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

998 '''

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

1000

1001 @deprecated_Property_RO

1002 def geodesic(self):

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

1004 return self.geodesicw

1005

1006 def geodesic_(self, exact=True):

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

1008

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

1010 otherwise a L{Geodesic}, L{GeodesicExact} or L{GeodesicSolve}

1011 instance for I{this} ellipsoid.

1012

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

1014

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

1016

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

1018 '''

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

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

1021 name=self.name)

1022 else:

1023 g = exact

1024 E = _xattr(g, ellipsoid=None)

1025 if not (E is self and isinstance(g, self._Geodesics)):

1026 raise _ValueError(exact=g, ellipsoid=E, txt=_not_(self.name))

1027 return g

1028

1029 @property_RO

1030 def _Geodesics(self):

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

1032 '''

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

1034 _MODS.geodesicx.GeodesicExact,

1035 _MODS.geodsolve.GeodesicSolve)

1036 return t

1037

1038 @property_RO

1039 def geodesicw(self):

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

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

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

1043 package is installed.

1044 '''

1045 # if not self.isEllipsoidal:

1046 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

1047 return _MODS.geodesicw.Geodesic(self)

1048

1049 @property_RO

1050 def geodesicx(self):

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

1052 '''

1053 # if not self.isEllipsoidal:

1054 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1056

1057 @property

1058 def geodsolve(self):

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

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

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

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

1063 '''

1064 # if not self.isEllipsoidal:

1065 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1067

1068 @geodsolve.setter # PYCHOK setter!

1069 def geodsolve(self, path):

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

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

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

1073 '''

1074 self._geodsolve = path

1075

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

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

1078 from a Point-Of-View in space.

1079

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

1081 or L{Vector3d}).

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

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

1084

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

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

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

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

1089

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

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

1092 outside this ellipsoid or points in an opposite

1093 direction.

1094

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

1096

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

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

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

1100 '''

1101 try:

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

1103 except Exception as x:

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

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

1106

1107 @Property_RO

1108 def _hash(self):

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

1110

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

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

1113 this ellipsoid's surface.

1114

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

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

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

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

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

1120

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

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

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

1124

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

1126

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

1128

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

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

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

1132 '''

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

1134 r = v.length

1135

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

1137 if r < EPS0: # EPS

1138 v = v.times(_0_0)

1139 h = -a

1140

1141 elif self.isSpherical:

1142 v = v.times(a / r)

1143 h = r - a

1144

1145 elif normal: # perpendicular to ellipsoid

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

1147 if x < EPS0: # PYCHOK no cover

1148 z = copysign0(b, v.z)

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

1150 h = y - b # polar

1151 elif y < EPS0: # PYCHOK no cover

1152 t = a / r

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

1154 h = x - a # equatorial

1155 else: # normal in 1st quadrant

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

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

1158 h = t.minus(v).length

1159

1160 else: # radial to ellipsoid's center

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

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

1163 v = v.times(t)

1164 h = r * (_1_0 - t)

1165

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

1167 name=self.height4.__name__)

1168

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

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

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

1172 '''

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

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

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

1176

1177 @Property_RO

1178 def isEllipsoidal(self):

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

1180 '''

1181 return self.f != 0

1182

1183 @Property_RO

1184 def isOblate(self):

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

1186 spherical otherwise.

1187 '''

1188 return self.f > 0

1189

1190 @Property_RO

1191 def isProlate(self):

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

1193 spherical otherwise.

1194 '''

1195 return self.f < 0

1196

1197 @Property_RO

1198 def isSpherical(self):

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

1200 '''

1201 return self.f == 0

1202

1203 def _Kseries(self, *AB8Ks):

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

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

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

1207

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

1209 coefficient tuples.

1210

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

1212

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

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

1215 '''

1216 k = self.KsOrder

1217 if self.n:

1218 ns = fpowers(self.n, k)

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

1220 else:

1221 ks = _0_0s(k)

1222 return ks

1223

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

1225 def KsOrder(self):

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

1227 '''

1228 return self._KsOrder

1229

1230 @KsOrder.setter # PYCHOK setter!

1231 def KsOrder(self, order):

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

1233

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

1235 '''

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

1237 raise _ValueError(order=order)

1238 if self._KsOrder != order:

1239 Ellipsoid.AlphaKs._update(self)

1240 Ellipsoid.BetaKs._update(self)

1241 self._KsOrder = order

1242

1243 @Property_RO

1244 def L(self):

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

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

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

1248 '''

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

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

1251

1252 def Llat(self, lat):

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

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

1255

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

1257

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

1259 hemisphere (C{meter}).

1260 '''

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

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

1263

1264 Lmeridian = Llat # meridional distance

1265

1266 @property_RO

1267 def _Lpd(self):

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

1269 '''

1270 return Meter(_Lpd=self.L / _90_0)

1271

1272 @property_RO

1273 def _Lpr(self):

1274 '''Get the I{quarter meridian} per radian (C{meter}), M{self.L / PI_2}.

1275 '''

1276 return Meter(_Lpr=self.L / PI_2)

1277

1278 @deprecated_Property_RO

1279 def majoradius(self): # PYCHOK no cover

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

1281 return self.a

1282

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

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

1285 along a parallel of (geodetic) latitude.

1286

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

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

1289

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

1291

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

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

1294

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

1296 '''

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

1298

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

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

1301 along a parallel of (geodetic) latitude.

1302

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

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

1305

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

1307

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

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

1310

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

1312 '''

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

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

1315

1316 @deprecated_Property_RO

1317 def minoradius(self): # PYCHOK no cover

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

1319 return self.b

1320

1321 @Property_RO

1322 def n(self):

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

1324 '''

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

1326

1327 flattening = f

1328 flattening1st = f

1329 flattening2nd = f2

1330 flattening3rd = n

1331

1332 polaradius = b # Rpolar

1333

1334# @Property_RO

1335# def Q(self):

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

1337#

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

1339#

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

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

1342# '''

1343# n = self.n

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

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

1346

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

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

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

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

1351

1352 @deprecated_Property_RO

1353 def quarteradius(self): # PYCHOK no cover

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

1355 return self.L

1356

1357 @Property_RO

1358 def R1(self):

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

1360

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

1362 and method C{Rgeometric}.

1363 '''

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

1365 return Radius(R1=r)

1366

1367 Rmean = R1

1368

1369 @Property_RO

1370 def R2(self):

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

1372

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

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

1375 '''

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

1377

1378 Rauthalic = R2

1379

1380# @Property_RO

1381# def R2(self):

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

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

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

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

1386

1387 @Property_RO

1388 def R2x(self):

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

1390

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

1392 '''

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

1394

1395 Rauthalicx = R2x

1396

1397 @Property_RO

1398 def R3(self):

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

1400

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

1402 '''

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

1404 return Radius(R3=r)

1405

1406 Rvolumetric = R3

1407

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

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

1410 along a parallel of (geodetic) latitude.

1411

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

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

1414

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

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

1417

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

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

1420

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

1422 '''

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

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

1425

1426 @Property_RO

1427 def Rbiaxial(self):

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

1429

1430 @see: C{Rtriaxial}

1431 '''

1432 a, b = self.a, self.b

1433 if b < a:

1434 b = sqrt(_0_5 + self.b2_a2 * _0_5) * a

1435 elif b > a:

1436 b *= sqrt(_0_5 + self.a2_b2 * _0_5)

1437 return Radius(Rbiaxial=b)

1438

1439 Requatorial = a # for consistent naming

1440

1441 def Rgeocentric(self, lat):

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

1443

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

1445

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

1447

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

1449

1450 @see: U{Geocentric Radius

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

1452 '''

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

1454 if a and self.f:

1455 if fabs(a) < PI_2:

1456 s2, c2 = _s2_c2(a)

1457 b2_a2_s2 = self.b2_a2 * s2

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

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

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

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

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

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

1464 else:

1465 r = self.b

1466 return Radius(Rgeocentric=r)

1467

1468 @Property_RO

1469 def Rgeometric(self):

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

1471

1472 @see: C{R1}.

1473 '''

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

1475 return Radius(Rgeometric=g)

1476

1477 def rhumb_(self, exact=True):

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

1479

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

1481 otherwise a L{Rhumb}, L{RhumbAux} or L{RhumbSolve} instance

1482 for I{this} ellipsoid.

1483

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

1485

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

1487

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

1489 '''

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

1491 r = _MODS.rhumb.ekx.Rhumb(self, exact=exact, name=self.name)

1492 else:

1493 r = exact

1494 E = _xattr(r, ellipsoid=None)

1495 if not (E is self and isinstance(r, self._Rhumbs)):

1496 raise _ValueError(exact=r, ellipsosid=E, txt=_not_(self.name))

1497 return r

1498

1499 @property_RO

1500 def rhumbaux(self):

1501 '''Get this ellipsoid's I{Auxiliary} C{rhumb.RhumbAux}.

1502 '''

1503 # if not self.isEllipsoidal:

1504 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

1505 return _MODS.rhumb.aux_.RhumbAux(self, name=self.name)

1506

1507 @property_RO

1508 def rhumbekx(self):

1509 '''Get this ellipsoid's I{Elliptic, Krüger} C{rhumb.Rhumb}.

1510 '''

1511 # if not self.isEllipsoidal:

1512 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

1513 return _MODS.rhumb.ekx.Rhumb(self, name=self.name)

1514

1515 @property_RO

1516 def _Rhumbs(self):

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

1518 '''

1519 p = _MODS.rhumb

1520 Ellipsoid._Rhumbs = t = (p.aux_.RhumbAux, # overwrite property_RO

1521 p.ekx.Rhumb, p.solve.RhumbSolve)

1522 return t

1523

1524 @property

1525 def rhumbsolve(self):

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

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

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

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

1530 '''

1531 # if not self.isEllipsoidal:

1532 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

1533 return _MODS.rhumb.solve.RhumbSolve(self, path=self._rhumbsolve, name=self.name)

1534

1535 @rhumbsolve.setter # PYCHOK setter!

1536 def rhumbsolve(self, path):

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

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

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

1540 '''

1541 self._rhumbsolve = path

1542

1543 @deprecated_property_RO

1544 def rhumbx(self):

1545 '''DEPRECATED on 2023.11.28, use property C{rhumbekx}. '''

1546 return self.rhumbekx

1547

1548 def Rlat(self, lat):

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

1550

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

1552

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

1554

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

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

1557

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

1559

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

1561

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

1563

1564 @see: Method C{circle4}.

1565 '''

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

1567 r = self.a

1568 if self.f and lat: # .isEllipsoidal

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

1570 r = Radius(Rlat=r)

1571 return r

1572

1573 Rpolar = b # for consistent naming

1574

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

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

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

1578

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

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

1581 to use an alternate formula.

1582

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

1584

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

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

1587

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

1589 '''

1590 if not self.f: # .isSpherical

1591 n = self.a

1592 elif ca is None:

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

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

1595 elif ca: # derived from EcefYou.forward

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

1597 n = self.a / h

1598 elif sa:

1599 n = self.a2_b / fabs(sa)

1600 else:

1601 n = self.a

1602 return n

1603

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

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

1606 radii of curvature of (geodetic) latitude.

1607

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

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

1610

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

1612 the radii of curvature.

1613

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

1615

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

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

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

1619 Earth_radius#Radii_of_curvature>}.

1620 '''

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

1622

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

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

1625 curvature of (geodetic) latitude.

1626

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

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

1629

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

1631 radii of curvature.

1632

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

1634

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

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

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

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

1639 '''

1640 a = fabs(Phi(phi))

1641 if self.f:

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

1643 if r > EPS02:

1644 n = self.a / sqrt(r)

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

1646 else:

1647 m = n = _0_0 # PYCHOK attr

1648 else:

1649 m = n = self.a

1650 if scaled and a:

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

1652 return Curvature2Tuple(Radius(rocMeridional=m),

1653 Radius(rocPrimeVertical=n))

1654

1655 def rocBearing(self, lat, bearing):

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

1657 latitude and compass direction.

1658

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

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

1661

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

1663

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

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

1666

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

1668

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

1670 '''

1671 if self.f:

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

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

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

1675 c2 *= n / m

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

1677 s2 *= m / n

1678 n = m

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

1680 else:

1681 b = self.b # == self.a

1682 return Radius(rocBearing=b)

1683

1684 @Property_RO

1685 def rocEquatorial2(self):

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

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

1688

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

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

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

1692 '''

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

1694 Radius(rocPrimeVertical0=self.a))

1695

1696 def rocGauss(self, lat):

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

1698

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

1700

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

1702

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

1704

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

1706 Earth_radius#Radii_of_curvature>}

1707 '''

1708 # using ...

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

1710 # return sqrt(m * n)

1711 # ... requires 1 or 2 sqrt

1712 g = self.b

1713 if self.f:

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

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

1716 return Radius(rocGauss=g)

1717

1718 def rocMean(self, lat):

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

1720

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

1722

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

1724

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

1726

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

1728 Earth_radius#Radii_of_curvature>}

1729 '''

1730 if self.f:

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

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

1733 else:

1734 m = self.a

1735 return Radius(rocMean=m)

1736

1737 def rocMeridional(self, lat):

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

1739

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

1741

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

1743

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

1745

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

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

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

1749 '''

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

1751 self.rocEquatorial2.meridional

1752

1753 rocPolar = a2_b # synonymous

1754

1755 def rocPrimeVertical(self, lat):

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

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

1758

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

1760

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

1762

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

1764

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

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

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

1768 Earth_radius#Radii_of_curvature>}.

1769 '''

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

1771 self.rocEquatorial2.prime_vertical

1772

1773 rocTransverse = rocPrimeVertical # synonymous

1774

1775 @deprecated_Property_RO

1776 def Rquadratic(self): # PYCHOK no cover

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

1778 return self.Rbiaxial

1779

1780 @deprecated_Property_RO

1781 def Rr(self): # PYCHOK no cover

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

1783 return self.Rrectifying

1784

1785 @Property_RO

1786 def Rrectifying(self):

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

1788

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

1790 '''

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

1792 return Radius(Rrectifying=r)

1793

1794 @deprecated_Property_RO

1795 def Rs(self): # PYCHOK no cover

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

1797 return self.Rgeometric

1798

1799 @Property_RO

1800 def Rtriaxial(self):

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

1802

1803 @see: C{Rbiaxial}

1804 '''

1805 a, b = self.a, self.b

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

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

1808 return Radius(Rtriaxial=q)

1809

1810 def toEllipsoid2(self, name=NN):

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

1812

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

1814

1815 @see: Property C{a_f}.

1816 '''

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

1818

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

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

1821

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

1823 @kwarg terse: Limit the number of items (C{int}, 0...18).

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

1825 this ellipsoid's name.

1826

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

1828 '''

1829 E = Ellipsoid

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

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

1832 if terse:

1833 t = t[:terse]

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

1835

1836 def toTriaxial(self, name=NN):

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

1838

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

1840 pointing east and C{Z} pointing north.

1841

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

1843 '''

1844 T = self._triaxial

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

1846

1847 @property_RO

1848 def _triaxial(self):

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

1850 '''

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

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

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

1854

1855 @Property_RO

1856 def volume(self):

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

1858

1859 @see: C{R3}.

1860 '''

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

1862

1863

1864 class Ellipsoid2(Ellipsoid):

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

1866 '''

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

1868 '''New L{Ellipsoid2}.

1869

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

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

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

1873

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

1875

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

1877

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

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

1880 '''

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

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

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

1884 else:

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

1886

1887

1888 def _spherical_a_b(a, b):

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

1890 '''

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

1892

1893

1894 def _spherical_f(f):

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

1896 '''

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

1898

1899

1900 def _spherical_f_(f_):

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

1902 '''

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

1904

1905

1906 def a_b2e(a, b):

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

1908

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

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

1911

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

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

1914

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

1916 '''

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

1918

1919

1920 def a_b2e2(a, b):

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

1922

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

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

1925

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

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

1928

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

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

1931 '''

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

1933

1934

1935 def a_b2e22(a, b):

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

1937

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

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

1940

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

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

1943

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

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

1946 '''

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

1948

1949

1950 def a_b2e32(a, b):

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

1952

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

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

1955

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

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

1958

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

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

1961 '''

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

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

1964

1965

1966 def a_b2f(a, b):

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

1968

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

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

1971

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

1973

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

1975 for I{near-spherical} ellipsoids.

1976 '''

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

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

1979

1980

1981 def a_b2f_(a, b):

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

1983

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

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

1986

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

1988

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

1990 for I{near-spherical} ellipsoids.

1991 '''

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

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

1994

1995

1996 def a_b2f2(a, b):

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

1998

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

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

2001

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

2003

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

2005 for I{near-spherical} ellipsoids.

2006 '''

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

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

2009

2010

2011 def a_b2n(a, b):

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

2013

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

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

2016

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

2018

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

2020 or C{0} for I{near-spherical} ellipsoids.

2021 '''

2022 t = 0 if _spherical_a_b(a, b) else float(a - b)

2023 return Float(n=_0_0 if fabs(t) < EPS0 else (t / (a + b)))

2024

2025

2026 def a_f2b(a, f):

2027 '''Return C{b}, the I{polar} radius for a given I{equatorial} radius and I{flattening}.

2028

2029 @arg a: Equatorial radius (C{scalar} > 0).

2030 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2031

2032 @return: The polar radius (C{float}), M{a * (1 - f)}.

2033 '''

2034 b = a if _spherical_f(f) else (a * (_1_0 - f))

2035 return Radius_(b=a if _spherical_a_b(a, b) else b)

2036

2037

2038 def a_f_2b(a, f_):

2039 '''Return C{b}, the I{polar} radius for a given I{equatorial} radius and I{inverse flattening}.

2040

2041 @arg a: Equatorial radius (C{scalar} > 0).

2042 @arg f_: Inverse flattening (C{scalar} >>> 1).

2043

2044 @return: The polar radius (C{float}), M{a * (f_ - 1) / f_}.

2045 '''

2046 b = a if _spherical_f_(f_) else (a * (f_ - _1_0) / f_)

2047 return Radius_(b=a if _spherical_a_b(a, b) else b)

2048

2049

2050 def b_f2a(b, f):

2051 '''Return C{a}, the I{equatorial} radius for a given I{polar} radius and I{flattening}.

2052

2053 @arg b: Polar radius (C{scalar} > 0).

2054 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2055

2056 @return: The equatorial radius (C{float}), M{b / (1 - f)}.

2057 '''

2058 t = _1_0 - f

2059 a = b if fabs(t) < EPS0 else (b / t)

2060 return Radius_(a=b if _spherical_a_b(a, b) else a)

2061

2062

2063 def b_f_2a(b, f_):

2064 '''Return C{a}, the I{equatorial} radius for a given I{polar} radius and I{inverse flattening}.

2065

2066 @arg b: Polar radius (C{scalar} > 0).

2067 @arg f_: Inverse flattening (C{scalar} >>> 1).

2068

2069 @return: The equatorial radius (C{float}), M{b * f_ / (f_ - 1)}.

2070 '''

2071 t = f_ - _1_0

2072 a = b if _spherical_f_(f_) or fabs(t - f_) < EPS0 \

2073 or fabs(t) < EPS0 else (b * f_ / t)

2074 return Radius_(a=b if _spherical_a_b(a, b) else a)

2075

2076

2077 def e2f(e):

2078 '''Return C{f}, the I{flattening} for a given I{1st eccentricity}.

2079

2080 @arg e: The (1st) eccentricity (0 <= C{float} < 1)

2081

2082 @return: The flattening (C{scalar} or C{0}).

2083

2084 @see: Function L{e22f}.

2085 '''

2086 return e22f(e**2)

2087

2088

2089 def e22f(e2):

2090 '''Return C{f}, the I{flattening} for a given I{1st eccentricity squared}.

2091

2092 @arg e2: The (1st) eccentricity I{squared}, I{signed} (L{NINF} < C{float} < 1)

2093

2094 @return: The flattening (C{float} or C{0}), M{e2 / (sqrt(e2 - 1) + 1)}.

2095 '''

2096 return Float(f=e2 / (sqrt(_1_0 - e2) + _1_0)) if e2 else _f_0_0

2097

2098

2099 def f2e2(f):

2100 '''Return C{e2}, the I{1st eccentricity squared} for a given I{flattening}.

2101

2102 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2103

2104 @return: The I{signed}, (1st) eccentricity I{squared} (C{float} < 1),

2105 M{f * (2 - f)}.

2106

2107 @note: The result is positive for I{oblate}, negative for I{prolate}

2108 or C{0} for I{near-spherical} ellipsoids.

2109

2110 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2111 C++/doc/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2112 <https://WikiPedia.org/wiki/Flattening>}.

2113 '''

2114 return Float(e2=_0_0 if _spherical_f(f) else (f * (_2_0 - f)))

2115

2116

2117 def f2e22(f):

2118 '''Return C{e22}, the I{2nd eccentricity squared} for a given I{flattening}.

2119

2120 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2121

2122 @return: The I{signed}, 2nd eccentricity I{squared} (C{float} > -1 or

2123 C{INF}), M{f * (2 - f) / (1 - f)**2}.

2124

2125 @note: The result is positive for I{oblate}, negative for I{prolate}

2126 or C{0} for near-spherical ellipsoids.

2127

2128 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2129 C++/doc/classGeographicLib_1_1Ellipsoid.html>}.

2130 '''

2131 # e2 / (1 - e2) == f * (2 - f) / (1 - f)**2

2132 t = (_1_0 - f)**2

2133 return Float(e22=INF if t < EPS0 else (f2e2(f) / t)) # PYCHOK type

2134

2135

2136 def f2e32(f):

2137 '''Return C{e32}, the I{3rd eccentricity squared} for a given I{flattening}.

2138

2139 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2140

2141 @return: The I{signed}, 3rd eccentricity I{squared} (C{float}),

2142 M{f * (2 - f) / (1 + (1 - f)**2)}.

2143

2144 @note: The result is positive for I{oblate}, negative for I{prolate}

2145 or C{0} for I{near-spherical} ellipsoids.

2146

2147 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2148 C++/doc/classGeographicLib_1_1Ellipsoid.html>}.

2149 '''

2150 # e2 / (2 - e2) == f * (2 - f) / (1 + (1 - f)**2)

2151 e2 = f2e2(f)

2152 return Float(e32=e2 / (_2_0 - e2))

2153

2154

2155 def f_2f(f_):

2156 '''Return C{f}, the I{flattening} for a given I{inverse flattening}.

2157

2158 @arg f_: Inverse flattening (C{scalar} >>> 1).

2159

2160 @return: The flattening (C{scalar} or C{0}), M{1 / f_}.

2161

2162 @note: The result is positive for I{oblate}, negative for I{prolate}

2163 or C{0} for I{near-spherical} ellipsoids.

2164 '''

2165 f = 0 if _spherical_f_(f_) else _1_0 / f_

2166 return _f_0_0 if _spherical_f(f) else Float(f=f) # PYCHOK type

2167

2168

2169 def f2f_(f):

2170 '''Return C{f_}, the I{inverse flattening} for a given I{flattening}.

2171

2172 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2173

2174 @return: The inverse flattening (C{scalar} or C{0}), M{1 / f}.

2175

2176 @note: The result is positive for I{oblate}, negative for I{prolate}

2177 or C{0} for I{near-spherical} ellipsoids.

2178 '''

2179 f_ = 0 if _spherical_f(f) else _1_0 / f

2180 return _f__0_0 if _spherical_f_(f_) else Float(f_=f_) # PYCHOK type

2181

2182

2183 def f2f2(f):

2184 '''Return C{f2}, the I{2nd flattening} for a given I{flattening}.

2185

2186 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2187

2188 @return: The I{signed}, 2nd flattening (C{scalar} or C{INF}), M{f / (1 - f)}.

2189

2190 @note: The result is positive for I{oblate}, negative for I{prolate}

2191 or C{0} for I{near-spherical} ellipsoids.

2192

2193 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2194 C++/doc/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2195 <https://WikiPedia.org/wiki/Flattening>}.

2196 '''

2197 t = _1_0 - f

2198 return Float(f2=_0_0 if _spherical_f(f) else (INF if fabs(t) < EPS

2199 else (f / t))) # PYCHOK type

2200

2201

2202 def f2n(f):

2203 '''Return C{n}, the I{3rd flattening} for a given I{flattening}.

2204

2205 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2206

2207 @return: The I{signed}, 3rd flattening (-1 <= C{float} < 1),

2208 M{f / (2 - f)}.

2209

2210 @note: The result is positive for I{oblate}, negative for I{prolate}

2211 or C{0} for I{near-spherical} ellipsoids.

2212

2213 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2214 C++/doc/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2215 <https://WikiPedia.org/wiki/Flattening>}.

2216 '''

2217 return Float(n=_0_0 if _spherical_f(f) else (f / float(_2_0 - f)))

2218

2219

2220 def n2e2(n):

2221 '''Return C{e2}, the I{1st eccentricity squared} for a given I{3rd flattening}.

2222

2223 @arg n: The 3rd flattening (-1 <= C{scalar} < 1).

2224

2225 @return: The I{signed}, (1st) eccentricity I{squared} (C{float} or NINF),

2226 M{4 * n / (1 + n)**2}.

2227

2228 @note: The result is positive for I{oblate}, negative for I{prolate}

2229 or C{0} for I{near-spherical} ellipsoids.

2230

2231 @see: U{Flattening<https://WikiPedia.org/wiki/Flattening>}.

2232 '''

2233 t = (n + _1_0)**2

2234 return Float(e2=_0_0 if fabs(n) < EPS0 else

2235 (NINF if t < EPS0 else (_4_0 * n / t)))

2236

2237

2238 def n2f(n):

2239 '''Return C{f}, the I{flattening} for a given I{3rd flattening}.

2240

2241 @arg n: The 3rd flattening (-1 <= C{scalar} < 1).

2242

2243 @return: The flattening (C{scalar} or NINF), M{2 * n / (1 + n)}.

2244

2245 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2246 C++/doc/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2247 <https://WikiPedia.org/wiki/Flattening>}.

2248 '''

2249 t = n + _1_0

2250 f = 0 if fabs(n) < EPS0 else (NINF if t < EPS0 else (_2_0 * n / t))

2251 return _f_0_0 if _spherical_f(f) else Float(f=f)

2252

2253

2254 def n2f_(n):

2255 '''Return C{f_}, the I{inverse flattening} for a given I{3rd flattening}.

2256

2257 @arg n: The 3rd flattening (-1 <= C{scalar} < 1).

2258

2259 @return: The inverse flattening (C{scalar} or C{0}), M{1 / f}.

2260

2261 @see: L{n2f} and L{f2f_}.

2262 '''

2263 return f2f_(n2f(n))

2264

2265

2266 def _normalTo3(px, py, E, eps=EPS): # in .height4 above

2267 '''(INTERNAL) Nearest point on a 2-D ellipse in 1st quadrant.

2268

2269 @see: Functions C{pygeodesy.triaxial._normalTo4} and C{-To5}.

2270 '''

2271 a, b, e0 = E.a, E.b, EPS0

2272 if min(px, py, a, b) < e0:

2273 raise _AssertionError(px=px, py=py, a=a, b=b, E=E)

2274

2275 a2 = a - b * E.b_a

2276 b2 = b - a * E.a_b

2277 tx = ty = _SQRT2_2

2278 _a, _h = fabs, hypot

2279 for i in range(16): # max 5

2280 ex = a2 * tx**3

2281 ey = b2 * ty**3

2282

2283 qx = px - ex

2284 qy = py - ey

2285 q = _h(qx, qy)

2286 if q < e0: # PYCHOK no cover

2287 break

2288 r = _h(ex - tx * a,

2289 ey - ty * b) / q

2290

2291 sx, tx = tx, min(_1_0, max(0, (ex + qx * r) / a))

2292 sy, ty = ty, min(_1_0, max(0, (ey + qy * r) / b))

2293 t = _h(ty, tx)

2294 if t < e0: # PYCHOK no cover

2295 break

2296 tx = tx / t # /= chokes PyChecker

2297 ty = ty / t

2298 if _a(sx - tx) < eps and \

2299 _a(sy - ty) < eps:

2300 break

2301

2302 tx *= a / px

2303 ty *= b / py

2304 return tx, ty, i # x and y as fractions

2305

2306

2307 class Ellipsoids(_NamedEnum):

2308 '''(INTERNAL) L{Ellipsoid} registry, I{must} be a sub-class

2309 to accommodate the L{_LazyNamedEnumItem} properties.

2310 '''

2311 def _Lazy(self, a, b, f_, **kwds):

2312 '''(INTERNAL) Instantiate the L{Ellipsoid}.

2313 '''

2314 return Ellipsoid(a, b=b, f_=f_, **kwds)

2315

2316 Ellipsoids = Ellipsoids(Ellipsoid) # PYCHOK singleton

2317'''Some pre-defined L{Ellipsoid}s, all I{lazily} instantiated.'''

2318# <https://www.GNU.org/software/gama/manual/html_node/Supported-ellipsoids.html>

2319# <https://GSSC.ESA.int/navipedia/index.php/Reference_Frames_in_GNSS>

2320# <https://kb.OSU.edu/dspace/handle/1811/77986>

2321# <https://www.IBM.com/docs/en/db2/11.5?topic=systems-supported-spheroids>

2322# <https://w3.Energistics.org/archive/Epicentre/Epicentre_v3.0/DataModel/LogicalDictionary/StandardValues/ellipsoid.html>

2323# <https://GitHub.com/locationtech/proj4j/blob/master/src/main/java/org/locationtech/proj4j/datum/Ellipsoid.java>

2324 Ellipsoids._assert( # <https://WikiPedia.org/wiki/Earth_ellipsoid>

2325 Airy1830 = _lazy(_Airy1830_, *_T(6377563.396, _0_0, 299.3249646)), # b=6356256.909

2326 AiryModified = _lazy(_AiryModified_, *_T(6377340.189, _0_0, 299.3249646)), # b=6356034.448

2327# APL4_9 = _lazy('APL4_9', *_T(6378137.0, _0_0, 298.24985392)), # Appl. Phys. Lab. 1965

2328# ANS = _lazy('ANS', *_T(6378160.0, _0_0, 298.25)), # Australian Nat. Spheroid

2329# AN_SA96 = _lazy('AN_SA96', *_T(6378160.0, _0_0, 298.24985392)), # Australian Nat. South America

2330 Australia1966 = _lazy('Australia1966', *_T(6378160.0, _0_0, 298.25)), # b=6356774.7192

2331 ATS1977 = _lazy('ATS1977', *_T(6378135.0, _0_0, 298.257)), # "Average Terrestrial System"

2332 Bessel1841 = _lazy(_Bessel1841_, *_T(6377397.155, 6356078.962818, 299.152812797)),

2333 BesselModified = _lazy('BesselModified', *_T(6377492.018, _0_0, 299.1528128)),

2334# BesselNamibia = _lazy('BesselNamibia', *_T(6377483.865, _0_0, 299.1528128)),

2335 CGCS2000 = _lazy('CGCS2000', *_T(R_MA, _0_0, 298.257222101)), # BeiDou Coord System (BDC)

2336# Clarke1858 = _lazy('Clarke1858', *_T(6378293.639, _0_0, 294.260676369)),

2337 Clarke1866 = _lazy(_Clarke1866_, *_T(6378206.4, 6356583.8, 294.978698214)),

2338 Clarke1880 = _lazy('Clarke1880', *_T(6378249.145, 6356514.86954978, 293.465)),

2339 Clarke1880IGN = _lazy(_Clarke1880IGN_, *_T(6378249.2, 6356515.0, 293.466021294)),

2340 Clarke1880Mod = _lazy('Clarke1880Mod', *_T(6378249.145, 6356514.96639549, 293.466307656)), # aka Clarke1880Arc

2341 CPM1799 = _lazy('CPM1799', *_T(6375738.7, 6356671.92557493, 334.39)), # Comm. des Poids et Mesures

2342 Delambre1810 = _lazy('Delambre1810', *_T(6376428.0, 6355957.92616372, 311.5)), # Belgium

2343 Engelis1985 = _lazy('Engelis1985', *_T(6378136.05, 6356751.32272154, 298.2566)),

2344# Everest1830 = _lazy('Everest1830', *_T(6377276.345, _0_0, 300.801699997)),

2345# Everest1948 = _lazy('Everest1948', *_T(6377304.063, _0_0, 300.801699997)),

2346# Everest1956 = _lazy('Everest1956', *_T(6377301.243, _0_0, 300.801699997)),

2347 Everest1969 = _lazy('Everest1969', *_T(6377295.664, 6356094.667915, 300.801699997)),

2348 Everest1975 = _lazy('Everest1975', *_T(6377299.151, 6356098.14512013, 300.8017255)),

2349 Fisher1968 = _lazy('Fisher1968', *_T(6378150.0, 6356768.33724438, 298.3)),

2350# Fisher1968Mod = _lazy('Fisher1968Mod', *_T(6378155.0, _0_0, 298.3)),

2351 GEM10C = _lazy('GEM10C', *_T(R_MA, 6356752.31424783, 298.2572236)),

2352 GPES = _lazy('GPES', *_T(6378135.0, 6356750.0, _0_0)), # "Gen. Purpose Earth Spheroid"

2353 GRS67 = _lazy('GRS67', *_T(6378160.0, _0_0, 298.247167427)), # Lucerne b=6356774.516

2354# GRS67Truncated = _lazy('GRS67Truncated', *_T(6378160.0, _0_0, 298.25)),

2355 GRS80 = _lazy(_GRS80_, *_T(R_MA, 6356752.314140347, 298.25722210088)), # IUGG, ITRS, ETRS89

2356# Hayford1924 = _lazy('Hayford1924', *_T(6378388.0, 6356911.94612795, None)), # aka Intl1924 f_=297

2357 Helmert1906 = _lazy('Helmert1906', *_T(6378200.0, 6356818.16962789, 298.3)),

2358# Hough1960 = _lazy('Hough1960', *_T(6378270.0, _0_0, 297.0)),

2359 IAU76 = _lazy('IAU76', *_T(6378140.0, _0_0, 298.257)), # Int'l Astronomical Union

2360 IERS1989 = _lazy('IERS1989', *_T(6378136.0, _0_0, 298.257)), # b=6356751.302

2361 IERS1992TOPEX = _lazy('IERS1992TOPEX', *_T(6378136.3, 6356751.61659215, 298.257223563)), # IERS/TOPEX/Poseidon/McCarthy

2362 IERS2003 = _lazy('IERS2003', *_T(6378136.6, 6356751.85797165, 298.25642)),

2363 Intl1924 = _lazy(_Intl1924_, *_T(6378388.0, _0_0, 297.0)), # aka Hayford b=6356911.9462795

2364 Intl1967 = _lazy('Intl1967', *_T(6378157.5, 6356772.2, 298.24961539)), # New Int'l

2365 Krassovski1940 = _lazy(_Krassovski1940_, *_T(6378245.0, 6356863.01877305, 298.3)), # spelling

2366 Krassowsky1940 = _lazy(_Krassowsky1940_, *_T(6378245.0, 6356863.01877305, 298.3)), # spelling

2367# Kaula = _lazy('Kaula', *_T(6378163.0, _0_0, 298.24)), # Kaula 1961

2368# Lerch = _lazy('Lerch', *_T(6378139.0, _0_0, 298.257)), # Lerch 1979

2369 Maupertuis1738 = _lazy('Maupertuis1738', *_T(6397300.0, 6363806.28272251, 191.0)), # France

2370 Mercury1960 = _lazy('Mercury1960', *_T(6378166.0, 6356784.28360711, 298.3)),

2371 Mercury1968Mod = _lazy('Mercury1968Mod', *_T(6378150.0, 6356768.33724438, 298.3)),

2372# MERIT = _lazy('MERIT', *_T(6378137.0, _0_0, 298.257)), # MERIT 1983

2373# NWL10D = _lazy('NWL10D', *_T(6378135.0, _0_0, 298.26)), # Naval Weapons Lab.

2374 NWL1965 = _lazy('NWL1965', *_T(6378145.0, 6356759.76948868, 298.25)), # Naval Weapons Lab.

2375# NWL9D = _lazy('NWL9D', *_T(6378145.0, 6356759.76948868, 298.25)), # NWL1965

2376 OSU86F = _lazy('OSU86F', *_T(6378136.2, 6356751.51693008, 298.2572236)),

2377 OSU91A = _lazy('OSU91A', *_T(6378136.3, 6356751.6165948, 298.2572236)),

2378# Plessis1817 = _lazy('Plessis1817', *_T(6397523.0, 6355863.0, 153.56512242)), # XXX incorrect?

2379 Plessis1817 = _lazy('Plessis1817', *_T(6376523.0, 6355862.93325557, 308.64)), # XXX IGN France 1972

2380# Prolate = _lazy('Prolate', *_T(6356752.3, R_MA, _0_0)),

2381 PZ90 = _lazy('PZ90', *_T(6378136.0, _0_0, 298.257839303)), # GLOSNASS PZ-90 and PZ-90.11

2382# SEAsia = _lazy('SEAsia', *_T(6378155.0, _0_0, 298.3)), # SouthEast Asia

2383 SGS85 = _lazy('SGS85', *_T(6378136.0, 6356751.30156878, 298.257)), # Soviet Geodetic System

2384 SoAmerican1969 = _lazy('SoAmerican1969', *_T(6378160.0, 6356774.71919531, 298.25)), # South American

2385 Sphere = _lazy(_Sphere_, *_T(R_M, R_M, _0_0)), # pseudo

2386 SphereAuthalic = _lazy('SphereAuthalic', *_T(R_FM, R_FM, _0_0)), # pseudo

2387 SpherePopular = _lazy('SpherePopular', *_T(R_MA, R_MA, _0_0)), # EPSG:3857 Spheroid

2388 Struve1860 = _lazy('Struve1860', *_T(6378298.3, 6356657.14266956, 294.73)),

2389# Walbeck = _lazy('Walbeck', *_T(6376896.0, _0_0, 302.78)),

2390# WarOffice = _lazy('WarOffice', *_T(6378300.0, _0_0, 296.0)),

2391 WGS60 = _lazy('WGS60', *_T(6378165.0, 6356783.28695944, 298.3)),

2392 WGS66 = _lazy('WGS66', *_T(6378145.0, 6356759.76948868, 298.25)),

2393 WGS72 = _lazy(_WGS72_, *_T(6378135.0, _0_0, 298.26)), # b=6356750.52

2394 WGS84 = _lazy(_WGS84_, *_T(R_MA, _0_0, _f__WGS84)), # GPS b=6356752.3142451793

2395# U{NOAA/NOS/NGS/inverse<https://GitHub.com/noaa-ngs/inverse/blob/main/invers3d.f>}

2396 WGS84_NGS = _lazy('WGS84_NGS', *_T(R_MA, _0_0, 298.257222100882711243162836600094))

2397)

2398

2399_EWGS84 = Ellipsoids.WGS84 # (INTERNAL) shared

2400

2401 if __name__ == '__main__':

2402

2403 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_

2404 from pygeodesy import nameof, printf

2405

2406 for E in (_EWGS84, Ellipsoids.GRS80, # NAD83,

2407 Ellipsoids.Sphere, Ellipsoids.SpherePopular,

2408 Ellipsoid(_EWGS84.b, _EWGS84.a, name='_Prolate')):

2409 e = f2n(E.f) - E.n

2410 printf('# %s: %s', _DOT_('Ellipsoids', E.name), E.toStr(prec=10), nl=1)

2411 printf('# e=%s, f_=%s, f=%s, n=%s (%s)', fstr(E.e, prec=13, fmt=Fmt.e),

2412 fstr(E.f_, prec=13, fmt=Fmt.e),

2413 fstr(E.f, prec=13, fmt=Fmt.e),

2414 fstr(E.n, prec=13, fmt=Fmt.e),

2415 fstr(e, prec=9, fmt=Fmt.e))

2416 printf('# %s %s', Ellipsoid.AlphaKs.name, fstr(E.AlphaKs, prec=20))

2417 printf('# %s %s', Ellipsoid.BetaKs.name, fstr(E.BetaKs, prec=20))

2418 printf('# %s %s', nameof(Ellipsoid.KsOrder), E.KsOrder) # property

2419

2420 # __doc__ of this file, force all into registry

2421 t = [NN] + Ellipsoids.toRepr(all=True, asorted=True).split(_NL_)

2422 printf(_NLATvar_.join(i.strip(_COMMA_) for i in t))

2423

2424# % python3 -m pygeodesy.ellipsoids

2425

2426# 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.7293127235, R1=6371008.7714150595, R2=6371007.1809184738, R3=6371000.7900091587, Rbiaxial=6367453.6345163295, Rtriaxial=6372797.5559594007

2427# e=8.1819190842622e-02, f_=2.98257223563e+02, f=3.3528106647475e-03, n=1.6792203863837e-03 (0.0e+00)

2428# AlphaKs 0.00083773182062446994, 0.00000076085277735725, 0.00000000119764550324, 0.00000000000242917068, 0.00000000000000571182, 0.0000000000000000148, 0.00000000000000000004, 0.0

2429# BetaKs 0.00083773216405794875, 0.0000000590587015222, 0.00000000016734826653, 0.00000000000021647981, 0.00000000000000037879, 0.00000000000000000072, 0.0, 0.0

2430# KsOrder 8

2431

2432# Ellipsoids.GRS80: name='GRS80', a=6378137, b=6356752.3141403468, f_=298.2572221009, f=0.0033528107, f2=0.0033640898, n=0.0016792204, e=0.081819191, e2=0.00669438, e21=0.99330562, e22=0.0067394968, e32=0.0033584313, A=6367449.1457710434, L=10001965.7292304561, R1=6371008.7713801153, R2=6371007.1808835147, R3=6371000.7899741363, Rbiaxial=6367453.6344640013, Rtriaxial=6372797.5559332585

2433# e=8.1819191042833e-02, f_=2.9825722210088e+02, f=3.3528106811837e-03, n=1.6792203946295e-03 (0.0e+00)

2434# AlphaKs 0.00083773182472890429, 0.00000076085278481561, 0.00000000119764552086, 0.00000000000242917073, 0.00000000000000571182, 0.0000000000000000148, 0.00000000000000000004, 0.0

2435# BetaKs 0.0008377321681623882, 0.00000005905870210374, 0.000000000167348269, 0.00000000000021647982, 0.00000000000000037879, 0.00000000000000000072, 0.0, 0.0

2436# KsOrder 8

2437

2438# 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

2439# e=0.0e+00, f_=0.0e+00, f=0.0e+00, n=0.0e+00 (0.0e+00)

2440# AlphaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2441# BetaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2442# KsOrder 8

2443

2444# 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

2445# e=0.0e+00, f_=0.0e+00, f=0.0e+00, n=0.0e+00 (0.0e+00)

2446# AlphaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2447# BetaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2448# KsOrder 8

2449

2450# 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.5204500314, R1=6363880.5428301189, R2=6363878.9413582645, R3=6363872.5644020075, Rbiaxial=6367453.6345163295, Rtriaxial=6362105.2243882557

2451# e=8.2094437949696e-02, f_=-2.97257223563e+02, f=-3.3640898209765e-03, n=-1.6792203863837e-03 (0.0e+00)

2452# AlphaKs -0.00084149152514366627, 0.00000076653480614871, -0.00000000120934503389, 0.0000000000024576225, -0.00000000000000578863, 0.00000000000000001502, -0.00000000000000000004, 0.0

2453# BetaKs -0.00084149187224351817, 0.00000005842735196773, -0.0000000001680487236, 0.00000000000021706261, -0.00000000000000038002, 0.00000000000000000073, -0.0, 0.0

2454# KsOrder 8

2455

2456# **) MIT License

2457#

2458# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.

2459#

2460# Permission is hereby granted, free of charge, to any person obtaining a

2461# copy of this software and associated documentation files (the "Software"),

2462# to deal in the Software without restriction, including without limitation

2463# the rights to use, copy, modify, merge, publish, distribute, sublicense,

2464# and/or sell copies of the Software, and to permit persons to whom the

2465# Software is furnished to do so, subject to the following conditions:

2466#

2467# The above copyright notice and this permission notice shall be included

2468# in all copies or substantial portions of the Software.

2469#

2470# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS

2471# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

2472# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL

2473# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR

2474# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,

2475# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR

2476# OTHER DEALINGS IN THE SOFTWARE.