1

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

3

4 u'''Ellipsoidal and spherical earth models.

5

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

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

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

9

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

11

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

13

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

62'''

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

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

65

66 from pygeodesy.basics import copysign0, 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_, property_RO

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

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

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

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

90

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

92

93__all__ = _ALL_LAZY.ellipsoids

94__version__ = '23.10.27'

95

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

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

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

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

100

101

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

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

104 '''

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

106

107

108 def _s2_c2(phi):

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

110 '''

111 if phi:

112 s2 = sin(phi)**2

113 if s2 > EPS:

114 c2 = _1_0 - s2

115 if c2 > EPS:

116 if c2 < EPS1:

117 return s2, c2

118 else:

119 return _1_0, _0_0 # phi == PI_2

120 return _0_0, _1_0 # phi == 0

121

122

123 class a_f2Tuple(_NamedTuple):

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

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

126

127 @see: Class L{Ellipsoid2}.

128 '''

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

130 _Units_ = (_Pass, _Pass)

131

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

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

134

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

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

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

138

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

140

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

142

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

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

145 '''

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

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

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

149 f = _f_0_0

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

151

152 @Property_RO

153 def b(self):

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

155 '''

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

157

158 def ellipsoid(self, name=NN):

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

160

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

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

163 '''

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

165

166 @Property_RO

167 def f_(self):

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

169 '''

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

171

172

173 class Circle4Tuple(_NamedTuple):

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

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

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

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

178

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

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

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

182 equals C{lat}.

183 '''

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

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

186

187

188 class Curvature2Tuple(_NamedTuple):

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

190 C{meter}, conventionally.

191 '''

192 _Names_ = (_meridional_, _prime_vertical_)

193 _Units_ = ( Meter, Meter)

194

195 @property_RO

196 def transverse(self):

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

198 '''

199 return self.prime_vertical

200

201

202 class Ellipsoid(_NamedEnumItem):

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

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

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

206 '''

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

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

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

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

211

212 _geodsolve = NN # means, use PYGEODESY_GEODSOLVE

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

214 _rhumbsolve = NN # means, use PYGEODESY_RHUMBSOLVE

215

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

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

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

219

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

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

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

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

224 spherical).

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

226

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

228

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

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

231

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

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

234 '''

235 ff_ = f, f_ # assertion below

236 try:

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

238 if not _isfinite(a):

239 raise ValueError(_SPACE_(_a_, _not_finite_))

240

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

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

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

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

245 elif f is not None:

246 f = Float(f=f)

247 b = a_f2b(a, f)

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

249 elif f_:

250 f_ = Float(f_=f_)

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

252 f = f_2f(f_)

253 else: # only a, spherical

254 f_ = f = 0

255 b = a # superfluous

256

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

258 raise ValueError(_SPACE_(_f_, _invalid_))

259 if not _isfinite(b):

260 raise ValueError(_SPACE_(_b_, _not_finite_))

261

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

263 b = a

264 f = _f_0_0

265 f_ = _f__0_0

266

267 except (TypeError, ValueError) as x:

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

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

270 raise _ValueError(t, cause=x)

271

272 self._a = a

273 self._b = b

274 self._f = f

275 self._f_ = f_

276

277 self._register(Ellipsoids, name)

278

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

280 d = dict(eps=_TOL)

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

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

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

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

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

286

287 def __eq__(self, other):

288 '''Compare this and an other ellipsoid.

289

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

291

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

293 '''

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

295 self.a == other.a and

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

297

298 def __hash__(self):

299 return self._hash # memoized

300

301 @Property_RO

302 def a(self):

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

304 '''

305 return self._a

306

307 equatoradius = a # = Requatorial

308

309 @Property_RO

310 def a2(self):

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

312 '''

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

314

315 @Property_RO

316 def a2_(self):

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

318 '''

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

320

321 @Property_RO

322 def a_b(self):

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

324 '''

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

326

327 @Property_RO

328 def a2_b(self):

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

330

331 @see: U{Radii of Curvature

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

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

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

335

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

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

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

339 '''

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

341

342 @Property_RO

343 def a2_b2(self):

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

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

346 '''

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

348

349 @Property_RO

350 def a_f(self):

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

352 '''

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

354

355 @Property_RO

356 def A(self):

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

358

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

360 '''

361 A, n = self.a, self.n

362 if n:

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

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

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

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

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

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

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

370 return A

371

372 @Property_RO

373 def _albersCyl(self):

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

375 '''

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

377

378 @Property_RO

379 def AlphaKs(self):

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

381 '''

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

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

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

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

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

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

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

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

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

391 _T(1424729850961/743921418240)) # PYCHOK unaligned

392

393 @Property_RO

394 def area(self):

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

396

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

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

399 '''

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

401

402 @Property_RO

403 def areax(self):

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

405 accurate for very I{oblate} ellipsoids.

406

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

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

409 '''

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

411

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

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

414 '''

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

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

417 t = (v, _vs_, val)

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

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

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

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

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

423 eps=eps, f0=f0, **name_value))

424

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

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

427

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

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

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

431

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

433

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

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

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

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

438 '''

439 if self.f:

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

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

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

443

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

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

446

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

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

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

450

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

452

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

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

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

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

457 '''

458 if self.f:

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

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

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

462

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

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

465

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

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

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

469

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

471

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

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

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

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

476 '''

477 if self.f:

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

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

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

481

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

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

484

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

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

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

488

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

490

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

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

493 thereof is the original geodetic latitude.

494

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

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

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

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

499 '''

500 if self.f:

501 r = Phi_(lat, clip=0)

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

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

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

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

506

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

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

509

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

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

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

513

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

515

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

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

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

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

520 '''

521 if self.f:

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

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

524

525 auxReduced = auxParametric # synonymous

526

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

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

529

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

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

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

533

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

535

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

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

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

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

540 '''

541 if self.f:

542 lat = Lat(lat)

543 if 0 < fabs(lat) < _90_0:

544 if inverse:

545 e = self._elliptic_e22

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

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

548 else:

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

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

551

552 @Property_RO

553 def b(self):

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

555 '''

556 return self._b

557

558 polaradius = b # = Rpolar

559

560 @Property_RO

561 def b_a(self):

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

563

564 @see: Property L{f1}.

565 '''

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

567

568 @Property_RO

569 def b2(self):

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

571 '''

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

573

574 @Property_RO

575 def b2_a(self):

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

577

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

579 '''

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

581

582 @Property_RO

583 def b2_a2(self):

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

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

586 '''

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

588

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

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

591 '''

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

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

594 return atan1d(s, c)

595

596 @Property_RO

597 def BetaKs(self):

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

599 '''

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

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

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

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

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

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

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

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

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

609 _T(191773887257/3719607091200)) # PYCHOK unaligned

610

611 @deprecated_Property_RO

612 def c(self): # PYCHOK no cover

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

614 return self.R2

615

616 @Property_RO

617 def c2(self):

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

619

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

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

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

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

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

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

626 '''

627 return self._c2f(False)

628

629 @Property_RO

630 def c2x(self):

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

632 ellipsoids.

633

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

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

636 '''

637 return self._c2f(True)

638

639 def _c2f(self, c2x):

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

641 '''

642 f, c2 = self.f, self.b2

643 if f:

644 e = self.e

645 if e > EPS0:

646 if f > 0: # .isOblate

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

648 elif f < 0: # .isProlate

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

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

651 return c2

652

653 def circle4(self, lat):

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

655

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

657

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

659 instance.

660

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

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

663

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

665

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

667

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

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

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

671 and method C{Rlat}.

672 '''

673 lat = Lat(lat)

674 if lat:

675 b = lat

676 if fabs(lat) < _90_0:

677 if self.f:

678 b = self._beta(lat)

679 z, r = sincos2d(b)

