1

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

3

4 u'''Ellipsoidal and spherical earth models.

5

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

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

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

9

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

11

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

13

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

63'''

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

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

66

67 from pygeodesy.basics import copysign0, isbool, isint

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

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

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

71 _0_25, _3_0 # PYCHOK used!

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

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

74 hypot1, hypot2, sqrt3, Fsum

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

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

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

78 _Intl1924_, _incompatible_, _invalid_, _Krassovski1940_, \

79 _Krassowsky1940_, _lat_, _meridional_, _negative_, _not_finite_, \

80 _prime_vertical_, _radius_, _Sphere_, _SPACE_, _vs_, _WGS72_, _WGS84_

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

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

83 _NamedEnumItem, _NamedTuple, _Pass, _ALL_LAZY, _MODS

84 from pygeodesy.namedTuples import Distance2Tuple, Vector3Tuple, Vector4Tuple

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

86 deprecated_property_RO, property_RO

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

88# from pygeodesy.triaxials import _hartzell3 # _MODS

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

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

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

92

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

94

95__all__ = _ALL_LAZY.ellipsoids

96__version__ = '24.05.28'

97

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

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

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

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

102

103

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

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

106 '''

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

108

109

110 def _s2_c2(phi):

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

112 '''

113 if phi:

114 s2 = sin(phi)**2

115 if s2 > EPS:

116 c2 = _1_0 - s2

117 if c2 > EPS:

118 if c2 < EPS1:

119 return s2, c2

120 else:

121 return _1_0, _0_0 # phi == PI_2

122 return _0_0, _1_0 # phi == 0

123

124

125 class a_f2Tuple(_NamedTuple):

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

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

128

129 @see: Class L{Ellipsoid2}.

130 '''

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

132 _Units_ = (_Pass, _Pass)

133

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

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

136

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

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

139 @kwarg name: Optional C{B{name}=NN} (C{str}).

140

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

142

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

144

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

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

147 '''

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

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

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

151 f = _f_0_0

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

153

154 @Property_RO

155 def b(self):

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

157 '''

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

159

160 def ellipsoid(self, **name):

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

162

163 @kwarg name: Optional C{B{name}=NN} (C{str}).

164

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

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

167 '''

168 return Ellipsoid(self.a, f=self.f, name=self._name__(name)) # PYCHOK .a and .f

169

170 @Property_RO

171 def f_(self):

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

173 '''

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

175

176

177 class Circle4Tuple(_NamedTuple):

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

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

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

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

182

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

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

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

186 equals C{lat}.

187 '''

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

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

190

191

192 class Curvature2Tuple(_NamedTuple):

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

194 C{meter}, conventionally.

195 '''

196 _Names_ = (_meridional_, _prime_vertical_)

197 _Units_ = ( Meter, Meter)

198

199 @property_RO

200 def transverse(self):

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

202 '''

203 return self.prime_vertical

204

205

206 class Ellipsoid(_NamedEnumItem):

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

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

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

210 '''

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

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

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

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

215

216 _geodsolve = NN # means, use PYGEODESY_GEODSOLVE

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

218 _rhumbsolve = NN # means, use PYGEODESY_RHUMBSOLVE

219

220 def __init__(self, a, b=None, f_=None, f=None, **name):

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

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

223

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

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

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

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

228 spherical).

229 @kwarg name: Optional, unique C{B{name}=NN} (C{str}).

230

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

232

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

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

235

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

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

238 '''

239 ff_ = f, f_ # assertion below

240 n = _name__(**name) if name else NN

241 try:

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

243 if not _isfinite(a):

244 raise ValueError(_SPACE_(_a_, _not_finite_))

245

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

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

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

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

250 elif f is not None:

251 f = Float(f=f)

252 b = a_f2b(a, f)

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

254 elif f_:

255 f_ = Float(f_=f_)

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

257 f = f_2f(f_)

258 else: # only a, spherical

259 f_ = f = 0

260 b = a # superfluous

261

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

263 raise ValueError(_SPACE_(_f_, _invalid_))

264 if not _isfinite(b):

265 raise ValueError(_SPACE_(_b_, _not_finite_))

266

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

268 b = a

269 f = _f_0_0

270 f_ = _f__0_0

271

272 except (TypeError, ValueError) as x:

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

274 t = instr(self, a=a, name=n, **d)

275 raise _ValueError(t, cause=x)

276

277 self._a = a

278 self._b = b

279 self._f = f

280 self._f_ = f_

281

282 self._register(Ellipsoids, n)

283

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

285 d = dict(eps=_TOL)

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

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

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

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

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

291

292 def __eq__(self, other):

293 '''Compare this and an other ellipsoid.

294

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

296

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

298 '''

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

300 self.a == other.a and

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

302

303 def __hash__(self):

304 return self._hash # memoized

305

306 @Property_RO

307 def a(self):

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

309 '''

310 return self._a

311

312 equatoradius = a # = Requatorial

313

314 @Property_RO

315 def a2(self):

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

317 '''

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

319

320 @Property_RO

321 def a2_(self):

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

323 '''

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

325

326 @Property_RO

327 def a_b(self):

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

329 '''

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

331

332 @Property_RO

333 def a2_b(self):

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

335

336 @see: U{Radii of Curvature

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

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

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

340

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

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

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

344 '''

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

346

347 @Property_RO

348 def a2_b2(self):

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

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

351 '''

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

353

354 @Property_RO

355 def a_f(self):

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

357 '''

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

359

360 @Property_RO

361 def A(self):

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

363

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

365 '''

366 A, n = self.a, self.n

367 if n:

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

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

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

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

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

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

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

375 return A

376

377 @Property_RO

378 def _albersCyl(self):

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

380 '''

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

382

383 @Property_RO

384 def AlphaKs(self):

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

386 '''

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

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

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

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

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

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

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

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

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

396 _T(1424729850961/743921418240)) # PYCHOK unaligned

397

398 @Property_RO

399 def area(self):

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

401

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

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

404 '''

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

406

407 @Property_RO

408 def areax(self):

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

410 accurate for very I{oblate} ellipsoids.

411

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

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

414 '''

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

416

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

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

419 '''

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

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

422 t = (v, _vs_, val)

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

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

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

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

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

428 eps=eps, f0=f0, **name_value))

429

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

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

432

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

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

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

436

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

438

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

440 C++/doc/classGeographicLib_1_1Ellipsoid.html>}, U{Authalic latitude

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

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

443 '''

444 if self.f:

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

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

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

448

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

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

451

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

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

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

455

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

457

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

459 C++/doc/classGeographicLib_1_1Ellipsoid.html>}, U{Conformal latitude

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

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

462 '''

463 if self.f:

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

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

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

467

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

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

470

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

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

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

474

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

476

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

478 C++/doc/classGeographicLib_1_1Ellipsoid.html>}, U{Geocentric latitude

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

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

481 '''

482 if self.f:

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

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

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

486

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

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

489

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

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

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

493

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

495

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

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

498 thereof is the original geodetic latitude.

499

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

501 C++/doc/classGeographicLib_1_1Ellipsoid.html>}, U{Isometric latitude

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

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

504 '''

505 if self.f:

506 r = Phi_(lat, clip=0)

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

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

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

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

511

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

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

514

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

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

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

518

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

520

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

522 C++/doc/classGeographicLib_1_1Ellipsoid.html>}, U{Parametric latitude

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

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

525 '''

526 if self.f:

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

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

529

530 auxReduced = auxParametric # synonymous

531

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

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

534

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

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

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

538

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

540

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

542 C++/doc/classGeographicLib_1_1Ellipsoid.html>}, U{Rectifying latitude

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

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

545 '''

546 if self.f:

547 lat = Lat(lat)

548 if 0 < fabs(lat) < _90_0:

549 if inverse:

550 e = self._elliptic_e22

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

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

553 else:

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

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

556

557 @Property_RO

558 def b(self):

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

560 '''

561 return self._b

562

563 polaradius = b # = Rpolar

564

565 @Property_RO

566 def b_a(self):

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

568

569 @see: Property L{f1}.

570 '''

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

572

573 @Property_RO

574 def b2(self):

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

576 '''

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

578

579 @Property_RO

580 def b2_a(self):

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

582

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

584 '''

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

586

587 @Property_RO

588 def b2_a2(self):

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

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

591 '''

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

593

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

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

596 '''

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

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

599 return atan1d(s, c)

600

601 @Property_RO

602 def BetaKs(self):

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

604 '''

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

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

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

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

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

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

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

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

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

614 _T(191773887257/3719607091200)) # PYCHOK unaligned

