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

95

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

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

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

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

100

101

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

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

104 '''

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

106

107

108 def _s2_c2(phi):

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

110 '''

111 if phi:

112 s2 = sin(phi)**2

113 if s2 > EPS:

114 c2 = _1_0 - s2

115 if c2 > EPS:

116 if c2 < EPS1:

117 return s2, c2

118 else:

119 return _1_0, _0_0 # phi == PI_2

120 return _0_0, _1_0 # phi == 0

121

122

123 class a_f2Tuple(_NamedTuple):

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

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

126

127 @see: Class L{Ellipsoid2}.

128 '''

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

130 _Units_ = (_Pass, _Pass)

131

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

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

134

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

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

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

138

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

140

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

142

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

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

145 '''

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

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

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

149 f = _f_0_0

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

151

152 @Property_RO

153 def b(self):

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

155 '''

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

157

158 def ellipsoid(self, name=NN):

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

160

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

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

163 '''

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

165

166 @Property_RO

167 def f_(self):

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

169 '''

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

171

172

173 class Circle4Tuple(_NamedTuple):

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

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

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

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

178

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

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

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

182 equals C{lat}.

183 '''

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

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

186

187

188 class Curvature2Tuple(_NamedTuple):

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

190 C{meter}, conventionally.

191 '''

192 _Names_ = (_meridional_, _prime_vertical_)

193 _Units_ = ( Meter, Meter)

194

195 @property_RO

196 def transverse(self):

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

198 '''

199 return self.prime_vertical

200

201

202 class Ellipsoid(_NamedEnumItem):

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

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

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

206 '''

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

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

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

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

211

212 _geodsolve = NN # means, use PYGEODESY_GEODSOLVE

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

214 _rhumbsolve = NN # means, use PYGEODESY_RHUMBSOLVE

215

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

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

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

219

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

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

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

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

224 spherical).

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

226

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

228

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

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

231

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

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

234 '''

235 ff_ = f, f_ # assertion below

236 try:

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

238 if not _isfinite(a):

239 raise ValueError(_SPACE_(_a_, _not_finite_))

240

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

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

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

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

245 elif f is not None:

246 f = Float(f=f)

247 b = a_f2b(a, f)

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

249 elif f_:

250 f_ = Float(f_=f_)

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

252 f = f_2f(f_)

253 else: # only a, spherical

254 f_ = f = 0

255 b = a # superfluous

256

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

258 raise ValueError(_SPACE_(_f_, _invalid_))

259 if not _isfinite(b):

260 raise ValueError(_SPACE_(_b_, _not_finite_))

261

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

263 b = a

264 f = _f_0_0

265 f_ = _f__0_0

266

267 except (TypeError, ValueError) as x:

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

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

270 raise _ValueError(t, cause=x)

271

272 self._a = a

273 self._b = b

274 self._f = f

275 self._f_ = f_

276

277 self._register(Ellipsoids, name)

278

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

280 d = dict(eps=_TOL)

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

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

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

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

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

286

287 def __eq__(self, other):

288 '''Compare this and an other ellipsoid.

289

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

291

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

293 '''

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

295 self.a == other.a and

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

297

298 def __hash__(self):

299 return self._hash # memoized

300

301 @Property_RO

302 def a(self):

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

304 '''

305 return self._a

306

307 equatoradius = a # = Requatorial

308

309 @Property_RO

310 def a2(self):

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

312 '''

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

314

315 @Property_RO

316 def a2_(self):

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

318 '''

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

320

321 @Property_RO

322 def a_b(self):

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

324 '''

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

326

327 @Property_RO

328 def a2_b(self):

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

330

331 @see: U{Radii of Curvature

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

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

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

335

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

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

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

339 '''

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

341

342 @Property_RO

343 def a2_b2(self):

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

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

346 '''

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

348

349 @Property_RO

350 def a_f(self):

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

352 '''

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

354

355 @Property_RO

356 def A(self):

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

358

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

360 '''

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

362 if n:

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

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

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

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

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

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

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

370 return A

371

372 @Property_RO

373 def _albersCyl(self):

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

375 '''

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

377

378 @Property_RO

379 def AlphaKs(self):

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

381 '''

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

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

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

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

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

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

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

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

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

391 _T(1424729850961/743921418240)) # PYCHOK unaligned

392

393 @Property_RO

394 def area(self):

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

396

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

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

399 '''

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

401

402 @Property_RO

403 def areax(self):

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

405 accurate for very I{oblate} ellipsoids.

406

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

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

409 '''

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

411

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

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

414 '''

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

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

417 t = (v, _vs_, val)

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

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

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

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

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

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

424

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

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

427

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

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

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

431

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

433

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

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

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

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

438 '''

439 if self.f:

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