680 r *= self.a

681 z *= self.b

682 else: # near-polar

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

684 else: # equator

685 r = self.a

686 z = lat = b = _0_0

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

688

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

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

691 along a parallel of (geodetic) latitude.

692

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

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

695

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

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

698

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

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

701

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

703 '''

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

705

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

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

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

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

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

711

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

713 and only between points not near-polar}.

714

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

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

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

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

719

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

721 in same units as this ellipsoid's axes.

722

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

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

725

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

727 '''

728 phi0 = Phi_(lat0=lat0)

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

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

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

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

733

734 @Property_RO

735 def e(self):

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

737

738 @see: Property L{es}.

739 '''

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

741

742 @deprecated_Property_RO

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

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

745 return self.e21

746

747# @Property_RO

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

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

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

751

752 @Property_RO

753 def e2(self):

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

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

756 '''

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

758

759 @Property_RO

760 def e2abs(self):

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

762 '''

763 return fabs(self.e2)

764

765 @Property_RO

766 def e21(self):

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

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

769 '''

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

771

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

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

774

775 @Property_RO

776 def e22(self):

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

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

779 '''

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

781

782 @Property_RO

783 def e22abs(self):

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

785 '''

786 return fabs(self.e22)

787

788 @Property_RO

789 def e32(self):

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

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

792 '''

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

794

795 @Property_RO

796 def e32abs(self):

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

798 '''

799 return fabs(self.e32)

800

801 @Property_RO

802 def e4(self):

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

804 '''

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

806

807 eccentricity = e # eccentricity

808# eccentricity2 = e2 # eccentricity squared

809 eccentricity1st2 = e2 # first eccentricity squared

810 eccentricity2nd2 = e22 # second eccentricity squared

811 eccentricity3rd2 = e32 # third eccentricity squared

812

813 def ecef(self, Ecef=None):

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

815

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

817

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

819

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

821

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

823 '''

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

825

826 @Property_RO

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

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

829 '''

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

831 return _MODS.elliptic.Elliptic(-e12)

832

833 @Property_RO

834 def _elliptic_e22(self): # aka ._elliptic_ep2

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

836 '''

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

838

839 equatoradius = a # Requatorial

840

841 def e2s(self, s):

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

843

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

845

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

847

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

849 '''

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

851

852 def e2s2(self, s):

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

854

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

856

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

858

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

860 '''

861 r = _1_0

862 if self.e2:

863 try:

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

865 if r < 0:

866 raise ValueError(_negative_)

867 except (TypeError, ValueError) as x:

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

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

870 return r

871

872 @Property_RO

873 def es(self):

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

875

876 @see: Property L{e}.

877 '''

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

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

880

881 def es_atanh(self, x):

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

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

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

885

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

887

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

889 html/classGeographicLib_1_1Math.html>}.

890 '''

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

892

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

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

895 '''

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

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

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

899 _0_0)) # .isSpherical

900

901 @Property_RO

902 def es_c(self):

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

904 '''

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

906

907 def es_tauf(self, taup):

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

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

910

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

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

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

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

915 '''

916 t = Scalar(taup=taup)

917 if self.f: # .isEllipsoidal

918 a = fabs(t)

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

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

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

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

923 if fabs(d) < s:

924 break

925 t = Scalar(tauf=T)

926 return t

927

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

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

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

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

932 '''

933 e = self._1_e21

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

935 _tf2 = self._es_taupf2

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

937 a, h = _tf2(T)

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

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

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

941 yield T, i, d

942

943 def es_taupf(self, tau):

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

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

946

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

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

949 '''

950 t = Scalar(tau=tau)

951 if self.f: # .isEllipsoidal

952 t, _ = self._es_taupf2(t)

953 t = Scalar(taupf=t)

954 return t

955

956 def _es_taupf2(self, tau):

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

958 '''

959 if _isfinite(tau):

960 h = hypot1(tau)

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

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

963 else:

964 a, h = tau, INF

965 return a, h

966

967 @Property_RO

968 def _exp_es_atanh_1(self):

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

970 '''

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

972

973 @Property_RO

974 def f(self):

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

976 '''

977 return self._f

978

979 @Property_RO

980 def f_(self):

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

982 '''

983 return self._f_

984

985 @Property_RO

986 def f1(self):

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

988

989 @see: Property L{b_a}.

990 '''

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

992

993 @Property_RO

994 def f2(self):

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

996 '''

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

998

999 @deprecated_Property_RO

1000 def geodesic(self):

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

1002 return self.geodesicw

1003

1004 def geodesic_(self, exact=True):

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

1006

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

1008 ...)}, otherwise a C{geodesic} for this ellipsoid

1009 (C{.geodesicw}, C{.geodesicx}, C{.geodsolve}) or an

1010 other L{Geodesic}, L{GeodesicExact} or L{GeodesicSolve}

1011 instance for I{this} ellipsoid.

1012

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

1014

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

1016

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

1018 '''

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

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

1021 name=self.name)

1022 else:

1023 g = exact

1024 E = _xattr(g, ellipsoid=None)

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

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

1027 return g

1028

1029 @property_RO

1030 def _Geodesics(self):

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

1032 '''

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

1034 _MODS.geodesicx.GeodesicExact,

1035 _MODS.geodsolve.GeodesicSolve)

1036 return t

1037

1038 @property_RO

1039 def geodesicw(self):

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

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

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

1043 package is installed.

1044 '''

1045 # if not self.isEllipsoidal:

1046 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

1047 return _MODS.geodesicw.Geodesic(self)

1048

1049 @property_RO

1050 def geodesicx(self):

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

1052 '''

1053 # if not self.isEllipsoidal:

1054 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1056

1057 @property

1058 def geodsolve(self):

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

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

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

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

1063 '''

1064 # if not self.isEllipsoidal:

1065 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1067

1068 @geodsolve.setter # PYCHOK setter!

1069 def geodsolve(self, path):

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

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

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

1073 '''

1074 self._geodsolve = path

1075

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

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

1078 from a Point-Of-View in space.

1079

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

1081 or L{Vector3d}).

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

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

1084

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

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

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

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

1089

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

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

1092 outside this ellipsoid or points in an opposite

1093 direction.

1094

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

1096

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

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

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

1100 '''

1101 try:

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

1103 except Exception as x:

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

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

1106

1107 @Property_RO

1108 def _hash(self):

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

1110

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

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

1113 this ellipsoid's surface.

1114

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

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

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

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

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

1120

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

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

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

1124

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

1126

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

1128

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

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

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

1132 '''

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

1134 r = v.length

1135

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

1137 if r < EPS0: # EPS

1138 v = v.times(_0_0)

1139 h = -a

1140

1141 elif self.isSpherical:

1142 v = v.times(a / r)

1143 h = r - a

1144

1145 elif normal: # perpendicular to ellipsoid

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

1147 if x < EPS0: # PYCHOK no cover

1148 z = copysign0(b, v.z)

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

1150 h = y - b # polar

1151 elif y < EPS0: # PYCHOK no cover

1152 t = a / r

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

1154 h = x - a # equatorial

1155 else: # normal in 1st quadrant

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

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

1158 h = t.minus(v).length

1159

1160 else: # radial to ellipsoid's center

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

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

1163 v = v.times(t)

1164 h = r * (_1_0 - t)

1165

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

1167 name=self.height4.__name__)

1168

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

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

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

1172 '''

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

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

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

1176

1177 @Property_RO

1178 def isEllipsoidal(self):

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

1180 '''

1181 return self.f != 0

1182

1183 @Property_RO

1184 def isOblate(self):

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

1186 spherical otherwise.

1187 '''

1188 return self.f > 0

1189

1190 @Property_RO

1191 def isProlate(self):

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

1193 spherical otherwise.