615

616 @deprecated_Property_RO

617 def c(self): # PYCHOK no cover

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

619 return self.R2

620

621 @Property_RO

622 def c2(self):

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

624

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

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

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

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

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

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

631 '''

632 return self._c2f(False)

633

634 @Property_RO

635 def c2x(self):

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

637 ellipsoids.

638

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

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

641 '''

642 return self._c2f(True)

643

644 def _c2f(self, c2x):

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

646 '''

647 f, c2 = self.f, self.b2

648 if f:

649 e = self.e

650 if e > EPS0:

651 if f > 0: # .isOblate

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

653 elif f < 0: # .isProlate

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

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

656 return c2

657

658 def circle4(self, lat):

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

660

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

662

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

664 instance.

665

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

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

668

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

670

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

672

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

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

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

676 and method C{Rlat}.

677 '''

678 lat = Lat(lat)

679 if lat:

680 b = lat

681 if fabs(lat) < _90_0:

682 if self.f:

683 b = self._beta(lat)

684 z, r = sincos2d(b)

685 r *= self.a

686 z *= self.b

687 else: # near-polar

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

689 else: # equator

690 r = self.a

691 z = lat = b = _0_0

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

693

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

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

696 along a parallel of (geodetic) latitude.

697

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

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

700

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

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

703

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

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

706

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

708 '''

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

710

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

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

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

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

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

716

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

718 and only between points not near-polar}.

719

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

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

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

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

724

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

726 in same units as this ellipsoid's axes.

727

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

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

730

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

732 '''

733 phi0 = Phi_(lat0=lat0)

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

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

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

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

738

739 @Property_RO

740 def e(self):

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

742

743 @see: Property L{es}.

744 '''

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

746

747 @deprecated_Property_RO

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

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

750 return self.e21

751

752# @Property_RO

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

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

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

756

757 @Property_RO

758 def e2(self):

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

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

761 '''

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

763

764 @Property_RO

765 def e2abs(self):

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

767 '''

768 return fabs(self.e2)

769

770 @Property_RO

771 def e21(self):

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

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

774 '''

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

776

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

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

779

780 @Property_RO

781 def e22(self):

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

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

784 '''

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

786

787 @Property_RO

788 def e22abs(self):

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

790 '''

791 return fabs(self.e22)

792

793 @Property_RO

794 def e32(self):

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

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

797 '''

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

799

800 @Property_RO

801 def e32abs(self):

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

803 '''

804 return fabs(self.e32)

805

806 @Property_RO

807 def e4(self):

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

809 '''

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

811

812 eccentricity = e # eccentricity

813# eccentricity2 = e2 # eccentricity squared

814 eccentricity1st2 = e2 # first eccentricity squared

815 eccentricity2nd2 = e22 # second eccentricity squared

816 eccentricity3rd2 = e32 # third eccentricity squared

817

818 def ecef(self, Ecef=None):

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

820

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

822

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

824

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

826

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

828 '''

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

830

831 @Property_RO

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

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

834 '''

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

836 return _MODS.elliptic.Elliptic(e12)

837

838 @Property_RO

839 def _elliptic_e22(self): # aka ._elliptic_ep2

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

841 '''

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

843

844 equatoradius = a # Requatorial

845

846 def e2s(self, s):

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

848

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

850

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

852

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

854 '''

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

856

857 def e2s2(self, s):

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

859

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

861

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

863

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

865 '''

866 r = _1_0

867 if self.e2:

868 try:

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

870 if r < 0:

871 raise ValueError(_negative_)

872 except (TypeError, ValueError) as x:

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

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

875 return r

876

877 @Property_RO

878 def es(self):

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

880

881 @see: Property L{e}.

882 '''

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

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

885

886 def es_atanh(self, x):

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

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

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

890

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

892

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

894 C++/doc/classGeographicLib_1_1Math.html>}.

895 '''

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

897

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

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

900 '''

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

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

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

904 _0_0)) # .isSpherical

905

906 @Property_RO

907 def es_c(self):

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

909 '''

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

911

912 def es_tauf(self, taup):

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

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

915

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

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

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

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

920 '''

921 t = Scalar(taup=taup)

922 if self.f: # .isEllipsoidal

923 a = fabs(t)

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

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

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

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

928 if fabs(d) < s:

929 break

930 t = Scalar(tauf=T)

931 return t

932

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

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

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

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

937 '''

938 e = self._1_e21

939 _F2_ = Fsum(T).fsum2f_ # τ0

940 _tf2 = self._es_taupf2

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

942 a, h = _tf2(T)

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

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

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

946 yield T, i, d

947

948 def es_taupf(self, tau):

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

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

951

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

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

954 '''

955 t = Scalar(tau=tau)

956 if self.f: # .isEllipsoidal

957 t, _ = self._es_taupf2(t)

958 t = Scalar(taupf=t)

959 return t

960

961 def _es_taupf2(self, tau):

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

963 '''

964 if _isfinite(tau):

965 h = hypot1(tau)

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

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

968 else:

969 a, h = tau, INF

970 return a, h

971

972 @Property_RO

973 def _exp_es_atanh_1(self):

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

975 '''

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

977

978 @Property_RO

979 def f(self):

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

981 '''

982 return self._f

983

984 @Property_RO

985 def f_(self):

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

987 '''

988 return self._f_

989

990 @Property_RO

991 def f1(self):

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

993

994 @see: Property L{b_a}.

995 '''

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

997

998 @Property_RO

999 def f2(self):

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

1001 '''

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

1003

1004 @deprecated_Property_RO

1005 def geodesic(self):

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

1007 return self.geodesicw

1008

1009 def geodesic_(self, exact=True):

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

1011

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

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

1014 instance for I{this} ellipsoid.

1015

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

1017

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

1019

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

1021 '''

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

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

1024 name=self.name)

1025 else:

1026 g = exact

1027 E = _xattr(g, ellipsoid=None)

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

1029 raise _ValueError(exact=g, ellipsoid=E, txt_not_=self.name)

1030 return g

1031

1032 @property_RO

1033 def _Geodesics(self):

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

1035 '''

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

1037 _MODS.geodesicx.GeodesicExact,

1038 _MODS.geodsolve.GeodesicSolve)

1039 return t

1040

1041 @property_RO

1042 def geodesicw(self):

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

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

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

1046 package is installed.

1047 '''

1048 # if not self.isEllipsoidal:

1049 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

1050 return _MODS.geodesicw.Geodesic(self)

1051

1052 @property_RO

1053 def geodesicx(self):

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

1055 '''

1056 # if not self.isEllipsoidal:

1057 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1059

1060 @property

1061 def geodsolve(self):

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

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

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

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

1066 '''

1067 # if not self.isEllipsoidal:

1068 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1070

1071 @geodsolve.setter # PYCHOK setter!

1072 def geodsolve(self, path):

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

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

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

1076 '''

1077 self._geodsolve = path

1078

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

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

1081 from a Point-Of-View in space.

1082

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

1084 or L{Vector3d}).

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

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

1087

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

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

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

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

1092

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

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

1095 outside this ellipsoid or points in an opposite

1096 direction.

1097

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

1099

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

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

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

1103 '''

1104 try:

1105 v, d, _ = _MODS.triaxials._hartzell3(pov, los, self._triaxial)

1106 except Exception as x:

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

1108 return Vector4Tuple(v.x, v.y, v.z, d, name__=self.hartzell4)

1109

1110 @Property_RO

1111 def _hash(self):

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

1113

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

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

