ecd@cs.umn.edu@ncs-med.UUCP (Elwood C. Downey) (03/11/90)
Posting-number: Volume 11, Issue 4
Submitted-by: ecd@cs.umn.edu@ncs-med.UUCP (Elwood C. Downey)
Archive-name: ephem4.12/part03
# This is the first line of a "shell archive" file.
# This means it contains several files that can be extracted into
# the current directory when run with the sh shell, as follows:
# sh < this_file_name
# This is file 3.
echo x mainmenu.c
sed -e 's/^X//' << 'EOFxEOF' > mainmenu.c
X/* printing routines for the main (upper) screen.
X */
X
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X#include "circum.h"
X#include "screen.h"
X
X/* #define PC_GRAPHICS */
X#ifdef PC_GRAPHICS
X#define JOINT 207
X#define VERT 179
X#define HORIZ 205
X#else
X#define JOINT '-'
X#define VERT '|'
X#define HORIZ '-'
X#endif
X
Xmm_borders()
X{
X char line[NC+1], *lp;
X register i;
X
X lp = line;
X for (i = 0; i < NC; i++)
X *lp++ = HORIZ;
X *lp = '\0';
X f_string (R_PLANTAB-1, 1, line);
X for (i = R_TOP; i < R_PLANTAB-1; i++)
X f_char (i, COL2-2, VERT);
X f_char (R_PLANTAB-1, COL2-2, JOINT);
X for (i = R_TOP; i < R_PLANTAB-1; i++)
X f_char (i, COL3-2, VERT);
X f_char (R_PLANTAB-1, COL3-2, JOINT);
X for (i = R_LST; i < R_PLANTAB-1; i++)
X f_char (i, COL4-2, VERT);
X f_char (R_PLANTAB-1, COL4-2, JOINT);
X}
X
X/* print the permanent labels, on the top menu and the planet names
X * in the bottom section.
X */
Xmm_labels()
X{
X f_string (R_TZN, C_TZN, "LT");
X f_string (R_UT, C_UT, "UTC");
X f_string (R_JD, C_JD, "JulianDate");
X f_string (R_WATCH, C_WATCH, "Watch");
X f_string (R_SRCH, C_SRCH, "Search");
X f_string (R_PLOT, C_PLOT, "Plot");
X f_string (R_ALTM, C_ALTM, "Menu");
X
X f_string (R_LST, C_LST, "LST");
X f_string (R_DAWN, C_DAWN, "Dawn");
X f_string (R_DUSK, C_DUSK, "Dusk");
X f_string (R_LON, C_LON, "NiteLn");
X f_string (R_NSTEP, C_NSTEP, "NStep");
X f_string (R_STPSZ, C_STPSZ, "StpSz");
X
X f_string (R_LAT, C_LAT, "Lat");
X f_string (R_LONG, C_LONG, "Long");
X f_string (R_HEIGHT, C_HEIGHT, "Elev");
X f_string (R_TEMP, C_TEMP, "Temp");
X f_string (R_PRES, C_PRES, "AtmPr");
X f_string (R_TZONE, C_TZONE, "TZ");
X f_string (R_EPOCH, C_EPOCH, "Epoch");
X
X f_string (R_PLANTAB, C_OBJ, "Ob");
X f_string (R_SUN, C_OBJ, "Su");
X f_string (R_MOON, C_OBJ, "Mo");
X f_string (R_MERCURY, C_OBJ, "Me");
X f_string (R_VENUS, C_OBJ, "Ve");
X f_string (R_MARS, C_OBJ, "Ma");
X f_string (R_JUPITER, C_OBJ, "Ju");
X f_string (R_SATURN, C_OBJ, "Sa");
X f_string (R_URANUS, C_OBJ, "Ur");
X f_string (R_NEPTUNE, C_OBJ, "Ne");
X f_string (R_PLUTO, C_OBJ, "Pl");
X f_string (R_OBJX, C_OBJ, "X");
X}
X
X/* print all the time/date/where related stuff: the Now structure.
X * print in a nice order, based on the field locations, as much as possible.
X */
Xmm_now (np, all)
XNow *np;
Xint all;
X{
X char buf[32];
X double lmjd = mjd - tz/24.0;
X double jd = mjd + 2415020L;
X double tmp;
X
X sprintf (buf, "%-3.3s", tznm);
X f_string (R_TZN, C_TZN, buf);
X f_time (R_LT, C_LT, mjd_hr(lmjd));
X f_date (R_LD, C_LD, mjd_day(lmjd));
X
X f_time (R_UT, C_UTV, mjd_hr(mjd));
X f_date (R_UD, C_UD, mjd_day(mjd));
X
X sprintf (buf, "%14.5f", jd);
X (void) flog_log (R_JD, C_JDV, jd);
X f_string (R_JD, C_JDV, buf);
X
X now_lst (np, &tmp);
X f_time (R_LST, C_LSTV, tmp);
X
X if (all) {
X f_gangle (R_LAT, C_LATV, lat);
X f_gangle (R_LONG, C_LONGV, -lng); /* + west */
X
X tmp = height * 2.093e7; /* want to see ft, not earth radii */
X sprintf (buf, "%5g ft", tmp);
X (void) flog_log (R_HEIGHT, C_HEIGHTV, tmp);
X f_string (R_HEIGHT, C_HEIGHTV, buf);
X
X tmp = 9./5.*temp + 32.0; /* want to see degrees F, not C */
X (void) flog_log (R_TEMP, C_TEMPV, tmp);
X sprintf (buf, "%6g F", tmp);
X f_string (R_TEMP, C_TEMPV, buf);
X
X tmp = pressure / 33.86; /* want to see in. Hg, not mBar */
X (void) flog_log (R_PRES, C_PRESV, tmp);
X sprintf (buf, "%5.2f in", tmp);
X f_string (R_PRES, C_PRESV, buf);
X
X f_signtime (R_TZONE, C_TZONEV, tz);
X
X if (epoch == EOD)
X f_string (R_EPOCH, C_EPOCHV, "(OfDate)");
X else {
X mjd_year (epoch, &tmp);
X f_double (R_EPOCH, C_EPOCHV, "%8.1f", tmp);
X }
X }
X
X /* print the calendar for local day, if new month/year. */
X mm_calendar (np, all > 1);
X}
X
X/* display dawn/dusk/length-of-night times.
X */
Xmm_twilight (np, force)
XNow *np;
Xint force;
X{
X double dusk, dawn;
X double tmp;
X int status;
X
X if (!twilight_cir (np, &dawn, &dusk, &status) && !force)
X return;
X
X if (status != 0) {
X f_blanks (R_DAWN, C_DAWNV, 5);
X f_blanks (R_DUSK, C_DUSKV, 5);
X f_string (R_LON, C_LONV, "-----");
X return;
X }
X
X f_mtime (R_DAWN, C_DAWNV, dawn);
X f_mtime (R_DUSK, C_DUSKV, dusk);
X tmp = dawn - dusk; range (&tmp, 24.0);
X f_mtime (R_LON, C_LONV, tmp);
X}
X
Xmm_newcir (y)
Xint y;
X{
X static char ncmsg[] = "NEW CIRCUMSTANCES";
X static char nomsg[] = " ";
X static int last_y = -1;
X
X if (y != last_y) {
X f_string (R_NEWCIR, C_NEWCIR, y ? ncmsg : nomsg);
X last_y = y;
X }
X}
X
Xstatic
Xmm_calendar (np, force)
XNow *np;
Xint force;
X{
X static char *mnames[] = {
X "January", "February", "March", "April", "May", "June",
X "July", "August", "September", "October", "November", "December"
X };
X static int last_m, last_y;
X static double last_tz;
X char str[64];
X int m, y;
X double d;
X int f, nd;
X int r;
X double jd0;
X
X /* get local m/d/y. do nothing if still same month and not forced. */
X mjd_cal (mjd_day(mjd-tz/24.0), &m, &d, &y);
X if (m == last_m && y == last_y && tz == last_tz && !force)
X return;
X last_m = m;
X last_y = y;
X last_tz = tz;
X
X /* find day of week of first day of month */
X cal_mjd (m, 1.0, y, &jd0);
X mjd_dow (jd0, &f);
X if (f < 0) {
X /* can't figure it out - too hard before Gregorian */
X int i;
X for (i = 8; --i >= 0; )
X f_string (R_CAL+i, C_CAL, " ");
X return;
X }
X
X /* print header */
X f_blanks (R_CAL, C_CAL, 20);
X sprintf (str, "%s %4d", mnames[m-1], y);
X f_string (R_CAL, C_CAL + (20 - (strlen(mnames[m-1]) + 5))/2, str);
X f_string (R_CAL+1, C_CAL, "Su Mo Tu We Th Fr Sa");
X
X /* find number of days in this month */
X mjd_dpm (jd0, &nd);
X
X /* print the calendar */
X for (r = 0; r < 6; r++) {
X char row[7*3+1], *rp = row;
X int c;
X for (c = 0; c < 7; c++) {
X int i = r*7+c;
X if (i < f || i >= f + nd)
X sprintf (rp, " ");
X else
X sprintf (rp, "%2d ", i-f+1);
X rp += 3;
X }
X row[sizeof(row)-2] = '\0'; /* don't print last blank; causes wrap*/
X f_string (R_CAL+2+r, C_CAL, row);
X }
X
X /* over print the new and full moons for this month.
X * TODO: don't really know which dates to use here (see moonnf())
X * so try several to be fairly safe. have to go back to 4/29/1988
X * to find the full moon on 5/1 for example.
X */
X mm_nfmoon (jd0-3, tz, m, f);
X mm_nfmoon (jd0+15, tz, m, f);
X}
X
Xstatic
Xmm_nfmoon (jd, tzone, m, f)
Xdouble jd, tzone;
Xint m, f;
X{
X static char nm[] = "NM", fm[] = "FM";
X double dm;
X int mm, ym;
X double jdn, jdf;
X int di;
X
X moonnf (jd, &jdn, &jdf);
X mjd_cal (jdn-tzone/24.0, &mm, &dm, &ym);
X if (m == mm) {
X di = dm + f - 1;
X f_string (R_CAL+2+di/7, C_CAL+3*(di%7), nm);
X }
X mjd_cal (jdf-tzone/24.0, &mm, &dm, &ym);
X if (m == mm) {
X di = dm + f - 1;
X f_string (R_CAL+2+di/7, C_CAL+3*(di%7), fm);
X }
X}
EOFxEOF
len=`wc -c < mainmenu.c`
if expr $len != 6662 > /dev/null
then echo Length of mainmenu.c is $len but it should be 6662.
fi
echo x moon.c
sed -e 's/^X//' << 'EOFxEOF' > moon.c
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X/* given the mjd, find the geocentric ecliptic longitude, lam, and latitude,
X * bet, and horizontal parallax, hp for the moon.
X * N.B. series for long and lat are good to about 10 and 3 arcseconds. however,
X * math errors cause up to 100 and 30 arcseconds error, even if use double.
X * why?? suspect highly sensitive nature of difference used to get m1..6.
X * N.B. still need to correct for nutation. then for topocentric location
X * further correct for parallax and refraction.
X */
Xmoon (mjd, lam, bet, hp)
Xdouble mjd;
Xdouble *lam, *bet, *hp;
X{
X double t, t2;
X double ld;
X double ms;
X double md;
X double de;
X double f;
X double n;
X double a, sa, sn, b, sb, c, sc, e, e2, l, g, w1, w2;
X double m1, m2, m3, m4, m5, m6;
X
X t = mjd/36525.;
X t2 = t*t;
X
X m1 = mjd/27.32158213;
X m1 = 360.0*(m1-(long)m1);
X m2 = mjd/365.2596407;
X m2 = 360.0*(m2-(long)m2);
X m3 = mjd/27.55455094;
X m3 = 360.0*(m3-(long)m3);
X m4 = mjd/29.53058868;
X m4 = 360.0*(m4-(long)m4);
X m5 = mjd/27.21222039;
X m5 = 360.0*(m5-(long)m5);
X m6 = mjd/6798.363307;
X m6 = 360.0*(m6-(long)m6);
X
X ld = 270.434164+m1-(.001133-.0000019*t)*t2;
X ms = 358.475833+m2-(.00015+.0000033*t)*t2;
X md = 296.104608+m3+(.009192+.0000144*t)*t2;
X de = 350.737486+m4-(.001436-.0000019*t)*t2;
X f = 11.250889+m5-(.003211+.0000003*t)*t2;
X n = 259.183275-m6+(.002078+.000022*t)*t2;
X
X a = degrad(51.2+20.2*t);
X sa = sin(a);
X sn = sin(degrad(n));
X b = 346.56+(132.87-.0091731*t)*t;
X sb = .003964*sin(degrad(b));
X c = degrad(n+275.05-2.3*t);
X sc = sin(c);
X ld = ld+.000233*sa+sb+.001964*sn;
X ms = ms-.001778*sa;
X md = md+.000817*sa+sb+.002541*sn;
X f = f+sb-.024691*sn-.004328*sc;
X de = de+.002011*sa+sb+.001964*sn;
X e = 1-(.002495+7.52e-06*t)*t;
X e2 = e*e;
X
X ld = degrad(ld);
X ms = degrad(ms);
X n = degrad(n);
X de = degrad(de);
X f = degrad(f);
X md = degrad(md);
X
X l = 6.28875*sin(md)+1.27402*sin(2*de-md)+.658309*sin(2*de)+
X .213616*sin(2*md)-e*.185596*sin(ms)-.114336*sin(2*f)+
X .058793*sin(2*(de-md))+.057212*e*sin(2*de-ms-md)+
X .05332*sin(2*de+md)+.045874*e*sin(2*de-ms)+.041024*e*sin(md-ms);
X l = l-.034718*sin(de)-e*.030465*sin(ms+md)+.015326*sin(2*(de-f))-
X .012528*sin(2*f+md)-.01098*sin(2*f-md)+.010674*sin(4*de-md)+
X .010034*sin(3*md)+.008548*sin(4*de-2*md)-e*.00791*sin(ms-md+2*de)-
X e*.006783*sin(2*de+ms);
X l = l+.005162*sin(md-de)+e*.005*sin(ms+de)+.003862*sin(4*de)+
X e*.004049*sin(md-ms+2*de)+.003996*sin(2*(md+de))+
X .003665*sin(2*de-3*md)+e*.002695*sin(2*md-ms)+
X .002602*sin(md-2*(f+de))+e*.002396*sin(2*(de-md)-ms)-
X .002349*sin(md+de);
X l = l+e2*.002249*sin(2*(de-ms))-e*.002125*sin(2*md+ms)-
X e2*.002079*sin(2*ms)+e2*.002059*sin(2*(de-ms)-md)-
X .001773*sin(md+2*(de-f))-.001595*sin(2*(f+de))+
X e*.00122*sin(4*de-ms-md)-.00111*sin(2*(md+f))+.000892*sin(md-3*de);
X l = l-e*.000811*sin(ms+md+2*de)+e*.000761*sin(4*de-ms-2*md)+
X e2*.000704*sin(md-2*(ms+de))+e*.000693*sin(ms-2*(md-de))+
X e*.000598*sin(2*(de-f)-ms)+.00055*sin(md+4*de)+.000538*sin(4*md)+
X e*.000521*sin(4*de-ms)+.000486*sin(2*md-de);
X l = l+e2*.000717*sin(md-2*ms);
X *lam = ld+degrad(l);
X range (lam, 2*PI);
X
X g = 5.12819*sin(f)+.280606*sin(md+f)+.277693*sin(md-f)+
X .173238*sin(2*de-f)+.055413*sin(2*de+f-md)+.046272*sin(2*de-f-md)+
X .032573*sin(2*de+f)+.017198*sin(2*md+f)+.009267*sin(2*de+md-f)+
X .008823*sin(2*md-f)+e*.008247*sin(2*de-ms-f);
X g = g+.004323*sin(2*(de-md)-f)+.0042*sin(2*de+f+md)+
X e*.003372*sin(f-ms-2*de)+e*.002472*sin(2*de+f-ms-md)+
X e*.002222*sin(2*de+f-ms)+e*.002072*sin(2*de-f-ms-md)+
X e*.001877*sin(f-ms+md)+.001828*sin(4*de-f-md)-e*.001803*sin(f+ms)-
X .00175*sin(3*f);
X g = g+e*.00157*sin(md-ms-f)-.001487*sin(f+de)-e*.001481*sin(f+ms+md)+
X e*.001417*sin(f-ms-md)+e*.00135*sin(f-ms)+.00133*sin(f-de)+
X .001106*sin(f+3*md)+.00102*sin(4*de-f)+.000833*sin(f+4*de-md)+
X .000781*sin(md-3*f)+.00067*sin(f+4*de-2*md);
X g = g+.000606*sin(2*de-3*f)+.000597*sin(2*(de+md)-f)+
X e*.000492*sin(2*de+md-ms-f)+.00045*sin(2*(md-de)-f)+
X .000439*sin(3*md-f)+.000423*sin(f+2*(de+md))+
X .000422*sin(2*de-f-3*md)-e*.000367*sin(ms+f+2*de-md)-
X e*.000353*sin(ms+f+2*de)+.000331*sin(f+4*de);
X g = g+e*.000317*sin(2*de+f-ms+md)+e2*.000306*sin(2*(de-ms)-f)-
X .000283*sin(md+3*f);
X w1 = .0004664*cos(n);
X w2 = .0000754*cos(c);
X *bet = degrad(g)*(1-w1-w2);
X
X *hp = .950724+.051818*cos(md)+.009531*cos(2*de-md)+.007843*cos(2*de)+
X .002824*cos(2*md)+.000857*cos(2*de+md)+e*.000533*cos(2*de-ms)+
X e*.000401*cos(2*de-md-ms)+e*.00032*cos(md-ms)-.000271*cos(de)-
X e*.000264*cos(ms+md)-.000198*cos(2*f-md);
X *hp = *hp+.000173*cos(3*md)+.000167*cos(4*de-md)-e*.000111*cos(ms)+
X .000103*cos(4*de-2*md)-.000084*cos(2*md-2*de)-
X e*.000083*cos(2*de+ms)+.000079*cos(2*de+2*md)+.000072*cos(4*de)+
X e*.000064*cos(2*de-ms+md)-e*.000063*cos(2*de+ms-md)+
X e*.000041*cos(ms+de);
X *hp = *hp+e*.000035*cos(2*md-ms)-.000033*cos(3*md-2*de)-
X .00003*cos(md+de)-.000029*cos(2*(f-de))-e*.000029*cos(2*md+ms)+
X e2*.000026*cos(2*(de-ms))-.000023*cos(2*(f-de)+md)+
X e*.000019*cos(4*de-ms-md);
X *hp = degrad(*hp);
X}
EOFxEOF
len=`wc -c < moon.c`
if expr $len != 5143 > /dev/null
then echo Length of moon.c is $len but it should be 5143.
fi
echo x moonnf.c
sed -e 's/^X//' << 'EOFxEOF' > moonnf.c
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X#define unw(w,z) ((w)-floor((w)/(z))*(z))
X
X/* given a modified Julian date, mjd, return the mjd of the new
X * and full moons about then, mjdn and mjdf.
X * TODO: exactly which ones does it find? eg:
X * 5/28/1988 yields 5/15 and 5/31
X * 5/29 6/14 6/29
X */
Xmoonnf (mjd, mjdn, mjdf)
Xdouble mjd;
Xdouble *mjdn, *mjdf;
X{
X int mo, yr;
X double dy;
X double mjd0;
X double k, tn, tf, t;
X
X mjd_cal (mjd, &mo, &dy, &yr);
X cal_mjd (1, 0., yr, &mjd0);
X k = (yr-1900+((mjd-mjd0)/365))*12.3685;
X k = floor(k+0.5);
X tn = k/1236.85;
X tf = (k+0.5)/1236.85;
X t = tn;
X m (t, k, mjdn);
X t = tf;
X k += 0.5;
X m (t, k, mjdf);
X}
X
Xstatic
Xm (t, k, mjd)
Xdouble t, k;
Xdouble *mjd;
X{
X double t2, a, a1, b, b1, c, ms, mm, f, ddjd;
X
X t2 = t*t;
X a = 29.53*k;
X c = degrad(166.56+(132.87-9.173e-3*t)*t);
X b = 5.8868e-4*k+(1.178e-4-1.55e-7*t)*t2+3.3e-4*sin(c)+7.5933E-1;
X ms = 359.2242+360*unw(k/1.236886e1,1)-(3.33e-5+3.47e-6*t)*t2;
X mm = 306.0253+360*unw(k/9.330851e-1,1)+(1.07306e-2+1.236e-5*t)*t2;
X f = 21.2964+360*unw(k/9.214926e-1,1)-(1.6528e-3+2.39e-6*t)*t2;
X ms = unw(ms,360);
X mm = unw(mm,360);
X f = unw(f,360);
X ms = degrad(ms);
X mm = degrad(mm);
X f = degrad(f);
X ddjd = (1.734e-1-3.93e-4*t)*sin(ms)+2.1e-3*sin(2*ms)
X -4.068e-1*sin(mm)+1.61e-2*sin(2*mm)-4e-4*sin(3*mm)
X +1.04e-2*sin(2*f)-5.1e-3*sin(ms+mm)-7.4e-3*sin(ms-mm)
X +4e-4*sin(2*f+ms)-4e-4*sin(2*f-ms)-6e-4*sin(2*f+mm)
X +1e-3*sin(2*f-mm)+5e-4*sin(ms+2*mm);
X a1 = (long)a;
X b = b+ddjd+(a-a1);
X b1 = (long)b;
X a = a1+b1;
X b = b-b1;
X *mjd = a + b;
X}
EOFxEOF
len=`wc -c < moonnf.c`
if expr $len != 1557 > /dev/null
then echo Length of moonnf.c is $len but it should be 1557.
fi
echo x nutation.c
sed -e 's/^X//' << 'EOFxEOF' > nutation.c
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X/* given the modified JD, mjd, find the nutation in obliquity, *deps, and
X * the nutation in longitude, *dpsi, each in radians.
X */
Xnutation (mjd, deps, dpsi)
Xdouble mjd;
Xdouble *deps, *dpsi;
X{
X static double lastmjd, lastdeps, lastdpsi;
X double ls, ld; /* sun's mean longitude, moon's mean longitude */
X double ms, md; /* sun's mean anomaly, moon's mean anomaly */
X double nm; /* longitude of moon's ascending node */
X double t, t2; /* number of Julian centuries of 36525 days since
X * Jan 0.5 1900.
X */
X double tls, tnm, tld; /* twice above */
X double a, b; /* temps */
X
X if (mjd == lastmjd) {
X *deps = lastdeps;
X *dpsi = lastdpsi;
X return;
X }
X
X t = mjd/36525.;
X t2 = t*t;
X
X a = 100.0021358*t;
X b = 360.*(a-(long)a);
X ls = 279.697+.000303*t2+b;
X
X a = 1336.855231*t;
X b = 360.*(a-(long)a);
X ld = 270.434-.001133*t2+b;
X
X a = 99.99736056000026*t;
X b = 360.*(a-(long)a);
X ms = 358.476-.00015*t2+b;
X
X a = 13255523.59*t;
X b = 360.*(a-(long)a);
X md = 296.105+.009192*t2+b;
X
X a = 5.372616667*t;
X b = 360.*(a-(long)a);
X nm = 259.183+.002078*t2-b;
X
X /* convert to radian forms for use with trig functions.
X */
X tls = 2*degrad(ls);
X nm = degrad(nm);
X tnm = 2*degrad(nm);
X ms = degrad(ms);
X tld = 2*degrad(ld);
X md = degrad(md);
X
X /* find delta psi and eps, in arcseconds.
X */
X lastdpsi = (-17.2327-.01737*t)*sin(nm)+(-1.2729-.00013*t)*sin(tls)
X +.2088*sin(tnm)-.2037*sin(tld)+(.1261-.00031*t)*sin(ms)
X +.0675*sin(md)-(.0497-.00012*t)*sin(tls+ms)
X -.0342*sin(tld-nm)-.0261*sin(tld+md)+.0214*sin(tls-ms)
X -.0149*sin(tls-tld+md)+.0124*sin(tls-nm)+.0114*sin(tld-md);
X lastdeps = (9.21+.00091*t)*cos(nm)+(.5522-.00029*t)*cos(tls)
X -.0904*cos(tnm)+.0884*cos(tld)+.0216*cos(tls+ms)
X +.0183*cos(tld-nm)+.0113*cos(tld+md)-.0093*cos(tls-ms)
X -.0066*cos(tls-nm);
X
X /* convert to radians.
X */
X lastdpsi = degrad(lastdpsi/3600);
X lastdeps = degrad(lastdeps/3600);
X
X lastmjd = mjd;
X *deps = lastdeps;
X *dpsi = lastdpsi;
X}
EOFxEOF
len=`wc -c < nutation.c`
if expr $len != 2010 > /dev/null
then echo Length of nutation.c is $len but it should be 2010.
fi
echo x objx.c
sed -e 's/^X//' << 'EOFxEOF' > objx.c
X/* functions to save the user-definable object.
X * this way, once defined, the object can be quieried for position just like the
X * other bodies, with objx_cir().
X */
X
X#include <math.h>
X#include "astro.h"
X#include "circum.h"
X#include "screen.h"
X
Xstatic int onflag; /* !=0 while object x is active */
X
Xstatic int xtype; /* see flags, below */
X#define FIXED 0 /* just simple ra/dec object */
X#define ELLIPTICAL 1 /* elliptical orbital elements */
X#define PARABOLIC 2 /* parabolic " */
X
X/* the defining info about object x.
X */
Xtypedef struct {
X double f_ra; /* ra, rads, at given epoch */
X double f_dec; /* dec, rads, at given epoch */
X double f_epoch; /* the given epoch, as an mjd */
X} ObjF; /* fixed object */
Xtypedef struct {
X double e_ep; /* epoch of perihelion, as an mjd */
X double e_p; /* orbital period, in (tropical?) years */
X double e_pl; /* perihelion longitude, degrees */
X double e_q; /* eccentricity */
X double e_inc; /* inclination, degrees */
X double e_om; /* longitude of ascending node, degrees */
X double e_sma; /* semi-major axis, AU */
X double e_m1, e_m2; /* magnitude model coefficients */
X} ObjE; /* object in heliocentric elliptical orbit */
Xtypedef struct {
X double p_ep; /* epoch of perihelion, as an mjd */
X double p_inc; /* inclination, rads */
X double p_qp; /* perihelion distance, AU */
X double p_ap; /* argument of perihelion, rads. */
X double p_om; /* longitude of ascending node, rads */
X double p_epoch; /* reference epoch, as an mjd */
X double p_m1, p_m2; /* magnitude model coefficients */
X} ObjP; /* object in heliocentric parabolic trajectory */
X
Xstatic ObjF objf;
Xstatic ObjE obje;
Xstatic ObjP objp;
X
X/* run when Objx is picked from menu.
X * let op define object and turn it on and off.
X */
Xobjx_setup()
X{
X static char *p[4] = {
X "Fixed coordinates", "Elliptical elements", "Parabolic elements"
X };
X int f;
X
X /* start pointing at selection for current object type */
X switch (xtype) {
X case FIXED: f = 0; break;
X case ELLIPTICAL: f = 1; break;
X case PARABOLIC: f = 2; break;
X }
X
X ask:
X p[3] = onflag ? "On" : "Off";
X switch (f = popup (p, f, 4)) {
X case 0: objx_dfixed((char**)0); goto ask;
X case 1: objx_delliptical((char**)0); goto ask;
X case 2: objx_dparabolic((char**)0); goto ask;
X case 3: onflag ^= 1; break;
X }
X}
X
X/* turn "on" or "off" but don't forget facts about object x.
X */
Xobjx_on ()
X{
X onflag = 1;
X}
X
Xobjx_off()
X{
X onflag = 0;
X}
X
X/* return true if object is now on, else 0.
X */
Xobjx_ison()
X{
X return (onflag);
X}
X
X/* fill in info about object x.
X * most arguments and conditions are the same as for plans().
X * only difference is that mag is already apparent, not absolute magnitude.
X * this is called by body_cir() for object x just like plans() is called
X * for the planets.
X */
Xobjx_cir (jd, lpd0, psi0, rp0, rho0, lam, bet, mag)
Xdouble jd; /* mjd now */
Xdouble *lpd0; /* heliocentric longitude, or NOHELIO */
Xdouble *psi0; /* heliocentric latitude, or 0 if *lpd0 set to NOHELIO */
Xdouble *rp0; /* distance from the sun, or 0 */
Xdouble *rho0; /* true distance from the Earth, or 0 */
Xdouble *lam; /* apparent geocentric ecliptic longitude */
Xdouble *bet; /* apparent geocentric ecliptic latitude */
Xdouble *mag; /* APPARENT magnitude, or NOMAG */
X{
X switch (xtype) {
X case FIXED: {
X double xr, xd;
X xr = objf.f_ra;
X xd = objf.f_dec;
X if (objf.f_epoch != jd)
X precess (objf.f_epoch, jd, &xr, &xd);
X eq_ecl (jd, xr, xd, bet, lam);
X
X *lpd0 = NOHELIO;
X *psi0 = *rp0 = *rho0 = 0.0;
X *mag = NOMAG;
X }
X break;
X case PARABOLIC: {
X double inc, ap, om;
X double lpd, psi, rp, rho;
X double dt;
X int pass;
X /* two passes to correct lam and bet for light travel time. */
X dt = 0.0;
X for (pass = 0; pass < 2; pass++) {
X reduce_elements (objp.p_epoch, jd-dt, objp.p_inc, objp.p_ap,
X objp.p_om, &inc, &ap, &om);
X comet (jd-dt, objp.p_ep, inc, ap, objp.p_qp, om,
X &lpd, &psi, &rp, &rho, lam, bet);
X if (pass == 0) {
X *lpd0 = lpd;
X *psi0 = psi;
X *rp0 = rp;
X *rho0 = rho;
X }
X dt = rho*5.775518e-3; /* au to light-days */
X }
X *mag = objp.p_m1 + 5*log10(*rho0) + objp.p_m2*log10(*rp0);
X }
X break;
X case ELLIPTICAL: {
X /* this is basically the same code as pelement() and plans()
X * combined and simplified for the special case of osculating
X * (unperturbed) elements.
X */
X double a0, a1;
X double t, aa, ml, dayl, dt, lg, lsn, rsn;
X double nu, ea;
X double ma, rp, lp, om, lo, slo, clo;
X double inc, psi, spsi, cpsi;
X double y, lpd, rpd, ll, rho, sll, cll;
X int pass;
X
X a0 = obje.e_pl - 360.0/36525.0*100.0*obje.e_ep/obje.e_p;
X a1 = 100.0/obje.e_p;
X
X t = jd/36525.0;
X aa = a1*t;
X ml = a0 + 360.0*(aa - (long)aa);
X range (&ml, 360.0);
X dayl = a1*9.856263e-3;
X
X dt = 0;
X sunpos (jd, &lsn, &rsn);
X lg = lsn + PI;
X
X for (pass = 0; pass < 2; pass++) {
X ma = degrad(ml-obje.e_pl-dt*dayl);
X anomaly (ma, obje.e_q, &nu, &ea);
X rp = obje.e_sma*(1-obje.e_q*obje.e_q)/(1+obje.e_q*cos(nu));
X lp = raddeg(nu)+obje.e_pl;
X lp = degrad(lp);
X om = degrad(obje.e_om);
X lo = lp-om;
X slo = sin(lo);
X clo = cos(lo);
X inc = degrad(obje.e_inc);
X spsi = slo*sin(inc);
X y = slo*cos(inc);
X psi = asin(spsi);
X lpd = atan(y/clo)+om;
X if (clo<0) lpd += PI;
X range (&lpd, 2*PI);
X cpsi = cos(psi);
X rpd = rp*cpsi;
X ll = lpd-lg;
X rho = sqrt(rsn*rsn+rp*rp-2*rsn*rp*cpsi*cos(ll));
X dt = rho*5.775518e-3;
X if (pass == 0) {
X *lpd0 = lpd;
X *psi0 = psi;
X *rp0 = rp;
X *rho0 = rho;
X }
X }
X
X sll = sin(ll);
X cll = cos(ll);
X *lam = atan(rsn*sll/(rpd-rsn*cll))+lpd;
X /* *lam = atan(-1*rpd*sll/(rsn-rpd*cll))+lg+PI; */
X range (lam, 2*PI);
X *bet = atan(rpd*spsi*sin(*lam-lpd)/(cpsi*rsn*sll));
X *mag = obje.e_m1 + 5*log10(*rho0) + obje.e_m2*log10(*rp0);
X }
X break;
X }
X}
X
X/* define objx based on the ephem.cfg line, s.
X * the "OBJX " keyword is already skipped.
X * format: xtype,[other fields, as per corresponding ObjX typedef]
X * N.B. we replace all ',' with '\0' within s IN PLACE.
X * return 0 if ok, else print reason why not with f_msg() and return -1.
X */
Xobjx_define (s)
Xchar *s;
X{
X#define MAXARGS 20
X char *av[MAXARGS]; /* point to each field for easy reference */
X char c;
X int ac;
X
X /* parse into comma separated fields */
X ac = 0;
X av[0] = s;
X do {
X c = *s++;
X if (c == ',' || c == '\0') {
X s[-1] = '\0';
X av[++ac] = s;
X }
X } while (c);
X
X if (ac < 1) {
X f_msg ("No type for OBJX");
X return (-1);
X }
X
X switch (av[0][0]) {
X case 'f':
X if (ac != 4) {
X f_msg ("Need ra,dec,epoch for \"fixed\" OBJX");
X return (-1);
X }
X objx_dfixed (av+1);
X break;
X case 'p':
X if (ac != 9) {
X f_msg ("Need ep,inc,ap,qp,om,epoch,abs,lum for \"parabolic\" OBJX");
X return (-1);
X }
X objx_dparabolic (av+1);
X break;
X case 'e':
X if (ac != 10) {
X f_msg ("Need ep,p,pl,ecc,inc,om,sma,abs,lum for \"elliptical\" OBJX");
X return (-1);
X }
X objx_delliptical (av+1);
X break;
X default:
X f_msg ("Unknown OBJX type");
X return (-1);
X }
X
X return (0);
X}
X
X/* define a fixed object.
X * args in av, in order, are ra, dec, end reference epoch.
X * if av then it is a list of strings to use for each parameter, else must
X * ask for each. the av option is for cracking the ephem.cfg line.
X * if asking show current settings and leave unchanged if hit RETURN.
X * END aborts without making any more changes.
X * xtype is set to FIXED.
X * N.B. we don't error check av in any way, not even for length.
X */
Xstatic
Xobjx_dfixed (av)
Xchar *av[];
X{
X char buf[NC];
X char *bp;
X int sts;
X
X xtype = FIXED;
X
X if (av) {
X bp = av[0];
X sts = 1;
X } else {
X static char p[] = "RA (h:m:s): (";
X f_prompt (p);
X f_ra (R_PROMPT, C_PROMPT+sizeof(p)-1, objf.f_ra);
X printf (") ");
X sts = read_line (buf, 8+1);
X if (sts < 0)
X return;
X bp = buf;
X }
X if (sts > 0) {
X int h, m, s;
X f_dec_sexsign (radhr(objf.f_ra), &h, &m, &s);
X f_sscansex (bp, &h, &m, &s);
X sex_dec (h, m, s, &objf.f_ra);
X objf.f_ra = hrrad(objf.f_ra);
X }
X
X if (av) {
X bp = av[1];
X sts = 1;
X } else {
X static char p[] = "Dec (d:m:s): (";
X f_prompt (p);
X f_angle (R_PROMPT, C_PROMPT+sizeof(p)-1, objf.f_dec);
X printf (") ");
X sts = read_line (buf, 9+1);
X if (sts < 0)
X return;
X bp = buf;
X }
X if (sts > 0) {
X int dg, m, s;
X f_dec_sexsign (raddeg(objf.f_dec), &dg, &m, &s);
X f_sscansex (bp, &dg, &m, &s);
X sex_dec (dg, m, s, &objf.f_dec);
X objf.f_dec = degrad(objf.f_dec);
X }
X
X if (av) {
X bp = av[2];
X sts = 1;
X } else {
X static char p[] = "Reference epoch (UT Date, m/d.d/y or year.d): ";
X double y;
X f_prompt (p);
X mjd_year (objf.f_epoch, &y);
X printf ("(%g) ", y);
X sts = read_line (buf, 8+1);
X if (sts < 0)
X return;
X bp = buf;
X }
X if (sts > 0)
X crack_year (bp, &objf.f_epoch);
X}
X
X/* define an object in an elliptical heliocentric orbit.
X * args in av, in order, are epoch of perihelion, orbital period, perihelion
X * longitude, eccentricity, inclination, longitude of ascending node, and
X * semi-major axis.
X * if av then it is a list of strings to use for each parameter, else must
X * ask for each. the av option is for cracking the ephem.cfg line.
X * if asking show current settings and leave unchanged if hit RETURN.
X * END aborts without making any more changes.
X * xtype is set to ELLIPTICAL.
X * N.B. we don't error check av in any way, not even for length.
X */
Xstatic
Xobjx_delliptical(av)
Xchar *av[];
X{
X char buf[NC];
X char *bp;
X int sts;
X
X xtype = ELLIPTICAL;
X
X if (av) {
X bp = av[0];
X sts = 1;
X } else {
X static char p[] =
X "Epoch of perihelion (UT Date, m/d.d/y or year.d): ";
X int m, y;
X double d;
X f_prompt (p);
X mjd_cal (obje.e_ep, &m, &d, &y);
X printf ("(%d/%g/%d) ", m, d, y);
X sts = read_line (buf, 15);
X if (sts < 0)
X return;
X bp = buf;
X }
X if (sts > 0)
X crack_year (bp, &obje.e_ep);
X
X if (av) {
X bp = av[1];
X sts = 1;
X } else {
X static char p[] = "Oribital period (years):";
X f_prompt (p);
X f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", obje.e_p);
X sts = read_line (buf, 8+1);
X if (sts < 0)
X return;
X bp = buf;
X }
X if (sts > 0)
X obje.e_p = atof(bp);
X
X if (av) {
X bp = av[2];
X sts = 1;
X } else {
X static char p[] = "Perihelion longitude (degs):";
X f_prompt (p);
X f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", obje.e_pl);
X sts = read_line (buf, 8+1);
X if (sts < 0)
X return;
X bp = buf;
X }
X if (sts > 0)
X obje.e_pl = atof(bp);
X
X if (av) {
X bp = av[3];
X sts = 1;
X } else {
X static char p[] = "Eccentricity:";
X f_prompt (p);
X f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", obje.e_q);
X sts = read_line (buf, 8+1);
X if (sts < 0)
X return;
X bp = buf;
X }
X if (sts > 0)
X obje.e_q = atof(bp);
X
X if (av) {
X bp = av[4];
X sts = 1;
X } else {
X static char p[] = "Inclination (degs):";
X f_prompt (p);
X f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", obje.e_inc);
X sts = read_line (buf, 8+1);
X if (sts < 0)
X return;
X bp = buf;
X }
X if (sts > 0)
X obje.e_inc = atof(bp);
X
X if (av) {
X bp = av[5];
X sts = 1;
X } else {
X static char p[] = "Longitude of ascending node (degs):";
X f_prompt (p);
X f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", obje.e_om);
X sts = read_line (buf, 8+1);
X if (sts < 0)
X return;
X bp = buf;
X }
X if (sts > 0)
X obje.e_om = atof(bp);
X
X if (av) {
X bp = av[6];
X sts = 1;
X } else {
X static char p[] = "Semi-major axis (AU):";
X f_prompt (p);
X f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", obje.e_sma);
X sts = read_line (buf, 8+1);
X if (sts < 0)
X return;
X bp = buf;
X }
X if (sts > 0)
X obje.e_sma = atof(bp);
X
X if (av) {
X bp = av[7];
X sts = 1;
X } else {
X static char p[] = "Absolute magnitude:";
X f_prompt (p);
X f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", obje.e_m1);
X sts = read_line (buf, 8+1);
X if (sts < 0)
X return;
X bp = buf;
X }
X if (sts > 0)
X obje.e_m1 = atof(bp);
X
X if (av) {
X bp = av[8];
X sts = 1;
X } else {
X static char p[] = "Luminosity index:";
X f_prompt (p);
X f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", obje.e_m2);
X sts = read_line (buf, 8+1);
X if (sts < 0)
X return;
X bp = buf;
X }
X if (sts > 0)
X obje.e_m2 = atof(bp);
X}
X
X/* define an object in heliocentric parabolic orbit.
X * args in av, in order, are epoch of perihelion, inclination, argument of
X * perihelion, perihelion distance, longitude of ascending node and
X * reference epoch.
X * if av then it is a list of strings to use for each parameter, else must
X * ask for each. the av option is for cracking the ephem.cfg line.
X * if asking show current settings and leave unchanged if hit RETURN.
X * END aborts without making any more changes.
X * xtype is set to PARABOLIC.
X * N.B. we don't error check av in any way, not even for length.
X */
Xstatic
Xobjx_dparabolic(av)
Xchar *av[];
X{
X char buf[NC];
X char *bp;
X int sts;
X
X xtype = PARABOLIC;
X
X if (av) {
X bp = av[0];
X sts = 1;
X } else {
X static char p[] =
X "Epoch of perihelion (UT Date, m/d.d/y or year.d): ";
X int m, y;
X double d;
X f_prompt (p);
X mjd_cal (objp.p_ep, &m, &d, &y);
X printf ("(%d/%g/%d) ", m, d, y);
X sts = read_line (buf, PW-sizeof(p));
X if (sts < 0)
X return;
X bp = buf;
X }
X if (sts > 0)
X crack_year (bp, &objp.p_ep);
X
X if (av) {
X bp = av[1];
X sts = 1;
X } else {
X static char p[] = "Inclination (degs):";
X f_prompt (p);
X f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ",raddeg(objp.p_inc));
X sts = read_line (buf, 8+1);
X if (sts < 0)
X return;
X bp = buf;
X }
X if (sts > 0)
X objp.p_inc = degrad(atof(bp));
X
X if (av) {
X bp = av[2];
X sts = 1;
X } else {
X static char p[] = "Argument of perihelion (degs):";
X f_prompt (p);
X f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", raddeg(objp.p_ap));
X sts = read_line (buf, 8+1);
X if (sts < 0)
X return;
X bp = buf;
X }
X if (sts > 0)
X objp.p_ap = degrad(atof(bp));
X
X if (av) {
X bp = av[3];
X sts = 1;
X } else {
X static char p[] = "Perihelion distance (AU):";
X f_prompt (p);
X f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", objp.p_qp);
X sts = read_line (buf, 8+1);
X if (sts < 0)
X return;
X bp = buf;
X }
X if (sts > 0)
X objp.p_qp = atof(bp);
X
X if (av) {
X bp = av[4];
X sts = 1;
X } else {
X static char p[] = "Longitude of ascending node (degs):";
X f_prompt (p);
X f_double(R_PROMPT,C_PROMPT+sizeof(p), "(%g) ", raddeg(objp.p_om));
X sts = read_line (buf, 8+1);
X if (sts < 0)
X return;
X bp = buf;
X }
X if (sts > 0)
X objp.p_om = degrad(atof(bp));
X
X if (av) {
X bp = av[5];
X sts = 1;
X } else {
X static char p[] = "Reference epoch (UT Date, m/d.d/y or year.d): ";
X double y;
X f_prompt (p);
X mjd_year (objp.p_epoch, &y);
X printf ("(%g) ", y);
X sts = read_line (buf, 8+1);
X if (sts < 0)
X return;
X bp = buf;
X }
X if (sts > 0)
X crack_year (bp, &objp.p_epoch);
X
X if (av) {
X bp = av[6];
X sts = 1;
X } else {
X static char p[] = "Absolute magnitude:";
X f_prompt (p);
X f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", objp.p_m1);
X sts = read_line (buf, 8+1);
X if (sts < 0)
X return;
X bp = buf;
X }
X if (sts > 0)
X objp.p_m1 = atof(bp);
X
X if (av) {
X bp = av[7];
X sts = 1;
X } else {
X static char p[] = "Luminosity index:";
X f_prompt (p);
X f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", objp.p_m2);
X sts = read_line (buf, 8+1);
X if (sts < 0)
X return;
X bp = buf;
X }
X if (sts > 0)
X objp.p_m2 = atof(bp);
X}
X
X/* given either a decimal year (xxxx. something) or a calendar (x/x/x)
X * convert it to an mjd and store it at *p;
X */
Xstatic
Xcrack_year (bp, p)
Xchar *bp;
Xdouble *p;
X{
X if (decimal_year(bp)) {
X double y = atof (bp);
X year_mjd (y, p);
X } else {
X int m, y;
X double d;
X mjd_cal (objp.p_ep, &m, &d, &y); /* init with current */
X f_sscandate (bp, &m, &d, &y);
X cal_mjd (m, d, y, p);
X }
X}
EOFxEOF
len=`wc -c < objx.c`
if expr $len != 16243 > /dev/null
then echo Length of objx.c is $len but it should be 16243.
fi
echo x obliq.c
sed -e 's/^X//' << 'EOFxEOF' > obliq.c
X#include <stdio.h>
X#include "astro.h"
X
X/* given the modified Julian date, mjd, find the obliquity of the
X * ecliptic, *eps, in radians.
X */
Xobliquity (mjd, eps)
Xdouble mjd;
Xdouble *eps;
X{
X static double lastmjd, lasteps;
X
X if (mjd != lastmjd) {
X double t;
X t = mjd/36525.;
X lasteps = degrad(2.345229444E1 - ((((-1.81E-3*t)+5.9E-3)*t+4.6845E1)*t)/3600.0);
X lastmjd = mjd;
X }
X *eps = lasteps;
X}
EOFxEOF
len=`wc -c < obliq.c`
if expr $len != 409 > /dev/null
then echo Length of obliq.c is $len but it should be 409.
fi
echo x parallax.c
sed -e 's/^X//' << 'EOFxEOF' > parallax.c
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X/* given true ha and dec, tha and tdec, the geographical latitude, phi, the
X * height above sea-level (as a fraction of the earths radius, 6378.16km),
X * ht, and the equatorial horizontal parallax, ehp, find the apparent
X * ha and dec, aha and adec allowing for parallax.
X * all angles in radians. ehp is the angle subtended at the body by the
X * earth's equator.
X */
Xta_par (tha, tdec, phi, ht, ehp, aha, adec)
Xdouble tha, tdec, phi, ht, ehp;
Xdouble *aha, *adec;
X{
X static double last_phi, last_ht, rsp, rcp;
X double rp; /* distance to object in Earth radii */
X double ctha;
X double stdec, ctdec;
X double tdtha, dtha;
X double caha;
X
X /* avoid calcs involving the same phi and ht */
X if (phi != last_phi || ht != last_ht) {
X double cphi, sphi, u;
X cphi = cos(phi);
X sphi = sin(phi);
X u = atan(9.96647e-1*sphi/cphi);
X rsp = (9.96647e-1*sin(u))+(ht*sphi);
X rcp = cos(u)+(ht*cphi);
X last_phi = phi;
X last_ht = ht;
X }
X
X rp = 1/sin(ehp);
X
X ctha = cos(tha);
X stdec = sin(tdec);
X ctdec = cos(tdec);
X tdtha = (rcp*sin(tha))/((rp*ctdec)-(rcp*ctha));
X dtha = atan(tdtha);
X *aha = tha+dtha;
X caha = cos(*aha);
X range (aha, 2*PI);
X *adec = atan(caha*(rp*stdec-rsp)/(rp*ctdec*ctha-rcp));
X}
X
X/* given the apparent ha and dec, aha and adec, the geographical latitude, phi,
X * the height above sea-level (as a fraction of the earths radius, 6378.16km),
X * ht, and the equatorial horizontal parallax, ehp, find the true ha and dec,
X * tha and tdec allowing for parallax.
X * all angles in radians. ehp is the angle subtended at the body by the
X * earth's equator.
X * uses ta_par() iteratively: find a set of true ha/dec that converts back
X * to the given apparent ha/dec.
X */
Xat_par (aha, adec, phi, ht, ehp, tha, tdec)
Xdouble aha, adec, phi, ht, ehp;
Xdouble *tha, *tdec;
X{
X double nha, ndec; /* ha/dec corres. to current true guesses */
X double eha, edec; /* error in ha/dec */
X
X /* first guess for true is just the apparent */
X *tha = aha;
X *tdec = adec;
X
X while (1) {
X ta_par (*tha, *tdec, phi, ht, ehp, &nha, &ndec);
X eha = aha - nha;
X edec = adec - ndec;
X if (fabs(eha)<1e-6 && fabs(edec)<1e-6)
X break;
X *tha += eha;
X *tdec += edec;
X }
X}
EOFxEOF
len=`wc -c < parallax.c`
if expr $len != 2280 > /dev/null
then echo Length of parallax.c is $len but it should be 2280.
fi
echo x pelement.c
sed -e 's/^X//' << 'EOFxEOF' > pelement.c
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X/* this array contains polynomial coefficients to find the various orbital
X * elements for the mean orbit at any instant in time for each major planet.
X * the first five elements are in the form a0 + a1*t + a2*t**2 + a3*t**3,
X * where t is the number of Julian centuries of 36525 Julian days since 1900
X * Jan 0.5. the last three elements are constants.
X *
X * the orbital element (column) indeces are:
X * [ 0- 3]: coefficients for mean longitude, in degrees;
X * [ 4- 7]: coefficients for longitude of the perihelion, in degrees;
X * [ 8-11]: coefficients for eccentricity;
X * [12-15]: coefficients for inclination, in degrees;
X * [16-19]: coefficients for longitude of the ascending node, in degrees;
X * [20]: semi-major axis, in AU;
X * [21]: angular diameter at 1 AU, in arcsec;
X * [22]: standard visual magnitude, ie, the visual magnitude of the planet
X * when at a distance of 1 AU from both the Sun and the Earth and
X * with zero phase angle.
X *
X * the planent (row) indeces are:
X * [0]: Mercury; [1]: Venus; [2]: Mars; [3]: Jupiter; [4]: Saturn;
X * [5]: Uranus; [6]: Neptune; [7]: Pluto.
X */
X#define NPELE (5*4 + 3) /* 4 coeffs for ea of 5 elems, + 3 constants */
Xstatic double elements[8][NPELE] = {
X
X { /* mercury... */
X
X 178.179078, 415.2057519, 3.011e-4, 0.0,
X 75.899697, 1.5554889, 2.947e-4, 0.0,
X .20561421, 2.046e-5, 3e-8, 0.0,
X 7.002881, 1.8608e-3, -1.83e-5, 0.0,
X 47.145944, 1.1852083, 1.739e-4, 0.0,
X .3870986, 6.74, -0.42
X },
X
X { /* venus... */
X
X 342.767053, 162.5533664, 3.097e-4, 0.0,
X 130.163833, 1.4080361, -9.764e-4, 0.0,
X 6.82069e-3, -4.774e-5, 9.1e-8, 0.0,
X 3.393631, 1.0058e-3, -1e-6, 0.0,
X 75.779647, .89985, 4.1e-4, 0.0,
X .7233316, 16.92, -4.4
X },
X
X { /* mars... */
X
X 293.737334, 53.17137642, 3.107e-4, 0.0,
X 3.34218203e2, 1.8407584, 1.299e-4, -1.19e-6,
X 9.33129e-2, 9.2064e-5, 7.7e-8, 0.0,
X 1.850333, -6.75e-4, 1.26e-5, 0.0,
X 48.786442, .7709917, -1.4e-6, -5.33e-6,
X 1.5236883, 9.36, -1.52
X },
X
X { /* jupiter... */
X
X 238.049257, 8.434172183, 3.347e-4, -1.65e-6,
X 1.2720972e1, 1.6099617, 1.05627e-3, -3.43e-6,
X 4.833475e-2, 1.6418e-4, -4.676e-7, -1.7e-9,
X 1.308736, -5.6961e-3, 3.9e-6, 0.0,
X 99.443414, 1.01053, 3.5222e-4, -8.51e-6,
X 5.202561, 196.74, -9.4
X },
X
X { /* saturn... */
X
X 266.564377, 3.398638567, 3.245e-4, -5.8e-6,
X 9.1098214e1, 1.9584158, 8.2636e-4, 4.61e-6,
X 5.589232e-2, -3.455e-4, -7.28e-7, 7.4e-10,
X 2.492519, -3.9189e-3, -1.549e-5, 4e-8,
X 112.790414, .8731951, -1.5218e-4, -5.31e-6,
X 9.554747, 165.6, -8.88
X },
X
X { /* uranus... */
X
X 244.19747, 1.194065406, 3.16e-4, -6e-7,
X 1.71548692e2, 1.4844328, 2.372e-4, -6.1e-7,
X 4.63444e-2, -2.658e-5, 7.7e-8, 0.0,
X .772464, 6.253e-4, 3.95e-5, 0.0,
X 73.477111, .4986678, 1.3117e-3, 0.0,
X 19.21814, 65.8, -7.19
X },
X
X { /* neptune... */
X
X 84.457994, .6107942056, 3.205e-4, -6e-7,
X 4.6727364e1, 1.4245744, 3.9082e-4, -6.05e-7,
X 8.99704e-3, 6.33e-6, -2e-9, 0.0,
X 1.779242, -9.5436e-3, -9.1e-6, 0.0,
X 130.681389, 1.098935, 2.4987e-4, -4.718e-6,
X 30.10957, 62.2, -6.87
X },
X
X { /* pluto...(osculating 1984 jan 21) */
X
X 95.3113544, .3980332167, 0.0, 0.0,
X 224.017, 0.0, 0.0, 0.0,
X .25515, 0.0, 0.0, 0.0,
X 17.1329, 0.0, 0.0, 0.0,
X 110.191, 0.0, 0.0, 0.0,
X 39.8151, 8.2, -1.0
X }
X};
X
X/* given a modified Julian date, mjd, return the elements for the mean orbit
X * at that instant of all the major planets, together with their
X * mean daily motions in longitude, angular diameter and standard visual
X * magnitude.
X * plan[i][j] contains all the values for all the planets at mjd, such that
X * i = 0..7: mercury, venus, mars, jupiter, saturn, unranus, neptune, pluto;
X * j = 0..8: mean longitude, mean daily motion in longitude, longitude of
X * the perihelion, eccentricity, inclination, longitude of the ascending
X * node, length of the semi-major axis, angular diameter from 1 AU, and
X * the standard visual magnitude (see elements[][] comment, above).
X */
Xpelement (mjd, plan)
Xdouble mjd;
Xdouble plan[8][9];
X{
X register double *ep, *pp;
X register double t = mjd/36525.;
X double aa;
X int planet, i;
X
X for (planet = 0; planet < 8; planet++) {
X ep = elements[planet];
X pp = plan[planet];
X aa = ep[1]*t;
X pp[0] = ep[0] + 360.*(aa-(long)aa) + (ep[3]*t + ep[2])*t*t;
X range (pp, 360.);
X pp[1] = (ep[1]*9.856263e-3) + (ep[2] + ep[3])/36525;
X
X for (i = 4; i < 20; i += 4)
X pp[i/4+1] = ((ep[i+3]*t + ep[i+2])*t + ep[i+1])*t + ep[i+0];
X
X pp[6] = ep[20];
X pp[7] = ep[21];
X pp[8] = ep[22];
X }
X}
EOFxEOF
len=`wc -c < pelement.c`
if expr $len != 4797 > /dev/null
then echo Length of pelement.c is $len but it should be 4797.
fi
echo x plans.c
sed -e 's/^X//' << 'EOFxEOF' > plans.c
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X#define TWOPI (2*PI)
X#define mod2PI(x) ((x) - (long)((x)/TWOPI)*TWOPI)
X
X/* given a modified Julian date, mjd, and a planet, p, find:
X * lpd0: heliocentric longitude,
X * psi0: heliocentric latitude,
X * rp0: distance from the sun to the planet,
X * rho0: distance from the Earth to the planet,
X * none corrected for light time, ie, they are the true values for the
X * given instant.
X * lam: geocentric ecliptic longitude,
X * bet: geocentric ecliptic latitude,
X * each corrected for light time, ie, they are the apparent values as
X * seen from the center of the Earth for the given instant.
X * dia: angular diameter in arcsec at 1 AU,
X * mag: visual magnitude when 1 AU from sun and earth at 0 phase angle.
X *
X * all angles are in radians, all distances in AU.
X * the mean orbital elements are found by calling pelement(), then mutual
X * perturbation corrections are applied as necessary.
X *
X * corrections for nutation and abberation must be made by the caller. The RA
X * and DEC calculated from the fully-corrected ecliptic coordinates are then
X * the apparent geocentric coordinates. Further corrections can be made, if
X * required, for atmospheric refraction and geocentric parallax although the
X * intrinsic error herein of about 10 arcseconds is usually the dominant
X * error at this stage.
X * TODO: combine the several intermediate expressions when get a good compiler.
X */
Xplans (mjd, p, lpd0, psi0, rp0, rho0, lam, bet, dia, mag)
Xdouble mjd;
Xint p;
Xdouble *lpd0, *psi0, *rp0, *rho0, *lam, *bet, *dia, *mag;
X{
X static double plan[8][9];
X static double lastmjd = -10000;
X double dl; /* perturbation correction for longitude */
X double dr; /* " orbital radius */
X double dml; /* " mean longitude */
X double ds; /* " eccentricity */
X double dm; /* " mean anomaly */
X double da; /* " semi-major axis */
X double dhl; /* " heliocentric longitude */
X double lsn, rsn; /* true geocentric longitude of sun and sun-earth rad */
X double mas; /* mean anomaly of the sun */
X double re; /* radius of earth's orbit */
X double lg; /* longitude of earth */
X double map[8]; /* array of mean anomalies for each planet */
X double lpd, psi, rp, rho;
X double ll, sll, cll;
X double t;
X double dt;
X int pass;
X int j;
X double s, ma;
X double nu, ea;
X double lp, om;
X double lo, slo, clo;
X double inc, y;
X double spsi, cpsi;
X double rpd;
X
X /* only need to fill in plan[] once for a given mjd */
X if (mjd != lastmjd) {
X pelement (mjd, plan);
X lastmjd = mjd;
X }
X
X dt = 0;
X t = mjd/36525.;
X sunpos (mjd, &lsn, &rsn);
X masun (mjd, &mas);
X re = rsn;
X lg = lsn+PI;
X
X /* first find the true position of the planet at mjd.
X * then repeat a second time for a slightly different time based
X * on the position found in the first pass to account for light-travel
X * time.
X */
X for (pass = 0; pass < 2; pass++) {
X
X for (j = 0; j < 8; j++)
X map[j] = degrad(plan[j][0]-plan[j][2]-dt*plan[j][1]);
X
X /* set initial corrections to 0.
X * then modify as necessary for the planet of interest.
X */
X dl = 0;
X dr = 0;
X dml = 0;
X ds = 0;
X dm = 0;
X da = 0;
X dhl = 0;
X
X switch (p) {
X
X case MERCURY:
X p_mercury (map, &dl, &dr);
X break;
X
X case VENUS:
X p_venus (t, mas, map, &dl, &dr, &dml, &dm);
X break;
X
X case MARS:
X p_mars (mas, map, &dl, &dr, &dml, &dm);
X break;
X
X case JUPITER:
X p_jupiter (t, plan[p][3], &dml, &ds, &dm, &da);
X break;
X
X case SATURN:
X p_saturn (t, plan[p][3], &dml, &ds, &dm, &da, &dhl);
X break;
X
X case URANUS:
X p_uranus (t, plan[p][3], &dl, &dr, &dml, &ds, &dm, &da, &dhl);
X break;
X
X case NEPTUNE:
X p_neptune (t, plan[p][3], &dl, &dr, &dml, &ds, &dm, &da, &dhl);
X break;
X
X case PLUTO:
X /* no perturbation theory for pluto */
X break;
X }
X
X s = plan[p][3]+ds;
X ma = map[p]+dm;
X anomaly (ma, s, &nu, &ea);
X rp = (plan[p][6]+da)*(1-s*s)/(1+s*cos(nu));
X lp = raddeg(nu)+plan[p][2]+raddeg(dml-dm);
X lp = degrad(lp);
X om = degrad(plan[p][5]);
X lo = lp-om;
X slo = sin(lo);
X clo = cos(lo);
X inc = degrad(plan[p][4]);
X rp = rp+dr;
X spsi = slo*sin(inc);
X y = slo*cos(inc);
X psi = asin(spsi)+dhl;
X spsi = sin(psi);
X lpd = atan(y/clo)+om+degrad(dl);
X if (clo<0) lpd += PI;
X range (&lpd, TWOPI);
X cpsi = cos(psi);
X rpd = rp*cpsi;
X ll = lpd-lg;
X rho = sqrt(re*re+rp*rp-2*re*rp*cpsi*cos(ll));
X
X /* when we view a planet we see it in the position it occupied
X * dt days ago, where rho is the distance between it and earth,
X * in AU. use this as the new time for the next pass.
X */
X dt = rho*5.775518e-3;
X
X if (pass == 0) {
X /* save heliocentric coordinates after first pass since, being
X * true, they are NOT to be corrected for light-travel time.
X */
X *lpd0 = lpd;
X range (lpd0, TWOPI);
X *psi0 = psi;
X *rp0 = rp;
X *rho0 = rho;
X }
X }
X
X sll = sin(ll);
X cll = cos(ll);
X if (p < MARS)
X *lam = atan(-1*rpd*sll/(re-rpd*cll))+lg+PI;
X else
X *lam = atan(re*sll/(rpd-re*cll))+lpd;
X range (lam, TWOPI);
X *bet = atan(rpd*spsi*sin(*lam-lpd)/(cpsi*re*sll));
X *dia = plan[p][7];
X *mag = plan[p][8];
X}
X
X/* set auxilliary variables used for jupiter, saturn, uranus, and neptune */
Xstatic
Xaux_jsun (t, x1, x2, x3, x4, x5, x6)
Xdouble t;
Xdouble *x1, *x2, *x3, *x4, *x5, *x6;
X{
X *x1 = t/5+0.1;
X *x2 = mod2PI(4.14473+5.29691e1*t);
X *x3 = mod2PI(4.641118+2.132991e1*t);
X *x4 = mod2PI(4.250177+7.478172*t);
X *x5 = 5 * *x3 - 2 * *x2;
X *x6 = 2 * *x2 - 6 * *x3 + 3 * *x4;
X}
X
X/* find the mean anomaly of the sun at mjd.
X * this is the same as that used in sun() but when it was converted to C it
X * was not known it would be required outside that routine.
X * TODO: add an argument to sun() to return mas and eliminate this routine.
X */
Xstatic
Xmasun (mjd, mas)
Xdouble mjd;
Xdouble *mas;
X{
X double t, t2;
X double a, b;
X
X t = mjd/36525;
X t2 = t*t;
X a = 9.999736042e1*t;
X b = 360.*(a-(long)a);
X *mas = degrad (3.5847583e2-(1.5e-4+3.3e-6*t)*t2+b);
X}
X
X/* perturbations for mercury */
Xstatic
Xp_mercury (map, dl, dr)
Xdouble map[];
Xdouble *dl, *dr;
X{
X *dl = 2.04e-3*cos(5*map[2-1]-2*map[1-1]+2.1328e-1)+
X 1.03e-3*cos(2*map[2-1]-map[1-1]-2.8046)+
X 9.1e-4*cos(2*map[3]-map[1-1]-6.4582e-1)+
X 7.8e-4*cos(5*map[2-1]-3*map[1-1]+1.7692e-1);
X
X *dr = 7.525e-6*cos(2*map[3]-map[1-1]+9.25251e-1)+
X 6.802e-6*cos(5*map[2-1]-3*map[1-1]-4.53642)+
X 5.457e-6*cos(2*map[2-1]-2*map[1-1]-1.24246)+
X 3.569e-6*cos(5*map[2-1]-map[1-1]-1.35699);
X}
X
X/* ....venus */
Xstatic
Xp_venus (t, mas, map, dl, dr, dml, dm)
Xdouble t, mas, map[];
Xdouble *dl, *dr, *dml, *dm;
X{
X *dml = degrad (7.7e-4*sin(4.1406+t*2.6227));
X *dm = *dml;
X
X *dl = 3.13e-3*cos(2*mas-2*map[2-1]-2.587)+
X 1.98e-3*cos(3*mas-3*map[2-1]+4.4768e-2)+
X 1.36e-3*cos(mas-map[2-1]-2.0788)+
X 9.6e-4*cos(3*mas-2*map[2-1]-2.3721)+
X 8.2e-4*cos(map[3]-map[2-1]-3.6318);
X
X *dr = 2.2501e-5*cos(2*mas-2*map[2-1]-1.01592)+
X 1.9045e-5*cos(3*mas-3*map[2-1]+1.61577)+
X 6.887e-6*cos(map[3]-map[2-1]-2.06106)+
X 5.172e-6*cos(mas-map[2-1]-5.08065e-1)+
X 3.62e-6*cos(5*mas-4*map[2-1]-1.81877)+
X 3.283e-6*cos(4*mas-4*map[2-1]+1.10851)+
X 3.074e-6*cos(2*map[3]-2*map[2-1]-9.62846e-1);
X}
X
X/* ....mars */
Xstatic
Xp_mars (mas, map, dl, dr, dml, dm)
Xdouble mas, map[];
Xdouble *dl, *dr, *dml, *dm;
X{
X double a;
X
X a = 3*map[3]-8*map[2]+4*mas;
X *dml = degrad (-1*(1.133e-2*sin(a)+9.33e-3*cos(a)));
X *dm = *dml;
X
X *dl = 7.05e-3*cos(map[3]-map[2]-8.5448e-1)+
X 6.07e-3*cos(2*map[3]-map[2]-3.2873)+
X 4.45e-3*cos(2*map[3]-2*map[2]-3.3492)+
X 3.88e-3*cos(mas-2*map[2]+3.5771e-1)+
X 2.38e-3*cos(mas-map[2]+6.1256e-1)+
X 2.04e-3*cos(2*mas-3*map[2]+2.7688)+
X 1.77e-3*cos(3*map[2]-map[2-1]-1.0053)+
X 1.36e-3*cos(2*mas-4*map[2]+2.6894)+
X 1.04e-3*cos(map[3]+3.0749e-1);
X
X *dr = 5.3227e-5*cos(map[3]-map[2]+7.17864e-1)+
X 5.0989e-5*cos(2*map[3]-2*map[2]-1.77997)+
X 3.8278e-5*cos(2*map[3]-map[2]-1.71617)+
X 1.5996e-5*cos(mas-map[2]-9.69618e-1)+
X 1.4764e-5*cos(2*mas-3*map[2]+1.19768)+
X 8.966e-6*cos(map[3]-2*map[2]+7.61225e-1);
X *dr += 7.914e-6*cos(3*map[3]-2*map[2]-2.43887)+
X 7.004e-6*cos(2*map[3]-3*map[2]-1.79573)+
X 6.62e-6*cos(mas-2*map[2]+1.97575)+
X 4.93e-6*cos(3*map[3]-3*map[2]-1.33069)+
X 4.693e-6*cos(3*mas-5*map[2]+3.32665)+
X 4.571e-6*cos(2*mas-4*map[2]+4.27086)+
X 4.409e-6*cos(3*map[3]-map[2]-2.02158);
X}
X
X/* ....jupiter */
Xstatic
Xp_jupiter (t, s, dml, ds, dm, da)
Xdouble t, s;
Xdouble *dml, *ds, *dm, *da;
X{
X double dp;
X double x1, x2, x3, x4, x5, x6, x7;
X double sx3, cx3, s2x3, c2x3;
X double sx5, cx5, s2x5;
X double sx6;
X double sx7, cx7, s2x7, c2x7, s3x7, c3x7, s4x7, c4x7, c5x7;
X
X aux_jsun (t, &x1, &x2, &x3, &x4, &x5, &x6);
X x7 = x3-x2;
X sx3 = sin(x3);
X cx3 = cos(x3);
X s2x3 = sin(2*x3);
X c2x3 = cos(2*x3);
X sx5 = sin(x5);
X cx5 = cos(x5);
X s2x5 = sin(2*x5);
X sx6 = sin(x6);
X sx7 = sin(x7);
X cx7 = cos(x7);
X s2x7 = sin(2*x7);
X c2x7 = cos(2*x7);
X s3x7 = sin(3*x7);
X c3x7 = cos(3*x7);
X s4x7 = sin(4*x7);
X c4x7 = cos(4*x7);
X c5x7 = cos(5*x7);
X
X *dml = (3.31364e-1-(1.0281e-2+4.692e-3*x1)*x1)*sx5+
X (3.228e-3-(6.4436e-2-2.075e-3*x1)*x1)*cx5-
X (3.083e-3+(2.75e-4-4.89e-4*x1)*x1)*s2x5+
X 2.472e-3*sx6+1.3619e-2*sx7+1.8472e-2*s2x7+6.717e-3*s3x7+
X 2.775e-3*s4x7+6.417e-3*s2x7*sx3+
X (7.275e-3-1.253e-3*x1)*sx7*sx3+
X 2.439e-3*s3x7*sx3-(3.5681e-2+1.208e-3*x1)*sx7*cx3;
X *dml += -3.767e-3*c2x7*sx3-(3.3839e-2+1.125e-3*x1)*cx7*sx3-
X 4.261e-3*s2x7*cx3+
X (1.161e-3*x1-6.333e-3)*cx7*cx3+
X 2.178e-3*cx3-6.675e-3*c2x7*cx3-2.664e-3*c3x7*cx3-
X 2.572e-3*sx7*s2x3-3.567e-3*s2x7*s2x3+2.094e-3*cx7*c2x3+
X 3.342e-3*c2x7*c2x3;
X *dml = degrad(*dml);
X
X *ds = (3606+(130-43*x1)*x1)*sx5+(1289-580*x1)*cx5-6764*sx7*sx3-
X 1110*s2x7*sx3-224*s3x7*sx3-204*sx3+(1284+116*x1)*cx7*sx3+
X 188*c2x7*sx3+(1460+130*x1)*sx7*cx3+224*s2x7*cx3-817*cx3+
X 6074*cx3*cx7+992*c2x7*cx3+
X 508*c3x7*cx3+230*c4x7*cx3+108*c5x7*cx3;
X *ds += -(956+73*x1)*sx7*s2x3+448*s2x7*s2x3+137*s3x7*s2x3+
X (108*x1-997)*cx7*s2x3+480*c2x7*s2x3+148*c3x7*s2x3+
X (99*x1-956)*sx7*c2x3+490*s2x7*c2x3+
X 158*s3x7*c2x3+179*c2x3+(1024+75*x1)*cx7*c2x3-
X 437*c2x7*c2x3-132*c3x7*c2x3;
X *ds *= 1e-7;
X
X dp = (7.192e-3-3.147e-3*x1)*sx5-4.344e-3*sx3+
X (x1*(1.97e-4*x1-6.75e-4)-2.0428e-2)*cx5+
X 3.4036e-2*cx7*sx3+(7.269e-3+6.72e-4*x1)*sx7*sx3+
X 5.614e-3*c2x7*sx3+2.964e-3*c3x7*sx3+3.7761e-2*sx7*cx3+
X 6.158e-3*s2x7*cx3-
X 6.603e-3*cx7*cx3-5.356e-3*sx7*s2x3+2.722e-3*s2x7*s2x3+
X 4.483e-3*cx7*s2x3-2.642e-3*c2x7*s2x3+4.403e-3*sx7*c2x3-
X 2.536e-3*s2x7*c2x3+5.547e-3*cx7*c2x3-2.689e-3*c2x7*c2x3;
X
X *dm = *dml-(degrad(dp)/s);
X
X *da = 205*cx7-263*cx5+693*c2x7+312*c3x7+147*c4x7+299*sx7*sx3+
X 181*c2x7*sx3+204*s2x7*cx3+111*s3x7*cx3-337*cx7*cx3-
X 111*c2x7*cx3;
X *da *= 1e-6;
X}
X
X/* ....saturn */
Xstatic
Xp_saturn (t, s, dml, ds, dm, da, dhl)
Xdouble t, s;
Xdouble *dml, *ds, *dm, *da, *dhl;
X{
X double dp;
X double x1, x2, x3, x4, x5, x6, x7, x8;
X double sx3, cx3, s2x3, c2x3, s3x3, c3x3, s4x3, c4x3;
X double sx5, cx5, s2x5, c2x5;
X double sx6;
X double sx7, cx7, s2x7, c2x7, s3x7, c3x7, s4x7, c4x7, c5x7, s5x7;
X double s2x8, c2x8, s3x8, c3x8;
X
X aux_jsun (t, &x1, &x2, &x3, &x4, &x5, &x6);
X x7 = x3-x2;
X sx3 = sin(x3);
X cx3 = cos(x3);
X s2x3 = sin(2*x3);
X c2x3 = cos(2*x3);
X sx5 = sin(x5);
X cx5 = cos(x5);
X s2x5 = sin(2*x5);
X sx6 = sin(x6);
X sx7 = sin(x7);
X cx7 = cos(x7);
X s2x7 = sin(2*x7);
X c2x7 = cos(2*x7);
X s3x7 = sin(3*x7);
X c3x7 = cos(3*x7);
X s4x7 = sin(4*x7);
X c4x7 = cos(4*x7);
X c5x7 = cos(5*x7);
X
X s3x3 = sin(3*x3);
X c3x3 = cos(3*x3);
X s4x3 = sin(4*x3);
X c4x3 = cos(4*x3);
X c2x5 = cos(2*x5);
X s5x7 = sin(5*x7);
X x8 = x4-x3;
X s2x8 = sin(2*x8);
X c2x8 = cos(2*x8);
X s3x8 = sin(3*x8);
X c3x8 = cos(3*x8);
X
X *dml = 7.581e-3*s2x5-7.986e-3*sx6-1.48811e-1*sx7-4.0786e-2*s2x7-
X (8.14181e-1-(1.815e-2-1.6714e-2*x1)*x1)*sx5-
X (1.0497e-2-(1.60906e-1-4.1e-3*x1)*x1)*cx5-1.5208e-2*s3x7-
X 6.339e-3*s4x7-6.244e-3*sx3-1.65e-2*s2x7*sx3+
X (8.931e-3+2.728e-3*x1)*sx7*sx3-5.775e-3*s3x7*sx3+
X (8.1344e-2+3.206e-3*x1)*cx7*sx3+1.5019e-2*c2x7*sx3;
X *dml += (8.5581e-2+2.494e-3*x1)*sx7*cx3+1.4394e-2*c2x7*cx3+
X (2.5328e-2-3.117e-3*x1)*cx7*cx3+
X 6.319e-3*c3x7*cx3+6.369e-3*sx7*s2x3+9.156e-3*s2x7*s2x3+
X 7.525e-3*s3x8*s2x3-5.236e-3*cx7*c2x3-7.736e-3*c2x7*c2x3-
X 7.528e-3*c3x8*c2x3;
X *dml = degrad(*dml);
X
X *ds = (-7927+(2548+91*x1)*x1)*sx5+(13381+(1226-253*x1)*x1)*cx5+
X (248-121*x1)*s2x5-(305+91*x1)*c2x5+412*s2x7+12415*sx3+
X (390-617*x1)*sx7*sx3+(165-204*x1)*s2x7*sx3+26599*cx7*sx3-
X 4687*c2x7*sx3-1870*c3x7*sx3-821*c4x7*sx3-
X 377*c5x7*sx3+497*c2x8*sx3+(163-611*x1)*cx3;
X *ds += -12696*sx7*cx3-4200*s2x7*cx3-1503*s3x7*cx3-619*s4x7*cx3-
X 268*s5x7*cx3-(282+1306*x1)*cx7*cx3+(-86+230*x1)*c2x7*cx3+
X 461*s2x8*cx3-350*s2x3+(2211-286*x1)*sx7*s2x3-
X 2208*s2x7*s2x3-568*s3x7*s2x3-346*s4x7*s2x3-
X (2780+222*x1)*cx7*s2x3+(2022+263*x1)*c2x7*s2x3+248*c3x7*s2x3+
X 242*s3x8*s2x3+467*c3x8*s2x3-490*c2x3-(2842+279*x1)*sx7*c2x3;
X *ds += (128+226*x1)*s2x7*c2x3+224*s3x7*c2x3+
X (-1594+282*x1)*cx7*c2x3+(2162-207*x1)*c2x7*c2x3+
X 561*c3x7*c2x3+343*c4x7*c2x3+469*s3x8*c2x3-242*c3x8*c2x3-
X 205*sx7*s3x3+262*s3x7*s3x3+208*cx7*c3x3-271*c3x7*c3x3-
X 382*c3x7*s4x3-376*s3x7*c4x3;
X *ds *= 1e-7;
X
X dp = (7.7108e-2+(7.186e-3-1.533e-3*x1)*x1)*sx5-7.075e-3*sx7+
X (4.5803e-2-(1.4766e-2+5.36e-4*x1)*x1)*cx5-7.2586e-2*cx3-
X 7.5825e-2*sx7*sx3-2.4839e-2*s2x7*sx3-8.631e-3*s3x7*sx3-
X 1.50383e-1*cx7*cx3+2.6897e-2*c2x7*cx3+1.0053e-2*c3x7*cx3-
X (1.3597e-2+1.719e-3*x1)*sx7*s2x3+1.1981e-2*s2x7*c2x3;
X dp += -(7.742e-3-1.517e-3*x1)*cx7*s2x3+
X (1.3586e-2-1.375e-3*x1)*c2x7*c2x3-
X (1.3667e-2-1.239e-3*x1)*sx7*c2x3+
X (1.4861e-2+1.136e-3*x1)*cx7*c2x3-
X (1.3064e-2+1.628e-3*x1)*c2x7*c2x3;
X
X *dm = *dml-(degrad(dp)/s);
X
X *da = 572*sx5-1590*s2x7*cx3+2933*cx5-647*s3x7*cx3+33629*cx7-
X 344*s4x7*cx3-3081*c2x7+2885*cx7*cx3-1423*c3x7+
X (2172+102*x1)*c2x7*cx3-671*c4x7+296*c3x7*cx3-320*c5x7-
X 267*s2x7*s2x3+1098*sx3-778*cx7*s2x3-2812*sx7*sx3;
X *da += 495*c2x7*s2x3+688*s2x7*sx3+250*c3x7*s2x3-393*s3x7*sx3-
X 856*sx7*c2x3-228*s4x7*sx3+441*s2x7*c2x3+2138*cx7*sx3+
X 296*c2x7*c2x3-999*c2x7*sx3+211*c3x7*c2x3-642*c3x7*sx3-
X 427*sx7*s3x3-325*c4x7*sx3+398*s3x7*s3x3-890*cx3+
X 344*cx7*c3x3+2206*sx7*cx3-427*c3x7*c3x3;
X *da *= 1e-6;
X
X *dhl = 7.47e-4*cx7*sx3+1.069e-3*cx7*cx3+2.108e-3*s2x7*s2x3+
X 1.261e-3*c2x7*s2x3+1.236e-3*s2x7*c2x3-2.075e-3*c2x7*c2x3;
X *dhl = degrad(*dhl);
X}
X
X/* ....uranus */
Xstatic
Xp_uranus (t, s, dl, dr, dml, ds, dm, da, dhl)
Xdouble t, s;
Xdouble *dl, *dr, *dml, *ds, *dm, *da, *dhl;
X{
X double dp;
X double x1, x2, x3, x4, x5, x6;
X double x8, x9, x10, x11, x12;
X double sx4, cx4, s2x4, c2x4;
X double sx9, cx9, s2x9, c2x9;
X double sx11, cx11;
X
X aux_jsun (t, &x1, &x2, &x3, &x4, &x5, &x6);
X
X x8 = mod2PI(1.46205+3.81337*t);
X x9 = 2*x8-x4;
X sx9 = sin(x9);
X cx9 = cos(x9);
X s2x9 = sin(2*x9);
X c2x9 = cos(2*x9);
X
X x10 = x4-x2;
X x11 = x4-x3;
X x12 = x8-x4;
X
X *dml = (8.64319e-1-1.583e-3*x1)*sx9+(8.2222e-2-6.833e-3*x1)*cx9+
X 3.6017e-2*s2x9-3.019e-3*c2x9+8.122e-3*sin(x6);
X *dml = degrad(*dml);
X
X dp = 1.20303e-1*sx9+6.197e-3*s2x9+(1.9472e-2-9.47e-4*x1)*cx9;
X *dm = *dml-(degrad(dp)/s);
X
X *ds = (163*x1-3349)*sx9+20981*cx9+1311*c2x9;
X *ds *= 1e-7;
X
X *da = -3.825e-3*cx9;
X
X *dl = (1.0122e-2-9.88e-4*x1)*sin(x4+x11)+
X (-3.8581e-2+(2.031e-3-1.91e-3*x1)*x1)*cos(x4+x11)+
X (3.4964e-2-(1.038e-3-8.68e-4*x1)*x1)*cos(2*x4+x11)+
X 5.594e-3*sin(x4+3*x12)-1.4808e-2*sin(x10)-
X 5.794e-3*sin(x11)+2.347e-3*cos(x11)+9.872e-3*sin(x12)+
X 8.803e-3*sin(2*x12)-4.308e-3*sin(3*x12);
X
X sx11 = sin(x11);
X cx11 = cos(x11);
X sx4 = sin(x4);
X cx4 = cos(x4);
X s2x4 = sin(2*x4);
X c2x4 = cos(2*x4);
X *dhl = (4.58e-4*sx11-6.42e-4*cx11-5.17e-4*cos(4*x12))*sx4-
X (3.47e-4*sx11+8.53e-4*cx11+5.17e-4*sin(4*x11))*cx4+
X 4.03e-4*(cos(2*x12)*s2x4+sin(2*x12)*c2x4);
X *dhl = degrad(*dhl);
X
X *dr = -25948+4985*cos(x10)-1230*cx4+3354*cos(x11)+904*cos(2*x12)+
X 894*(cos(x12)-cos(3*x12))+(5795*cx4-1165*sx4+1388*c2x4)*sx11+
X (1351*cx4+5702*sx4+1388*s2x4)*cos(x11);
X *dr *= 1e-6;
X}
X
X/* ....neptune */
Xstatic
Xp_neptune (t, s, dl, dr, dml, ds, dm, da, dhl)
Xdouble t, s;
Xdouble *dl, *dr, *dml, *ds, *dm, *da, *dhl;
X{
X double dp;
X double x1, x2, x3, x4, x5, x6;
X double x8, x9, x10, x11, x12;
X double sx8, cx8;
X double sx9, cx9, s2x9, c2x9;
X double s2x12, c2x12;
X
X aux_jsun (t, &x1, &x2, &x3, &x4, &x5, &x6);
X
X x8 = mod2PI(1.46205+3.81337*t);
X x9 = 2*x8-x4;
X sx9 = sin(x9);
X cx9 = cos(x9);
X s2x9 = sin(2*x9);
X c2x9 = cos(2*x9);
X
X x10 = x8-x2;
X x11 = x8-x3;
X x12 = x8-x4;
X
X *dml = (1.089e-3*x1-5.89833e-1)*sx9+(4.658e-3*x1-5.6094e-2)*cx9-
X 2.4286e-2*s2x9;
X *dml = degrad(*dml);
X
X dp = 2.4039e-2*sx9-2.5303e-2*cx9+6.206e-3*s2x9-5.992e-3*c2x9;
X
X *dm = *dml-(degrad(dp)/s);
X
X *ds = 4389*sx9+1129*s2x9+4262*cx9+1089*c2x9;
X *ds *= 1e-7;
X
X *da = 8189*cx9-817*sx9+781*c2x9;
X *da *= 1e-6;
X
X s2x12 = sin(2*x12);
X c2x12 = cos(2*x12);
X sx8 = sin(x8);
X cx8 = cos(x8);
X *dl = -9.556e-3*sin(x10)-5.178e-3*sin(x11)+2.572e-3*s2x12-
X 2.972e-3*c2x12*sx8-2.833e-3*s2x12*cx8;
X
X *dhl = 3.36e-4*c2x12*sx8+3.64e-4*s2x12*cx8;
X *dhl = degrad(*dhl);
X
X *dr = -40596+4992*cos(x10)+2744*cos(x11)+2044*cos(x12)+1051*c2x12;
X *dr *= 1e-6;
X}
X
EOFxEOF
len=`wc -c < plans.c`
if expr $len != 17648 > /dev/null
then echo Length of plans.c is $len but it should be 17648.
fi