1194 '''

1195 return self.f < 0

1196

1197 @Property_RO

1198 def isSpherical(self):

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

1200 '''

1201 return self.f == 0

1202

1203 def _Kseries(self, *AB8Ks):

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

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

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

1207

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

1209 coefficient tuples.

1210

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

1212

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

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

1215 '''

1216 k = self.KsOrder

1217 if self.n:

1218 ns = fpowers(self.n, k)

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

1220 else:

1221 ks = _0_0s(k)

1222 return ks

1223

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

1225 def KsOrder(self):

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

1227 '''

1228 return self._KsOrder

1229

1230 @KsOrder.setter # PYCHOK setter!

1231 def KsOrder(self, order):

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

1233

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

1235 '''

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

1237 raise _ValueError(order=order)

1238 if self._KsOrder != order:

1239 Ellipsoid.AlphaKs._update(self)

1240 Ellipsoid.BetaKs._update(self)

1241 self._KsOrder = order

1242

1243 @Property_RO

1244 def L(self):

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

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

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

1248 '''

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

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

1251

1252 def Llat(self, lat):

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

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

1255

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

1257

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

1259 hemisphere (C{meter}).

1260 '''

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

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

1263

1264 Lmeridian = Llat # meridional distance

1265

1266 @property_RO

1267 def _Lpd(self):

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

1269 '''

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

1271

1272 @property_RO

1273 def _Lpr(self):

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

1275 '''

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

1277

1278 @deprecated_Property_RO

1279 def majoradius(self): # PYCHOK no cover

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

1281 return self.a

1282

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

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

1285 along a parallel of (geodetic) latitude.

1286

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

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

1289

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

1291

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

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

1294

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

1296 '''

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

1298

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

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

1301 along a parallel of (geodetic) latitude.

1302

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

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

1305

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

1307

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

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

1310

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

1312 '''

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

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

1315

1316 @deprecated_Property_RO

1317 def minoradius(self): # PYCHOK no cover

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

1319 return self.b

1320

1321 @Property_RO

1322 def n(self):

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

1324 '''

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

1326

1327 flattening = f

1328 flattening1st = f

1329 flattening2nd = f2

1330 flattening3rd = n

1331

1332 polaradius = b # Rpolar

1333

1334# @Property_RO

1335# def Q(self):

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

1337#

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

1339#

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

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

1342# '''

1343# n = self.n

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

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

1346

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

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

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

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

1351

1352 @deprecated_Property_RO

1353 def quarteradius(self): # PYCHOK no cover

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

1355 return self.L

1356

1357 @Property_RO

1358 def R1(self):

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

1360

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

1362 and method C{Rgeometric}.

1363 '''

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

1365 return Radius(R1=r)

1366

1367 Rmean = R1

1368

1369 @Property_RO

1370 def R2(self):

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

1372

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

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

1375 '''

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

1377

1378 Rauthalic = R2

1379

1380# @Property_RO

1381# def R2(self):

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

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

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

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

1386

1387 @Property_RO

1388 def R2x(self):

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

1390

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

1392 '''

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

1394

1395 Rauthalicx = R2x

1396

1397 @Property_RO

1398 def R3(self):

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

1400

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

1402 '''

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

1404 return Radius(R3=r)

1405

1406 Rvolumetric = R3

1407

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

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

1410 along a parallel of (geodetic) latitude.

1411

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

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

1414

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

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

1417

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

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

1420

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

1422 '''

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

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

1425

1426 @Property_RO

1427 def Rbiaxial(self):

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

1429

1430 @see: C{Rtriaxial}

1431 '''

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

1433 if b < a:

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

1435 elif b > a:

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

1437 return Radius(Rbiaxial=b)

1438

1439 Requatorial = a # for consistent naming

1440

1441 def Rgeocentric(self, lat):

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

1443

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

1445

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

1447

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

1449

1450 @see: U{Geocentric Radius

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

1452 '''

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

1454 if a and self.f:

1455 if fabs(a) < PI_2:

1456 s2, c2 = _s2_c2(a)

1457 b2_a2_s2 = self.b2_a2 * s2

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

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

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

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

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

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

1464 else:

1465 r = self.b

1466 return Radius(Rgeocentric=r)

1467

1468 @Property_RO

1469 def Rgeometric(self):

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

1471

1472 @see: C{R1}.

1473 '''

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

1475 return Radius(Rgeometric=g)

1476

1477 def rhumb_(self, exact=True):

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

1479

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

1481 ...)}, otherwise a C{rhumb} for this ellipsoid (C{.rhumbaux},

1482 C{.rhumbsolve}, C{.rhumbx}) or an other L{Rhumb}, L{RhumbAux}

1483 or L{RhumbSolve} instance for I{this} ellipsoid.

1484

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

1486

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

1488

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

1490 '''

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

1492 r = _MODS.rhumbx.Rhumb(self, exact=exact, name=self.name)

1493 else:

1494 r = exact

1495 E = _xattr(r, ellipsoid=None)

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

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

1498 return r

1499

1500 @property_RO

1501 def rhumbaux(self):

1502 '''Get this ellipsoid's I{Auxiliary} C{rhumbaux.RhumbAux}.

1503 '''

1504 # if not self.isEllipsoidal:

1505 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

1506 return _MODS.rhumbaux.RhumbAux(self, name=self.name)

1507

1508 @property_RO

1509 def _Rhumbs(self):

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

1511 '''

1512 Ellipsoid._Rhumbs = t = (_MODS.rhumbaux.RhumbAux, # overwrite property_RO

1513 _MODS.rhumbsolve.RhumbSolve,

1514 _MODS.rhumbx.Rhumb)

1515 return t

1516

1517 @property

1518 def rhumbsolve(self):

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

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

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

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

1523 '''

1524 # if not self.isEllipsoidal:

1525 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

1526 return _MODS.rhumbsolve.RhumbSolve(self, path=self._rhumbsolve, name=self.name)

1527

1528 @rhumbsolve.setter # PYCHOK setter!

1529 def rhumbsolve(self, path):

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

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

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

1533 '''

1534 self._rhumbsolve = path

1535

1536 @property_RO

1537 def rhumbx(self):

1538 '''Get this ellipsoid's I{Elliptic} C{rhumbx.Rhumb}.

1539 '''

1540 # if not self.isEllipsoidal:

1541 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1543

1544 def Rlat(self, lat):

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

1546

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

1548

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

1550

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

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

1553

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

1555

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

1557

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

1559

1560 @see: Method C{circle4}.

1561 '''

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

1563 r = self.a

1564 if self.f and lat: # .isEllipsoidal

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

1566 r = Radius(Rlat=r)

1567 return r

1568

1569 Rpolar = b # for consistent naming

1570

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

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

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

1574

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

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

1577 to use an alternate formula.

1578

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

1580

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

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

1583

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

1585 '''

1586 if not self.f: # .isSpherical

1587 n = self.a

1588 elif ca is None:

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

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

1591 elif ca: # derived from EcefYou.forward

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

1593 n = self.a / h

1594 elif sa:

1595 n = self.a2_b / fabs(sa)

1596 else:

1597 n = self.a

1598 return n

1599

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

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

1602 radii of curvature of (geodetic) latitude.

1603

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

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

1606

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

1608 the radii of curvature.

1609

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

1611

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

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

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

1615 Earth_radius#Radii_of_curvature>}.

1616 '''

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

1618

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

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

1621 curvature of (geodetic) latitude.

1622

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

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

1625

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

1627 radii of curvature.

1628

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

1630

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

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

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

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

1635 '''

1636 a = fabs(Phi(phi))

1637 if self.f:

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

1639 if r > EPS02:

1640 n = self.a / sqrt(r)

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

1642 else:

1643 m = n = _0_0 # PYCHOK attr

