Coverage for pygeodesy/ellipsoids.py: 98%
692 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-11 14:35 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-11 14:35 -0400
2# -*- coding: utf-8 -*-
4u'''Ellipsoidal and spherical earth models.
6Classes L{a_f2Tuple}, L{Ellipsoid} and L{Ellipsoid2}, an L{Ellipsoids} registry and
72 dozen functions to convert I{equatorial} radius, I{polar} radius, I{eccentricities},
8I{flattenings} and I{inverse flattening}.
10See module L{datums} for L{Datum} and L{Transform} information and other details.
12Following is the list of predefined L{Ellipsoid}s, all instantiated lazily.
14@var Ellipsoids.Airy1830: Ellipsoid(name='Airy1830', a=6377563.396, b=6356256.90923729, f_=299.3249646, f=0.00334085, f2=0.00335205, n=0.00167322, e=0.08167337, e2=0.00667054, e22=0.00671533, e32=0.00334643, A=6366914.60892522, L=10001126.0807165, R1=6370461.23374576, R2=6370459.65470808, R3=6370453.30994572, Rbiaxial=6366919.065224, Rtriaxial=6372243.45317691)
15@var Ellipsoids.AiryModified: Ellipsoid(name='AiryModified', a=6377340.189, b=6356034.44793853, f_=299.3249646, f=0.00334085, f2=0.00335205, n=0.00167322, e=0.08167337, e2=0.00667054, e22=0.00671533, e32=0.00334643, A=6366691.77461988, L=10000776.05340819, R1=6370238.27531284, R2=6370236.69633043, R3=6370230.35179013, Rbiaxial=6366696.2307627, Rtriaxial=6372020.43236847)
16@var Ellipsoids.ATS1977: Ellipsoid(name='ATS1977', a=6378135, b=6356750.30492159, f_=298.257, f=0.00335281, f2=0.00336409, n=0.00167922, e=0.08181922, e2=0.00669438, e22=0.0067395, e32=0.00335843, A=6367447.14116695, L=10001962.58040571, R1=6371006.7683072, R2=6371005.17780873, R3=6370998.78689182, Rbiaxial=6367451.62986519, Rtriaxial=6372795.55363648)
17@var Ellipsoids.Australia1966: Ellipsoid(name='Australia1966', a=6378160, b=6356774.71919531, f_=298.25, f=0.00335289, f2=0.00336417, n=0.00167926, e=0.08182018, e2=0.00669454, e22=0.00673966, e32=0.00335851, A=6367471.84853228, L=10002001.39064442, R1=6371031.5730651, R2=6371029.9824858, R3=6371023.59124344, Rbiaxial=6367476.337459, Rtriaxial=6372820.40754721)
18@var Ellipsoids.Bessel1841: Ellipsoid(name='Bessel1841', a=6377397.155, b=6356078.962818, f_=299.1528128, f=0.00334277, f2=0.00335398, n=0.00167418, e=0.08169683, e2=0.00667437, e22=0.00671922, e32=0.00334836, A=6366742.52023395, L=10000855.76443237, R1=6370291.09093933, R2=6370289.51012659, R3=6370283.15821523, Rbiaxial=6366746.98155108, Rtriaxial=6372074.29334012)
19@var Ellipsoids.BesselModified: Ellipsoid(name='BesselModified', a=6377492.018, b=6356173.5087127, f_=299.1528128, f=0.00334277, f2=0.00335398, n=0.00167418, e=0.08169683, e2=0.00667437, e22=0.00671922, e32=0.00334836, A=6366837.22474766, L=10001004.52593463, R1=6370385.84823756, R2=6370384.26740131, R3=6370377.91539546, Rbiaxial=6366841.68613115, Rtriaxial=6372169.07716325)
20@var Ellipsoids.CGCS2000: Ellipsoid(name='CGCS2000', a=6378137, b=6356752.31414036, f_=298.2572221, f=0.00335281, f2=0.00336409, n=0.00167922, e=0.08181919, e2=0.00669438, e22=0.0067395, e32=0.00335843, A=6367449.14577105, L=10001965.72923046, R1=6371008.77138012, R2=6371007.18088352, R3=6371000.78997414, Rbiaxial=6367453.63446401, Rtriaxial=6372797.55593326)
21@var Ellipsoids.Clarke1866: Ellipsoid(name='Clarke1866', a=6378206.4, b=6356583.8, f_=294.97869821, f=0.00339008, f2=0.00340161, n=0.00169792, e=0.08227185, e2=0.00676866, e22=0.00681478, e32=0.00339582, A=6367399.68916978, L=10001888.04298286, R1=6370998.86666667, R2=6370997.240633, R3=6370990.70659881, Rbiaxial=6367404.2783313, Rtriaxial=6372807.62791066)
22@var Ellipsoids.Clarke1880: Ellipsoid(name='Clarke1880', a=6378249.145, b=6356514.86954978, f_=293.465, f=0.00340756, f2=0.00341921, n=0.00170669, e=0.0824834, e2=0.00680351, e22=0.00685012, e32=0.00341337, A=6367386.64398051, L=10001867.55164747, R1=6371004.38651659, R2=6371002.74366963, R3=6370996.1419165, Rbiaxial=6367391.2806777, Rtriaxial=6372822.52526083)
23@var Ellipsoids.Clarke1880IGN: Ellipsoid(name='Clarke1880IGN', a=6378249.2, b=6356515, f_=293.46602129, f=0.00340755, f2=0.0034192, n=0.00170668, e=0.08248326, e2=0.00680349, e22=0.00685009, e32=0.00341336, A=6367386.73667336, L=10001867.69724907, R1=6371004.46666667, R2=6371002.82383112, R3=6370996.22212395, Rbiaxial=6367391.37333829, Rtriaxial=6372822.59907505)
24@var Ellipsoids.Clarke1880Mod: Ellipsoid(name='Clarke1880Mod', a=6378249.145, b=6356514.96639549, f_=293.46630766, f=0.00340755, f2=0.0034192, n=0.00170668, e=0.08248322, e2=0.00680348, e22=0.00685009, e32=0.00341335, A=6367386.69236201, L=10001867.62764495, R1=6371004.4187985, R2=6371002.77596616, R3=6370996.17427195, Rbiaxial=6367391.32901784, Rtriaxial=6372822.5494103)
25@var Ellipsoids.CPM1799: Ellipsoid(name='CPM1799', a=6375738.7, b=6356671.92557493, f_=334.39, f=0.00299052, f2=0.00299949, n=0.0014975, e=0.07727934, e2=0.0059721, e22=0.00600798, e32=0.00299499, A=6366208.88184784, L=10000017.52721564, R1=6369383.10852498, R2=6369381.8434158, R3=6369376.76247022, Rbiaxial=6366212.45090321, Rtriaxial=6370977.3559758)
26@var Ellipsoids.Delambre1810: Ellipsoid(name='Delambre1810', a=6376428, b=6355957.92616372, f_=311.5, f=0.00321027, f2=0.00322061, n=0.00160772, e=0.08006397, e2=0.00641024, e22=0.0064516, e32=0.00321543, A=6366197.07684334, L=9999998.98395793, R1=6369604.64205457, R2=6369603.18419749, R3=6369597.32739068, Rbiaxial=6366201.19059818, Rtriaxial=6371316.64722284)
27@var Ellipsoids.Engelis1985: Ellipsoid(name='Engelis1985', a=6378136.05, b=6356751.32272154, f_=298.2566, f=0.00335282, f2=0.0033641, n=0.00167922, e=0.08181928, e2=0.00669439, e22=0.00673951, e32=0.00335844, A=6367448.17507971, L=10001964.20447208, R1=6371007.80757385, R2=6371006.21707085, R3=6370999.82613573, Rbiaxial=6367452.66379074, Rtriaxial=6372796.59560563)
28@var Ellipsoids.Everest1969: Ellipsoid(name='Everest1969', a=6377295.664, b=6356094.667915, f_=300.8017, f=0.00332445, f2=0.00333554, n=0.00166499, e=0.08147298, e2=0.00663785, e22=0.0066822, e32=0.00332998, A=6366699.57839501, L=10000788.3115495, R1=6370228.665305, R2=6370227.10178537, R3=6370220.81951618, Rbiaxial=6366703.99082487, Rtriaxial=6372002.02812501)
29@var Ellipsoids.Everest1975: Ellipsoid(name='Everest1975', a=6377299.151, b=6356098.14512013, f_=300.8017255, f=0.00332445, f2=0.00333554, n=0.00166499, e=0.08147298, e2=0.00663785, e22=0.0066822, e32=0.00332997, A=6366703.06049924, L=10000793.78122603, R1=6370232.14904004, R2=6370230.58551983, R3=6370224.30324826, Rbiaxial=6366707.47293076, Rtriaxial=6372005.51267879)
30@var Ellipsoids.Fisher1968: Ellipsoid(name='Fisher1968', a=6378150, b=6356768.33724438, f_=298.3, f=0.00335233, f2=0.00336361, n=0.00167898, e=0.08181333, e2=0.00669342, e22=0.00673853, e32=0.00335795, A=6367463.65604381, L=10001988.52191361, R1=6371022.77908146, R2=6371021.18903735, R3=6371014.79995035, Rbiaxial=6367468.14345752, Rtriaxial=6372811.30979281)
31@var Ellipsoids.GEM10C: Ellipsoid(name='GEM10C', a=6378137, b=6356752.31424783, f_=298.2572236, f=0.00335281, f2=0.00336409, n=0.00167922, e=0.08181919, e2=0.00669438, e22=0.0067395, e32=0.00335843, A=6367449.14582474, L=10001965.7293148, R1=6371008.77141594, R2=6371007.18091936, R3=6371000.79001005, Rbiaxial=6367453.63451765, Rtriaxial=6372797.55596006)
32@var Ellipsoids.GPES: Ellipsoid(name='GPES', a=6378135, b=6378135, f_=0, f=0, f2=0, n=0, e=0, e2=0, e22=0, e32=0, A=6378135, L=10018751.02980197, R1=6378135, R2=6378135, R3=6378135, Rbiaxial=6378135, Rtriaxial=6378135)
33@var Ellipsoids.GRS67: Ellipsoid(name='GRS67', a=6378160, b=6356774.51609071, f_=298.24716743, f=0.00335292, f2=0.0033642, n=0.00167928, e=0.08182057, e2=0.00669461, e22=0.00673973, e32=0.00335854, A=6367471.74706533, L=10002001.2312605, R1=6371031.50536357, R2=6371029.91475409, R3=6371023.52339015, Rbiaxial=6367476.23607738, Rtriaxial=6372820.3568989)
34@var Ellipsoids.GRS80: Ellipsoid(name='GRS80', a=6378137, b=6356752.31414035, f_=298.2572221, f=0.00335281, f2=0.00336409, n=0.00167922, e=0.08181919, e2=0.00669438, e22=0.0067395, e32=0.00335843, A=6367449.14577104, L=10001965.72923046, R1=6371008.77138012, R2=6371007.18088351, R3=6371000.78997414, Rbiaxial=6367453.634464, Rtriaxial=6372797.55593326)
35@var Ellipsoids.Helmert1906: Ellipsoid(name='Helmert1906', a=6378200, b=6356818.16962789, f_=298.3, f=0.00335233, f2=0.00336361, n=0.00167898, e=0.08181333, e2=0.00669342, e22=0.00673853, e32=0.00335795, A=6367513.57227074, L=10002066.93013953, R1=6371072.7232093, R2=6371071.13315272, R3=6371064.74401563, Rbiaxial=6367518.05971963, Rtriaxial=6372861.26794141)
36@var Ellipsoids.IAU76: Ellipsoid(name='IAU76', a=6378140, b=6356755.28815753, f_=298.257, f=0.00335281, f2=0.00336409, n=0.00167922, e=0.08181922, e2=0.00669438, e22=0.0067395, e32=0.00335843, A=6367452.13278844, L=10001970.42122641, R1=6371011.76271918, R2=6371010.17221946, R3=6371003.78129754, Rbiaxial=6367456.6214902, Rtriaxial=6372800.54945074)
37@var Ellipsoids.IERS1989: Ellipsoid(name='IERS1989', a=6378136, b=6356751.30156878, f_=298.257, f=0.00335281, f2=0.00336409, n=0.00167922, e=0.08181922, e2=0.00669438, e22=0.0067395, e32=0.00335843, A=6367448.13949125, L=10001964.14856985, R1=6371007.76718959, R2=6371006.17669088, R3=6370999.78577297, Rbiaxial=6367452.62819019, Rtriaxial=6372796.55279934)
38@var Ellipsoids.IERS1992TOPEX: Ellipsoid(name='IERS1992TOPEX', a=6378136.3, b=6356751.61659215, f_=298.25722356, f=0.00335281, f2=0.00336409, n=0.00167922, e=0.08181919, e2=0.00669438, e22=0.0067395, e32=0.00335843, A=6367448.44699641, L=10001964.63159782, R1=6371008.07219738, R2=6371006.48170097, R3=6371000.09079236, Rbiaxial=6367452.93568883, Rtriaxial=6372796.85654541)
39@var Ellipsoids.IERS2003: Ellipsoid(name='IERS2003', a=6378136.6, b=6356751.85797165, f_=298.25642, f=0.00335282, f2=0.0033641, n=0.00167922, e=0.0818193, e2=0.0066944, e22=0.00673951, e32=0.00335844, A=6367448.71771058, L=10001965.05683465, R1=6371008.35265722, R2=6371006.76215217, R3=6371000.37120877, Rbiaxial=6367453.20642742, Rtriaxial=6372797.14192686)
40@var Ellipsoids.Intl1924: Ellipsoid(name='Intl1924', a=6378388, b=6356911.94612795, f_=297, f=0.003367, f2=0.00337838, n=0.00168634, e=0.08199189, e2=0.00672267, e22=0.00676817, e32=0.00337267, A=6367654.50005758, L=10002288.29898944, R1=6371229.31537598, R2=6371227.71133444, R3=6371221.26587487, Rbiaxial=6367659.02704315, Rtriaxial=6373025.77129687)
41@var Ellipsoids.Intl1967: Ellipsoid(name='Intl1967', a=6378157.5, b=6356772.2, f_=298.24961539, f=0.0033529, f2=0.00336418, n=0.00167926, e=0.08182023, e2=0.00669455, e22=0.00673967, e32=0.00335852, A=6367469.33894446, L=10001997.44859308, R1=6371029.06666667, R2=6371027.47608389, R3=6371021.08482752, Rbiaxial=6367473.827881, Rtriaxial=6372817.9027631)
42@var Ellipsoids.Krassovski1940: Ellipsoid(name='Krassovski1940', a=6378245, b=6356863.01877305, f_=298.3, f=0.00335233, f2=0.00336361, n=0.00167898, e=0.08181333, e2=0.00669342, e22=0.00673853, e32=0.00335795, A=6367558.49687498, L=10002137.49754285, R1=6371117.67292435, R2=6371116.08285656, R3=6371109.69367439, Rbiaxial=6367562.98435553, Rtriaxial=6372906.23027515)
43@var Ellipsoids.Krassowsky1940: Ellipsoid(name='Krassowsky1940', a=6378245, b=6356863.01877305, f_=298.3, f=0.00335233, f2=0.00336361, n=0.00167898, e=0.08181333, e2=0.00669342, e22=0.00673853, e32=0.00335795, A=6367558.49687498, L=10002137.49754285, R1=6371117.67292435, R2=6371116.08285656, R3=6371109.69367439, Rbiaxial=6367562.98435553, Rtriaxial=6372906.23027515)
44@var Ellipsoids.Maupertuis1738: Ellipsoid(name='Maupertuis1738', a=6397300, b=6363806.28272251, f_=191, f=0.0052356, f2=0.00526316, n=0.00262467, e=0.10219488, e2=0.01044379, e22=0.01055402, e32=0.00524931, A=6380564.13011837, L=10022566.69846922, R1=6386135.42757417, R2=6386131.54144847, R3=6386115.8862823, Rbiaxial=6380575.11882818, Rtriaxial=6388943.03218495)
45@var Ellipsoids.Mercury1960: Ellipsoid(name='Mercury1960', a=6378166, b=6356784.28360711, f_=298.3, f=0.00335233, f2=0.00336361, n=0.00167898, e=0.08181333, e2=0.00669342, e22=0.00673853, e32=0.00335795, A=6367479.62923643, L=10002013.61254591, R1=6371038.76120237, R2=6371037.17115427, R3=6371030.78205124, Rbiaxial=6367484.1166614, Rtriaxial=6372827.29640037)
46@var Ellipsoids.Mercury1968Mod: Ellipsoid(name='Mercury1968Mod', a=6378150, b=6356768.33724438, f_=298.3, f=0.00335233, f2=0.00336361, n=0.00167898, e=0.08181333, e2=0.00669342, e22=0.00673853, e32=0.00335795, A=6367463.65604381, L=10001988.52191361, R1=6371022.77908146, R2=6371021.18903735, R3=6371014.79995035, Rbiaxial=6367468.14345752, Rtriaxial=6372811.30979281)
47@var Ellipsoids.NWL1965: Ellipsoid(name='NWL1965', a=6378145, b=6356759.76948868, f_=298.25, f=0.00335289, f2=0.00336417, n=0.00167926, e=0.08182018, e2=0.00669454, e22=0.00673966, e32=0.00335851, A=6367456.87366841, L=10001977.86818326, R1=6371016.58982956, R2=6371014.999254, R3=6371008.60802667, Rbiaxial=6367461.36258457, Rtriaxial=6372805.42010473)
48@var Ellipsoids.OSU86F: Ellipsoid(name='OSU86F', a=6378136.2, b=6356751.51693008, f_=298.2572236, f=0.00335281, f2=0.00336409, n=0.00167922, e=0.08181919, e2=0.00669438, e22=0.0067395, e32=0.00335843, A=6367448.3471653, L=10001964.47478349, R1=6371007.97231003, R2=6371006.38181364, R3=6370999.99090513, Rbiaxial=6367452.83585765, Rtriaxial=6372796.75662978)
49@var Ellipsoids.OSU91A: Ellipsoid(name='OSU91A', a=6378136.3, b=6356751.6165948, f_=298.2572236, f=0.00335281, f2=0.00336409, n=0.00167922, e=0.08181919, e2=0.00669438, e22=0.0067395, e32=0.00335843, A=6367448.44699773, L=10001964.63159991, R1=6371008.07219827, R2=6371006.48170186, R3=6371000.09079324, Rbiaxial=6367452.93569015, Rtriaxial=6372796.85654607)
50@var Ellipsoids.Plessis1817: Ellipsoid(name='Plessis1817', a=6376523, b=6355862.93325557, f_=308.64, f=0.00324002, f2=0.00325055, n=0.00162264, e=0.08043347, e2=0.00646954, e22=0.00651167, e32=0.00324527, A=6366197.15710739, L=9999999.11003639, R1=6369636.31108519, R2=6369634.82608583, R3=6369628.85999668, Rbiaxial=6366201.34758009, Rtriaxial=6371364.26393357)
51@var Ellipsoids.PZ90: Ellipsoid(name='PZ90', a=6378136, b=6356751.36174571, f_=298.2578393, f=0.0033528, f2=0.00336408, n=0.00167922, e=0.08181911, e2=0.00669437, e22=0.00673948, e32=0.00335842, A=6367448.16955443, L=10001964.19579298, R1=6371007.78724857, R2=6371006.1967588, R3=6370999.80587691, Rbiaxial=6367452.65822809, Rtriaxial=6372796.56780569)
52@var Ellipsoids.SGS85: Ellipsoid(name='SGS85', a=6378136, b=6356751.30156878, f_=298.257, f=0.00335281, f2=0.00336409, n=0.00167922, e=0.08181922, e2=0.00669438, e22=0.0067395, e32=0.00335843, A=6367448.13949125, L=10001964.14856985, R1=6371007.76718959, R2=6371006.17669087, R3=6370999.78577297, Rbiaxial=6367452.62819019, Rtriaxial=6372796.55279934)
53@var Ellipsoids.SoAmerican1969: Ellipsoid(name='SoAmerican1969', a=6378160, b=6356774.71919531, f_=298.25, f=0.00335289, f2=0.00336417, n=0.00167926, e=0.08182018, e2=0.00669454, e22=0.00673966, e32=0.00335851, A=6367471.84853228, L=10002001.39064442, R1=6371031.5730651, R2=6371029.98248581, R3=6371023.59124344, Rbiaxial=6367476.337459, Rtriaxial=6372820.40754721)
54@var Ellipsoids.Sphere: Ellipsoid(name='Sphere', a=6371008.771415, b=6371008.771415, f_=0, f=0, f2=0, n=0, e=0, e2=0, e22=0, e32=0, A=6371008.771415, L=10007557.17611675, R1=6371008.771415, R2=6371008.771415, R3=6371008.771415, Rbiaxial=6371008.771415, Rtriaxial=6371008.771415)
55@var Ellipsoids.SphereAuthalic: Ellipsoid(name='SphereAuthalic', a=6371000, b=6371000, f_=0, f=0, f2=0, n=0, e=0, e2=0, e22=0, e32=0, A=6371000, L=10007543.39801029, R1=6371000, R2=6371000, R3=6371000, Rbiaxial=6371000, Rtriaxial=6371000)
56@var Ellipsoids.SpherePopular: Ellipsoid(name='SpherePopular', a=6378137, b=6378137, f_=0, f=0, f2=0, n=0, e=0, e2=0, e22=0, e32=0, A=6378137, L=10018754.17139462, R1=6378137, R2=6378137, R3=6378137, Rbiaxial=6378137, Rtriaxial=6378137)
57@var Ellipsoids.Struve1860: Ellipsoid(name='Struve1860', a=6378298.3, b=6356657.14266956, f_=294.73, f=0.00339294, f2=0.00340449, n=0.00169935, e=0.0823065, e2=0.00677436, e22=0.00682056, e32=0.00339869, A=6367482.31832549, L=10002017.83655714, R1=6371084.58088985, R2=6371082.95208988, R3=6371076.40691418, Rbiaxial=6367486.91530791, Rtriaxial=6372894.90029454)
58@var Ellipsoids.WGS60: Ellipsoid(name='WGS60', a=6378165, b=6356783.28695944, f_=298.3, f=0.00335233, f2=0.00336361, n=0.00167898, e=0.08181333, e2=0.00669342, e22=0.00673853, e32=0.00335795, A=6367478.63091189, L=10002012.0443814, R1=6371037.76231981, R2=6371036.17227197, R3=6371029.78316994, Rbiaxial=6367483.11833616, Rtriaxial=6372826.29723739)
59@var Ellipsoids.WGS66: Ellipsoid(name='WGS66', a=6378145, b=6356759.76948868, f_=298.25, f=0.00335289, f2=0.00336417, n=0.00167926, e=0.08182018, e2=0.00669454, e22=0.00673966, e32=0.00335851, A=6367456.87366841, L=10001977.86818326, R1=6371016.58982956, R2=6371014.999254, R3=6371008.60802667, Rbiaxial=6367461.36258457, Rtriaxial=6372805.42010473)
60@var Ellipsoids.WGS72: Ellipsoid(name='WGS72', a=6378135, b=6356750.52001609, f_=298.26, f=0.00335278, f2=0.00336406, n=0.0016792, e=0.08181881, e2=0.00669432, e22=0.00673943, e32=0.0033584, A=6367447.24862383, L=10001962.74919858, R1=6371006.84000536, R2=6371005.24953886, R3=6370998.8587507, Rbiaxial=6367451.7372317, Rtriaxial=6372795.60727472)
61@var Ellipsoids.WGS84: Ellipsoid(name='WGS84', a=6378137, b=6356752.31424518, f_=298.25722356, f=0.00335281, f2=0.00336409, n=0.00167922, e=0.08181919, e2=0.00669438, e22=0.0067395, e32=0.00335843, A=6367449.14582341, L=10001965.72931272, R1=6371008.77141506, R2=6371007.18091847, R3=6371000.79000916, Rbiaxial=6367453.63451633, Rtriaxial=6372797.5559594)
62'''
63# make sure int/int division yields float quotient, see .basics
64from __future__ import division as _; del _ # PYCHOK semicolon
66from pygeodesy.basics import copysign0, isint
67from pygeodesy.constants import EPS, EPS0, EPS02, EPS1, INF, NINF, PI4, PI_2, PI_3, R_M, R_MA, R_FM, \
68 _EPSqrt, _EPStol as _TOL, _floatuple as _T, _isfinite, _SQRT2_2, \
69 _0_0s, _0_0, _0_5, _1_0, _1_EPS, _2_0, _4_0, _90_0, \
70 _0_25, _3_0 # PYCHOK used!
71from pygeodesy.errors import _AssertionError, _ValueError, _xkwds_not
72from pygeodesy.fmath import cbrt, cbrt2, fdot, Fhorner, fpowers, Fsum, hypot, hypot_, \
73 hypot1, hypot2, sqrt3
74# from pygeodesy.fsums import Fsum # from .fmath
75from pygeodesy.interns import NN, _a_, _Airy1830_, _AiryModified_, _b_, _Bessel1841_, _beta_, \
76 _Clarke1866_, _Clarke1880IGN_, _DOT_, _f_, _GRS80_, _height_, \
77 _Intl1924_, _incompatible_, _invalid_, _Krassovski1940_, \
78 _Krassowsky1940_, _meridional_, _lat_, _negative_, _not_finite_, \
79 _vs_, _prime_vertical_, _radius_, _Sphere_, _SPACE_, _WGS72_, _WGS84_
80from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS
81from pygeodesy.named import _lazyNamedEnumItem as _lazy, _NamedEnum, \
82 _NamedEnumItem, _NamedTuple, _Pass
83from pygeodesy.namedTuples import Distance2Tuple, Vector3Tuple, Vector4Tuple
84from pygeodesy.props import deprecated_Property_RO, Property_RO, property_doc_, property_RO
85from pygeodesy.streprs import Fmt, fstr, instr, strs, unstr
86from pygeodesy.units import Bearing_, Distance, Float, Float_, Height, Lam_, Lat, Meter, \
87 Meter2, Meter3, Phi, Phi_, Radius, Radius_, Scalar
88from pygeodesy.utily import atand, atan2b, atan2d, degrees90, m2radians, radians2m, sincos2d
90from math import asinh, atan, atanh, cos, degrees, exp, fabs, radians, sin, sinh, sqrt, tan
92__all__ = _ALL_LAZY.ellipsoids
93__version__ = '23.04.11'
95_f_0_0 = Float(f =_0_0) # zero flattening
96_f__0_0 = Float(f_=_0_0) # zero inverse flattening
97# see U{WGS84_f<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Constants.html>}
98_f__WGS84 = Float(f_=_1_0 / (1000000000 / 298257223563)) # 298.25722356299997 vs 298.257223563
101def _aux(lat, inverse, auxLat, clip=90):
102 '''Return a named auxiliary latitude in C{degrees}.
103 '''
104 return Lat(lat, clip=clip, name=_lat_ if inverse else auxLat.__name__)
107def _s2_c2(phi):
108 '''(INTERNAL) Return 2-tuple C{(sin(B{phi})**2, cos(B{phi})**2)}.
109 '''
110 if phi:
111 s2 = sin(phi)**2
112 if s2 > EPS:
113 c2 = _1_0 - s2
114 if c2 > EPS:
115 if c2 < EPS1:
116 return s2, c2
117 else:
118 return _1_0, _0_0 # phi == PI_2
119 return _0_0, _1_0 # phi == 0
122class a_f2Tuple(_NamedTuple):
123 '''2-Tuple C{(a, f)} specifying an ellipsoid by I{equatorial}
124 radius C{a} in C{meter} and scalar I{flattening} C{f}.
126 @see: Class L{Ellipsoid2}.
127 '''
128 _Names_ = (_a_, _f_) # name 'f' not 'f_'
129 _Units_ = (_Pass, _Pass)
131 def __new__(cls, a, f, **name):
132 '''New L{a_f2Tuple} ellipsoid specification.
134 @arg a: Equatorial radius (C{scalar} > 0).
135 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).
136 @kwarg name: Optional name (C{str}).
138 @return: An L{a_f2Tuple}C{(a, f)} instance.
140 @raise UnitError: Invalid B{C{a}} or B{C{f}}.
142 @note: C{abs(B{f}) < EPS} is forced to C{B{f}=0}, I{spherical}.
143 Negative C{B{f}} produces a I{prolate} ellipsoid.
144 '''
145 a = Radius_(a=a) # low=EPS, high=None
146 f = Float_( f=f, low=None, high=EPS1)
147 if fabs(f) < EPS: # force spherical
148 f = _f_0_0
149 return _NamedTuple.__new__(cls, a, f, **name)
151 @Property_RO
152 def b(self):
153 '''Get the I{polar} radius (C{meter}), M{a * (1 - f)}.
154 '''
155 return a_f2b(self.a, self.f) # PYCHOK .a and .f
157 def ellipsoid(self, name=NN):
158 '''Return an L{Ellipsoid} for this 2-tuple C{(a, f)}.
160 @raise NameError: A registered C{ellipsoid} with the
161 same B{C{name}} already exists.
162 '''
163 return Ellipsoid(self.a, f=self.f, name=name or self.name) # PYCHOK .a and .f
165 @Property_RO
166 def f_(self):
167 '''Get the I{inverse} flattening (C{float}), M{1 / f} == M{a / (a - b)}.
168 '''
169 return f2f_(self.f) # PYCHOK .f
172class Circle4Tuple(_NamedTuple):
173 '''4-Tuple C{(radius, height, lat, beta)} of the C{radius} and C{height},
174 both conventionally in C{meter} of a parallel I{circle of latitude} at
175 (geodetic) latitude C{lat} and the I{parametric (or reduced) auxiliary
176 latitude} C{beta}, both in C{degrees90}.
178 The C{height} is the (signed) distance along the z-axis between the
179 parallel and the equator. At near-polar C{lat}s, the C{radius} is C{0},
180 the C{height} is the ellipsoid's (signed) polar radius and C{beta}
181 equals C{lat}.
182 '''
183 _Names_ = (_radius_, _height_, _lat_, _beta_)
184 _Units_ = ( Radius, Height, Lat, Lat)
187class Curvature2Tuple(_NamedTuple):
188 '''2-Tuple C{(meridional, prime_vertical)} of radii of curvature, both in
189 C{meter}, conventionally.
190 '''
191 _Names_ = (_meridional_, _prime_vertical_)
192 _Units_ = ( Meter, Meter)
194 @property_RO
195 def transverse(self):
196 '''Get this I{prime_vertical}, aka I{transverse} radius of curvature.
197 '''
198 return self.prime_vertical
201class Ellipsoid(_NamedEnumItem):
202 '''Ellipsoid with I{equatorial} and I{polar} radii, I{flattening}, I{inverse
203 flattening} and other, often used, I{cached} attributes, supporting
204 I{oblate} and I{prolate} ellipsoidal and I{spherical} earth models.
205 '''
206 _a = 0 # equatorial radius, semi-axis (C{meter})
207 _b = 0 # polar radius, semi-axis (C{meter}): a * (f - 1) / f
208 _f = 0 # (1st) flattening: (a - b) / a
209 _f_ = 0 # inverse flattening: 1 / f = a / (a - b)
211 _KsOrder = 8 # Krüger series order (4, 6 or 8)
213 def __init__(self, a, b=None, f_=None, f=None, name=NN):
214 '''New L{Ellipsoid} from the I{equatorial} radius and either the
215 I{polar} radius or I{inverse flattening} or I{flattening}.
217 @arg a: Equatorial radius, semi-axis (C{meter}).
218 @arg b: Optional polar radius, semi-axis (C{meter}).
219 @arg f_: Inverse flattening: M{a / (a - b)} (C{float} >>> 1.0).
220 @arg f: Flattening: M{(a - b) / a} (C{float}, near zero for
221 spherical).
222 @kwarg name: Optional, unique name (C{str}).
224 @raise NameError: Ellipsoid with the same B{C{name}} already exists.
226 @raise ValueError: Invalid B{C{a}}, B{C{b}}, B{C{f_}} or B{C{f}} or
227 B{C{f_}} and B{C{f}} are incompatible.
229 @note: M{abs(f_) > 1 / EPS} or M{abs(1 / f_) < EPS} is forced
230 to M{1 / f_ = 0}, spherical.
231 '''
232 ff_ = f, f_ # assertion below
233 try:
234 a = Radius_(a=a) # low=EPS
235 if not _isfinite(a):
236 raise ValueError(_SPACE_(_a_, _not_finite_))
238 if b: # not in (_0_0, None)
239 b = Radius_(b=b) # low=EPS
240 f = a_b2f(a, b) if f is None else Float(f=f)
241 f_ = f2f_(f) if f_ is None else Float(f_=f_)
242 elif f is not None:
243 f = Float(f=f)
244 b = a_f2b(a, f)
245 f_ = f2f_(f) if f_ is None else Float(f_=f_)
246 elif f_:
247 f_ = Float(f_=f_)
248 b = a_f_2b(a, f_) # a * (f_ - 1) / f_
249 f = f_2f(f_)
250 else: # only a, spherical
251 f_ = f = 0
252 b = a # superfluous
254 if not f < _1_0: # sanity check, see .ecef.Ecef.__init__
255 raise ValueError(_SPACE_(_f_, _invalid_))
256 if not _isfinite(b):
257 raise ValueError(_SPACE_(_b_, _not_finite_))
259 if fabs(f) < EPS or a == b or not f_: # spherical
260 b = a
261 f = _f_0_0
262 f_ = _f__0_0
264 except (TypeError, ValueError) as x:
265 d = _xkwds_not(None, b=b, f_=f_, f=f)
266 t = instr(self, a=a, name=name, **d)
267 raise _ValueError(t, cause=x)
269 self._a = a
270 self._b = b
271 self._f = f
272 self._f_ = f_
274 self._register(Ellipsoids, name)
276 if f and f_: # see .test/testEllipsoidal.py
277 d = dict(eps=_TOL)
278 if None in ff_: # both f_ and f given
279 d.update(Error=_ValueError, txt=_incompatible_)
280 self._assert(_1_0 / f, f_=f_, **d)
281 self._assert(_1_0 / f_, f =f, **d)
282 self._assert(self.b2_a2, e21=self.e21, eps=EPS)
284 def __eq__(self, other):
285 '''Compare this and an other ellipsoid.
287 @arg other: The other ellipsoid (L{Ellipsoid} or L{Ellipsoid2}).
289 @return: C{True} if equal, C{False} otherwise.
290 '''
291 return self is other or (isinstance(other, Ellipsoid) and
292 self.a == other.a and
293 (self.f == other.f or self.b == other.b))
295 @Property_RO
296 def a(self):
297 '''Get the I{equatorial} radius, semi-axis (C{meter}).
298 '''
299 return self._a
301 equatoradius = a # = Requatorial
303 @Property_RO
304 def a2(self):
305 '''Get the I{equatorial} radius I{squared} (C{meter} I{squared}), M{a**2}.
306 '''
307 return Meter2(a2=self.a**2)
309 @Property_RO
310 def a2_(self):
311 '''Get the inverse of the I{equatorial} radius I{squared} (C{meter} I{squared}), M{1 / a**2}.
312 '''
313 return Float(a2_=_1_0 / self.a2)
315 @Property_RO
316 def a_b(self):
317 '''Get the ratio I{equatorial} over I{polar} radius (C{float}), M{a / b} == M{1 / (1 - f)}.
318 '''
319 return Float(a_b=self.a / self.b if self.f else _1_0)
321 @Property_RO
322 def a2_b(self):
323 '''Get the I{polar} meridional (or polar) radius of curvature (C{meter}), M{a**2 / b}.
325 @see: U{Radii of Curvature
326 <https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>}
327 and U{Moritz, H. (1980), Geodetic Reference System 1980
328 <https://WikiPedia.org/wiki/Earth_radius#cite_note-Moritz-2>}.
330 @note: Symbol C{c} is used by IUGG and IERS for the U{polar radius of curvature
331 <https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>}, see L{c2}
332 and L{R2} or L{Rauthalic}.
333 '''
334 return Radius(a2_b=self.a2 / self.b if self.f else self.a) # = rocPolar
336 @Property_RO
337 def a2_b2(self):
338 '''Get the ratio I{equatorial} over I{polar} radius I{squared} (C{float}),
339 M{(a / b)**2} == M{1 / (1 - e**2)} == M{1 / (1 - e2)} == M{1 / e21}.
340 '''
341 return Float(a2_b2=self.a_b**2 if self.f else _1_0)
343 @Property_RO
344 def a_f(self):
345 '''Get the I{equatorial} radius and I{flattening} (L{a_f2Tuple}), see method C{toEllipsoid2}.
346 '''
347 return a_f2Tuple(self.a, self.f, name=self.name)
349 @Property_RO
350 def A(self):
351 '''Get the UTM I{meridional (or rectifying)} radius (C{meter}).
353 @see: I{Meridian arc unit} U{Q<https://StudyLib.net/doc/7443565/>}.
354 '''
355 A, n = self.a, self.n
356 if n:
357 d = (n + _1_0) * 1048576 / A
358 if d: # use 6 n**2 terms, half-way between the _KsOrder's 4, 6, 8
359 # <https://GeographicLib.SourceForge.io/C++/doc/tmseries30.html>
360 # <https://GeographicLib.SourceForge.io/C++/doc/transversemercator.html> and
361 # <https://www.MyGeodesy.id.AU/documents/Karney-Krueger%20equations.pdf> (3)
362 # A *= fhorner(n**2, 1048576, 262144, 16384, 4096, 1600, 784, 441) / 1048576) / (1 + n)
363 A = Radius(A=Fhorner(n**2, 1048576, 262144, 16384, 4096, 1600, 784, 441).fover(d))
364 return A
366 @Property_RO
367 def _albersCyl(self):
368 '''(INTERNAL) Helper for C{auxAuthalic}.
369 '''
370 return _MODS.albers.AlbersEqualAreaCylindrical(datum=self, name=self.name)
372 @Property_RO
373 def AlphaKs(self):
374 '''Get the I{Krüger} U{Alpha series coefficients<https://GeographicLib.SourceForge.io/C++/doc/tmseries30.html>} (C{KsOrder}C{-tuple}).
375 '''
376 return self._Kseries( # XXX int/int quotients may require from __future__ import division as _; del _ # PYCHOK semicolon
377 # n n**2 n**3 n**4 n**5 n**6 n**7 n**8
378 _T(1/2, -2/3, 5/16, 41/180, -127/288, 7891/37800, 72161/387072, -18975107/50803200),
379 _T(13/48, -3/5, 557/1440, 281/630, -1983433/1935360, 13769/28800, 148003883/174182400), # PYCHOK unaligned
380 _T(61/240, -103/140, 15061/26880, 167603/181440, -67102379/29030400, 79682431/79833600), # PYCHOK unaligned
381 _T(49561/161280, -179/168, 6601661/7257600, 97445/49896, -40176129013/7664025600), # PYCHOK unaligned
382 _T(34729/80640, -3418889/1995840, 14644087/9123840, 2605413599/622702080), # PYCHOK unaligned
383 _T(212378941/319334400, -30705481/10378368, 175214326799/58118860800), # PYCHOK unaligned
384 _T(1522256789/1383782400, -16759934899/3113510400), # PYCHOK unaligned
385 _T(1424729850961/743921418240)) # PYCHOK unaligned
387 @Property_RO
388 def area(self):
389 '''Get the ellipsoid's surface area (C{meter} I{squared}), M{4 * PI * c2}.
391 @see: Properties L{areax}, L{c2} and L{R2} and functions
392 L{ellipsoidalExact.areaOf} and L{ellipsoidalKarney.areaOf}.
393 '''
394 return Meter2(area=self.c2 * PI4)
396 @Property_RO
397 def areax(self):
398 '''Get the ellipsoid's surface area (C{meter} I{squared}), M{4 * PI * c2x}, more
399 accurate for very I{oblate} ellipsoids.
401 @see: Properties L{area}, L{c2x} and L{R2x}, class L{GeodesicExact} and
402 functions L{ellipsoidalExact.areaOf} and L{ellipsoidalKarney.areaOf}.
403 '''
404 return Meter2(areax=self.c2x * PI4)
406 def _assert(self, val, eps=_TOL, f0=_0_0, Error=_AssertionError, txt=NN, **name_value):
407 '''(INTERNAL) Assert a C{name=value} vs C{val}.
408 '''
409 for n, v in name_value.items():
410 if fabs(v - val) > eps: # PYCHOK no cover
411 t = (v, _vs_, val)
412 t = _SPACE_.join(strs(t, prec=12, fmt=Fmt.g))
413 t = Fmt.EQUAL(self._DOT_(n), t)
414 raise Error(t, txt=txt or Fmt.exceeds_eps(eps))
415 return Float(v if self.f else f0, name=n)
416 raise Error(unstr(self._DOT_(self._assert.__name__), val,
417 eps=eps, f0=f0, **name_value))
419 def auxAuthalic(self, lat, inverse=False):
420 '''Compute the I{authalic} auxiliary latitude or the I{inverse} thereof.
422 @arg lat: The geodetic (or I{authalic}) latitude (C{degrees90}).
423 @kwarg inverse: If C{True}, B{C{lat}} is the I{authalic} and
424 return the geodetic latitude (C{bool}).
426 @return: The I{authalic} (or geodetic) latitude in C{degrees90}.
428 @see: U{Inverse-/AuthalicLatitude<https://GeographicLib.SourceForge.io/
429 html/classGeographicLib_1_1Ellipsoid.html>}, U{Authalic latitude
430 <https://WikiPedia.org/wiki/Latitude#Authalic_latitude>}, and
431 U{Snyder<https://Pubs.USGS.gov/pp/1395/report.pdf>}, p 16.
432 '''
433 if self.f:
434 f = self._albersCyl._tanf if inverse else self._albersCyl._txif # PYCHOK attr
435 lat = atand(f(tan(Phi_(lat)))) # PYCHOK attr
436 return _aux(lat, inverse, Ellipsoid.auxAuthalic)
438 def auxConformal(self, lat, inverse=False):
439 '''Compute the I{conformal} auxiliary latitude or the I{inverse} thereof.
441 @arg lat: The geodetic (or I{conformal}) latitude (C{degrees90}).
442 @kwarg inverse: If C{True}, B{C{lat}} is the I{conformal} and
443 return the geodetic latitude (C{bool}).
445 @return: The I{conformal} (or geodetic) latitude in C{degrees90}.
447 @see: U{Inverse-/ConformalLatitude<https://GeographicLib.SourceForge.io/
448 html/classGeographicLib_1_1Ellipsoid.html>}, U{Conformal latitude
449 <https://WikiPedia.org/wiki/Latitude#Conformal_latitude>}, and
450 U{Snyder<https://Pubs.USGS.gov/pp/1395/report.pdf>}, pp 15-16.
451 '''
452 if self.f:
453 f = self.es_tauf if inverse else self.es_taupf # PYCHOK attr
454 lat = atand(f(tan(Phi_(lat)))) # PYCHOK attr
455 return _aux(lat, inverse, Ellipsoid.auxConformal)
457 def auxGeocentric(self, lat, inverse=False):
458 '''Compute the I{geocentric} auxiliary latitude or the I{inverse} thereof.
460 @arg lat: The geodetic (or I{geocentric}) latitude (C{degrees90}).
461 @kwarg inverse: If C{True}, B{C{lat}} is the geocentric and
462 return the I{geocentric} latitude (C{bool}).
464 @return: The I{geocentric} (or geodetic) latitude in C{degrees90}.
466 @see: U{Inverse-/GeocentricLatitude<https://GeographicLib.SourceForge.io/
467 html/classGeographicLib_1_1Ellipsoid.html>}, U{Geocentric latitude
468 <https://WikiPedia.org/wiki/Latitude#Geocentric_latitude>}, and
469 U{Snyder<https://Pubs.USGS.gov/pp/1395/report.pdf>}, pp 17-18.
470 '''
471 if self.f:
472 f = self.a2_b2 if inverse else self.b2_a2
473 lat = atand(f * tan(Phi_(lat)))
474 return _aux(lat, inverse, Ellipsoid.auxGeocentric)
476 def auxIsometric(self, lat, inverse=False):
477 '''Compute the I{isometric} auxiliary latitude or the I{inverse} thereof.
479 @arg lat: The geodetic (or I{isometric}) latitude (C{degrees}).
480 @kwarg inverse: If C{True}, B{C{lat}} is the I{isometric} and
481 return the geodetic latitude (C{bool}).
483 @return: The I{isometric} (or geodetic) latitude in C{degrees}.
485 @note: The I{isometric} latitude for geodetic C{+/-90} is far
486 outside the C{[-90..+90]} range but the inverse
487 thereof is the original geodetic latitude.
489 @see: U{Inverse-/IsometricLatitude<https://GeographicLib.SourceForge.io/
490 html/classGeographicLib_1_1Ellipsoid.html>}, U{Isometric latitude
491 <https://WikiPedia.org/wiki/Latitude#Isometric_latitude>}, and
492 U{Snyder<https://Pubs.USGS.gov/pp/1395/report.pdf>}, pp 15-16.
493 '''
494 if self.f:
495 r = Phi_(lat, clip=0)
496 lat = degrees(atan(self.es_tauf(sinh(r))) if inverse else
497 asinh(self.es_taupf(tan(r))))
498 # clip=0, since auxIsometric(+/-90) is far outside [-90..+90]
499 return _aux(lat, inverse, Ellipsoid.auxIsometric, clip=0)
501 def auxParametric(self, lat, inverse=False):
502 '''Compute the I{parametric} auxiliary latitude or the I{inverse} thereof.
504 @arg lat: The geodetic (or I{parametric}) latitude (C{degrees90}).
505 @kwarg inverse: If C{True}, B{C{lat}} is the I{parametric} and
506 return the geodetic latitude (C{bool}).
508 @return: The I{parametric} (or geodetic) latitude in C{degrees90}.
510 @see: U{Inverse-/ParametricLatitude<https://GeographicLib.SourceForge.io/
511 html/classGeographicLib_1_1Ellipsoid.html>}, U{Parametric latitude
512 <https://WikiPedia.org/wiki/Latitude#Parametric_(or_reduced)_latitude>},
513 and U{Snyder<https://Pubs.USGS.gov/pp/1395/report.pdf>}, p 18.
514 '''
515 if self.f:
516 lat = self._beta(Lat(lat), inverse=inverse)
517 return _aux(lat, inverse, Ellipsoid.auxParametric)
519 auxReduced = auxParametric # synonymous
521 def auxRectifying(self, lat, inverse=False):
522 '''Compute the I{rectifying} auxiliary latitude or the I{inverse} thereof.
524 @arg lat: The geodetic (or I{rectifying}) latitude (C{degrees90}).
525 @kwarg inverse: If C{True}, B{C{lat}} is the I{rectifying} and
526 return the geodetic latitude (C{bool}).
528 @return: The I{rectifying} (or geodetic) latitude in C{degrees90}.
530 @see: U{Inverse-/RectifyingLatitude<https://GeographicLib.SourceForge.io/
531 html/classGeographicLib_1_1Ellipsoid.html>}, U{Rectifying latitude
532 <https://WikiPedia.org/wiki/Latitude#Rectifying_latitude>}, and
533 U{Snyder<https://Pubs.USGS.gov/pp/1395/report.pdf>}, pp 16-17.
534 '''
535 if self.f:
536 lat = Lat(lat)
537 if 0 < fabs(lat) < _90_0:
538 if inverse:
539 e = self._elliptic_e22
540 d = degrees90(e.fEinv(e.cE * lat / _90_0))
541 lat = self.auxParametric(d, inverse=True)
542 else:
543 lat = _90_0 * self.Llat(lat) / self.L
544 return _aux(lat, inverse, Ellipsoid.auxRectifying)
546 @Property_RO
547 def b(self):
548 '''Get the I{polar} radius, semi-axis (C{meter}).
549 '''
550 return self._b
552 polaradius = b # = Rpolar
554 @Property_RO
555 def b_a(self):
556 '''Get the ratio I{polar} over I{equatorial} radius (C{float}), M{b / a == f1 == 1 - f}.
558 @see: Property L{f1}.
559 '''
560 return self._assert(self.b / self.a, b_a=self.f1, f0=_1_0)
562 @Property_RO
563 def b2(self):
564 '''Get the I{polar} radius I{squared} (C{float}), M{b**2}.
565 '''
566 return Meter2(b2=self.b**2)
568 @Property_RO
569 def b2_a(self):
570 '''Get the I{equatorial} meridional radius of curvature (C{meter}), M{b**2 / a}, see C{rocMeridional}C{(0)}.
572 @see: U{Radii of Curvature<https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>}.
573 '''
574 return Radius(b2_a=self.b2 / self.a if self.f else self.b)
576 @Property_RO
577 def b2_a2(self):
578 '''Get the ratio I{polar} over I{equatorial} radius I{squared} (C{float}), M{(b / a)**2}
579 == M{(1 - f)**2} == M{1 - e**2} == C{e21}.
580 '''
581 return Float(b2_a2=self.b_a**2 if self.f else _1_0)
583 def _beta(self, lat, inverse=False):
584 '''(INTERNAL) Get the I{parametric (or reduced) auxiliary latitude} or inverse thereof.
585 '''
586 s, c = sincos2d(lat) # like Karney's tand(lat)
587 s *= self.a_b if inverse else self.b_a
588 return atan2d(s, c) # == atand(s / c) if c else copysign0(_90_0, lat)
590 @Property_RO
591 def BetaKs(self):
592 '''Get the I{Krüger} U{Beta series coefficients<https://GeographicLib.SourceForge.io/C++/doc/tmseries30.html>} (C{KsOrder}C{-tuple}).
593 '''
594 return self._Kseries( # XXX int/int quotients may require from __future__ import division as _; del _ # PYCHOK semicolon
595 # n n**2 n**3 n**4 n**5 n**6 n**7 n**8
596 _T(1/2, -2/3, 37/96, -1/360, -81/512, 96199/604800, -5406467/38707200, 7944359/67737600),
597 _T(1/48, 1/15, -437/1440, 46/105, -1118711/3870720, 51841/1209600, 24749483/348364800), # PYCHOK unaligned
598 _T(17/480, -37/840, -209/4480, 5569/90720, 9261899/58060800, -6457463/17740800), # PYCHOK unaligned
599 _T(4397/161280, -11/504, -830251/7257600, 466511/2494800, 324154477/7664025600), # PYCHOK unaligned
600 _T(4583/161280, -108847/3991680, -8005831/63866880, 22894433/124540416), # PYCHOK unaligned
601 _T(20648693/638668800, -16363163/518918400, -2204645983/12915302400), # PYCHOK unaligne
602 _T(219941297/5535129600, -497323811/12454041600), # PYCHOK unaligned
603 _T(191773887257/3719607091200)) # PYCHOK unaligned
605 @deprecated_Property_RO
606 def c(self): # PYCHOK no cover
607 '''DEPRECATED, use property C{R2} or C{Rauthalic}.'''
608 return self.R2
610 @Property_RO
611 def c2(self):
612 '''Get the I{authalic} earth radius I{squared} (C{meter} I{squared}).
614 @see: Properties L{c2x}, L{area}, L{R2}, L{Rauthalic}, I{Karney's} U{equation (60)
615 <https://Link.Springer.com/article/10.1007%2Fs00190-012-0578-z>} and C++ U{Ellipsoid.Area
616 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Ellipsoid.html>},
617 U{Authalic radius<https://WikiPedia.org/wiki/Earth_radius#Authalic_radius>}, U{Surface area
618 <https://WikiPedia.org/wiki/Ellipsoid>} and U{surface area
619 <https://www.Numericana.com/answer/geometry.htm#oblate>}.
620 '''
621 return self._c2f(False)
623 @Property_RO
624 def c2x(self):
625 '''Get the I{authalic} earth radius I{squared} (C{meter} I{squared}), more accurate for very I{oblate}
626 ellipsoids.
628 @see: Properties L{c2}, L{areax}, L{R2x}, L{Rauthalicx}, class L{GeodesicExact} and I{Karney}'s comments at C++
629 attribute U{GeodesicExact._c2<https://GeographicLib.SourceForge.io/C++/doc/GeodesicExact_8cpp_source.html>}.
630 '''
631 return self._c2f(True)
633 def _c2f(self, c2x):
634 '''(INTERNAL) Helper for C{.c2} and C{.c2x}.
635 '''
636 f, c2 = self.f, self.b2
637 if f:
638 e = self.e
639 if e > EPS0:
640 if f > 0: # .isOblate
641 c2 *= (asinh(sqrt(self.e22abs)) if c2x else atanh(e)) / e
642 elif f < 0: # .isProlate
643 c2 *= atan(e) / e # XXX asin?
644 c2 = Meter2(c2=(self.a2 + c2) * _0_5)
645 return c2
647 def circle4(self, lat):
648 '''Get the equatorial or a parallel I{circle of latitude}.
650 @arg lat: Geodetic latitude (C{degrees90}, C{str}).
652 @return: A L{Circle4Tuple}C{(radius, height, lat, beta)}
653 instance.
655 @raise RangeError: Latitude B{C{lat}} outside valid range and
656 L{pygeodesy.rangerrors} set to C{True}.
658 @raise TypeError: Invalid B{C{lat}}.
660 @raise ValueError: Invalid B{C{lat}}.
662 @see: Definition of U{I{p} and I{z} under B{Parametric (or reduced) latitude}
663 <https://WikiPedia.org/wiki/Latitude>}, I{Karney's} C++ U{CircleRadius and CircleHeight
664 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Ellipsoid.html>}
665 and method C{Rlat}.
666 '''
667 lat = Lat(lat)
668 if lat:
669 b = lat
670 if fabs(lat) < _90_0:
671 if self.f:
672 b = self._beta(lat)
673 z, r = sincos2d(b)
674 r *= self.a
675 z *= self.b
676 else: # near-polar
677 r, z = _0_0, copysign0(self.b, lat)
678 else: # equator
679 r = self.a
680 z = lat = b = _0_0
681 return Circle4Tuple(r, z, lat, b)
683 def degrees2m(self, deg, lat=0):
684 '''Convert an angle to the distance along the equator or
685 along a parallel of (geodetic) latitude.
687 @arg deg: The angle (C{degrees}).
688 @kwarg lat: Parallel latitude (C{degrees90}, C{str}).
690 @return: Distance (C{meter}, same units as the equatorial
691 and polar radii) or C{0} for near-polar B{C{lat}}.
693 @raise RangeError: Latitude B{C{lat}} outside valid range and
694 L{pygeodesy.rangerrors} set to C{True}.
696 @raise ValueError: Invalid B{C{deg}} or B{C{lat}}.
697 '''
698 return self.radians2m(radians(deg), lat=lat)
700 def distance2(self, lat0, lon0, lat1, lon1):
701 '''I{Approximate} the distance and (initial) bearing between
702 two points based on the U{local, flat earth approximation
703 <https://www.EdWilliams.org/avform.htm#flat>} aka U{Hubeny
704 <https://www.OVG.AT/de/vgi/files/pdf/3781/>} formula.
706 I{Suitable only for distances of several hundred Km or Miles
707 and only between points not near-polar}.
709 @arg lat0: From latitude (C{degrees}).
710 @arg lon0: From longitude (C{degrees}).
711 @arg lat1: To latitude (C{degrees}).
712 @arg lon1: To longitude (C{degrees}).
714 @return: A L{Distance2Tuple}C{(distance, initial)} with C{distance}
715 in same units as this ellipsoid's axes.
717 @note: The meridional and prime_vertical radii of curvature are
718 taken and scaled I{at the initial latitude}, see C{roc2}.
720 @see: Function L{pygeodesy.flatLocal}/L{pygeodesy.hubeny}.
721 '''
722 phi0 = Phi_(lat0=lat0)
723 m, n = self.roc2_(phi0, scaled=True)
724 m *= Phi_(lat1=lat1) - phi0
725 n *= Lam_(lon1=lon1) - Lam_(lon0=lon0)
726 return Distance2Tuple(hypot(m, n), atan2b(n, m))
728 @Property_RO
729 def e(self):
730 '''Get the I{unsigned, (1st) eccentricity} (C{float}), M{sqrt(1 - (b / a)**2))}, see C{a_b2e}.
732 @see: Property L{es}.
733 '''
734 return Float(e=sqrt(self.e2abs) if self.e2 else _0_0)
736 @deprecated_Property_RO
737 def e12(self): # see property ._e12
738 '''DEPRECATED, use property C{e21}.'''
739 return self.e21
741# @Property_RO
742# def _e12(self): # see property ._elliptic_e12
743# # (INTERNAL) until e12 above can be replaced with e21.
744# return self.e2 / (_1_0 - self.e2) # see I{Karney}'s Ellipsoid._e12 = e2 / (1 - e2)
746 @Property_RO
747 def e2(self):
748 '''Get the I{signed, (1st) eccentricity squared} (C{float}), M{f * (2 - f)
749 == 1 - (b / a)**2}, see C{a_b2e2}.
750 '''
751 return self._assert(a_b2e2(self.a, self.b), e2=f2e2(self.f))
753 @Property_RO
754 def e2abs(self):
755 '''Get the I{unsigned, (1st) eccentricity squared} (C{float}).
756 '''
757 return fabs(self.e2)
759 @Property_RO
760 def e21(self):
761 '''Get 1 less I{1st eccentricity squared} (C{float}), M{1 - e**2}
762 == M{1 - e2} == M{(1 - f)**2} == M{b**2 / a**2}, see C{b2_a2}.
763 '''
764 return self._assert((_1_0 - self.f)**2, e21=_1_0 - self.e2, f0=_1_0)
766# _e2m = e21 # see I{Karney}'s Ellipsoid._e2m = 1 - _e2
767 _1_e21 = a2_b2 # == M{1 / e21} == M{1 / (1 - e**2)}
769 @Property_RO
770 def e22(self):
771 '''Get the I{signed, 2nd eccentricity squared} (C{float}), M{e2 / (1 - e2)
772 == e2 / (1 - f)**2 == (a / b)**2 - 1}, see C{a_b2e22}.
773 '''
774 return self._assert(a_b2e22(self.a, self.b), e22=f2e22(self.f))
776 @Property_RO
777 def e22abs(self):
778 '''Get the I{unsigned, 2nd eccentricity squared} (C{float}).
779 '''
780 return fabs(self.e22)
782 @Property_RO
783 def e32(self):
784 '''Get the I{signed, 3rd eccentricity squared} (C{float}), M{e2 / (2 - e2)
785 == (a**2 - b**2) / (a**2 + b**2)}, see C{a_b2e32}.
786 '''
787 return self._assert(a_b2e32(self.a, self.b), e32=f2e32(self.f))
789 @Property_RO
790 def e32abs(self):
791 '''Get the I{unsigned, 3rd eccentricity squared} (C{float}).
792 '''
793 return fabs(self.e32)
795 @Property_RO
796 def e4(self):
797 '''Get the I{unsignd, (1st) eccentricity} to 4th power (C{float}), M{e**4 == e2**2}.
798 '''
799 return Float(e4=self.e2**2 if self.e2 else _0_0)
801 eccentricity = e # eccentricity
802# eccentricity2 = e2 # eccentricity squared
803 eccentricity1st2 = e2 # first eccentricity squared
804 eccentricity2nd2 = e22 # second eccentricity squared
805 eccentricity3rd2 = e32 # third eccentricity squared
807 def ecef(self, Ecef=None):
808 '''Return U{ECEF<https://WikiPedia.org/wiki/ECEF>} converter.
810 @kwarg Ecef: ECEF class to use, default L{EcefKarney}.
812 @return: An ECEF converter for this C{ellipsoid}.
814 @raise TypeError: Invalid B{C{Ecef}}.
816 @see: Module L{pygeodesy.ecef}.
817 '''
818 return _MODS.ecef._4Ecef(self, Ecef)
820 @Property_RO
821 def _elliptic_e12(self): # see I{Karney}'s Ellipsoid._e12
822 '''(INTERNAL) Elliptic helper for C{rhumbx.Rhumb}.
823 '''
824 e12 = self.e2 / (_1_0 - self.e2) # NOT DEPRECATED .e12!
825 return _MODS.elliptic.Elliptic(-e12)
827 @Property_RO
828 def _elliptic_e22(self): # aka ._elliptic_ep2
829 '''(INTERNAL) Elliptic helper for C{auxRectifying}, C{L}, C{Llat}.
830 '''
831 return _MODS.elliptic.Elliptic(-self.e22abs) # complex
833 equatoradius = a # Requatorial
835 def e2s(self, s):
836 '''Compute norm M{sqrt(1 - e2 * s**2)}.
838 @arg s: Sine value (C{scalar}).
840 @return: Norm (C{float}).
842 @raise ValueError: Invalid B{C{s}}.
843 '''
844 return sqrt(self.e2s2(s)) if self.e2 else _1_0
846 def e2s2(self, s):
847 '''Compute M{1 - e2 * s**2}.
849 @arg s: Sine value (C{scalar}).
851 @return: Result (C{float}).
853 @raise ValueError: Invalid B{C{s}}.
854 '''
855 r = _1_0
856 if self.e2:
857 try:
858 r -= self.e2 * Scalar(s=s)**2
859 if r < 0:
860 raise ValueError(_negative_)
861 except (TypeError, ValueError) as x:
862 t = self._DOT_(Ellipsoid.e2s2.__name__)
863 raise _ValueError(t, s, cause=x)
864 return r
866 @Property_RO
867 def es(self):
868 '''Get the I{signed (1st) eccentricity} (C{float}).
870 @see: Property L{e}.
871 '''
872 # note, self.e is always non-negative
873 return Float(es=copysign0(self.e, self.f)) # see .ups
875 def es_atanh(self, x):
876 '''Compute M{es * atanh(es * x)} or M{-es * atan(es * x)}
877 for I{oblate} respectively I{prolate} ellipsoids where
878 I{es} is the I{signed} (1st) eccentricity.
880 @raise ValueError: Invalid B{C{x}}.
882 @see: Function U{Math::eatanhe<https://GeographicLib.SourceForge.io/
883 html/classGeographicLib_1_1Math.html>}.
884 '''
885 return self._es_atanh(Scalar(x=x)) if self.f else _0_0
887 def _es_atanh(self, x):
888 '''(INTERNAL) Helper for .es_atanh, ._es_taupf2 and ._exp_es_atanh.
889 '''
890 es = self.es # signOf(es) == signOf(f)
891 return es * (atanh(es * x) if es > 0 else # .isOblate
892 (-atan(es * x) if es < 0 else # .isProlate
893 _0_0)) # .isSpherical
895 @Property_RO
896 def es_c(self):
897 '''Get M{(1 - f) * exp(es_atanh(1))} (C{float}), M{b_a * exp(es_atanh(1))}.
898 '''
899 return Float(es_c=(self._exp_es_atanh_1 * self.b_a) if self.f else _1_0)
901 def es_tauf(self, taup):
902 '''Compute I{Karney}'s U{equations (19), (20) and (21)
903 <https://ArXiv.org/abs/1002.1417>}.
905 @see: I{Karney}'s C++ method U{Math::tauf<https://GeographicLib.
906 SourceForge.io/html/classGeographicLib_1_1Math.html>} and
907 and I{Veness}' JavaScript method U{toLatLon<https://www.
908 Movable-Type.co.UK/scripts/latlong-utm-mgrs.html>}.
909 '''
910 t = Scalar(taup=taup)
911 if self.f: # .isEllipsoidal
912 a = fabs(t)
913 T = t * (self._exp_es_atanh_1 if a > 70 else self._1_e21)
914 if fabs(T * _EPSqrt) < _2_0: # handles +/- INF and NAN
915 s = (a * _TOL) if a > _1_0 else _TOL
916 for T, _, d in self._es_tauf3(t, T): # max 2
917 if fabs(d) < s:
918 break
919 t = Scalar(tauf=T)
920 return t
922 def _es_tauf3(self, taup, T, N=9): # in .utm.Utm._toLLEB
923 '''(INTERNAL) Yield a 3-tuple C{(τi, iteration, delta)} for at most
924 B{C{N}} Newton iterations, converging rapidly except when C{delta}
925 toggles on +/-1.12e-16 or +/-4.47e-16, see C{.utm.Utm._toLLEB}.
926 '''
927 e = self._1_e21
928 _F2_ = Fsum(T).fsum2_ # τ0
929 _tf2 = self._es_taupf2
930 for i in range(1, N + 1):
931 a, h = _tf2(T)
932 d = (taup - a) * (e + T**2) / (hypot1(a) * h)
933 # = (taup - a) / hypot1(a) / ((e + T**2) / h)
934 T, d = _F2_(d) # τi, (τi - τi-1)
935 yield T, i, d
937 def es_taupf(self, tau):
938 '''Compute I{Karney}'s U{equations (7), (8) and (9)
939 <https://ArXiv.org/abs/1002.1417>}.
941 @see: I{Karney}'s C++ method U{Math::taupf<https://GeographicLib.
942 SourceForge.io/html/classGeographicLib_1_1Math.html>}.
943 '''
944 t = Scalar(tau=tau)
945 if self.f: # .isEllipsoidal
946 t, _ = self._es_taupf2(t)
947 t = Scalar(taupf=t)
948 return t
950 def _es_taupf2(self, tau):
951 '''(INTERNAL) Return 2-tuple C{(es_taupf(tau), hypot1(tau))}.
952 '''
953 if _isfinite(tau):
954 h = hypot1(tau)
955 s = sinh(self._es_atanh(tau / h))
956 a = hypot1(s) * tau - h * s
957 else:
958 a, h = tau, INF
959 return a, h
961 @Property_RO
962 def _exp_es_atanh_1(self):
963 '''(INTERNAL) Helper for .es_c and .es_tauf.
964 '''
965 return exp(self._es_atanh(_1_0)) if self.es else _1_0
967 @Property_RO
968 def f(self):
969 '''Get the I{flattening} (C{float}), M{(a - b) / a}, C{0} for spherical, negative for prolate.
970 '''
971 return self._f
973 @Property_RO
974 def f_(self):
975 '''Get the I{inverse flattening} (C{float}), M{1 / f} == M{a / (a - b)}, C{0} for spherical, see C{a_b2f_}.
976 '''
977 return self._f_
979 @Property_RO
980 def f1(self):
981 '''Get the I{1 - flattening} (C{float}), M{f1 == 1 - f == b / a}.
983 @see: Property L{b_a}.
984 '''
985 return Float(f1=_1_0 - self.f)
987 @Property_RO
988 def f2(self):
989 '''Get the I{2nd flattening} (C{float}), M{(a - b) / b == f / (1 - f)}, C{0} for spherical, see C{a_b2f2}.
990 '''
991 return self._assert(self.a_b - _1_0, f2=f2f2(self.f))
993 @Property_RO
994 def geodesic(self):
995 '''Get this ellipsoid's I{wrapped} U{geodesic.Geodesic
996 <https://GeographicLib.SourceForge.io/Python/doc/code.html>}, provided
997 I{Karney}'s U{geographiclib<https://PyPI.org/project/geographiclib>}
998 package is installed.
999 '''
1000 # if not self.isEllipsoidal:
1001 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)
1002 return _MODS.geodesicw.Geodesic(self)
1004 @Property_RO
1005 def geodesicx(self):
1006 '''Get this ellipsoid's I{exact} L{GeodesicExact}.
1007 '''
1008 # if not self.isEllipsoidal:
1009 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)
1010 return _MODS.geodesicx.GeodesicExact(self, name=self.name)
1012 @Property_RO
1013 def geodsolve(self):
1014 '''Get this ellipsoid's L{GeodesicSolve}, the I{wrapper} around utility
1015 U{GeodSolve<https://GeographicLib.SourceForge.io/C++/doc/GeodSolve.1.html>},
1016 provided the path to the C{GeodSolve} executable is specified with env
1017 variable C{PYGEODESY_GEODSOLVE}.
1018 '''
1019 # if not self.isEllipsoidal:
1020 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)
1021 return _MODS.geodsolve.GeodesicSolve(self, name=self.name)
1023 def height4(self, xyz, normal=True):
1024 '''Compute the projection on and the height of a cartesian above or below
1025 this ellipsoid's surface.
1027 @arg xyz: The cartesian (C{Cartesian}, L{Ecef9Tuple}, L{Vector3d},
1028 L{Vector3Tuple} or L{Vector4Tuple}).
1029 @kwarg normal: If C{True} the projection is perpendicular to (the nearest
1030 point on) this ellipsoid's surface, otherwise the C{radial}
1031 line to this ellipsoid's center (C{bool}).
1033 @return: L{Vector4Tuple}C{(x, y, z, h)} with the cartesian coordinates
1034 C{x}, C{y} and C{z} of the projection on or the intersection
1035 with and with the height C{h} above the ellipsoid's surface
1036 in C{meter}, conventionally.
1038 @raise ValueError: Null B{C{xyz}}.
1040 @raise TypeError: Non-cartesian B{C{xyz}}.
1042 @see: U{Distance to<https://StackOverflow.com/questions/22959698/distance-from-given-point-to-given-ellipse>}
1043 and U{intersection with<https://MathWorld.wolfram.com/Ellipse-LineIntersection.html>} an ellipse and
1044 method L{Triaxial.height4}.
1045 '''
1046 v = _MODS.vector3d._otherV3d(xyz=xyz)
1047 r = v.length
1049 a, b, i = self.a, self.b, None
1050 if r < EPS0: # EPS
1051 v = v.times(_0_0)
1052 h = -a
1054 elif self.isSpherical:
1055 v = v.times(a / r)
1056 h = r - a
1058 elif normal: # perpendicular to ellipsoid
1059 x, y = hypot(v.x, v.y), fabs(v.z)
1060 if x < EPS0: # PYCHOK no cover
1061 z = copysign0(b, v.z)
1062 v = Vector3Tuple(v.x, v.y, z)
1063 h = y - b # polar
1064 elif y < EPS0: # PYCHOK no cover
1065 t = a / r
1066 v = v.times_(t, t, 0) # force z=0.0
1067 h = x - a # equatorial
1068 else: # normal in 1st quadrant
1069 x, y, i = _normalTo3(x, y, self)
1070 t, v = v, v.times_(x, x, y)
1071 h = t.minus(v).length
1073 else: # radial to ellipsoid's center
1074 h = hypot_(a * v.z, b * v.x, b * v.y)
1075 t = (a * b / h) if h > EPS0 else _0_0 # EPS
1076 v = v.times(t)
1077 h = r * (_1_0 - t)
1079 return Vector4Tuple(v.x, v.y, v.z, h, iteration=i,
1080 name=self.height4.__name__)
1082 def _hubeny_2(self, phi2, phi1, lam21):
1083 '''(INTERNAL) like function C{pygeodesy.flatLocal_}/C{pygeodesy.hubeny_}
1084 but returning the I{angular} distance in C{radians squared}.
1085 '''
1086 m, n = self.roc2_((phi2 + phi1) * _0_5, scaled=True)
1087 return hypot2(m * (phi2 - phi1), n * lam21) * self.a2_
1089 @Property_RO
1090 def isEllipsoidal(self):
1091 '''Is this model I{ellipsoidal} (C{bool})?
1092 '''
1093 return self.f != 0
1095 @Property_RO
1096 def isOblate(self):
1097 '''Is this ellipsoid I{oblate} (C{bool})? I{Prolate} or
1098 spherical otherwise.
1099 '''
1100 return self.f > 0
1102 @Property_RO
1103 def isProlate(self):
1104 '''Is this ellipsoid I{prolate} (C{bool})? I{Oblate} or
1105 spherical otherwise.
1106 '''
1107 return self.f < 0
1109 @Property_RO
1110 def isSpherical(self):
1111 '''Is this ellipsoid I{spherical} (C{bool})?
1112 '''
1113 return self.f == 0
1115 def _Kseries(self, *AB8Ks):
1116 '''(INTERNAL) Compute the 4-, 6- or 8-th order I{Krüger} Alpha
1117 or Beta series coefficients per I{Karney}'s U{equations (35)
1118 and (36)<https://ArXiv.org/pdf/1002.1417v3.pdf>}.
1120 @arg AB8Ks: 8-Tuple of 8-th order I{Krüger} Alpha or Beta series
1121 coefficient tuples.
1123 @return: I{Krüger} series coefficients (L{KsOrder}C{-tuple}).
1125 @see: I{Karney}'s 30-th order U{TMseries30
1126 <https://GeographicLib.SourceForge.io/C++/doc/tmseries30.html>}.
1127 '''
1128 k = self.KsOrder
1129 if self.n:
1130 ns = fpowers(self.n, k)
1131 ks = tuple(fdot(AB8Ks[i][:k-i], *ns[i:]) for i in range(k))
1132 else:
1133 ks = _0_0s(k)
1134 return ks
1136 @property_doc_(''' the I{Krüger} series' order (C{int}), see properties C{AlphaKs}, C{BetaKs}.''')
1137 def KsOrder(self):
1138 '''Get the I{Krüger} series' order (C{int} 4, 6 or 8).
1139 '''
1140 return self._KsOrder
1142 @KsOrder.setter # PYCHOK setter!
1143 def KsOrder(self, order):
1144 '''Set the I{Krüger} series' order (C{int} 4, 6 or 8).
1146 @raise ValueError: Invalid B{C{order}}.
1147 '''
1148 if not (isint(order) and order in (4, 6, 8)):
1149 raise _ValueError(order=order)
1150 if order != self._KsOrder:
1151 Ellipsoid.AlphaKs._update(self)
1152 Ellipsoid.BetaKs._update(self)
1153 self._KsOrder = order
1155 @Property_RO
1156 def L(self):
1157 '''Get the I{quarter meridian} C{L}, aka the C{polar distance}
1158 along a meridian between the equator and a pole (C{meter}),
1159 M{b * Elliptic(-e2 / (1 - e2)).E} or M{b * PI / 2}.
1160 '''
1161 r = self._elliptic_e22.cE if self.f else PI_2
1162 return Distance(L=self.b * r)
1164 @Property_RO
1165 def _L_90(self):
1166 '''Get the I{quarter meridian} per degree (C{meter}), M{self.L / 90}.
1167 '''
1168 return self.L / _90_0
1170 def Llat(self, lat):
1171 '''Return the I{meridional length}, the distance along a meridian
1172 between the equator and a (geodetic) latitude, see C{L}.
1174 @arg lat: Geodetic latitude (C{degrees90}).
1176 @return: The meridional length at B{C{lat}}, negative on southern
1177 hemisphere (C{meter}).
1178 '''
1179 r = self._elliptic_e22.fEd(self.auxParametric(lat)) if self.f else Phi_(lat)
1180 return Distance(Llat=self.b * r)
1182 Lmeridian = Llat # meridional distance
1184 @deprecated_Property_RO
1185 def majoradius(self): # PYCHOK no cover
1186 '''DEPRECATED, use property C{a} or C{Requatorial}.'''
1187 return self.a
1189 def m2degrees(self, distance, lat=0):
1190 '''Convert a distance to an angle along the equator or
1191 along a parallel of (geodetic) latitude.
1193 @arg distance: Distance (C{meter}).
1194 @kwarg lat: Parallel latitude (C{degrees90}, C{str}).
1196 @return: Angle (C{degrees}) or C{INF} for near-polar B{C{lat}}.
1198 @raise RangeError: Latitude B{C{lat}} outside valid range and
1199 L{pygeodesy.rangerrors} set to C{True}.
1201 @raise ValueError: Invalid B{C{distance}} or B{C{lat}}.
1202 '''
1203 return degrees(self.m2radians(distance, lat=lat))
1205 def m2radians(self, distance, lat=0):
1206 '''Convert a distance to an angle along the equator or
1207 along a parallel of (geodetic) latitude.
1209 @arg distance: Distance (C{meter}).
1210 @kwarg lat: Parallel latitude (C{degrees90}, C{str}).
1212 @return: Angle (C{radians}) or C{INF} for near-polar B{C{lat}}.
1214 @raise RangeError: Latitude B{C{lat}} outside valid range and
1215 L{pygeodesy.rangerrors} set to C{True}.
1217 @raise ValueError: Invalid B{C{distance}} or B{C{lat}}.
1218 '''
1219 r = self.circle4(lat).radius if lat else self.a
1220 return m2radians(distance, radius=r, lat=0)
1222 @deprecated_Property_RO
1223 def minoradius(self): # PYCHOK no cover
1224 '''DEPRECATED, use property C{b} or C{Rpolar}.'''
1225 return self.b
1227 @Property_RO
1228 def n(self):
1229 '''Get the I{3rd flattening} (C{float}), M{f / (2 - f) == (a - b) / (a + b)}, see C{a_b2n}.
1230 '''
1231 return self._assert(a_b2n(self.a, self.b), n=f2n(self.f))
1233 flattening = f
1234 flattening1st = f
1235 flattening2nd = f2
1236 flattening3rd = n
1238 polaradius = b # Rpolar
1240# @Property_RO
1241# def Q(self):
1242# '''Get the I{meridian arc unit} C{Q}, the mean, meridional length I{per radian} C({float}).
1243#
1244# @note: C{Q * PI / 2} ≈ C{L}, the I{quarter meridian}.
1245#
1246# @see: Property C{A} and U{Engsager, K., Poder, K.<https://StudyLib.net/doc/7443565/
1247# a-highly-accurate-world-wide-algorithm-for-the-transverse...>}.
1248# '''
1249# n = self.n
1250# d = (n + _1_0) / self.a
1251# return Float(Q=Fhorner(n**2, _1_0, _0_25, _1_16th, _0_25).fover(d) if d else self.b)
1253# # Moritz, H. <https://Geodesy.Geology.Ohio-State.EDU/course/refpapers/00740128.pdf>
1254# # Q = (1 - 3/4 * e'2 + 45/64 * e'4 - 175/256 * e'6 + 11025/16384 * e'8) * rocPolar
1255# # = (4 + e'2 * (-3 + e'2 * (45/16 + e'2 * (-175/64 + e'2 * 11025/4096)))) * rocPolar / 4
1256# return Fhorner(self.e22, 4, -3, 45 / 16, -175 / 64, 11025 / 4096).fover(4 / self.rocPolar)
1258 @deprecated_Property_RO
1259 def quarteradius(self): # PYCHOK no cover
1260 '''DEPRECATED, use property C{L} or method C{Llat}.'''
1261 return self.L
1263 @Property_RO
1264 def R1(self):
1265 '''Get the I{mean} earth radius per I{IUGG} (C{meter}), M{(2 * a + b) / 3 == a * (1 - f / 3)}.
1267 @see: U{Earth radius<https://WikiPedia.org/wiki/Earth_radius>}
1268 and method C{Rgeometric}.
1269 '''
1270 r = Fsum(self.a, self.a, self.b).fover(_3_0) if self.f else self.a
1271 return Radius(R1=r)
1273 Rmean = R1
1275 @Property_RO
1276 def R2(self):
1277 '''Get the I{authalic} earth radius (C{meter}), M{sqrt(c2)}.
1279 @see: C{R2x}, C{c2}, C{area} and U{Earth radius
1280 <https://WikiPedia.org/wiki/Earth_radius>}.
1281 '''
1282 return Radius(R2=sqrt(self.c2) if self.f else self.a)
1284 Rauthalic = R2
1286# @Property_RO
1287# def R2(self):
1288# # Moritz, H. <https://Geodesy.Geology.Ohio-State.EDU/course/refpapers/00740128.pdf>
1289# # R2 = (1 - 2/3 * e'2 + 26/45 * e'4 - 100/189 * e'6 + 7034/14175 * e'8) * rocPolar
1290# # = (3 + e'2 * (-2 + e'2 * (26/15 + e'2 * (-100/63 + e'2 * 7034/4725)))) * rocPolar / 3
1291# return Fhorner(self.e22, 3, -2, 26 / 15, -100 / 63, 7034 / 4725).fover(3 / self.rocPolar)
1293 @Property_RO
1294 def R2x(self):
1295 '''Get the I{authalic} earth radius (C{meter}), M{sqrt(c2x)}.
1297 @see: C{R2}, C{c2x} and C{areax}.
1298 '''
1299 return Radius(R2x=sqrt(self.c2x) if self.f else self.a)
1301 Rauthalicx = R2x
1303 @Property_RO
1304 def R3(self):
1305 '''Get the I{volumetric} earth radius (C{meter}), M{(a * a * b)**(1/3)}.
1307 @see: U{Earth radius<https://WikiPedia.org/wiki/Earth_radius>} and C{volume}.
1308 '''
1309 r = (cbrt(self.b_a) * self.a) if self.f else self.a
1310 return Radius(R3=r)
1312 Rvolumetric = R3
1314 def radians2m(self, rad, lat=0):
1315 '''Convert an angle to the distance along the equator or
1316 along a parallel of (geodetic) latitude.
1318 @arg rad: The angle (C{radians}).
1319 @kwarg lat: Parallel latitude (C{degrees90}, C{str}).
1321 @return: Distance (C{meter}, same units as the equatorial
1322 and polar radii) or C{0} for near-polar B{C{lat}}.
1324 @raise RangeError: Latitude B{C{lat}} outside valid range and
1325 L{pygeodesy.rangerrors} set to C{True}.
1327 @raise ValueError: Invalid B{C{rad}} or B{C{lat}}.
1328 '''
1329 r = self.circle4(lat).radius if lat else self.a
1330 return radians2m(rad, radius=r, lat=0)
1332 @Property_RO
1333 def Rbiaxial(self):
1334 '''Get the I{biaxial, quadratic} mean earth radius (C{meter}), M{sqrt((a**2 + b**2) / 2)}.
1336 @see: C{Rtriaxial}
1337 '''
1338 a, b = self.a, self.b
1339 q = (sqrt((_1_0 + self.b2_a2) * _0_5) * a) if a > b else (
1340 (sqrt((_1_0 + self.a2_b2) * _0_5) * b) if a < b else a)
1341 return Radius(Rbiaxial=q)
1343 Requatorial = a # for consistent naming
1345 def Rgeocentric(self, lat):
1346 '''Compute the I{geocentric} earth radius of (geodetic) latitude.
1348 @arg lat: Latitude (C{degrees90}).
1350 @return: Geocentric earth radius (C{meter}).
1352 @raise ValueError: Invalid B{C{lat}}.
1354 @see: U{Geocentric Radius
1355 <https://WikiPedia.org/wiki/Earth_radius#Geocentric_radius>}
1356 '''
1357 r, a = self.a, Phi_(lat)
1358 if a and self.f:
1359 if fabs(a) < PI_2:
1360 s2, c2 = _s2_c2(a)
1361 b2_a2_s2 = self.b2_a2 * s2
1362 # R == sqrt((a2**2 * c2 + b2**2 * s2) / (a2 * c2 + b2 * s2))
1363 # == sqrt(a2**2 * (c2 + (b2 / a2)**2 * s2) / (a2 * (c2 + b2 / a2 * s2)))
1364 # == sqrt(a2 * (c2 + (b2 / a2)**2 * s2) / (c2 + (b2 / a2) * s2))
1365 # == a * sqrt((c2 + b2_a2 * b2_a2 * s2) / (c2 + b2_a2 * s2))
1366 # == a * sqrt((c2 + b2_a2 * b2_a2_s2) / (c2 + b2_a2_s2))
1367 r *= sqrt((c2 + b2_a2_s2 * self.b2_a2) / (c2 + b2_a2_s2))
1368 else:
1369 r = self.b
1370 return Radius(Rgeocentric=r)
1372 @Property_RO
1373 def Rgeometric(self):
1374 '''Get the I{geometric} mean earth radius (C{meter}), M{sqrt(a * b)}.
1376 @see: C{R1}.
1377 '''
1378 g = sqrt(self.a * self.b) if self.f else self.a
1379 return Radius(Rgeometric=g)
1381 @Property_RO
1382 def rhumbsolve(self):
1383 '''Get this ellipsoid's L{RhumbSolve}, the I{wrapper} around utility
1384 U{RhumbSolve<https://GeographicLib.SourceForge.io/C++/doc/GeodSolve.1.html>},
1385 provided the path to the C{RhumbSolve} executable is specified with env
1386 variable C{PYGEODESY_RHUMBSOLVE}.
1387 '''
1388 # if not self.isEllipsoidal:
1389 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)
1390 return _MODS.rhumbsolve.RhumbSolve(self, name=self.name)
1392 @Property_RO
1393 def rhumbx(self):
1394 '''Get this ellipsoid's L{Rhumb}.
1395 '''
1396 # if not self.isEllipsoidal:
1397 # raise _IsnotError(_ellipsoidal_, ellipsoid=self)
1398 return _MODS.rhumbx.Rhumb(self, name=self.name)
1400 def Rlat(self, lat):
1401 '''I{Approximate} the earth radius of (geodetic) latitude.
1403 @arg lat: Latitude (C{degrees90}).
1405 @return: Approximate earth radius (C{meter}).
1407 @raise RangeError: Latitude B{C{lat}} outside valid range and
1408 L{pygeodesy.rangerrors} set to C{True}.
1410 @raise TypeError: Invalid B{C{lat}}.
1412 @raise ValueError: Invalid B{C{lat}}.
1414 @note: C{Rlat(B{90})} equals C{Rpolar}.
1416 @see: Method C{circle4}.
1417 '''
1418 # r = a - (a - b) * |lat| / 90
1419 r = self.a
1420 if self.f and lat: # .isEllipsoidal
1421 r -= (r - self.b) * fabs(Lat(lat)) / _90_0
1422 r = Radius(Rlat=r)
1423 return r
1425 Rpolar = b # for consistent naming
1427 @deprecated_Property_RO
1428 def Rquadratic(self): # PYCHOK no cover
1429 '''DEPRECATED, use property C{Rbiaxial} or C{Rtriaxial}.'''
1430 return self.Rbiaxial
1432 @deprecated_Property_RO
1433 def Rr(self): # PYCHOK no cover
1434 '''DEPRECATED, use property C{Rrectifying}.'''
1435 return self.Rrectifying
1437 @Property_RO
1438 def Rrectifying(self):
1439 '''Get the I{rectifying} earth radius (C{meter}), M{((a**(3/2) + b**(3/2)) / 2)**(2/3)}.
1441 @see: U{Earth radius<https://WikiPedia.org/wiki/Earth_radius>}.
1442 '''
1443 r = (cbrt2((_1_0 + sqrt3(self.b_a)) * _0_5) * self.a) if self.f else self.a
1444 return Radius(Rrectifying=r)
1446 def roc1_(self, sa, ca=None):
1447 '''Compute the I{prime-vertical}, I{normal} radius of curvature
1448 of (geodetic) latitude, I{unscaled}.
1450 @arg sa: Sine of the latitude (C{float}, [-1.0..+1.0]).
1451 @kwarg ca: Optional cosine of the latitude (C{float}, [-1.0..+1.0])
1452 to use an alternate formula.
1454 @return: The prime-vertical radius of curvature (C{float}).
1456 @note: The delta between both formulae with C{Ellipsoids.WGS84}
1457 is less than 2 nanometer over the entire latitude range.
1459 @see: Method L{roc2_} and class L{EcefYou}.
1460 '''
1461 if not self.f: # .isSpherical
1462 n = self.a
1463 elif ca is None:
1464 r = self.e2s2(sa) # see .roc2_ and _EcefBase._forward
1465 n = (self.a / sqrt(r)) if r > EPS02 else _0_0
1466 elif ca: # derived from EcefYou.forward
1467 h = hypot(ca, self.b_a * sa) if sa else fabs(ca)
1468 n = self.a / h
1469 elif sa:
1470 n = self.a2_b / fabs(sa)
1471 else:
1472 n = self.a
1473 return n
1475 def roc2(self, lat, scaled=False):
1476 '''Compute the I{meridional} and I{prime-vertical}, I{normal}
1477 radii of curvature of (geodetic) latitude.
1479 @arg lat: Latitude (C{degrees90}).
1480 @kwarg scaled: Scale prime_vertical by C{cos(radians(B{lat}))} (C{bool}).
1482 @return: An L{Curvature2Tuple}C{(meridional, prime_vertical)} with
1483 the radii of curvature.
1485 @raise ValueError: Invalid B{C{lat}}.
1487 @see: Methods L{roc2_} and L{roc1_}, U{Local, flat earth approximation
1488 <https://www.EdWilliams.org/avform.htm#flat>} and meridional and
1489 prime vertical U{Radii of Curvature<https://WikiPedia.org/wiki/
1490 Earth_radius#Radii_of_curvature>}.
1491 '''
1492 return self.roc2_(Phi_(lat), scaled=scaled)
1494 def roc2_(self, phi, scaled=False):
1495 '''Compute the I{meridional} and I{prime-vertical}, I{normal} radii of
1496 curvature of (geodetic) latitude.
1498 @arg phi: Latitude (C{radians}).
1499 @kwarg scaled: Scale prime_vertical by C{cos(B{phi})} (C{bool}).
1501 @return: An L{Curvature2Tuple}C{(meridional, prime_vertical)} with the
1502 radii of curvature.
1504 @raise ValueError: Invalid B{C{phi}}.
1506 @see: Methods L{roc2} and L{roc1_}, property L{rocEquatorial2}, U{Local,
1507 flat earth approximation<https://www.EdWilliams.org/avform.htm#flat>}
1508 and the meridional and prime vertical U{Radii of Curvature
1509 <https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>}.
1510 '''
1511 a = fabs(Phi(phi))
1512 if self.f:
1513 r = self.e2s2(sin(a))
1514 if r > EPS02:
1515 n = self.a / sqrt(r)
1516 m = n * self.e21 / r # PYCHOK attr
1517 else:
1518 m = n = _0_0 # PYCHOK attr
1519 else:
1520 m = n = self.a
1521 if scaled and a:
1522 n *= cos(a) if a < PI_2 else _0_0
1523 return Curvature2Tuple(Radius(rocMeridional=m),
1524 Radius(rocPrimeVertical=n))
1526 def rocBearing(self, lat, bearing):
1527 '''Compute the I{directional} radius of curvature of (geodetic)
1528 latitude and compass direction.
1530 @arg lat: Latitude (C{degrees90}).
1531 @arg bearing: Direction (compass C{degrees360}).
1533 @return: Directional radius of curvature (C{meter}).
1535 @raise RangeError: Latitude B{C{lat}} outside valid range and
1536 L{pygeodesy.rangerrors} set to C{True}.
1538 @raise ValueError: Invalid B{C{lat}} or B{C{bearing}}.
1540 @see: U{Radii of Curvature<https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>}
1541 '''
1542 if self.f:
1543 s2, c2 = _s2_c2(Bearing_(bearing))
1544 m, n = self.roc2_(Phi_(lat))
1545 if n < m: # == n / (c2 * n / m + s2)
1546 c2 *= n / m
1547 elif m < n: # == m / (c2 + s2 * m / n)
1548 s2 *= m / n
1549 n = m
1550 b = n / (c2 + s2) # == 1 / (c2 / m + s2 / n)
1551 else:
1552 b = self.b # == self.a
1553 return Radius(rocBearing=b)
1555 @Property_RO
1556 def rocEquatorial2(self):
1557 '''Get the I{meridional} and I{prime-vertical}, I{normal} radii of curvature
1558 at the equator as L{Curvature2Tuple}C{(meridional, prime_vertical)}.
1560 @see: Methods L{rocMeridional} and L{rocPrimeVertical}, properties L{b2_a},
1561 L{a2_b}, C{rocPolar} and polar and equatorial U{Radii of Curvature
1562 <https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>}.
1563 '''
1564 return Curvature2Tuple(Radius(rocMeridional0=self.b2_a if self.f else self.a),
1565 Radius(rocPrimeVertical0=self.a))
1567 def rocGauss(self, lat):
1568 '''Compute the I{Gaussian} radius of curvature of (geodetic) latitude.
1570 @arg lat: Latitude (C{degrees90}).
1572 @return: Gaussian radius of curvature (C{meter}).
1574 @raise ValueError: Invalid B{C{lat}}.
1576 @see: Non-directional U{Radii of Curvature<https://WikiPedia.org/wiki/
1577 Earth_radius#Radii_of_curvature>}
1578 '''
1579 # using ...
1580 # m, n = self.roc2_(Phi_(lat))
1581 # return sqrt(m * n)
1582 # ... requires 1 or 2 sqrt
1583 g = self.b
1584 if self.f:
1585 s2, c2 = _s2_c2(Phi_(lat))
1586 g = g / (c2 + self.b2_a2 * s2)
1587 return Radius(rocGauss=g)
1589 def rocMean(self, lat):
1590 '''Compute the I{mean} radius of curvature of (geodetic) latitude.
1592 @arg lat: Latitude (C{degrees90}).
1594 @return: Mean radius of curvature (C{meter}).
1596 @raise ValueError: Invalid B{C{lat}}.
1598 @see: Non-directional U{Radii of Curvature<https://WikiPedia.org/wiki/
1599 Earth_radius#Radii_of_curvature>}
1600 '''
1601 if self.f:
1602 m, n = self.roc2_(Phi_(lat))
1603 m *= n * _2_0 / (m + n) # == 2 / (1 / m + 1 / n)
1604 else:
1605 m = self.a
1606 return Radius(rocMean=m)
1608 def rocMeridional(self, lat):
1609 '''Compute the I{meridional} radius of curvature of (geodetic) latitude.
1611 @arg lat: Latitude (C{degrees90}).
1613 @return: Meridional radius of curvature (C{meter}).
1615 @raise ValueError: Invalid B{C{lat}}.
1617 @see: Methods L{roc2} and L{roc2_}, U{Local, flat earth approximation
1618 <https://www.EdWilliams.org/avform.htm#flat>} and U{Radii of
1619 Curvature<https://WikiPedia.org/wiki/Earth_radius#Radii_of_curvature>}.
1620 '''
1621 return self.roc2_(Phi_(lat)).meridional if lat else \
1622 self.rocEquatorial2.meridional
1624 rocPolar = a2_b # synonymous
1626 def rocPrimeVertical(self, lat):
1627 '''Compute the I{prime-vertical}, I{normal} radius of curvature of
1628 (geodetic) latitude, aka the I{transverse} radius of curvature.
1630 @arg lat: Latitude (C{degrees90}).
1632 @return: Prime-vertical radius of curvature (C{meter}).
1634 @raise ValueError: Invalid B{C{lat}}.
1636 @see: Methods L{roc2}, L{roc2_} and L{roc1_}, U{Local, flat earth
1637 approximation<https://www.EdWilliams.org/avform.htm#flat>} and
1638 U{Radii of Curvature<https://WikiPedia.org/wiki/
1639 Earth_radius#Radii_of_curvature>}.
1640 '''
1641 return self.roc2_(Phi_(lat)).prime_vertical if lat else \
1642 self.rocEquatorial2.prime_vertical
1644 rocTransverse = rocPrimeVertical # synonymous
1646 @deprecated_Property_RO
1647 def Rs(self): # PYCHOK no cover
1648 '''DEPRECATED, use property C{Rgeometric}.'''
1649 return self.Rgeometric
1651 @Property_RO
1652 def Rtriaxial(self):
1653 '''Get the I{triaxial, quadratic} mean earth radius (C{meter}), M{sqrt((3 * a**2 + b**2) / 4)}.
1655 @see: C{Rbiaxial}
1656 '''
1657 a, b = self.a, self.b
1658 q = (sqrt((_3_0 + self.b2_a2) * _0_25) * a) if a > b else (
1659 (sqrt((_3_0 * self.a2_b2 + _1_0) * _0_25) * b) if a < b else a)
1660 return Radius(Rtriaxial=q)
1662 def toEllipsoid2(self, name=NN):
1663 '''Get a copy of this ellipsoid as an L{Ellipsoid2}.
1665 @kwarg name: Optional, unique name (C{str}).
1667 @see: Property C{a_f}.
1668 '''
1669 return Ellipsoid2(self, None, name=name)
1671 def toStr(self, prec=8, name=NN, **unused): # PYCHOK expected
1672 '''Return this ellipsoid as a text string.
1674 @kwarg prec: Number of decimal digits, unstripped (C{int}).
1675 @kwarg name: Override name (C{str}) or C{None} to exclude
1676 this ellipsoid's name.
1678 @return: This C{Ellipsoid}'s attributes (C{str}).
1679 '''
1680 E = Ellipsoid
1681 t = E.a, E.b, E.f_, E.f, E.f2, E.n, E.e, E.e2, E.e22, E.e32, \
1682 E.A, E.L, E.R1, E.R2, E.R3, E.Rbiaxial, E.Rtriaxial
1683 return self._instr(name, prec, props=t)
1685 def toTriaxial(self, name=NN):
1686 '''Convert this ellipsoid to a L{Triaxial_}.
1688 @return: A L{Triaxial_} or L{Triaxial} with the C{X} axis
1689 pointing east and C{Z} pointing north.
1691 @see: Method L{Triaxial_.toEllipsoid}.
1692 '''
1693 a, b, n = self.a, self.b, (name or self.name)
1694 return _MODS.triaxials.Triaxial_(a, a, b, name=n) if a < b else \
1695 _MODS.triaxials.Triaxial( a, a, b, name=n)
1697 @Property_RO
1698 def volume(self):
1699 '''Get the ellipsoid's I{volume} (C{meter**3}), M{4 / 3 * PI * R3**3}.
1701 @see: C{R3}.
1702 '''
1703 return Meter3(volume=self.a2 * self.b * PI_3 * _4_0)
1706class Ellipsoid2(Ellipsoid):
1707 '''An L{Ellipsoid} specified by I{equatorial} radius and I{flattening}.
1708 '''
1709 def __init__(self, a, f, name=NN):
1710 '''New L{Ellipsoid2}.
1712 @arg a: Equatorial radius, semi-axis (C{meter}).
1713 @arg f: Flattening: (C{float} < 1.0, negative for I{prolate}).
1714 @kwarg name: Optional, unique name (C{str}).
1716 @raise NameError: Ellipsoid with that B{C{name}} already exists.
1718 @raise ValueError: Invalid B{C{a}} or B{C{f}}.
1720 @note: C{abs(B{f}) < EPS} is forced to C{B{f}=0}, I{spherical}.
1721 Negative C{B{f}} produces a I{prolate} ellipsoid.
1722 '''
1723 if f is None and isinstance(a, Ellipsoid):
1724 Ellipsoid.__init__(self, a.a, f =a.f,
1725 b=a.b, f_=a.f_, name=name)
1726 else:
1727 Ellipsoid.__init__(self, a, f=f, name=name)
1730def _spherical_a_b(a, b):
1731 '''(INTERNAL) C{True} for spherical or invalid C{a} or C{b}.
1732 '''
1733 return a < EPS or b < EPS or fabs(a - b) < EPS
1736def _spherical_f(f):
1737 '''(INTERNAL) C{True} for spherical or invalid C{f}.
1738 '''
1739 return fabs(f) < EPS or f > EPS1
1742def _spherical_f_(f_):
1743 '''(INTERNAL) C{True} for spherical or invalid C{f_}.
1744 '''
1745 return fabs(f_) < EPS or fabs(f_) > _1_EPS
1748def a_b2e(a, b):
1749 '''Return C{e}, the I{1st eccentricity} for a given I{equatorial} and I{polar} radius.
1751 @arg a: Equatorial radius (C{scalar} > 0).
1752 @arg b: Polar radius (C{scalar} > 0).
1754 @return: The I{unsigned}, (1st) eccentricity (C{float} or C{0}),
1755 M{sqrt(1 - (b / a)**2)}.
1757 @note: The result is always I{non-negative} and C{0} for I{near-spherical} ellipsoids.
1758 '''
1759 return Float(e=sqrt(fabs(a_b2e2(a, b))))
1762def a_b2e2(a, b):
1763 '''Return C{e2}, the I{1st eccentricity squared} for a given I{equatorial} and I{polar} radius.
1765 @arg a: Equatorial radius (C{scalar} > 0).
1766 @arg b: Polar radius (C{scalar} > 0).
1768 @return: The I{signed}, (1st) eccentricity I{squared} (C{float} or C{0}),
1769 M{1 - (b / a)**2}.
1771 @note: The result is positive for I{oblate}, negative for I{prolate}
1772 or C{0} for I{near-spherical} ellipsoids.
1773 '''
1774 return Float(e2=_0_0 if _spherical_a_b(a, b) else (_1_0 - (b / a)**2))
1777def a_b2e22(a, b):
1778 '''Return C{e22}, the I{2nd eccentricity squared} for a given I{equatorial} and I{polar} radius.
1780 @arg a: Equatorial radius (C{scalar} > 0).
1781 @arg b: Polar radius (C{scalar} > 0).
1783 @return: The I{signed}, 2nd eccentricity I{squared} (C{float} or C{0}),
1784 M{(a / b)**2 - 1}.
1786 @note: The result is positive for I{oblate}, negative for I{prolate}
1787 or C{0} for I{near-spherical} ellipsoids.
1788 '''
1789 return Float(e22=_0_0 if _spherical_a_b(a, b) else ((a / b)**2 - _1_0))
1792def a_b2e32(a, b):
1793 '''Return C{e32}, the I{3rd eccentricity squared} for a given I{equatorial} and I{polar} radius.
1795 @arg a: Equatorial radius (C{scalar} > 0).
1796 @arg b: Polar radius (C{scalar} > 0).
1798 @return: The I{signed}, 3rd eccentricity I{squared} (C{float} or C{0}),
1799 M{(a**2 - b**2) / (a**2 + b**2)}.
1801 @note: The result is positive for I{oblate}, negative for I{prolate}
1802 or C{0} for I{near-spherical} ellipsoids.
1803 '''
1804 a2, b2 = a**2, b**2
1805 return Float(e32=_0_0 if _spherical_a_b(a2, b2) else ((a2 - b2) / (a2 + b2)))
1808def a_b2f(a, b):
1809 '''Return C{f}, the I{flattening} for a given I{equatorial} and I{polar} radius.
1811 @arg a: Equatorial radius (C{scalar} > 0).
1812 @arg b: Polar radius (C{scalar} > 0).
1814 @return: The flattening (C{float} or C{0}), M{(a - b) / a}.
1816 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0}
1817 for I{near-spherical} ellipsoids.
1818 '''
1819 f = 0 if _spherical_a_b(a, b) else ((a - b) / a)
1820 return _f_0_0 if _spherical_f(f) else Float(f=f)
1823def a_b2f_(a, b):
1824 '''Return C{f_}, the I{inverse flattening} for a given I{equatorial} and I{polar} radius.
1826 @arg a: Equatorial radius (C{scalar} > 0).
1827 @arg b: Polar radius (C{scalar} > 0).
1829 @return: The inverse flattening (C{float} or C{0}), M{a / (a - b)}.
1831 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0}
1832 for I{near-spherical} ellipsoids.
1833 '''
1834 f_ = 0 if _spherical_a_b(a, b) else (a / float(a - b))
1835 return _f__0_0 if _spherical_f_(f_) else Float(f_=f_)
1838def a_b2f2(a, b):
1839 '''Return C{f2}, the I{2nd flattening} for a given I{equatorial} and I{polar} radius.
1841 @arg a: Equatorial radius (C{scalar} > 0).
1842 @arg b: Polar radius (C{scalar} > 0).
1844 @return: The I{signed}, 2nd flattening (C{float} or C{0}), M{(a - b) / b}.
1846 @note: The result is positive for I{oblate}, negative for I{prolate} or C{0}
1847 for I{near-spherical} ellipsoids.
1848 '''
1849 t = 0 if _spherical_a_b(a, b) else float(a - b)
1850 return Float(f2=_0_0 if fabs(t) < EPS else (t / b))
1853def a_b2n(a, b):
1854 '''Return C{n}, the I{3rd flattening} for a given I{equatorial} and I{polar} radius.
1856 @arg a: Equatorial radius (C{scalar} > 0).
1857 @arg b: Polar radius (C{scalar} > 0).
1859 @return: The I{signed}, 3rd flattening (C{float} or C{0}), M{(a - b) / (a + b)}.
1861 @note: The result is positive for I{oblate}, negative for I{prolate}
1862 or C{0} for I{near-spherical} ellipsoids.
1863 '''
1864 t = 0 if _spherical_a_b(a, b) else float(a - b)
1865 return Float(n=_0_0 if fabs(t) < EPS else (t / (a + b)))
1868def a_f2b(a, f):
1869 '''Return C{b}, the I{polar} radius for a given I{equatorial} radius and I{flattening}.
1871 @arg a: Equatorial radius (C{scalar} > 0).
1872 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).
1874 @return: The polar radius (C{float}), M{a * (1 - f)}.
1875 '''
1876 b = a if _spherical_f(f) else (a * (_1_0 - f))
1877 return Radius_(b=a if _spherical_a_b(a, b) else b)
1880def a_f_2b(a, f_):
1881 '''Return C{b}, the I{polar} radius for a given I{equatorial} radius and I{inverse flattening}.
1883 @arg a: Equatorial radius (C{scalar} > 0).
1884 @arg f_: Inverse flattening (C{scalar} >>> 1).
1886 @return: The polar radius (C{float}), M{a * (f_ - 1) / f_}.
1887 '''
1888 b = a if _spherical_f_(f_) else (a * (f_ - _1_0) / f_)
1889 return Radius_(b=a if _spherical_a_b(a, b) else b)
1892def b_f2a(b, f):
1893 '''Return C{a}, the I{equatorial} radius for a given I{polar} radius and I{flattening}.
1895 @arg b: Polar radius (C{scalar} > 0).
1896 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).
1898 @return: The equatorial radius (C{float}), M{b / (1 - f)}.
1899 '''
1900 t = _1_0 - f
1901 a = b if fabs(t < EPS) else (b / t)
1902 return Radius_(a=b if _spherical_a_b(a, b) else a)
1905def b_f_2a(b, f_):
1906 '''Return C{a}, the I{equatorial} radius for a given I{polar} radius and I{inverse flattening}.
1908 @arg b: Polar radius (C{scalar} > 0).
1909 @arg f_: Inverse flattening (C{scalar} >>> 1).
1911 @return: The equatorial radius (C{float}), M{b * f_ / (f_ - 1)}.
1912 '''
1913 t = f_ - _1_0
1914 a = b if _spherical_f_(f_) or fabs(t - f_) < EPS \
1915 or fabs(t) < EPS else (b * f_ / t)
1916 return Radius_(a=b if _spherical_a_b(a, b) else a)
1919def e2f(e):
1920 '''Return C{f}, the I{flattening} for a given I{1st eccentricity}.
1922 @arg e: The (1st) eccentricity (0 <= C{float} < 1)
1924 @return: The flattening (C{float} or C{0}).
1926 @see: Function L{e22f}.
1927 '''
1928 return e22f(e**2)
1931def e22f(e2):
1932 '''Return C{f}, the I{flattening} for a given I{1st eccentricity squared}.
1934 @arg e2: The (1st) eccentricity I{squared}, I{signed} (L{NINF} < C{float} < 1)
1936 @return: The flattening (C{float} or C{0}), M{e2 / (sqrt(e2 - 1) + 1)}.
1937 '''
1938 return Float(f=e2 / (sqrt(_1_0 - e2) + _1_0)) if e2 else _f_0_0
1941def f2e2(f):
1942 '''Return C{e2}, the I{1st eccentricity squared} for a given I{flattening}.
1944 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).
1946 @return: The I{signed}, (1st) eccentricity I{squared} (C{float} < 1),
1947 M{f * (2 - f)}.
1949 @note: The result is positive for I{oblate}, negative for I{prolate}
1950 or C{0} for I{near-spherical} ellipsoids.
1952 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/
1953 html/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening
1954 <https://WikiPedia.org/wiki/Flattening>}.
1955 '''
1956 return Float(e2=_0_0 if _spherical_f(f) else (f * (_2_0 - f)))
1959def f2e22(f):
1960 '''Return C{e22}, the I{2nd eccentricity squared} for a given I{flattening}.
1962 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).
1964 @return: The I{signed}, 2nd eccentricity I{squared} (C{float} > -1 or
1965 C{INF}), M{f * (2 - f) / (1 - f)**2}.
1967 @note: The result is positive for I{oblate}, negative for I{prolate}
1968 or C{0} for near-spherical ellipsoids.
1970 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/
1971 html/classGeographicLib_1_1Ellipsoid.html>}.
1972 '''
1973 # e2 / (1 - e2) == f * (2 - f) / (1 - f)**2
1974 t = (_1_0 - f)**2
1975 return Float(e22=INF if t < EPS else (f2e2(f) / t)) # PYCHOK type
1978def f2e32(f):
1979 '''Return C{e32}, the I{3rd eccentricity squared} for a given I{flattening}.
1981 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).
1983 @return: The I{signed}, 3rd eccentricity I{squared} (C{float}),
1984 M{f * (2 - f) / (1 + (1 - f)**2)}.
1986 @note: The result is positive for I{oblate}, negative for I{prolate}
1987 or C{0} for I{near-spherical} ellipsoids.
1989 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/
1990 html/classGeographicLib_1_1Ellipsoid.html>}.
1991 '''
1992 # e2 / (2 - e2) == f * (2 - f) / (1 + (1 - f)**2)
1993 e2 = f2e2(f)
1994 return Float(e32=e2 / (_2_0 - e2))
1997def f_2f(f_):
1998 '''Return C{f}, the I{flattening} for a given I{inverse flattening}.
2000 @arg f_: Inverse flattening (C{scalar} >>> 1).
2002 @return: The flattening (C{float} or C{0}), M{1 / f_}.
2004 @note: The result is positive for I{oblate}, negative for I{prolate}
2005 or C{0} for I{near-spherical} ellipsoids.
2006 '''
2007 f = 0 if _spherical_f_(f_) else _1_0 / f_
2008 return _f_0_0 if _spherical_f(f) else Float(f=f) # PYCHOK type
2011def f2f_(f):
2012 '''Return C{f_}, the I{inverse flattening} for a given I{flattening}.
2014 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).
2016 @return: The inverse flattening (C{float} or C{0}), M{1 / f}.
2018 @note: The result is positive for I{oblate}, negative for I{prolate}
2019 or C{0} for I{near-spherical} ellipsoids.
2020 '''
2021 f_ = 0 if _spherical_f(f) else _1_0 / f
2022 return _f__0_0 if _spherical_f_(f_) else Float(f_=f_) # PYCHOK type
2025def f2f2(f):
2026 '''Return C{f2}, the I{2nd flattening} for a given I{flattening}.
2028 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).
2030 @return: The I{signed}, 2nd flattening (C{float} or C{INF}), M{f / (1 - f)}.
2032 @note: The result is positive for I{oblate}, negative for I{prolate}
2033 or C{0} for I{near-spherical} ellipsoids.
2035 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/
2036 html/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening
2037 <https://WikiPedia.org/wiki/Flattening>}.
2038 '''
2039 t = _1_0 - f
2040 return Float(f2=_0_0 if _spherical_f(f) else (INF if fabs(t) < EPS
2041 else (f / t))) # PYCHOK type
2044def f2n(f):
2045 '''Return C{n}, the I{3rd flattening} for a given I{flattening}.
2047 @arg f: Flattening (C{scalar} < 1, negative for I{prolate}).
2049 @return: The I{signed}, 3rd flattening (-1 <= C{float} < 1),
2050 M{f / (2 - f)}.
2052 @note: The result is positive for I{oblate}, negative for I{prolate}
2053 or C{0} for I{near-spherical} ellipsoids.
2055 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/
2056 html/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening
2057 <https://WikiPedia.org/wiki/Flattening>}.
2058 '''
2059 return Float(n=_0_0 if _spherical_f(f) else (f / float(_2_0 - f)))
2062def n2e2(n):
2063 '''Return C{e2}, the I{1st eccentricity squared} for a given I{3rd flattening}.
2065 @arg n: The 3rd flattening (-1 <= C{scalar} < 1).
2067 @return: The I{signed}, (1st) eccentricity I{squared} (C{float} or NINF),
2068 M{4 * n / (1 + n)**2}.
2070 @note: The result is positive for I{oblate}, negative for I{prolate}
2071 or C{0} for I{near-spherical} ellipsoids.
2073 @see: U{Flattening<https://WikiPedia.org/wiki/Flattening>}.
2074 '''
2075 t = (n + _1_0)**2
2076 return Float(e2=_0_0 if fabs(n) < EPS else
2077 (NINF if t < EPS else (_4_0 * n / t)))
2080def n2f(n):
2081 '''Return C{f}, the I{flattening} for a given I{3rd flattening}.
2083 @arg n: The 3rd flattening (-1 <= C{scalar} < 1).
2085 @return: The flattening (C{float} or NINF), M{2 * n / (1 + n)}.
2087 @see: U{Eccentricity conversions<https://GeographicLib.SourceForge.io/
2088 html/classGeographicLib_1_1Ellipsoid.html>} and U{Flattening
2089 <https://WikiPedia.org/wiki/Flattening>}.
2090 '''
2091 t = n + _1_0
2092 f = 0 if fabs(n) < EPS else (NINF if t < EPS else (_2_0 * n / t))
2093 return _f_0_0 if _spherical_f(f) else Float(f=f)
2096def n2f_(n):
2097 '''Return C{f_}, the I{inverse flattening} for a given I{3rd flattening}.
2099 @arg n: The 3rd flattening (-1 <= C{scalar} < 1).
2101 @return: The inverse flattening (C{float} or C{0}), M{1 / f}.
2103 @see: L{n2f} and L{f2f_}.
2104 '''
2105 return f2f_(n2f(n))
2108def _normalTo3(px, py, E): # in .height4 above
2109 '''(INTERNAL) Nearest point on a 2-D ellipse in 1st quadrant.
2111 @see: Functions C{.triaxial._normalTo4} and C{-To5}.
2112 '''
2113 a, b = E.a, E.b
2114 if min(px, py, a, b) < EPS0:
2115 raise _AssertionError(px=px, py=py, a=a, b=b, E=E)
2117 a2 = a - b * E.b_a
2118 b2 = b - a * E.a_b
2119 tx = ty = _SQRT2_2
2120 for i in range(10): # max 5
2121 ex = a2 * tx**3
2122 ey = b2 * ty**3
2124 qx = px - ex
2125 qy = py - ey
2126 q = hypot(qx, qy)
2127 if q < EPS0: # PYCHOK no cover
2128 break
2129 r = hypot(ex - tx * a, ey - ty * b) / q
2131 sx, tx = tx, min(_1_0, max(0, (ex + qx * r) / a))
2132 sy, ty = ty, min(_1_0, max(0, (ey + qy * r) / b))
2133 t = hypot(ty, tx)
2134 if t < EPS0: # PYCHOK no cover
2135 break
2136 tx = tx / t # /= chokes PyChecker
2137 ty = ty / t
2138 if max(fabs(sx - tx), fabs(sy - ty)) < EPS:
2139 break
2141 tx *= a / px
2142 ty *= b / py
2143 return tx, ty, i # x and y as fractions
2146class Ellipsoids(_NamedEnum):
2147 '''(INTERNAL) L{Ellipsoid} registry, I{must} be a sub-class
2148 to accommodate the L{_LazyNamedEnumItem} properties.
2149 '''
2150 def _Lazy(self, a, b, f_, **kwds):
2151 '''(INTERNAL) Instantiate the L{Ellipsoid}.
2152 '''
2153 return Ellipsoid(a, b=b, f_=f_, **kwds)
2155Ellipsoids = Ellipsoids(Ellipsoid) # PYCHOK singleton
2156'''Some pre-defined L{Ellipsoid}s, all I{lazily} instantiated.'''
2157# <https://www.GNU.org/software/gama/manual/html_node/Supported-ellipsoids.html>
2158# <https://GSSC.ESA.int/navipedia/index.php/Reference_Frames_in_GNSS>
2159# <https://kb.OSU.edu/dspace/handle/1811/77986>
2160# <https://www.IBM.com/docs/en/db2/11.5?topic=systems-supported-spheroids>
2161# <https://w3.Energistics.org/archive/Epicentre/Epicentre_v3.0/DataModel/LogicalDictionary/StandardValues/ellipsoid.html>
2162# <https://GitHub.com/locationtech/proj4j/blob/master/src/main/java/org/locationtech/proj4j/datum/Ellipsoid.java>
2163Ellipsoids._assert( # <https://WikiPedia.org/wiki/Earth_ellipsoid>
2164 Airy1830 = _lazy(_Airy1830_, *_T(6377563.396, _0_0, 299.3249646)), # b=6356256.909
2165 AiryModified = _lazy(_AiryModified_, *_T(6377340.189, _0_0, 299.3249646)), # b=6356034.448
2166# APL4_9 = _lazy('APL4_9', *_T(6378137.0, _0_0, 298.24985392)), # Appl. Phys. Lab. 1965
2167# ANS = _lazy('ANS', *_T(6378160.0, _0_0, 298.25)), # Australian Nat. Spheroid
2168# AN_SA96 = _lazy('AN_SA96', *_T(6378160.0, _0_0, 298.24985392)), # Australian Nat. South America
2169 Australia1966 = _lazy('Australia1966', *_T(6378160.0, _0_0, 298.25)), # b=6356774.7192
2170 ATS1977 = _lazy('ATS1977', *_T(6378135.0, _0_0, 298.257)), # "Average Terrestrial System"
2171 Bessel1841 = _lazy(_Bessel1841_, *_T(6377397.155, 6356078.962818, 299.152812797)),
2172 BesselModified = _lazy('BesselModified', *_T(6377492.018, _0_0, 299.1528128)),
2173# BesselNamibia = _lazy('BesselNamibia', *_T(6377483.865, _0_0, 299.1528128)),
2174 CGCS2000 = _lazy('CGCS2000', *_T(R_MA, _0_0, 298.257222101)), # BeiDou Coord System (BDC)
2175# Clarke1858 = _lazy('Clarke1858', *_T(6378293.639, _0_0, 294.260676369)),
2176 Clarke1866 = _lazy(_Clarke1866_, *_T(6378206.4, 6356583.8, 294.978698214)),
2177 Clarke1880 = _lazy('Clarke1880', *_T(6378249.145, 6356514.86954978, 293.465)),
2178 Clarke1880IGN = _lazy(_Clarke1880IGN_, *_T(6378249.2, 6356515.0, 293.466021294)),
2179 Clarke1880Mod = _lazy('Clarke1880Mod', *_T(6378249.145, 6356514.96639549, 293.466307656)), # aka Clarke1880Arc
2180 CPM1799 = _lazy('CPM1799', *_T(6375738.7, 6356671.92557493, 334.39)), # Comm. des Poids et Mesures
2181 Delambre1810 = _lazy('Delambre1810', *_T(6376428.0, 6355957.92616372, 311.5)), # Belgium
2182 Engelis1985 = _lazy('Engelis1985', *_T(6378136.05, 6356751.32272154, 298.2566)),
2183# Everest1830 = _lazy('Everest1830', *_T(6377276.345, _0_0, 300.801699997)),
2184# Everest1948 = _lazy('Everest1948', *_T(6377304.063, _0_0, 300.801699997)),
2185# Everest1956 = _lazy('Everest1956', *_T(6377301.243, _0_0, 300.801699997)),
2186 Everest1969 = _lazy('Everest1969', *_T(6377295.664, 6356094.667915, 300.801699997)),
2187 Everest1975 = _lazy('Everest1975', *_T(6377299.151, 6356098.14512013, 300.8017255)),
2188 Fisher1968 = _lazy('Fisher1968', *_T(6378150.0, 6356768.33724438, 298.3)),
2189# Fisher1968Mod = _lazy('Fisher1968Mod', *_T(6378155.0, _0_0, 298.3)),
2190 GEM10C = _lazy('GEM10C', *_T(R_MA, 6356752.31424783, 298.2572236)),
2191 GPES = _lazy('GPES', *_T(6378135.0, 6356750.0, _0_0)), # "Gen. Purpose Earth Spheroid"
2192 GRS67 = _lazy('GRS67', *_T(6378160.0, _0_0, 298.247167427)), # Lucerne b=6356774.516
2193# GRS67Truncated = _lazy('GRS67Truncated', *_T(6378160.0, _0_0, 298.25)),
2194 GRS80 = _lazy(_GRS80_, *_T(R_MA, 6356752.314140347, 298.257222101)), # IUGG, ITRS, ETRS89
2195# Hayford1924 = _lazy('Hayford1924', *_T(6378388.0, 6356911.94612795, None)), # aka Intl1924 f_=297
2196 Helmert1906 = _lazy('Helmert1906', *_T(6378200.0, 6356818.16962789, 298.3)),
2197# Hough1960 = _lazy('Hough1960', *_T(6378270.0, _0_0, 297.0)),
2198 IAU76 = _lazy('IAU76', *_T(6378140.0, _0_0, 298.257)), # Int'l Astronomical Union
2199 IERS1989 = _lazy('IERS1989', *_T(6378136.0, _0_0, 298.257)), # b=6356751.302
2200 IERS1992TOPEX = _lazy('IERS1992TOPEX', *_T(6378136.3, 6356751.61659215, 298.257223563)), # IERS/TOPEX/Poseidon/McCarthy
2201 IERS2003 = _lazy('IERS2003', *_T(6378136.6, 6356751.85797165, 298.25642)),
2202 Intl1924 = _lazy(_Intl1924_, *_T(6378388.0, _0_0, 297.0)), # aka Hayford b=6356911.9462795
2203 Intl1967 = _lazy('Intl1967', *_T(6378157.5, 6356772.2, 298.24961539)), # New Int'l
2204 Krassovski1940 = _lazy(_Krassovski1940_, *_T(6378245.0, 6356863.01877305, 298.3)), # spelling
2205 Krassowsky1940 = _lazy(_Krassowsky1940_, *_T(6378245.0, 6356863.01877305, 298.3)), # spelling
2206# Kaula = _lazy('Kaula', *_T(6378163.0, _0_0, 298.24)), # Kaula 1961
2207# Lerch = _lazy('Lerch', *_T(6378139.0, _0_0, 298.257)), # Lerch 1979
2208 Maupertuis1738 = _lazy('Maupertuis1738', *_T(6397300.0, 6363806.28272251, 191.0)), # France
2209 Mercury1960 = _lazy('Mercury1960', *_T(6378166.0, 6356784.28360711, 298.3)),
2210 Mercury1968Mod = _lazy('Mercury1968Mod', *_T(6378150.0, 6356768.33724438, 298.3)),
2211# MERIT = _lazy('MERIT', *_T(6378137.0, _0_0, 298.257)), # MERIT 1983
2212# NWL10D = _lazy('NWL10D', *_T(6378135.0, _0_0, 298.26)), # Naval Weapons Lab.
2213 NWL1965 = _lazy('NWL1965', *_T(6378145.0, 6356759.76948868, 298.25)), # Naval Weapons Lab.
2214# NWL9D = _lazy('NWL9D', *_T(6378145.0, 6356759.76948868, 298.25)), # NWL1965
2215 OSU86F = _lazy('OSU86F', *_T(6378136.2, 6356751.51693008, 298.2572236)),
2216 OSU91A = _lazy('OSU91A', *_T(6378136.3, 6356751.6165948, 298.2572236)),
2217# Plessis1817 = _lazy('Plessis1817', *_T(6397523.0, 6355863.0, 153.56512242)), # XXX incorrect?
2218 Plessis1817 = _lazy('Plessis1817', *_T(6376523.0, 6355862.93325557, 308.64)), # XXX IGN France 1972
2219 PZ90 = _lazy('PZ90', *_T(6378136.0, _0_0, 298.257839303)), # GLOSNASS PZ-90 and PZ-90.11
2220# SEAsia = _lazy('SEAsia', *_T(6378155.0, _0_0, 298.3)), # SouthEast Asia
2221 SGS85 = _lazy('SGS85', *_T(6378136.0, 6356751.30156878, 298.257)), # Soviet Geodetic System
2222 SoAmerican1969 = _lazy('SoAmerican1969', *_T(6378160.0, 6356774.71919531, 298.25)), # South American
2223 Struve1860 = _lazy('Struve1860', *_T(6378298.3, 6356657.14266956, 294.73)),
2224# Walbeck = _lazy('Walbeck', *_T(6376896.0, _0_0, 302.78)),
2225# WarOffice = _lazy('WarOffice', *_T(6378300.0, _0_0, 296.0)),
2226 WGS60 = _lazy('WGS60', *_T(6378165.0, 6356783.28695944, 298.3)),
2227 WGS66 = _lazy('WGS66', *_T(6378145.0, 6356759.76948868, 298.25)),
2228 WGS72 = _lazy(_WGS72_, *_T(6378135.0, _0_0, 298.26)), # b=6356750.52
2229 WGS84 = _lazy(_WGS84_, *_T(R_MA, _0_0, _f__WGS84)), # GPS b=6356752.3142451793
2230# Prolate = _lazy('Prolate', *_T(6356752.3, R_MA, _0_0)),
2231 Sphere = _lazy(_Sphere_, *_T(R_M, R_M, _0_0)), # pseudo
2232 SphereAuthalic = _lazy('SphereAuthalic', *_T(R_FM, R_FM, _0_0)), # pseudo
2233 SpherePopular = _lazy('SpherePopular', *_T(R_MA, R_MA, _0_0)) # EPSG:3857 Spheroid
2234)
2237if __name__ == '__main__':
2239 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_
2240 from pygeodesy.lazily import printf
2241 from pygeodesy.named import nameof
2243 for E in (Ellipsoids.WGS84, Ellipsoids.GRS80, # NAD83,
2244 Ellipsoids.Sphere, Ellipsoids.SpherePopular,
2245 Ellipsoid(Ellipsoids.WGS84.b, Ellipsoids.WGS84.a, name='_Prolate')):
2246 e = f2n(E.f) - E.n
2247 printf('# %s: %s', _DOT_('Ellipsoids', E.name), E.toStr(prec=10), nl=1)
2248 printf('# e=%s, f_=%s, f=%s, n=%s (%s)', fstr(E.e, prec=13, fmt=Fmt.e),
2249 fstr(E.f_, prec=13, fmt=Fmt.e),
2250 fstr(E.f, prec=13, fmt=Fmt.e),
2251 fstr(E.n, prec=13, fmt=Fmt.e),
2252 fstr(e, prec=9, fmt=Fmt.e))
2253 printf('# %s %s', Ellipsoid.AlphaKs.name, fstr(E.AlphaKs, prec=20))
2254 printf('# %s %s', Ellipsoid.BetaKs.name, fstr(E.BetaKs, prec=20))
2255 printf('# %s %s', nameof(Ellipsoid.KsOrder), E.KsOrder) # property
2257 # __doc__ of this file, force all into registry
2258 t = [NN] + Ellipsoids.toRepr(all=True, asorted=True).split(_NL_)
2259 printf(_NLATvar_.join(i.strip(_COMMA_) for i in t))
2261# **) MIT License
2262#
2263# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
2264#
2265# Permission is hereby granted, free of charge, to any person obtaining a
2266# copy of this software and associated documentation files (the "Software"),
2267# to deal in the Software without restriction, including without limitation
2268# the rights to use, copy, modify, merge, publish, distribute, sublicense,
2269# and/or sell copies of the Software, and to permit persons to whom the
2270# Software is furnished to do so, subject to the following conditions:
2271#
2272# The above copyright notice and this permission notice shall be included
2273# in all copies or substantial portions of the Software.
2274#
2275# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
2276# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
2277# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
2278# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
2279# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
2280# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
2281# OTHER DEALINGS IN THE SOFTWARE.
2283# % python3 -m pygeodesy.ellipsoids
2285# Ellipsoids.WGS84: name='WGS84', a=6378137, b=6356752.3142451793, f_=298.257223563, f=0.0033528107, f2=0.0033640898, n=0.0016792204, e=0.0818191908, e2=0.00669438, e22=0.0067394967, e32=0.0033584313, A=6367449.1458234144, L=10001965.7293127216, R1=6371008.7714150595, R2=6371007.1809184738, R3=6371000.7900091587, Rbiaxial=6367453.6345163295, Rtriaxial=6372797.5559594007
2286# e=8.1819190842622e-02, f_=2.98257223563e+02, f=3.3528106647475e-03, n=1.6792203863837e-03 (0.0e+00)
2287# AlphaKs 0.00083773182062446994, 0.00000076085277735725, 0.00000000119764550324, 0.00000000000242917068, 0.00000000000000571182, 0.0000000000000000148, 0.00000000000000000004, 0.0
2288# BetaKs 0.00083773216405794875, 0.0000000590587015222, 0.00000000016734826653, 0.00000000000021647981, 0.00000000000000037879, 0.00000000000000000072, 0.0, 0.0
2289# KsOrder 8
2291# Ellipsoids.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
2292# e=8.1819191042833e-02, f_=2.98257222101e+02, f=3.3528106811837e-03, n=1.6792203946295e-03 (0.0e+00)
2293# AlphaKs 0.00083773182472890429, 0.00000076085278481561, 0.00000000119764552086, 0.00000000000242917073, 0.00000000000000571182, 0.0000000000000000148, 0.00000000000000000004, 0.0
2294# BetaKs 0.0008377321681623882, 0.00000005905870210374, 0.000000000167348269, 0.00000000000021647982, 0.00000000000000037879, 0.00000000000000000072, 0.0, 0.0
2295# KsOrder 8
2297# Ellipsoids.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
2298# e=0.0e+00, f_=0.0e+00, f=0.0e+00, n=0.0e+00 (0.0e+00)
2299# AlphaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
2300# BetaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
2301# KsOrder 8
2303# Ellipsoids.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
2304# e=0.0e+00, f_=0.0e+00, f=0.0e+00, n=0.0e+00 (0.0e+00)
2305# AlphaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
2306# BetaKs 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
2307# KsOrder 8
2309# Ellipsoids._Prolate: name='_Prolate', a=6356752.3142451793, b=6378137, f_=-297.257223563, f=-0.0033640898, f2=-0.0033528107, n=-0.0016792204, e=0.0820944379, e2=-0.0067394967, e22=-0.00669438, e32=-0.0033584313, A=6367449.1458234144, L=10035500.5204500332, R1=6363880.5428301189, R2=6363878.9413582645, R3=6363872.5644020075, Rbiaxial=6367453.6345163295, Rtriaxial=6362105.2243882557
2310# e=8.2094437949696e-02, f_=-2.97257223563e+02, f=-3.3640898209765e-03, n=-1.6792203863837e-03 (0.0e+00)
2311# AlphaKs -0.00084149152514366627, 0.00000076653480614871, -0.00000000120934503389, 0.0000000000024576225, -0.00000000000000578863, 0.00000000000000001502, -0.00000000000000000004, 0.0
2312# BetaKs -0.00084149187224351817, 0.00000005842735196773, -0.0000000001680487236, 0.00000000000021706261, -0.00000000000000038002, 0.00000000000000000073, -0.0, 0.0
2313# KsOrder 8