diff options
Diffstat (limited to 'ast/erfa/plan94.c')
| -rw-r--r-- | ast/erfa/plan94.c | 523 |
1 files changed, 0 insertions, 523 deletions
diff --git a/ast/erfa/plan94.c b/ast/erfa/plan94.c deleted file mode 100644 index ef878dd..0000000 --- a/ast/erfa/plan94.c +++ /dev/null @@ -1,523 +0,0 @@ -#include "erfa.h" - -int eraPlan94(double date1, double date2, int np, double pv[2][3]) -/* -** - - - - - - - - - - -** e r a P l a n 9 4 -** - - - - - - - - - - -** -** Approximate heliocentric position and velocity of a nominated major -** planet: Mercury, Venus, EMB, Mars, Jupiter, Saturn, Uranus or -** Neptune (but not the Earth itself). -** -** Given: -** date1 double TDB date part A (Note 1) -** date2 double TDB date part B (Note 1) -** np int planet (1=Mercury, 2=Venus, 3=EMB, 4=Mars, -** 5=Jupiter, 6=Saturn, 7=Uranus, 8=Neptune) -** -** Returned (argument): -** pv double[2][3] planet p,v (heliocentric, J2000.0, AU,AU/d) -** -** Returned (function value): -** int status: -1 = illegal NP (outside 1-8) -** 0 = OK -** +1 = warning: year outside 1000-3000 -** +2 = warning: failed to converge -** -** Notes: -** -** 1) The date date1+date2 is in the TDB time scale (in practice TT can -** be used) and is a Julian Date, apportioned in any convenient way -** between the two arguments. For example, JD(TDB)=2450123.7 could -** be expressed in any of these ways, among others: -** -** date1 date2 -** -** 2450123.7 0.0 (JD method) -** 2451545.0 -1421.3 (J2000 method) -** 2400000.5 50123.2 (MJD method) -** 2450123.5 0.2 (date & time method) -** -** The JD method is the most natural and convenient to use in cases -** where the loss of several decimal digits of resolution is -** acceptable. The J2000 method is best matched to the way the -** argument is handled internally and will deliver the optimum -** resolution. The MJD method and the date & time methods are both -** good compromises between resolution and convenience. The limited -** accuracy of the present algorithm is such that any of the methods -** is satisfactory. -** -** 2) If an np value outside the range 1-8 is supplied, an error status -** (function value -1) is returned and the pv vector set to zeroes. -** -** 3) For np=3 the result is for the Earth-Moon Barycenter. To obtain -** the heliocentric position and velocity of the Earth, use instead -** the ERFA function eraEpv00. -** -** 4) On successful return, the array pv contains the following: -** -** pv[0][0] x } -** pv[0][1] y } heliocentric position, AU -** pv[0][2] z } -** -** pv[1][0] xdot } -** pv[1][1] ydot } heliocentric velocity, AU/d -** pv[1][2] zdot } -** -** The reference frame is equatorial and is with respect to the -** mean equator and equinox of epoch J2000.0. -** -** 5) The algorithm is due to J.L. Simon, P. Bretagnon, J. Chapront, -** M. Chapront-Touze, G. Francou and J. Laskar (Bureau des -** Longitudes, Paris, France). From comparisons with JPL -** ephemeris DE102, they quote the following maximum errors -** over the interval 1800-2050: -** -** L (arcsec) B (arcsec) R (km) -** -** Mercury 4 1 300 -** Venus 5 1 800 -** EMB 6 1 1000 -** Mars 17 1 7700 -** Jupiter 71 5 76000 -** Saturn 81 13 267000 -** Uranus 86 7 712000 -** Neptune 11 1 253000 -** -** Over the interval 1000-3000, they report that the accuracy is no -** worse than 1.5 times that over 1800-2050. Outside 1000-3000 the -** accuracy declines. -** -** Comparisons of the present function with the JPL DE200 ephemeris -** give the following RMS errors over the interval 1960-2025: -** -** position (km) velocity (m/s) -** -** Mercury 334 0.437 -** Venus 1060 0.855 -** EMB 2010 0.815 -** Mars 7690 1.98 -** Jupiter 71700 7.70 -** Saturn 199000 19.4 -** Uranus 564000 16.4 -** Neptune 158000 14.4 -** -** Comparisons against DE200 over the interval 1800-2100 gave the -** following maximum absolute differences. (The results using -** DE406 were essentially the same.) -** -** L (arcsec) B (arcsec) R (km) Rdot (m/s) -** -** Mercury 7 1 500 0.7 -** Venus 7 1 1100 0.9 -** EMB 9 1 1300 1.0 -** Mars 26 1 9000 2.5 -** Jupiter 78 6 82000 8.2 -** Saturn 87 14 263000 24.6 -** Uranus 86 7 661000 27.4 -** Neptune 11 2 248000 21.4 -** -** 6) The present ERFA re-implementation of the original Simon et al. -** Fortran code differs from the original in the following respects: -** -** * C instead of Fortran. -** -** * The date is supplied in two parts. -** -** * The result is returned only in equatorial Cartesian form; -** the ecliptic longitude, latitude and radius vector are not -** returned. -** -** * The result is in the J2000.0 equatorial frame, not ecliptic. -** -** * More is done in-line: there are fewer calls to subroutines. -** -** * Different error/warning status values are used. -** -** * A different Kepler's-equation-solver is used (avoiding -** use of double precision complex). -** -** * Polynomials in t are nested to minimize rounding errors. -** -** * Explicit double constants are used to avoid mixed-mode -** expressions. -** -** None of the above changes affects the result significantly. -** -** 7) The returned status indicates the most serious condition -** encountered during execution of the function. Illegal np is -** considered the most serious, overriding failure to converge, -** which in turn takes precedence over the remote date warning. -** -** Called: -** eraAnp normalize angle into range 0 to 2pi -** -** Reference: Simon, J.L, Bretagnon, P., Chapront, J., -** Chapront-Touze, M., Francou, G., and Laskar, J., -** Astron. Astrophys. 282, 663 (1994). -** -** Copyright (C) 2013-2016, NumFOCUS Foundation. -** Derived, with permission, from the SOFA library. See notes at end of file. -*/ -{ -/* Gaussian constant */ - static const double GK = 0.017202098950; - -/* Sin and cos of J2000.0 mean obliquity (IAU 1976) */ - static const double SINEPS = 0.3977771559319137; - static const double COSEPS = 0.9174820620691818; - -/* Maximum number of iterations allowed to solve Kepler's equation */ - static const int KMAX = 10; - - int jstat, i, k; - double t, da, dl, de, dp, di, dom, dmu, arga, argl, am, - ae, dae, ae2, at, r, v, si2, xq, xp, tl, xsw, - xcw, xm2, xf, ci2, xms, xmc, xpxq2, x, y, z; - -/* Planetary inverse masses */ - static const double amas[] = { 6023600.0, /* Mercury */ - 408523.5, /* Venus */ - 328900.5, /* EMB */ - 3098710.0, /* Mars */ - 1047.355, /* Jupiter */ - 3498.5, /* Saturn */ - 22869.0, /* Uranus */ - 19314.0 }; /* Neptune */ - -/* -** Tables giving the mean Keplerian elements, limited to t^2 terms: -** -** a semi-major axis (AU) -** dlm mean longitude (degree and arcsecond) -** e eccentricity -** pi longitude of the perihelion (degree and arcsecond) -** dinc inclination (degree and arcsecond) -** omega longitude of the ascending node (degree and arcsecond) -*/ - - static const double a[][3] = { - { 0.3870983098, 0.0, 0.0 }, /* Mercury */ - { 0.7233298200, 0.0, 0.0 }, /* Venus */ - { 1.0000010178, 0.0, 0.0 }, /* EMB */ - { 1.5236793419, 3e-10, 0.0 }, /* Mars */ - { 5.2026032092, 19132e-10, -39e-10 }, /* Jupiter */ - { 9.5549091915, -0.0000213896, 444e-10 }, /* Saturn */ - { 19.2184460618, -3716e-10, 979e-10 }, /* Uranus */ - { 30.1103868694, -16635e-10, 686e-10 } /* Neptune */ - }; - - static const double dlm[][3] = { - { 252.25090552, 5381016286.88982, -1.92789 }, - { 181.97980085, 2106641364.33548, 0.59381 }, - { 100.46645683, 1295977422.83429, -2.04411 }, - { 355.43299958, 689050774.93988, 0.94264 }, - { 34.35151874, 109256603.77991, -30.60378 }, - { 50.07744430, 43996098.55732, 75.61614 }, - { 314.05500511, 15424811.93933, -1.75083 }, - { 304.34866548, 7865503.20744, 0.21103 } - }; - - static const double e[][3] = { - { 0.2056317526, 0.0002040653, -28349e-10 }, - { 0.0067719164, -0.0004776521, 98127e-10 }, - { 0.0167086342, -0.0004203654, -0.0000126734 }, - { 0.0934006477, 0.0009048438, -80641e-10 }, - { 0.0484979255, 0.0016322542, -0.0000471366 }, - { 0.0555481426, -0.0034664062, -0.0000643639 }, - { 0.0463812221, -0.0002729293, 0.0000078913 }, - { 0.0094557470, 0.0000603263, 0.0 } - }; - - static const double pi[][3] = { - { 77.45611904, 5719.11590, -4.83016 }, - { 131.56370300, 175.48640, -498.48184 }, - { 102.93734808, 11612.35290, 53.27577 }, - { 336.06023395, 15980.45908, -62.32800 }, - { 14.33120687, 7758.75163, 259.95938 }, - { 93.05723748, 20395.49439, 190.25952 }, - { 173.00529106, 3215.56238, -34.09288 }, - { 48.12027554, 1050.71912, 27.39717 } - }; - - static const double dinc[][3] = { - { 7.00498625, -214.25629, 0.28977 }, - { 3.39466189, -30.84437, -11.67836 }, - { 0.0, 469.97289, -3.35053 }, - { 1.84972648, -293.31722, -8.11830 }, - { 1.30326698, -71.55890, 11.95297 }, - { 2.48887878, 91.85195, -17.66225 }, - { 0.77319689, -60.72723, 1.25759 }, - { 1.76995259, 8.12333, 0.08135 } - }; - - static const double omega[][3] = { - { 48.33089304, -4515.21727, -31.79892 }, - { 76.67992019, -10008.48154, -51.32614 }, - { 174.87317577, -8679.27034, 15.34191 }, - { 49.55809321, -10620.90088, -230.57416 }, - { 100.46440702, 6362.03561, 326.52178 }, - { 113.66550252, -9240.19942, -66.23743 }, - { 74.00595701, 2669.15033, 145.93964 }, - { 131.78405702, -221.94322, -0.78728 } - }; - -/* Tables for trigonometric terms to be added to the mean elements of */ -/* the semi-major axes */ - - static const double kp[][9] = { - { 69613, 75645, 88306, 59899, 15746, 71087, 142173, 3086, 0 }, - { 21863, 32794, 26934, 10931, 26250, 43725, 53867, 28939, 0 }, - { 16002, 21863, 32004, 10931, 14529, 16368, 15318, 32794, 0 }, - { 6345, 7818, 15636, 7077, 8184, 14163, 1107, 4872, 0 }, - { 1760, 1454, 1167, 880, 287, 2640, 19, 2047, 1454 }, - { 574, 0, 880, 287, 19, 1760, 1167, 306, 574 }, - { 204, 0, 177, 1265, 4, 385, 200, 208, 204 }, - { 0, 102, 106, 4, 98, 1367, 487, 204, 0 } - }; - - static const double ca[][9] = { - { 4, -13, 11, -9, -9, -3, -1, 4, 0 }, - { -156, 59, -42, 6, 19, -20, -10, -12, 0 }, - { 64, -152, 62, -8, 32, -41, 19, -11, 0 }, - { 124, 621, -145, 208, 54, -57, 30, 15, 0 }, - { -23437, -2634, 6601, 6259, -1507,-1821, 2620, -2115, -1489 }, - { 62911,-119919, 79336,17814,-24241,12068, 8306, -4893, 8902 }, - { 389061,-262125,-44088, 8387,-22976,-2093, -615, -9720, 6633 }, - { -412235,-157046,-31430,37817, -9740, -13, -7449, 9644, 0 } - }; - - static const double sa[][9] = { - { -29, -1, 9, 6, -6, 5, 4, 0, 0 }, - { -48, -125, -26, -37, 18, -13, -20, -2, 0 }, - { -150, -46, 68, 54, 14, 24, -28, 22, 0 }, - { -621, 532, -694, -20, 192, -94, 71, -73, 0 }, - { -14614,-19828, -5869, 1881, -4372, -2255, 782, 930, 913 }, - { 139737, 0, 24667, 51123, -5102, 7429, -4095, -1976, -9566 }, - { -138081, 0, 37205,-49039,-41901,-33872,-27037,-12474, 18797 }, - { 0, 28492,133236, 69654, 52322,-49577,-26430, -3593, 0 } - }; - -/* Tables giving the trigonometric terms to be added to the mean */ -/* elements of the mean longitudes */ - - static const double kq[][10] = { - { 3086,15746,69613,59899,75645,88306, 12661, 2658, 0, 0 }, - { 21863,32794,10931, 73, 4387,26934, 1473, 2157, 0, 0 }, - { 10,16002,21863,10931, 1473,32004, 4387, 73, 0, 0 }, - { 10, 6345, 7818, 1107,15636, 7077, 8184, 532, 10, 0 }, - { 19, 1760, 1454, 287, 1167, 880, 574, 2640, 19, 1454 }, - { 19, 574, 287, 306, 1760, 12, 31, 38, 19, 574 }, - { 4, 204, 177, 8, 31, 200, 1265, 102, 4, 204 }, - { 4, 102, 106, 8, 98, 1367, 487, 204, 4, 102 } - }; - - static const double cl[][10] = { - { 21, -95, -157, 41, -5, 42, 23, 30, 0, 0 }, - { -160, -313, -235, 60, -74, -76, -27, 34, 0, 0 }, - { -325, -322, -79, 232, -52, 97, 55, -41, 0, 0 }, - { 2268, -979, 802, 602, -668, -33, 345, 201, -55, 0 }, - { 7610, -4997,-7689,-5841,-2617, 1115,-748,-607, 6074, 354 }, - { -18549, 30125,20012, -730, 824, 23,1289,-352, -14767, -2062 }, - { -135245,-14594, 4197,-4030,-5630,-2898,2540,-306, 2939, 1986 }, - { 89948, 2103, 8963, 2695, 3682, 1648, 866,-154, -1963, -283 } - }; - - static const double sl[][10] = { - { -342, 136, -23, 62, 66, -52, -33, 17, 0, 0 }, - { 524, -149, -35, 117, 151, 122, -71, -62, 0, 0 }, - { -105, -137, 258, 35, -116, -88,-112, -80, 0, 0 }, - { 854, -205, -936, -240, 140, -341, -97, -232, 536, 0 }, - { -56980, 8016, 1012, 1448,-3024,-3710, 318, 503, 3767, 577 }, - { 138606,-13478,-4964, 1441,-1319,-1482, 427, 1236, -9167, -1918 }, - { 71234,-41116, 5334,-4935,-1848, 66, 434, -1748, 3780, -701 }, - { -47645, 11647, 2166, 3194, 679, 0,-244, -419, -2531, 48 } - }; - -/*--------------------------------------------------------------------*/ - -/* Validate the planet number. */ - if ((np < 1) || (np > 8)) { - jstat = -1; - - /* Reset the result in case of failure. */ - for (k = 0; k < 2; k++) { - for (i = 0; i < 3; i++) { - pv[k][i] = 0.0; - } - } - - } else { - - /* Decrement the planet number to start at zero. */ - np--; - - /* Time: Julian millennia since J2000.0. */ - t = ((date1 - ERFA_DJ00) + date2) / ERFA_DJM; - - /* OK status unless remote date. */ - jstat = fabs(t) <= 1.0 ? 0 : 1; - - /* Compute the mean elements. */ - da = a[np][0] + - (a[np][1] + - a[np][2] * t) * t; - dl = (3600.0 * dlm[np][0] + - (dlm[np][1] + - dlm[np][2] * t) * t) * ERFA_DAS2R; - de = e[np][0] + - ( e[np][1] + - e[np][2] * t) * t; - dp = eraAnpm((3600.0 * pi[np][0] + - (pi[np][1] + - pi[np][2] * t) * t) * ERFA_DAS2R); - di = (3600.0 * dinc[np][0] + - (dinc[np][1] + - dinc[np][2] * t) * t) * ERFA_DAS2R; - dom = eraAnpm((3600.0 * omega[np][0] + - (omega[np][1] + - omega[np][2] * t) * t) * ERFA_DAS2R); - - /* Apply the trigonometric terms. */ - dmu = 0.35953620 * t; - for (k = 0; k < 8; k++) { - arga = kp[np][k] * dmu; - argl = kq[np][k] * dmu; - da += (ca[np][k] * cos(arga) + - sa[np][k] * sin(arga)) * 1e-7; - dl += (cl[np][k] * cos(argl) + - sl[np][k] * sin(argl)) * 1e-7; - } - arga = kp[np][8] * dmu; - da += t * (ca[np][8] * cos(arga) + - sa[np][8] * sin(arga)) * 1e-7; - for (k = 8; k < 10; k++) { - argl = kq[np][k] * dmu; - dl += t * (cl[np][k] * cos(argl) + - sl[np][k] * sin(argl)) * 1e-7; - } - dl = fmod(dl, ERFA_D2PI); - - /* Iterative soln. of Kepler's equation to get eccentric anomaly. */ - am = dl - dp; - ae = am + de * sin(am); - k = 0; - dae = 1.0; - while (k < KMAX && fabs(dae) > 1e-12) { - dae = (am - ae + de * sin(ae)) / (1.0 - de * cos(ae)); - ae += dae; - k++; - if (k == KMAX-1) jstat = 2; - } - - /* True anomaly. */ - ae2 = ae / 2.0; - at = 2.0 * atan2(sqrt((1.0 + de) / (1.0 - de)) * sin(ae2), - cos(ae2)); - - /* Distance (AU) and speed (radians per day). */ - r = da * (1.0 - de * cos(ae)); - v = GK * sqrt((1.0 + 1.0 / amas[np]) / (da * da * da)); - - si2 = sin(di / 2.0); - xq = si2 * cos(dom); - xp = si2 * sin(dom); - tl = at + dp; - xsw = sin(tl); - xcw = cos(tl); - xm2 = 2.0 * (xp * xcw - xq * xsw); - xf = da / sqrt(1 - de * de); - ci2 = cos(di / 2.0); - xms = (de * sin(dp) + xsw) * xf; - xmc = (de * cos(dp) + xcw) * xf; - xpxq2 = 2 * xp * xq; - - /* Position (J2000.0 ecliptic x,y,z in AU). */ - x = r * (xcw - xm2 * xp); - y = r * (xsw + xm2 * xq); - z = r * (-xm2 * ci2); - - /* Rotate to equatorial. */ - pv[0][0] = x; - pv[0][1] = y * COSEPS - z * SINEPS; - pv[0][2] = y * SINEPS + z * COSEPS; - - /* Velocity (J2000.0 ecliptic xdot,ydot,zdot in AU/d). */ - x = v * (( -1.0 + 2.0 * xp * xp) * xms + xpxq2 * xmc); - y = v * (( 1.0 - 2.0 * xq * xq) * xmc - xpxq2 * xms); - z = v * (2.0 * ci2 * (xp * xms + xq * xmc)); - - /* Rotate to equatorial. */ - pv[1][0] = x; - pv[1][1] = y * COSEPS - z * SINEPS; - pv[1][2] = y * SINEPS + z * COSEPS; - - } - -/* Return the status. */ - return jstat; - -} -/*---------------------------------------------------------------------- -** -** -** Copyright (C) 2013-2016, NumFOCUS Foundation. -** All rights reserved. -** -** This library is derived, with permission, from the International -** Astronomical Union's "Standards of Fundamental Astronomy" library, -** available from http://www.iausofa.org. -** -** The ERFA version is intended to retain identical functionality to -** the SOFA library, but made distinct through different function and -** file names, as set out in the SOFA license conditions. The SOFA -** original has a role as a reference standard for the IAU and IERS, -** and consequently redistribution is permitted only in its unaltered -** state. The ERFA version is not subject to this restriction and -** therefore can be included in distributions which do not support the -** concept of "read only" software. -** -** Although the intent is to replicate the SOFA API (other than -** replacement of prefix names) and results (with the exception of -** bugs; any that are discovered will be fixed), SOFA is not -** responsible for any errors found in this version of the library. -** -** If you wish to acknowledge the SOFA heritage, please acknowledge -** that you are using a library derived from SOFA, rather than SOFA -** itself. -** -** -** TERMS AND CONDITIONS -** -** Redistribution and use in source and binary forms, with or without -** modification, are permitted provided that the following conditions -** are met: -** -** 1 Redistributions of source code must retain the above copyright -** notice, this list of conditions and the following disclaimer. -** -** 2 Redistributions in binary form must reproduce the above copyright -** notice, this list of conditions and the following disclaimer in -** the documentation and/or other materials provided with the -** distribution. -** -** 3 Neither the name of the Standards Of Fundamental Astronomy Board, -** the International Astronomical Union nor the names of its -** contributors may be used to endorse or promote products derived -** from this software without specific prior written permission. -** -** THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -** "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -** LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -** FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -** COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, -** INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, -** BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; -** LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER -** CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT -** LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN -** ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -** POSSIBILITY OF SUCH DAMAGE. -** -*/ |