1116 this ellipsoid's surface.

1117

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

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

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

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

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

1123

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

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

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

1127

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

1129

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

1131

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

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

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

1135 '''

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

1137 r = v.length

1138

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

1140 if r < EPS0: # EPS

1141 v = v.times(_0_0)

1142 h = -a

1143

1144 elif self.isSpherical:

1145 v = v.times(a / r)

1146 h = r - a

1147

1148 elif normal: # perpendicular to ellipsoid

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

1150 if x < EPS0: # PYCHOK no cover

1151 z = copysign0(b, v.z)

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

1153 h = y - b # polar

1154 elif y < EPS0: # PYCHOK no cover

1155 t = a / r

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

1157 h = x - a # equatorial

1158 else: # normal in 1st quadrant

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

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

1161 h = t.minus(v).length

1162

1163 else: # radial to ellipsoid's center

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

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

1166 v = v.times(t)

1167 h = r * (_1_0 - t)

1168

1169 return Vector4Tuple(v.x, v.y, v.z, h, iteration=i, name__=self.height4)

1170

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

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

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

1174 '''

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

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

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

1178

1179 @Property_RO

1180 def isEllipsoidal(self):

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

1182 '''

1183 return self.f != 0

1184

1185 @Property_RO

1186 def isOblate(self):

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

1188 spherical otherwise.

1189 '''

1190 return self.f > 0

1191

1192 @Property_RO

1193 def isProlate(self):

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

1195 spherical otherwise.

1196 '''

1197 return self.f < 0

1198

1199 @Property_RO

1200 def isSpherical(self):

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

1202 '''

1203 return self.f == 0

1204

1205 def _Kseries(self, *AB8Ks):

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

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

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

1209

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

1211 coefficient tuples.

1212

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

1214

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

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

1217 '''

1218 k = self.KsOrder

1219 if self.n:

1220 ns = fpowers(self.n, k)

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

1222 else:

1223 ks = _0_0s(k)

1224 return ks

1225

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

1227 def KsOrder(self):

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

1229 '''

1230 return self._KsOrder

1231

1232 @KsOrder.setter # PYCHOK setter!

1233 def KsOrder(self, order):

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

1235

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

1237 '''

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

1239 raise _ValueError(order=order)

1240 if self._KsOrder != order:

1241 Ellipsoid.AlphaKs._update(self)

1242 Ellipsoid.BetaKs._update(self)

1243 self._KsOrder = order

1244

1245 @Property_RO

1246 def L(self):

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

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

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

1250 '''

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

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

1253

1254 def Llat(self, lat):

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

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

1257

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

1259

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

1261 hemisphere (C{meter}).

1262 '''

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

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

1265

1266 Lmeridian = Llat # meridional distance

1267

1268 @property_RO

1269 def _Lpd(self):

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

1271 '''

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

1273

1274 @property_RO

1275 def _Lpr(self):

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

1277 '''

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

1279

1280 @deprecated_Property_RO

1281 def majoradius(self): # PYCHOK no cover

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

1283 return self.a

1284

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

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

1287 along a parallel of (geodetic) latitude.

1288

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

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

1291

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

1293

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

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

1296

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

1298 '''

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

1300

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

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

1303 along a parallel of (geodetic) latitude.

1304

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

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

1307

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

1309

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

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

1312

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

1314 '''

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

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

1317

1318 @deprecated_Property_RO

1319 def minoradius(self): # PYCHOK no cover

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

1321 return self.b

1322

1323 @Property_RO

1324 def n(self):

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

1326 '''

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

1328

1329 flattening = f

1330 flattening1st = f

1331 flattening2nd = f2

1332 flattening3rd = n

1333

1334 polaradius = b # Rpolar

1335

1336# @Property_RO

1337# def Q(self):

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

1339#

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

1341#

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

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

1344# '''

1345# n = self.n

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

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

1348

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

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

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

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

1353

1354 @deprecated_Property_RO

1355 def quarteradius(self): # PYCHOK no cover

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

1357 return self.L

1358

1359 @Property_RO

1360 def R1(self):

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

1362

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

1364 and method C{Rgeometric}.

1365 '''

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

1367 return Radius(R1=r)

1368

1369 Rmean = R1

1370

1371 @Property_RO

1372 def R2(self):

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

1374

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

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

1377 '''

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

1379

1380 Rauthalic = R2

1381

1382# @Property_RO

1383# def R2(self):

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

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

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

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

1388

1389 @Property_RO

1390 def R2x(self):

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

1392

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

1394 '''

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

1396

1397 Rauthalicx = R2x

1398

1399 @Property_RO

1400 def R3(self):

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

1402

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

1404 '''

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

1406 return Radius(R3=r)

1407

1408 Rvolumetric = R3

1409

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

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

1412 along a parallel of (geodetic) latitude.

1413

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

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

1416

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

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

1419

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

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

1422

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

1424 '''

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

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

1427

1428 @Property_RO

1429 def Rbiaxial(self):

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

1431

1432 @see: C{Rtriaxial}

1433 '''

1434 a, b = self.a, self.b

1435 if b < a:

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

1437 elif b > a:

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

1439 return Radius(Rbiaxial=b)

1440

1441 Requatorial = a # for consistent naming

1442

1443 def Rgeocentric(self, lat):

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

1445

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

1447

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

1449

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

1451

1452 @see: U{Geocentric Radius

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

1454 '''

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

1456 if a and self.f:

1457 if fabs(a) < PI_2:

1458 s2, c2 = _s2_c2(a)

1459 b2_a2_s2 = self.b2_a2 * s2

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

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

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

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

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

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

1466 else:

1467 r = self.b

1468 return Radius(Rgeocentric=r)

1469

1470 @Property_RO

1471 def Rgeometric(self):

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

1473

1474 @see: C{R1}.

1475 '''

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

1477 return Radius(Rgeometric=g)

1478

1479 def rhumb_(self, exact=True):

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

1481

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

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

1484 for I{this} ellipsoid.

1485

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

1487

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

1489

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

1491 '''

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

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

1494 else:

1495 r = exact

1496 E = _xattr(r, ellipsoid=None)

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

1498 raise _ValueError(exact=r, ellipsosid=E, txt_not_=self.name)

1499 return r

1500

1501 @property_RO

1502 def rhumbaux(self):

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

1504 '''

1505 # if not self.isEllipsoidal:

1506 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1508

1509 @property_RO

1510 def rhumbekx(self):

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

1512 '''

1513 # if not self.isEllipsoidal:

1514 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1516

1517 @property_RO

1518 def _Rhumbs(self):

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

1520 '''

1521 p = _MODS.rhumb

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

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

1524 return t

1525

1526 @property

1527 def rhumbsolve(self):

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

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

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

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

1532 '''

1533 # if not self.isEllipsoidal:

1534 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)

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

1536

1537 @rhumbsolve.setter # PYCHOK setter!

1538 def rhumbsolve(self, path):

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

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

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

1542 '''

1543 self._rhumbsolve = path

1544

1545 @deprecated_property_RO

1546 def rhumbx(self):

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

1548 return self.rhumbekx

1549

1550 def Rlat(self, lat):

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

1552

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

1554

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

1556

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

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

1559

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

1561

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

1563

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

1565

1566 @see: Method C{circle4}.

1567 '''

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

1569 r = self.a

1570 if self.f and lat: # .isEllipsoidal

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

1572 r = Radius(Rlat=r)

1573 return r

1574

1575 Rpolar = b # for consistent naming

1576

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

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

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

1580

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

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

1583 to use an alternate formula.

1584

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

1586

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

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

1589

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

1591 '''

1592 if not self.f: # .isSpherical

1593 n = self.a

1594 elif ca is None:

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

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

1597 elif ca: # derived from EcefYou.forward

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

1599 n = self.a / h

1600 elif sa:

1601 n = self.a2_b / fabs(sa)

1602 else:

1603 n = self.a

1604 return n

1605

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

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

1608 radii of curvature of (geodetic) latitude.

1609

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

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

1612

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

1614 the radii of curvature.

1615

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

1617

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

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

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

1621 Earth_radius#Radii_of_curvature>}.

1622 '''

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

1624

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

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

1627 curvature of (geodetic) latitude.

1628

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

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

1631

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

1633 radii of curvature.

1634

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

1636

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

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

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

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

1641 '''

1642 a = fabs(Phi(phi))

1643 if self.f:

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

1645 if r > EPS02:

1646 n = self.a / sqrt(r)

1647 m = n * self.e21 / r

1648 else:

1649 m = n = _0_0

1650 else:

1651 m = n = self.a

1652 if scaled and a:

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

1654 return Curvature2Tuple(Radius(rocMeridional=m),

1655 Radius(rocPrimeVertical=n))

1656

1657 def rocBearing(self, lat, bearing):

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

1659 latitude and compass direction.

1660

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

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

1663

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

1665

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

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

1668

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

1670

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

1672 '''

1673 if self.f:

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

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

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

1677 c2 *= n / m

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

1679 s2 *= m / n

1680 n = m

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

1682 else:

1683 b = self.b # == self.a

1684 return Radius(rocBearing=b)

1685

1686 @Property_RO

1687 def rocEquatorial2(self):

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

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

1690

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

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

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

1694 '''

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

1696 Radius(rocPrimeVertical0=self.a))

