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.09'

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 html/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 html/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 html/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 html/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 html/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 html/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 html/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, 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 name: Override name (C{str}) or C{None} to exclude

1824 this ellipsoid's name.

1825

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

1827 '''

1828 E = Ellipsoid

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

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

1831 if Fmt.form: # terse form[at]

1832 t = t[:4]

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

1834

1835 def toTriaxial(self, name=NN):

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

1837

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

1839 pointing east and C{Z} pointing north.

1840

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

1842 '''

1843 T = self._triaxial

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

1845

1846 @property_RO

1847 def _triaxial(self):

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

1849 '''

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

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

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

1853

1854 @Property_RO

1855 def volume(self):

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

1857

1858 @see: C{R3}.

1859 '''

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

1861

1862

1863 class Ellipsoid2(Ellipsoid):

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

1865 '''

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

1867 '''New L{Ellipsoid2}.

1868

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

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

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

1872

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

1874

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

1876

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

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

1879 '''

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

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

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

1883 else:

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

1885

1886

1887 def _spherical_a_b(a, b):

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

1889 '''

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

1891

1892

1893 def _spherical_f(f):

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

1895 '''

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

1897

1898

1899 def _spherical_f_(f_):

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

1901 '''

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

1903

1904

1905 def a_b2e(a, b):

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

1907

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

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

1910

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

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

1913

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

1915 '''

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

1917

1918

1919 def a_b2e2(a, b):

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

1921

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

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

1924

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

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

1927

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

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

1930 '''

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

1932

1933

1934 def a_b2e22(a, b):

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

1936

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

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

1939

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

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

1942

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

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

1945 '''

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

1947

1948

1949 def a_b2e32(a, b):

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

1951

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

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

1954

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

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

1957

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

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

1960 '''

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

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

1963

1964

1965 def a_b2f(a, b):

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

1967

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

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

1970

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

1972

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

1974 for I{near-spherical} ellipsoids.

1975 '''

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

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

1978

1979

1980 def a_b2f_(a, b):

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

1982

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

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

1985

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

1987

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

1989 for I{near-spherical} ellipsoids.

1990 '''

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

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

1993

1994

1995 def a_b2f2(a, b):

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

1997

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

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

2000

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

2002

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

2004 for I{near-spherical} ellipsoids.

2005 '''

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

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

2008

2009

2010 def a_b2n(a, b):

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

2012

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

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

2015

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

2017

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

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

2020 '''

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

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

2023

2024

2025 def a_f2b(a, f):

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

2027

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

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

2030

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

2032 '''

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

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

2035

2036

2037 def a_f_2b(a, f_):

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

2039

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

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

2042

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

2044 '''

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

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

2047

2048

2049 def b_f2a(b, f):

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

2051

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

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

2054

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

2056 '''

2057 t = _1_0 - f

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

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

2060

2061

2062 def b_f_2a(b, f_):

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

2064

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

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

2067

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

2069 '''

2070 t = f_ - _1_0

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

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

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

2074

2075

2076 def e2f(e):

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

2078

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

2080

2081 @return: The flattening (C{scalar} or C{0}).

2082

2083 @see: Function L{e22f}.

2084 '''

2085 return e22f(e**2)

2086

2087

2088 def e22f(e2):

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

2090

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

2092

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

2094 '''

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

2096

2097

2098 def f2e2(f):

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

2100

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

2102

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

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

2105

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

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

2108

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

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

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

2112 '''

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

2114

2115

2116 def f2e22(f):

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

2118

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

2120

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

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

2123

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

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

2126

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

2128 html/classGeographicLib_1_1Ellipsoid.html>}.

2129 '''

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

2131 t = (_1_0 - f)**2

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

2133

2134

2135 def f2e32(f):

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

2137

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

2139

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

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

2142

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

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

2145

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

2147 html/classGeographicLib_1_1Ellipsoid.html>}.

2148 '''

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

2150 e2 = f2e2(f)

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

2152

2153

2154 def f_2f(f_):

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

2156

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

2158

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

2160

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

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

2163 '''

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

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

2166

2167

2168 def f2f_(f):

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

2170

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

2172

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

2174

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

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

2177 '''

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

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

2180

2181

2182 def f2f2(f):

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

2184

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

2186

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

2188

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

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

2191

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

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

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

2195 '''

2196 t = _1_0 - f

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

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

2199

2200

2201 def f2n(f):

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

2203

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

2205

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

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

2208

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

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

2211

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

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

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

2215 '''

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

2217

2218

2219 def n2e2(n):

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

2221

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

2223

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

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

2226

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

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

2229

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

2231 '''

2232 t = (n + _1_0)**2

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

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

2235

2236

2237 def n2f(n):

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

2239

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

2241

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

2243

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

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

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

2247 '''

2248 t = n + _1_0

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

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

2251

2252

2253 def n2f_(n):

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

2255

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

2257

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

2259

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

2261 '''

2262 return f2f_(n2f(n))

2263

2264

2265 def _normalTo3(px, py, E, eps=EPS): # in .height4 above

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

2267

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

2269 '''

2270 a, b, e0 = E.a, E.b, EPS0

2271 if min(px, py, a, b) < e0:

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

2273

2274 a2 = a - b * E.b_a

2275 b2 = b - a * E.a_b

2276 tx = ty = _SQRT2_2

2277 _a, _h = fabs, hypot

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

2279 ex = a2 * tx**3

2280 ey = b2 * ty**3

2281

2282 qx = px - ex

2283 qy = py - ey

2284 q = _h(qx, qy)

2285 if q < e0: # PYCHOK no cover

2286 break

2287 r = _h(ex - tx * a,

2288 ey - ty * b) / q

2289

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

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

2292 t = _h(ty, tx)

2293 if t < e0: # PYCHOK no cover

2294 break

2295 tx = tx / t # /= chokes PyChecker

2296 ty = ty / t

2297 if _a(sx - tx) < eps and \

2298 _a(sy - ty) < eps:

2299 break

2300

2301 tx *= a / px

2302 ty *= b / py

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

2304

2305

2306 class Ellipsoids(_NamedEnum):

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

2308 to accommodate the L{_LazyNamedEnumItem} properties.

2309 '''

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

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

2312 '''

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

2314

2315 Ellipsoids = Ellipsoids(Ellipsoid) # PYCHOK singleton

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2394# U{NOAA/NOS/NGS/inverse<https://GitHub.com/noaa-ngs/inverse/blob/main/invers3d.f>}

2395 WGS84_NGS = _lazy('WGS84_NGS', *_T(R_MA, _0_0, 298.257222100882711243162836600094))

2396)

2397

2398_EWGS84 = Ellipsoids.WGS84 # (INTERNAL) shared

2399

2400 if __name__ == '__main__':

2401

2402 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_

2403 from pygeodesy import nameof, printf

2404

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

2406 Ellipsoids.Sphere, Ellipsoids.SpherePopular,

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

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

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

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

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

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

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

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

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

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

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

2418

2419 # __doc__ of this file, force all into registry

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

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

2422

2423# % python3 -m pygeodesy.ellipsoids

2424

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

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

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

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

2429# KsOrder 8

2430

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

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

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

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

2435# KsOrder 8

2436

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

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

2439# AlphaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2440# BetaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2441# KsOrder 8

2442

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

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

2445# AlphaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2446# BetaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2447# KsOrder 8

2448

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

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

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

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

2453# KsOrder 8

2454

2455# **) MIT License

2456#

2457# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.

2458#

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

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

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

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

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

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

2465#

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

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

2468#

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

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

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

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

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

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

2475# OTHER DEALINGS IN THE SOFTWARE.