441 lat = 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{scalar}), M{(a - b) / a}, C{0} for spherical, negative for prolate.

976 '''

977 return self._f

978

979 @Property_RO

980 def f_(self):

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

982 '''

983 return self._f_

984

985 @Property_RO

986 def f1(self):

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

988

989 @see: Property L{b_a}.

990 '''

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

992

993 @Property_RO

994 def f2(self):

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

996 '''

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

998

999 @deprecated_Property_RO

1000 def geodesic(self):

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

1002 return self.geodesicw

1003

1004 def geodesic_(self, exact=True):

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

1006

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

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

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

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

1011 instance for I{this} ellipsoid.

1012

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

1014

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

1016

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

1018 '''

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

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

1021 name=self.name)

1022 else:

1023 g = exact

1024 _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 self._KsOrder != order:

1239 Ellipsoid.AlphaKs._update(self)

1240 Ellipsoid.BetaKs._update(self)

1241 self._KsOrder = order

1242

1243 @Property_RO

1244 def L(self):

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

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

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

1248 '''

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

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

1251

1252 @Property_RO

1253 def _Lpd(self):

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

1255 '''

1256 return Distance(_Lpd=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 if b < a:

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

1429 elif b > a:

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

1431 return Radius(Rbiaxial=b)

1432

1433 Requatorial = a # for consistent naming

1434

1435 def Rgeocentric(self, lat):

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

1437

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

1439

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

1441

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

1443

1444 @see: U{Geocentric Radius

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

1446 '''

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

1448 if a and self.f:

1449 if fabs(a) < PI_2:

1450 s2, c2 = _s2_c2(a)

1451 b2_a2_s2 = self.b2_a2 * s2

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

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

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

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

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

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

1458 else:

1459 r = self.b

1460 return Radius(Rgeocentric=r)

1461

1462 @Property_RO

1463 def Rgeometric(self):

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

1465

1466 @see: C{R1}.

1467 '''

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

1469 return Radius(Rgeometric=g)

1470

1471 def rhumb_(self, exact=True):

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

1473

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

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

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

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

1478

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

1480

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

1482

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

1484 '''

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

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

1487 else:

1488 r = exact

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

1490 if r.ellipsoid != self:

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

1492 return r

1493

1494 @Property_RO

1495 def rhumbaux(self):

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

1497 '''

1498 # if not self.isEllipsoidal:

1499 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1501

1502 @property_RO

1503 def _Rhumbs(self):

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

1505 '''

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

1507 _MODS.rhumbsolve.RhumbSolve,

1508 _MODS.rhumbx.Rhumb)

1509 return t

1510

1511 @property

1512 def rhumbsolve(self):

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

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

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

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

1517 '''

1518 # if not self.isEllipsoidal:

1519 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1521

1522 @rhumbsolve.setter # PYCHOK setter!

1523 def rhumbsolve(self, path):

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

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

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

1527 '''

1528 self._rhumbsolve = path

1529

1530 @Property_RO

1531 def rhumbx(self):

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

1533 '''

1534 # if not self.isEllipsoidal:

1535 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1537

1538 def Rlat(self, lat):

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

1540

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

1542

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

1544

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

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

1547

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

1549

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

1551

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

1553

1554 @see: Method C{circle4}.

1555 '''

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

1557 r = self.a

1558 if self.f and lat: # .isEllipsoidal

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

1560 r = Radius(Rlat=r)

1561 return r

1562

1563 Rpolar = b # for consistent naming

1564

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

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

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

1568

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

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

1571 to use an alternate formula.

1572

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

1574

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

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

1577

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

1579 '''

1580 if not self.f: # .isSpherical

1581 n = self.a

1582 elif ca is None:

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

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

1585 elif ca: # derived from EcefYou.forward

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

1587 n = self.a / h

1588 elif sa:

1589 n = self.a2_b / fabs(sa)

1590 else:

1591 n = self.a

1592 return n

1593

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

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

1596 radii of curvature of (geodetic) latitude.

1597

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

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

1600

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

1602 the radii of curvature.

1603

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

1605

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

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

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

1609 Earth_radius#Radii_of_curvature>}.

1610 '''

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

1612

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

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

1615 curvature of (geodetic) latitude.

1616

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

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

1619

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

1621 radii of curvature.

1622

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

1624

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

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

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

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

1629 '''

1630 a = fabs(Phi(phi))

1631 if self.f:

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

1633 if r > EPS02:

1634 n = self.a / sqrt(r)

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

1636 else:

1637 m = n = _0_0 # PYCHOK attr

1638 else:

1639 m = n = self.a

1640 if scaled and a:

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

1642 return Curvature2Tuple(Radius(rocMeridional=m),

1643 Radius(rocPrimeVertical=n))

1644

1645 def rocBearing(self, lat, bearing):

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

1647 latitude and compass direction.

1648

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

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

1651

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

1653

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

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

1656

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

1658

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

1660 '''