1697

1698 def rocGauss(self, lat):

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

1700

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

1702

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

1704

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

1706

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

1708 Earth_radius#Radii_of_curvature>}

1709 '''

1710 # using ...

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

1712 # return sqrt(m * n)

1713 # ... requires 1 or 2 sqrt

1714 g = self.b

1715 if self.f:

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

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

1718 return Radius(rocGauss=g)

1719

1720 def rocMean(self, lat):

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

1722

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

1724

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

1726

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

1728

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

1730 Earth_radius#Radii_of_curvature>}

1731 '''

1732 if self.f:

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

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

1735 else:

1736 m = self.a

1737 return Radius(rocMean=m)

1738

1739 def rocMeridional(self, lat):

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

1741

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

1743

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

1745

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

1747

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

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

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

1751 '''

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

1753 self.rocEquatorial2.meridional

1754

1755 rocPolar = a2_b # synonymous

1756

1757 def rocPrimeVertical(self, lat):

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

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

1760

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

1762

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

1764

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

1766

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

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

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

1770 Earth_radius#Radii_of_curvature>}.

1771 '''

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

1773 self.rocEquatorial2.prime_vertical

1774

1775 rocTransverse = rocPrimeVertical # synonymous

1776

1777 @deprecated_Property_RO

1778 def Rquadratic(self): # PYCHOK no cover

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

1780 return self.Rbiaxial

1781

1782 @deprecated_Property_RO

1783 def Rr(self): # PYCHOK no cover

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

1785 return self.Rrectifying

1786

1787 @Property_RO

1788 def Rrectifying(self):

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

1790

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

1792 '''

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

1794 return Radius(Rrectifying=r)

1795

1796 @deprecated_Property_RO

1797 def Rs(self): # PYCHOK no cover

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

1799 return self.Rgeometric

1800

1801 @Property_RO

1802 def Rtriaxial(self):

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

1804

1805 @see: C{Rbiaxial}

1806 '''

1807 a, b = self.a, self.b

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

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

1810 return Radius(Rtriaxial=q)

1811

1812 def toEllipsoid2(self, **name):

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

1814

1815 @kwarg name: Optional, unique C{B{name}=NN} (C{str}).

1816

1817 @see: Property C{a_f}.

1818 '''

1819 return Ellipsoid2(self, None, **name)

1820

1821 def toStr(self, prec=8, terse=0, **name): # PYCHOK expected

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

1823

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

1825 @kwarg terse: Limit the number of items (C{int}, 0...18).

1826 @kwarg name: Optional C{B{name}=NN} (C{str}) or C{None} to

1827 exclude this ellipsoid's name.

1828

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

1830 '''

1831 E = Ellipsoid

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

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

1834 if terse:

1835 t = t[:terse]

1836 n, _ = _name2__(**name) # name=None

1837 return self._instr(n, prec, props=t)

1838

1839 def toTriaxial(self, **name):

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

1841

1842 @kwarg name: Optional C{B{name}=NN} (C{str}).

1843

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

1845 pointing east and C{Z} pointing north.

1846

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

1848 '''

1849 T = self._triaxial

1850 return T.copy(**name) if name else T

1851

1852 @property_RO

1853 def _triaxial(self):

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

1855 '''

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

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

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

1859

1860 @Property_RO

1861 def volume(self):

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

1863

1864 @see: C{R3}.

1865 '''

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

1867

1868

1869 class Ellipsoid2(Ellipsoid):

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

1871 '''

1872 def __init__(self, a, f=None, **name):

1873 '''New L{Ellipsoid2}.

1874

1875 @arg a: Equatorial radius, semi-axis (C{meter}) or a previous

1876 L{Ellipsoid} instance.

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

1878 if B{C{a}} is in C{meter}.

1879 @kwarg name: Optional, unique C{B{name}=NN} (C{str}).

1880

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

1882

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

1884

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

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

1887 '''

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

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

1890 b=a.b, f_=a.f_, **name)

1891 else:

1892 Ellipsoid.__init__(self, a, f=f, **name)

1893

1894

1895 def _spherical_a_b(a, b):

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

1897 '''

1898 return a < EPS0 or b < EPS0 or fabs(a - b) < EPS0

1899

1900

1901 def _spherical_f(f):

1902 '''(INTERNAL) C{True} for spherical or invalid C{f}.

1903 '''

1904 return fabs(f) < EPS or f > EPS1

1905

1906

1907 def _spherical_f_(f_):

1908 '''(INTERNAL) C{True} for spherical or invalid C{f_}.

1909 '''

1910 return fabs(f_) < EPS or fabs(f_) > _1_EPS

1911

1912

1913 def a_b2e(a, b):

1914 '''Return C{e}, the I{1st eccentricity} for a given I{equatorial} and I{polar} radius.

1915

1916 @arg a: Equatorial radius (C{scalar} > 0).

1917 @arg b: Polar radius (C{scalar} > 0).

1918

1919 @return: The I{unsigned}, (1st) eccentricity (C{float} or C{0}),

1920 M{sqrt(1 - (b / a)**2)}.

1921

1922 @note: The result is always I{non-negative} and C{0} for I{near-spherical} ellipsoids.

1923 '''

1924 return Float(e=sqrt(fabs(a_b2e2(a, b)))) # == sqrt(fabs(a - b) * (a + b)) / a)

1925

1926

1927 def a_b2e2(a, b):

1928 '''Return C{e2}, the I{1st eccentricity squared} for a given I{equatorial} and I{polar} radius.

1929

1930 @arg a: Equatorial radius (C{scalar} > 0).

1931 @arg b: Polar radius (C{scalar} > 0).

1932

1933 @return: The I{signed}, (1st) eccentricity I{squared} (C{float} or C{0}),

1934 M{1 - (b / a)**2}.

1935

1936 @note: The result is positive for I{oblate}, negative for I{prolate}

1937 or C{0} for I{near-spherical} ellipsoids.

1938 '''

1939 return Float(e2=_0_0 if _spherical_a_b(a, b) else ((a - b) * (a + b) / a**2))

1940

1941

1942 def a_b2e22(a, b):

1943 '''Return C{e22}, the I{2nd eccentricity squared} for a given I{equatorial} and I{polar} radius.

1944

1945 @arg a: Equatorial radius (C{scalar} > 0).

1946 @arg b: Polar radius (C{scalar} > 0).

1947

1948 @return: The I{signed}, 2nd eccentricity I{squared} (C{float} or C{0}),

1949 M{(a / b)**2 - 1}.

1950

1951 @note: The result is positive for I{oblate}, negative for I{prolate}

1952 or C{0} for I{near-spherical} ellipsoids.