1644 else:

1645 m = n = self.a

1646 if scaled and a:

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

1648 return Curvature2Tuple(Radius(rocMeridional=m),

1649 Radius(rocPrimeVertical=n))

1650

1651 def rocBearing(self, lat, bearing):

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

1653 latitude and compass direction.

1654

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

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

1657

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

1659

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

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

1662

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

1664

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

1666 '''

1667 if self.f:

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

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

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

1671 c2 *= n / m

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

1673 s2 *= m / n

1674 n = m

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

1676 else:

1677 b = self.b # == self.a

1678 return Radius(rocBearing=b)

1679

1680 @Property_RO

1681 def rocEquatorial2(self):

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

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

1684

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

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

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

1688 '''

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

1690 Radius(rocPrimeVertical0=self.a))

1691

1692 def rocGauss(self, lat):

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

1694

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

1696

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

1698

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

1700

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

1702 Earth_radius#Radii_of_curvature>}

1703 '''

1704 # using ...

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

1706 # return sqrt(m * n)

1707 # ... requires 1 or 2 sqrt

1708 g = self.b

1709 if self.f:

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

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

1712 return Radius(rocGauss=g)

1713

1714 def rocMean(self, lat):

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

1716

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

1718

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

1720

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

1722

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

1724 Earth_radius#Radii_of_curvature>}

1725 '''

1726 if self.f:

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

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

1729 else:

1730 m = self.a

1731 return Radius(rocMean=m)

1732

1733 def rocMeridional(self, lat):

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

1735

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

1737

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

1739

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

1741

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

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

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

1745 '''

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

1747 self.rocEquatorial2.meridional

1748

1749 rocPolar = a2_b # synonymous

1750

1751 def rocPrimeVertical(self, lat):

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

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

1754

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

1756

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

1758

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

1760

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

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

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

1764 Earth_radius#Radii_of_curvature>}.

1765 '''

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

1767 self.rocEquatorial2.prime_vertical

1768

1769 rocTransverse = rocPrimeVertical # synonymous

1770

1771 @deprecated_Property_RO

1772 def Rquadratic(self): # PYCHOK no cover

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

1774 return self.Rbiaxial

1775

1776 @deprecated_Property_RO

1777 def Rr(self): # PYCHOK no cover

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

1779 return self.Rrectifying

1780

1781 @Property_RO

1782 def Rrectifying(self):

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

1784

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

1786 '''

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

1788 return Radius(Rrectifying=r)

1789

1790 @deprecated_Property_RO

1791 def Rs(self): # PYCHOK no cover

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

1793 return self.Rgeometric

1794

1795 @Property_RO

1796 def Rtriaxial(self):

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

1798

1799 @see: C{Rbiaxial}

1800 '''

1801 a, b = self.a, self.b

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

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

1804 return Radius(Rtriaxial=q)

1805

1806 def toEllipsoid2(self, name=NN):

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

1808

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

1810

1811 @see: Property C{a_f}.

1812 '''

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

1814

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

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

1817

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

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

1820 this ellipsoid's name.

1821

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

1823 '''

1824 E = Ellipsoid

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

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

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

1828 t = t[:4]

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

1830

1831 def toTriaxial(self, name=NN):

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

1833

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

1835 pointing east and C{Z} pointing north.

1836

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

1838 '''

1839 T = self._triaxial

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

1841

1842 @property_RO

1843 def _triaxial(self):

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

1845 '''

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

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

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

1849

1850 @Property_RO

1851 def volume(self):

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

1853

1854 @see: C{R3}.

1855 '''

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

1857

1858

1859 class Ellipsoid2(Ellipsoid):

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

1861 '''

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

1863 '''New L{Ellipsoid2}.

1864

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

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

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

1868

1869 @raise NameError: Ellipsoid with that B{C{name}} already exists.

1870

1871 @raise ValueError: Invalid B{C{a}} or B{C{f}}.

1872

1873 @note: C{abs(B{f}) < EPS} is forced to C{B{f}=0}, I{spherical}.

1874 Negative C{B{f}} produces a I{prolate} ellipsoid.

1875 '''

1876 if f is None and isinstance(a, Ellipsoid):

1877 Ellipsoid.__init__(self, a.a, f =a.f,

1878 b=a.b, f_=a.f_, name=name)

1879 else:

1880 Ellipsoid.__init__(self, a, f=f, name=name)

1881

1882

1883 def _spherical_a_b(a, b):

1884 '''(INTERNAL) C{True} for spherical or invalid C{a} or C{b}.

1885 '''

1886 return a < EPS0 or b < EPS0 or fabs(a - b) < EPS0

1887

1888

1889 def _spherical_f(f):

1890 '''(INTERNAL) C{True} for spherical or invalid C{f}.

1891 '''

1892 return fabs(f) < EPS or f > EPS1

1893

1894

1895 def _spherical_f_(f_):

1896 '''(INTERNAL) C{True} for spherical or invalid C{f_}.

1897 '''

1898 return fabs(f_) < EPS or fabs(f_) > _1_EPS

1899

1900

1901 def a_b2e(a, b):

1902 '''Return C{e}, the I{1st eccentricity} for a given I{equatorial} and I{polar} radius.

1903

1904 @arg a: Equatorial radius (C{scalar} > 0).

1905 @arg b: Polar radius (C{scalar} > 0).

1906

1907 @return: The I{unsigned}, (1st) eccentricity (C{float} or C{0}),

1908 M{sqrt(1 - (b / a)**2)}.

1909

1910 @note: The result is always I{non-negative} and C{0} for I{near-spherical} ellipsoids.

1911 '''

1912 return Float(e=sqrt(fabs(a_b2e2(a, b)))) # == sqrt(fabs(a - b) * (a + b)) / a)

1913

1914

1915 def a_b2e2(a, b):

1916 '''Return C{e2}, the I{1st eccentricity squared} for a given I{equatorial} and I{polar} radius.

1917

1918 @arg a: Equatorial radius (C{scalar} > 0).

1919 @arg b: Polar radius (C{scalar} > 0).

1920

1921 @return: The I{signed}, (1st) eccentricity I{squared} (C{float} or C{0}),

1922 M{1 - (b / a)**2}.

1923

1924 @note: The result is positive for I{oblate}, negative for I{prolate}

1925 or C{0} for I{near-spherical} ellipsoids.

1926 '''

1927 return Float(e2=_0_0 if _spherical_a_b(a, b) else ((a - b) * (a + b) / a**2))

1928

1929

1930 def a_b2e22(a, b):

1931 '''Return C{e22}, the I{2nd eccentricity squared} for a given I{equatorial} and I{polar} radius.

1932

1933 @arg a: Equatorial radius (C{scalar} > 0).

1934 @arg b: Polar radius (C{scalar} > 0).

1935

1936 @return: The I{signed}, 2nd eccentricity I{squared} (C{float} or C{0}),

1937 M{(a / b)**2 - 1}.

1938

1939 @note: The result is positive for I{oblate}, negative for I{prolate}

1940 or C{0} for I{near-spherical} ellipsoids.

1941 '''

1942 return Float(e22=_0_0 if _spherical_a_b(a, b) else ((a - b) * (a + b) / b**2))

1943

1944

1945 def a_b2e32(a, b):

1946 '''Return C{e32}, the I{3rd eccentricity squared} for a given I{equatorial} and I{polar} radius.

1947

1948 @arg a: Equatorial radius (C{scalar} > 0).

1949 @arg b: Polar radius (C{scalar} > 0).

1950

1951 @return: The I{signed}, 3rd eccentricity I{squared} (C{float} or C{0}),

1952 M{(a**2 - b**2) / (a**2 + b**2)}.

1953

1954 @note: The result is positive for I{oblate}, negative for I{prolate}

1955 or C{0} for I{near-spherical} ellipsoids.

1956 '''

