home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / x / volume19 / xephem / part18 / circum.c next >
Encoding:
C/C++ Source or Header  |  1993-05-15  |  19.5 KB  |  644 lines
  1. /* given a Now and an Obj with the object definition portion filled in,
  2.  * fill in the sky position (s_*) portions.
  3.  */
  4.  
  5. #include <stdio.h>
  6. #include <math.h>
  7. #if defined(__STDC__)
  8. #include <stdlib.h>
  9. #endif
  10. #include "astro.h"
  11. #include "circum.h"
  12. #include "preferences.h"
  13.  
  14. #if defined(__STDC__) || defined(__cplusplus)
  15. #define P_(s) s
  16. #else
  17. #define P_(s) ()
  18. #endif
  19.  
  20. extern void anomaly P_((double ma, double s, double *nu, double *ea));
  21. extern void comet P_((double Mjd, double ep, double inc, double ap, double qp, double om, double *lpd, double *psi, double *rp, double *rho, double *lam, double *bet));
  22. extern void ecl_eq P_((double Mjd, double lati, double lngi, double *ra, double *dec));
  23. extern void eq_ecl P_((double Mjd, double ra, double dec, double *lati, double *lngi));
  24. extern void gk_mag P_((double g, double k, double rp, double rho, double *mp));
  25. extern void hadec_aa P_((double lati, double ha, double dec, double *alt, double *az));
  26. extern void hg_mag P_((double h, double g, double rp, double rho, double rsn, double *mp));
  27. extern void moon P_((double Mjd, double *lam, double *bet, double *hp));
  28. extern void now_lst P_((Now *np, double *lst));
  29. extern void nutation P_((double Mjd, double *deps, double *dpsi));
  30. extern void plans P_((double Mjd, int p, double *lpd0, double *psi0, double *rp0, double *rho0, double *lam, double *bet, double *dia, double *mag));
  31. extern void precess P_((double mjd1, double mjd2, double *ra, double *dec));
  32. extern void range P_((double *v, double r));
  33. extern void reduce_elements P_((double mjd0, double Mjd, double inc0, double ap0, double om0, double *inc, double *ap, double *om));
  34. extern void refract P_((double pr, double tr, double ta, double *aa));
  35. extern void sunpos P_((double Mjd, double *lsn, double *rsn));
  36. extern void ta_par P_((double tha, double tdec, double phi, double ht, double ehp, double *aha, double *adec));
  37.  
  38. int obj_cir P_((Now *np, Obj *op));
  39. static obj_planet P_((Now *np, Obj *op));
  40. static obj_fixed P_((Now *np, Obj *op));
  41. static obj_elliptical P_((Now *np, Obj *op));
  42. static obj_hyperbolic P_((Now *np, Obj *op));
  43. static obj_parabolic P_((Now *np, Obj *op));
  44. static sun_cir P_((Now *np, Obj *op));
  45. static moon_cir P_((Now *np, Obj *op));
  46. static void cir_sky P_((Now *np, double lpd, double psi, double rp, double rho, double lam, double bet, double lsn, double rsn, Obj *op));
  47. static void elongation P_((double lam, double bet, double lsn, double *el));
  48.  
  49. #undef P_
  50.  
  51. /* given a Now and an Obj, fill in the approprirate s_* fields within Obj.
  52.  * return 0 if all ok, else -1.
  53.  */
  54. obj_cir (np, op)
  55. Now *np;
  56. Obj *op;
  57. {
  58.     switch (op->type) {
  59.     case FIXED:     return (obj_fixed (np, op));
  60.     case ELLIPTICAL: return (obj_elliptical (np, op));
  61.     case HYPERBOLIC: return (obj_hyperbolic (np, op));
  62.     case PARABOLIC:  return (obj_parabolic (np, op));
  63.     case PLANET:     return (obj_planet (np, op));
  64.     default:
  65.         printf ("obj_cir() called with type %d\n", op->type);
  66.         exit(1);
  67.         return (-1);    /* just for lint */
  68.     }
  69. }
  70.  
  71. static
  72. obj_planet (np, op)
  73. Now *np;
  74. Obj *op;
  75. {
  76.     double lsn, rsn;    /* true geoc lng of sun; dist from sn to earth*/
  77.     double lpd, psi;    /* heliocentric ecliptic long and lat */
  78.     double rp;        /* dist from sun */
  79.     double rho;        /* dist from earth */
  80.     double lam, bet;    /* geocentric ecliptic long and lat */
  81.     double dia, mag;    /* angular diameter at 1 AU and magnitude */
  82.     double f;        /* fractional phase from earth */
  83.     int p;
  84.  
  85.     /* validate code and check for a few special cases */
  86.     p = op->pl.code;
  87.     if (p < 0 || p > MOON) {
  88.         printf ("unknown planet code: %d \n", p);
  89.         exit(1);
  90.     }
  91.     else if (p == SUN)
  92.         return (sun_cir (np, op));
  93.     else if (p == MOON)
  94.         return (moon_cir (np, op));
  95.  
  96.  
  97.     /* find solar ecliptical longitude and distance to sun from earth */
  98.     sunpos (mjd, &lsn, &rsn);
  99.  
  100.     /* find helio long/lat; sun/planet and earth/plant dist; ecliptic
  101.      * long/lat; diameter and mag.
  102.      */
  103.     plans(mjd, p, &lpd, &psi, &rp, &rho, &lam, &bet, &dia, &mag);
  104.  
  105.     /* fill in all of op->s_* stuff except s_size and s_mag */
  106.     cir_sky (np, lpd, psi, rp, rho, lam, bet, lsn, rsn, op);
  107.  
  108.     /* compute magnitude and angular size */
  109.     f = 0.25 * (((rp+rho)*(rp+rho) - rsn*rsn)/(rp*rho));
  110.     op->s_mag = MAGSCALE * (mag + 5.0*log(rp*rho/sqrt(f))/log(10.0));
  111.     op->s_size = dia/rho + 0.5;
  112.  
  113.     return (0);
  114. }
  115.  
  116. static
  117. obj_fixed (np, op)
  118. Now *np;
  119. Obj *op;
  120. {
  121.     double lsn, rsn;    /* true geoc lng of sun, dist from sn to earth*/
  122.     double lam, bet;    /* geocentric ecliptic long and lat */
  123.     double lst;        /* local sidereal time */
  124.     double ha;        /* local hour angle */
  125.     double el;        /* elongation */
  126.     double alt, az;        /* current alt, az */
  127.     double ra, dec;        /* ra and dec at epoch of date */
  128.  
  129.     /* set ra/dec to their values at epoch of date */
  130.     ra = op->f_RA;
  131.     dec = op->f_dec;
  132.     precess (op->f_epoch, mjd, &ra, &dec);
  133.  
  134.     /* set s_ra/dec at desired epoch.  */
  135.     if (epoch == EOD) {
  136.         op->s_ra = ra;
  137.         op->s_dec = dec;
  138.     } else {
  139.         /* want a certain epoch -- if it's not what the database is at
  140.          * we change the original to save time next time assuming the
  141.          * user is likely to stick with this for a while.
  142.          */
  143.         if ((float)epoch != op->f_epoch) {
  144.         double tra = op->f_RA, tdec = op->f_dec;
  145.         float tepoch = epoch;    /* compare to float precision */
  146.         precess (op->f_epoch, tepoch, &tra, &tdec);
  147.         op->f_epoch = tepoch;
  148.         op->f_RA = tra;
  149.         op->f_dec = tdec;
  150.         }
  151.         op->s_ra = op->f_RA;
  152.         op->s_dec = op->f_dec;
  153.     }
  154.  
  155.     /* convert equitorial ra/dec to geocentric ecliptic lat/long */
  156.     eq_ecl (mjd, ra, dec, &bet, &lam);
  157.  
  158.     /* find solar ecliptical longitude and distance to sun from earth */
  159.     sunpos (mjd, &lsn, &rsn);
  160.  
  161.     /* compute elongation from ecliptic long/lat and sun geocentric long */
  162.     elongation (lam, bet, lsn, &el);
  163.     el = raddeg(el);
  164.     op->s_elong = el;
  165.  
  166.     /* these are really the same fields ...
  167.     op->s_mag = op->f_mag;
  168.     op->s_size = op->f_size;
  169.     */
  170.  
  171.     /* alt, az: correct for refraction; use eod ra/dec. */
  172.     now_lst (np, &lst);
  173.     ha = hrrad(lst) - ra;
  174.     hadec_aa (lat, ha, dec, &alt, &az);
  175.     refract (pressure, temp, alt, &alt);
  176.     op->s_alt = alt;
  177.     op->s_az = az;
  178.  
  179.     return (0);
  180. }
  181.  
  182. /* compute sky circumstances of an object in heliocentric elliptic orbit at *np.
  183.  */
  184. static
  185. obj_elliptical (np, op)
  186. Now *np;
  187. Obj *op;
  188. {
  189.     double lsn, rsn;    /* true geoc lng of sun; dist from sn to earth*/
  190.     double dt;        /* light travel time to object */
  191.     double lg;        /* helio long of earth */
  192.     double nu, ea;        /* true anomaly and eccentric anomaly */
  193.     double ma;        /* mean anomaly */
  194.     double rp;        /* distance from the sun */
  195.     double lo, slo, clo;    /* angle from ascending node */
  196.     double inc;        /* inclination */
  197.     double psi, spsi, cpsi;    /* heliocentric latitude */
  198.     double lpd;         /* heliocentric longitude */
  199.     double rho;        /* distance from the Earth */
  200.     double om;        /* arg of perihelion */
  201.     double Om;        /* long of ascending node. */
  202.     double lam;            /* geocentric ecliptic longitude */
  203.     double bet;            /* geocentric ecliptic latitude */
  204.     double e;        /* fast eccentricity */
  205.     double ll, sll, cll;    /* helio angle between object and earth */
  206.     double mag;        /* magnitude */
  207.     double rpd;
  208.     double y;
  209.     int pass;
  210.  
  211.     /* find location of earth from sun now */
  212.     sunpos (mjd, &lsn, &rsn);
  213.     lg = lsn + PI;
  214.  
  215.     /* faster access to eccentricty */
  216.     e = op->e_e;
  217.  
  218.     /* mean daily motion is optional -- fill in from period if 0 */
  219.     if (op->e_n == 0.0) {
  220.         double a = op->e_a;
  221.         op->e_n = 0.9856076686/sqrt(a*a*a);
  222.     }
  223.  
  224.     /* correct for light time by computing position at time mjd, then
  225.      *   again at mjd-dt, where
  226.      *   dt = time it takes light to travel earth-object distance.
  227.      * this is basically the same code as pelement() and plans()
  228.      *   combined and simplified for the special case of osculating
  229.      *   (unperturbed) elements. we have added reduction of elements using
  230.      *   reduce_elements().
  231.      */
  232.     dt = 0;
  233.     for (pass = pref_get(PREF_ALGO)==PREF_ACCURATE ? 0 : 1; pass<2; pass++){
  234.  
  235.         reduce_elements (op->e_epoch, mjd-dt, degrad(op->e_inc),
  236.                 degrad (op->e_om), degrad (op->e_Om),
  237.                 &inc, &om, &Om);
  238.  
  239.         ma = degrad (op->e_M + (mjd - op->e_cepoch - dt) * op->e_n);
  240.         anomaly (ma, e, &nu, &ea);
  241.         rp = op->e_a * (1-e*e) / (1+e*cos(nu));
  242.         lo = nu + om;
  243.         slo = sin(lo);
  244.         clo = cos(lo);
  245.         spsi = slo*sin(inc);
  246.         y = slo*cos(inc);
  247.         psi = asin(spsi);
  248.         lpd = atan(y/clo)+Om;
  249.         if (clo<0) lpd += PI;
  250.         range (&lpd, 2*PI);
  251.         cpsi = cos(psi);
  252.         rpd = rp*cpsi;
  253.         ll = lpd-lg;
  254.         rho = sqrt(rsn*rsn+rp*rp-2*rsn*rp*cpsi*cos(ll));
  255.  
  256.         dt = rho*LTAU/3600.0/24.0;    /* light travel time, in days / AU */
  257.     }
  258.  
  259.     /* compute sin and cos of ll */
  260.     sll = sin(ll);
  261.     cll = cos(ll);
  262.  
  263.     /* find geocentric ecliptic longitude and latitude */
  264.     if (rpd < rsn)
  265.         lam = atan(-1*rpd*sll/(rsn-rpd*cll))+lg+PI;
  266.     else
  267.         lam = atan(rsn*sll/(rpd-rsn*cll))+lpd;
  268.     range (&lam, 2*PI);
  269.     bet = atan(rpd*spsi*sin(lam-lpd)/(cpsi*rsn*sll));
  270.  
  271.     /* fill in all of op->s_* stuff except s_size and s_mag */
  272.     cir_sky (np, lpd, psi, rp, rho, lam, bet, lsn, rsn, op);
  273.  
  274.     /* compute magnitude and size */
  275.     if (op->e_mag.whichm == MAG_HG) {
  276.         /* the H and G parameters from the Astro. Almanac.
  277.          */
  278.         hg_mag (op->e_mag.m1, op->e_mag.m2, rp, rho, rsn, &mag);
  279.     } else {
  280.         /* the g/k model of comets */
  281.         gk_mag (op->e_mag.m1, op->e_mag.m2, rp, rho, &mag);
  282.     }
  283.     op->s_mag = mag * MAGSCALE;
  284.     op->s_size = op->e_size / rho;
  285.  
  286.     return (0);
  287. }
  288.  
  289. /* compute sky circumstances of an object in heliocentric hyperbolic orbit.
  290.  */
  291. static
  292. obj_hyperbolic (np, op)
  293. Now *np;
  294. Obj *op;
  295. {
  296.     double lsn, rsn;    /* true geoc lng of sun; dist from sn to earth*/
  297.     double dt;        /* light travel time to object */
  298.     double lg;        /* helio long of earth */
  299.     double nu, ea;        /* true anomaly and eccentric anomaly */
  300.     double ma;        /* mean anomaly */
  301.     double rp;        /* distance from the sun */
  302.     double lo, slo, clo;    /* angle from ascending node */
  303.     double inc;        /* inclination */
  304.     double psi, spsi, cpsi;    /* heliocentric latitude */
  305.     double lpd;         /* heliocentric longitude */
  306.     double rho;        /* distance from the Earth */
  307.     double om;        /* arg of perihelion */
  308.     double Om;        /* long of ascending node. */
  309.     double lam;            /* geocentric ecliptic longitude */
  310.     double bet;            /* geocentric ecliptic latitude */
  311.     double e;        /* fast eccentricity */
  312.     double ll, sll, cll;    /* helio angle between object and earth */
  313.     double n;        /* mean daily motion */
  314.     double mag;        /* magnitude */
  315.     double a;        /* mean distance */
  316.     double rpd;
  317.     double y;
  318.     int pass;
  319.  
  320.     /* find solar ecliptical longitude and distance to sun from earth */
  321.     sunpos (mjd, &lsn, &rsn);
  322.  
  323.     lg = lsn + PI;
  324.     e = op->h_e;
  325.     a = op->h_qp/(e - 1.0);
  326.     n = .98563/sqrt(a*a*a);
  327.  
  328.     /* correct for light time by computing position at time mjd, then
  329.      *   again at mjd-dt, where
  330.      *   dt = time it takes light to travel earth-object distance.
  331.      */
  332.     dt = 0;
  333.     for (pass = pref_get(PREF_ALGO)==PREF_ACCURATE ? 0 : 1; pass<2; pass++){
  334.  
  335.         reduce_elements (op->h_epoch, mjd-dt, degrad(op->h_inc),
  336.                 degrad (op->h_om), degrad (op->h_Om),
  337.                 &inc, &om, &Om);
  338.  
  339.         ma = degrad ((mjd - op->h_ep - dt) * n);
  340.         anomaly (ma, e, &nu, &ea);
  341.         rp = a * (e*e-1.0) / (1.0+e*cos(nu));
  342.         lo = nu + om;
  343.         slo = sin(lo);
  344.         clo = cos(lo);
  345.         spsi = slo*sin(inc);
  346.         y = slo*cos(inc);
  347.         psi = asin(spsi);
  348.         lpd = atan(y/clo)+Om;
  349.         if (clo<0) lpd += PI;
  350.         range (&lpd, 2*PI);
  351.         cpsi = cos(psi);
  352.         rpd = rp*cpsi;
  353.         ll = lpd-lg;
  354.         rho = sqrt(rsn*rsn+rp*rp-2*rsn*rp*cpsi*cos(ll));
  355.  
  356.         dt = rho*5.775518e-3;    /* light travel time, in days */
  357.     }
  358.  
  359.     /* compute sin and cos of ll */
  360.     sll = sin(ll);
  361.     cll = cos(ll);
  362.  
  363.     /* find geocentric ecliptic longitude and latitude */
  364.     if (rpd < rsn)
  365.         lam = atan(-1*rpd*sll/(rsn-rpd*cll))+lg+PI;
  366.     else
  367.         lam = atan(rsn*sll/(rpd-rsn*cll))+lpd;
  368.     range (&lam, 2*PI);
  369.     bet = atan(rpd*spsi*sin(lam-lpd)/(cpsi*rsn*sll));
  370.  
  371.     /* fill in all of op->s_* stuff except s_size and s_mag */
  372.     cir_sky (np, lpd, psi, rp, rho, lam, bet, lsn, rsn, op);
  373.  
  374.     /* compute magnitude and size */
  375.     gk_mag (op->h_g, op->h_k, rp, rho, &mag);
  376.     op->s_mag = mag * MAGSCALE;
  377.     op->s_size = op->h_size / rho;
  378.  
  379.     return (0);
  380. }
  381.  
  382. /* compute sky circumstances of an object in heliocentric hyperbolic orbit.
  383.  */
  384. static
  385. obj_parabolic (np, op)
  386. Now *np;
  387. Obj *op;
  388. {
  389.     double lsn, rsn;    /* true geoc lng of sun; dist from sn to earth*/
  390.     double lam;            /* geocentric ecliptic longitude */
  391.     double bet;            /* geocentric ecliptic latitude */
  392.     double mag;        /* magnitude */
  393.     double inc, om, Om;
  394.     double lpd, psi, rp, rho;
  395.     double dt;
  396.     int pass;
  397.  
  398.     /* find solar ecliptical longitude and distance to sun from earth */
  399.     sunpos (mjd, &lsn, &rsn);
  400.  
  401.     /* two passes to correct lam and bet for light travel time.e_ */
  402.     dt = 0.0;
  403.     for (pass = pref_get(PREF_ALGO)==PREF_ACCURATE ? 0 : 1; pass<2; pass++){
  404.         reduce_elements (op->p_epoch, mjd-dt, degrad(op->p_inc),
  405.         degrad(op->p_om), degrad(op->p_Om), &inc, &om, &Om);
  406.         comet (mjd-dt, op->p_ep, inc, om, op->p_qp, Om,
  407.                     &lpd, &psi, &rp, &rho, &lam, &bet);
  408.         dt = rho*LTAU/3600.0/24.0;    /* light travel time, in days / AU */
  409.     }
  410.  
  411.     /* fill in all of op->s_* stuff except s_size and s_mag */
  412.     cir_sky (np, lpd, psi, rp, rho, lam, bet, lsn, rsn, op);
  413.  
  414.     /* compute magnitude and size */
  415.     gk_mag (op->p_g, op->p_k, rp, rho, &mag);
  416.     op->s_mag = mag * MAGSCALE;
  417.     op->s_size = op->p_size / rho;
  418.  
  419.     return (0);
  420. }
  421.  
  422. /* find sun's circumstances now.
  423.  */
  424. static
  425. sun_cir (np, op)
  426. Now *np;
  427. Obj *op;
  428. {
  429.     double lsn, rsn;    /* true geoc lng of sun; dist from sn to earth*/
  430.     double deps, dpsi;    /* nutation on obliquity and longitude */
  431.     double lst;        /* local sidereal time */
  432.     double ehp;        /* angular diamter of earth from object */
  433.     double ha;        /* hour angle */
  434.     double ra, dec;        /* ra and dec now */
  435.     double alt, az;        /* alt and az */
  436.     double dhlong;
  437.  
  438.     sunpos (mjd, &lsn, &rsn);    /* sun's true ecliptic long and dist */
  439.     nutation (mjd, &deps, &dpsi);    /* correct for nutation ... */
  440.     lsn += dpsi;
  441.     lsn -= degrad(20.4/3600);    /* and aberation */
  442.  
  443.     op->s_edist = rsn;
  444.     op->s_sdist = 0.0;
  445.     op->s_elong = 0.0;
  446.     op->s_phase = 100.0;
  447.     op->s_size = raddeg(4.65242e-3/rsn)*3600*2;
  448.     op->s_mag = MAGSCALE * -26.8;    /* TODO */
  449.     dhlong = lsn-PI;    /* geo- to helio- centric */
  450.     range (&dhlong, 2*PI);
  451.     op->s_hlong = dhlong;
  452.     op->s_hlat = 0.0;
  453.  
  454.  
  455.     /* convert geocentric ecliptic lat/long to equitorial ra/dec of date */
  456.     ecl_eq (mjd, 0.0, lsn, &ra, &dec);
  457.  
  458.     /* find s_ra/dec at desired epoch */
  459.     if (epoch == EOD) {
  460.         op->s_ra = ra;
  461.         op->s_dec = dec;
  462.     } else {
  463.         double tra = ra, tdec = dec;
  464.         precess (mjd, epoch, &tra, &tdec);
  465.         op->s_ra = tra;
  466.         op->s_dec = tdec;
  467.     }
  468.  
  469.     /* find alt/az based on unprecessed ra/dec */
  470.     now_lst (np, &lst);
  471.     ha = hrrad(lst) - ra;
  472.     ehp = (2.0 * 6378.0 / 146.0e6) / op->s_edist;
  473.     ta_par (ha, dec, lat, elev, ehp, &ha, &dec);
  474.     hadec_aa (lat, ha, dec, &alt, &az);
  475.     refract (pressure, temp, alt, &alt);
  476.     op->s_alt = alt;
  477.     op->s_az = az;
  478.     return (0);
  479. }
  480.  
  481. /* find moon's circumstances now.
  482.  */
  483. static
  484. moon_cir (np, op)
  485. Now *np;
  486. Obj *op;
  487. {
  488.     double lsn, rsn;    /* true geoc lng of sun; dist from sn to earth*/
  489.     double deps, dpsi;    /* nutation on obliquity and longitude */
  490.     double lst;        /* local sidereal time */
  491.     double ehp;        /* angular diamter of earth from object */
  492.     double ha;        /* hour angle */
  493.     double ra, dec;        /* ra and dec now */
  494.     double alt, az;        /* alt and az */
  495.     double lam;            /* geocentric ecliptic longitude */
  496.     double bet;            /* geocentric ecliptic latitude */
  497.     double edistau;        /* earth-moon dist, in au */
  498.     double el;        /* elongation, rads east */
  499.  
  500.     moon (mjd, &lam, &bet, &ehp);    /* moon's true ecliptic loc */
  501.     nutation (mjd, &deps, &dpsi);    /* correct for nutation */
  502.     lam += dpsi;
  503.     range (&lam, 2*PI);
  504.  
  505.     op->s_edist = 6378.14/sin(ehp);    /* earth-moon dist, want km */
  506.     op->s_size = 3600*31.22512*sin(ehp);/* moon angular dia, seconds */
  507.     op->s_hlong = lam;        /* save geo in helio fields */
  508.     op->s_hlat = bet;
  509.  
  510.     /* convert geocentric ecliptic lat/long to equitorial ra/dec of date */
  511.     ecl_eq (mjd, bet, lam, &ra, &dec);
  512.  
  513.     /* find s_ra/dec at desired epoch */
  514.     if (epoch == EOD) {
  515.         op->s_ra = ra;
  516.         op->s_dec = dec;
  517.     } else {
  518.         double tra = ra, tdec = dec;
  519.         precess (mjd, epoch, &tra, &tdec);
  520.         op->s_ra = tra;
  521.         op->s_dec = tdec;
  522.     }
  523.  
  524.     sunpos (mjd, &lsn, &rsn);
  525.     range (&lsn, 2*PI);
  526.     elongation (lam, bet, lsn, &el);
  527.  
  528.     /* solve triangle of earth, sun, and elongation for moon-sun dist */
  529.     edistau = op->s_edist/1.495979e8;    /* km -> au */
  530.     op->s_sdist= sqrt (edistau*edistau + rsn*rsn - 2.0*edistau*rsn*cos(el));
  531.  
  532.     /* TODO: improve mag; this is based on a flat moon model. */
  533.     op->s_mag = MAGSCALE * 
  534.         (-12.7 + 2.5*(log10(PI) - log10(PI/2*(1+1.e-6-cos(el)))));
  535.  
  536.     op->s_elong = raddeg(el);        /* want degrees */
  537.     op->s_phase = fabs(el)/PI*100.0;    /* want non-negative % */
  538.  
  539.     /* show topocentric alt/az by correcting ra/dec for parallax 
  540.      * as well as refraction.
  541.      * use unprecessed ra/dec.
  542.      */
  543.     now_lst (np, &lst);
  544.     ha = hrrad(lst) - ra;
  545.     ta_par (ha, dec, lat, elev, ehp, &ha, &dec);
  546.     hadec_aa (lat, ha, dec, &alt, &az);
  547.     refract (pressure, temp, alt, &alt);
  548.     op->s_alt = alt;
  549.     op->s_az = az;
  550.  
  551.     return (0);
  552. }
  553.  
  554. /* fill in all of op->s_* stuff except s_size and s_mag */
  555. static void
  556. cir_sky (np, lpd, psi, rp, rho, lam, bet, lsn, rsn, op)
  557. Now *np;
  558. double lpd, psi;    /* heliocentric ecliptic long and lat */
  559. double rp;        /* dist from sun */
  560. double rho;        /* dist from earth */
  561. double lam, bet;    /* true geocentric ecliptic long and lat */
  562. double lsn, rsn;    /* true geoc lng of sun; dist from sn to earth*/
  563. Obj *op;
  564. {
  565.     double a;        /* geoc angle between object and sun */
  566.     double ra, dec;        /* ra and dec at epoch of date */
  567.     double el;        /* elongation */
  568.     double lst;        /* local sidereal time */
  569.     double ehp;        /* ehp: angular dia of earth from body */
  570.     double ha;        /* local hour angle */
  571.     double f;        /* fractional phase from earth */
  572.     double alt, az;        /* current alt, az */
  573.     double deps, dpsi;    /* nutation on obliquity and longitude */
  574.  
  575.     /* correct for nutation ... */
  576.     nutation (mjd, &deps, &dpsi);
  577.     lam += dpsi;
  578.  
  579.     /* and correct for 20.4" aberation */
  580.     a = lsn - lam;
  581.     lam -= degrad(20.4/3600)*cos(a)/cos(bet);
  582.     bet -= degrad(20.4/3600)*sin(a)*sin(bet);
  583.  
  584.     /* convert geocentric ecliptic lat/long to equitorial ra/dec of date */
  585.     ecl_eq (mjd, bet, lam, &ra, &dec);
  586.  
  587.     /* set s_ra/s_dec to that of desired display epoch. */
  588.     if (epoch == EOD) {
  589.         op->s_ra = ra;
  590.         op->s_dec = dec;
  591.     } else {
  592.         double tra = ra, tdec = dec;
  593.         precess (mjd, epoch, &tra, &tdec);
  594.         op->s_ra = tra;
  595.         op->s_dec = tdec;
  596.     }
  597.  
  598.     /* set earth/planet and sun/planet distance */
  599.     op->s_edist = rho;
  600.     op->s_sdist = rp;
  601.  
  602.     /* compute elongation and phase */
  603.     elongation (lam, bet, lsn, &el);
  604.     el = raddeg(el);
  605.     op->s_elong = el;
  606.     f = 0.25 * ((rp+rho)*(rp+rho) - rsn*rsn)/(rp*rho);
  607.     op->s_phase = f*100.0; /* percent */
  608.  
  609.     /* set heliocentric long/lat */
  610.     op->s_hlong = lpd;
  611.     op->s_hlat = psi;
  612.  
  613.     /* alt, az: correct for parallax and refraction; use eod ra/dec. */
  614.     now_lst (np, &lst);
  615.     ha = hrrad(lst) - ra;
  616.     ehp = (2.0*6378.0/146.0e6) / rho;
  617.     ta_par (ha, dec, lat, elev, ehp, &ha, &dec);
  618.     hadec_aa (lat, ha, dec, &alt, &az);
  619.     refract (pressure, temp, alt, &alt);
  620.     op->s_alt = alt;
  621.     op->s_az = az;
  622. }
  623.  
  624. /* given geocentric ecliptic longitude and latitude, lam and bet, of some object
  625.  * and the longitude of the sun, lsn, find the elongation, el. this is the
  626.  * actual angular separation of the object from the sun, not just the difference
  627.  * in the longitude. the sign, however, IS set simply as a test on longitude
  628.  * such that el will be >0 for an evening object <0 for a morning object.
  629.  * to understand the test for el sign, draw a graph with lam going from 0-2*PI
  630.  *   down the vertical axis, lsn going from 0-2*PI across the hor axis. then
  631.  *   define the diagonal regions bounded by the lines lam=lsn+PI, lam=lsn and
  632.  *   lam=lsn-PI. the "morning" regions are any values to the lower left of the
  633.  *   first line and bounded within the second pair of lines.
  634.  * all angles in radians.
  635.  */
  636. static void
  637. elongation (lam, bet, lsn, el)
  638. double lam, bet, lsn;
  639. double *el;
  640. {
  641.     *el = acos(cos(bet)*cos(lam-lsn));
  642.     if (lam>lsn+PI || lam>lsn-PI && lam<lsn) *el = - *el;
  643. }
  644.