1953 '''

1954 return Float(e22=_0_0 if _spherical_a_b(a, b) else ((a - b) * (a + b) / b**2))

1955

1956

1957 def a_b2e32(a, b):

1958 '''Return C{e32}, the I{3rd eccentricity squared} for a given I{equatorial} and I{polar} radius.

1959

1960 @arg a: Equatorial radius (C{scalar} > 0).

1961 @arg b: Polar radius (C{scalar} > 0).

1962

1963 @return: The I{signed}, 3rd eccentricity I{squared} (C{float} or C{0}),

1964 M{(a**2 - b**2) / (a**2 + b**2)}.

1965

1966 @note: The result is positive for I{oblate}, negative for I{prolate}

1967 or C{0} for I{near-spherical} ellipsoids.

1968 '''

1969 a2, b2 = a**2, b**2

1970 return Float(e32=_0_0 if _spherical_a_b(a2, b2) else ((a2 - b2) / (a2 + b2)))

1971

1972

1973 def a_b2f(a, b):

1974 '''Return C{f}, the I{flattening} for a given I{equatorial} and I{polar} radius.

1975

1976 @arg a: Equatorial radius (C{scalar} > 0).

1977 @arg b: Polar radius (C{scalar} > 0).

1978

1979 @return: The flattening (C{scalar} or C{0}), M{(a - b) / a}.

1980

1981 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0}

1982 for I{near-spherical} ellipsoids.

1983 '''

1984 f = 0 if _spherical_a_b(a, b) else ((a - b) / a)

1985 return _f_0_0 if _spherical_f(f) else Float(f=f)

1986

1987

1988 def a_b2f_(a, b):

1989 '''Return C{f_}, the I{inverse flattening} for a given I{equatorial} and I{polar} radius.

1990

1991 @arg a: Equatorial radius (C{scalar} > 0).

1992 @arg b: Polar radius (C{scalar} > 0).

1993

1994 @return: The inverse flattening (C{scalar} or C{0}), M{a / (a - b)}.

1995

1996 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0}

1997 for I{near-spherical} ellipsoids.

1998 '''

1999 f_ = 0 if _spherical_a_b(a, b) else (a / float(a - b))

2000 return _f__0_0 if _spherical_f_(f_) else Float(f_=f_)

2001

2002

2003 def a_b2f2(a, b):

2004 '''Return C{f2}, the I{2nd flattening} for a given I{equatorial} and I{polar} radius.

2005

2006 @arg a: Equatorial radius (C{scalar} > 0).

2007 @arg b: Polar radius (C{scalar} > 0).

2008

2009 @return: The I{signed}, 2nd flattening (C{scalar} or C{0}), M{(a - b) / b}.

2010

2011 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0}

2012 for I{near-spherical} ellipsoids.

2013 '''

2014 t = 0 if _spherical_a_b(a, b) else float(a - b)

2015 return Float(f2=_0_0 if fabs(t) < EPS0 else (t / b))

2016

2017

2018 def a_b2n(a, b):

2019 '''Return C{n}, the I{3rd flattening} for a given I{equatorial} and I{polar} radius.

2020

2021 @arg a: Equatorial radius (C{scalar} > 0).

2022 @arg b: Polar radius (C{scalar} > 0).

2023

2024 @return: The I{signed}, 3rd flattening (C{scalar} or C{0}), M{(a - b) / (a + b)}.

2025

2026 @note: The result is positive for I{oblate}, negative for I{prolate}

2027 or C{0} for I{near-spherical} ellipsoids.

2028 '''

2029 t = 0 if _spherical_a_b(a, b) else float(a - b)

2030 return Float(n=_0_0 if fabs(t) < EPS0 else (t / (a + b)))

2031

2032

2033 def a_f2b(a, f):

2034 '''Return C{b}, the I{polar} radius for a given I{equatorial} radius and I{flattening}.

2035

2036 @arg a: Equatorial radius (C{scalar} > 0).

2037 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2038

2039 @return: The polar radius (C{float}), M{a * (1 - f)}.

2040 '''

2041 b = a if _spherical_f(f) else (a * (_1_0 - f))

2042 return Radius_(b=a if _spherical_a_b(a, b) else b)

2043

2044

2045 def a_f_2b(a, f_):

2046 '''Return C{b}, the I{polar} radius for a given I{equatorial} radius and I{inverse flattening}.

2047

2048 @arg a: Equatorial radius (C{scalar} > 0).

2049 @arg f_: Inverse flattening (C{scalar} >>> 1).

2050

2051 @return: The polar radius (C{float}), M{a * (f_ - 1) / f_}.

2052 '''

2053 b = a if _spherical_f_(f_) else (a * (f_ - _1_0) / f_)

2054 return Radius_(b=a if _spherical_a_b(a, b) else b)

2055

2056

2057 def b_f2a(b, f):

2058 '''Return C{a}, the I{equatorial} radius for a given I{polar} radius and I{flattening}.

2059

2060 @arg b: Polar radius (C{scalar} > 0).

2061 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2062

2063 @return: The equatorial radius (C{float}), M{b / (1 - f)}.

2064 '''

2065 t = _1_0 - f

2066 a = b if fabs(t) < EPS0 else (b / t)

2067 return Radius_(a=b if _spherical_a_b(a, b) else a)

2068

2069

2070 def b_f_2a(b, f_):

2071 '''Return C{a}, the I{equatorial} radius for a given I{polar} radius and I{inverse flattening}.

2072

2073 @arg b: Polar radius (C{scalar} > 0).

2074 @arg f_: Inverse flattening (C{scalar} >>> 1).

2075

2076 @return: The equatorial radius (C{float}), M{b * f_ / (f_ - 1)}.

2077 '''

2078 t = f_ - _1_0

2079 a = b if _spherical_f_(f_) or fabs(t - f_) < EPS0 \

2080 or fabs(t) < EPS0 else (b * f_ / t)

2081 return Radius_(a=b if _spherical_a_b(a, b) else a)

2082

2083

2084 def e2f(e):

2085 '''Return C{f}, the I{flattening} for a given I{1st eccentricity}.

2086

2087 @arg e: The (1st) eccentricity (0 <= C{float} < 1)

2088

2089 @return: The flattening (C{scalar} or C{0}).

2090

2091 @see: Function L{e22f}.

2092 '''

2093 return e22f(e**2)

2094

2095

2096 def e22f(e2):

2097 '''Return C{f}, the I{flattening} for a given I{1st eccentricity squared}.

2098

2099 @arg e2: The (1st) eccentricity I{squared}, I{signed} (L{NINF} < C{float} < 1)

2100

2101 @return: The flattening (C{float} or C{0}), M{e2 / (sqrt(e2 - 1) + 1)}.

2102 '''

2103 return Float(f=e2 / (sqrt(_1_0 - e2) + _1_0)) if e2 else _f_0_0

2104

2105

2106 def f2e2(f):

2107 '''Return C{e2}, the I{1st eccentricity squared} for a given I{flattening}.

2108

2109 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2110

2111 @return: The I{signed}, (1st) eccentricity I{squared} (C{float} < 1),

2112 M{f * (2 - f)}.

2113

2114 @note: The result is positive for I{oblate}, negative for I{prolate}

2115 or C{0} for I{near-spherical} ellipsoids.

2116

2117 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2118 C++/doc/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2119 <https://WikiPedia.org/wiki/Flattening>}.

2120 '''

2121 return Float(e2=_0_0 if _spherical_f(f) else (f * (_2_0 - f)))

2122

2123

2124 def f2e22(f):

2125 '''Return C{e22}, the I{2nd eccentricity squared} for a given I{flattening}.

2126

2127 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2128

2129 @return: The I{signed}, 2nd eccentricity I{squared} (C{float} > -1 or

2130 C{INF}), M{f * (2 - f) / (1 - f)**2}.

2131

2132 @note: The result is positive for I{oblate}, negative for I{prolate}

2133 or C{0} for near-spherical ellipsoids.

2134

2135 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2136 C++/doc/classGeographicLib_1_1Ellipsoid.html>}.