1661 if self.f:

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

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

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

1665 c2 *= n / m

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

1667 s2 *= m / n

1668 n = m

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

1670 else:

1671 b = self.b # == self.a

1672 return Radius(rocBearing=b)

1673

1674 @Property_RO

1675 def rocEquatorial2(self):

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

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

1678

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

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

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

1682 '''

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

1684 Radius(rocPrimeVertical0=self.a))

1685

1686 def rocGauss(self, lat):

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

1688

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

1690

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

1692

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

1694

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

1696 Earth_radius#Radii_of_curvature>}

1697 '''

1698 # using ...

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

1700 # return sqrt(m * n)

1701 # ... requires 1 or 2 sqrt

1702 g = self.b

1703 if self.f:

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

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

1706 return Radius(rocGauss=g)

1707

1708 def rocMean(self, lat):

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

1710

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

1712

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

1714

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

1716

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

1718 Earth_radius#Radii_of_curvature>}

1719 '''

1720 if self.f:

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

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

1723 else:

1724 m = self.a

1725 return Radius(rocMean=m)

1726

1727 def rocMeridional(self, lat):

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

1729

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

1731

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

1733

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

1735

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

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

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

1739 '''

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

1741 self.rocEquatorial2.meridional

1742

1743 rocPolar = a2_b # synonymous

1744

1745 def rocPrimeVertical(self, lat):

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

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

1748

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

1750

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

1752

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

1754

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

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

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

1758 Earth_radius#Radii_of_curvature>}.

1759 '''

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

1761 self.rocEquatorial2.prime_vertical

1762

1763 rocTransverse = rocPrimeVertical # synonymous

1764

1765 @deprecated_Property_RO

1766 def Rquadratic(self): # PYCHOK no cover

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

1768 return self.Rbiaxial

1769

1770 @deprecated_Property_RO

1771 def Rr(self): # PYCHOK no cover

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

1773 return self.Rrectifying

1774

1775 @Property_RO

1776 def Rrectifying(self):

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

1778

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

1780 '''

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

1782 return Radius(Rrectifying=r)

1783

1784 @deprecated_Property_RO

1785 def Rs(self): # PYCHOK no cover

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

1787 return self.Rgeometric

1788

1789 @Property_RO

1790 def Rtriaxial(self):

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

1792

1793 @see: C{Rbiaxial}

1794 '''

1795 a, b = self.a, self.b

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

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

1798 return Radius(Rtriaxial=q)

1799

1800 def toEllipsoid2(self, name=NN):

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

1802

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

1804

1805 @see: Property C{a_f}.

1806 '''

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

1808

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

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

1811

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

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

1814 this ellipsoid's name.

1815

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

1817 '''

1818 E = Ellipsoid

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

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

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

1822

1823 def toTriaxial(self, name=NN):

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

1825

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

1827 pointing east and C{Z} pointing north.

1828

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

1830 '''

1831 T = self._triaxial

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

1833

1834 @Property_RO

1835 def _triaxial(self):

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

1837 '''

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

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

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

1841

1842 @Property_RO

1843 def volume(self):

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

1845

1846 @see: C{R3}.

1847 '''

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

1849

1850

1851 class Ellipsoid2(Ellipsoid):

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

1853 '''

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

1855 '''New L{Ellipsoid2}.

1856

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

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

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

1860

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

1862

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

1864

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

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

1867 '''

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

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

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

1871 else:

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

1873

1874

1875 def _spherical_a_b(a, b):

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

1877 '''

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

1879

1880

1881 def _spherical_f(f):

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

1883 '''

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

1885

1886

1887 def _spherical_f_(f_):

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

1889 '''

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

1891

1892

1893 def a_b2e(a, b):

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

1895

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

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

1898

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

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

1901

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

1903 '''

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

1905

1906

1907 def a_b2e2(a, b):

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

1909

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

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

1912

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

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

1915

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

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

1918 '''

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

1920

1921

1922 def a_b2e22(a, b):

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

1924

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

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

1927

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

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

1930

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

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

1933 '''

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

1935

1936

1937 def a_b2e32(a, b):

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

1939

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

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

1942

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

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

1945

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

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

1948 '''

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

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

1951

1952

1953 def a_b2f(a, b):

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

1955

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

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

1958

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

1960

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

1962 for I{near-spherical} ellipsoids.

1963 '''

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

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

1966

1967

1968 def a_b2f_(a, b):

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

1970

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

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

1973

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

1975

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

1977 for I{near-spherical} ellipsoids.

1978 '''

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

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

1981

1982

1983 def a_b2f2(a, b):

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

1985

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

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

1988

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

1990

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

1992 for I{near-spherical} ellipsoids.

