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(). : 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; solar system barycenter = 0, center of mass of the Sun = 1. OUTPUT ARGUMENTS: pos[3] (double) Position vector of 'body' at 'tjd'; equatorial rectangular coordinates in AU referred to the mean equator and equinox of J2000.0. vel[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. 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. 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 data for the four largest planets. Masses are DE405 values; elements are from Explanatory Supplement, p. 316). These data are used for barycenter computations only. */ const double pm[4] = {1047.349, 3497.898, 22903.0, 19412.2}; const double pa[4] = {5.203363, 9.537070, 19.191264, 30.068963}; const double pl[4] = {0.600470, 0.871693, 5.466933, 5.321160}; 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. */ const double obl = 23.43929111; static double tlast = 0.0; static double sine, cose, tmass, pbary[3], vbary[3]; double oblr, qjd, ras, decs, diss, pos1[3], p[3][3], dlon, sinl, cosl, x, y, z, xdot, ydot, zdot, f; /* Initialize constants. */ if (tlast == 0.0) { oblr = obl * TWOPI / 360.0; sine = sin (oblr); cose = cos (oblr); tmass = 1.0; for (i = 0; i < 4; i++) tmass += 1.0 / pm[i]; tlast = 1.0; } /* Check if input Julian date is within range. */ if ((tjd < 2340000.5) || (tjd > 2560000.5)) return (ierr = 1); /* Form helicentric coordinates of the Sun or Earth, depending on 'body'. */ if ((body == 0) || (body == 1) || (body == 10)) for (i = 0; i < 3; i++) pos[i] = vel[i] = 0.0; else if ((body == 2) || (body == 3)) { 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, pos); p[i][0] = -pos[0]; p[i][1] = -pos[1]; p[i][2] = -pos[2]; } for (i = 0; i < 3; i++) { pos[i] = p[1][i]; vel[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 rough 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. 'sinl' and 'cosl' are the sine and cosine of the planet's mean longitude. */ for (i = 0; i < 4; i++) { dlon = pl[i] + pn[i] * (tjd - T0); dlon = fmod (dlon, TWOPI); sinl = sin (dlon); cosl = cos (dlon); x = pa[i] * cosl; y = pa[i] * sinl * cose; z = pa[i] * sinl * sine; xdot = -pa[i] * pn[i] * sinl; ydot = pa[i] * pn[i] * cosl * cose; zdot = pa[i] * pn[i] * cosl * sine; f = 1.0 / (pm[i] * tmass); pbary[0] += x * f; pbary[1] += y * f; pbary[2] += z * f; vbary[0] += xdot * f; vbary[1] += ydot * f; vbary[2] += zdot * f; } tlast = tjd; } for (i = 0; i < 3; i++) { pos[i] -= pbary[i]; vel[i] -= vbary[i]; } } return (ierr); }

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

Definition at line 249 of file solsys3.c. References dis, 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 } }

---