1 2# -*- coding: utf-8 -*- 3 4u'''Ellipsoidal and spherical earth models. 5 6Classes L{a_f2Tuple}, L{Ellipsoid} and L{Ellipsoid2}, an L{Ellipsoids} registry and 7a dozen functions to convert I{equatorial} radius, I{polar} radius, I{eccentricities}, 8I{flattenings} and I{inverse flattening}. 9 10See module L{datums} for more information and other details. 11 12@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, 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) 13@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, 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) 14@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, 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) 15@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, 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) 16@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, 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) 17@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, 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) 18@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, 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) 19@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, 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) 20@var Ellipsoids.Clarke1880Mod: Ellipsoid(name='Clarke1880Mod', a=6378249.145, b=6356514.96582849, f_=293.4663, f=0.00340755, f2=0.0034192, n=0.00170668, e=0.08248322, e2=0.00680348, e22=0.00685009, e32=0.00341335, A=6367386.69207875, L=10001867.62720001, R1=6371004.4186095, R2=6371002.77577708, R3=6370996.17408252, Rbiaxial=6367391.32873482, Rtriaxial=6372822.54926891) 21@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, 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) 22@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, 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) 23@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, 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) 24@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, 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) 25@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, 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) 26@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, 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) 27@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, 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) 28@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, 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) 29@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, 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) 30@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, 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) 31@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, 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) 32@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, 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) 33@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, 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) 34@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, 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) 35@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, 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) 36@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, 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) 37@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, 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) 38@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, 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) 39@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, 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) 40@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, 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) 41@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, e22=0.0067395, e32=0.00335843, A=6367448.44699773, L=10001964.6315999, R1=6371008.07219827, R2=6371006.48170186, R3=6371000.09079324, Rbiaxial=6367452.93569015, Rtriaxial=6372796.85654607) 42@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, 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) 43@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, 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) 44@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, 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) 45@var Ellipsoids.Sphere: Ellipsoid(name='Sphere', a=6371008.771415, b=6371008.771415, f_=0, f=0, f2=0, n=0, e=0, e2=0, 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) 46@var Ellipsoids.SphereAuthalic: Ellipsoid(name='SphereAuthalic', a=6371000, b=6371000, f_=0, f=0, f2=0, n=0, e=0, e2=0, e22=0, e32=0, A=6371000, L=10007543.39801029, R1=6371000, R2=6371000, R3=6371000, Rbiaxial=6371000, Rtriaxial=6371000) 47@var Ellipsoids.SpherePopular: Ellipsoid(name='SpherePopular', a=6378137, b=6378137, f_=0, f=0, f2=0, n=0, e=0, e2=0, e22=0, e32=0, A=6378137, L=10018754.17139462, R1=6378137, R2=6378137, R3=6378137, Rbiaxial=6378137, Rtriaxial=6378137) 48@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, 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) 49@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, e22=0.00673853, e32=0.00335795, A=6367478.63091189, L=10002012.0443814, R1=6371037.76231981, R2=6371036.17227197, R3=6371029.78316994, Rbiaxial=6367483.11833616, Rtriaxial=6372826.29723739) 50@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, 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) 51@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, 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) 52@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, 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) 53''' 54# make sure int/int division yields float quotient, see .basics 55from __future__ import division 56 57from pygeodesy.basics import copysign0, isfinite, isint, _xinstanceof 58from pygeodesy.errors import _AssertionError, _ValueError 59from pygeodesy.fmath import cbrt, cbrt2, fdot, fhorner, fpowers, Fsum, fsum_, \ 60 hypot, hypot_, hypot1, hypot2, sqrt3 61from pygeodesy.interns import EPS, EPS0, EPS02, EPS1, INF, NN, PI4, PI_2, R_M, _a_, \ 62 _Airy1830_, _AiryModified_, _Bessel1841_, _Clarke1866_, \ 63 _Clarke1880IGN_, _DOT_, _1_EPS, _EPStol as _TOL, _f_, \ 64 _finite_, _float as _F, _floatuple as _T, _GRS80_, _height_, \ 65 _Intl1924_, _Krassovski1940_, _Krassowsky1940_, _lat_, \ 66 _meridional_, _negative_, _not_, _null_, _prime_vertical_, \ 67 _radius_, _Sphere_, _SPACE_, _vs_, _WGS72_, _WGS84_, \ 68 _0_0, _0_5, _1_0, _2_0, _4_0, _90_0 69from pygeodesy.interns import _0_25, _3_0 # PYCHOK used! 70from pygeodesy.lazily import _ALL_LAZY 71from pygeodesy.named import _lazyNamedEnumItem as _lazy, _NamedEnum, \ 72 _NamedEnumItem, _NamedTuple, _Pass 73from pygeodesy.namedTuples import Distance2Tuple, Vector3Tuple, Vector4Tuple 74from pygeodesy.props import deprecated_Property_RO, Property_RO, property_doc_ 75from pygeodesy.streprs import Fmt, fstr, instr, strs, unstr 76from pygeodesy.units import Bearing_, Distance, Float, Float_, Height, Lam_, Lat, Meter, \ 77 Meter2, Meter3, Phi, Phi_, Radius, Radius_, Scalar 78from pygeodesy.utily import atand, atan2b, atan2d, degrees90, m2km, m2NM, m2SM, \ 79 m2radians, radians2m, sincos2d 80 81from math import asinh, atan, atanh, cos, degrees, exp, radians, sin, sinh, sqrt, tan 82 83R_M = Radius(R_M =R_M) # mean (spherical) earth radius (C{meter}) 84R_MA = Radius(R_MA=_F(6378137.0)) # equatorial earth radius (C{meter}), WGS84, EPSG:3785 85R_MB = Radius(R_MB=_F(6356752.3)) # polar earth radius (C{meter}), WGS84, EPSG:3785 86R_KM = Radius(R_KM=_F(m2km(R_M))) # mean (spherical) earth radius (C{KM}, kilo meter) 87R_NM = Radius(R_NM=_F(m2NM(R_M))) # mean (spherical) earth radius (C{NM}, nautical miles) 88R_SM = Radius(R_SM=_F(m2SM(R_M))) # mean (spherical) earth radius (C{SM}, statute miles) 89# See <https://www.EdWilliams.org/avform.htm>, 90# <https://www.DTIC.mil/dtic/tr/fulltext/u2/a216843.pdf> and 91# <https://GitHub.com/NASA/MultiDop/blob/master/src/share/man/man3/geog_lib.3> 92# based on International Standard Nautical Mile of 1,852 meter (1' latitude) 93R_FM = Radius(R_FM=_F(6371000.0)) # former FAI Sphere earth radius (C{meter}) 94R_GM = Radius(R_GM=_F(6371230.0)) # Avg. radius, distance to geoid surface (C{meter}) 95R_VM = Radius(R_VM=_F(6366707.0194937)) # Aviation/Navigation earth radius (C{meter}) 96# R_ = Radius(R_ =_F(6372797.560856)) # XXX some other earth radius??? 97 98__all__ = _ALL_LAZY.ellipsoids 99__version__ = '21.09.09' 100 101_f_0_0 = Float(f =_0_0) 102_f__0_0 = Float(f_=_0_0) 103# like U{WGS84_f()<https://GeographicLib.SourceForge.io/html/classGeographicLib_1_1Constants.html>} 104_f_WGS84 = Float(f = 1 / (1000000000 / 298257223563)) # 10000000000 / 2982572235629 105 106 107def _aux(lat, inverse, auxLat, clip=90): 108 '''Return named auxiliary latitude in C{degrees}. 109 ''' 110 return Lat(lat, clip=clip, name=_lat_ if inverse else auxLat) 111 112 113def _4Ecef(this, Ecef): 114 '''Return an ECEF converter. 115 ''' 116 from pygeodesy.ecef import EcefKarney, EcefVeness, EcefYou 117 118 if Ecef is None: 119 Ecef = EcefKarney 120 else: 121 _xinstanceof(EcefKarney, EcefVeness, EcefYou, Ecef=Ecef) 122 return Ecef(this, name=this.name) # datum or ellipsoid 123 124 125def _s2_c2(phi): 126 '''(INTERNAL) Return 2-tuple C{(sin(B{phi})**2, cos(B{phi})**2)}. 127 ''' 128 if phi: 129 s2 = sin(phi)**2 130 if s2 > 0: 131 c2 = _1_0 - s2 132 if c2 > 0: 133 if c2 < _1_0: 134 return s2, c2 135 else: 136 return _1_0, _0_0 137 return _0_0, _1_0 138 139 140class a_f2Tuple(_NamedTuple): 141 '''2-Tuple C{(a, f)} specifying an ellipsoid by I{equatorial} 142 radius C{a} in C{meter} and scalar I{flattening} C{f}. 143 144 @see: Class L{Ellipsoid2}. 145 ''' 146 _Names_ = (_a_, _f_) # name 'f' not 'f_' 147 _Units_ = (_Pass, _Pass) 148 149 def __new__(cls, a, f, **name): 150 '''New L{a_f2Tuple} ellipsoid specification. 151 152 @arg a: Equatorial radius (C{scalar} > 0). 153 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}). 154 155 @return: An L{a_f2Tuple}C{(a, f)} instance. 156 157 @raise UnitError: Invalid B{C{a}} or B{C{f}}. 158 159 @note: C{abs(B{f}) < EPS} is forced to C{B{f}=0}, I{spherical}. 160 Negative C{B{f}} produces a I{prolate} ellipsoid. 161 ''' 162 a = Radius_(a=a) 163 f = Float_( f=f, low=None, high=EPS1) 164 if abs(f) < EPS: # force spherical 165 f = _f_0_0 166 return _NamedTuple.__new__(cls, a, f, **name) 167 168 @Property_RO 169 def b(self): 170 '''Get the I{polar} radius (C{meter}), M{a * (1 - f)}. 171 ''' 172 return a_f2b(self.a, self.f) # PYCHOK .a and .f 173 174 @Property_RO 175 def f_(self): 176 '''Get the I{inverse} flattening (C{float}), M{1 / f} == M{a / (a - b)}. 177 ''' 178 return f2f_(self.f) # PYCHOK .f 179 180 181class Circle4Tuple(_NamedTuple): 182 '''4-Tuple C{(radius, height, lat, beta)} of the C{radius} and C{height}, 183 both conventionally in C{meter} of a parallel I{circle of latitude} at 184 (geodetic) latitude C{lat} and the I{parametric (or reduced) auxiliary 185 latitude} C{beta}, both in C{degrees90}. 186 187 The C{height} is the (signed) distance along the z-axis between the 188 parallel and the equator. At near-polar C{lat}s, the C{radius} is C{0}, 189 the C{height} is the ellipsoid's (signed) polar radius and C{beta} 190 equals C{lat}. 191 ''' 192 _Names_ = (_radius_, _height_, _lat_, 'beta') 193 _Units_ = ( Radius, Height, Lat, Lat) 194 195 196class Curvature2Tuple(_NamedTuple): 197 '''2-Tuple C{(meridional, prime_vertical)} of radii of curvature, 198 both in C{meter}, conventionally. 199 ''' 200 _Names_ = (_meridional_, _prime_vertical_) 201 _Units_ = ( Meter, Meter) 202 203 204class Ellipsoid(_NamedEnumItem): 205 '''Ellipsoid with I{equatorial} and I{polar} radii, I{flattening}, I{inverse 206 flattening} and other, often used, I{cached} attributes, supporting 207 I{oblate} and I{prolate} ellipsoidal and I{spherical} earth models. 208 ''' 209 _a = 0 # equatorial radius, semi-axis (C{meter}) 210 _b = 0 # polar radius, semi-axis (C{meter}): a * (f - 1) / f 211 _f = 0 # (1st) flattening: (a - b) / a 212 _f_ = 0 # inverse flattening: 1 / f = a / (a - b) 213 214 _KsOrder = 8 # Krüger series order (4, 6 or 8) 215 _Math = None # cached karney._wrapped.Math module 216 217 def __init__(self, a, b=None, f_=None, name=NN): 218 '''New L{Ellipsoid} from I{equatorial} and I{polar} radius or 219 I{equatorial} radius and I{inverse flattening}. 220 221 @arg a: Equatorial radius, semi-axis (C{meter}). 222 @arg b: Optional, polar radius, semi-axis (C{meter}). 223 @arg f_: Inverse flattening: M{a / (a - b)} (C{float} >>> 1.0). 224 @kwarg name: Optional, unique name (C{str}). 225 226 @raise NameError: Ellipsoid with that B{C{name}} already exists. 227 228 @raise ValueError: Invalid B{C{a}}, B{C{b}} or B{C{f_}}. 229 230 @note: M{abs(f_) > 1 / EPS} or M{abs(1 / f_) < EPS} is forced 231 to M{1 / f_ = 0}, spherical. 232 ''' 233 try: 234 a = Radius_(a=a) # low=EPS 235 if not isfinite(a): 236 raise ValueError(_not_(_finite_)) 237 238 if b: 239 b = Radius_(b=b) # low=EPS 240 f = a_b2f(a, b) 241 f_ = f2f_(f) if f_ is None else Float(f_=f_) 242 elif f_: 243 f_ = Float(f_=f_) 244 b = a_f_2b(a, f_) # a * (f_ - 1) / f_ 245 f = a_b2f(a, b) 246 else: # only a, spherical 247 f = f_ = 0 248 b = a # superfluous 249 250 if not isfinite(b): 251 raise ValueError(_not_(_finite_)) 252 253 if abs(f) < EPS or a == b or not f_: # spherical 254 b = a 255 f = _f_0_0 256 f_ = _f__0_0 257 elif f > EPS1: # sanity check, see .ecef.Ecef.__init__ 258 raise _ValueError(f=f) 259 260 except (TypeError, ValueError) as x: 261 t = instr(self, a=a, b=b, f_=f_, name=name) 262 raise _ValueError(t, txt=str(x)) 263 264 self._a = a 265 self._b = b 266 self._f = f 267 self._f_ = f_ 268 269 self._register(Ellipsoids, name) 270 271 if f and f_: # see .test/testEllipsoidal.py 272 self._assert(_1_0 / f, f_=f_, eps=_TOL) 273 self._assert(_1_0 / f_, f =f, eps=_TOL) 274 self._assert(self.b2_a2, e12=self.e12, eps=EPS) 275 276 def __eq__(self, other): 277 '''Compare this and an other ellipsoid. 278 279 @arg other: The other ellipsoid (L{Ellipsoid} or L{Ellipsoid2}). 280 281 @return: C{True} if equal, C{False} otherwise. 282 ''' 283 return self is other or (isinstance(other, Ellipsoid) and 284 self.a == other.a and 285 (self.b == other.b or self.f == other.f)) 286 287 @Property_RO 288 def a(self): 289 '''Get the I{equatorial} radius, semi-axis (C{meter}). 290 ''' 291 return self._a 292 293 equatoradius = a # = Requatorial 294 295 @Property_RO 296 def a2(self): 297 '''Get the I{equatorial} radius I{squared} (C{meter**2}), M{a**2}. 298 ''' 299 return Meter2(a2=self.a**2) 300 301 @Property_RO 302 def a2_(self): 303 '''Get the inverse of the I{equatorial} radius I{squared} (C{meter**2}), M{1 / a**2}. 304 ''' 305 return Float(a2_=_1_0 / self.a2) 306 307 @Property_RO 308 def a_b(self): 309 '''Get the ratio I{equatorial} over I{polar} radius (C{float}), M{a / b} == M{1 / (1 - f)}. 310 ''' 311 return Float(a_b=self.a / self.b if self.f else _1_0) 312 313 @Property_RO 314 def a2_b(self): 315 '''Get the I{polar} meridional radius of curvature (C{meter}), M{a**2 / b}. 316 317 @see: U{Radii of Curvature 318 <https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>} 319 and U{Moritz, H. (1980), Geodetic Reference System 1980 320 <https://WikiPedia.org/wiki/Earth_radius#cite_note-Moritz-2>}. 321 322 @note: Symbol C{c} is used by IUGG and IERS for the U{polar radius of curvature 323 <https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>}, see L{c2} and 324 L{R2} or L{Rauthalic}. 325 ''' 326 return Radius(a2_b=self.a2 / self.b if self.f else self.a) # = rocPolar 327 328 @Property_RO 329 def a2_b2(self): 330 '''Get the ratio I{equatorial} over I{polar} radius I{squared} (C{float}), M{(a / b)**2} 331 == M{1 / (1 - e**2)} == M{1 / (1 - e2)} == M{1 / e12}. 332 ''' 333 return Float(a2_b2=self.a_b**2 if self.f else _1_0) 334 335 @Property_RO 336 def a_f(self): 337 '''Get the I{equatorial} radius and I{flattening} (L{a_f2Tuple}). 338 ''' 339 return a_f2Tuple(self.a, self.f, name=self.name) 340 341 @Property_RO 342 def A(self): 343 '''Get the UTM I{meridional (or rectifying)} radius (C{meter}). 344 ''' 345 A, n = self.a, self.n 346 if n: 347 n1 = _1_0 + n 348 if n1: # use 6 n**2 terms, half-way between the _KsOrder's 4, 6, 8 349 # <https://GeographicLib.SourceForge.io/html/tmseries30.html> 350 # <https://GeographicLib.SourceForge.io/html/transversemercator.html> and 351 # <https://www.MyGeodesy.id.AU/documents/Karney-Krueger%20equations.pdf> (3) 352 A = Radius(A=A / n1 * (fhorner(n**2, 1048576, 262144, 16384, 4096, 1600, 784, 441) / 1048576)) 353 return A 354 355 @Property_RO 356 def _albersCyl(self): 357 '''(INTERNAL) Helper for C{auxAuthalic}. 358 ''' 359 from pygeodesy.albers import AlbersEqualAreaCylindrical 360 return AlbersEqualAreaCylindrical(datum=self, name=self.name) 361 362 @Property_RO 363 def AlphaKs(self): 364 '''Get the I{Krüger} U{Alpha series coefficients<https://GeographicLib.SourceForge.io/html/tmseries30.html>} (C{KsOrder}C{-tuple}). 365 ''' 366 return self._Kseries( # XXX int/int quotients may require from __future__ import division 367 # n n**2 n**3 n**4 n**5 n**6 n**7 n**8 368 _T(1/2, -2/3, 5/16, 41/180, -127/288, 7891/37800, 72161/387072, -18975107/50803200), 369 _T(13/48, -3/5, 557/1440, 281/630, -1983433/1935360, 13769/28800, 148003883/174182400), # PYCHOK unaligned 370 _T(61/240, -103/140, 15061/26880, 167603/181440, -67102379/29030400, 79682431/79833600), # PYCHOK unaligned 371 _T(49561/161280, -179/168, 6601661/7257600, 97445/49896, -40176129013/7664025600), # PYCHOK unaligned 372 _T(34729/80640, -3418889/1995840, 14644087/9123840, 2605413599/622702080), # PYCHOK unaligned 373 _T(212378941/319334400, -30705481/10378368, 175214326799/58118860800), # PYCHOK unaligned 374 _T(1522256789/1383782400, -16759934899/3113510400), # PYCHOK unaligned 375 _T(1424729850961/743921418240)) # PYCHOK unaligned 376 377 @Property_RO 378 def area(self): 379 '''Get the ellipsoid's surface area (C{meter**2}), M{4 * PI * c2}. 380 381 @see: L{areax}, L{c2} and L{R2}. 382 ''' 383 return Meter2(area=self.c2 * PI4) 384 385 @Property_RO 386 def areax(self): 387 '''Get the ellipsoid's surface area (C{meter**2}), M{4 * PI * c2x}, 388 more accurate for very I{oblate} ellipsoids. 389 390 @see: L{area}, L{c2x}, L{R2x} and L{GeodesicExact}. 391 ''' 392 return Meter2(areax=self.c2x * PI4) 393 394 def _assert(self, val, eps=_TOL, f0=_0_0, **name_value): 395 '''(INTERNAL) Assert a C{name=value} vs C{val}. 396 ''' 397 for n, v in name_value.items(): 398 if abs(v - val) > eps: 399 t = (v, _vs_, val) 400 t = _SPACE_.join(strs(t, prec=12, fmt=Fmt.g)) 401 t = Fmt.EQUAL(self._DOT_(n), t) 402 raise _AssertionError(t, txt=Fmt.exceeds_eps(eps)) 403 return Float(v if self.f else f0, name=n) 404 raise _AssertionError(unstr(self._DOT_(self._assert.__name__), val, 405 eps=eps, f0=f0, **name_value)) 406 407 def auxAuthalic(self, lat, inverse=False): 408 '''Compute the I{authalic} auxiliary latitude or inverse thereof. 409 410 @arg lat: The geodetic (or I{authalic}) latitude (C{degrees90}). 411 @kwarg inverse: If C{True}, B{C{lat}} is the I{authalic} and 412 return the geodetic latitude (C{bool}). 413 414 @return: The I{authalic} (or geodetic) latitude in C{degrees90}. 415 416 @see: U{Inverse-/AuthalicLatitude<https://GeographicLib.SourceForge.io/ 417 html/classGeographicLib_1_1Ellipsoid.html>}, U{Authalic latitude 418 <https://WikiPedia.org/wiki/Latitude#Authalic_latitude>}, and 419 U{Snyder<https://Pubs.USGS.gov/pp/1395/report.pdf>}, p 16. 420 421 ''' 422 if self.f: 423 f = self._albersCyl._tanf if inverse else self._albersCyl._txif # PYCHOK attr 424 lat = atand(f(tan(Phi_(lat)))) # PYCHOK attr 425 return _aux(lat, inverse, Ellipsoid.auxAuthalic.__name__) 426 427 def auxConformal(self, lat, inverse=False): 428 '''Compute the I{conformal} auxiliary latitude or inverse thereof. 429 430 @arg lat: The geodetic (or I{conformal}) latitude (C{degrees90}). 431 @kwarg inverse: If C{True}, B{C{lat}} is the I{conformal} and 432 return the geodetic latitude (C{bool}). 433 434 @return: The I{conformal} (or geodetic) latitude in C{degrees90}. 435 436 @see: U{Inverse-/ConformalLatitude<https://GeographicLib.SourceForge.io/ 437 html/classGeographicLib_1_1Ellipsoid.html>}, U{Conformal latitude 438 <https://WikiPedia.org/wiki/Latitude#Conformal_latitude>}, and 439 U{Snyder<https://Pubs.USGS.gov/pp/1395/report.pdf>}, pp 15-16. 440 ''' 441 if self.f: 442 f = self.es_tauf if inverse else self.es_taupf # PYCHOK attr 443 lat = atand(f(tan(Phi_(lat)))) # PYCHOK attr 444 return _aux(lat, inverse, Ellipsoid.auxConformal.__name__) 445 446 def auxGeocentric(self, lat, inverse=False): 447 '''Compute the I{geocentric} auxiliary latitude or inverse thereof. 448 449 @arg lat: The geodetic (or I{geocentric}) latitude (C{degrees90}). 450 @kwarg inverse: If C{True}, B{C{lat}} is the geocentric and 451 return the I{geocentric} latitude (C{bool}). 452 453 @return: The I{geocentric} (or geodetic) latitude in C{degrees90}. 454 455 @see: U{Inverse-/GeocentricLatitude<https://GeographicLib.SourceForge.io/ 456 html/classGeographicLib_1_1Ellipsoid.html>}, U{Geocentric latitude 457 <https://WikiPedia.org/wiki/Latitude#Geocentric_latitude>}, and 458 U{Snyder<<https://Pubs.USGS.gov/pp/1395/report.pdf>}, pp 17-18. 459 ''' 460 if self.f: 461 f = self.a2_b2 if inverse else self.b2_a2 462 lat = atand(f * tan(Phi_(lat))) 463 return _aux(lat, inverse, Ellipsoid.auxGeocentric.__name__) 464 465 def auxIsometric(self, lat, inverse=False): 466 '''Compute the I{isometric} auxiliary latitude or inverse thereof. 467 468 @arg lat: The geodetic (or I{isometric}) latitude (C{degrees}). 469 @kwarg inverse: If C{True}, B{C{lat}} is the I{isometric} and 470 return the geodetic latitude (C{bool}). 471 472 @return: The I{isometric} (or geodetic) latitude in C{degrees}. 473 474 @note: The I{isometric} latitude for geodetic C{+/-90} is far 475 outside the C{[-90..+90]} range but the inverse 476 thereof is the original geodetic latitude. 477 478 @see: U{Inverse-/IsometricLatitude<https://GeographicLib.SourceForge.io/ 479 html/classGeographicLib_1_1Ellipsoid.html>}, U{Isometric latitude 480 <https://WikiPedia.org/wiki/Latitude#Isometric_latitude>}, and 481 U{Snyder<https://Pubs.USGS.gov/pp/1395/report.pdf>}, pp 15-16. 482 ''' 483 if self.f: 484 r = Phi_(lat, clip=0) 485 lat = degrees( atan(self.es_tauf(sinh(r))) if inverse else 486 asinh(self.es_taupf(tan(r)))) 487 # clip=0, since auxIsometric(+/-90) is far outside [-90..+90] 488 return _aux(lat, inverse, Ellipsoid.auxIsometric.__name__, clip=0) 489 490 def auxParametric(self, lat, inverse=False): 491 '''Compute the I{parametric} auxiliary latitude or inverse thereof. 492 493 @arg lat: The geodetic (or I{parametric}) latitude (C{degrees90}). 494 @kwarg inverse: If C{True}, B{C{lat}} is the I{parametric} and 495 return the geodetic latitude (C{bool}). 496 497 @return: The I{parametric} (or geodetic) latitude in C{degrees90}. 498 499 @see: U{Inverse-/ParametricLatitude<https://GeographicLib.SourceForge.io/ 500 html/classGeographicLib_1_1Ellipsoid.html>}, U{Parametric latitude 501 <https://WikiPedia.org/wiki/Latitude#Parametric_(or_reduced)_latitude>}, 502 and U{Snyder<https://Pubs.USGS.gov/pp/1395/report.pdf>}, p 18. 503 ''' 504 if self.f: 505 lat = self._beta(Lat(lat), inverse=inverse) 506 return _aux(lat, inverse, Ellipsoid.auxParametric.__name__) 507 508 auxReduced = auxParametric # synonyms 509 510 def auxRectifying(self, lat, inverse=False): 511 '''Compute the I{rectifying} auxiliary latitude or inverse thereof. 512 513 @arg lat: The geodetic (or I{rectifying}) latitude (C{degrees90}). 514 @kwarg inverse: If C{True}, B{C{lat}} is the I{rectifying} and 515 return the geodetic latitude (C{bool}). 516 517 @return: The I{rectifying} (or geodetic) latitude in C{degrees90}. 518 519 @see: U{Inverse-/RectifyingLatitude<https://GeographicLib.SourceForge.io/ 520 html/classGeographicLib_1_1Ellipsoid.html>}, U{Rectifying latitude 521 <https://WikiPedia.org/wiki/Latitude#Rectifying_latitude>}, and 522 U{Snyder<https://Pubs.USGS.gov/pp/1395/report.pdf>}, pp 16-17. 523 ''' 524 if self.f: 525 lat = Lat(lat) 526 if 0 < abs(lat) < _90_0: 527 if inverse: 528 e = self._elliptic_e22 529 lat = degrees90(e.fEinv(e.cE * lat / _90_0)) 530 lat = self.auxParametric(lat, inverse=True) 531 else: 532 lat = _90_0 * self.Llat(lat) / self.L 533 return _aux(lat, inverse, Ellipsoid.auxRectifying.__name__) 534 535 @Property_RO 536 def b(self): 537 '''Get the I{polar} radius, semi-axis (C{meter}). 538 ''' 539 return self._b 540 541 polaradius = b # = Rpolar 542 543 @Property_RO 544 def b_a(self): 545 '''Get the ratio I{polar} over I{equatorial} radius (C{float}), M{b / a} == M{1 - f}. 546 ''' 547 return self._assert(self.b / self.a, b_a=_1_0 - self.f, f0=_1_0) 548 549 @Property_RO 550 def b2(self): 551 '''Get the I{polar} radius I{squared} (C{float}), M{b**2}. 552 ''' 553 return Meter2(b2=self.b**2) 554 555 @Property_RO 556 def b2_a(self): 557 '''Get the I{equatorial} meridional radius of curvature (C{meter}), M{b**2 / a}, see C{rocMeridional}C{(0)}. 558 559 @see: U{Radii of Curvature 560 <https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>}. 561 ''' 562 return Radius(b2_a=self.b2 / self.a if self.f else self.b) 563 564 @Property_RO 565 def b2_a2(self): 566 '''Get the ratio I{polar} over I{equatorial} radius I{squared} (C{float}), M{(b / a)**2} 567 == M{(1 - f)**2} == M{1 - e**2} == C{e12}. 568 ''' 569 return Float(b2_a2=self.b_a**2 if self.f else _1_0) 570 571 def _beta(self, lat, inverse=False): 572 '''(INTERNAL) Get the I{parametric (or reduced) auxiliary latitude} or inverse thereof. 573 ''' 574 s, c = sincos2d(lat) # like Karney's tand(lat) 575 s *= self.a_b if inverse else self.b_a 576 return atan2d(s, c) # == atand(s / c) if c else copysign0(_90_0, lat) 577 578 @Property_RO 579 def BetaKs(self): 580 '''Get the I{Krüger} U{Beta series coefficients<https://GeographicLib.SourceForge.io/html/tmseries30.html>} (C{KsOrder}C{-tuple}). 581 ''' 582 return self._Kseries( # XXX int/int quotients may require from __future__ import division 583 # n n**2 n**3 n**4 n**5 n**6 n**7 n**8 584 _T(1/2, -2/3, 37/96, -1/360, -81/512, 96199/604800, -5406467/38707200, 7944359/67737600), 585 _T(1/48, 1/15, -437/1440, 46/105, -1118711/3870720, 51841/1209600, 24749483/348364800), # PYCHOK unaligned 586 _T(17/480, -37/840, -209/4480, 5569/90720, 9261899/58060800, -6457463/17740800), # PYCHOK unaligned 587 _T(4397/161280, -11/504, -830251/7257600, 466511/2494800, 324154477/7664025600), # PYCHOK unaligned 588 _T(4583/161280, -108847/3991680, -8005831/63866880, 22894433/124540416), # PYCHOK unaligned 589 _T(20648693/638668800, -16363163/518918400, -2204645983/12915302400), # PYCHOK unaligne 590 _T(219941297/5535129600, -497323811/12454041600), # PYCHOK unaligned 591 _T(191773887257/3719607091200)) # PYCHOK unaligned 592 593 @deprecated_Property_RO 594 def c(self): 595 '''DEPRECATED, use property C{R2} or C{Rauthalic}.''' 596 return self.R2 597 598 @Property_RO 599 def c2(self): 600 '''Get the I{authalic} earth radius I{squared} (C{meter**2}). 601 602 603 @see: L{c2x}, L{area}, L{R2}, L{Rauthalic}, I{Karney's} U{equation 60 604 <https://Link.Springer.com/article/10.1007%2Fs00190-012-0578-z>} and C++ U{Ellipsoid.Area() 605 <https://GeographicLib.SourceForge.io/html/classGeographicLib_1_1Ellipsoid.html>}, 606 U{Authalic radius<https://WikiPedia.org/wiki/Earth_radius#Authalic_radius>}, U{Surface area 607 <https://WikiPedia.org/wiki/Ellipsoid>} and U{surface area 608 <https://www.Numericana.com/answer/geometry.htm#oblate>}. 609 ''' 610 return self._c2f(False) 611 612 @Property_RO 613 def c2x(self): 614 '''Get the I{authalic} earth radius I{squared} (C{meter**2}), more accurate for very I{oblate} 615 ellipsoids. 616 617 @see: L{c2}, L{areax}, L{R2x}, L{Rauthalicx}, L{GeodesicExact} and I{Karney}'s comments at C++ 618 attribute U{GeodesicExact._c2<https://GeographicLib.SourceForge.io/html/GeodesicExact_8cpp_source.html>}. 619 ''' 620 return self._c2f(True) 621 622 def _c2f(self, c2x): 623 '''(INTERNAL) Helper for C{.c2} and C{.c2x}. 624 ''' 625 f = self.f 626 if f: 627 c2, e = self.b2, self.e 628 if e > EPS0: 629 if f > 0: # .isOblate 630 c2 *= (asinh(sqrt(self.e22abs)) if c2x else atanh(e)) / e 631 elif f < 0: # .isProlate 632 c2 *= atan(e) / e # XXX asin? 633 c2 = Meter2(c2=(self.a2 + c2) * _0_5) 634 else: 635 c2 = self.a2 636 return c2 637 638 def circle4(self, lat): 639 '''Get the equatorial or a parallel I{circle of latitude}. 640 641 @arg lat: Geodetic latitude (C{degrees90}, C{str}). 642 643 @return: A L{Circle4Tuple}C{(radius, height, lat, beta)} 644 instance. 645 646 @raise RangeError: Latitude B{C{lat}} outside valid range 647 and L{rangerrors} set to C{True}. 648 649 @raise TypeError: Invalid B{C{lat}}. 650 651 @raise ValueError: Invalid B{C{lat}}. 652 653 @see: Definition of U{I{p} and I{z} under B{Parametric (or 654 reduced) latitude}<https://WikiPedia.org/wiki/Latitude>} 655 and I{Karney's} C++ U{CircleRadius and CircleHeight 656 <https://GeographicLib.SourceForge.io/html/classGeographicLib_1_1Ellipsoid.html>}. 657 ''' 658 lat = Lat(lat) 659 if lat: 660 b = lat 661 if abs(lat) < _90_0: 662 if self.f: 663 b = self._beta(lat) 664 z, r = sincos2d(b) 665 r *= self.a 666 z *= self.b 667 else: # near-polar 668 r, z = _0_0, copysign0(self.b, lat) 669 else: # equator 670 r = self.a 671 z = lat = b = _0_0 672 return Circle4Tuple(r, z, lat, b) 673 674 def degrees2m(self, deg, lat=0): 675 '''Convert an angle to the distance along the equator or 676 along a parallel of (geodetic) latitude. 677 678 @arg deg: The angle (C{degrees}). 679 @kwarg lat: Parallel latitude (C{degrees90}, C{str}). 680 681 @return: Distance (C{meter}, same units as the equatorial 682 and polar radii) or C{0} for near-polar B{C{lat}}. 683 684 @raise RangeError: Latitude B{C{lat}} outside valid range 685 and L{rangerrors} set to C{True}. 686 687 @raise ValueError: Invalid B{C{deg}} or B{C{lat}}. 688 ''' 689 return self.radians2m(radians(deg), lat=lat) 690 691 def distance2(self, lat0, lon0, lat1, lon1): 692 '''I{Approximate} the distance and (initial) bearing between 693 two points based on the U{local, flat earth approximation 694 <https://www.EdWilliams.org/avform.htm#flat>} aka U{Hubeny 695 <https://www.OVG.AT/de/vgi/files/pdf/3781/>} formula. 696 697 I{Suitable only for distances of several hundred Km or Miles 698 and only between points not near-polar}. 699 700 @arg lat0: From latitude (C{degrees}). 701 @arg lon0: From longitude (C{degrees}). 702 @arg lat1: To latitude (C{degrees}). 703 @arg lon1: To longitude (C{degrees}). 704 705 @return: A L{Distance2Tuple}C{(distance, initial)} with C{distance} 706 in same units as this ellipsoid's axes. 707 708 @note: The meridional and prime_vertical radii of curvature are 709 taken and scaled I{at the initial latitude}, see C{roc2}. 710 711 @see: Function L{flatLocal}/L{hubeny}. 712 ''' 713 phi0 = Phi_(lat0=lat0) 714 m, n = self.roc2_(phi0, scaled=True) 715 m *= Phi_(lat1=lat1) - phi0 716 n *= Lam_(lon1=lon1) - Lam_(lon0=lon0) 717 return Distance2Tuple(hypot(m, n), atan2b(n, m)) 718 719 @Property_RO 720 def e(self): 721 '''Get the I{(1st) eccentricity} (C{float}), M{sqrt(1 - (b / a)**2))}, see C{a_b2e}. 722 ''' 723 return Float(e=sqrt(self.e2abs) if self.e2 else _0_0) 724 725 eccentricity = e # eccentricity 726 727 @Property_RO 728 def e2(self): 729 '''Get the I{(1st) eccentricity squared} (C{float}), M{f * (2 - f) 730 == 1 - (b / a)**2}, see C{a_b2e2}. 731 ''' 732 return self._assert(a_b2e2(self.a, self.b), e2=f2e2(self.f)) 733 734 eccentricity2 = e2 # eccentricity squared 735 736 @Property_RO 737 def e2abs(self): 738 '''Get C{abs} value of the I{(1st) eccentricity squared} (C{float}). 739 ''' 740 return abs(self.e2) 741 742 @Property_RO 743 def e12(self): 744 '''Get 1 less I{(1st) eccentricity squared} (C{float}), M{1 - e**2 745 == (1 - f)**2} == M{b**2 / a**2}, see C{b2_a2}. 746 ''' 747 return self._assert((_1_0 - self.f)**2, e12=_1_0 - self.e2, f0=_1_0) # 1 - e2 748 749 @Property_RO 750 def e22(self): 751 '''Get the I{2nd eccentricity squared} (C{float}), M{e2 / (1 - e2) 752 == e2 / (1 - f)**2 == (a / b)**2 - 1}, see C{a_b2e22}. 753 ''' 754 return self._assert(a_b2e22(self.a, self.b), e22=f2e22(self.f)) 755 756 eccentricity2nd2 = e22 # second eccentricity squared 757 758 @Property_RO 759 def e22abs(self): 760 '''Get C{abs} value of the I{(2nd) eccentricity squared} (C{float}). 761 ''' 762 return abs(self.e22) 763 764 @Property_RO 765 def e32(self): 766 '''Get the I{3rd eccentricity squared} (C{float}), M{e2 / (2 - e2) 767 == (a**2 - b**2) / (a**2 + b**2)}, see C{a_b2e32}. 768 ''' 769 return self._assert(a_b2e32(self.a, self.b), e32=f2e32(self.f)) 770 771 eccentricity3rd2 = e32 # third eccentricity squared 772 773 @Property_RO 774 def e32abs(self): 775 '''Get C{abs} value of the I{(3rd) eccentricity squared} (C{float}). 776 ''' 777 return abs(self.e32) 778 779 @Property_RO 780 def e4(self): 781 '''Get the I{(1st) eccentricity} to 4th power (C{float}), M{e**4 == e2**2}. 782 ''' 783 return Float(e4=self.e2**2 if self.e2 else _0_0) 784 785 def ecef(self, Ecef=None): 786 '''Return U{ECEF<https://WikiPedia.org/wiki/ECEF>} converter. 787 788 @kwarg Ecef: ECEF class to use (L{EcefKarney}, L{EcefVeness} 789 or L{EcefYou}). 790 791 @return: An ECEF converter for this C{ellipsoid} (L{EcefKarney}, 792 L{EcefVeness} or L{EcefYou}). 793 794 @raise TypeError: Invalid B{C{Ecef}}. 795 ''' 796 return _4Ecef(self, Ecef) 797 798 @Property_RO 799 def _elliptic_e22(self): 800 '''(INTERNAL) Elliptic helper for C{auxRectifying}, C{L}, C{Llat}. 801 ''' 802 from pygeodesy.elliptic import Elliptic 803 return Elliptic(-abs(self.e22)) 804 805 def e2s(self, s): 806 '''Compute norm M{sqrt(1 - e2 * s**2)}. 807 808 @arg s: Sine value (C{scalar}). 809 810 @return: Norm (C{float}). 811 812 @raise ValueError: Invalid B{C{s}}. 813 ''' 814 return sqrt(self.e2s2(s)) if self.e2 else _1_0 815 816 def e2s2(self, s): 817 '''Compute M{1 - e2 * s**2}. 818 819 @arg s: S value (C{scalar}). 820 821 @return: Result (C{float}). 822 823 @raise ValueError: Invalid B{C{s}}. 824 ''' 825 r = _1_0 826 if self.e2: 827 try: 828 r -= self.e2 * Scalar(s=s)**2 829 if r < 0: 830 raise ValueError(_negative_) 831 except (TypeError, ValueError) as x: 832 t = self._DOT_(Ellipsoid.e2s2.__name__) 833 raise _ValueError(t, s, txt=str(x)) 834 return r 835 836 @Property_RO 837 def es(self): 838 '''Get the I{(1st) eccentricity signed} (C{float}). 839 ''' 840 # note, self.e is always non-negative 841 return Float(es=copysign0(self.e, self.f)) # see .ups 842 843 def es_atanh(self, x): 844 '''Compute M{es * atanh(es * x)} where I{es} is the I{signed} 845 (1st) eccentricity. 846 847 @raise ValueError: Invalid B{C{x}}. 848 849 @see: Function U{Math::eatanhe<https://GeographicLib.SourceForge.io/ 850 html/classGeographicLib_1_1Math.html>}. 851 ''' 852 # note, self.e is always non-negative 853 if self.f > 0: # .isOblate 854 r = atanh(self.e * Scalar(x=x)) * self.e 855 elif self.f < 0: # .isProlate 856 r = -atan(self.e * Scalar(x=x)) * self.e 857 else: # .isSpherical 858 r = _0_0 859 return r 860 861 @Property_RO 862 def es_c(self): 863 '''Get M{(1 - f) * exp(es_atanh(1))} (C{float}), M{b_a * exp(es_atanh(1))}. 864 ''' 865 return Float(es_c=(exp(self.es_atanh(_1_0)) * self.b_a) if self.f else _1_0) 866 867 def es_tauf(self, taup): 868 '''Compute I{Karney}'s U{equations (19), (20) and (21) 869 <https://ArXiv.org/abs/1002.1417>}. 870 871 @see: U{Math::tauf<https://GeographicLib.SourceForge.io/ 872 html/classGeographicLib_1_1Math.html>}. 873 ''' 874 tp = Scalar(taup=taup) 875 if not self.f: # .isSpherical 876 return tp 877 tol = max(abs(tp), _1_0) * _TOL 878 # To lowest order in e^2, taup = (1 - e^2) * tau, so use starting 879 # guess tau = taup / (1 - e^2) = taup * a^2 / b^2. This starting 880 # guess is the geocentric latitude which, to first order in the 881 # flattening, is equal to the conformal auxiliary latitude. 882 e = self.a2_b2 # == _1_0 / self.e12 883 t = tp * e 884 T = Fsum(t) 885 for _ in range(9): 886 a, h = self._es_taupf2(t) 887 d = (tp - a) * (e + t**2) / (h * hypot1(a)) 888 t, d = T.fsum2_(d) 889 if abs(d) < tol: 890 break 891 return t 892 893 def es_taupf(self, tau): 894 '''Compute I{Karney}'s U{equations (7), (8) and (9) 895 <https://ArXiv.org/abs/1002.1417>}. 896 897 @see: U{Math::taupf<https://GeographicLib.SourceForge.io/ 898 html/classGeographicLib_1_1Math.html>}. 899 ''' 900 t = Scalar(tau=tau) 901 if self.f: 902 t, _ = self._es_taupf2(t) 903 return t 904 905 def _es_taupf2(self, t): 906 '''(INTERNAL) Return C{(es_taupf(t), hypot1(t))}. 907 ''' 908 h = hypot1(t) 909 s = sinh(self.es_atanh(t / h)) 910 a = hypot1(s) * t - h * s 911 return a, h 912 913 @Property_RO 914 def f(self): 915 '''Get the I{flattening} (C{float}), M{(a - b) / a}, C{0} for spherical, negative for prolate. 916 ''' 917 return self._f 918 919 flattening = f 920 921 @Property_RO 922 def f_(self): 923 '''Get the I{inverse flattening} (C{float}), M{1 / f} == M{a / (a - b)}, C{0} for spherical, see C{a_b2f_}. 924 ''' 925 return self._f_ 926 927 @Property_RO 928 def f2(self): 929 '''Get the I{2nd flattening} (C{float}), M{(a - b) / b == f / (1 - f)}, C{0} for spherical, see C{a_b2f2}. 930 ''' 931 return self._assert(self.a_b - _1_0, f2=f2f2(self.f)) 932 933 flattening2nd = f2 934 935 @Property_RO 936 def geodesic(self): 937 '''Get this ellipsoid's I{wrapped Karney} U{Geodesic 938 <https://GeographicLib.SourceForge.io/html/python/code.html>}, 939 provided the U{geographiclib<https://PyPI.org/project/geographiclib>} 940 package is installed. 941 ''' 942 # if not self.isEllipsoidal: 943 # raise _IsnotError(_ellipsoidal_, ellipsoid=self) 944 from pygeodesy.karney import _wrapped 945 g = _wrapped.Geodesic(self) 946 Ellipsoid._Math = _wrapped.Math # hold 947 return g 948 949 @Property_RO 950 def _geodesic_Math2(self): 951 '''(INTERNAL) Get this ellipsoid's I{wrapped Karney} C{Geodesic} 952 and I{Karney}'s C{Math} class, see L{geodesic}. 953 ''' 954 g = self.geodesic 955 return g, Ellipsoid._Math 956 957 @Property_RO 958 def geodesicx(self): 959 '''Get this ellipsoid's L{GeodesicExact}. 960 ''' 961 # if not self.isEllipsoidal: 962 # raise _IsnotError(_ellipsoidal_, ellipsoid=self) 963 from pygeodesy.geodesicx import GeodesicExact 964 return GeodesicExact(self, name=self.name) 965 966 @Property_RO 967 def geodsolve(self): 968 '''Get this ellipsoid's L{GeodesicSolve}, the I{wrapper} around utility 969 U{GeodSolve<https://GeographicLib.SourceForge.io/html/GeodSolve.1.html>}, 970 provided the path to the C{GeodSolve} executable is specified with env 971 variable C{PYGEODESY_GEODSOLVE}. 972 ''' 973 # if not self.isEllipsoidal: 974 # raise _IsnotError(_ellipsoidal_, ellipsoid=self) 975 from pygeodesy.geodsolve import GeodesicSolve 976 return GeodesicSolve(self, name=self.name) 977 978 def height4(self, xyz, normal=True): 979 '''Compute the height of a cartesian above or below and the projection 980 on this ellipsoid's surface. 981 982 @arg xyz: The cartesian (C{Cartesian}, L{Ecef9Tuple}, L{Vector3d}, 983 L{Vector3Tuple} or L{Vector4Tuple}). 984 @kwarg normal: If C{True} the projection is the nearest point on the 985 ellipsoid's surface, otherwise the intersection of the 986 radial line to the center and the ellipsoid's surface. 987 988 @return: L{Vector4Tuple}C{(x, y, z, h)} with the intersection C{x}, C{y} 989 and C{z} coordinates and height C{h} in C{meter}, conventionally. 990 991 @raise ValueError: Null B{C{xyz}}. 992 993 @raise TypeError: Non-cartesian B{C{xyz}}. 994 995 @see: U{Distance to<https://StackOverflow.com/questions/22959698/distance-from-given-point-to-given-ellipse>} 996 and U{intersection with<https://MathWorld.wolfram.com/Ellipse-LineIntersection.html>} an ellipse. 997 ''' 998 from pygeodesy.vector3d import _otherV3d 999 v = _otherV3d(xyz=xyz) 1000 r = v.length 1001 if r < EPS0: # EPS 1002 raise _ValueError(xyz=xyz, txt=_null_) 1003 1004 a, b = self.a, self.b 1005 if self.isSpherical: 1006 v = v.times(a / r) 1007 h = r - a 1008 1009 elif normal: # perpendicular to ellipsoid 1010 x, y = hypot(v.x, v.y), abs(v.z) 1011 if x < EPS0: # polar 1012 z = copysign0(b, v.z) 1013 v = Vector3Tuple(v.x, v.y, z) 1014 h = y - b 1015 elif y < EPS0: # equatorial 1016 t = a / r 1017 v = v.times_(t, t, 0) # force z=0 1018 h = x - a 1019 else: # normal in 1st quadrant 1020 x, y = _normal2(x, y, self) 1021 t, v = v, v.times_(x, x, y) 1022 h = t.minus(v).length 1023 1024 else: # radial to ellipsoid center 1025 t = hypot_(a * v.z, b * v.x, b * v.y) 1026 if t < EPS0: # EPS 1027 raise _ValueError(xyz=xyz, txt=_null_) 1028 t = a * b / t 1029 v = v.times(t) 1030 h = r * (_1_0 - t) 1031 1032 return Vector4Tuple(v.x, v.y, v.z, h, name=self.height4.__name__) 1033 1034 def _hubeny2_(self, phi2, phi1, lam21): 1035 '''(INTERNAL) like function L{flatLocal_}/L{hubeny_} but 1036 returning the I{angular} distance in C{radians squared}. 1037 ''' 1038 m, n = self.roc2_((phi2 + phi1) * _0_5, scaled=True) 1039 return hypot2(m * (phi2 - phi1), n * lam21) * self.a2_ 1040 1041 @Property_RO 1042 def isEllipsoidal(self): 1043 '''Is this model I{ellipsoidal} (C{bool})? 1044 ''' 1045 return self.f != 0 1046 1047 @Property_RO 1048 def isOblate(self): 1049 '''Is this ellipsoid I{oblate} (C{bool})? I{Prolate} or 1050 spherical otherwise. 1051 ''' 1052 return self.f > 0 1053 1054 @Property_RO 1055 def isProlate(self): 1056 '''Is this ellipsoid I{prolate} (C{bool})? I{Oblate} or 1057 spherical otherwise. 1058 ''' 1059 return self.f < 0 1060 1061 @Property_RO 1062 def isSpherical(self): 1063 '''Is this model I{spherical} (C{bool})? 1064 ''' 1065 return self.f == 0 1066 1067 def _Kseries(self, *AB8Ks): 1068 '''(INTERNAL) Compute the 4-, 6- or 8-th order I{Krüger} Alpha 1069 or Beta series coefficients per I{Karney}'s U{equations 35 1070 and 36<https://Arxiv.org/pdf/1002.1417v3.pdf>}. 1071 1072 @arg AB8Ks: 8-Tuple of 8-th order I{Krüger} Alpha or Beta series 1073 coefficient tuples. 1074 1075 @return: I{Krüger} series coefficients (L{KsOrder}C{-tuple}). 1076 1077 @see: I{Karney}s 30-th order U{TMseries30 1078 <https://GeographicLib.SourceForge.io/html/tmseries30.html>}. 1079 ''' 1080 k = self.KsOrder 1081 if self.n: 1082 ns = fpowers(self.n, k) 1083 ks = tuple(fdot(AB8Ks[i][:k-i], *ns[i:]) for i in range(k)) 1084 else: 1085 ks = (_0_0,) * k 1086 return ks 1087 1088 @property_doc_(''' the I{Krüger} series' order (C{int}), see properties C{AlphaKs}, C{BetaKs}.''') 1089 def KsOrder(self): 1090 '''Get the I{Krüger} series' order (C{int} 4, 6 or 8). 1091 ''' 1092 return self._KsOrder 1093 1094 @KsOrder.setter # PYCHOK setter! 1095 def KsOrder(self, order): 1096 '''Set the I{Krüger} series' order. 1097 1098 @arg order: New I{Krüger} series' order (C{int} 4, 6 or 8). 1099 1100 @raise ValueError: Invalid B{C{order}}. 1101 ''' 1102 if not (isint(order) and order in (4, 6, 8)): 1103 raise _ValueError(order=order) 1104 if order != self._KsOrder: 1105 Ellipsoid.AlphaKs._update(self) 1106 Ellipsoid.BetaKs._update(self) 1107 self._KsOrder = order 1108 1109 @Property_RO 1110 def L(self): 1111 '''Get the I{quarter meridian} C{L}, aka C{polar distance}, the distance 1112 along a meridian between the equator and a pole (C{meter}), 1113 M{b * Elliptic(-e2 / (1 - e2)).E} or M{a * PI / 2}. 1114 ''' 1115 if self.f: # complete integral 2nd ... 1116 # kind: Elliptic(-e2 / (1 - e2)).E 1117 r = self._elliptic_e22.cE 1118 else: # spherical 1119 r = PI_2 1120 return Distance(L=self.b * r) 1121 1122 def Llat(self, lat): 1123 '''Return the I{meridional length}, the distance along a meridian 1124 between the equator and a (geodetic) latitude, see C{L}. 1125 1126 @arg lat: Geodetic latitude (C{degrees90}). 1127 1128 @return: The meridional length at B{C{lat}}, negative on southern 1129 hemisphere (C{meter}). 1130 ''' 1131 r = self._elliptic_e22.fEd(self.auxParametric(lat)) if self.f else Phi_(lat) 1132 return Distance(Llat=self.b * r) 1133 1134 Lmeridian = Llat # meridional distance 1135 1136 @Property_RO 1137 def Mabcd(self): 1138 '''Get the OSGR meridional coefficients (C{4-Tuple}), C{Airy130} only. 1139 ''' 1140 if self.n: 1141 n1, n2, n3 = fpowers(self.n, 3) # PYCHOK false! 1142 Mabcd = (fsum_(4, 4 * n1, 5 * n2, 5 * n3) / 4, 1143 fsum_( 24 * n1, 24 * n2, 21 * n3) / 8, 1144 fsum_( 15 * n2, 15 * n3) / 8, 1145 (35 * n3) / 24) 1146 else: 1147 Mabcd = _1_0, _0_0, _0_0, _0_0 1148 return Mabcd 1149 1150 @deprecated_Property_RO 1151 def majoradius(self): # PYCHOK no cover 1152 '''DEPRECATED, use property C{a} or C{Requatorial}.''' 1153 return self.a 1154 1155 def m2degrees(self, distance, lat=0): 1156 '''Convert a distance to an angle along the equator or 1157 along a parallel of (geodetic) latitude. 1158 1159 @arg distance: Distance (C{meter}). 1160 @kwarg lat: Parallel latitude (C{degrees90}, C{str}). 1161 1162 @return: Angle (C{degrees}) or C{INF} for near-polar B{C{lat}}. 1163 1164 @raise RangeError: Latitude B{C{lat}} outside valid range 1165 and L{rangerrors} set to C{True}. 1166 1167 @raise ValueError: Invalid B{C{distance}} or B{C{lat}}. 1168 ''' 1169 return degrees(self.m2radians(distance, lat=lat)) 1170 1171 def m2radians(self, distance, lat=0): 1172 '''Convert a distance to an angle along the equator or 1173 along a parallel of (geodetic) latitude. 1174 1175 @arg distance: Distance (C{meter}). 1176 @kwarg lat: Parallel latitude (C{degrees90}, C{str}). 1177 1178 @return: Angle (C{radians}) or C{INF} for near-polar B{C{lat}}. 1179 1180 @raise RangeError: Latitude B{C{lat}} outside valid range 1181 and L{rangerrors} set to C{True}. 1182 1183 @raise ValueError: Invalid B{C{distance}} or B{C{lat}}. 1184 ''' 1185 r = self.circle4(lat).radius if lat else self.a 1186 return m2radians(distance, radius=r, lat=0) 1187 1188 @deprecated_Property_RO 1189 def minoradius(self): # PYCHOK no cover 1190 '''DEPRECATED, use property C{b} or C{Rpolar}.''' 1191 return self.b 1192 1193 @Property_RO 1194 def n(self): 1195 '''Get the I{3rd flattening} (C{float}), M{f / (2 - f) == (a - b) / (a + b)}, see C{a_b2n}. 1196 ''' 1197 return self._assert(a_b2n(self.a, self.b), n=f2n(self.f)) 1198 1199 flattening3rd = n 1200 1201 @deprecated_Property_RO 1202 def quarteradius(self): # PYCHOK no cover 1203 '''DEPRECATED, use property C{L} or method C{Llat}.''' 1204 return self.L 1205 1206 @Property_RO 1207 def R1(self): 1208 '''Get the I{mean} earth radius per I{IUGG} (C{meter}), M{(2 * a + b) / 3}. 1209 1210 @see: U{Earth radius<https://WikiPedia.org/wiki/Earth_radius>} 1211 and method C{Rgeometric}. 1212 ''' 1213 r = (fsum_(self.a, self.a, self.b) / _3_0) if self.f else self.a 1214 return Radius(R1=r) 1215 1216 Rmean = R1 1217 1218 @Property_RO 1219 def R2(self): 1220 '''Get the I{authalic} earth radius (C{meter}), M{sqrt(c2)}. 1221 1222 @see: C{R2x}, C{c2}, C{area} and U{Earth radius 1223 <https://WikiPedia.org/wiki/Earth_radius>}. 1224 ''' 1225 return Radius(R2=sqrt(self.c2) if self.f else self.a) 1226 1227 Rauthalic = R2 1228 1229 @Property_RO 1230 def R2x(self): 1231 '''Get the I{authalic} earth radius (C{meter}), M{sqrt(c2x)}. 1232 1233 @see: C{R2}, C{c2x} and C{areax}. 1234 ''' 1235 return Radius(R2x=sqrt(self.c2x) if self.f else self.a) 1236 1237 Rauthalicx = R2x 1238 1239 @Property_RO 1240 def R3(self): 1241 '''Get the I{volumetric} earth radius (C{meter}), M{(a * a * b)**(1/3)}. 1242 1243 @see: U{Earth radius<https://WikiPedia.org/wiki/Earth_radius>} and C{volume}. 1244 ''' 1245 r = (cbrt(self.b_a) * self.a) if self.f else self.a 1246 return Radius(R3=r) 1247 1248 Rvolumetric = R3 1249 1250 def radians2m(self, rad, lat=0): 1251 '''Convert an angle to the distance along the equator or 1252 along a parallel of (geodetic) latitude. 1253 1254 @arg rad: The angle (C{radians}). 1255 @kwarg lat: Parallel latitude (C{degrees90}, C{str}). 1256 1257 @return: Distance (C{meter}, same units as the equatorial 1258 and polar radii) or C{0} for near-polar B{C{lat}}. 1259 1260 @raise RangeError: Latitude B{C{lat}} outside valid range 1261 and L{rangerrors} set to C{True}. 1262 1263 @raise ValueError: Invalid B{C{rad}} or B{C{lat}}. 1264 ''' 1265 r = self.circle4(lat).radius if lat else self.a 1266 return radians2m(rad, radius=r, lat=0) 1267 1268 @Property_RO 1269 def Rbiaxial(self): 1270 '''Get the I{biaxial, quadratic} mean earth radius (C{meter}), M{sqrt((a**2 + b**2) / 2)}. 1271 1272 @see: C{Rtriaxial} 1273 ''' 1274 b = (sqrt((_1_0 + self.b2_a2) * _0_5) * self.a) if self.f else self.a 1275 return Radius(Rbiaxial=b) 1276 1277 Requatorial = a # for consistent naming 1278 1279 def Rgeocentric(self, lat): 1280 '''Compute the I{geocentric} earth radius of (geodetic) latitude. 1281 1282 @arg lat: Latitude (C{degrees90}). 1283 1284 @return: Geocentric earth radius (C{meter}). 1285 1286 @raise ValueError: Invalid B{C{lat}}. 1287 1288 @see: U{Geocentric Radius 1289 <https://WikiPedia.org/wiki/Earth_radius#Geocentric_radius>} 1290 ''' 1291 r, a = self.a, Phi_(lat) 1292 if a and self.f: 1293 if abs(a) < PI_2: 1294 s2, c2 = _s2_c2(a) 1295 b2_a2_s2 = self.b2_a2 * s2 1296 # R == sqrt((a2**2 * c2 + b2**2 * s2) / (a2 * c2 + b2 * s2)) 1297 # == sqrt(a2**2 * (c2 + (b2 / a2)**2 * s2) / (a2 * (c2 + b2 / a2 * s2))) 1298 # == sqrt(a2 * (c2 + (b2 / a2)**2 * s2) / (c2 + (b2 / a2) * s2)) 1299 # == a * sqrt((c2 + b2_a2 * b2_a2 * s2) / (c2 + b2_a2 * s2)) 1300 # == a * sqrt((c2 + b2_a2 * b2_a2_s2) / (c2 + b2_a2_s2)) 1301 r *= sqrt((c2 + b2_a2_s2 * self.b2_a2) / (c2 + b2_a2_s2)) 1302 else: 1303 r = self.b 1304 return Radius(Rgeocentric=r) 1305 1306 @Property_RO 1307 def Rgeometric(self): 1308 '''Get the I{geometric} mean earth radius (C{meter}), M{sqrt(a * b)}. 1309 1310 @see: C{R1}. 1311 ''' 1312 g = sqrt(self.a * self.b) if self.f else self.a 1313 return Radius(Rgeometric=g) 1314 1315 def Rlat(self, lat): 1316 '''I{Approximate} the earth radius of (geodetic) latitude. 1317 1318 @arg lat: Latitude (C{degrees90}). 1319 1320 @return: Approximate earth radius (C{meter}). 1321 1322 @raise RangeError: Latitude B{C{lat}} outside valid range 1323 and L{rangerrors} set to C{True}. 1324 1325 @raise TypeError: Invalid B{C{lat}}. 1326 1327 @raise ValueError: Invalid B{C{lat}}. 1328 1329 @note: C{Rlat(B{90})} equals C{Rpolar}. 1330 1331 @see: Method C{Rparallel}. 1332 ''' 1333 # r = a - (a - b) * |lat| / 90 1334 r = self.a 1335 if self.f and lat: # .isEllipsoidal 1336 r -= (r - self.b) * abs(Lat(lat)) / _90_0 1337 r = Radius(Rlat=r) 1338 return r 1339 1340 Rpolar = b # for consistent naming 1341 1342 @deprecated_Property_RO 1343 def Rquadratic(self): # PYCHOK no cover 1344 '''DEPRECATED, use property C{Rbiaxial} or C{Rtriaxial}.''' 1345 return self.Rbiaxial 1346 1347 @deprecated_Property_RO 1348 def Rr(self): # PYCHOK no cover 1349 '''DEPRECATED, use property C{Rrectifying}.''' 1350 return self.Rrectifying 1351 1352 @Property_RO 1353 def Rrectifying(self): 1354 '''Get the I{rectifying} earth radius (C{meter}), M{((a**(3/2) + b**(3/2)) / 2)**(2/3)}. 1355 1356 @see: U{Earth radius<https://WikiPedia.org/wiki/Earth_radius>}. 1357 ''' 1358 r = (cbrt2((_1_0 + sqrt3(self.b_a)) * _0_5) * self.a) if self.f else self.a 1359 return Radius(Rrectifying=r) 1360 1361 def roc1_(self, sa, ca=None): 1362 '''Compute the I{prime-vertical}, I{normal} radius of curvature 1363 of (geodetic) latitude, I{unscaled}. 1364 1365 @arg sa: Sine of the latitude (C{float}, [-1.0..+1.0]). 1366 @kwarg ca: Optional cosine of the latitude (C{float}, [-1.0..+1.0]) 1367 to use an alternate formula. 1368 1369 @return: The prime-vertical radius of curvature (C{float}). 1370 1371 @note: The delta between both formulae with C{Ellipsoids.WGS84} 1372 is less than 2 nanometer over the entire latitude range. 1373 1374 @see: Method L{roc2_} and class L{EcefYou}. 1375 ''' 1376 if not self.f: # .isSpherical 1377 n = self.a 1378 elif ca is None: 1379 r = self.e2s2(sa) # see .roc2_ and _EcefBase._forward 1380 n = (self.a / sqrt(r)) if r > EPS02 else _0_0 1381 elif ca: # derived from EcefYou.forward 1382 h = hypot(ca, self.b_a * sa) if sa else abs(ca) 1383 n = self.a / h 1384 elif sa: 1385 n = self.a2_b / abs(sa) 1386 else: 1387 n = self.a 1388 return n 1389 1390 def roc2(self, lat, scaled=False): 1391 '''Compute the I{meridional} and I{prime-vertical}, I{normal} 1392 radii of curvature of (geodetic) latitude. 1393 1394 @arg lat: Latitude (C{degrees90}). 1395 @kwarg scaled: Scale prime_vertical by C{cos(radians(B{lat}))} (C{bool}). 1396 1397 @return: An L{Curvature2Tuple}C{(meridional, prime_vertical)} with 1398 the radii of curvature. 1399 1400 @raise ValueError: Invalid B{C{lat}}. 1401 1402 @see: Methods L{roc2_} and L{roc1_} and U{Local, flat earth 1403 approximation<https://www.EdWilliams.org/avform.htm#flat>} 1404 and meridional and prime vertical U{Radii of Curvature 1405 <https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>}. 1406 ''' 1407 return self.roc2_(Phi_(lat), scaled=scaled) 1408 1409 def roc2_(self, phi, scaled=False): 1410 '''Compute the I{meridional} and I{prime-vertical}, I{normal} 1411 radii of curvature of (geodetic) latitude. 1412 1413 @arg phi: Latitude (C{radians}). 1414 @kwarg scaled: Scale prime_vertical by C{cos(B{phi})} (C{bool}). 1415 1416 @return: An L{Curvature2Tuple}C{(meridional, prime_vertical)} with 1417 the radii of curvature. 1418 1419 @raise ValueError: Invalid B{C{phi}}. 1420 1421 @see: Methods L{roc2} and L{roc1_} and U{Local, flat earth 1422 approximation<https://www.EdWilliams.org/avform.htm#flat>} 1423 and the meridional and prime vertical U{Radii of Curvature 1424 <https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>}. 1425 ''' 1426 a = abs(Phi(phi)) 1427 if self.f: 1428 r = self.e2s2(sin(a)) 1429 if r > EPS02: 1430 n = self.a / sqrt(r) 1431 m = n * self.e12 / r # PYCHOK attr 1432 else: 1433 m = n = _0_0 # PYCHOK attr 1434 else: 1435 m = n = self.a 1436 if scaled and a: 1437 n *= cos(a) if a < PI_2 else _0_0 1438 return Curvature2Tuple(Radius(rocMeridional=m), 1439 Radius(rocPrimeVertical=n)) 1440 1441 def rocBearing(self, lat, bearing): 1442 '''Compute the I{directional} radius of curvature 1443 of (geodetic) latitude and compass direction. 1444 1445 @arg lat: Latitude (C{degrees90}). 1446 @arg bearing: Direction (compass C{degrees360}). 1447 1448 @return: Directional radius of curvature (C{meter}). 1449 1450 @raise RangeError: Latitude B{C{lat}} outside valid range 1451 and L{rangerrors} set to C{True}. 1452 1453 @raise ValueError: Invalid B{C{lat}} or B{C{bearing}}. 1454 1455 @see: U{Radii of Curvature 1456 <https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>} 1457 ''' 1458 if self.f: 1459 s2, c2 = _s2_c2(Bearing_(bearing)) 1460 m, n = self.roc2_(Phi_(lat)) 1461 if n < m: # == n / (c2 * n / m + s2) 1462 c2 *= n / m 1463 elif m < n: # == m / (c2 + s2 * m / n) 1464 s2 *= m / n 1465 n = m 1466 b = n / (c2 + s2) # == 1 / (c2 / m + s2 / n) 1467 else: 1468 b = self.b # == self.a 1469 return Radius(rocBearing=b) 1470 1471 def rocGauss(self, lat): 1472 '''Compute the I{Gaussian} radius of curvature of (geodetic) latitude. 1473 1474 @arg lat: Latitude (C{degrees90}). 1475 1476 @return: Gaussian radius of curvature (C{meter}). 1477 1478 @raise ValueError: Invalid B{C{lat}}. 1479 1480 @see: U{Radii of Curvature 1481 <https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>} 1482 ''' 1483 # using ... 1484 # m, n = self.roc2_(Phi_(lat)) 1485 # return sqrt(m * n) 1486 # ... requires 1 or 2 sqrt 1487 if self.f: 1488 s2, c2 = _s2_c2(Phi_(lat)) 1489 g = self.b / (c2 + self.b2_a2 * s2) 1490 else: 1491 g = self.b 1492 return Radius(rocGauss=g) 1493 1494 def rocMean(self, lat): 1495 '''Compute the I{mean} radius of curvature of (geodetic) latitude. 1496 1497 @arg lat: Latitude (C{degrees90}). 1498 1499 @return: Mean radius of curvature (C{meter}). 1500 1501 @raise ValueError: Invalid B{C{lat}}. 1502 1503 @see: U{Radii of Curvature 1504 <https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>} 1505 ''' 1506 if self.f: 1507 m, n = self.roc2_(Phi_(lat)) 1508 m *= 2 * n / (m + n) # == 2 / (1 / m + 1 / n) 1509 else: 1510 m = self.a 1511 return Radius(rocMean=m) 1512 1513 def rocMeridional(self, lat): 1514 '''Compute the I{meridional} radius of curvature of (geodetic) latitude. 1515 1516 @arg lat: Latitude (C{degrees90}). 1517 1518 @return: Meridional radius of curvature (C{meter}). 1519 1520 @raise ValueError: Invalid B{C{lat}}. 1521 1522 @see: Methods L{roc2} and L{roc2_} and U{Local, flat earth 1523 approximation<https://www.EdWilliams.org/avform.htm#flat>} 1524 and U{Radii of Curvature 1525 <https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>}. 1526 ''' 1527 return self.roc2_(Phi_(lat)).meridional 1528 1529 rocPolar = a2_b # synonymous 1530 1531 def rocPrimeVertical(self, lat): 1532 '''Compute the I{prime-vertical}, I{normal} radius of curvature 1533 of (geodetic) latitude, aka the transverse radius of curvature. 1534 1535 @arg lat: Latitude (C{degrees90}). 1536 1537 @return: Prime-vertical radius of curvature (C{meter}). 1538 1539 @raise ValueError: Invalid B{C{lat}}. 1540 1541 @see: Methods L{roc2}, L{roc2_} and L{roc1_} and U{Local, flat earth 1542 approximation<https://www.EdWilliams.org/avform.htm#flat>} 1543 and U{Radii of Curvature 1544 <https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>}. 1545 ''' 1546 return self.roc2_(Phi_(lat)).prime_vertical 1547 1548 rocTransverse = rocPrimeVertical # synonyms 1549 1550 @deprecated_Property_RO 1551 def Rs(self): # PYCHOK no cover 1552 '''DEPRECATED, use property C{Rgeometric}.''' 1553 return self.Rgeometric 1554 1555 @Property_RO 1556 def Rtriaxial(self): 1557 '''Get the I{triaxial, quadratic} mean earth radius (C{meter}), M{sqrt((3 * a**2 + b**2) / 4)}. 1558 1559 @see: C{Rbiaxial} 1560 ''' 1561 t = (sqrt((_3_0 + self.b2_a2) * _0_25) * self.a) if self.f else self.a 1562 return Radius(Rtriaxial=t) 1563 1564 def toStr(self, prec=8): # PYCHOK expected 1565 '''Return this ellipsoid as a text string. 1566 1567 @kwarg prec: Optional number of decimals, unstripped (C{int}). 1568 1569 @return: Ellipsoid attributes (C{str}). 1570 ''' 1571 E = Ellipsoid 1572 return self._instr(prec, _a_, E.b.name, E.f_.name, _f_, E.f2.name, E.n.name, 1573 E.e.name, E.e2.name, E.e22.name, E.e32.name, 1574 E.A.name, E.L.name, E.R1.name, E.R2.name, E.R3.name, 1575 E.Rbiaxial.name, E.Rtriaxial.name) 1576 1577 @Property_RO 1578 def volume(self): 1579 '''Get the ellipsoid's I{volume} (C{meter**3}), M{4 / 3 * PI * R3**3}. 1580 1581 @see: C{R3}. 1582 ''' 1583 return Meter3(volume=self.a2 * self.b * (PI4 / _3_0)) 1584 1585 1586class Ellipsoid2(Ellipsoid): 1587 '''An L{Ellipsoid} specified by I{equatorial} radius and I{flattening}. 1588 ''' 1589 def __init__(self, a, f, name=NN): 1590 '''New L{Ellipsoid2}. 1591 1592 @arg a: Equatorial radius, semi-axis (C{meter}). 1593 @arg f: Flattening: (C{float} < 1.0, negative for I{prolate}). 1594 @kwarg name: Optional, unique name (C{str}). 1595 1596 @raise NameError: Ellipsoid with that B{C{name}} already exists. 1597 1598 @raise ValueError: Invalid B{C{a}} or B{C{f}}. 1599 1600 @note: C{abs(B{f}) < EPS} is forced to C{B{f}=0}, I{spherical}. 1601 Negative C{B{f}} produces a I{prolate} ellipsoid. 1602 ''' 1603 try: 1604 a = Radius_(a=a, low=EPS, high=None) # like Radius_ 1605 f = Float_( f=f, low=None, high=EPS1) 1606 except (TypeError, ValueError) as x: 1607 t = instr(self, a=a, f=f, name=name) # PYCHOK no cover 1608 raise _ValueError(t, txt=str(x)) 1609 Ellipsoid.__init__(self, a, a_f2b(a, f), name=name) 1610 1611 1612def _spherical_a_b(a, b): 1613 '''(INTERNAL) C{True} for spherical or invalid C{a} or C{b}. 1614 ''' 1615 return a < EPS or b < EPS or abs(a - b) < EPS 1616 1617 1618def _spherical_f(f): 1619 '''(INTERNAL) C{True} for spherical or invalid C{f}. 1620 ''' 1621 return abs(f) < EPS or f > EPS1 1622 1623 1624def _spherical_f_(f_): 1625 '''(INTERNAL) C{True} for spherical or invalid C{f_}. 1626 ''' 1627 return abs(f_) < EPS or abs(f_) > _1_EPS 1628 1629 1630def a_b2e(a, b): 1631 '''Return C{e}, the I{1st eccentricity} for a given I{equatorial} and I{polar} radius. 1632 1633 @arg a: Equatorial radius (C{scalar} > 0). 1634 @arg b: Polar radius (C{scalar} > 0). 1635 1636 @return: The (1st) eccentricity (C{float} or C{0}), M{sqrt(1 - (b / a)**2)}. 1637 1638 @note: The result is always non-negative and C{0} for I{near-spherical} ellipsoids. 1639 ''' 1640 return Float(e=sqrt(abs(a_b2e2(a, b)))) 1641 1642 1643def a_b2e2(a, b): 1644 '''Return C{e2}, the I{1st eccentricity squared} for a given I{equatorial} and I{polar} radius. 1645 1646 @arg a: Equatorial radius (C{scalar} > 0). 1647 @arg b: Polar radius (C{scalar} > 0). 1648 1649 @return: The (1st) eccentricity I{squared} (C{float} or C{0}), M{1 - (b / a)**2}. 1650 1651 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0} 1652 for I{near-spherical} ellipsoids. 1653 ''' 1654 return Float(e2=_0_0 if _spherical_a_b(a, b) else (_1_0 - (b / a)**2)) 1655 1656 1657def a_b2e22(a, b): 1658 '''Return C{e22}, the I{2nd eccentricity squared} for a given I{equatorial} and I{polar} radius. 1659 1660 @arg a: Equatorial radius (C{scalar} > 0). 1661 @arg b: Polar radius (C{scalar} > 0). 1662 1663 @return: The 2nd eccentricity I{squared} (C{float} or C{0}), M{(a / b)**2 - 1}. 1664 1665 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0} 1666 for I{near-spherical} ellipsoids. 1667 ''' 1668 return Float(e22=_0_0 if _spherical_a_b(a, b) else ((a / b)**2 - _1_0)) 1669 1670 1671def a_b2e32(a, b): 1672 '''Return C{e32}, the I{3rd eccentricity squared} for a given I{equatorial} and I{polar} radius. 1673 1674 @arg a: Equatorial radius (C{scalar} > 0). 1675 @arg b: Polar radius (C{scalar} > 0). 1676 1677 @return: The 3rd eccentricity I{squared} (C{float} or C{0}), M{(a**2 - b**2) / (a**2 + b**2)}. 1678 1679 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0} 1680 for I{near-spherical} ellipsoids. 1681 ''' 1682 a2, b2 = a**2, b**2 1683 return Float(e32=_0_0 if _spherical_a_b(a2, b2) else ((a2 - b2) / (a2 + b2))) 1684 1685 1686def a_b2f(a, b): 1687 '''Return C{f}, the I{flattening} for a given I{equatorial} and I{polar} radius. 1688 1689 @arg a: Equatorial radius (C{scalar} > 0). 1690 @arg b: Polar radius (C{scalar} > 0). 1691 1692 @return: The flattening (C{float} or C{0}), M{(a - b) / a}. 1693 1694 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0} 1695 for I{near-spherical} ellipsoids. 1696 ''' 1697 f = 0 if _spherical_a_b(a, b) else ((a - b) / a) 1698 return _f_0_0 if _spherical_f(f) else Float(f=f) 1699 1700 1701def a_b2f_(a, b): 1702 '''Return C{f_}, the I{inverse flattening} for a given I{equatorial} and I{polar} radius. 1703 1704 @arg a: Equatorial radius (C{scalar} > 0). 1705 @arg b: Polar radius (C{scalar} > 0). 1706 1707 @return: The inverse flattening (C{float} or C{0}), M{a / (a - b)}. 1708 1709 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0} 1710 for I{near-spherical} ellipsoids. 1711 ''' 1712 f_ = 0 if _spherical_a_b(a, b) else (a / float(a - b)) 1713 return _f__0_0 if _spherical_f_(f_) else Float(f_=f_) 1714 1715 1716def a_b2f2(a, b): 1717 '''Return C{f2}, the I{2nd flattening} for a given I{equatorial} and I{polar} radius. 1718 1719 @arg a: Equatorial radius (C{scalar} > 0). 1720 @arg b: Polar radius (C{scalar} > 0). 1721 1722 @return: The 2nd flattening (C{float} or C{0}), M{(a - b) / b}. 1723 1724 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0} 1725 for I{near-spherical} ellipsoids. 1726 ''' 1727 t = 0 if _spherical_a_b(a, b) else float(a - b) 1728 return Float(f2=_0_0 if abs(t) < EPS else (t / b)) 1729 1730 1731def a_b2n(a, b): 1732 '''Return C{n}, the I{3rd flattening} for a given I{equatorial} and I{polar} radius. 1733 1734 @arg a: Equatorial radius (C{scalar} > 0). 1735 @arg b: Polar radius (C{scalar} > 0). 1736 1737 @return: The 3rd flattening (C{float} or C{0}), M{(a - b) / (a + b)}. 1738 1739 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0} 1740 for I{near-spherical} ellipsoids. 1741 ''' 1742 t = 0 if _spherical_a_b(a, b) else float(a - b) 1743 return Float(n=_0_0 if abs(t) < EPS else (t / (a + b))) 1744 1745 1746def a_f2b(a, f): 1747 '''Return C{b}, the I{polar} radius for a given I{equatorial} radius and I{flattening}. 1748 1749 @arg a: Equatorial radius (C{scalar} > 0). 1750 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}). 1751 1752 @return: The polar radius (C{float}), M{a * (1 - f)}. 1753 ''' 1754 b = a if _spherical_f(f) else (a * (_1_0 - f)) 1755 return Radius_(b=a if _spherical_a_b(a, b) else b) 1756 1757 1758def a_f_2b(a, f_): 1759 '''Return C{b}, the I{polar} radius for a given I{equatorial} radius and I{inverse flattening}. 1760 1761 @arg a: Equatorial radius (C{scalar} > 0). 1762 @arg f_: Inverse flattening (C{scalar} >>> 1). 1763 1764 @return: The polar radius (C{float}), M{a * (f_ - 1) / f_}. 1765 ''' 1766 b = a if _spherical_f_(f_) else (a * (f_ - _1_0) / f_) 1767 return Radius_(b=a if _spherical_a_b(a, b) else b) 1768 1769 1770def b_f2a(b, f): 1771 '''Return C{a}, the I{equatorial} radius for a given I{polar} radius and I{flattening}. 1772 1773 @arg b: Polar radius (C{scalar} > 0). 1774 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}). 1775 1776 @return: The equatorial radius (C{float}), M{b / (1 - f)}. 1777 ''' 1778 t = _1_0 - f 1779 a = b if abs(t < EPS) else (b / t) 1780 return Radius_(a=b if _spherical_a_b(a, b) else a) 1781 1782 1783def b_f_2a(b, f_): 1784 '''Return C{a}, the I{equatorial} radius for a given I{polar} radius and I{inverse flattening}. 1785 1786 @arg b: Polar radius (C{scalar} > 0). 1787 @arg f_: Inverse flattening (C{scalar} >>> 1). 1788 1789 @return: The equatorial radius (C{float}), M{b * f_ / (f_ - 1)}. 1790 ''' 1791 t = f_ - _1_0 1792 a = b if _spherical_f_(f_) or abs(t - f_) < EPS \ 1793 or abs(t) < EPS else (b * f_ / t) 1794 return Radius_(a=b if _spherical_a_b(a, b) else a) 1795 1796 1797def f2e2(f): 1798 '''Return C{e2}, the I{1st eccentricity squared} for a given I{flattening}. 1799 1800 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}). 1801 1802 @return: The (1st) eccentricity I{squared} (C{float} < 1), M{f * (2 - f)}. 1803 1804 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0} 1805 for I{near-spherical} ellipsoids. 1806 1807 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/ 1808 html/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening 1809 <https://WikiPedia.org/wiki/Flattening>}. 1810 ''' 1811 return Float(e2=_0_0 if _spherical_f(f) else (f * (_2_0 - f))) 1812 1813 1814def f2e22(f): 1815 '''Return C{e22}, the I{2nd eccentricity squared} for a given I{flattening}. 1816 1817 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}). 1818 1819 @return: The 2nd eccentricity I{squared} (C{float} > -1 or C{INF}), 1820 M{f * (2 - f) / (1 - f)**2}. 1821 1822 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0} 1823 for near-spherical ellipsoids. 1824 1825 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/ 1826 html/classGeographicLib_1_1Ellipsoid.html>}. 1827 ''' 1828 # e2 / (1 - e2) == f * (2 - f) / (1 - f)**2 1829 t = (_1_0 - f)**2 1830 return Float(e22=INF if t < EPS else (f2e2(f) / t)) # PYCHOK type 1831 1832 1833def f2e32(f): 1834 '''Return C{e32}, the I{3rd eccentricity squared} for a given I{flattening}. 1835 1836 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}). 1837 1838 @return: The 3rd eccentricity I{squared} (C{float}), M{f * (2 - f) / (1 + (1 - f)**2)}. 1839 1840 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0} 1841 for I{near-spherical} ellipsoids. 1842 1843 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/ 1844 html/classGeographicLib_1_1Ellipsoid.html>}. 1845 ''' 1846 # e2 / (2 - e2) == f * (2 - f) / (1 + (1 - f)**2) 1847 e2 = f2e2(f) 1848 return Float(e32=e2 / (_2_0 - e2)) 1849 1850 1851def f_2f(f_): 1852 '''Return C{f}, the I{flattening} for a given I{inverse flattening}. 1853 1854 @arg f_: Inverse flattening (C{scalar} >>> 1). 1855 1856 @return: The flattening (C{float} or C{0}), M{1 / f_}. 1857 1858 @note: The result is positive for I{oblate}, negative for I{prolate} 1859 or C{0} for I{near-spherical} ellipsoids. 1860 ''' 1861 f = 0 if _spherical_f_(f_) else _1_0 / f_ 1862 return _f_0_0 if _spherical_f(f) else Float(f=f) # PYCHOK type 1863 1864 1865def f2f_(f): 1866 '''Return C{f_}, the I{inverse flattening} for a given I{flattening}. 1867 1868 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}). 1869 1870 @return: The inverse flattening (C{float} or C{0}), M{1 / f}. 1871 1872 @note: The result is positive for I{oblate}, negative for I{prolate} 1873 or C{0} for I{near-spherical} ellipsoids. 1874 ''' 1875 f_ = 0 if _spherical_f(f) else _1_0 / f 1876 return _f__0_0 if _spherical_f_(f_) else Float(f_=f_) # PYCHOK type 1877 1878 1879def f2f2(f): 1880 '''Return C{f2}, the I{2nd flattening} for a given I{flattening}. 1881 1882 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}). 1883 1884 @return: The 2nd flattening (C{float} or C{INF}), M{f / (1 - f)}. 1885 1886 @note: The result is positive for I{oblate}, negative for I{prolate} 1887 or C{0} for I{near-spherical} ellipsoids. 1888 1889 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/ 1890 html/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening 1891 <https://WikiPedia.org/wiki/Flattening>}. 1892 ''' 1893 t = _1_0 - f 1894 return Float(f2=_0_0 if _spherical_f(f) else 1895 (INF if abs(t) < EPS else (f / t))) # PYCHOK type 1896 1897 1898def f2n(f): 1899 '''Return C{n}, the I{3rd flattening} for a given I{flattening}. 1900 1901 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}). 1902 1903 @return: The 3rd flattening (-1 <= C{float} < 1), M{f / (2 - f)}. 1904 1905 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0} 1906 for I{near-spherical} ellipsoids. 1907 1908 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/ 1909 html/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening 1910 <https://WikiPedia.org/wiki/Flattening>}. 1911 ''' 1912 return Float(n=_0_0 if _spherical_f(f) else (f / float(_2_0 - f))) 1913 1914 1915def n2e2(n): 1916 '''Return C{e2}, the I{1st eccentricity squared} for a given I{3rd flattening}. 1917 1918 @arg n: The 3rd flattening (-1 <= C{scalar} < 1). 1919 1920 @return: The (1st) eccentricity I{squared} (C{float} or -INF), M{4 * n / (1 + n)**2}. 1921 1922 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0} 1923 for I{near-spherical} ellipsoids. 1924 1925 @see: U{Flattening<https://WikiPedia.org/wiki/Flattening>}. 1926 ''' 1927 t = (_1_0 + n)**2 1928 return Float(e2=_0_0 if abs(n) < EPS else 1929 (-INF if t < EPS else (_4_0 * n / t))) 1930 1931 1932def _normal2(px, py, E): 1933 '''(INTERNAL) Nearest point on a 2-D ellipse. 1934 ''' 1935 a, b = E.a, E.b 1936 if min(px, py, a, b) < EPS0: 1937 raise _AssertionError(px=px, py=py, a=a, b=b, E=E) 1938 1939 a2 = a - b * E.b_a 1940 b2 = b - a * E.a_b 1941 tx = ty = 0.70710678118 1942 for _ in range(9): # 4..5 trips max 1943 ex = a2 * tx**3 1944 ey = b2 * ty**3 1945 1946 qx = px - ex 1947 qy = py - ey 1948 q = hypot(qx, qy) 1949 if q < EPS0: 1950 break 1951 r = hypot(ex - tx * a, ey - ty * b) / q 1952 1953 sx, tx = tx, min(_1_0, max(0, (ex + qx * r) / a)) 1954 sy, ty = ty, min(_1_0, max(0, (ey + qy * r) / b)) 1955 t = hypot(ty, tx) 1956 if t < EPS0: 1957 break 1958 tx = tx / t # /= t chokes PyChecker 1959 ty = ty / t 1960 if max(abs(sx - tx), abs(sy - ty)) < EPS: 1961 break 1962 1963 tx *= a / px 1964 ty *= b / py 1965 return tx, ty # x and y as fractions 1966 1967 1968def n2f(n): 1969 '''Return C{f}, the I{flattening} for a given I{3rd flattening}. 1970 1971 @arg n: The 3rd flattening (-1 <= C{scalar} < 1). 1972 1973 @return: The flattening (C{float} or -INF), M{2 * n / (1 + n)}. 1974 1975 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/ 1976 html/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening 1977 <https://WikiPedia.org/wiki/Flattening>}. 1978 ''' 1979 t = n + _1_0 1980 f = 0 if abs(n) < EPS else (-INF if t < EPS else _2_0 * n / t) 1981 return _f_0_0 if _spherical_f(f) else Float(f=f) 1982 1983 1984def n2f_(n): 1985 '''Return C{f_}, the I{inverse flattening} for a given I{3rd flattening}. 1986 1987 @arg n: The 3rd flattening (-1 <= C{scalar} < 1). 1988 1989 @return: The inverse flattening (C{float} or C{0}), M{1 / f}. 1990 1991 @see: L{n2f} and L{f2f_}. 1992 ''' 1993 return f2f_(n2f(n)) 1994 1995 1996class Ellipsoids(_NamedEnum): 1997 '''(INTERNAL) L{Ellipsoid} registry, I{must} be a sub-class 1998 to accommodate the L{_LazyNamedEnumItem} properties. 1999 ''' 2000 def _Lazy(self, a, b, f_, **kwds): 2001 '''(INTERNAL) Instantiate the L{Ellipsoid}. 2002 ''' 2003 return Ellipsoid(a, b, f_, **kwds) 2004 2005Ellipsoids = Ellipsoids(Ellipsoid) # PYCHOK singleton 2006'''Some pre-defined L{Ellipsoid}s, all I{lazily} instantiated.''' 2007# <https://www.GNU.org/software/gama/manual/html_node/Supported-ellipsoids.html> 2008# <https://w3.Energistics.org/archive/Epicentre/Epicentre_v3.0/DataModel/ 2009# LogicalDictionary/StandardValues/ellipsoid.html> 2010# <https://kb.OSU.edu/dspace/handle/1811/77986> 2011Ellipsoids._assert( # <https://WikiPedia.org/wiki/Earth_ellipsoid> 2012 Airy1830 = _lazy(_Airy1830_, *_T(6377563.396, _0_0, 299.3249646)), # b=6356256.909 2013 AiryModified = _lazy(_AiryModified_, *_T(6377340.189, _0_0, 299.3249646)), # b=6356034.448 2014# ANS = _lazy('ANS', *_T(6378160.0, _0_0, 298.25)), # Australian Nat. Spheroid 2015 Australia1966 = _lazy('Australia1966', *_T(6378160.0, _0_0, 298.25)), # b=6356774.7192 2016# Bessel1841 = _lazy(_Bessel1841_, *_T(6377397.155, 6356078.963, 299.152815351)), 2017 Bessel1841 = _lazy(_Bessel1841_, *_T(6377397.155, 6356078.962818, 299.152812797)), 2018 Clarke1866 = _lazy(_Clarke1866_, *_T(6378206.4, 6356583.8, 294.978698214)), 2019 Clarke1880 = _lazy('Clarke1880', *_T(6378249.145, 6356514.86954978, 293.465)), 2020 Clarke1880IGN = _lazy(_Clarke1880IGN_, *_T(6378249.2, 6356515.0, 293.466021294)), 2021 Clarke1880Mod = _lazy('Clarke1880Mod', *_T(6378249.145, 6356514.96582849, 293.4663)), 2022 CPM1799 = _lazy('CPM1799', *_T(6375738.7, 6356671.92557493, 334.39)), # Comm. des Poids et Mesures 2023 Delambre1810 = _lazy('Delambre1810', *_T(6376428.0, 6355957.92616372, 311.5)), # Belgium 2024 Engelis1985 = _lazy('Engelis1985', *_T(6378136.05, 6356751.32272154, 298.2566)), 2025 Everest1969 = _lazy('Everest1969', *_T(6377295.664, 6356094.667915, 300.801699997)), 2026 Fisher1968 = _lazy('Fisher1968', *_T(6378150.0, 6356768.33724438, 298.3)), 2027 GEM10C = _lazy('GEM10C', *_T(R_MA, 6356752.31424783, 298.2572236)), 2028# GPES = _lazy('GPES', *_T(6378135.0, 6356750.0, _0_0)), # "Gen. Purpose Earth Spheroid" 2029 GRS67 = _lazy('GRS67', *_T(6378160.0, _0_0, 298.247167427)), # Lucerne b=6356774.516 2030 GRS80 = _lazy(_GRS80_, *_T(R_MA, 6356752.314140347, 298.257222101)), # ITRS, ETRS89 2031# Hayford1924 = _lazy('Hayford1924', *_T(6378388.0, 6356911.94612795, _0_0)), # aka Intl1924 f_=297 2032 Helmert1906 = _lazy('Helmert1906', *_T(6378200.0, 6356818.16962789, 298.3)), 2033 IERS1989 = _lazy('IERS1989', *_T(6378136.0, _0_0, 298.257)), # b=6356751.302 2034 IERS1992TOPEX = _lazy('IERS1992TOPEX', *_T(6378136.3, 6356751.61659215, 298.257223563)), # IERS/TOPEX/Poseidon/McCarthy 2035 IERS2003 = _lazy('IERS2003', *_T(6378136.6, 6356751.85797165, 298.25642)), 2036 Intl1924 = _lazy(_Intl1924_, *_T(6378388.0, _0_0, 297.0)), # aka Hayford b=6356911.9462795 2037 Intl1967 = _lazy('Intl1967', *_T(6378157.5, 6356772.2, 298.24961539)), # New Int'l 2038 Krassovski1940 = _lazy(_Krassovski1940_, *_T(6378245.0, 6356863.01877305, 298.3)), # spelling 2039 Krassowsky1940 = _lazy(_Krassowsky1940_, *_T(6378245.0, 6356863.01877305, 298.3)), # spelling 2040 Maupertuis1738 = _lazy('Maupertuis1738', *_T(6397300.0, 6363806.28272251, 191.0)), # France 2041 Mercury1960 = _lazy('Mercury1960', *_T(6378166.0, 6356784.28360711, 298.3)), 2042 Mercury1968Mod = _lazy('Mercury1968Mod', *_T(6378150.0, 6356768.33724438, 298.3)), 2043 NWL1965 = _lazy('NWL1965', *_T(6378145.0, 6356759.76948868, 298.25)), # Naval Weapons Lab. 2044 OSU86F = _lazy('OSU86F', *_T(6378136.2, 6356751.51693008, 298.2572236)), 2045 OSU91A = _lazy('OSU91A', *_T(6378136.3, 6356751.6165948, 298.2572236)), 2046# Plessis1817 = _lazy('Plessis1817', *_T(6397523.0, 6355863.0, 153.56512242)), # XXX incorrect? 2047 Plessis1817 = _lazy('Plessis1817', *_T(6376523.0, 6355862.93325557, 308.64)), # XXX IGN France 1972 2048 SGS85 = _lazy('SGS85', *_T(6378136.0, 6356751.30156878, 298.257)), # Soviet Geodetic System 2049 SoAmerican1969 = _lazy('SoAmerican1969', *_T(6378160.0, 6356774.71919531, 298.25)), # South American 2050 Struve1860 = _lazy('Struve1860', *_T(6378298.3, 6356657.14266956, 294.73)), 2051 WGS60 = _lazy('WGS60', *_T(6378165.0, 6356783.28695944, 298.3)), 2052 WGS66 = _lazy('WGS66', *_T(6378145.0, 6356759.76948868, 298.25)), 2053 WGS72 = _lazy(_WGS72_, *_T(6378135.0, _0_0, 298.26)), # b=6356750.52 2054 WGS84 = _lazy(_WGS84_, *_T(R_MA, _0_0, _f_WGS84)), # GPS b=6356752.3142451793 2055# Prolate = _lazy('Prolate', *_T(6356752.3, R_MA, _0_0)), 2056 Sphere = _lazy(_Sphere_, *_T(R_M, R_M, _0_0)), # pseudo 2057 SphereAuthalic = _lazy('SphereAuthalic', *_T(R_FM, R_FM, _0_0)), # pseudo 2058 SpherePopular = _lazy('SpherePopular', *_T(R_MA, R_MA, _0_0)) # EPSG:3857 Spheroid 2059) 2060 2061 2062if __name__ == '__main__': 2063 2064 from pygeodesy.interns import _COMMA_, _NL_, _NL_hash_, _NL_var_ 2065 from pygeodesy.named import nameof 2066 2067 for E in (Ellipsoids.WGS84, Ellipsoids.GRS80, # NAD83, 2068 Ellipsoids.Sphere, Ellipsoids.SpherePopular, 2069 Ellipsoid(Ellipsoids.WGS84.b, Ellipsoids.WGS84.a, name='_Prolate')): 2070 e = f2n(E.f) - E.n 2071 t = NN(_COMMA_, _NL_hash_, _SPACE_)(E.toStr(prec=10), # re-callable 2072 'e=%s, f_=%s, f=%s, n=%s (%s)' % (fstr(E.e, prec=13, fmt=Fmt.e), 2073 fstr(E.f_, prec=13, fmt=Fmt.e), 2074 fstr(E.f, prec=13, fmt=Fmt.e), 2075 fstr(E.n, prec=13, fmt=Fmt.e), 2076 fstr(e, prec=9, fmt=Fmt.e),), 2077 '%s=(%s)' % (Ellipsoid.AlphaKs.name, fstr(E.AlphaKs, prec=20),), 2078 '%s= (%s)' % (Ellipsoid.BetaKs.name, fstr(E.BetaKs, prec=20),), 2079 '%s= %s' % (nameof(Ellipsoid.KsOrder), E.KsOrder), # property 2080 '%s= (%s)' % (Ellipsoid.Mabcd.name, fstr(E.Mabcd, prec=20),)) 2081 print('%s%s: %s' % (_NL_hash_, _DOT_(E.classname, E.name), t)) 2082 2083 # __doc__ of this file, force all into registry 2084 t = [NN] + Ellipsoids.toRepr(all=True).split(_NL_) 2085 print(_NL_var_.join(i.strip(_COMMA_) for i in t)) 2086 2087# **) MIT License 2088# 2089# Copyright (C) 2016-2021 -- mrJean1 at Gmail -- All Rights Reserved. 2090# 2091# Permission is hereby granted, free of charge, to any person obtaining a 2092# copy of this software and associated documentation files (the "Software"), 2093# to deal in the Software without restriction, including without limitation 2094# the rights to use, copy, modify, merge, publish, distribute, sublicense, 2095# and/or sell copies of the Software, and to permit persons to whom the 2096# Software is furnished to do so, subject to the following conditions: 2097# 2098# The above copyright notice and this permission notice shall be included 2099# in all copies or substantial portions of the Software. 2100# 2101# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 2102# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 2103# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 2104# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR 2105# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, 2106# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR 2107# OTHER DEALINGS IN THE SOFTWARE. 2108 2109# % python3 -m pygeodesy.ellipsoids 2110 2111# Ellipsoid.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, e22=0.0067394967, e32=0.0033584313, A=6367449.1458234144, L=10001965.7293127235, R1=6371008.7714150595, R2=6371007.1809184747, R3=6371000.7900091587, Rbiaxial=6367453.6345163295, Rtriaxial=6372797.5559594007, 2112# e=8.1819190842622e-02, f_=2.98257223563e+02, f=3.3528106647475e-03, n=1.6792203863837e-03 (0.0e+00), 2113# AlphaKs=(0.00083773182062447786, 0.00000076085277735726, 0.00000000119764550324, 0.00000000000242917068, 0.00000000000000571182, 0.0000000000000000148, 0.00000000000000000004, 0.0), 2114# BetaKs= (0.00083773216405795667, 0.0000000590587015222, 0.00000000016734826653, 0.00000000000021647981, 0.00000000000000037879, 0.00000000000000000072, 0.0, 0.0), 2115# KsOrder= 8, 2116# Mabcd= (1.00168275103155868244, 0.00504613293193333871, 0.00000529596776243457, 0.00000000690525779769) 2117 2118# Ellipsoid.GRS80: name='GRS80', a=6378137, b=6356752.3141403468, f_=298.257222101, f=0.0033528107, f2=0.0033640898, n=0.0016792204, e=0.081819191, e2=0.00669438, e22=0.0067394968, e32=0.0033584313, A=6367449.1457710434, L=10001965.7292304579, R1=6371008.7713801153, R2=6371007.1808835147, R3=6371000.7899741363, Rbiaxial=6367453.6344640013, Rtriaxial=6372797.5559332585, 2119# e=8.1819191042833e-02, f_=2.98257222101e+02, f=3.3528106811837e-03, n=1.6792203946295e-03 (0.0e+00), 2120# AlphaKs=(0.00083773182472890429, 0.00000076085278481561, 0.00000000119764552086, 0.00000000000242917073, 0.00000000000000571182, 0.0000000000000000148, 0.00000000000000000004, 0.0), 2121# BetaKs= (0.0008377321681623882, 0.00000005905870210374, 0.000000000167348269, 0.00000000000021647982, 0.00000000000000037879, 0.00000000000000000072, 0.0, 0.0), 2122# KsOrder= 8, 2123# Mabcd= (1.00168275103983916985, 0.0050461329567537995, 0.00000529596781448937, 0.00000000690525789941) 2124 2125# Ellipsoid.Sphere: name='Sphere', a=6371008.7714149999, b=6371008.7714149999, f_=0, f=0, f2=0, n=0, e=0, e2=0, 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, 2126# e=0.0e+00, f_=0.0e+00, f=0.0e+00, n=0.0e+00 (0.0e+00), 2127# AlphaKs=(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0), 2128# BetaKs= (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0), 2129# KsOrder= 8, 2130# Mabcd= (1.0, 0.0, 0.0, 0.0) 2131 2132# Ellipsoid.SpherePopular: name='SpherePopular', a=6378137, b=6378137, f_=0, f=0, f2=0, n=0, e=0, e2=0, e22=0, e32=0, A=6378137, L=10018754.171394622, R1=6378137, R2=6378137, R3=6378137, Rbiaxial=6378137, Rtriaxial=6378137, 2133# e=0.0e+00, f_=0.0e+00, f=0.0e+00, n=0.0e+00 (0.0e+00), 2134# AlphaKs=(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0), 2135# BetaKs= (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0), 2136# KsOrder= 8, 2137# Mabcd= (1.0, 0.0, 0.0, 0.0) 2138 2139# Ellipsoid._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, e22=-0.00669438, e32=-0.0033584313, A=6367449.1458234154, L=10035500.5204500332, R1=6363880.5428301198, R2=6363878.9413582645, R3=6363872.5644020075, Rbiaxial=6367453.6345163304, Rtriaxial=6362105.2243882548, 2140# e=8.2094437949696e-02, f_=-2.97257223563e+02, f=-3.3640898209765e-03, n=-1.6792203863837e-03 (0.0e+00), 2141# AlphaKs=(-0.00084149152514366627, 0.00000076653480614871, -0.00000000120934503389, 0.0000000000024576225, -0.00000000000000578863, 0.00000000000000001502, -0.00000000000000000004, 0.0), 2142# BetaKs= (-0.00084149187224351817, 0.00000005842735196773, -0.0000000001680487236, 0.00000000000021706261, -0.00000000000000038002, 0.00000000000000000073, -0.0, 0.0), 2143# KsOrder= 8, 2144# Mabcd= (0.99832429842120640195, -0.00502921424529705757, 0.00000527821138524052, -0.00000000690525779769) 2145