1993 '''

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

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

1996

1997

1998 def a_b2n(a, b):

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

2000

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

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

2003

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

2005

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

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

2008 '''

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

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

2011

2012

2013 def a_f2b(a, f):

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

2015

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

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

2018

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

2020 '''

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

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

2023

2024

2025 def a_f_2b(a, f_):

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

2027

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

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

2030

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

2032 '''

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

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

2035

2036

2037 def b_f2a(b, f):

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

2039

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

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

2042

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

2044 '''

2045 t = _1_0 - f

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

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

2048

2049

2050 def b_f_2a(b, f_):

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

2052

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

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

2055

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

2057 '''

2058 t = f_ - _1_0

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

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

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

2062

2063

2064 def e2f(e):

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

2066

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

2068

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

2070

2071 @see: Function L{e22f}.

2072 '''

2073 return e22f(e**2)

2074

2075

2076 def e22f(e2):

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

2078

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

2080

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

2082 '''

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

2084

2085

2086 def f2e2(f):

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

2088

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

2090

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

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

2093

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

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

2096

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

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

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

2100 '''

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

2102

2103

2104 def f2e22(f):

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

2106

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

2108

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

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

2111

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

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

2114

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

2116 html/classGeographicLib_1_1Ellipsoid.html>}.

2117 '''

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

2119 t = (_1_0 - f)**2

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

2121

2122

2123 def f2e32(f):

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

2125

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

2127

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

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

2130

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

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

2133

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

2135 html/classGeographicLib_1_1Ellipsoid.html>}.

2136 '''

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

2138 e2 = f2e2(f)

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

2140

2141

2142 def f_2f(f_):

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

2144

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

2146

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

2148

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

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

2151 '''

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

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

2154

2155

2156 def f2f_(f):

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

2158

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

2160

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

2162

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

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

2165 '''

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

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

2168

2169

2170 def f2f2(f):

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

2172

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

2174

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

2176

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

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

2179

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

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

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

2183 '''

2184 t = _1_0 - f

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

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

2187

2188

2189 def f2n(f):

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

2191

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

2193

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

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

2196

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

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

2199

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

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

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

2203 '''

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

2205

2206

2207 def n2e2(n):

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

2209

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

2211

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

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

2214

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

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

2217

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

2219 '''

2220 t = (n + _1_0)**2

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

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

2223

2224

2225 def n2f(n):

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

2227

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

2229

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

2231

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

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

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

2235 '''

2236 t = n + _1_0

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

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

2239

2240

2241 def n2f_(n):

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

2243

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

2245

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

2247

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

2249 '''

2250 return f2f_(n2f(n))

2251

2252

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

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

2255

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

2257 '''

2258 a, b = E.a, E.b

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

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

2261

2262 a2 = a - b * E.b_a

2263 b2 = b - a * E.a_b

2264 tx = ty = _SQRT2_2

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

2266 ex = a2 * tx**3

2267 ey = b2 * ty**3

2268

2269 qx = px - ex

2270 qy = py - ey

2271 q = hypot(qx, qy)

2272 if q < EPS0: # PYCHOK no cover

2273 break

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

2275

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

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

2278 t = hypot(ty, tx)

2279 if t < EPS0: # PYCHOK no cover

2280 break

2281 tx = tx / t # /= chokes PyChecker

2282 ty = ty / t

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

2284 fabs(sy - ty) < EPS:

2285 break

2286

2287 tx *= a / px

2288 ty *= b / py

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

2290

2291

2292 class Ellipsoids(_NamedEnum):

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

2294 to accommodate the L{_LazyNamedEnumItem} properties.

2295 '''

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

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

2298 '''

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

2300

2301 Ellipsoids = Ellipsoids(Ellipsoid) # PYCHOK singleton

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2380)

2381

2382_EWGS84 = Ellipsoids.WGS84 # (INTERNAL) shared

2383

2384 if __name__ == '__main__':

2385

2386 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_

2387 from pygeodesy import nameof, printf

2388

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

2390 Ellipsoids.Sphere, Ellipsoids.SpherePopular,

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

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

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

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

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

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

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

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

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

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

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

2402

2403 # __doc__ of this file, force all into registry

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

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

2406

2407# % python3 -m pygeodesy.ellipsoids

2408

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

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

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

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

2413# KsOrder 8

2414

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

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

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

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

2419# KsOrder 8

2420

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

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

2423# AlphaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2424# BetaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2425# KsOrder 8

2426

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

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

2429# AlphaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2430# BetaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2431# KsOrder 8

2432

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

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

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

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

2437# KsOrder 8

2438

2439# **) MIT License

2440#

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

2442#

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

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

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

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

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

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

2449#

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

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

2452#

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

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

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

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

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

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

2459# OTHER DEALINGS IN THE SOFTWARE.