1957 a2, b2 = a**2, b**2

1958 return Float(e32=_0_0 if _spherical_a_b(a2, b2) else ((a2 - b2) / (a2 + b2)))

1959

1960

1961 def a_b2f(a, b):

1962 '''Return C{f}, the I{flattening} for a given I{equatorial} and I{polar} radius.

1963

1964 @arg a: Equatorial radius (C{scalar} > 0).

1965 @arg b: Polar radius (C{scalar} > 0).

1966

1967 @return: The flattening (C{scalar} or C{0}), M{(a - b) / a}.

1968

1969 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0}

1970 for I{near-spherical} ellipsoids.

1971 '''

1972 f = 0 if _spherical_a_b(a, b) else ((a - b) / a)

1973 return _f_0_0 if _spherical_f(f) else Float(f=f)

1974

1975

1976 def a_b2f_(a, b):

1977 '''Return C{f_}, the I{inverse flattening} for a given I{equatorial} and I{polar} radius.

1978

1979 @arg a: Equatorial radius (C{scalar} > 0).

1980 @arg b: Polar radius (C{scalar} > 0).

1981

1982 @return: The inverse flattening (C{scalar} or C{0}), M{a / (a - b)}.

1983

1984 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0}

1985 for I{near-spherical} ellipsoids.

1986 '''

1987 f_ = 0 if _spherical_a_b(a, b) else (a / float(a - b))

1988 return _f__0_0 if _spherical_f_(f_) else Float(f_=f_)

1989

1990

1991 def a_b2f2(a, b):

1992 '''Return C{f2}, the I{2nd flattening} for a given I{equatorial} and I{polar} radius.

1993

1994 @arg a: Equatorial radius (C{scalar} > 0).

1995 @arg b: Polar radius (C{scalar} > 0).

1996

1997 @return: The I{signed}, 2nd flattening (C{scalar} or C{0}), M{(a - b) / b}.

1998

1999 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0}

2000 for I{near-spherical} ellipsoids.

2001 '''

2002 t = 0 if _spherical_a_b(a, b) else float(a - b)

2003 return Float(f2=_0_0 if fabs(t) < EPS0 else (t / b))

2004

2005

2006 def a_b2n(a, b):

2007 '''Return C{n}, the I{3rd flattening} for a given I{equatorial} and I{polar} radius.

2008

2009 @arg a: Equatorial radius (C{scalar} > 0).

2010 @arg b: Polar radius (C{scalar} > 0).

2011

2012 @return: The I{signed}, 3rd flattening (C{scalar} or C{0}), M{(a - b) / (a + b)}.

2013

2014 @note: The result is positive for I{oblate}, negative for I{prolate}

2015 or C{0} for I{near-spherical} ellipsoids.

2016 '''

2017 t = 0 if _spherical_a_b(a, b) else float(a - b)

2018 return Float(n=_0_0 if fabs(t) < EPS0 else (t / (a + b)))

2019

2020

2021 def a_f2b(a, f):

2022 '''Return C{b}, the I{polar} radius for a given I{equatorial} radius and I{flattening}.

2023

2024 @arg a: Equatorial radius (C{scalar} > 0).

2025 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2026

2027 @return: The polar radius (C{float}), M{a * (1 - f)}.

2028 '''

2029 b = a if _spherical_f(f) else (a * (_1_0 - f))

2030 return Radius_(b=a if _spherical_a_b(a, b) else b)

2031

2032

2033 def a_f_2b(a, f_):

2034 '''Return C{b}, the I{polar} radius for a given I{equatorial} radius and I{inverse flattening}.

2035

2036 @arg a: Equatorial radius (C{scalar} > 0).

2037 @arg f_: Inverse flattening (C{scalar} >>> 1).

2038

2039 @return: The polar radius (C{float}), M{a * (f_ - 1) / f_}.

2040 '''

2041 b = a if _spherical_f_(f_) else (a * (f_ - _1_0) / f_)

2042 return Radius_(b=a if _spherical_a_b(a, b) else b)

2043

2044

2045 def b_f2a(b, f):

2046 '''Return C{a}, the I{equatorial} radius for a given I{polar} radius and I{flattening}.

2047

2048 @arg b: Polar radius (C{scalar} > 0).

2049 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2050

2051 @return: The equatorial radius (C{float}), M{b / (1 - f)}.

2052 '''

2053 t = _1_0 - f

2054 a = b if fabs(t) < EPS0 else (b / t)

2055 return Radius_(a=b if _spherical_a_b(a, b) else a)

2056

2057

2058 def b_f_2a(b, f_):

2059 '''Return C{a}, the I{equatorial} radius for a given I{polar} radius and I{inverse flattening}.

2060

2061 @arg b: Polar radius (C{scalar} > 0).

2062 @arg f_: Inverse flattening (C{scalar} >>> 1).

2063

2064 @return: The equatorial radius (C{float}), M{b * f_ / (f_ - 1)}.

2065 '''

2066 t = f_ - _1_0

2067 a = b if _spherical_f_(f_) or fabs(t - f_) < EPS0 \

2068 or fabs(t) < EPS0 else (b * f_ / t)

2069 return Radius_(a=b if _spherical_a_b(a, b) else a)

2070

2071

2072 def e2f(e):

2073 '''Return C{f}, the I{flattening} for a given I{1st eccentricity}.

2074

2075 @arg e: The (1st) eccentricity (0 <= C{float} < 1)

2076

2077 @return: The flattening (C{scalar} or C{0}).

2078

2079 @see: Function L{e22f}.

2080 '''

2081 return e22f(e**2)

2082

2083

2084 def e22f(e2):

2085 '''Return C{f}, the I{flattening} for a given I{1st eccentricity squared}.

2086

2087 @arg e2: The (1st) eccentricity I{squared}, I{signed} (L{NINF} < C{float} < 1)

2088

2089 @return: The flattening (C{float} or C{0}), M{e2 / (sqrt(e2 - 1) + 1)}.

2090 '''

2091 return Float(f=e2 / (sqrt(_1_0 - e2) + _1_0)) if e2 else _f_0_0

2092

2093

2094 def f2e2(f):

2095 '''Return C{e2}, the I{1st eccentricity squared} for a given I{flattening}.

2096

2097 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2098

2099 @return: The I{signed}, (1st) eccentricity I{squared} (C{float} < 1),

2100 M{f * (2 - f)}.

2101

2102 @note: The result is positive for I{oblate}, negative for I{prolate}

2103 or C{0} for I{near-spherical} ellipsoids.

2104

2105 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2106 html/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2107 <https://WikiPedia.org/wiki/Flattening>}.

2108 '''

2109 return Float(e2=_0_0 if _spherical_f(f) else (f * (_2_0 - f)))

2110

2111

2112 def f2e22(f):

2113 '''Return C{e22}, the I{2nd eccentricity squared} for a given I{flattening}.

2114

2115 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2116

2117 @return: The I{signed}, 2nd eccentricity I{squared} (C{float} > -1 or

2118 C{INF}), M{f * (2 - f) / (1 - f)**2}.

2119

2120 @note: The result is positive for I{oblate}, negative for I{prolate}

2121 or C{0} for near-spherical ellipsoids.

2122

2123 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2124 html/classGeographicLib_1_1Ellipsoid.html>}.

2125 '''

2126 # e2 / (1 - e2) == f * (2 - f) / (1 - f)**2

2127 t = (_1_0 - f)**2

2128 return Float(e22=INF if t < EPS0 else (f2e2(f) / t)) # PYCHOK type

2129

2130

2131 def f2e32(f):

