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

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

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

65

66 from pygeodesy.basics import copysign0, isbool, isint

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

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

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

70 _0_25, _3_0 # PYCHOK used!

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

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

73 hypot1, hypot2, sqrt3

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

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

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

77 _Intl1924_, _incompatible_, _invalid_, _Krassovski1940_, \

78 _Krassowsky1940_, _meridional_, _lat_, _negative_, _not_, \

79 _not_finite_, _prime_vertical_, _radius_, _Sphere_, _SPACE_, \

80 _vs_, _WGS72_, _WGS84_

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

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

83 _NamedTuple, _Pass, _ALL_LAZY, _MODS

84 from pygeodesy.namedTuples import Distance2Tuple, Vector3Tuple, Vector4Tuple

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

86 deprecated_property_RO, property_RO

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

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

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

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

91

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

93

94__all__ = _ALL_LAZY.ellipsoids

95__version__ = '24.02.27'

96

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

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

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

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

101

102

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

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

105 '''

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

107

108

109 def _s2_c2(phi):

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

111 '''

112 if phi:

113 s2 = sin(phi)**2

114 if s2 > EPS:

115 c2 = _1_0 - s2

116 if c2 > EPS:

117 if c2 < EPS1:

118 return s2, c2

119 else:

120 return _1_0, _0_0 # phi == PI_2

121 return _0_0, _1_0 # phi == 0

122

123

124 class a_f2Tuple(_NamedTuple):

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

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

127

128 @see: Class L{Ellipsoid2}.

129 '''

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

131 _Units_ = (_Pass, _Pass)

132

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

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

135

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

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

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

139

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

141

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

143

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

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

146 '''

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

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

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

150 f = _f_0_0

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

152

153 @Property_RO

154 def b(self):

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

156 '''

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

158

159 def ellipsoid(self, name=NN):

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

161

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

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

164 '''

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

166

167 @Property_RO

168 def f_(self):

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

170 '''

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

172

173

174 class Circle4Tuple(_NamedTuple):

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

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

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

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

179

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

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

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

183 equals C{lat}.

184 '''

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

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

187

188

189 class Curvature2Tuple(_NamedTuple):

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

191 C{meter}, conventionally.

192 '''

193 _Names_ = (_meridional_, _prime_vertical_)

194 _Units_ = ( Meter, Meter)

195

196 @property_RO

197 def transverse(self):

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

199 '''

200 return self.prime_vertical

201

202

203 class Ellipsoid(_NamedEnumItem):

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

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

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

207 '''

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

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

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

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

212

213 _geodsolve = NN # means, use PYGEODESY_GEODSOLVE

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

215 _rhumbsolve = NN # means, use PYGEODESY_RHUMBSOLVE

216

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

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

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

220

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

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

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

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

225 spherical).

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

227

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

229

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

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

232

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

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

235 '''

236 ff_ = f, f_ # assertion below

237 try:

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

239 if not _isfinite(a):

240 raise ValueError(_SPACE_(_a_, _not_finite_))

241

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

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

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

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

246 elif f is not None:

247 f = Float(f=f)

248 b = a_f2b(a, f)

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

250 elif f_:

251 f_ = Float(f_=f_)

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

253 f = f_2f(f_)

254 else: # only a, spherical

255 f_ = f = 0

256 b = a # superfluous

257

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

259 raise ValueError(_SPACE_(_f_, _invalid_))

260 if not _isfinite(b):

261 raise ValueError(_SPACE_(_b_, _not_finite_))

262

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

264 b = a

265 f = _f_0_0

266 f_ = _f__0_0

267

268 except (TypeError, ValueError) as x:

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

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

271 raise _ValueError(t, cause=x)

272

273 self._a = a

274 self._b = b

275 self._f = f

276 self._f_ = f_

277

278 self._register(Ellipsoids, name)

279

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

281 d = dict(eps=_TOL)

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

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

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

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

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

287

288 def __eq__(self, other):

289 '''Compare this and an other ellipsoid.

290

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

292

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

294 '''

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

296 self.a == other.a and

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

298

299 def __hash__(self):

300 return self._hash # memoized

301

302 @Property_RO

303 def a(self):

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

305 '''

306 return self._a

307

308 equatoradius = a # = Requatorial

309

310 @Property_RO

311 def a2(self):

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

313 '''

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

315

316 @Property_RO

317 def a2_(self):

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

319 '''

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

321

322 @Property_RO

323 def a_b(self):

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

325 '''

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

327

328 @Property_RO

329 def a2_b(self):

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

331

332 @see: U{Radii of Curvature

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

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

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

336

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

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

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

340 '''

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

342

343 @Property_RO

344 def a2_b2(self):

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

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

347 '''

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

349

350 @Property_RO

351 def a_f(self):

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

353 '''

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

355

356 @Property_RO

357 def A(self):

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

359

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

361 '''

362 A, n = self.a, self.n

363 if n:

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

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

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

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

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

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

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

371 return A

372

373 @Property_RO

374 def _albersCyl(self):

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

376 '''

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

378

379 @Property_RO

380 def AlphaKs(self):

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

382 '''

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

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

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

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

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

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

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

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

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

392 _T(1424729850961/743921418240)) # PYCHOK unaligned

393

394 @Property_RO

395 def area(self):

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

397

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

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

400 '''

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

402

403 @Property_RO

404 def areax(self):

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

406 accurate for very I{oblate} ellipsoids.

407

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

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

410 '''

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

412

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

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

415 '''

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

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

418 t = (v, _vs_, val)

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

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

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

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

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

424 eps=eps, f0=f0, **name_value))

425

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

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

428

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

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

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

432

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

434

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

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

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

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

439 '''

440 if self.f:

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

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

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

444

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

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

447

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

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

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

451

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

453

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

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

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

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

458 '''

459 if self.f:

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

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

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

463

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

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

466

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

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

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

470

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

472

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

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

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

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

477 '''

478 if self.f:

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

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

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

482

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

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

485

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

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

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

489

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

491

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

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

494 thereof is the original geodetic latitude.

495

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

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

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

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

500 '''

501 if self.f:

502 r = Phi_(lat, clip=0)

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

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

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

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

507

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

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

510

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

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

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

514

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

516

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

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

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

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

521 '''

522 if self.f:

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

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

525

526 auxReduced = auxParametric # synonymous

527

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

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

530

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

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

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

534

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

536

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

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

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

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

541 '''

542 if self.f:

543 lat = Lat(lat)

544 if 0 < fabs(lat) < _90_0:

545 if inverse:

546 e = self._elliptic_e22

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

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

549 else:

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

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

552

553 @Property_RO

554 def b(self):

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

556 '''

557 return self._b

558

559 polaradius = b # = Rpolar

560

561 @Property_RO

562 def b_a(self):

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

564

565 @see: Property L{f1}.

566 '''

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

568

569 @Property_RO

570 def b2(self):

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

572 '''

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

574

575 @Property_RO

576 def b2_a(self):

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

578

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

580 '''

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

582

583 @Property_RO

584 def b2_a2(self):

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

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

587 '''

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

589

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

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

592 '''

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

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

595 return atan1d(s, c)

596

597 @Property_RO

598 def BetaKs(self):

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

600 '''

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

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

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

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

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

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

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

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

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

