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, _xinstanceof

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

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

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

70 _0_25, _3_0 # PYCHOK used!

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

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

73 hypot1, hypot2, sqrt3

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

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

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

77 _Intl1924_, _incompatible_, _invalid_, _Krassovski1940_, \

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

79 _prime_vertical_, _radius_, _Sphere_, _SPACE_, _vs_, _WGS72_, \

80 _WGS84_

81 from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS

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

83 _NamedEnumItem, _NamedTuple, _Pass

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 atand, atan2b, atan2d, 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.08.20'

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{float}), 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 and either the

218 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{float}, 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 = atand(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 = atand(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 = atand(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(atan(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 atan2d(s, c) # == atand(s / c) if c else copysign0(_90_0, lat)

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 *= atan(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{float}), 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{float}), 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 _xinstanceof(*self._Geodesics, exact=g)

1025 if g.ellipsoid != self:

1026 raise _ValueError(exact=g, ellipsosid=g.ellipsoid)

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.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._hartzell3d2(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 order != self._KsOrder:

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)).E} 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 @Property_RO

1253 def _L_90(self):

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

1255 '''

1256 return self.L / _90_0

1257

1258 def Llat(self, lat):

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

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

1261

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

1263

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

1265 hemisphere (C{meter}).

1266 '''

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

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

1269

1270 Lmeridian = Llat # meridional distance

1271

1272 @deprecated_Property_RO

1273 def majoradius(self): # PYCHOK no cover

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

1275 return self.a

1276

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

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

1279 along a parallel of (geodetic) latitude.

1280

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

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

1283

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

1285

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

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

1288

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

1290 '''

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

1292

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

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

1295 along a parallel of (geodetic) latitude.

1296

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

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

1299

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

1301

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

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

1304

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

1306 '''

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

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

1309

1310 @deprecated_Property_RO

1311 def minoradius(self): # PYCHOK no cover

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

1313 return self.b

1314

1315 @Property_RO

1316 def n(self):

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

1318 '''

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

1320

1321 flattening = f

1322 flattening1st = f

1323 flattening2nd = f2

1324 flattening3rd = n

1325

1326 polaradius = b # Rpolar

1327

1328# @Property_RO

1329# def Q(self):

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

1331#

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

1333#

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

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

1336# '''

1337# n = self.n

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

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

1340

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

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

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

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

1345

1346 @deprecated_Property_RO

1347 def quarteradius(self): # PYCHOK no cover

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

1349 return self.L

1350

1351 @Property_RO

1352 def R1(self):

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

1354

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

1356 and method C{Rgeometric}.

1357 '''

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

1359 return Radius(R1=r)

1360

1361 Rmean = R1

1362

1363 @Property_RO

1364 def R2(self):

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

1366

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

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

1369 '''

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

1371

1372 Rauthalic = R2

1373

1374# @Property_RO

1375# def R2(self):

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

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

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

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

1380

1381 @Property_RO

1382 def R2x(self):

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

1384

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

1386 '''

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

1388

1389 Rauthalicx = R2x

1390

1391 @Property_RO

1392 def R3(self):

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

1394

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

1396 '''

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

1398 return Radius(R3=r)

1399

1400 Rvolumetric = R3

1401

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

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

1404 along a parallel of (geodetic) latitude.

1405

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

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

1408

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

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

1411

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

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

1414

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

1416 '''

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

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

1419

1420 @Property_RO

1421 def Rbiaxial(self):

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

1423

1424 @see: C{Rtriaxial}

1425 '''

1426 a, b = self.a, self.b

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

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

1429 return Radius(Rbiaxial=q)

1430

1431 Requatorial = a # for consistent naming

1432

1433 def Rgeocentric(self, lat):

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

1435

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

1437

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

1439

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

1441

1442 @see: U{Geocentric Radius

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

1444 '''

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

1446 if a and self.f:

1447 if fabs(a) < PI_2:

1448 s2, c2 = _s2_c2(a)

1449 b2_a2_s2 = self.b2_a2 * s2

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

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

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

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

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

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

1456 else:

1457 r = self.b

1458 return Radius(Rgeocentric=r)

1459

1460 @Property_RO

1461 def Rgeometric(self):

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

1463

1464 @see: C{R1}.

1465 '''

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

1467 return Radius(Rgeometric=g)

1468

1469 def rhumb_(self, exact=True):

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

1471

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

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

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

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

1476

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

1478

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

1480

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

1482 '''

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

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

1485 else:

1486 r = exact

1487 _xinstanceof(*self._Rhumbs, exact=r)

1488 if r.ellipsoid != self:

1489 raise _ValueError(exact=r, ellipsosid=r.ellipsoid)

1490 return r

1491

1492 @Property_RO

1493 def rhumbaux(self):

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

1495 '''

1496 # if not self.isEllipsoidal:

1497 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1499

1500 @property_RO

1501 def _Rhumbs(self):

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

1503 '''

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

1505 _MODS.rhumbsolve.RhumbSolve,

1506 _MODS.rhumbx.Rhumb)

1507 return t

1508

1509 @property

1510 def rhumbsolve(self):

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

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

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

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

1515 '''

1516 # if not self.isEllipsoidal:

1517 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1519

1520 @rhumbsolve.setter # PYCHOK setter!

1521 def rhumbsolve(self, path):

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

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

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

1525 '''

1526 self._rhumbsolve = path

1527

1528 @Property_RO

1529 def rhumbx(self):

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

1531 '''

1532 # if not self.isEllipsoidal:

1533 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1535

1536 def Rlat(self, lat):

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

1538

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

1540

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

1542

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

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

1545

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

1547

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

1549

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

1551

1552 @see: Method C{circle4}.

1553 '''

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

1555 r = self.a

1556 if self.f and lat: # .isEllipsoidal

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

1558 r = Radius(Rlat=r)

1559 return r

1560

1561 Rpolar = b # for consistent naming

1562

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

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

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

1566

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

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

1569 to use an alternate formula.

1570

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

1572

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

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

1575

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

1577 '''

1578 if not self.f: # .isSpherical

1579 n = self.a

1580 elif ca is None:

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

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

1583 elif ca: # derived from EcefYou.forward

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

1585 n = self.a / h

1586 elif sa:

1587 n = self.a2_b / fabs(sa)

1588 else:

1589 n = self.a

1590 return n

1591

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

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

1594 radii of curvature of (geodetic) latitude.

1595

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

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

1598

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

1600 the radii of curvature.

1601

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

1603

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

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

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

1607 Earth_radius#Radii_of_curvature>}.

1608 '''

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

1610

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

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

1613 curvature of (geodetic) latitude.

1614

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

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

1617

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

1619 radii of curvature.

1620

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

1622

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

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

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

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

1627 '''

1628 a = fabs(Phi(phi))

1629 if self.f:

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

1631 if r > EPS02:

1632 n = self.a / sqrt(r)

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

1634 else:

1635 m = n = _0_0 # PYCHOK attr

1636 else:

1637 m = n = self.a

1638 if scaled and a:

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

1640 return Curvature2Tuple(Radius(rocMeridional=m),

1641 Radius(rocPrimeVertical=n))

1642

1643 def rocBearing(self, lat, bearing):

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

1645 latitude and compass direction.

1646

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

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

1649

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

1651

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

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

1654

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

1656

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

1658 '''

1659 if self.f:

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

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

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

1663 c2 *= n / m

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

1665 s2 *= m / n

1666 n = m

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

1668 else:

1669 b = self.b # == self.a

1670 return Radius(rocBearing=b)

1671

1672 @Property_RO

1673 def rocEquatorial2(self):

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

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

1676

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

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

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

1680 '''

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

1682 Radius(rocPrimeVertical0=self.a))

1683

1684 def rocGauss(self, lat):

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

1686

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

1688

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

1690

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

1692

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

1694 Earth_radius#Radii_of_curvature>}

1695 '''

1696 # using ...

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

1698 # return sqrt(m * n)

1699 # ... requires 1 or 2 sqrt

1700 g = self.b

1701 if self.f:

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

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

1704 return Radius(rocGauss=g)

1705

1706 def rocMean(self, lat):

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

1708

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

1710

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

1712

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

1714

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

1716 Earth_radius#Radii_of_curvature>}

1717 '''

1718 if self.f:

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

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

1721 else:

1722 m = self.a

1723 return Radius(rocMean=m)

1724

1725 def rocMeridional(self, lat):

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

1727

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

1729

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

1731

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

1733

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

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

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

1737 '''

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

1739 self.rocEquatorial2.meridional

1740

1741 rocPolar = a2_b # synonymous

1742

1743 def rocPrimeVertical(self, lat):

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

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

1746

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

1748

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

1750

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

1752

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

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

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

1756 Earth_radius#Radii_of_curvature>}.

1757 '''

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

1759 self.rocEquatorial2.prime_vertical

1760

1761 rocTransverse = rocPrimeVertical # synonymous

1762

1763 @deprecated_Property_RO

1764 def Rquadratic(self): # PYCHOK no cover

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

1766 return self.Rbiaxial

1767

1768 @deprecated_Property_RO

1769 def Rr(self): # PYCHOK no cover

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

1771 return self.Rrectifying

1772

1773 @Property_RO

1774 def Rrectifying(self):

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

1776

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

1778 '''

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

1780 return Radius(Rrectifying=r)

1781

1782 @deprecated_Property_RO

1783 def Rs(self): # PYCHOK no cover

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

1785 return self.Rgeometric

1786

1787 @Property_RO

1788 def Rtriaxial(self):

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

1790

1791 @see: C{Rbiaxial}

1792 '''

1793 a, b = self.a, self.b

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

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

1796 return Radius(Rtriaxial=q)

1797

1798 def toEllipsoid2(self, name=NN):

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

1800

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

1802

1803 @see: Property C{a_f}.

1804 '''

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

1806

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

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

1809

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

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

1812 this ellipsoid's name.

1813

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

1815 '''

1816 E = Ellipsoid

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

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

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

1820

1821 def toTriaxial(self, name=NN):

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

1823

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

1825 pointing east and C{Z} pointing north.

1826

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

1828 '''

1829 T = self._triaxial

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

1831

1832 @Property_RO

1833 def _triaxial(self):

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

1835 '''

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

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

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

1839

1840 @Property_RO

1841 def volume(self):

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

1843

1844 @see: C{R3}.

1845 '''

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

1847

1848

1849 class Ellipsoid2(Ellipsoid):

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

1851 '''

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

1853 '''New L{Ellipsoid2}.

1854

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

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

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

1858

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

1860

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

1862

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

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

1865 '''

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

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

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

1869 else:

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

1871

1872

1873 def _spherical_a_b(a, b):

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

1875 '''

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

1877

1878

1879 def _spherical_f(f):

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

1881 '''

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

1883

1884

1885 def _spherical_f_(f_):

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

1887 '''

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

1889

1890

1891 def a_b2e(a, b):

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

1893

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

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

1896

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

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

1899

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

1901 '''

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

1903

1904

1905 def a_b2e2(a, b):

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

1907

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

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

1910

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

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

1913

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

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

1916 '''

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

1918

1919

1920 def a_b2e22(a, b):

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

1922

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

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

1925

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

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

1928

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

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

1931 '''

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

1933

1934

1935 def a_b2e32(a, b):

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

1937

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

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

1940

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

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

1943

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

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

1946 '''

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

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

1949

1950

1951 def a_b2f(a, b):

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

1953

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

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

1956

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

1958

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

1960 for I{near-spherical} ellipsoids.

1961 '''

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

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

1964

1965

1966 def a_b2f_(a, b):

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

1968

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

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

1971

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

1973

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

1975 for I{near-spherical} ellipsoids.

1976 '''

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

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

1979

1980

1981 def a_b2f2(a, b):

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

1983

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

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

1986

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

1988

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

1990 for I{near-spherical} ellipsoids.

1991 '''

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

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

1994

1995

1996 def a_b2n(a, b):

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

1998

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

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

2001

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

2003

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

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

2006 '''

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

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

2009

2010

2011 def a_f2b(a, f):

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

2013

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

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

2016

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

2018 '''

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

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

2021

2022

2023 def a_f_2b(a, f_):

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

2025

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

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

2028

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

2030 '''

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

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

2033

2034

2035 def b_f2a(b, f):

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

2037

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

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

2040

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

2042 '''

2043 t = _1_0 - f

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

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

2046

2047

2048 def b_f_2a(b, f_):

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

2050

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

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

2053

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

2055 '''

2056 t = f_ - _1_0

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

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

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

2060

2061

2062 def e2f(e):

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

2064

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

2066

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

2068

2069 @see: Function L{e22f}.

2070 '''

2071 return e22f(e**2)

2072

2073

2074 def e22f(e2):

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

2076

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

2078

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

2080 '''

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

2082

2083

2084 def f2e2(f):

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

2086

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

2088

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

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

2091

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

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

2094

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

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

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

2098 '''

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

2100

2101

2102 def f2e22(f):

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

2104

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

2106

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

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

2109

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

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

2112

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

2114 html/classGeographicLib_1_1Ellipsoid.html>}.

2115 '''

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

2117 t = (_1_0 - f)**2

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

2119

2120

2121 def f2e32(f):

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

2123

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

2125

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

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

2128

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

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

2131

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

2133 html/classGeographicLib_1_1Ellipsoid.html>}.

2134 '''

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

2136 e2 = f2e2(f)

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

2138

2139

2140 def f_2f(f_):

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

2142

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

2144

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

2146

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

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

2149 '''

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

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

2152

2153

2154 def f2f_(f):

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

2156

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

2158

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

2160

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

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

2163 '''

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

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

2166

2167

2168 def f2f2(f):

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

2170

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

2172

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

2174

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

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

2177

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

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

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

2181 '''

2182 t = _1_0 - f

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

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

2185

2186

2187 def f2n(f):

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

2189

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

2191

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

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

2194

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

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

2197

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

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

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

2201 '''

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

2203

2204

2205 def n2e2(n):

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

2207

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

2209

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

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

2212

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

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

2215

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

2217 '''

2218 t = (n + _1_0)**2

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

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

2221

2222

2223 def n2f(n):

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

2225

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

2227

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

2229

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

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

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

2233 '''

2234 t = n + _1_0

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

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

2237

2238

2239 def n2f_(n):

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

2241

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

2243

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

2245

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

2247 '''

2248 return f2f_(n2f(n))

2249

2250

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

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

2253

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

2255 '''

2256 a, b = E.a, E.b

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

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

2259

2260 a2 = a - b * E.b_a

2261 b2 = b - a * E.a_b

2262 tx = ty = _SQRT2_2

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

2264 ex = a2 * tx**3

2265 ey = b2 * ty**3

2266

2267 qx = px - ex

2268 qy = py - ey

2269 q = hypot(qx, qy)

2270 if q < EPS0: # PYCHOK no cover

2271 break

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

2273

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

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

2276 t = hypot(ty, tx)

2277 if t < EPS0: # PYCHOK no cover

2278 break

2279 tx = tx / t # /= chokes PyChecker

2280 ty = ty / t

2281 if fabs(sx - tx) < EPS and \

2282 fabs(sy - ty) < EPS:

2283 break

2284

2285 tx *= a / px

2286 ty *= b / py

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

2288

2289

2290 class Ellipsoids(_NamedEnum):

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

2292 to accommodate the L{_LazyNamedEnumItem} properties.

2293 '''

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

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

2296 '''

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

2298

2299 Ellipsoids = Ellipsoids(Ellipsoid) # PYCHOK singleton

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2378)

2379

2380_EWGS84 = Ellipsoids.WGS84 # (INTERNAL) shared

2381

2382 if __name__ == '__main__':

2383

2384 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_

2385 from pygeodesy import nameof, printf

2386

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

2388 Ellipsoids.Sphere, Ellipsoids.SpherePopular,

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

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

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

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

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

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

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

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

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

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

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

2400

2401 # __doc__ of this file, force all into registry

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

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

2404

2405# % python3 -m pygeodesy.ellipsoids

2406

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

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

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

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

2411# KsOrder 8

2412

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

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

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

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

2417# KsOrder 8

2418

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

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

2421# AlphaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2422# BetaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2423# KsOrder 8

2424

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

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

2427# AlphaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2428# BetaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2429# KsOrder 8

2430

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

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

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

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

2435# KsOrder 8

2436

2437# **) MIT License

2438#

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

2440#

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

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

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

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

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

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

2447#

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

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

2450#

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

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

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

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

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

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

2457# OTHER DEALINGS IN THE SOFTWARE.