home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS - Coast to Coast / simteldosarchivecoasttocoast2.iso / astrnomy / ephem421.zip / RISET_C.C < prev    next >
C/C++ Source or Header  |  1990-09-13  |  9KB  |  276 lines
  1. /* find rise and set circumstances, ie, riset_cir() and related functions. */
  2.  
  3. #include <stdio.h>
  4. #include <math.h>
  5. #include "astro.h"
  6. #include "circum.h"
  7. #include "screen.h"    /* just for SUN and MOON */
  8.  
  9. #define    TRACE(x)    {FILE *fp = fopen("trace","a"); fprintf x; fclose(fp);}
  10.  
  11. #define    STDREF    degrad(34./60.)    /* nominal horizon refraction amount */
  12. #define    TWIREF    degrad(18.)    /* twilight horizon displacement */
  13. #define    TMACC    (15./3600.)    /* convergence accuracy, hours */
  14.  
  15. /* find where and when a body, p, will rise and set and
  16.  *   its transit circumstances. all times are local, angles rads e of n.
  17.  * return 0 if just returned same stuff as previous call, else 1 if new.
  18.  * status is set from the RS_* #defines in circum.h.
  19.  * also used to find astro twilight by calling with dis of 18 degrees.
  20.  */
  21. riset_cir (p, np, force, hzn, ltr, lts, ltt, azr, azs, altt, status)
  22. int p;        /* one of the body defines in astro.h or screen.h */
  23. Now *np;
  24. int force;    /* set !=0 to force computations */
  25. int hzn;    /* STDHZN or ADPHZN or TWILIGHT */
  26. double *ltr, *lts; /* local rise and set times */
  27. double *ltt;    /* local transit time */
  28. double *azr, *azs; /* local rise and set azimuths, rads e of n */
  29. double *altt;    /* local altitude at transit */
  30. int *status;    /* one or more of the RS_* defines */
  31. {
  32.     typedef struct {
  33.         Now l_now;
  34.         double l_ltr, l_lts, l_ltt, l_azr, l_azs, l_altt;
  35.         int l_hzn;
  36.         int l_status;
  37.     } Last;
  38.     /* must be in same order as the astro.h/screen.h #define's */
  39.     static Last last[NOBJ] = {
  40.         {NOMJD}, {NOMJD}, {NOMJD}, {NOMJD}, {NOMJD}, {NOMJD},
  41.         {NOMJD}, {NOMJD}, {NOMJD}, {NOMJD}, {NOMJD}, {NOMJD}
  42.     };
  43.     Last *lp;
  44.     int new;
  45.  
  46.     lp = last + p;
  47.     if (!force && same_cir (np, &lp->l_now) && same_lday (np, &lp->l_now)
  48.                         && lp->l_hzn == hzn) {
  49.         *ltr = lp->l_ltr;
  50.         *lts = lp->l_lts;
  51.         *ltt = lp->l_ltt;
  52.         *azr = lp->l_azr;
  53.         *azs = lp->l_azs;
  54.         *altt = lp->l_altt;
  55.         *status = lp->l_status;
  56.         new = 0;
  57.     } else {
  58.         *status = 0;
  59.         iterative_riset (p, np, hzn, ltr, lts, ltt, azr, azs, altt, status);
  60.         lp->l_ltr = *ltr;
  61.         lp->l_lts = *lts;
  62.         lp->l_ltt = *ltt;
  63.         lp->l_azr = *azr;
  64.         lp->l_azs = *azs;
  65.         lp->l_altt = *altt;
  66.         lp->l_status = *status;
  67.         lp->l_hzn = hzn;
  68.         lp->l_now = *np;
  69.         new = 1;
  70.     }
  71.     return (new);
  72. }
  73.  
  74. static
  75. iterative_riset (p, np, hzn, ltr, lts, ltt, azr, azs, altt, status)
  76. int p;
  77. Now *np;
  78. int hzn;
  79. double *ltr, *lts, *ltt;    /* local times of rise, set and transit */
  80. double *azr, *azs, *altt;/* local azimuths of rise, set and transit altitude */
  81. int *status;
  82. {
  83. #define    MAXPASSES    6
  84.     double lstr, lsts, lstt; /* local sidereal times of rising/setting */
  85.     double mjd0;        /* mjd estimates of rise/set event */
  86.     double lnoon;        /* mjd of local noon */
  87.     double x;        /* discarded tmp value */
  88.     Now n;            /* just used to call now_lst() */
  89.     double lst;        /* lst at local noon */
  90.     double diff, lastdiff;    /* iterative improvement to mjd0 */
  91.     int pass;
  92.     int rss;
  93.  
  94.     /* first approximation is to find rise/set times of a fixed object
  95.      * in its position at local noon.
  96.      */
  97.     lnoon = mjd_day(mjd - tz/24.0) + (12.0 + tz)/24.0; /*mjd of local noon*/
  98.     n.n_mjd = lnoon;
  99.     n.n_lng = lng;
  100.     now_lst (&n, &lst);    /* lst at local noon */
  101.     mjd0 = lnoon;
  102.     stationary_riset (p,mjd0,np,hzn,&lstr,&lsts,&lstt,&x,&x,&x,&rss);
  103.     chkrss:
  104.     switch (rss) {
  105.     case  0:  break;
  106.     case  1: *status = RS_NEVERUP; return;
  107.     case -1: *status = RS_CIRCUMPOLAR; goto transit;
  108.     default: *status = RS_ERROR; return;
  109.     }
  110.  
  111.     /* find a better approximation to the rising circumstances based on
  112.      * more passes, each using a "fixed" object at the location at
  113.      * previous approximation of the rise time.
  114.      */
  115.     lastdiff = 1000.0;
  116.     for (pass = 1; pass < MAXPASSES; pass++) {
  117.         diff = (lstr - lst)*SIDRATE; /* next guess at rise time wrt noon */
  118.         if (diff > 12.0)
  119.         diff -= 24.0*SIDRATE;    /* not tomorrow, today */
  120.         else if (diff < -12.0)
  121.         diff += 24.0*SIDRATE;    /* not yesterday, today */
  122.         mjd0 = lnoon + diff/24.0;    /* next guess at mjd of rise */
  123.         stationary_riset (p,mjd0,np,hzn,&lstr,&x,&x,azr,&x,&x,&rss);
  124.         if (rss != 0) goto chkrss;
  125.         if (fabs (diff - lastdiff) < TMACC)
  126.         break;
  127.         lastdiff = diff;
  128.     }
  129.     if (pass == MAXPASSES)
  130.         *status |= RS_NORISE;    /* didn't converge - no rise today */
  131.     else {
  132.         *ltr = 12.0 + diff;
  133.         if (p != MOON &&
  134.             (*ltr <= 24.0*(1.0-SIDRATE) || *ltr >= 24.0*SIDRATE))
  135.         *status |= RS_2RISES;
  136.     }
  137.  
  138.     /* find a better approximation to the setting circumstances based on
  139.      * more passes, each using a "fixed" object at the location at
  140.      * previous approximation of the set time.
  141.      */
  142.     lastdiff = 1000.0;
  143.     for (pass = 1; pass < MAXPASSES; pass++) {
  144.         diff = (lsts - lst)*SIDRATE; /* next guess at set time wrt noon */
  145.         if (diff > 12.0)
  146.         diff -= 24.0*SIDRATE;    /* not tomorrow, today */
  147.         else if (diff < -12.0)
  148.         diff += 24.0*SIDRATE;    /* not yesterday, today */
  149.         mjd0 = lnoon + diff/24.0;    /* next guess at mjd of set */
  150.         stationary_riset (p,mjd0,np,hzn,&x,&lsts,&x,&x,azs,&x,&rss);
  151.         if (rss != 0) goto chkrss;
  152.         if (fabs (diff - lastdiff) < TMACC)
  153.         break;
  154.         lastdiff = diff;
  155.     }
  156.     if (pass == MAXPASSES)
  157.         *status |= RS_NOSET;    /* didn't converge - no set today */
  158.     else {
  159.         *lts = 12.0 + diff;
  160.         if (p != MOON &&
  161.             (*lts <= 24.0*(1.0-SIDRATE) || *lts >= 24.0*SIDRATE))
  162.         *status |= RS_2SETS;
  163.     }
  164.  
  165.     transit:
  166.     /* find a better approximation to the transit circumstances based on
  167.      * more passes, each using a "fixed" object at the location at
  168.      * previous approximation of the transit time.
  169.      */
  170.     lastdiff = 1000.0;
  171.     for (pass = 1; pass < MAXPASSES; pass++) {
  172.         diff = (lstt - lst)*SIDRATE; /*next guess at transit time wrt noon*/
  173.         if (diff > 12.0)
  174.         diff -= 24.0*SIDRATE;    /* not tomorrow, today */
  175.         else if (diff < -12.0)
  176.         diff += 24.0*SIDRATE;    /* not yesterday, today */
  177.         mjd0 = lnoon + diff/24.0;    /* next guess at mjd of transit */
  178.         stationary_riset (p,mjd0,np,hzn,&x,&x,&lstt,&x,&x,altt,&rss);
  179.         if (fabs (diff - lastdiff) < TMACC)
  180.         break;
  181.         lastdiff = diff;
  182.     }
  183.     if (pass == MAXPASSES)
  184.         *status |= RS_NOTRANS;    /* didn't converge - no transit today */
  185.     else {
  186.         *ltt = 12.0 + diff;
  187.         if (p != MOON &&
  188.             (*ltt <= 24.0*(1.0-SIDRATE) || *ltt >= 24.0*SIDRATE))
  189.         *status |= RS_2TRANS;
  190.     }
  191. }
  192.  
  193. static
  194. stationary_riset (p, mjd0, np, hzn, lstr, lsts, lstt, azr, azs, altt, status)
  195. int p;
  196. double mjd0;
  197. Now *np;
  198. int hzn;
  199. double *lstr, *lsts, *lstt;
  200. double *azr, *azs, *altt;
  201. int *status;
  202. {
  203.     extern void bye();
  204.     double dis;
  205.     Now n;
  206.     Sky s;
  207.  
  208.     /* find object p's topocentric ra/dec at mjd0
  209.      * (this must include parallax)
  210.      */
  211.     n = *np;
  212.     n.n_mjd = mjd0;
  213.     (void) body_cir (p, 0.0, &n, &s);
  214.     if (epoch != EOD)
  215.         precess (epoch, mjd0, &s.s_ra, &s.s_dec);
  216.     if (s.s_edist > 0) {
  217.         /* parallax, if we can */
  218.         double ehp, lst, ha;
  219.         if (p == MOON)
  220.         ehp = asin (6378.14/s.s_edist);
  221.         else
  222.         ehp = (2.*6378./146e6)/s.s_edist;
  223.         now_lst (&n, &lst);
  224.         ha = hrrad(lst) - s.s_ra;
  225.         ta_par (ha, s.s_dec, lat, height, ehp, &ha, &s.s_dec);
  226.         s.s_ra = hrrad(lst) - ha;
  227.         range (&s.s_ra, 2*PI);
  228.     }
  229.  
  230.     switch (hzn) {
  231.     case STDHZN:
  232.         /* nominal atmospheric refraction.
  233.          * then add nominal moon or sun semi-diameter, as appropriate.
  234.          * other objects assumes to be negligibly small.
  235.          */
  236.         dis = STDREF;
  237.         if (p == MOON || p == SUN)
  238.         dis += degrad (32./60./2.);
  239.         break;
  240.     case TWILIGHT:
  241.         if (p != SUN) {
  242.         f_msg ("Non-sun twilight bug!");
  243.         bye();
  244.         }
  245.         dis = TWIREF;
  246.         break;
  247.     case ADPHZN:
  248.         /* adaptive includes actual refraction conditions and also
  249.          * includes object's semi-diameter.
  250.          */
  251.         unrefract (pressure, temp, 0.0, &dis);
  252.         dis = -dis;
  253.         dis += degrad(s.s_size/3600./2.0);
  254.         break;
  255.     }
  256.  
  257.     riset (s.s_ra, s.s_dec, lat, dis, lstr, lsts, azr, azs, status);
  258.     transit (s.s_ra, s.s_dec, np, lstt, altt);
  259. }
  260.  
  261.  
  262. /* find when and how hi object at (r,d) is when it transits. */
  263. static
  264. transit (r, d, np, lstt, altt)
  265. double r, d;    /* ra and dec, rads */
  266. Now *np;    /* for refraction info */
  267. double *lstt;    /* local sidereal time of transit */
  268. double *altt;    /* local, refracted, altitude at time of transit */
  269. {
  270.     *lstt = radhr(r);
  271.     *altt = PI/2 - lat + d;
  272.     if (*altt > PI/2)
  273.         *altt = PI - *altt;
  274.     refract (pressure, temp, *altt, altt);
  275. }
  276.