610 _T(191773887257/3719607091200)) # PYCHOK unaligned

611

612 @deprecated_Property_RO

613 def c(self): # PYCHOK no cover

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

615 return self.R2

616

617 @Property_RO

618 def c2(self):

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

620

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

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

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

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

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

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

627 '''

628 return self._c2f(False)

629

630 @Property_RO

631 def c2x(self):

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

633 ellipsoids.

634

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

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

637 '''

638 return self._c2f(True)

639

640 def _c2f(self, c2x):

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

642 '''

643 f, c2 = self.f, self.b2

644 if f:

645 e = self.e

646 if e > EPS0:

647 if f > 0: # .isOblate

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

649 elif f < 0: # .isProlate

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

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

652 return c2

653

654 def circle4(self, lat):

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

656

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

658

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

660 instance.

661

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

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

664

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

666

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

668

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

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

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

672 and method C{Rlat}.

673 '''

674 lat = Lat(lat)

675 if lat:

676 b = lat

677 if fabs(lat) < _90_0:

678 if self.f:

679 b = self._beta(lat)

680 z, r = sincos2d(b)

681 r *= self.a

682 z *= self.b

683 else: # near-polar

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

685 else: # equator

686 r = self.a

687 z = lat = b = _0_0

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

689

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

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

692 along a parallel of (geodetic) latitude.

693

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

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

696

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

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

699

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

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

702

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

704 '''

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

706

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

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

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

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

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

712

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

714 and only between points not near-polar}.

715

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

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

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

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

720

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

722 in same units as this ellipsoid's axes.

723

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

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

726

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

728 '''

729 phi0 = Phi_(lat0=lat0)

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

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

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

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

734

735 @Property_RO

736 def e(self):

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

738

739 @see: Property L{es}.

740 '''

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

742

743 @deprecated_Property_RO

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

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

746 return self.e21

747

748# @Property_RO

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

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

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

752

753 @Property_RO

754 def e2(self):

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

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

757 '''

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

759

760 @Property_RO

761 def e2abs(self):

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

763 '''

764 return fabs(self.e2)

765

766 @Property_RO

767 def e21(self):

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

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

770 '''

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

772

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

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

775

776 @Property_RO

777 def e22(self):

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

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

780 '''

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

782

783 @Property_RO

784 def e22abs(self):

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

786 '''

787 return fabs(self.e22)

788

789 @Property_RO

790 def e32(self):

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

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

793 '''

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

795

796 @Property_RO

797 def e32abs(self):

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

799 '''

800 return fabs(self.e32)

801

802 @Property_RO

803 def e4(self):

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

805 '''

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

807

808 eccentricity = e # eccentricity

809# eccentricity2 = e2 # eccentricity squared

810 eccentricity1st2 = e2 # first eccentricity squared

811 eccentricity2nd2 = e22 # second eccentricity squared

812 eccentricity3rd2 = e32 # third eccentricity squared

813

814 def ecef(self, Ecef=None):

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

816

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

818

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

820

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

822

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

824 '''

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

826

827 @Property_RO

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

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

830 '''

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

832 return _MODS.elliptic.Elliptic(e12)

833

834 @Property_RO

835 def _elliptic_e22(self): # aka ._elliptic_ep2

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

837 '''

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

839

840 equatoradius = a # Requatorial

841

842 def e2s(self, s):

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

844

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

846

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

848

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

850 '''

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

852

853 def e2s2(self, s):

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

855

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

857

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

859

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

861 '''

862 r = _1_0

863 if self.e2:

864 try:

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

866 if r < 0:

867 raise ValueError(_negative_)

868 except (TypeError, ValueError) as x:

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

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

871 return r

872

873 @Property_RO

874 def es(self):

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

876

877 @see: Property L{e}.

878 '''

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

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

881

882 def es_atanh(self, x):

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

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

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

886

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

888

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

890 html/classGeographicLib_1_1Math.html>}.

891 '''

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

893

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

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

896 '''

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

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

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

900 _0_0)) # .isSpherical

901

902 @Property_RO

903 def es_c(self):

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

905 '''

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

907

908 def es_tauf(self, taup):

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

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

911

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

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

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

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

916 '''

917 t = Scalar(taup=taup)

918 if self.f: # .isEllipsoidal

919 a = fabs(t)

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

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

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

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

924 if fabs(d) < s:

925 break

926 t = Scalar(tauf=T)

927 return t

928

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

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

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

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

933 '''

934 e = self._1_e21

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

936 _tf2 = self._es_taupf2

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

938 a, h = _tf2(T)

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

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

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

942 yield T, i, d

943

944 def es_taupf(self, tau):

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

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

947

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

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

950 '''

951 t = Scalar(tau=tau)

952 if self.f: # .isEllipsoidal

953 t, _ = self._es_taupf2(t)

954 t = Scalar(taupf=t)

955 return t

956

957 def _es_taupf2(self, tau):

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

959 '''

960 if _isfinite(tau):

961 h = hypot1(tau)

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

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

964 else:

965 a, h = tau, INF

966 return a, h

967

968 @Property_RO

969 def _exp_es_atanh_1(self):

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

971 '''

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

973

974 @Property_RO

975 def f(self):

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

977 '''

978 return self._f

979

980 @Property_RO

981 def f_(self):

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

983 '''

984 return self._f_

985

986 @Property_RO

987 def f1(self):

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

989

990 @see: Property L{b_a}.

991 '''

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

993

994 @Property_RO

995 def f2(self):

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

997 '''

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

999

1000 @deprecated_Property_RO

1001 def geodesic(self):

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

1003 return self.geodesicw

1004

1005 def geodesic_(self, exact=True):

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

1007

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

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

1010 instance for I{this} ellipsoid.

1011

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

1013

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

1015

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

1017 '''

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

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

1020 name=self.name)

1021 else:

1022 g = exact

1023 E = _xattr(g, ellipsoid=None)

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

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

1026 return g

1027

1028 @property_RO

1029 def _Geodesics(self):

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

1031 '''

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

1033 _MODS.geodesicx.GeodesicExact,

1034 _MODS.geodsolve.GeodesicSolve)

1035 return t

1036

1037 @property_RO

1038 def geodesicw(self):

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

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

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

1042 package is installed.

1043 '''

1044 # if not self.isEllipsoidal:

1045 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

1046 return _MODS.geodesicw.Geodesic(self)

1047

1048 @property_RO

1049 def geodesicx(self):

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

1051 '''

1052 # if not self.isEllipsoidal:

1053 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1055

1056 @property

1057 def geodsolve(self):

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

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

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

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

1062 '''

1063 # if not self.isEllipsoidal:

1064 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1066

1067 @geodsolve.setter # PYCHOK setter!

1068 def geodsolve(self, path):

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

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

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

1072 '''

1073 self._geodsolve = path

1074

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

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

1077 from a Point-Of-View in space.

1078

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

1080 or L{Vector3d}).

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

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

1083

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

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

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

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

1088

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

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

1091 outside this ellipsoid or points in an opposite

1092 direction.

1093

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

1095

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

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

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

1099 '''

1100 try:

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

1102 except Exception as x:

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

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

1105

1106 @Property_RO

1107 def _hash(self):

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

1109

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

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

1112 this ellipsoid's surface.