2137 '''

2138 # e2 / (1 - e2) == f * (2 - f) / (1 - f)**2

2139 t = (_1_0 - f)**2

2140 return Float(e22=INF if t < EPS0 else (f2e2(f) / t)) # PYCHOK type

2141

2142

2143 def f2e32(f):

2144 '''Return C{e32}, the I{3rd eccentricity squared} for a given I{flattening}.

2145

2146 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2147

2148 @return: The I{signed}, 3rd eccentricity I{squared} (C{float}),

2149 M{f * (2 - f) / (1 + (1 - f)**2)}.

2150

2151 @note: The result is positive for I{oblate}, negative for I{prolate}

2152 or C{0} for I{near-spherical} ellipsoids.

2153

2154 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2155 C++/doc/classGeographicLib_1_1Ellipsoid.html>}.

2156 '''

2157 # e2 / (2 - e2) == f * (2 - f) / (1 + (1 - f)**2)

2158 e2 = f2e2(f)

2159 return Float(e32=e2 / (_2_0 - e2))

2160

2161

2162 def f_2f(f_):

2163 '''Return C{f}, the I{flattening} for a given I{inverse flattening}.

2164

2165 @arg f_: Inverse flattening (C{scalar} >>> 1).

2166

2167 @return: The flattening (C{scalar} or C{0}), M{1 / f_}.

2168

2169 @note: The result is positive for I{oblate}, negative for I{prolate}

2170 or C{0} for I{near-spherical} ellipsoids.

2171 '''

2172 f = 0 if _spherical_f_(f_) else _1_0 / f_

2173 return _f_0_0 if _spherical_f(f) else Float(f=f) # PYCHOK type

2174

2175

2176 def f2f_(f):

2177 '''Return C{f_}, the I{inverse flattening} for a given I{flattening}.

2178

2179 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2180

2181 @return: The inverse flattening (C{scalar} or C{0}), M{1 / f}.

2182

2183 @note: The result is positive for I{oblate}, negative for I{prolate}

2184 or C{0} for I{near-spherical} ellipsoids.

2185 '''

2186 f_ = 0 if _spherical_f(f) else _1_0 / f

2187 return _f__0_0 if _spherical_f_(f_) else Float(f_=f_) # PYCHOK type

2188

2189

2190 def f2f2(f):

2191 '''Return C{f2}, the I{2nd flattening} for a given I{flattening}.

2192

2193 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2194

