home *** CD-ROM | disk | FTP | other *** search
/ rtsi.com / 2014.01.www.rtsi.com.tar / www.rtsi.com / OS9 / OSK / EFFO / pd8.lzh / SRC / riset_c.c < prev    next >
Text File  |  1990-04-13  |  11KB  |  352 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.  
  14. /* find where and when a body, p, will rise and set and
  15.  *   its transit circumstances. all times are local, angles rads e of n.
  16.  * return 0 if just returned same stuff as previous call, else 1 if new.
  17.  * status is set from the RS_* #defines in circum.h.
  18.  * also used to find astro twilight by calling with dis of 18 degrees.
  19.  */
  20. riset_cir (p, np, force, hzn, ltr, lts, ltt, azr, azs, altt, status)
  21. int p;        /* one of the body defines in astro.h or screen.h */
  22. Now *np;
  23. int force;    /* set !=0 to force computations */
  24. int hzn;    /* STDHZN or ADPHZN or TWILIGHT */
  25. double *ltr, *lts; /* local rise and set times */
  26. double *ltt;    /* local transit time */
  27. double *azr, *azs; /* local rise and set azimuths, rads e of n */
  28. double *altt;    /* local altitude at transit */
  29. int *status;    /* one or more of the RS_* defines */
  30. {
  31.     typedef struct {
  32.         Now l_now;
  33.         double l_ltr, l_lts, l_ltt, l_azr, l_azs, l_altt;
  34.         int l_hzn;
  35.         int l_status;
  36.     } Last;
  37.     /* must be in same order as the astro.h/screen.h #define's */
  38.     static Last last[NOBJ];
  39.     Last *lp;
  40.     int new;
  41.  
  42.     lp = last + p;
  43.     if (!force && same_cir (np, &lp->l_now) && same_lday (np, &lp->l_now)
  44.                         && lp->l_hzn == hzn) {
  45.         *ltr = lp->l_ltr;
  46.         *lts = lp->l_lts;
  47.         *ltt = lp->l_ltt;
  48.         *azr = lp->l_azr;
  49.         *azs = lp->l_azs;
  50.         *altt = lp->l_altt;
  51.         *status = lp->l_status;
  52.         new = 0;
  53.     } else {
  54.         *status = 0;
  55.         iterative_riset (p, np, hzn, ltr, lts, ltt, azr, azs, altt, status);
  56.         lp->l_ltr = *ltr;
  57.         lp->l_lts = *lts;
  58.         lp->l_ltt = *ltt;
  59.         lp->l_azr = *azr;
  60.         lp->l_azs = *azs;
  61.         lp->l_altt = *altt;
  62.         lp->l_status = *status;
  63.         lp->l_hzn = hzn;
  64.         lp->l_now = *np;
  65.         new = 1;
  66.     }
  67.     return (new);
  68. }
  69.  
  70. static
  71. iterative_riset (p, np, hzn, ltr, lts, ltt, azr, azs, altt, status)
  72. int p;
  73. Now *np;
  74. int hzn;
  75. double *ltr, *lts, *ltt;    /* local times of rise, set and transit */
  76. double *azr, *azs, *altt;/* local azimuths of rise, set and transit altitude */
  77. int *status;
  78. {
  79. #define    MAXPASSES    6
  80.     double lstr, lsts, lstt; /* local sidereal times of rising/setting */
  81.     double mjd0;        /* mjd estimates of rise/set event */
  82.     double lnoon;        /* mjd of local noon */
  83.     double x;        /* discarded tmp value */
  84.     Now n;            /* just used to call now_lst() */
  85.     double lst;        /* lst at local noon */
  86.     double diff, lastdiff;    /* iterative improvement to mjd0 */
  87.     int pass;
  88.     int rss;
  89.  
  90.     /* first approximation is to find rise/set times of a fixed object
  91.      * in its position at local noon.
  92.      */
  93.     lnoon = mjd_day(mjd - tz/24.0) + (12.0 + tz)/24.0; /*mjd of local noon*/
  94.     n.n_mjd = lnoon;
  95.     n.n_lng = lng;
  96.     now_lst (&n, &lst);    /* lst at local noon */
  97.     mjd0 = lnoon;
  98.     stationary_riset (p,mjd0,np,hzn,&lstr,&lsts,&lstt,&x,&x,&x,&rss);
  99.     chkrss:
  100.     switch (rss) {
  101.     case  0:  break;
  102.     case  1: *status = RS_NEVERUP; return;
  103.     case -1: *status = RS_CIRCUMPOLAR; goto transit;
  104.     default: *status = RS_ERROR; return;
  105.     }
  106.  
  107.     /* find a better approximation to the rising circumstances based on
  108.      * more passes, each using a "fixed" object at the location at
  109.      * previous approximation of the rise time.
  110.      */
  111.     lastdiff = 1000.0;
  112.     for (pass = 1; pass < MAXPASSES; pass++) {
  113.         diff = (lstr - lst)*SIDRATE; /* next guess at rise time wrt noon */
  114.         if (diff > 12.0)
  115.         diff -= 24.0*SIDRATE;    /* not tomorrow, today */
  116.         else if (diff < -12.0)
  117.         diff += 24.0*SIDRATE;    /* not yesterday, today */
  118.         mjd0 = lnoon + diff/24.0;    /* next guess at mjd of rise */
  119.         stationary_riset (p,mjd0,np,hzn,&lstr,&x,&x,azr,&x,&x,&rss);
  120.         if (rss != 0) goto chkrss;
  121.         if (fabs (diff - lastdiff) < 0.25/60.0)
  122.         break;             /* converged to better than 15 secs */
  123.         lastdiff = diff;
  124.     }
  125.     if (pass == MAXPASSES)
  126.         *status |= RS_NORISE;    /* didn't converge - no rise today */
  127.     else {
  128.         *ltr = 12.0 + diff;
  129.         if (p != MOON &&
  130.             (*ltr <= 24.0*(1.0-SIDRATE) || *ltr >= 24.0*SIDRATE))
  131.         *status |= RS_2RISES;
  132.     }
  133.  
  134.     /* find a better approximation to the setting circumstances based on
  135.      * more passes, each using a "fixed" object at the location at
  136.      * previous approximation of the set time.
  137.      */
  138.     lastdiff = 1000.0;
  139.     for (pass = 1; pass < MAXPASSES; pass++) {
  140.         diff = (lsts - lst)*SIDRATE; /* next guess at set time wrt noon */
  141.         if (diff > 12.0)
  142.         diff -= 24.0*SIDRATE;    /* not tomorrow, today */
  143.         else if (diff < -12.0)
  144.         diff += 24.0*SIDRATE;    /* not yesterday, today */
  145.         mjd0 = lnoon + diff/24.0;    /* next guess at mjd of set */
  146.         stationary_riset (p,mjd0,np,hzn,&x,&lsts,&x,&x,azs,&x,&rss);
  147.         if (rss != 0) goto chkrss;
  148.         if (fabs (diff - lastdiff) < 0.25/60.0)
  149.         break;             /* converged to better than 15 secs */
  150.         lastdiff = diff;
  151.     }
  152.     if (pass == MAXPASSES)
  153.         *status |= RS_NOSET;    /* didn't converge - no set today */
  154.     else {
  155.         *lts = 12.0 + diff;
  156.         if (p != MOON &&
  157.             (*lts <= 24.0*(1.0-SIDRATE) || *lts >= 24.0*SIDRATE))
  158.         *status |= RS_2SETS;
  159.     }
  160.  
  161.     transit:
  162.     /* find a better approximation to the transit circumstances based on
  163.      * more passes, each using a "fixed" object at the location at
  164.      * previous approximation of the transit time.
  165.      */
  166.     lastdiff = 1000.0;
  167.     for (pass = 1; pass < MAXPASSES; pass++) {
  168.         diff = (lstt - lst)*SIDRATE; /*next guess at transit time wrt noon*/
  169.         if (diff > 12.0)
  170.         diff -= 24.0*SIDRATE;    /* not tomorrow, today */
  171.         else if (diff < -12.0)
  172.         diff += 24.0*SIDRATE;    /* not yesterday, today */
  173.         mjd0 = lnoon + diff/24.0;    /* next guess at mjd of set */
  174.         stationary_riset (p,mjd0,np,hzn,&x,&x,&lstt,&x,&x,altt,&rss);
  175.         if (fabs (diff - lastdiff) < 0.25/60.0)
  176.         break;             /* converged to better than 15 secs */
  177.         lastdiff = diff;
  178.     }
  179.     if (pass == MAXPASSES)
  180.         *status |= RS_NOTRANS;    /* didn't converge - no transit today */
  181.     else {
  182.         *ltt = 12.0 + diff;
  183.         if (p != MOON &&
  184.             (*ltt <= 24.0*(1.0-SIDRATE) || *ltt >= 24.0*SIDRATE))
  185.         *status |= RS_2TRANS;
  186.     }
  187. }
  188.  
  189. static
  190. stationary_riset (p, mjd0, np, hzn, lstr, lsts, lstt, azr, azs, altt, status)
  191. int p;
  192. double mjd0;
  193. Now *np;
  194. int hzn;
  195. double *lstr, *lsts, *lstt;
  196. double *azr, *azs, *altt;
  197. int *status;
  198. {
  199.     extern void bye();
  200.     double dis;
  201.     Now n;
  202.     Sky s;
  203.  
  204.     switch (hzn) {
  205.     case STDHZN:
  206.         /* nominal atmospheric refraction.
  207.          * then add nominal moon or sun semi-diameter, as appropriate.
  208.          */
  209.         dis = STDREF;
  210.         break;
  211.     case TWILIGHT:
  212.         if (p != SUN) {
  213.         f_msg ("Non-sun twilight bug!");
  214.         bye();
  215.         }
  216.         dis = TWIREF;
  217.         break;
  218.     default:
  219.         /* horizon displacement is refraction plus object's semi-diameter */
  220.         unrefract (pressure, temp, 0.0, &dis);
  221.         dis = -dis;
  222.         n = *np;
  223.         n.n_mjd = mjd0;
  224.         (void) body_cir (p, 0.0, &n, &s);
  225.         dis += degrad(s.s_size/3600./2.0);
  226.     }
  227.  
  228.     switch (p) {
  229.     case SUN:
  230.         if (hzn == STDHZN)
  231.         dis += degrad (32./60./2.);
  232.         fixedsunriset (mjd0,np,dis,lstr,lsts,lstt,azr,azs,altt,status);
  233.         break;
  234.     case MOON:
  235.         if (hzn == STDHZN)
  236.         dis += degrad (32./60./2.);
  237.         fixedmoonriset (mjd0,np,dis,lstr,lsts,lstt,azr,azs,altt,status);
  238.         break;
  239.     default:
  240.         if (hzn == STDHZN) {
  241.         /* assume planets have negligible diameters.
  242.          * also, need to set s if hzn was STDHZN.
  243.          */
  244.         n = *np;
  245.         n.n_mjd = mjd0;
  246.         (void) body_cir (p, 0.0, &n, &s);
  247.         }
  248.         riset (s.s_ra, s.s_dec, lat, dis, lstr, lsts, azr, azs, status);
  249.         transit (s.s_ra, s.s_dec, np, lstt, altt);
  250.     }
  251. }
  252.  
  253. /* find the local rise/set sidereal times and azimuths of an object fixed
  254.  * at the ra/dec of the sun on the given mjd time as seen from lat.
  255.  * times are for the upper limb. dis should account for refraction and
  256.  *   sun's semi-diameter; we add corrections for nutation and aberration.
  257.  */
  258. static
  259. fixedsunriset (mjd0, np, dis, lstr, lsts, lstt, azr, azs, altt, status)
  260. double mjd0;
  261. Now *np;
  262. double dis;
  263. double *lstr, *lsts, *lstt;
  264. double *azr, *azs, *altt;
  265. int *status;
  266. {
  267.     double lsn, rsn;
  268.     double deps, dpsi;
  269.     double r, d;
  270.  
  271.     /* find ecliptic position of sun at mjd */
  272.     sunpos (mjd0, &lsn, &rsn);
  273.  
  274.     /* allow for nutation and 20.4" light travel time */
  275.     nutation (mjd0, &deps, &dpsi);
  276.         lsn += dpsi;
  277.         lsn -= degrad(20.4/3600.0);
  278.  
  279.     /* convert ecliptic to equatorial coords */
  280.     ecl_eq (mjd0, 0.0, lsn, &r, &d);
  281.  
  282.     /* find circumstances for given horizon displacement */
  283.     riset (r, d, lat, dis, lstr, lsts, azr, azs, status);
  284.     transit (r, d, np, lstt, altt);
  285. }
  286.  
  287. /* find the local sidereal rise/set times and azimuths of an object fixed
  288.  * at the ra/dec of the moon on the given mjd time as seen from lat.
  289.  * times are for the upper limb. dis should account for refraction and moon's
  290.  *    semi-diameter; we add corrections for parallax, and nutation.
  291.  * accuracy is to nearest minute.
  292.  */
  293. static
  294. fixedmoonriset (mjd0, np, dis, lstr, lsts, lstt, azr, azs, altt, status)
  295. double mjd0;
  296. Now *np;
  297. double dis;
  298. double *lstr, *lsts, *lstt;
  299. double *azr, *azs, *altt;
  300. int *status;
  301. {
  302.     double lam, bet, hp;
  303.     double deps, dpsi;
  304.     double r, d;
  305.     double ha;
  306.  
  307.     /* find geocentric ecliptic location and equatorial horizontal parallax
  308.      * of moon at mjd.
  309.      */
  310.     moon (mjd0, &lam, &bet, &hp);
  311.  
  312.     /* allow for nutation */
  313.     nutation (mjd0, &deps, &dpsi);
  314.         lam += dpsi;
  315.  
  316.     /* convert ecliptic to equatorial coords */
  317.     ecl_eq (mjd0, bet, lam, &r, &d);
  318.  
  319.     /* find local sidereal times of rise/set times/azimuths for given
  320.      * equatorial coords, allowing for
  321.      * .27249*sin(hp)   parallax (hp is radius of earth from moon;
  322.      *                  .27249 is radius of moon from earth where
  323.      *                  the ratio is the dia_moon/dia_earth).
  324.      * hp               nominal angular earth radius. subtract because
  325.      *                  tangential line-of-sight makes moon appear lower
  326.      *                  as move out from center of earth.
  327.      */
  328.     dis += .27249*sin(hp) - hp;
  329.     riset (r, d, lat, dis, lstr, lsts, azr, azs, status);
  330.  
  331.     /* account for parallax on the meridian (0 hour-angle).
  332.      * TODO: is this really right? other stuff too? better way?? help!
  333.      */
  334.     ta_par (0.0, d, lat, height, hp, &ha, &d);
  335.     transit (r+ha, d, np, lstt, altt);
  336. }
  337.  
  338. /* find when and how hi object at (r,d) is when it transits. */
  339. static
  340. transit (r, d, np, lstt, altt)
  341. double r, d;    /* ra and dec, rads */
  342. Now *np;    /* for refraction info */
  343. double *lstt;    /* local sidereal time of transit */
  344. double *altt;    /* local, refracted, altitude at time of transit */
  345. {
  346.     *lstt = radhr(r);
  347.     *altt = PI/2 - lat + d;
  348.     if (*altt > PI/2)
  349.         *altt = PI - *altt;
  350.     refract (pressure, temp, *altt, altt);
  351. }
  352.