1113

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

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

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

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

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

1119

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

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

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

1123

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

1125

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

1127

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

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

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

1131 '''

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

1133 r = v.length

1134

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

1136 if r < EPS0: # EPS

1137 v = v.times(_0_0)

1138 h = -a

1139

1140 elif self.isSpherical:

1141 v = v.times(a / r)

1142 h = r - a

1143

1144 elif normal: # perpendicular to ellipsoid

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

1146 if x < EPS0: # PYCHOK no cover

1147 z = copysign0(b, v.z)

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

1149 h = y - b # polar

1150 elif y < EPS0: # PYCHOK no cover

1151 t = a / r

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

1153 h = x - a # equatorial

1154 else: # normal in 1st quadrant

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

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

1157 h = t.minus(v).length

1158

1159 else: # radial to ellipsoid's center

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

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

1162 v = v.times(t)

1163 h = r * (_1_0 - t)

1164

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

1166 name=self.height4.__name__)

1167

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

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

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

1171 '''

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

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

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

1175

1176 @Property_RO

1177 def isEllipsoidal(self):

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

1179 '''

1180 return self.f != 0

1181

1182 @Property_RO

1183 def isOblate(self):

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

1185 spherical otherwise.

1186 '''

1187 return self.f > 0

1188

1189 @Property_RO

1190 def isProlate(self):

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

1192 spherical otherwise.

1193 '''

1194 return self.f < 0

1195

1196 @Property_RO

1197 def isSpherical(self):

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

1199 '''

1200 return self.f == 0

1201

1202 def _Kseries(self, *AB8Ks):

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

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

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

1206

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

1208 coefficient tuples.

1209

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

1211

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

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

1214 '''

1215 k = self.KsOrder

1216 if self.n:

1217 ns = fpowers(self.n, k)

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

1219 else:

1220 ks = _0_0s(k)

1221 return ks

1222

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

1224 def KsOrder(self):

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

1226 '''

1227 return self._KsOrder

1228

1229 @KsOrder.setter # PYCHOK setter!

1230 def KsOrder(self, order):

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

1232

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

1234 '''

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

1236 raise _ValueError(order=order)

1237 if self._KsOrder != order:

1238 Ellipsoid.AlphaKs._update(self)

1239 Ellipsoid.BetaKs._update(self)

1240 self._KsOrder = order

1241

1242 @Property_RO

1243 def L(self):

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

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

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

1247 '''

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

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

1250

1251 def Llat(self, lat):

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

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

1254

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

1256

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

1258 hemisphere (C{meter}).

1259 '''

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

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

1262

1263 Lmeridian = Llat # meridional distance

1264

1265 @property_RO

1266 def _Lpd(self):

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

1268 '''

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

1270

1271 @property_RO

1272 def _Lpr(self):

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

1274 '''

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

1276

1277 @deprecated_Property_RO

1278 def majoradius(self): # PYCHOK no cover

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

1280 return self.a

1281

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

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

1284 along a parallel of (geodetic) latitude.

1285

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

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

1288

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

1290

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

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

1293

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

1295 '''

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

1297

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

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

1300 along a parallel of (geodetic) latitude.

1301

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

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

1304

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

1306

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

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

1309

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

1311 '''

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

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

1314

1315 @deprecated_Property_RO

1316 def minoradius(self): # PYCHOK no cover

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

1318 return self.b

1319

1320 @Property_RO

1321 def n(self):

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

1323 '''

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

1325

1326 flattening = f

1327 flattening1st = f

1328 flattening2nd = f2

1329 flattening3rd = n

1330

1331 polaradius = b # Rpolar

1332

1333# @Property_RO

1334# def Q(self):

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

1336#

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

1338#

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

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

1341# '''

1342# n = self.n

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

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

1345

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

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

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

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

1350

1351 @deprecated_Property_RO

1352 def quarteradius(self): # PYCHOK no cover

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

1354 return self.L

1355

1356 @Property_RO

1357 def R1(self):

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

1359

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

1361 and method C{Rgeometric}.

1362 '''

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

1364 return Radius(R1=r)

1365

1366 Rmean = R1

1367

1368 @Property_RO

1369 def R2(self):

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

1371

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

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

1374 '''

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

1376

1377 Rauthalic = R2

1378

1379# @Property_RO

1380# def R2(self):

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

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

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

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

1385

1386 @Property_RO

1387 def R2x(self):

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

1389

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

1391 '''

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

1393

1394 Rauthalicx = R2x

1395

1396 @Property_RO

1397 def R3(self):

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

1399

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

1401 '''

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

1403 return Radius(R3=r)

1404

1405 Rvolumetric = R3

1406

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

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

1409 along a parallel of (geodetic) latitude.

1410

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

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

1413

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

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

1416

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

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

1419

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

1421 '''

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

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

1424

1425 @Property_RO

1426 def Rbiaxial(self):

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

1428

1429 @see: C{Rtriaxial}

1430 '''

1431 a, b = self.a, self.b

1432 if b < a:

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

1434 elif b > a:

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

1436 return Radius(Rbiaxial=b)

1437

1438 Requatorial = a # for consistent naming

1439

1440 def Rgeocentric(self, lat):

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

1442

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

1444

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

1446

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

1448

1449 @see: U{Geocentric Radius

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

1451 '''

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

1453 if a and self.f:

1454 if fabs(a) < PI_2:

1455 s2, c2 = _s2_c2(a)

1456 b2_a2_s2 = self.b2_a2 * s2

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

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

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

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

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

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

1463 else:

1464 r = self.b

1465 return Radius(Rgeocentric=r)

1466

1467 @Property_RO

1468 def Rgeometric(self):

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

1470

1471 @see: C{R1}.

1472 '''

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

1474 return Radius(Rgeometric=g)

1475

1476 def rhumb_(self, exact=True):

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

1478

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

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

1481 for I{this} ellipsoid.

1482

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

1484

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

1486

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

1488 '''

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

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

1491 else:

1492 r = exact

1493 E = _xattr(r, ellipsoid=None)

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

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

1496 return r

1497

1498 @property_RO

1499 def rhumbaux(self):

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

1501 '''

1502 # if not self.isEllipsoidal:

1503 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1505

1506 @property_RO

1507 def rhumbekx(self):

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

1509 '''

1510 # if not self.isEllipsoidal:

1511 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1513

1514 @property_RO

1515 def _Rhumbs(self):

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

1517 '''

1518 p = _MODS.rhumb

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

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

1521 return t

1522

1523 @property

1524 def rhumbsolve(self):

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

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

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

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

1529 '''

1530 # if not self.isEllipsoidal:

1531 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1533

1534 @rhumbsolve.setter # PYCHOK setter!

1535 def rhumbsolve(self, path):

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

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

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

1539 '''

1540 self._rhumbsolve = path

1541

1542 @deprecated_property_RO

1543 def rhumbx(self):

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

1545 return self.rhumbekx

1546

1547 def Rlat(self, lat):

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

1549

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

1551

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

1553

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

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

1556

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

1558

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

1560

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

1562

1563 @see: Method C{circle4}.

1564 '''

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

1566 r = self.a

1567 if self.f and lat: # .isEllipsoidal

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

1569 r = Radius(Rlat=r)

1570 return r

1571

1572 Rpolar = b # for consistent naming

1573

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

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

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

1577

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

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

1580 to use an alternate formula.

1581

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

1583

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

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

1586

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

1588 '''

1589 if not self.f: # .isSpherical

1590 n = self.a

1591 elif ca is None:

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

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

1594 elif ca: # derived from EcefYou.forward

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

1596 n = self.a / h

1597 elif sa:

1598 n = self.a2_b / fabs(sa)

1599 else:

1600 n = self.a

1601 return n

1602

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

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

1605 radii of curvature of (geodetic) latitude.

1606

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

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

1609

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

1611 the radii of curvature.

1612

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

1614

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

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

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

1618 Earth_radius#Radii_of_curvature>}.

1619 '''

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

1621

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

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

1624 curvature of (geodetic) latitude.

1625

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

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

1628

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

1630 radii of curvature.

1631

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

1633

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

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

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

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

1638 '''

1639 a = fabs(Phi(phi))

1640 if self.f:

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

1642 if r > EPS02:

1643 n = self.a / sqrt(r)

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

1645 else:

1646 m = n = _0_0 # PYCHOK attr

1647 else:

1648 m = n = self.a

1649 if scaled and a:

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

1651 return Curvature2Tuple(Radius(rocMeridional=m),

1652 Radius(rocPrimeVertical=n))

1653

1654 def rocBearing(self, lat, bearing):

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

1656 latitude and compass direction.

1657

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

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

1660

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

1662

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

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

1665

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

1667

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

1669 '''

1670 if self.f:

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

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

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

1674 c2 *= n / m

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

1676 s2 *= m / n

1677 n = m

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

1679 else:

1680 b = self.b # == self.a

1681 return Radius(rocBearing=b)

1682

1683 @Property_RO

1684 def rocEquatorial2(self):

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

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

1687

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

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

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

1691 '''

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

1693 Radius(rocPrimeVertical0=self.a))