2132 '''Return C{e32}, the I{3rd eccentricity squared} for a given I{flattening}.

2133

2134 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2135

2136 @return: The I{signed}, 3rd eccentricity I{squared} (C{float}),

2137 M{f * (2 - f) / (1 + (1 - f)**2)}.

2138

2139 @note: The result is positive for I{oblate}, negative for I{prolate}

2140 or C{0} for I{near-spherical} ellipsoids.

2141

2142 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2143 html/classGeographicLib_1_1Ellipsoid.html>}.

2144 '''

2145 # e2 / (2 - e2) == f * (2 - f) / (1 + (1 - f)**2)

2146 e2 = f2e2(f)

2147 return Float(e32=e2 / (_2_0 - e2))

2148

2149

2150 def f_2f(f_):

2151 '''Return C{f}, the I{flattening} for a given I{inverse flattening}.

2152

2153 @arg f_: Inverse flattening (C{scalar} >>> 1).

2154

2155 @return: The flattening (C{scalar} or C{0}), M{1 / f_}.

2156

2157 @note: The result is positive for I{oblate}, negative for I{prolate}

2158 or C{0} for I{near-spherical} ellipsoids.

2159 '''

2160 f = 0 if _spherical_f_(f_) else _1_0 / f_

2161 return _f_0_0 if _spherical_f(f) else Float(f=f) # PYCHOK type

2162

2163

2164 def f2f_(f):

2165 '''Return C{f_}, the I{inverse flattening} for a given I{flattening}.

2166

2167 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2168

2169 @return: The inverse flattening (C{scalar} or C{0}), M{1 / f}.

2170

2171 @note: The result is positive for I{oblate}, negative for I{prolate}

2172 or C{0} for I{near-spherical} ellipsoids.

2173 '''

2174 f_ = 0 if _spherical_f(f) else _1_0 / f

2175 return _f__0_0 if _spherical_f_(f_) else Float(f_=f_) # PYCHOK type

2176

2177

2178 def f2f2(f):

2179 '''Return C{f2}, the I{2nd flattening} for a given I{flattening}.

2180

2181 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2182

2183 @return: The I{signed}, 2nd flattening (C{scalar} or C{INF}), M{f / (1 - f)}.

2184

2185 @note: The result is positive for I{oblate}, negative for I{prolate}

2186 or C{0} for I{near-spherical} ellipsoids.

2187

2188 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2189 html/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2190 <https://WikiPedia.org/wiki/Flattening>}.

2191 '''

2192 t = _1_0 - f

2193 return Float(f2=_0_0 if _spherical_f(f) else (INF if fabs(t) < EPS

2194 else (f / t))) # PYCHOK type

2195

2196

2197 def f2n(f):

2198 '''Return C{n}, the I{3rd flattening} for a given I{flattening}.

2199

2200 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2201

2202 @return: The I{signed}, 3rd flattening (-1 <= C{float} < 1),

2203 M{f / (2 - f)}.

2204

2205 @note: The result is positive for I{oblate}, negative for I{prolate}

2206 or C{0} for I{near-spherical} ellipsoids.

2207

2208 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2209 html/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2210 <https://WikiPedia.org/wiki/Flattening>}.

2211 '''

2212 return Float(n=_0_0 if _spherical_f(f) else (f / float(_2_0 - f)))

2213

2214

2215 def n2e2(n):

2216 '''Return C{e2}, the I{1st eccentricity squared} for a given I{3rd flattening}.

2217

2218 @arg n: The 3rd flattening (-1 <= C{scalar} < 1).

2219

2220 @return: The I{signed}, (1st) eccentricity I{squared} (C{float} or NINF),

2221 M{4 * n / (1 + n)**2}.

2222

2223 @note: The result is positive for I{oblate}, negative for I{prolate}

2224 or C{0} for I{near-spherical} ellipsoids.

2225

2226 @see: U{Flattening<https://WikiPedia.org/wiki/Flattening>}.

2227 '''

2228 t = (n + _1_0)**2

2229 return Float(e2=_0_0 if fabs(n) < EPS0 else

2230 (NINF if t < EPS0 else (_4_0 * n / t)))

2231

2232

2233 def n2f(n):

2234 '''Return C{f}, the I{flattening} for a given I{3rd flattening}.

2235

2236 @arg n: The 3rd flattening (-1 <= C{scalar} < 1).

2237

2238 @return: The flattening (C{scalar} or NINF), M{2 * n / (1 + n)}.

2239

2240 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2241 html/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2242 <https://WikiPedia.org/wiki/Flattening>}.

2243 '''

2244 t = n + _1_0

2245 f = 0 if fabs(n) < EPS0 else (NINF if t < EPS0 else (_2_0 * n / t))

2246 return _f_0_0 if _spherical_f(f) else Float(f=f)

2247

2248

2249 def n2f_(n):

2250 '''Return C{f_}, the I{inverse flattening} for a given I{3rd flattening}.

2251

2252 @arg n: The 3rd flattening (-1 <= C{scalar} < 1).

2253

2254 @return: The inverse flattening (C{scalar} or C{0}), M{1 / f}.

2255

2256 @see: L{n2f} and L{f2f_}.

2257 '''

2258 return f2f_(n2f(n))

2259

2260

2261 def _normalTo3(px, py, E): # in .height4 above

2262 '''(INTERNAL) Nearest point on a 2-D ellipse in 1st quadrant.

2263

2264 @see: Functions C{.triaxial._normalTo4} and C{-To5}.

2265 '''

2266 a, b = E.a, E.b

2267 if min(px, py, a, b) < EPS0:

2268 raise _AssertionError(px=px, py=py, a=a, b=b, E=E)

2269

2270 a2 = a - b * E.b_a

2271 b2 = b - a * E.a_b

2272 tx = ty = _SQRT2_2

2273 for i in range(16): # max 5

2274 ex = a2 * tx**3

2275 ey = b2 * ty**3

2276

2277 qx = px - ex

2278 qy = py - ey

2279 q = hypot(qx, qy)

2280 if q < EPS0: # PYCHOK no cover

2281 break

2282 r = hypot(ex - tx * a, ey - ty * b) / q

2283

2284 sx, tx = tx, min(_1_0, max(0, (ex + qx * r) / a))

2285 sy, ty = ty, min(_1_0, max(0, (ey + qy * r) / b))

2286 t = hypot(ty, tx)

2287 if t < EPS0: # PYCHOK no cover

2288 break

2289 tx = tx / t # /= chokes PyChecker

2290 ty = ty / t

2291 if fabs(sx - tx) < EPS and \

2292 fabs(sy - ty) < EPS:

2293 break

2294

2295 tx *= a / px

2296 ty *= b / py

2297 return tx, ty, i # x and y as fractions

2298

2299

2300 class Ellipsoids(_NamedEnum):

2301 '''(INTERNAL) L{Ellipsoid} registry, I{must} be a sub-class

2302 to accommodate the L{_LazyNamedEnumItem} properties.

2303 '''

2304 def _Lazy(self, a, b, f_, **kwds):

2305 '''(INTERNAL) Instantiate the L{Ellipsoid}.

2306 '''

2307 return Ellipsoid(a, b=b, f_=f_, **kwds)

2308

2309 Ellipsoids = Ellipsoids(Ellipsoid) # PYCHOK singleton

2310'''Some pre-defined L{Ellipsoid}s, all I{lazily} instantiated.'''

2311# <https://www.GNU.org/software/gama/manual/html_node/Supported-ellipsoids.html>

2312# <https://GSSC.ESA.int/navipedia/index.php/Reference_Frames_in_GNSS>

2313# <https://kb.OSU.edu/dspace/handle/1811/77986>

2314# <https://www.IBM.com/docs/en/db2/11.5?topic=systems-supported-spheroids>

2315# <https://w3.Energistics.org/archive/Epicentre/Epicentre_v3.0/DataModel/LogicalDictionary/StandardValues/ellipsoid.html>