2195 @return: The I{signed}, 2nd flattening (C{scalar} or C{INF}), M{f / (1 - 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 C++/doc/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2202 <https://WikiPedia.org/wiki/Flattening>}.

2203 '''

2204 t = _1_0 - f

2205 return Float(f2=_0_0 if _spherical_f(f) else (INF if fabs(t) < EPS

2206 else (f / t))) # PYCHOK type

2207

2208

2209 def f2n(f):

2210 '''Return C{n}, the I{3rd flattening} for a given I{flattening}.

2211

2212 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).

2213

2214 @return: The I{signed}, 3rd flattening (-1 <= C{float} < 1),

2215 M{f / (2 - f)}.

2216

2217 @note: The result is positive for I{oblate}, negative for I{prolate}

2218 or C{0} for I{near-spherical} ellipsoids.

2219

2220 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2221 C++/doc/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2222 <https://WikiPedia.org/wiki/Flattening>}.

2223 '''

2224 return Float(n=_0_0 if _spherical_f(f) else (f / float(_2_0 - f)))

2225

2226

2227 def n2e2(n):

2228 '''Return C{e2}, the I{1st eccentricity squared} for a given I{3rd flattening}.

2229

2230 @arg n: The 3rd flattening (-1 <= C{scalar} < 1).

2231

2232 @return: The I{signed}, (1st) eccentricity I{squared} (C{float} or NINF),

2233 M{4 * n / (1 + n)**2}.

2234

2235 @note: The result is positive for I{oblate}, negative for I{prolate}

2236 or C{0} for I{near-spherical} ellipsoids.

2237

2238 @see: U{Flattening<https://WikiPedia.org/wiki/Flattening>}.

2239 '''

2240 t = (n + _1_0)**2

2241 return Float(e2=_0_0 if fabs(n) < EPS0 else

2242 (NINF if t < EPS0 else (_4_0 * n / t)))

2243

2244

2245 def n2f(n):

2246 '''Return C{f}, the I{flattening} for a given I{3rd flattening}.

2247

2248 @arg n: The 3rd flattening (-1 <= C{scalar} < 1).

2249

2250 @return: The flattening (C{scalar} or NINF), M{2 * n / (1 + n)}.

2251

2252 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/

2253 C++/doc/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening

2254 <https://WikiPedia.org/wiki/Flattening>}.

2255 '''

2256 t = n + _1_0

2257 f = 0 if fabs(n) < EPS0 else (NINF if t < EPS0 else (_2_0 * n / t))

2258 return _f_0_0 if _spherical_f(f) else Float(f=f)

2259

2260

2261 def n2f_(n):

2262 '''Return C{f_}, the I{inverse flattening} for a given I{3rd flattening}.

2263

2264 @arg n: The 3rd flattening (-1 <= C{scalar} < 1).

2265

2266 @return: The inverse flattening (C{scalar} or C{0}), M{1 / f}.

2267

2268 @see: L{n2f} and L{f2f_}.

2269 '''

2270 return f2f_(n2f(n))

2271

2272

2273 def _normalTo3(px, py, E, eps=EPS): # in .height4 above

2274 '''(INTERNAL) Nearest point on a 2-D ellipse in 1st quadrant.

2275

2276 @see: Functions C{pygeodesy.triaxial._normalTo4} and C{-To5}.

2277 '''

2278 a, b, e0 = E.a, E.b, EPS0

2279 if min(px, py, a, b) < e0:

2280 raise _AssertionError(px=px, py=py, a=a, b=b, E=E)

2281

2282 a2 = a - b * E.b_a

2283 b2 = b - a * E.a_b

2284 tx = ty = _SQRT2_2

2285 _a, _h = fabs, hypot

2286 for i in range(16): # max 5

2287 ex = a2 * tx**3

2288 ey = b2 * ty**3

2289

2290 qx = px - ex

2291 qy = py - ey

2292 q = _h(qx, qy)

2293 if q < e0: # PYCHOK no cover

2294 break

2295 r = _h(ex - tx * a,

2296 ey - ty * b) / q

2297

2298 sx, tx = tx, min(_1_0, max(0, (ex + qx * r) / a))

2299 sy, ty = ty, min(_1_0, max(0, (ey + qy * r) / b))

2300 t = _h(ty, tx)

2301 if t < e0: # PYCHOK no cover

2302 break

2303 tx = tx / t # /= chokes PyChecker

2304 ty = ty / t

2305 if _a(sx - tx) < eps and \

2306 _a(sy - ty) < eps:

2307 break

2308

2309 tx *= a / px

2310 ty *= b / py

2311 return tx, ty, i # x and y as fractions

2312

2313

2314 class Ellipsoids(_NamedEnum):

2315 '''(INTERNAL) L{Ellipsoid} registry, I{must} be a sub-class

2316 to accommodate the L{_LazyNamedEnumItem} properties.

2317 '''

2318 def _Lazy(self, a, b, f_, **kwds):

2319 '''(INTERNAL) Instantiate the L{Ellipsoid}.

2320 '''

2321 return Ellipsoid(a, b=b, f_=f_, **kwds)

2322

2323 Ellipsoids = Ellipsoids(Ellipsoid) # PYCHOK singleton

2324'''Some pre-defined L{Ellipsoid}s, all I{lazily} instantiated.'''

2325# <https://www.GNU.org/software/gama/manual/html_node/Supported-ellipsoids.html>

2326# <https://GSSC.ESA.int/navipedia/index.php/Reference_Frames_in_GNSS>

2327# <https://kb.OSU.edu/dspace/handle/1811/77986>

2328# <https://www.IBM.com/docs/en/db2/11.5?topic=systems-supported-spheroids>

2329# <https://w3.Energistics.org/archive/Epicentre/Epicentre_v3.0/DataModel/LogicalDictionary/StandardValues/ellipsoid.html>

2330# <https://GitHub.com/locationtech/proj4j/blob/master/src/main/java/org/locationtech/proj4j/datum/Ellipsoid.java>

2331 Ellipsoids._assert( # <https://WikiPedia.org/wiki/Earth_ellipsoid>

2332 Airy1830 = _lazy(_Airy1830_, *_T(6377563.396, _0_0, 299.3249646)), # b=6356256.909

2333 AiryModified = _lazy(_AiryModified_, *_T(6377340.189, _0_0, 299.3249646)), # b=6356034.448

2334# APL4_9 = _lazy('APL4_9', *_T(6378137.0, _0_0, 298.24985392)), # Appl. Phys. Lab. 1965

2335# ANS = _lazy('ANS', *_T(6378160.0, _0_0, 298.25)), # Australian Nat. Spheroid

2336# AN_SA96 = _lazy('AN_SA96', *_T(6378160.0, _0_0, 298.24985392)), # Australian Nat. South America

2337 Australia1966 = _lazy('Australia1966', *_T(6378160.0, _0_0, 298.25)), # b=6356774.7192

2338 ATS1977 = _lazy('ATS1977', *_T(6378135.0, _0_0, 298.257)), # "Average Terrestrial System"

2339 Bessel1841 = _lazy(_Bessel1841_, *_T(6377397.155, 6356078.962818, 299.152812797)),

2340 BesselModified = _lazy('BesselModified', *_T(6377492.018, _0_0, 299.1528128)),

2341# BesselNamibia = _lazy('BesselNamibia', *_T(6377483.865, _0_0, 299.1528128)),

2342 CGCS2000 = _lazy('CGCS2000', *_T(R_MA, _0_0, 298.257222101)), # BeiDou Coord System (BDC)

2343# Clarke1858 = _lazy('Clarke1858', *_T(6378293.639, _0_0, 294.260676369)),

2344 Clarke1866 = _lazy(_Clarke1866_, *_T(6378206.4, 6356583.8, 294.978698214)),

2345 Clarke1880 = _lazy('Clarke1880', *_T(6378249.145, 6356514.86954978, 293.465)),

2346 Clarke1880IGN = _lazy(_Clarke1880IGN_, *_T(6378249.2, 6356515.0, 293.466021294)),

2347 Clarke1880Mod = _lazy('Clarke1880Mod', *_T(6378249.145, 6356514.96639549, 293.466307656)), # aka Clarke1880Arc

2348 CPM1799 = _lazy('CPM1799', *_T(6375738.7, 6356671.92557493, 334.39)), # Comm. des Poids et Mesures

2349 Delambre1810 = _lazy('Delambre1810', *_T(6376428.0, 6355957.92616372, 311.5)), # Belgium

2350 Engelis1985 = _lazy('Engelis1985', *_T(6378136.05, 6356751.32272154, 298.2566)),

2351# Everest1830 = _lazy('Everest1830', *_T(6377276.345, _0_0, 300.801699997)),

2352# Everest1948 = _lazy('Everest1948', *_T(6377304.063, _0_0, 300.801699997)),

2353# Everest1956 = _lazy('Everest1956', *_T(6377301.243, _0_0, 300.801699997)),

2354 Everest1969 = _lazy('Everest1969', *_T(6377295.664, 6356094.667915, 300.801699997)),

2355 Everest1975 = _lazy('Everest1975', *_T(6377299.151, 6356098.14512013, 300.8017255)),

2356 Fisher1968 = _lazy('Fisher1968', *_T(6378150.0, 6356768.33724438, 298.3)),

2357# Fisher1968Mod = _lazy('Fisher1968Mod', *_T(6378155.0, _0_0, 298.3)),

2358 GEM10C = _lazy('GEM10C', *_T(R_MA, 6356752.31424783, 298.2572236)),

2359 GPES = _lazy('GPES', *_T(6378135.0, 6356750.0, _0_0)), # "Gen. Purpose Earth Spheroid"

2360 GRS67 = _lazy('GRS67', *_T(6378160.0, _0_0, 298.247167427)), # Lucerne b=6356774.516

2361# GRS67Truncated = _lazy('GRS67Truncated', *_T(6378160.0, _0_0, 298.25)),

2362 GRS80 = _lazy(_GRS80_, *_T(R_MA, 6356752.314140347, 298.25722210088)), # IUGG, ITRS, ETRS89

2363# Hayford1924 = _lazy('Hayford1924', *_T(6378388.0, 6356911.94612795, None)), # aka Intl1924 f_=297

2364 Helmert1906 = _lazy('Helmert1906', *_T(6378200.0, 6356818.16962789, 298.3)),

2365# Hough1960 = _lazy('Hough1960', *_T(6378270.0, _0_0, 297.0)),

2366 IAU76 = _lazy('IAU76', *_T(6378140.0, _0_0, 298.257)), # Int'l Astronomical Union

2367 IERS1989 = _lazy('IERS1989', *_T(6378136.0, _0_0, 298.257)), # b=6356751.302

2368 IERS1992TOPEX = _lazy('IERS1992TOPEX', *_T(6378136.3, 6356751.61659215, 298.257223563)), # IERS/TOPEX/Poseidon/McCarthy

2369 IERS2003 = _lazy('IERS2003', *_T(6378136.6, 6356751.85797165, 298.25642)),

2370 Intl1924 = _lazy(_Intl1924_, *_T(6378388.0, _0_0, 297.0)), # aka Hayford b=6356911.9462795

2371 Intl1967 = _lazy('Intl1967', *_T(6378157.5, 6356772.2, 298.24961539)), # New Int'l

2372 Krassovski1940 = _lazy(_Krassovski1940_, *_T(6378245.0, 6356863.01877305, 298.3)), # spelling

2373 Krassowsky1940 = _lazy(_Krassowsky1940_, *_T(6378245.0, 6356863.01877305, 298.3)), # spelling

2374# Kaula = _lazy('Kaula', *_T(6378163.0, _0_0, 298.24)), # Kaula 1961

2375# Lerch = _lazy('Lerch', *_T(6378139.0, _0_0, 298.257)), # Lerch 1979

2376 Maupertuis1738 = _lazy('Maupertuis1738', *_T(6397300.0, 6363806.28272251, 191.0)), # France

2377 Mercury1960 = _lazy('Mercury1960', *_T(6378166.0, 6356784.28360711, 298.3)),

2378 Mercury1968Mod = _lazy('Mercury1968Mod', *_T(6378150.0, 6356768.33724438, 298.3)),

2379# MERIT = _lazy('MERIT', *_T(6378137.0, _0_0, 298.257)), # MERIT 1983

2380# NWL10D = _lazy('NWL10D', *_T(6378135.0, _0_0, 298.26)), # Naval Weapons Lab.

2381 NWL1965 = _lazy('NWL1965', *_T(6378145.0, 6356759.76948868, 298.25)), # Naval Weapons Lab.

2382# NWL9D = _lazy('NWL9D', *_T(6378145.0, 6356759.76948868, 298.25)), # NWL1965

2383 OSU86F = _lazy('OSU86F', *_T(6378136.2, 6356751.51693008, 298.2572236)),

2384 OSU91A = _lazy('OSU91A', *_T(6378136.3, 6356751.6165948, 298.2572236)),

2385# Plessis1817 = _lazy('Plessis1817', *_T(6397523.0, 6355863.0, 153.56512242)), # XXX incorrect?

2386 Plessis1817 = _lazy('Plessis1817', *_T(6376523.0, 6355862.93325557, 308.64)), # XXX IGN France 1972

2387# Prolate = _lazy('Prolate', *_T(6356752.3, R_MA, _0_0)),

2388 PZ90 = _lazy('PZ90', *_T(6378136.0, _0_0, 298.257839303)), # GLOSNASS PZ-90 and PZ-90.11

2389# SEAsia = _lazy('SEAsia', *_T(6378155.0, _0_0, 298.3)), # SouthEast Asia

2390 SGS85 = _lazy('SGS85', *_T(6378136.0, 6356751.30156878, 298.257)), # Soviet Geodetic System

2391 SoAmerican1969 = _lazy('SoAmerican1969', *_T(6378160.0, 6356774.71919531, 298.25)), # South American

2392 Sphere = _lazy(_Sphere_, *_T(R_M, R_M, _0_0)), # pseudo

2393 SphereAuthalic = _lazy('SphereAuthalic', *_T(R_FM, R_FM, _0_0)), # pseudo

2394 SpherePopular = _lazy('SpherePopular', *_T(R_MA, R_MA, _0_0)), # EPSG:3857 Spheroid

2395 Struve1860 = _lazy('Struve1860', *_T(6378298.3, 6356657.14266956, 294.73)),

2396# Walbeck = _lazy('Walbeck', *_T(6376896.0, _0_0, 302.78)),

2397# WarOffice = _lazy('WarOffice', *_T(6378300.0, _0_0, 296.0)),

2398 WGS60 = _lazy('WGS60', *_T(6378165.0, 6356783.28695944, 298.3)),

2399 WGS66 = _lazy('WGS66', *_T(6378145.0, 6356759.76948868, 298.25)),

2400 WGS72 = _lazy(_WGS72_, *_T(6378135.0, _0_0, 298.26)), # b=6356750.52

2401 WGS84 = _lazy(_WGS84_, *_T(R_MA, _0_0, _f__WGS84)), # GPS b=6356752.3142451793

2402# U{NOAA/NOS/NGS/inverse<https://GitHub.com/noaa-ngs/inverse/blob/main/invers3d.f>}

2403 WGS84_NGS = _lazy('WGS84_NGS', *_T(R_MA, _0_0, 298.257222100882711243162836600094))

2404)

2405

2406_EWGS84 = Ellipsoids.WGS84 # (INTERNAL) shared

2407

2408 if __name__ == '__main__':

2409

2410 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_

2411 from pygeodesy import nameof, printf

2412

2413 for E in (_EWGS84, Ellipsoids.GRS80, # NAD83,

2414 Ellipsoids.Sphere, Ellipsoids.SpherePopular,

2415 Ellipsoid(_EWGS84.b, _EWGS84.a, name='_Prolate')):

2416 e = f2n(E.f) - E.n

2417 printf('# %s: %s', _DOT_('Ellipsoids', E.name), E.toStr(prec=10), nl=1)

2418 printf('# e=%s, f_=%s, f=%s, n=%s (%s)', fstr(E.e, prec=13, fmt=Fmt.e),

2419 fstr(E.f_, prec=13, fmt=Fmt.e),

2420 fstr(E.f, prec=13, fmt=Fmt.e),

2421 fstr(E.n, prec=13, fmt=Fmt.e),

2422 fstr(e, prec=9, fmt=Fmt.e))

2423 printf('# %s %s', Ellipsoid.AlphaKs.name, fstr(E.AlphaKs, prec=20))

2424 printf('# %s %s', Ellipsoid.BetaKs.name, fstr(E.BetaKs, prec=20))

2425 printf('# %s %s', nameof(Ellipsoid.KsOrder), E.KsOrder) # property

2426

2427 # __doc__ of this file, force all into registry

2428 t = [NN] + Ellipsoids.toRepr(all=True, asorted=True).split(_NL_)

2429 printf(_NLATvar_.join(i.strip(_COMMA_) for i in t))

2430

2431# % python3 -m pygeodesy.ellipsoids

2432

2433# Ellipsoids.WGS84: name='WGS84', a=6378137, b=6356752.3142451793, f_=298.257223563, f=0.0033528107, f2=0.0033640898, n=0.0016792204, e=0.0818191908, e2=0.00669438, e21=0.99330562, e22=0.0067394967, e32=0.0033584313, A=6367449.1458234144, L=10001965.7293127235, R1=6371008.7714150595, R2=6371007.1809184738, R3=6371000.7900091587, Rbiaxial=6367453.6345163295, Rtriaxial=6372797.5559594007

2434# e=8.1819190842622e-02, f_=2.98257223563e+02, f=3.3528106647475e-03, n=1.6792203863837e-03 (0.0e+00)

2435# AlphaKs 0.00083773182062446994, 0.00000076085277735725, 0.00000000119764550324, 0.00000000000242917068, 0.00000000000000571182, 0.0000000000000000148, 0.00000000000000000004, 0.0

2436# BetaKs 0.00083773216405794875, 0.0000000590587015222, 0.00000000016734826653, 0.00000000000021647981, 0.00000000000000037879, 0.00000000000000000072, 0.0, 0.0

2437# KsOrder 8

2438

2439# Ellipsoids.GRS80: name='GRS80', a=6378137, b=6356752.3141403468, f_=298.2572221009, f=0.0033528107, f2=0.0033640898, n=0.0016792204, e=0.081819191, e2=0.00669438, e21=0.99330562, e22=0.0067394968, e32=0.0033584313, A=6367449.1457710434, L=10001965.7292304561, R1=6371008.7713801153, R2=6371007.1808835147, R3=6371000.7899741363, Rbiaxial=6367453.6344640013, Rtriaxial=6372797.5559332585

2440# e=8.1819191042833e-02, f_=2.9825722210088e+02, f=3.3528106811837e-03, n=1.6792203946295e-03 (0.0e+00)

2441# AlphaKs 0.00083773182472890429, 0.00000076085278481561, 0.00000000119764552086, 0.00000000000242917073, 0.00000000000000571182, 0.0000000000000000148, 0.00000000000000000004, 0.0

2442# BetaKs 0.0008377321681623882, 0.00000005905870210374, 0.000000000167348269, 0.00000000000021647982, 0.00000000000000037879, 0.00000000000000000072, 0.0, 0.0

2443# KsOrder 8

2444

2445# 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

2446# e=0.0e+00, f_=0.0e+00, f=0.0e+00, n=0.0e+00 (0.0e+00)

2447# AlphaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2448# BetaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2449# KsOrder 8

2450

2451# 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

2452# e=0.0e+00, f_=0.0e+00, f=0.0e+00, n=0.0e+00 (0.0e+00)

2453# AlphaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2454# BetaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

2455# KsOrder 8

2456

2457# Ellipsoids._Prolate: name='_Prolate', a=6356752.3142451793, b=6378137, f_=-297.257223563, f=-0.0033640898, f2=-0.0033528107, n=-0.0016792204, e=0.0820944379, e2=-0.0067394967, e21=1.0067394967, e22=-0.00669438, e32=-0.0033584313, A=6367449.1458234144, L=10035500.5204500314, R1=6363880.5428301189, R2=6363878.9413582645, R3=6363872.5644020075, Rbiaxial=6367453.6345163295, Rtriaxial=6362105.2243882557

2458# e=8.2094437949696e-02, f_=-2.97257223563e+02, f=-3.3640898209765e-03, n=-1.6792203863837e-03 (0.0e+00)

2459# AlphaKs -0.00084149152514366627, 0.00000076653480614871, -0.00000000120934503389, 0.0000000000024576225, -0.00000000000000578863, 0.00000000000000001502, -0.00000000000000000004, 0.0

2460# BetaKs -0.00084149187224351817, 0.00000005842735196773, -0.0000000001680487236, 0.00000000000021706261, -0.00000000000000038002, 0.00000000000000000073, -0.0, 0.0

2461# KsOrder 8

2462

2463# **) MIT License

2464#

2465# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.

2466#

2467# Permission is hereby granted, free of charge, to any person obtaining a

2468# copy of this software and associated documentation files (the "Software"),

2469# to deal in the Software without restriction, including without limitation

2470# the rights to use, copy, modify, merge, publish, distribute, sublicense,

2471# and/or sell copies of the Software, and to permit persons to whom the

2472# Software is furnished to do so, subject to the following conditions:

2473#

2474# The above copyright notice and this permission notice shall be included

2475# in all copies or substantial portions of the Software.

2476#

2477# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS

2478# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

2479# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL

2480# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR

2481# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,

2482# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR

2483# OTHER DEALINGS IN THE SOFTWARE.