1694

1695 def rocGauss(self, lat):

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

1697

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

1699

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

1701

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

1703

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

1705 Earth_radius#Radii_of_curvature>}

1706 '''

1707 # using ...

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

1709 # return sqrt(m * n)

1710 # ... requires 1 or 2 sqrt

1711 g = self.b

1712 if self.f:

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

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

1715 return Radius(rocGauss=g)

1716

1717 def rocMean(self, lat):

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

1719

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

1721

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

1723

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

1725

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

1727 Earth_radius#Radii_of_curvature>}

1728 '''

1729 if self.f:

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

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

1732 else:

1733 m = self.a

1734 return Radius(rocMean=m)

1735

1736 def rocMeridional(self, lat):

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

1738

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

1740

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

1742

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

1744

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

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

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

1748 '''

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

1750 self.rocEquatorial2.meridional

1751

1752 rocPolar = a2_b # synonymous

1753

1754 def rocPrimeVertical(self, lat):

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

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

1757

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

1759

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

1761

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

1763

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

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

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

1767 Earth_radius#Radii_of_curvature>}.

1768 '''

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

1770 self.rocEquatorial2.prime_vertical

1771

1772 rocTransverse = rocPrimeVertical # synonymous

1773

1774 @deprecated_Property_RO

1775 def Rquadratic(self): # PYCHOK no cover

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

1777 return self.Rbiaxial

1778

1779 @deprecated_Property_RO

1780 def Rr(self): # PYCHOK no cover

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

1782 return self.Rrectifying

1783

1784 @Property_RO

1785 def Rrectifying(self):

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

1787

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

1789 '''

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

1791 return Radius(Rrectifying=r)

1792

1793 @deprecated_Property_RO

1794 def Rs(self): # PYCHOK no cover

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

1796 return self.Rgeometric

1797

1798 @Property_RO

1799 def Rtriaxial(self):

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

1801

1802 @see: C{Rbiaxial}

1803 '''

1804 a, b = self.a, self.b

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

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

1807 return Radius(Rtriaxial=q)

1808

1809 def toEllipsoid2(self, name=NN):

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

1811

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

1813

1814 @see: Property C{a_f}.

1815 '''

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

1817

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

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

1820

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

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

1823 this ellipsoid's name.

1824

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

1826 '''

1827 E = Ellipsoid

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

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

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

1831 t = t[:4]

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

1833

1834 def toTriaxial(self, name=NN):

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

1836

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

1838 pointing east and C{Z} pointing north.

1839

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

1841 '''

1842 T = self._triaxial

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

1844

1845 @property_RO

1846 def _triaxial(self):

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

1848 '''

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

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

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

1852

1853 @Property_RO

1854 def volume(self):

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

1856

1857 @see: C{R3}.

1858 '''

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

1860

1861

1862 class Ellipsoid2(Ellipsoid):

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

1864 '''

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

1866 '''New L{Ellipsoid2}.

1867

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

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

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

1871

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

1873

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

1875

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

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

1878 '''

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

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

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

1882 else:

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

1884

1885

1886 def _spherical_a_b(a, b):

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

1888 '''

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

1890

1891

1892 def _spherical_f(f):

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

1894 '''

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

1896

1897

1898 def _spherical_f_(f_):

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

1900 '''

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

1902

1903

1904 def a_b2e(a, b):

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

1906

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

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

1909

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

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

1912

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

1914 '''

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

1916

1917

1918 def a_b2e2(a, b):

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

1920

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

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

1923

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

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

1926

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

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

1929 '''

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

1931

1932

1933 def a_b2e22(a, b):

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

1935

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

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

1938

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

1940 M{(a / b)**2 - 1}.

1941

1942 @note: The result is positive for I{oblate}, negative for I{prolate}

1943 or C{0} for I{near-spherical} ellipsoids.

1944 '''

1945 return Float(e22=_0_0 if _spherical_a_b(a, b) else ((a - b) * (a + b) / b**2))

1946

1947

1948 def a_b2e32(a, b):

1949 '''Return C{e32}, the I{3rd eccentricity squared} for a given I{equatorial} and I{polar} radius.

1950

1951 @arg a: Equatorial radius (C{scalar} > 0).

1952 @arg b: Polar radius (C{scalar} > 0).

1953

1954 @return: The I{signed}, 3rd eccentricity I{squared} (C{float} or C{0}),

1955 M{(a**2 - b**2) / (a**2 + b**2)}.

1956

1957 @note: The result is positive for I{oblate}, negative for I{prolate}

1958 or C{0} for I{near-spherical} ellipsoids.

1959 '''

1960 a2, b2 = a**2, b**2

1961 return Float(e32=_0_0 if _spherical_a_b(a2, b2) else ((a2 - b2) / (a2 + b2)))

1962

1963

1964 def a_b2f(a, b):

1965 '''Return C{f}, the I{flattening} for a given I{equatorial} and I{polar} radius.

1966

1967 @arg a: Equatorial radius (C{scalar} > 0).

1968 @arg b: Polar radius (C{scalar} > 0).

1969

1970 @return: The flattening (C{scalar} or C{0}), M{(a - b) / a}.

1971

1972 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0}

1973 for I{near-spherical} ellipsoids.

1974 '''

1975 f = 0 if _spherical_a_b(a, b) else ((a - b) / a)

1976 return _f_0_0 if _spherical_f(f) else Float(f=f)