2316# <https://GitHub.com/locationtech/proj4j/blob/master/src/main/java/org/locationtech/proj4j/datum/Ellipsoid.java>

2317 Ellipsoids._assert( # <https://WikiPedia.org/wiki/Earth_ellipsoid>

2318 Airy1830 = _lazy(_Airy1830_, *_T(6377563.396, _0_0, 299.3249646)), # b=6356256.909

2319 AiryModified = _lazy(_AiryModified_, *_T(6377340.189, _0_0, 299.3249646)), # b=6356034.448

2320# APL4_9 = _lazy('APL4_9', *_T(6378137.0, _0_0, 298.24985392)), # Appl. Phys. Lab. 1965

2321# ANS = _lazy('ANS', *_T(6378160.0, _0_0, 298.25)), # Australian Nat. Spheroid

2322# AN_SA96 = _lazy('AN_SA96', *_T(6378160.0, _0_0, 298.24985392)), # Australian Nat. South America

2323 Australia1966 = _lazy('Australia1966', *_T(6378160.0, _0_0, 298.25)), # b=6356774.7192

2324 ATS1977 = _lazy('ATS1977', *_T(6378135.0, _0_0, 298.257)), # "Average Terrestrial System"

2325 Bessel1841 = _lazy(_Bessel1841_, *_T(6377397.155, 6356078.962818, 299.152812797)),

2326 BesselModified = _lazy('BesselModified', *_T(6377492.018, _0_0, 299.1528128)),

2327# BesselNamibia = _lazy('BesselNamibia', *_T(6377483.865, _0_0, 299.1528128)),

2328 CGCS2000 = _lazy('CGCS2000', *_T(R_MA, _0_0, 298.257222101)), # BeiDou Coord System (BDC)

2329# Clarke1858 = _lazy('Clarke1858', *_T(6378293.639, _0_0, 294.260676369)),

2330 Clarke1866 = _lazy(_Clarke1866_, *_T(6378206.4, 6356583.8, 294.978698214)),

2331 Clarke1880 = _lazy('Clarke1880', *_T(6378249.145, 6356514.86954978, 293.465)),

2332 Clarke1880IGN = _lazy(_Clarke1880IGN_, *_T(6378249.2, 6356515.0, 293.466021294)),

2333 Clarke1880Mod = _lazy('Clarke1880Mod', *_T(6378249.145, 6356514.96639549, 293.466307656)), # aka Clarke1880Arc

2334 CPM1799 = _lazy('CPM1799', *_T(6375738.7, 6356671.92557493, 334.39)), # Comm. des Poids et Mesures

2335 Delambre1810 = _lazy('Delambre1810', *_T(6376428.0, 6355957.92616372, 311.5)), # Belgium

2336 Engelis1985 = _lazy('Engelis1985', *_T(6378136.05, 6356751.32272154, 298.2566)),

2337# Everest1830 = _lazy('Everest1830', *_T(6377276.345, _0_0, 300.801699997)),

2338# Everest1948 = _lazy('Everest1948', *_T(6377304.063, _0_0, 300.801699997)),

2339# Everest1956 = _lazy('Everest1956', *_T(6377301.243, _0_0, 300.801699997)),

2340 Everest1969 = _lazy('Everest1969', *_T(6377295.664, 6356094.667915, 300.801699997)),

2341 Everest1975 = _lazy('Everest1975', *_T(6377299.151, 6356098.14512013, 300.8017255)),

2342 Fisher1968 = _lazy('Fisher1968', *_T(6378150.0, 6356768.33724438, 298.3)),

2343# Fisher1968Mod = _lazy('Fisher1968Mod', *_T(6378155.0, _0_0, 298.3)),

2344 GEM10C = _lazy('GEM10C', *_T(R_MA, 6356752.31424783, 298.2572236)),

2345 GPES = _lazy('GPES', *_T(6378135.0, 6356750.0, _0_0)), # "Gen. Purpose Earth Spheroid"

2346 GRS67 = _lazy('GRS67', *_T(6378160.0, _0_0, 298.247167427)), # Lucerne b=6356774.516

2347# GRS67Truncated = _lazy('GRS67Truncated', *_T(6378160.0, _0_0, 298.25)),

2348 GRS80 = _lazy(_GRS80_, *_T(R_MA, 6356752.314140347, 298.25722210088)), # IUGG, ITRS, ETRS89

2349# Hayford1924 = _lazy('Hayford1924', *_T(6378388.0, 6356911.94612795, None)), # aka Intl1924 f_=297

2350 Helmert1906 = _lazy('Helmert1906', *_T(6378200.0, 6356818.16962789, 298.3)),

2351# Hough1960 = _lazy('Hough1960', *_T(6378270.0, _0_0, 297.0)),

2352 IAU76 = _lazy('IAU76', *_T(6378140.0, _0_0, 298.257)), # Int'l Astronomical Union

2353 IERS1989 = _lazy('IERS1989', *_T(6378136.0, _0_0, 298.257)), # b=6356751.302

2354 IERS1992TOPEX = _lazy('IERS1992TOPEX', *_T(6378136.3, 6356751.61659215, 298.257223563)), # IERS/TOPEX/Poseidon/McCarthy

2355 IERS2003 = _lazy('IERS2003', *_T(6378136.6, 6356751.85797165, 298.25642)),

2356 Intl1924 = _lazy(_Intl1924_, *_T(6378388.0, _0_0, 297.0)), # aka Hayford b=6356911.9462795

2357 Intl1967 = _lazy('Intl1967', *_T(6378157.5, 6356772.2, 298.24961539)), # New Int'l

2358 Krassovski1940 = _lazy(_Krassovski1940_, *_T(6378245.0, 6356863.01877305, 298.3)), # spelling

2359 Krassowsky1940 = _lazy(_Krassowsky1940_, *_T(6378245.0, 6356863.01877305, 298.3)), # spelling

2360# Kaula = _lazy('Kaula', *_T(6378163.0, _0_0, 298.24)), # Kaula 1961

2361# Lerch = _lazy('Lerch', *_T(6378139.0, _0_0, 298.257)), # Lerch 1979

2362 Maupertuis1738 = _lazy('Maupertuis1738', *_T(6397300.0, 6363806.28272251, 191.0)), # France

2363 Mercury1960 = _lazy('Mercury1960', *_T(6378166.0, 6356784.28360711, 298.3)),

2364 Mercury1968Mod = _lazy('Mercury1968Mod', *_T(6378150.0, 6356768.33724438, 298.3)),

2365# MERIT = _lazy('MERIT', *_T(6378137.0, _0_0, 298.257)), # MERIT 1983

2366# NWL10D = _lazy('NWL10D', *_T(6378135.0, _0_0, 298.26)), # Naval Weapons Lab.

2367 NWL1965 = _lazy('NWL1965', *_T(6378145.0, 6356759.76948868, 298.25)), # Naval Weapons Lab.

2368# NWL9D = _lazy('NWL9D', *_T(6378145.0, 6356759.76948868, 298.25)), # NWL1965

2369 OSU86F = _lazy('OSU86F', *_T(6378136.2, 6356751.51693008, 298.2572236)),

2370 OSU91A = _lazy('OSU91A', *_T(6378136.3, 6356751.6165948, 298.2572236)),

2371# Plessis1817 = _lazy('Plessis1817', *_T(6397523.0, 6355863.0, 153.56512242)), # XXX incorrect?

2372 Plessis1817 = _lazy('Plessis1817', *_T(6376523.0, 6355862.93325557, 308.64)), # XXX IGN France 1972

2373 PZ90 = _lazy('PZ90', *_T(6378136.0, _0_0, 298.257839303)), # GLOSNASS PZ-90 and PZ-90.11

2374# SEAsia = _lazy('SEAsia', *_T(6378155.0, _0_0, 298.3)), # SouthEast Asia

2375 SGS85 = _lazy('SGS85', *_T(6378136.0, 6356751.30156878, 298.257)), # Soviet Geodetic System

2376 SoAmerican1969 = _lazy('SoAmerican1969', *_T(6378160.0, 6356774.71919531, 298.25)), # South American

2377 Struve1860 = _lazy('Struve1860', *_T(6378298.3, 6356657.14266956, 294.73)),

2378# Walbeck = _lazy('Walbeck', *_T(6376896.0, _0_0, 302.78)),

2379# WarOffice = _lazy('WarOffice', *_T(6378300.0, _0_0, 296.0)),

2380 WGS60 = _lazy('WGS60', *_T(6378165.0, 6356783.28695944, 298.3)),

2381 WGS66 = _lazy('WGS66', *_T(6378145.0, 6356759.76948868, 298.25)),

2382 WGS72 = _lazy(_WGS72_, *_T(6378135.0, _0_0, 298.26)), # b=6356750.52

2383 WGS84 = _lazy(_WGS84_, *_T(R_MA, _0_0, _f__WGS84)), # GPS b=6356752.3142451793

2384# Prolate = _lazy('Prolate', *_T(6356752.3, R_MA, _0_0)),

2385 Sphere = _lazy(_Sphere_, *_T(R_M, R_M, _0_0)), # pseudo

2386 SphereAuthalic = _lazy('SphereAuthalic', *_T(R_FM, R_FM, _0_0)), # pseudo

2387 SpherePopular = _lazy('SpherePopular', *_T(R_MA, R_MA, _0_0)) # EPSG:3857 Spheroid

2388)

