/* Naval Observatory Vector Astrometry Software (NOVAS) C Edition, Version 3.1 solsys3.c: Self-contained Sun-Earth ephemeris U. S. Naval Observatory Astronomical Applications Dept. Washington, DC http://www.usno.navy.mil/USNO/astronomical-applications */ #ifndef _NOVAS_ #include "novas.h" #endif #include /* Additional function prototype. */ void sun_eph (double jd, double *ra, double *dec, double *dis); /********solarsystem */ short int solarsystem (double tjd, short int body, short int origin, double *position, double *velocity) /* ------------------------------------------------------------------------ PURPOSE: Provides the position and velocity of the Earth at epoch 'tjd' by evaluating a closed-form theory without reference to an external file. This function can also provide the position and velocity of the Sun. REFERENCES: Kaplan, G. H. "NOVAS: Naval Observatory Vector Astrometry Subroutines"; USNO internal document dated 20 Oct 1988; revised 15 Mar 1990. Explanatory Supplement to The Astronomical Almanac (1992). INPUT ARGUMENTS: tjd (double) TDB Julian date. body (short int) Body identification number. Set 'body' = 0 or 'body' = 1 or 'body' = 10 for the Sun. Set 'body' = 2 or 'body' = 3 for the Earth. origin (short int) Origin code = 0 ... solar system barycenter = 1 ... center of mass of the Sun OUTPUT ARGUMENTS: position[3] (double) Position vector of 'body' at 'tjd'; equatorial rectangular coordinates in AU referred to the mean equator and equinox of J2000.0. velocity[3] (double) Velocity vector of 'body' at 'tjd'; equatorial rectangular system referred to the mean equator and equinox of J2000.0, in AU/Day. RETURNED VALUE: (short int) 0...Everything OK. 1...Input Julian date ('tjd') out of range. 2...Invalid value of 'body'. GLOBALS USED: T0, TWOPI. novascon.c FUNCTIONS CALLED: sun_eph solsys3.c radec2vector novas.c precession novas.c sin math.h cos math.h fabs math.h fmod math.h VER./DATE/ PROGRAMMER: V1.0/05-96/JAB (USNO/AA) Convert to C; substitute new theory of Sun. V1.1/06-98/JAB (USNO/AA) Updated planetary masses & mean elements. V1.2/04-09/JAB (USNO/AA) Updated value of obliquity at J2000.0. V1.3/11-09/WKP (USNO/AA) Updated output argument names to 'position' and 'velocity'. V1.4/11-09/JAB (USNO/AA) Update barycenter computation. V1.5/02-11/JLB (USNO/AA) Reformatted description of origin for consistency with other documentation. NOTES: 1. This function is the "C" version of Fortran NOVAS routine 'solsys' version 3. ------------------------------------------------------------------------ */ { short int ierr = 0; short int i; /* The arrays below contain masses and orbital elements for the four largest planets -- Jupiter, Saturn, Uranus, and Neptune -- (see Explanatory Supplement (1992), p. 316) with angles in radians. These data are used for barycenter computations only. */ static const double pm[4] = {1047.349, 3497.898, 22903.0, 19412.2}; static const double pa[4] = {5.203363, 9.537070, 19.191264, 30.068963}; static const double pe[4] = {0.048393, 0.054151, 0.047168, 0.008586}; static const double pj[4] = {0.022782, 0.043362, 0.013437, 0.030878}; static const double po[4] = {1.755036, 1.984702, 1.295556, 2.298977}; static const double pw[4] = {0.257503, 1.613242, 2.983889, 0.784898}; static const double pl[4] = {0.600470, 0.871693, 5.466933, 5.321160}; static const double pn[4] = {1.450138e-3, 5.841727e-4, 2.047497e-4, 1.043891e-4}; /* 'obl' is the obliquity of ecliptic at epoch J2000.0 in degrees. */ static const double obl = 23.4392794444; static double tlast = 0.0; static double tmass, a[3][4], b[3][4], pbary[3], vbary[3]; double oblr, se, ce, si, ci, sn, cn, sw, cw, p1, p2, p3, q1, q2, q3, roote, qjd, ras, decs, diss, pos1[3], p[3][3], e, mlon, ma, u, sinu, cosu, anr, pplan[3], vplan[3], f; /* Initialize constants. Initial value of 'tmass' is mass of Sun plus four inner planets. */ if (tlast < 1.0) { tmass = 1.0 + 5.977e-6; oblr = obl * TWOPI / 360.0; se = sin (oblr); ce = cos (oblr); for (i = 0; i < 4; i++) { tmass += 1.0 / pm[i]; /* Compute sine and cosine of orbital angles. */ si = sin (pj[i]); ci = cos (pj[i]); sn = sin (po[i]); cn = cos (po[i]); sw = sin (pw[i] - po[i]); cw = cos (pw[i] - po[i]); /* Compute p and q vectors (see Brouwer & Clemence (1961), Methods of Celestial Mechanics, pp. 35-36.) */ p1 = cw * cn - sw * sn * ci; p2 = ( cw * sn + sw * cn * ci) * ce - sw * si * se; p3 = ( cw * sn + sw * cn * ci) * se + sw * si * ce; q1 = -sw * cn - cw * sn * ci; q2 = (-sw * sn + cw * cn * ci) * ce - cw * si * se; q3 = (-sw * sn + cw * cn * ci) * se + cw * si * ce; roote = sqrt (1.0 - pe[i] * pe[i]); a[0][i] = pa[i] * p1; a[1][i] = pa[i] * p2; a[2][i] = pa[i] * p3; b[0][i] = pa[i] * roote * q1; b[1][i] = pa[i] * roote * q2; b[2][i] = pa[i] * roote * q3; } tlast = 1.0; } /* Check if input Julian date is within range (within 3 centuries of J2000). */ if ((tjd < 2340000.5) || (tjd > 2560000.5)) return (ierr = 1); /* Form heliocentric coordinates of the Sun or Earth, depending on 'body'. Velocities are obtained from crude numerical differentiation. */ if ((body == 0) || (body == 1) || (body == 10)) /* Sun */ { for (i = 0; i < 3; i++) position[i] = velocity[i] = 0.0; } else if ((body == 2) || (body == 3)) /* Earth */ { for (i = 0; i < 3; i++) { qjd = tjd + (double) (i - 1) * 0.1; sun_eph (qjd, &ras,&decs,&diss); radec2vector (ras,decs,diss, pos1); precession (qjd,pos1,T0, position); p[i][0] = -position[0]; p[i][1] = -position[1]; p[i][2] = -position[2]; } for (i = 0; i < 3; i++) { position[i] = p[1][i]; velocity[i] = (p[2][i] - p[0][i]) / 0.2; } } else return (ierr = 2); /* If 'origin' = 0, move origin to solar system barycenter. Solar system barycenter coordinates are computed from Keplerian approximations of the coordinates of the four largest planets. */ if (origin == 0) { if (fabs (tjd - tlast) >= 1.0e-06) { for (i = 0; i < 3; i++) pbary[i] = vbary[i] = 0.0; /* The following loop cycles once for each of the four planets. */ for (i = 0; i < 4; i++) { /* Compute mean longitude, mean anomaly, and eccentric anomaly. */ e = pe[i]; mlon = pl[i] + pn[i] * (tjd - T0); ma = fmod ((mlon - pw[i]), TWOPI); u = ma + e * sin (ma) + 0.5 * e * e * sin (2.0 * ma); sinu = sin (u); cosu = cos (u); /* Compute velocity factor. */ anr = pn[i] / (1.0 - e * cosu); /* Compute planet's position and velocity wrt eq & eq J2000. */ pplan[0] = a[0][i] * (cosu - e) + b[0][i] * sinu; pplan[1] = a[1][i] * (cosu - e) + b[1][i] * sinu; pplan[2] = a[2][i] * (cosu - e) + b[2][i] * sinu; vplan[0] = anr * (-a[0][i] * sinu + b[0][i] * cosu); vplan[1] = anr * (-a[1][i] * sinu + b[1][i] * cosu); vplan[2] = anr * (-a[2][i] * sinu + b[2][i] * cosu); /* Compute mass factor and add in to total displacement. */ f = 1.0 / (pm[i] * tmass); pbary[0] += pplan[0] * f; pbary[1] += pplan[1] * f; pbary[2] += pplan[2] * f; vbary[0] += vplan[0] * f; vbary[1] += vplan[1] * f; vbary[2] += vplan[2] * f; } tlast = tjd; } for (i = 0; i < 3; i++) { position[i] -= pbary[i]; velocity[i] -= vbary[i]; } } return (ierr); } /********solarsystem_hp */ short int solarsystem_hp (double tjd[2], short int body, short int origin, double *position, double *velocity) /* ------------------------------------------------------------------------ PURPOSE: Provides the position and velocity of the Earth at epoch 'tjd' to highest precision. However, this is a dummy routine especially for use with the NOVAS-C file 'solsys3.c'. ** See notes below.** REFERENCES: None. INPUT ARGUMENTS: tjd[2] (double) Two-element array containing the Julian date, which may be split any way (although the first element is usually the "integer" part, and the second element is the "fractional" part). Julian date is on the TDB or "T_eph" time scale. body (short int) Body identification number. Set 'body' = 0 or 'body' = 1 or 'body' = 10 for the Sun. Set 'body' = 2 or 'body' = 3 for the Earth. origin (short int) Origin code = 0 ... solar system barycenter = 1 ... center of mass of the Sun OUTPUT ARGUMENTS: position[3] (double) Position vector of 'body' at 'tjd'; equatorial rectangular coordinates in AU referred to the mean equator and equinox of J2000.0. velocity[3] (double) Velocity vector of 'body' at 'tjd'; equatorial rectangular system referred to the mean equator and equinox of J2000.0, in AU/Day. RETURNED VALUE: (short int) 0...Everything OK. 1...Input Julian date ('tjd') out of range. 2...Invalid value of 'body'. 3...This version of 'solarsystem' not valid for use with NOVAS-C. GLOBALS USED: None. FUNCTIONS CALLED: solarsystem (solsys3.c) VER./DATE/ PROGRAMMER: V1.0/01-07/JAB (USNO/AA) V1.1/11-09/WKP (USNO/AA) Updated output argument names to 'position' and 'velocity'. v1.2/07-10/JLB (USNO/AA) Revised note 3 per Help desk msg V1.3/02-11/JLB (USNO/AA) Reformatted description of origin for consistency with other documentation. NOTES: 1. Function 'solarsystem' in file 'solsys3.c' provides the position and velocity of the Earth and Sun by evaluating a low-precison, closed-form theory without reference to an external file. NOVAS-C does not provide a high-precision counterpart of this function that works without requiring an external file. Thus, this function ('solarsystem_hp') is essentially a dummy function that acts according to the value of variable 'action', which is set below. 2. If 'action' = 1 (the default), this function returns an error code of '3' indicating that the function does not provide high-precision position and velocity of the Earth and Sun. This is the appropriate action. 3. If 'action' = 2, this function simply calls function 'solarsystem' and returns the low-precision position and velocity. An error code of '0' (no error) is also returned. This action may be useful for code testing purposes, but is not appropriate or recommended for normal use of NOVAS-C. To access this option, the code must be changed manually to set 'action' = 2. 4. Use alternate versions of 'solarsystem' (located in the various 'solsysn.c' files, where n is an integer) when highest precision is needed. ------------------------------------------------------------------------ */ { /* Set the value of 'action' according to the rules explained under NOTES in the prolog. */ short int action = 1; short int error = 0; short int i; double jd, p[3], v[3]; switch (action) { case 2: /* Return low-precision position and velocity */ jd = tjd[0] + tjd[1]; if ((error = solarsystem (jd,body,origin, p,v)) != 0) { return (error); } else { for (i = 0; i < 3; i++) { position[i] = p[i]; velocity[i] = v[i]; } } break; case 1: /* Set an error flag and return. */ default: return (error = 3); } return (error); } /********sun_eph */ void sun_eph (double jd, double *ra, double *dec, double *dis) /* ------------------------------------------------------------------------ PURPOSE: To compute equatorial spherical coordinates of Sun referred to the mean equator and equinox of date. REFERENCES: Bretagnon, P. and Simon, J.L. (1986). Planetary Programs and Tables from -4000 to + 2800. (Richmond, VA: Willmann-Bell). Kaplan, G.H. (2005). US Naval Observatory Circular 179. INPUT ARGUMENTS: jd (double) Julian date on TDT or ET time scale. OUTPUT ARGUMENTS: ra (double) Right ascension referred to mean equator and equinox of date (hours). dec (double) Declination referred to mean equator and equinox of date (degrees). dis (double) Geocentric distance (AU). RETURNED VALUE: None. GLOBALS USED: T0, TWOPI, ASEC2RAD FUNCTIONS CALLED: sin math.h cos math.h asin math.h atan2 math.h VER./DATE/ PROGRAMMER: V1.0/08-94/JAB (USNO/AA) V1.1/05-96/JAB (USNO/AA): Compute mean coordinates instead of apparent. V1.2/01-07/JAB (USNO/AA): Use 'ASEC2RAD' instead of 'RAD2SEC'. V1.3/04-09/JAB (USNO/AA): Update the equation for mean obliquity of the ecliptic, and correct longitude based on a linear fit to DE405 in the interval 1900-2100 (see notes). NOTES: 1. Quoted accuracy is 2.0 + 0.03 * T^2 arcsec, where T is measured in units of 1000 years from J2000.0. See reference. 2. The obliquity equation is updated to equation 5.12 of the second reference. 3. The linear fit to DE405 primarily corrects for the difference between "old" (Lieske) and "new" (IAU 2006) precession. The difference, new - old, is -0.3004 arcsec/cy. ------------------------------------------------------------------------ */ { short int i; double sum_lon = 0.0; double sum_r = 0.0; const double factor = 1.0e-07; double u, arg, lon, lat, t, emean, sin_lon; struct sun_con { double l; double r; double alpha; double nu; }; static const struct sun_con con[50] = {{403406.0, 0.0, 4.721964, 1.621043}, {195207.0, -97597.0, 5.937458, 62830.348067}, {119433.0, -59715.0, 1.115589, 62830.821524}, {112392.0, -56188.0, 5.781616, 62829.634302}, { 3891.0, -1556.0, 5.5474 , 125660.5691 }, { 2819.0, -1126.0, 1.5120 , 125660.9845 }, { 1721.0, -861.0, 4.1897 , 62832.4766 }, { 0.0, 941.0, 1.163 , 0.813 }, { 660.0, -264.0, 5.415 , 125659.310 }, { 350.0, -163.0, 4.315 , 57533.850 }, { 334.0, 0.0, 4.553 , -33.931 }, { 314.0, 309.0, 5.198 , 777137.715 }, { 268.0, -158.0, 5.989 , 78604.191 }, { 242.0, 0.0, 2.911 , 5.412 }, { 234.0, -54.0, 1.423 , 39302.098 }, { 158.0, 0.0, 0.061 , -34.861 }, { 132.0, -93.0, 2.317 , 115067.698 }, { 129.0, -20.0, 3.193 , 15774.337 }, { 114.0, 0.0, 2.828 , 5296.670 }, { 99.0, -47.0, 0.52 , 58849.27 }, { 93.0, 0.0, 4.65 , 5296.11 }, { 86.0, 0.0, 4.35 , -3980.70 }, { 78.0, -33.0, 2.75 , 52237.69 }, { 72.0, -32.0, 4.50 , 55076.47 }, { 68.0, 0.0, 3.23 , 261.08 }, { 64.0, -10.0, 1.22 , 15773.85 }, { 46.0, -16.0, 0.14 , 188491.03 }, { 38.0, 0.0, 3.44 , -7756.55 }, { 37.0, 0.0, 4.37 , 264.89 }, { 32.0, -24.0, 1.14 , 117906.27 }, { 29.0, -13.0, 2.84 , 55075.75 }, { 28.0, 0.0, 5.96 , -7961.39 }, { 27.0, -9.0, 5.09 , 188489.81 }, { 27.0, 0.0, 1.72 , 2132.19 }, { 25.0, -17.0, 2.56 , 109771.03 }, { 24.0, -11.0, 1.92 , 54868.56 }, { 21.0, 0.0, 0.09 , 25443.93 }, { 21.0, 31.0, 5.98 , -55731.43 }, { 20.0, -10.0, 4.03 , 60697.74 }, { 18.0, 0.0, 4.27 , 2132.79 }, { 17.0, -12.0, 0.79 , 109771.63 }, { 14.0, 0.0, 4.24 , -7752.82 }, { 13.0, -5.0, 2.01 , 188491.91 }, { 13.0, 0.0, 2.65 , 207.81 }, { 13.0, 0.0, 4.98 , 29424.63 }, { 12.0, 0.0, 0.93 , -7.99 }, { 10.0, 0.0, 2.21 , 46941.14 }, { 10.0, 0.0, 3.59 , -68.29 }, { 10.0, 0.0, 1.50 , 21463.25 }, { 10.0, -9.0, 2.55 , 157208.40 }}; /* Define the time units 'u', measured in units of 10000 Julian years from J2000.0, and 't', measured in Julian centuries from J2000.0. */ u = (jd - T0) / 3652500.0; t = u * 100.0; /* Compute longitude and distance terms from the series. */ for (i = 0; i < 50; i++) { arg = con[i].alpha + con[i].nu * u; sum_lon += con[i].l * sin (arg); sum_r += con[i].r * cos (arg); } /* Compute longitude, latitude, and distance referred to mean equinox and ecliptic of date. Apply correction to longitude based on a linear fit to DE405 in the interval 1900-2100. */ lon = 4.9353929 + 62833.1961680 * u + factor * sum_lon; lon += ((-0.1371679461 - 0.2918293271 * t) * ASEC2RAD); lon = fmod (lon, TWOPI); if (lon < 0.0) lon += TWOPI; lat = 0.0; *dis = 1.0001026 + factor * sum_r; /* Compute mean obliquity of the ecliptic. */ emean = (84381.406 + (-46.836769 + (-0.0001831 + 0.00200340 * t) * t) * t) * ASEC2RAD; /* Compute equatorial spherical coordinates referred to the mean equator and equinox of date. */ sin_lon = sin (lon); *ra = atan2 ((cos (emean) * sin_lon), cos (lon)) * RAD2DEG; *ra = fmod (*ra, 360.0); if (*ra < 0.0) *ra += 360.0; *ra = *ra / 15.0; *dec = asin (sin (emean) * sin_lon) * RAD2DEG; return; }