1977

1978

1979 def a_b2f_(a, b):

1980 '''Return C{f_}, the I{inverse flattening} for a given I{equatorial} and I{polar} radius.

1981

1982 @arg a: Equatorial radius (C{scalar} > 0).

1983 @arg b: Polar radius (C{scalar} > 0).

1984

1985 @return: The inverse flattening (C{scalar} or C{0}), M{a / (a - b)}.

1986

1987 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0}

1988 for I{near-spherical} ellipsoids.

1989 '''

1990 f_ = 0 if _spherical_a_b(a, b) else (a / float(a - b))

1991 return _f__0_0 if _spherical_f_(f_) else Float(f_=f_)

1992

1993

1994 def a_b2f2(a, b):

1995 '''Return C{f2}, the I{2nd flattening} for a given I{equatorial} and I{polar} radius.

1996

1997 @arg a: Equatorial radius (C{scalar} > 0).

1998 @arg b: Polar radius (C{scalar} > 0).

1999

2000 @return: The I{signed}, 2nd flattening (C{scalar} or C{0}), M{(a - b) / b}.

2001

2002 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0}

2003 for I{near-spherical} ellipsoids.

2004 '''

2005 t = 0 if _spherical_a_b(a, b) else float(a - b)

2006 return Float(f2=_0_0 if fabs(t) < EPS0 else (t / b))

2007

2008

2009 def a_b2n(a, b):

2010 '''Return C{n}, the I{3rd flattening} for a given I{equatorial} and I{polar} radius.

2011

2012 @arg a: Equatorial radius (C{scalar} > 0).

2013 @arg b: Polar radius (C{scalar} > 0).

2014

2015 @return: The I{signed}, 3rd flattening (C{scalar} or C{0}), M{(a - b) / (a + b)}.

2016

2017 @note: The result is positive for I{oblate}, negative for I{prolate}

2018 or C{0} for I{near-spherical} ellipsoids.

2019 '''

2020 t = 0 if _spherical_a_b(a, b) else float(a - b)

2021 return Float(n=_0_0 if fabs(t) < EPS0 else (t / (a + b)))

2022

2023

2024 def a_f2b(a, f):

2025 '''Return C{b}, the I{polar} radius for a given I{equatorial} radius and I{flattening}.

2026

2027 @arg a: Equatorial radius (C{scalar} > 0).

2028 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2029

2030 @return: The polar radius (C{float}), M{a * (1 - f)}.

2031 '''

2032 b = a if _spherical_f(f) else (a * (_1_0 - f))

2033 return Radius_(b=a if _spherical_a_b(a, b) else b)

2034

2035

2036 def a_f_2b(a, f_):

2037 '''Return C{b}, the I{polar} radius for a given I{equatorial} radius and I{inverse flattening}.

2038

2039 @arg a: Equatorial radius (C{scalar} > 0).

2040 @arg f_: Inverse flattening (C{scalar} >>> 1).

2041

2042 @return: The polar radius (C{float}), M{a * (f_ - 1) / f_}.

2043 '''

2044 b = a if _spherical_f_(f_) else (a * (f_ - _1_0) / f_)

2045 return Radius_(b=a if _spherical_a_b(a, b) else b)

2046

2047

2048 def b_f2a(b, f):

2049 '''Return C{a}, the I{equatorial} radius for a given I{polar} radius and I{flattening}.

2050

2051 @arg b: Polar radius (C{scalar} > 0).

2052 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2053

2054 @return: The equatorial radius (C{float}), M{b / (1 - f)}.

2055 '''

2056 t = _1_0 - f

2057 a = b if fabs(t) < EPS0 else (b / t)

2058 return Radius_(a=b if _spherical_a_b(a, b) else a)

2059

2060

2061 def b_f_2a(b, f_):

2062 '''Return C{a}, the I{equatorial} radius for a given I{polar} radius and I{inverse flattening}.

2063

2064 @arg b: Polar radius (C{scalar} > 0).

2065 @arg f_: Inverse flattening (C{scalar} >>> 1).

2066

2067 @return: The equatorial radius (C{float}), M{b * f_ / (f_ - 1)}.

2068 '''

2069 t = f_ - _1_0

2070 a = b if _spherical_f_(f_) or fabs(t - f_) < EPS0 \

2071 or fabs(t) < EPS0 else (b * f_ / t)

2072 return Radius_(a=b if _spherical_a_b(a, b) else a)

2073

2074

2075 def e2f(e):

2076 '''Return C{f}, the I{flattening} for a given I{1st eccentricity}.

2077

2078 @arg e: The (1st) eccentricity (0 <= C{float} < 1)

2079

2080 @return: The flattening (C{scalar} or C{0}).

2081

2082 @see: Function L{e22f}.

2083 '''

2084 return e22f(e**2)

2085

2086

2087 def e22f(e2):

2088 '''Return C{f}, the I{flattening} for a given I{1st eccentricity squared}.

2089

2090 @arg e2: The (1st) eccentricity I{squared}, I{signed} (L{NINF} < C{float} < 1)

2091

2092 @return: The flattening (C{float} or C{0}), M{e2 / (sqrt(e2 - 1) + 1)}.

2093 '''

2094 return Float(f=e2 / (sqrt(_1_0 - e2) + _1_0)) if e2 else _f_0_0

2095

2096

2097 def f2e2(f):

2098 '''Return C{e2}, the I{1st eccentricity squared} for a given I{flattening}.

2099

2100 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2101

2102 @return: The I{signed}, (1st) eccentricity I{squared} (C{float} < 1),

2103 M{f * (2 - f)}.

2104

2105 @note: The result is positive for I{oblate}, negative for I{prolate}

2106 or C{0} for I{near-spherical} ellipsoids.

2107

2108 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2109 html/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2110 <https://WikiPedia.org/wiki/Flattening>}.

2111 '''

2112 return Float(e2=_0_0 if _spherical_f(f) else (f * (_2_0 - f)))

2113

2114

2115 def f2e22(f):

2116 '''Return C{e22}, the I{2nd eccentricity squared} for a given I{flattening}.

2117

2118 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2119

2120 @return: The I{signed}, 2nd eccentricity I{squared} (C{float} > -1 or

2121 C{INF}), M{f * (2 - f) / (1 - f)**2}.

2122

2123 @note: The result is positive for I{oblate}, negative for I{prolate}

2124 or C{0} for near-spherical ellipsoids.

2125

2126 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2127 html/classGeographicLib_1_1Ellipsoid.html>}.

2128 '''

2129 # e2 / (1 - e2) == f * (2 - f) / (1 - f)**2

2130 t = (_1_0 - f)**2

2131 return Float(e22=INF if t < EPS0 else (f2e2(f) / t)) # PYCHOK type

2132

2133

2134 def f2e32(f):

2135 '''Return C{e32}, the I{3rd eccentricity squared} for a given I{flattening}.

2136

2137 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2138

2139 @return: The I{signed}, 3rd eccentricity I{squared} (C{float}),

2140 M{f * (2 - f) / (1 + (1 - f)**2)}.

2141

2142 @note: The result is positive for I{oblate}, negative for I{prolate}

2143 or C{0} for I{near-spherical} ellipsoids.

2144

2145 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2146 html/classGeographicLib_1_1Ellipsoid.html>}.

2147 '''