2389

2390_EWGS84 = Ellipsoids.WGS84 # (INTERNAL) shared

2391

2392 if __name__ == '__main__':

2393

2394 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_

2395 from pygeodesy import nameof, printf

2396

2397 for E in (_EWGS84, Ellipsoids.GRS80, # NAD83,

2398 Ellipsoids.Sphere, Ellipsoids.SpherePopular,

2399 Ellipsoid(_EWGS84.b, _EWGS84.a, name='_Prolate')):

2400 e = f2n(E.f) - E.n

2401 printf('# %s: %s', _DOT_('Ellipsoids', E.name), E.toStr(prec=10), nl=1)

2402 printf('# e=%s, f_=%s, f=%s, n=%s (%s)', fstr(E.e, prec=13, fmt=Fmt.e),

2403 fstr(E.f_, prec=13, fmt=Fmt.e),

2404 fstr(E.f, prec=13, fmt=Fmt.e),

2405 fstr(E.n, prec=13, fmt=Fmt.e),

2406 fstr(e, prec=9, fmt=Fmt.e))

2407 printf('# %s %s', Ellipsoid.AlphaKs.name, fstr(E.AlphaKs, prec=20))

2408 printf('# %s %s', Ellipsoid.BetaKs.name, fstr(E.BetaKs, prec=20))

2409 printf('# %s %s', nameof(Ellipsoid.KsOrder), E.KsOrder) # property

2410

2411 # __doc__ of this file, force all into registry

2412 t = [NN] + Ellipsoids.toRepr(all=True, asorted=True).split(_NL_)

2413 printf(_NLATvar_.join(i.strip(_COMMA_) for i in t))

2414

2415# % python3 -m pygeodesy.ellipsoids

2416

2417# Ellipsoids.WGS84: name='WGS84', a=6378137, b=6356752.3142451793, f_=298.257223563, f=0.0033528107, f2=0.0033640898, n=0.0016792204, e=0.0818191908, e2=0.00669438, e21=0.99330562, e22=0.0067394967, e32=0.0033584313, A=6367449.1458234144, L=10001965.7293127216, R1=6371008.7714150595, R2=6371007.1809184738, R3=6371000.7900091587, Rbiaxial=6367453.6345163295, Rtriaxial=6372797.5559594007

2418# e=8.1819190842622e-02, f_=2.98257223563e+02, f=3.3528106647475e-03, n=1.6792203863837e-03 (0.0e+00)

2419# AlphaKs 0.00083773182062446994, 0.00000076085277735725, 0.00000000119764550324, 0.00000000000242917068, 0.00000000000000571182, 0.0000000000000000148, 0.00000000000000000004, 0.0

2420# BetaKs 0.00083773216405794875, 0.0000000590587015222, 0.00000000016734826653, 0.00000000000021647981, 0.00000000000000037879, 0.00000000000000000072, 0.0, 0.0

2421# KsOrder 8

2422

2423# 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.7292304579, R1=6371008.7713801153, R2=6371007.1808835147, R3=6371000.7899741363, Rbiaxial=6367453.6344640013, Rtriaxial=6372797.5559332585

2424# e=8.1819191042833e-02, f_=2.9825722210088e+02, f=3.3528106811837e-03, n=1.6792203946295e-03 (0.0e+00)

2425# AlphaKs 0.00083773182472890429, 0.00000076085278481561, 0.00000000119764552086, 0.00000000000242917073, 0.00000000000000571182, 0.0000000000000000148, 0.00000000000000000004, 0.0

2426# BetaKs 0.0008377321681623882, 0.00000005905870210374, 0.000000000167348269, 0.00000000000021647982, 0.00000000000000037879, 0.00000000000000000072, 0.0, 0.0

2427# KsOrder 8

2428

2429# 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

2430# e=0.0e+00, f_=0.0e+00, f=0.0e+00, n=0.0e+00 (0.0e+00)

2431# AlphaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2432# BetaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2433# KsOrder 8

2434

2435# 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

2436# e=0.0e+00, f_=0.0e+00, f=0.0e+00, n=0.0e+00 (0.0e+00)

2437# AlphaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2438# BetaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2439# KsOrder 8

2440

2441# Ellipsoids._Prolate: name='_Prolate', a=6356752.3142451793, b=6378137, f_=-297.257223563, f=-0.0033640898, f2=-0.0033528107, n=-0.0016792204, e=0.0820944379, e2=-0.0067394967, e21=1.0067394967, e22=-0.00669438, e32=-0.0033584313, A=6367449.1458234144, L=10035500.5204500332, R1=6363880.5428301189, R2=6363878.9413582645, R3=6363872.5644020075, Rbiaxial=6367453.6345163295, Rtriaxial=6362105.2243882557

2442# e=8.2094437949696e-02, f_=-2.97257223563e+02, f=-3.3640898209765e-03, n=-1.6792203863837e-03 (0.0e+00)

2443# AlphaKs -0.00084149152514366627, 0.00000076653480614871, -0.00000000120934503389, 0.0000000000024576225, -0.00000000000000578863, 0.00000000000000001502, -0.00000000000000000004, 0.0

2444# BetaKs -0.00084149187224351817, 0.00000005842735196773, -0.0000000001680487236, 0.00000000000021706261, -0.00000000000000038002, 0.00000000000000000073, -0.0, 0.0

2445# KsOrder 8

2446

2447# **) MIT License

2448#

2449# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.

2450#

2451# Permission is hereby granted, free of charge, to any person obtaining a

2452# copy of this software and associated documentation files (the "Software"),

2453# to deal in the Software without restriction, including without limitation

2454# the rights to use, copy, modify, merge, publish, distribute, sublicense,

2455# and/or sell copies of the Software, and to permit persons to whom the

2456# Software is furnished to do so, subject to the following conditions:

2457#

2458# The above copyright notice and this permission notice shall be included

2459# in all copies or substantial portions of the Software.

2460#

2461# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS

2462# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

2463# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL

2464# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR

2465# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,

2466# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR

2467# OTHER DEALINGS IN THE SOFTWARE.