loon: solarsystem.h File Reference

00040           :    
00041       Provides the position and velocity of the Earth at epoch 'tjd'
00042       by evaluating a closed-form theory without reference to an 
00043       external file.  This function can also provide the position
00044       and velocity of the Sun.
00045 
00046    REFERENCES: 
00047       Kaplan, G. H. "NOVAS: Naval Observatory Vector Astrometry
00048          Subroutines"; USNO internal document dated 20 Oct 1988;
00049          revised 15 Mar 1990.
00050       Explanatory Supplement to The Astronomical Almanac (1992).
00051 
00052    INPUT
00053    ARGUMENTS:
00054       tjd (double)
00055          TDB Julian date.
00056       body (short int)
00057          Body identification number.
00058          Set 'body' = 0 or 'body' = 1 or 'body' = 10 for the Sun.
00059          Set 'body' = 2 or 'body' = 3 for the Earth.
00060       origin (short int)
00061          Origin code; solar system barycenter   = 0,
00062                       center of mass of the Sun = 1.
00063 
00064    OUTPUT
00065    ARGUMENTS:
00066       pos[3] (double)
00067          Position vector of 'body' at 'tjd'; equatorial rectangular
00068          coordinates in AU referred to the mean equator and equinox
00069          of J2000.0.
00070       vel[3] (double)
00071          Velocity vector of 'body' at 'tjd'; equatorial rectangular
00072          system referred to the mean equator and equinox of J2000.0,
00073          in AU/Day.
00074 
00075    RETURNED
00076    VALUE:
00077       (short int)
00078          0...Everything OK.
00079          1...Input Julian date ('tjd') out of range.
00080          2...Invalid value of 'body'.
00081 
00082    GLOBALS
00083    USED:
00084       T0, TWOPI.
00085 
00086    FUNCTIONS
00087    CALLED:
00088       sun_eph          solsys3.c
00089       radec2vector     novas.c
00090       precession       novas.c
00091       sin              math.h
00092       cos              math.h
00093       fabs             math.h
00094       fmod             math.h
00095 
00096    VER./DATE/
00097    PROGRAMMER:
00098       V1.0/05-96/JAB (USNO/AA) Convert to C; substitute new theory of
00099                                Sun.
00100       V1.1/06-98/JAB (USNO/AA) Updated planetary masses & mean elements.
00101 
00102    NOTES:
00103       1. This function is the "C" version of Fortran NOVAS routine
00104       'solsys' version 3.
00105 
00106 ------------------------------------------------------------------------
00107 */
00108 {
00109    short int ierr = 0;
00110    short int i;
00111 
00112 /*
00113    The arrays below contain data for the four largest planets.  Masses
00114    are DE405 values; elements are from Explanatory Supplement, p. 316). 
00115    These data are used for barycenter computations only.
00116 */
00117 
00118    const double pm[4] = {1047.349, 3497.898, 22903.0, 19412.2};
00119    const double pa[4] = {5.203363, 9.537070, 19.191264, 30.068963};
00120    const double pl[4] = {0.600470, 0.871693, 5.466933, 5.321160};
00121    const double pn[4] = {1.450138e-3, 5.841727e-4, 2.047497e-4, 
00122                          1.043891e-4};
00123 
00124 /*
00125    'obl' is the obliquity of ecliptic at epoch J2000.0 in degrees.
00126 */
00127 
00128    const double obl = 23.43929111;
00129 
00130    static double tlast = 0.0;
00131    static double sine, cose, tmass, pbary[3], vbary[3];
00132 
00133    double oblr, qjd, ras, decs, diss, pos1[3], p[3][3], dlon, sinl,
00134       cosl, x, y, z, xdot, ydot, zdot, f;
00135 
00136 /*
00137    Initialize constants.
00138 */
00139 
00140    if (tlast == 0.0)
00141    {
00142       oblr = obl * TWOPI / 360.0;
00143       sine = sin (oblr);
00144       cose = cos (oblr);
00145       tmass = 1.0;
00146       for (i = 0; i < 4; i++)
00147          tmass += 1.0 / pm[i];
00148       tlast = 1.0;
00149    }
00150 
00151 /*
00152    Check if input Julian date is within range.
00153 */
00154 
00155    if ((tjd < 2340000.5) || (tjd > 2560000.5))
00156       return (ierr = 1);
00157 
00158 /*
00159    Form helicentric coordinates of the Sun or Earth, depending on
00160    'body'.
00161 */
00162 
00163    if ((body == 0) || (body == 1) || (body == 10))
00164       for (i = 0; i < 3; i++)
00165          pos[i] = vel[i] = 0.0;
00166 
00167     else if ((body == 2) || (body == 3))
00168     {
00169       for (i = 0; i < 3; i++)
00170       {
00171          qjd = tjd + (double) (i - 1) * 0.1;
00172          sun_eph (qjd, &ras,&decs,&diss);
00173          radec2vector (ras,decs,diss, pos1);
00174          precession (qjd,pos1,T0, pos);
00175          p[i][0] = -pos[0];
00176          p[i][1] = -pos[1];
00177          p[i][2] = -pos[2];
00178       }
00179       for (i = 0; i < 3; i++)
00180       {
00181          pos[i] = p[1][i];
00182          vel[i] = (p[2][i] - p[0][i]) / 0.2;
00183       }
00184     }
00185 
00186     else
00187       return (ierr = 2);
00188            
00189 /*
00190    If 'origin' = 0, move origin to solar system barycenter.
00191 
00192    Solar system barycenter coordinates are computed from rough
00193    approximations of the coordinates of the four largest planets.
00194 */
00195 
00196    if (origin == 0)
00197    {
00198       if (fabs (tjd - tlast) >= 1.0e-06)
00199       {
00200          for (i = 0; i < 3; i++)
00201             pbary[i] = vbary[i] = 0.0;
00202 
00203 /*
00204    The following loop cycles once for each of the four planets.
00205 
00206    'sinl' and 'cosl' are the sine and cosine of the planet's mean
00207    longitude.
00208 */
00209 
00210          for (i = 0; i < 4; i++)
00211          {
00212             dlon = pl[i] + pn[i] * (tjd - T0);
00213             dlon = fmod (dlon, TWOPI);
00214             sinl = sin (dlon);
00215             cosl = cos (dlon);
00216 
00217             x =  pa[i] * cosl;
00218             y =  pa[i] * sinl * cose;
00219             z =  pa[i] * sinl * sine;
00220             xdot = -pa[i] * pn[i] * sinl;
00221             ydot =  pa[i] * pn[i] * cosl * cose;
00222             zdot =  pa[i] * pn[i] * cosl * sine;
00223 
00224             f = 1.0 / (pm[i] * tmass);
00225 
00226             pbary[0] += x * f;
00227             pbary[1] += y * f;
00228             pbary[2] += z * f;
00229             vbary[0] += xdot * f;
00230             vbary[1] += ydot * f;
00231             vbary[2] += zdot * f;
00232          }
00233 
00234          tlast = tjd;
00235       }
00236 
00237       for (i = 0; i < 3; i++)
00238       {
00239          pos[i] -= pbary[i];
00240          vel[i] -= vbary[i];
00241       }
00242    }
00243 
00244    return (ierr);
00245 }
solarsystem.h File Reference

Functions

Function Documentation


Functions
short int	solarsystem (double tjd, short int body, short int origin, double pos, double vel)