2148 # e2 / (2 - e2) == f * (2 - f) / (1 + (1 - f)**2)

2149 e2 = f2e2(f)

2150 return Float(e32=e2 / (_2_0 - e2))

2151

2152

2153 def f_2f(f_):

2154 '''Return C{f}, the I{flattening} for a given I{inverse flattening}.

2155

2156 @arg f_: Inverse flattening (C{scalar} >>> 1).

2157

2158 @return: The flattening (C{scalar} or C{0}), M{1 / f_}.

2159

2160 @note: The result is positive for I{oblate}, negative for I{prolate}

2161 or C{0} for I{near-spherical} ellipsoids.

2162 '''

2163 f = 0 if _spherical_f_(f_) else _1_0 / f_

2164 return _f_0_0 if _spherical_f(f) else Float(f=f) # PYCHOK type

2165

2166

2167 def f2f_(f):

2168 '''Return C{f_}, the I{inverse flattening} for a given I{flattening}.

2169

2170 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2171

2172 @return: The inverse flattening (C{scalar} or C{0}), M{1 / f}.

2173

2174 @note: The result is positive for I{oblate}, negative for I{prolate}

2175 or C{0} for I{near-spherical} ellipsoids.

2176 '''

2177 f_ = 0 if _spherical_f(f) else _1_0 / f

2178 return _f__0_0 if _spherical_f_(f_) else Float(f_=f_) # PYCHOK type

2179

2180

2181 def f2f2(f):

2182 '''Return C{f2}, the I{2nd flattening} for a given I{flattening}.

2183

2184 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2185

2186 @return: The I{signed}, 2nd flattening (C{scalar} or C{INF}), M{f / (1 - f)}.

2187

2188 @note: The result is positive for I{oblate}, negative for I{prolate}

2189 or C{0} for I{near-spherical} ellipsoids.

2190

2191 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2192 html/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2193 <https://WikiPedia.org/wiki/Flattening>}.

2194 '''

2195 t = _1_0 - f

2196 return Float(f2=_0_0 if _spherical_f(f) else (INF if fabs(t) < EPS

2197 else (f / t))) # PYCHOK type

2198

2199

2200 def f2n(f):

2201 '''Return C{n}, the I{3rd flattening} for a given I{flattening}.

2202

2203 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2204

2205 @return: The I{signed}, 3rd flattening (-1 <= C{float} < 1),

2206 M{f / (2 - f)}.

2207

2208 @note: The result is positive for I{oblate}, negative for I{prolate}

2209 or C{0} for I{near-spherical} ellipsoids.

2210

2211 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2212 html/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2213 <https://WikiPedia.org/wiki/Flattening>}.

2214 '''

2215 return Float(n=_0_0 if _spherical_f(f) else (f / float(_2_0 - f)))

2216

2217

2218 def n2e2(n):

2219 '''Return C{e2}, the I{1st eccentricity squared} for a given I{3rd flattening}.

2220

2221 @arg n: The 3rd flattening (-1 <= C{scalar} < 1).

2222

2223 @return: The I{signed}, (1st) eccentricity I{squared} (C{float} or NINF),

2224 M{4 * n / (1 + n)**2}.

2225

2226 @note: The result is positive for I{oblate}, negative for I{prolate}

2227 or C{0} for I{near-spherical} ellipsoids.

2228

2229 @see: U{Flattening<https://WikiPedia.org/wiki/Flattening>}.

2230 '''

2231 t = (n + _1_0)**2

2232 return Float(e2=_0_0 if fabs(n) < EPS0 else

2233 (NINF if t < EPS0 else (_4_0 * n / t)))

2234

2235

2236 def n2f(n):

2237 '''Return C{f}, the I{flattening} for a given I{3rd flattening}.

2238

2239 @arg n: The 3rd flattening (-1 <= C{scalar} < 1).

2240

2241 @return: The flattening (C{scalar} or NINF), M{2 * n / (1 + n)}.

2242

2243 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2244 html/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2245 <https://WikiPedia.org/wiki/Flattening>}.

2246 '''

2247 t = n + _1_0

2248 f = 0 if fabs(n) < EPS0 else (NINF if t < EPS0 else (_2_0 * n / t))

2249 return _f_0_0 if _spherical_f(f) else Float(f=f)

2250

2251

2252 def n2f_(n):

2253 '''Return C{f_}, the I{inverse flattening} for a given I{3rd flattening}.

2254

2255 @arg n: The 3rd flattening (-1 <= C{scalar} < 1).

2256

2257 @return: The inverse flattening (C{scalar} or C{0}), M{1 / f}.

2258

2259 @see: L{n2f} and L{f2f_}.

2260 '''

2261 return f2f_(n2f(n))

2262

2263

2264 def _normalTo3(px, py, E, eps=EPS): # in .height4 above

2265 '''(INTERNAL) Nearest point on a 2-D ellipse in 1st quadrant.

2266

2267 @see: Functions C{pygeodesy.triaxial._normalTo4} and C{-To5}.

2268 '''

2269 a, b, e0 = E.a, E.b, EPS0

2270 if min(px, py, a, b) < e0:

2271 raise _AssertionError(px=px, py=py, a=a, b=b, E=E)

2272

2273 a2 = a - b * E.b_a

2274 b2 = b - a * E.a_b

2275 tx = ty = _SQRT2_2

2276 _a, _h = fabs, hypot

2277 for i in range(16): # max 5

2278 ex = a2 * tx**3

2279 ey = b2 * ty**3

2280

2281 qx = px - ex

2282 qy = py - ey

2283 q = _h(qx, qy)

2284 if q < e0: # PYCHOK no cover

2285 break

2286 r = _h(ex - tx * a,

2287 ey - ty * b) / q

2288

2289 sx, tx = tx, min(_1_0, max(0, (ex + qx * r) / a))

2290 sy, ty = ty, min(_1_0, max(0, (ey + qy * r) / b))

2291 t = _h(ty, tx)

2292 if t < e0: # PYCHOK no cover

2293 break

2294 tx = tx / t # /= chokes PyChecker

2295 ty = ty / t

2296 if _a(sx - tx) < eps and \

2297 _a(sy - ty) < eps:

2298 break

2299

2300 tx *= a / px

2301 ty *= b / py

2302 return tx, ty, i # x and y as fractions

2303

2304

2305 class Ellipsoids(_NamedEnum):

2306 '''(INTERNAL) L{Ellipsoid} registry, I{must} be a sub-class

2307 to accommodate the L{_LazyNamedEnumItem} properties.

2308 '''

2309 def _Lazy(self, a, b, f_, **kwds):

2310 '''(INTERNAL) Instantiate the L{Ellipsoid}.

2311 '''

2312 return Ellipsoid(a, b=b, f_=f_, **kwds)

2313

2314 Ellipsoids = Ellipsoids(Ellipsoid) # PYCHOK singleton

2315'''Some pre-defined L{Ellipsoid}s, all I{lazily} instantiated.'''

2316# <https://www.GNU.org/software/gama/manual/html_node/Supported-ellipsoids.html>

2317# <https://GSSC.ESA.int/navipedia/index.php/Reference_Frames_in_GNSS>

2318# <https://kb.OSU.edu/dspace/handle/1811/77986>

