solsys3.c File Reference

#include "novas.h"
#include <math.h>

Go to the source code of this file.

Functions

void sun_eph (double jd, double *ra, double *dec, double *dis)
short int solarsystem (double tjd, short int body, short int origin, double *pos, double *vel)

Function Documentation

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

Definition at line 34 of file solsys3.c.

References precession(), radec2vector(), sun_eph(), T0, and TWOPI.

Referenced by ephemeris(), and get_earth().

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 }

void sun_eph ( double  jd,
double *  ra,
double *  dec,
double *  dis 
)

Definition at line 249 of file solsys3.c.

References RAD2DEG, RAD2SEC, T0, and TWOPI.

Referenced by solarsystem().

00255           :
00256       To compute equatorial spherical coordinates of Sun referred to
00257       the mean equator and equinox of date.
00258 
00259    REFERENCES:
00260       Bretagnon, P. and Simon, J.L. (1986).  Planetary Programs and
00261          Tables from -4000 to + 2800. (Richmond, VA: Willmann-Bell).
00262 
00263    INPUT
00264    ARGUMENTS:
00265       jd (double)
00266          Julian date on TDT or ET time scale.
00267 
00268    OUTPUT
00269    ARGUMENTS:
00270       ra (double)
00271          Right ascension referred to mean equator and equinox of date
00272          (hours).
00273       dec (double)
00274          Declination referred to mean equator and equinox of date 
00275          (degrees).
00276       dis (double)
00277          Geocentric distance (AU).
00278 
00279    RETURNED
00280    VALUE:
00281       None.
00282 
00283    GLOBALS
00284    USED:
00285       T0
00286       TWOPI
00287       RAD2DEG
00288 
00289    FUNCTIONS
00290    CALLED:
00291       sin           math.h
00292       cos           math.h
00293       asin          math.h
00294       atan2         math.h
00295 
00296    VER./DATE/
00297    PROGRAMMER:
00298       V1.0/08-94/JAB (USNO/AA)
00299       V1.1/05-96/JAB (USNO/AA): Compute mean coordinates instead of
00300                                 apparent.
00301 
00302    NOTES:
00303       1. Quoted accuracy is 2.0 + 0.03 * T^2 arcsec, where T is
00304       measured in units of 1000 years from J2000.0.  See reference.
00305 
00306 ------------------------------------------------------------------------
00307 */
00308 {
00309    short int i;
00310 
00311    double sum_lon = 0.0;
00312    double sum_r = 0.0;
00313    const double factor = 1.0e-07;
00314    double u, arg, lon, lat, t, t2, emean, sin_lon;
00315 
00316    struct sun_con
00317    {
00318    double l;
00319    double r;
00320    double alpha;
00321    double nu;
00322    };
00323 
00324    static const struct sun_con con[50] =
00325       {{403406.0,      0.0, 4.721964,     1.621043},
00326        {195207.0, -97597.0, 5.937458, 62830.348067}, 
00327        {119433.0, -59715.0, 1.115589, 62830.821524}, 
00328        {112392.0, -56188.0, 5.781616, 62829.634302}, 
00329        {  3891.0,  -1556.0, 5.5474  , 125660.5691 }, 
00330        {  2819.0,  -1126.0, 1.5120  , 125660.9845 }, 
00331        {  1721.0,   -861.0, 4.1897  ,  62832.4766 }, 
00332        {     0.0,    941.0, 1.163   ,      0.813  }, 
00333        {   660.0,   -264.0, 5.415   , 125659.310  }, 
00334        {   350.0,   -163.0, 4.315   ,  57533.850  }, 
00335        {   334.0,      0.0, 4.553   ,    -33.931  }, 
00336        {   314.0,    309.0, 5.198   , 777137.715  }, 
00337        {   268.0,   -158.0, 5.989   ,  78604.191  }, 
00338        {   242.0,      0.0, 2.911   ,      5.412  }, 
00339        {   234.0,    -54.0, 1.423   ,  39302.098  }, 
00340        {   158.0,      0.0, 0.061   ,    -34.861  }, 
00341        {   132.0,    -93.0, 2.317   , 115067.698  }, 
00342        {   129.0,    -20.0, 3.193   ,  15774.337  }, 
00343        {   114.0,      0.0, 2.828   ,   5296.670  }, 
00344        {    99.0,    -47.0, 0.52    ,  58849.27   }, 
00345        {    93.0,      0.0, 4.65    ,   5296.11   }, 
00346        {    86.0,      0.0, 4.35    ,  -3980.70   }, 
00347        {    78.0,    -33.0, 2.75    ,  52237.69   }, 
00348        {    72.0,    -32.0, 4.50    ,  55076.47   }, 
00349        {    68.0,      0.0, 3.23    ,    261.08   }, 
00350        {    64.0,    -10.0, 1.22    ,  15773.85   }, 
00351        {    46.0,    -16.0, 0.14    ,  188491.03  }, 
00352        {    38.0,      0.0, 3.44    ,   -7756.55  }, 
00353        {    37.0,      0.0, 4.37    ,     264.89  }, 
00354        {    32.0,    -24.0, 1.14    ,  117906.27  }, 
00355        {    29.0,    -13.0, 2.84    ,   55075.75  }, 
00356        {    28.0,      0.0, 5.96    ,   -7961.39  }, 
00357        {    27.0,     -9.0, 5.09    ,  188489.81  }, 
00358        {    27.0,      0.0, 1.72    ,    2132.19  }, 
00359        {    25.0,    -17.0, 2.56    ,  109771.03  }, 
00360        {    24.0,    -11.0, 1.92    ,   54868.56  }, 
00361        {    21.0,      0.0, 0.09    ,   25443.93  }, 
00362        {    21.0,     31.0, 5.98    ,  -55731.43  }, 
00363        {    20.0,    -10.0, 4.03    ,   60697.74  }, 
00364        {    18.0,      0.0, 4.27    ,    2132.79  }, 
00365        {    17.0,    -12.0, 0.79    ,  109771.63  }, 
00366        {    14.0,      0.0, 4.24    ,   -7752.82  }, 
00367        {    13.0,     -5.0, 2.01    ,  188491.91  }, 
00368        {    13.0,      0.0, 2.65    ,     207.81  }, 
00369        {    13.0,      0.0, 4.98    ,   29424.63  }, 
00370        {    12.0,      0.0, 0.93    ,      -7.99  }, 
00371        {    10.0,      0.0, 2.21    ,   46941.14  }, 
00372        {    10.0,      0.0, 3.59    ,     -68.29  }, 
00373        {    10.0,      0.0, 1.50    ,   21463.25  }, 
00374        {    10.0,     -9.0, 2.55    ,  157208.40  }};
00375 
00376 /*
00377    Define the time unit 'u', measured in units of 10000 Julian years
00378    from J2000.0.
00379 */
00380 
00381    u = (jd - T0) / 3652500.0;
00382    
00383 /*
00384    Compute longitude and distance terms from the series.
00385 */
00386 
00387    for (i = 0; i < 50; i++)
00388    {
00389       arg = con[i].alpha + con[i].nu * u;
00390       sum_lon += con[i].l * sin (arg);
00391       sum_r += con[i].r * cos (arg);
00392    }
00393 
00394 /*
00395    Compute longitude, latitude, and distance referred to mean equinox
00396    and ecliptic of date.
00397 */
00398 
00399    lon = 4.9353929 + 62833.1961680 * u + factor * sum_lon;
00400 
00401    lon = fmod (lon, TWOPI);
00402    if (lon < 0.0)
00403       lon += TWOPI;
00404 
00405    lat = 0.0;
00406 
00407    *dis = 1.0001026 + factor * sum_r;
00408 
00409 /*
00410    Compute mean obliquity of the ecliptic.
00411 */
00412 
00413    t = u * 100.0;
00414    t2 = t * t;
00415    emean = (0.001813 * t2 * t - 0.00059 * t2 - 46.8150 * t +
00416       84381.448) / RAD2SEC;
00417 
00418 /*
00419    Compute equatorial spherical coordinates referred to the mean equator 
00420    and equinox of date.
00421 */
00422 
00423    sin_lon = sin (lon);
00424    *ra = atan2 ((cos (emean) * sin_lon), cos (lon)) * RAD2DEG;
00425    *ra = fmod (*ra, 360.0);
00426    if (*ra < 0.0)
00427       *ra += 360.0;
00428    *ra = *ra / 15.0;
00429 
00430    *dec = asin (sin (emean) * sin_lon) * RAD2DEG;
00431    
00432    return;
00433 }
}


Generated on 8 Jul 2019 for loon by  doxygen 1.6.1