home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / x / volume19 / xephem / part15 / riset_cir.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-05-15  |  7.7 KB  |  252 lines
  1. /* find rise and set circumstances, ie, riset_cir() and related functions. */
  2.  
  3. #include <stdio.h>
  4. #include <math.h>
  5. #if defined(__STDC__)
  6. #include <stdlib.h>
  7. #endif
  8. #include "astro.h"
  9. #include "circum.h"
  10.  
  11. #if defined(__STDC__) || defined(__cplusplus)
  12. #define P_(s) s
  13. #else
  14. #define P_(s) ()
  15. #endif
  16.  
  17. extern double mjd_day P_((double jd));
  18. extern int obj_cir P_((Now *np, Obj *op));
  19. extern void now_lst P_((Now *np, double *lst));
  20. extern void range P_((double *v, double r));
  21. extern void refract P_((double pr, double tr, double ta, double *aa));
  22. extern void riset P_((double ra, double dec, double lati, double dis, double *lstr, double *lsts, double *azr, double *azs, int *status));
  23. extern void riset_cir P_((Now *np, Obj *op, double dis, RiseSet *rp));
  24. extern void ta_par P_((double tha, double tdec, double phi, double ht, double ehp, double *aha, double *adec));
  25. extern void unrefract P_((double pr, double tr, double aa, double *ta));
  26.  
  27. void riset_cir P_((Now *np, Obj *op, double dis, RiseSet *rp));
  28. void twilight_cir P_((Now *np, double dis, double *dawn, double *dusk, int *status));
  29. static void stationary_riset P_((Obj *op, Now *np, double dis, double *lstr, double *lsts, double *lstt, double *azr, double *azs, double *altt, int *status));
  30. static void transit P_((double r, double d, Now *np, double *lstt, double *altt));
  31.  
  32. #undef P_
  33.  
  34.  
  35. #define    TMACC    (15./3600.)    /* convergence accuracy, hours */
  36.  
  37. /* find where and when an object, op, will rise and set and
  38.  *   it's transit circumstances. all times are local, angles rads e of n.
  39.  * dis is the angle down from an ideal horizon, in rads (see riset()).
  40.  */
  41. void
  42. riset_cir (np, op, dis, rp)
  43. Now *np;
  44. Obj *op;
  45. double dis;
  46. RiseSet *rp;
  47. {
  48. #define    MAXPASSES    6
  49.     double lstr, lsts, lstt; /* local sidereal times of rising/setting */
  50.     double lnoon;        /* mjd of local noon */
  51.     double x;        /* discarded tmp value */
  52.     Now n;            /* copy to move time around */
  53.     Obj o;            /* copy to get circumstances at n */
  54.     double lst;        /* lst at local noon */
  55.     double diff, lastdiff;    /* iterative improvement to mjd0 */
  56.     int pass;
  57.     int rss;
  58.  
  59.     /* find mjd of local noon -- all iterations start then */
  60.     lnoon = mjd_day(mjd - tz/24.0) + (12.0 + tz)/24.0;
  61.  
  62.     /* work with local copies so we can move the time around */
  63.     n = *np;
  64.     o = *op;
  65.  
  66.     /* assume no problems initially */
  67.     rp->rs_flags = 0;
  68.  
  69.     /* first approximation is to find rise/set times of a fixed object
  70.      * at the current epoch in its position at local noon.
  71.      */
  72.     n.n_mjd = lnoon;
  73.     n.n_epoch = EOD;
  74.     now_lst (&n, &lst);    /* lst at local noon */
  75.     stationary_riset (&o, &n, dis, &lstr, &lsts, &lstt, &x, &x, &x, &rss);
  76.  
  77.     chkrss:    /* come here to check rss status and control next step */
  78.     switch (rss) {
  79.     case  0:  break;
  80.     case  1: rp->rs_flags = RS_NEVERUP; return;
  81.     case -1: rp->rs_flags = RS_CIRCUMPOLAR; goto dotransit;
  82.     default: rp->rs_flags = RS_ERROR; return;
  83.     }
  84.  
  85.     /* find a better approximation to the rising circumstances based on
  86.      * more passes, each using a "fixed" object at the location at
  87.      * previous approximation of the rise time.
  88.      */
  89.     lastdiff = 1000.0;
  90.     for (pass = 1; pass < MAXPASSES; pass++) {
  91.         diff = (lstr - lst)*SIDRATE; /* next guess at rise time wrt noon */
  92.         if (diff > 12.0)
  93.         diff -= 24.0*SIDRATE;    /* not tomorrow, today */
  94.         else if (diff < -12.0)
  95.         diff += 24.0*SIDRATE;    /* not yesterday, today */
  96.         n.n_mjd = lnoon + diff/24.0;    /* next guess at mjd of rise */
  97.         stationary_riset (&o, &n, dis, &lstr, &x, &x, &rp->rs_riseaz,
  98.                                 &x, &x, &rss);
  99.         if (rss != 0) goto chkrss;
  100.         if (fabs (diff - lastdiff) < TMACC)
  101.         break;
  102.         lastdiff = diff;
  103.     }
  104.     if (pass == MAXPASSES)
  105.         rp->rs_flags |= RS_NORISE;    /* didn't converge - no rise today */
  106.     else {
  107.         rp->rs_risetm = 12.0 + diff;
  108.         if (!is_planet(op, MOON) &&
  109.             (rp->rs_risetm <= 24.0*(1.0-SIDRATE)
  110.             || rp->rs_risetm >= 24.0*SIDRATE))
  111.         rp->rs_flags |= RS_2RISES;
  112.     }
  113.  
  114.     /* find a better approximation to the setting circumstances based on
  115.      * more passes, each using a "fixed" object at the location at
  116.      * previous approximation of the set time.
  117.      */
  118.     lastdiff = 1000.0;
  119.     for (pass = 1; pass < MAXPASSES; pass++) {
  120.         diff = (lsts - lst)*SIDRATE; /* next guess at set time wrt noon */
  121.         if (diff > 12.0)
  122.         diff -= 24.0*SIDRATE;    /* not tomorrow, today */
  123.         else if (diff < -12.0)
  124.         diff += 24.0*SIDRATE;    /* not yesterday, today */
  125.         n.n_mjd = lnoon + diff/24.0;    /* next guess at mjd of set */
  126.         stationary_riset (&o, &n, dis, &x, &lsts, &x, &x,
  127.                             &rp->rs_setaz, &x, &rss);
  128.         if (rss != 0) goto chkrss;
  129.         if (fabs (diff - lastdiff) < TMACC)
  130.         break;
  131.         lastdiff = diff;
  132.     }
  133.     if (pass == MAXPASSES)
  134.         rp->rs_flags |= RS_NOSET;    /* didn't converge - no set today */
  135.     else {
  136.         rp->rs_settm = 12.0 + diff;
  137.         if (!is_planet(op, MOON) &&
  138.             (rp->rs_settm <= 24.0*(1.0-SIDRATE)
  139.             || rp->rs_settm >= 24.0*SIDRATE))
  140.         rp->rs_flags |= RS_2SETS;
  141.     }
  142.  
  143.     dotransit:
  144.     /* find a better approximation to the transit circumstances based on
  145.      * more passes, each using a "fixed" object at the location at
  146.      * previous approximation of the transit time.
  147.      */
  148.     lastdiff = 1000.0;
  149.     for (pass = 1; pass < MAXPASSES; pass++) {
  150.         diff = (lstt - lst)*SIDRATE; /*next guess at transit time wrt noon*/
  151.         if (diff > 12.0)
  152.         diff -= 24.0*SIDRATE;    /* not tomorrow, today */
  153.         else if (diff < -12.0)
  154.         diff += 24.0*SIDRATE;    /* not yesterday, today */
  155.         n.n_mjd = lnoon + diff/24.0; /* next guess at mjd of transit */
  156.         stationary_riset (&o, &n, dis, &x, &x, &lstt, &x, &x,
  157.                             &rp->rs_tranalt, &rss);
  158.         if (fabs (diff - lastdiff) < TMACC)
  159.         break;
  160.         lastdiff = diff;
  161.     }
  162.     if (pass == MAXPASSES)
  163.         rp->rs_flags |= RS_NOTRANS;    /* didn't converge - no transit today */
  164.     else {
  165.         rp->rs_trantm = 12.0 + diff;
  166.         if (!is_planet(op, MOON) &&
  167.             (rp->rs_trantm <= 24.0*(1.0-SIDRATE)
  168.             || rp->rs_trantm >= 24.0*SIDRATE))
  169.         rp->rs_flags |= RS_2TRANS;
  170.     }
  171. }
  172.  
  173. /* find local times when sun is dis rads below horizon.
  174.  */
  175. void
  176. twilight_cir (np, dis, dawn, dusk, status)
  177. Now *np;
  178. double dis;
  179. double *dawn, *dusk;
  180. int *status;
  181. {
  182.     RiseSet rs;
  183.     Obj o;
  184.  
  185.     o.type = PLANET;
  186.     o.pl.code = SUN;
  187.     riset_cir (np, &o, dis, &rs);
  188.     *dawn = rs.rs_risetm;
  189.     *dusk = rs.rs_settm;
  190.     *status = rs.rs_flags;
  191. }
  192.  
  193. /* find the local rise/set/transit circumstances of a fixed object, *op.
  194.  * use lp to decide the day and the location circumstances.
  195.  * fill *status is have any problems.
  196.  */
  197. static void
  198. stationary_riset (op, np, dis, lstr, lsts, lstt, azr, azs, altt, status)
  199. Obj *op;
  200. Now *np;
  201. double dis;
  202. double *lstr, *lsts, *lstt;
  203. double *azr, *azs, *altt;
  204. int *status;
  205. {
  206.     /* find object op's topocentric ra/dec at np..
  207.      * (this must include parallax if it's in the solar system).
  208.      */
  209.     if (obj_cir (np, op) < 0) {
  210.         printf ("stationary_riset: can't get object loc\n");
  211.         exit (1);
  212.     }
  213.     
  214.     if (is_ssobj(op)) {
  215.         /* obj_cir() gives geocentric ra/dec; we need to account for
  216.          * parallax in solar system objects to get topocentric ra/dec.
  217.          */
  218.         double tra, tdec;
  219.         double ehp, lst, ha;
  220.         if (is_planet(op, MOON))
  221.         ehp = asin (6378.14/op->s_edist);
  222.         else
  223.         ehp = (2.*6378./146e6)/op->s_edist;
  224.         now_lst (np, &lst);
  225.         ha = hrrad(lst) - op->s_ra;
  226.         ta_par (ha, op->s_dec, lat, elev, ehp, &ha, &tdec);
  227.         tra = hrrad(lst) - ha;
  228.         range (&tra, 2*PI);
  229.         op->s_ra = tra;
  230.         op->s_dec = tdec;
  231.     }
  232.  
  233.     riset (op->s_ra, op->s_dec, lat, dis, lstr, lsts, azr, azs, status);
  234.     transit (op->s_ra, op->s_dec, np, lstt, altt);
  235. }
  236.  
  237.  
  238. /* find when and how hi object at (r,d) is when it transits. */
  239. static void
  240. transit (r, d, np, lstt, altt)
  241. double r, d;    /* ra and dec, rads */
  242. Now *np;    /* for refraction info */
  243. double *lstt;    /* local sidereal time of transit */
  244. double *altt;    /* local, refracted, altitude at time of transit */
  245. {
  246.     *lstt = radhr(r);
  247.     *altt = PI/2 - lat + d;
  248.     if (*altt > PI/2)
  249.         *altt = PI - *altt;
  250.     refract (pressure, temp, *altt, altt);
  251. }
  252.