2319# <https://www.IBM.com/docs/en/db2/11.5?topic=systems-supported-spheroids>

2320# <https://w3.Energistics.org/archive/Epicentre/Epicentre_v3.0/DataModel/LogicalDictionary/StandardValues/ellipsoid.html>

2321# <https://GitHub.com/locationtech/proj4j/blob/master/src/main/java/org/locationtech/proj4j/datum/Ellipsoid.java>

2322 Ellipsoids._assert( # <https://WikiPedia.org/wiki/Earth_ellipsoid>

2323 Airy1830 = _lazy(_Airy1830_, *_T(6377563.396, _0_0, 299.3249646)), # b=6356256.909

2324 AiryModified = _lazy(_AiryModified_, *_T(6377340.189, _0_0, 299.3249646)), # b=6356034.448

2325# APL4_9 = _lazy('APL4_9', *_T(6378137.0, _0_0, 298.24985392)), # Appl. Phys. Lab. 1965

2326# ANS = _lazy('ANS', *_T(6378160.0, _0_0, 298.25)), # Australian Nat. Spheroid

2327# AN_SA96 = _lazy('AN_SA96', *_T(6378160.0, _0_0, 298.24985392)), # Australian Nat. South America

2328 Australia1966 = _lazy('Australia1966', *_T(6378160.0, _0_0, 298.25)), # b=6356774.7192

2329 ATS1977 = _lazy('ATS1977', *_T(6378135.0, _0_0, 298.257)), # "Average Terrestrial System"

2330 Bessel1841 = _lazy(_Bessel1841_, *_T(6377397.155, 6356078.962818, 299.152812797)),

2331 BesselModified = _lazy('BesselModified', *_T(6377492.018, _0_0, 299.1528128)),

2332# BesselNamibia = _lazy('BesselNamibia', *_T(6377483.865, _0_0, 299.1528128)),

2333 CGCS2000 = _lazy('CGCS2000', *_T(R_MA, _0_0, 298.257222101)), # BeiDou Coord System (BDC)

2334# Clarke1858 = _lazy('Clarke1858', *_T(6378293.639, _0_0, 294.260676369)),

2335 Clarke1866 = _lazy(_Clarke1866_, *_T(6378206.4, 6356583.8, 294.978698214)),

2336 Clarke1880 = _lazy('Clarke1880', *_T(6378249.145, 6356514.86954978, 293.465)),

2337 Clarke1880IGN = _lazy(_Clarke1880IGN_, *_T(6378249.2, 6356515.0, 293.466021294)),

2338 Clarke1880Mod = _lazy('Clarke1880Mod', *_T(6378249.145, 6356514.96639549, 293.466307656)), # aka Clarke1880Arc

2339 CPM1799 = _lazy('CPM1799', *_T(6375738.7, 6356671.92557493, 334.39)), # Comm. des Poids et Mesures

2340 Delambre1810 = _lazy('Delambre1810', *_T(6376428.0, 6355957.92616372, 311.5)), # Belgium

2341 Engelis1985 = _lazy('Engelis1985', *_T(6378136.05, 6356751.32272154, 298.2566)),

2342# Everest1830 = _lazy('Everest1830', *_T(6377276.345, _0_0, 300.801699997)),

2343# Everest1948 = _lazy('Everest1948', *_T(6377304.063, _0_0, 300.801699997)),

2344# Everest1956 = _lazy('Everest1956', *_T(6377301.243, _0_0, 300.801699997)),

2345 Everest1969 = _lazy('Everest1969', *_T(6377295.664, 6356094.667915, 300.801699997)),

2346 Everest1975 = _lazy('Everest1975', *_T(6377299.151, 6356098.14512013, 300.8017255)),

2347 Fisher1968 = _lazy('Fisher1968', *_T(6378150.0, 6356768.33724438, 298.3)),

2348# Fisher1968Mod = _lazy('Fisher1968Mod', *_T(6378155.0, _0_0, 298.3)),

2349 GEM10C = _lazy('GEM10C', *_T(R_MA, 6356752.31424783, 298.2572236)),

2350 GPES = _lazy('GPES', *_T(6378135.0, 6356750.0, _0_0)), # "Gen. Purpose Earth Spheroid"

2351 GRS67 = _lazy('GRS67', *_T(6378160.0, _0_0, 298.247167427)), # Lucerne b=6356774.516

2352# GRS67Truncated = _lazy('GRS67Truncated', *_T(6378160.0, _0_0, 298.25)),

2353 GRS80 = _lazy(_GRS80_, *_T(R_MA, 6356752.314140347, 298.25722210088)), # IUGG, ITRS, ETRS89

2354# Hayford1924 = _lazy('Hayford1924', *_T(6378388.0, 6356911.94612795, None)), # aka Intl1924 f_=297

2355 Helmert1906 = _lazy('Helmert1906', *_T(6378200.0, 6356818.16962789, 298.3)),

2356# Hough1960 = _lazy('Hough1960', *_T(6378270.0, _0_0, 297.0)),

2357 IAU76 = _lazy('IAU76', *_T(6378140.0, _0_0, 298.257)), # Int'l Astronomical Union

2358 IERS1989 = _lazy('IERS1989', *_T(6378136.0, _0_0, 298.257)), # b=6356751.302

2359 IERS1992TOPEX = _lazy('IERS1992TOPEX', *_T(6378136.3, 6356751.61659215, 298.257223563)), # IERS/TOPEX/Poseidon/McCarthy

2360 IERS2003 = _lazy('IERS2003', *_T(6378136.6, 6356751.85797165, 298.25642)),

2361 Intl1924 = _lazy(_Intl1924_, *_T(6378388.0, _0_0, 297.0)), # aka Hayford b=6356911.9462795

2362 Intl1967 = _lazy('Intl1967', *_T(6378157.5, 6356772.2, 298.24961539)), # New Int'l

2363 Krassovski1940 = _lazy(_Krassovski1940_, *_T(6378245.0, 6356863.01877305, 298.3)), # spelling

2364 Krassowsky1940 = _lazy(_Krassowsky1940_, *_T(6378245.0, 6356863.01877305, 298.3)), # spelling

2365# Kaula = _lazy('Kaula', *_T(6378163.0, _0_0, 298.24)), # Kaula 1961

2366# Lerch = _lazy('Lerch', *_T(6378139.0, _0_0, 298.257)), # Lerch 1979

2367 Maupertuis1738 = _lazy('Maupertuis1738', *_T(6397300.0, 6363806.28272251, 191.0)), # France

2368 Mercury1960 = _lazy('Mercury1960', *_T(6378166.0, 6356784.28360711, 298.3)),

2369 Mercury1968Mod = _lazy('Mercury1968Mod', *_T(6378150.0, 6356768.33724438, 298.3)),

2370# MERIT = _lazy('MERIT', *_T(6378137.0, _0_0, 298.257)), # MERIT 1983

2371# NWL10D = _lazy('NWL10D', *_T(6378135.0, _0_0, 298.26)), # Naval Weapons Lab.

2372 NWL1965 = _lazy('NWL1965', *_T(6378145.0, 6356759.76948868, 298.25)), # Naval Weapons Lab.

2373# NWL9D = _lazy('NWL9D', *_T(6378145.0, 6356759.76948868, 298.25)), # NWL1965

2374 OSU86F = _lazy('OSU86F', *_T(6378136.2, 6356751.51693008, 298.2572236)),

2375 OSU91A = _lazy('OSU91A', *_T(6378136.3, 6356751.6165948, 298.2572236)),

2376# Plessis1817 = _lazy('Plessis1817', *_T(6397523.0, 6355863.0, 153.56512242)), # XXX incorrect?

2377 Plessis1817 = _lazy('Plessis1817', *_T(6376523.0, 6355862.93325557, 308.64)), # XXX IGN France 1972

2378# Prolate = _lazy('Prolate', *_T(6356752.3, R_MA, _0_0)),

2379 PZ90 = _lazy('PZ90', *_T(6378136.0, _0_0, 298.257839303)), # GLOSNASS PZ-90 and PZ-90.11

2380# SEAsia = _lazy('SEAsia', *_T(6378155.0, _0_0, 298.3)), # SouthEast Asia

2381 SGS85 = _lazy('SGS85', *_T(6378136.0, 6356751.30156878, 298.257)), # Soviet Geodetic System

2382 SoAmerican1969 = _lazy('SoAmerican1969', *_T(6378160.0, 6356774.71919531, 298.25)), # South American

2383 Sphere = _lazy(_Sphere_, *_T(R_M, R_M, _0_0)), # pseudo

2384 SphereAuthalic = _lazy('SphereAuthalic', *_T(R_FM, R_FM, _0_0)), # pseudo

2385 SpherePopular = _lazy('SpherePopular', *_T(R_MA, R_MA, _0_0)), # EPSG:3857 Spheroid

2386 Struve1860 = _lazy('Struve1860', *_T(6378298.3, 6356657.14266956, 294.73)),

2387# Walbeck = _lazy('Walbeck', *_T(6376896.0, _0_0, 302.78)),

2388# WarOffice = _lazy('WarOffice', *_T(6378300.0, _0_0, 296.0)),

2389 WGS60 = _lazy('WGS60', *_T(6378165.0, 6356783.28695944, 298.3)),

2390 WGS66 = _lazy('WGS66', *_T(6378145.0, 6356759.76948868, 298.25)),

2391 WGS72 = _lazy(_WGS72_, *_T(6378135.0, _0_0, 298.26)), # b=6356750.52

2392 WGS84 = _lazy(_WGS84_, *_T(R_MA, _0_0, _f__WGS84)) # GPS b=6356752.3142451793

2393)

2394

2395_EWGS84 = Ellipsoids.WGS84 # (INTERNAL) shared

2396

2397 if __name__ == '__main__':

2398

2399 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_

2400 from pygeodesy import nameof, printf

2401

2402 for E in (_EWGS84, Ellipsoids.GRS80, # NAD83,

2403 Ellipsoids.Sphere, Ellipsoids.SpherePopular,

2404 Ellipsoid(_EWGS84.b, _EWGS84.a, name='_Prolate')):

2405 e = f2n(E.f) - E.n

2406 printf('# %s: %s', _DOT_('Ellipsoids', E.name), E.toStr(prec=10), nl=1)

2407 printf('# e=%s, f_=%s, f=%s, n=%s (%s)', fstr(E.e, prec=13, fmt=Fmt.e),

2408 fstr(E.f_, prec=13, fmt=Fmt.e),

2409 fstr(E.f, prec=13, fmt=Fmt.e),

2410 fstr(E.n, prec=13, fmt=Fmt.e),

2411 fstr(e, prec=9, fmt=Fmt.e))

2412 printf('# %s %s', Ellipsoid.AlphaKs.name, fstr(E.AlphaKs, prec=20))

2413 printf('# %s %s', Ellipsoid.BetaKs.name, fstr(E.BetaKs, prec=20))

2414 printf('# %s %s', nameof(Ellipsoid.KsOrder), E.KsOrder) # property

2415

2416 # __doc__ of this file, force all into registry

2417 t = [NN] + Ellipsoids.toRepr(all=True, asorted=True).split(_NL_)

2418 printf(_NLATvar_.join(i.strip(_COMMA_) for i in t))

2419

2420# % python3 -m pygeodesy.ellipsoids

2421

2422# 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

2423# e=8.1819190842622e-02, f_=2.98257223563e+02, f=3.3528106647475e-03, n=1.6792203863837e-03 (0.0e+00)

2424# AlphaKs 0.00083773182062446994, 0.00000076085277735725, 0.00000000119764550324, 0.00000000000242917068, 0.00000000000000571182, 0.0000000000000000148, 0.00000000000000000004, 0.0

2425# BetaKs 0.00083773216405794875, 0.0000000590587015222, 0.00000000016734826653, 0.00000000000021647981, 0.00000000000000037879, 0.00000000000000000072, 0.0, 0.0

2426# KsOrder 8

2427

2428# 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

2429# e=8.1819191042833e-02, f_=2.9825722210088e+02, f=3.3528106811837e-03, n=1.6792203946295e-03 (0.0e+00)

2430# AlphaKs 0.00083773182472890429, 0.00000076085278481561, 0.00000000119764552086, 0.00000000000242917073, 0.00000000000000571182, 0.0000000000000000148, 0.00000000000000000004, 0.0

2431# BetaKs 0.0008377321681623882, 0.00000005905870210374, 0.000000000167348269, 0.00000000000021647982, 0.00000000000000037879, 0.00000000000000000072, 0.0, 0.0

2432# KsOrder 8

2433

2434# 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

2435# e=0.0e+00, f_=0.0e+00, f=0.0e+00, n=0.0e+00 (0.0e+00)

2436# AlphaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2437# BetaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2438# KsOrder 8

2439

2440# 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

2441# e=0.0e+00, f_=0.0e+00, f=0.0e+00, n=0.0e+00 (0.0e+00)

2442# AlphaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2443# BetaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2444# KsOrder 8

2445

2446# 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

2447# e=8.2094437949696e-02, f_=-2.97257223563e+02, f=-3.3640898209765e-03, n=-1.6792203863837e-03 (0.0e+00)

2448# AlphaKs -0.00084149152514366627, 0.00000076653480614871, -0.00000000120934503389, 0.0000000000024576225, -0.00000000000000578863, 0.00000000000000001502, -0.00000000000000000004, 0.0

2449# BetaKs -0.00084149187224351817, 0.00000005842735196773, -0.0000000001680487236, 0.00000000000021706261, -0.00000000000000038002, 0.00000000000000000073, -0.0, 0.0

2450# KsOrder 8

2451

2452# **) MIT License

2453#

2454# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.

2455#

2456# Permission is hereby granted, free of charge, to any person obtaining a

2457# copy of this software and associated documentation files (the "Software"),

2458# to deal in the Software without restriction, including without limitation

2459# the rights to use, copy, modify, merge, publish, distribute, sublicense,

2460# and/or sell copies of the Software, and to permit persons to whom the

2461# Software is furnished to do so, subject to the following conditions:

2462#

2463# The above copyright notice and this permission notice shall be included

2464# in all copies or substantial portions of the Software.

2465#

2466# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS

2467# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

2468# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL

2469# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR

2470# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,

2471# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR

2472# OTHER DEALINGS IN THE SOFTWARE.