home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / misc / volume14 / ephem-4.21 / part03 / circum.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-08-30  |  10.2 KB  |  373 lines
  1. /* fill in a Sky struct with all we know about each object.
  2.  *(the user defined objects are in obj.c)
  3.  */
  4.  
  5. #include <stdio.h>
  6. #include <math.h>
  7. #include "astro.h"
  8. #include "circum.h"
  9. #include "screen.h"    /* just for SUN and MOON */
  10.  
  11. /* find body p's circumstances now.
  12.  * to save some time the caller may specify a desired accuracy, in arc seconds.
  13.  * if, based on its mean motion, it would not have moved this much since the
  14.  * last time we were called we only recompute altitude and azimuth and avoid
  15.  * recomputing the planet's heliocentric position. use 0.0 for best possible.
  16.  * we always recompute the user-defined objects' position regardless.
  17.  * return 0 if only alt/az changes, else 1 if all other stuff updated too.
  18.  * TODO: beware of fast retrograde motions.
  19.  */
  20. body_cir (p, as, np, sp)
  21. int p;
  22. double as;
  23. Now *np;
  24. Sky *sp;
  25. {
  26.     typedef struct {
  27.         double l_dpas;    /* mean days per arc second */
  28.         Now l_now;        /* when l_sky was found */
  29.         double l_ra, l_dec;    /* the eod, ie, unprecessed, ra/dec values */
  30.         Sky l_sky;
  31.     } Last;
  32.     /* must be in same order as the astro.h object #define's */
  33.     static Last last[8] = {
  34.         {.000068, {NOMJD}},
  35.         {.00017, {NOMJD}},
  36.         {.00053, {NOMJD}},
  37.         {.0034, {NOMJD}},
  38.         {.0092, {NOMJD}},
  39.         {.027, {NOMJD}},
  40.         {.046, {NOMJD}},
  41.         {.069, {NOMJD}}
  42.     };
  43.     Last objxlast, objylast;
  44.     double lst, alt, az;
  45.     double ehp, ha, dec;    /* ehp: angular dia of earth from body */
  46.     Last *lp;
  47.     int new;
  48.  
  49.     switch (p) {
  50.     case SUN: return (sun_cir (as, np, sp));
  51.     case MOON: return (moon_cir (as, np, sp));
  52.     case OBJX: lp = &objxlast; break;
  53.     case OBJY: lp = &objylast; break;
  54.     default: lp = last + p; break;
  55.     }
  56.  
  57.     /* if less than l_every days from last time for this planet
  58.      * just redo alt/az.
  59.      * ALWAYS redo objects x and y.
  60.      */
  61.     if (p != OBJX && p != OBJY && same_cir (np, &lp->l_now)
  62.               && about_now (np, &lp->l_now, as*lp->l_dpas)) {
  63.         *sp = lp->l_sky;
  64.         new = 0;
  65.     } else {
  66.         double lpd0, psi0;    /* heliocentric ecliptic long and lat */
  67.         double rp0;        /* dist from sun */
  68.         double rho0;    /* dist from earth */
  69.         double lam, bet;    /* geocentric ecliptic long and lat */
  70.         double dia, mag;    /* angular diameter at 1 AU and magnitude */
  71.         double lsn, rsn;    /* true geoc lng of sun, dist from sn to earth*/
  72.         double el;    /* elongation */
  73.         double f;   /* phase from earth */
  74.  
  75.         lp->l_now = *np;
  76.         sunpos (mjd, &lsn, &rsn);
  77.         if (p == OBJX || p == OBJY)
  78.         obj_cir(mjd,p,&lpd0,&psi0,&rp0,&rho0,&lam,&bet,&sp->s_mag);
  79.         else {
  80.         double deps, dpsi;
  81.         double a;
  82.         plans(mjd, p, &lpd0, &psi0, &rp0, &rho0, &lam, &bet, &dia,&mag);
  83.         nutation (mjd, &deps, &dpsi);    /* correct for nutation */
  84.         lam += dpsi;
  85.         a = lsn-lam;            /* and 20.4" aberation */
  86.         lam -= degrad(20.4/3600)*cos(a)/cos(bet);
  87.         bet -= degrad(20.4/3600)*sin(a)*sin(bet);
  88.         }
  89.  
  90.         ecl_eq (mjd, bet, lam, &lp->l_ra, &lp->l_dec);
  91.  
  92.         sp->s_ra = lp->l_ra;
  93.         sp->s_dec = lp->l_dec;
  94.         if (epoch != EOD)
  95.         precess (mjd, epoch, &sp->s_ra, &sp->s_dec);
  96.         sp->s_edist = rho0;
  97.         sp->s_sdist = rp0;
  98.         elongation (lam, bet, lsn, &el);
  99.         el = raddeg(el);
  100.         sp->s_elong = el;
  101.         f = (rp0 > 0.0)
  102.         ? 0.25 * (((rp0+rho0)*(rp0+rho0) - rsn*rsn)/(rp0*rho0)) : 0.0;
  103.         sp->s_phase = f*100.0; /* percent */
  104.         if (p != OBJX && p != OBJY) {
  105.         sp->s_size = dia/rho0;
  106.         sp->s_mag = mag + 5.0*log(rp0*rho0/sqrt(f))/log(10.0);
  107.         }
  108.         sp->s_hlong = lpd0;
  109.         sp->s_hlat = psi0;
  110.         new = 1;
  111.     }
  112.  
  113.     /* alt, az; correct for parallax and refraction; use eod ra/dec */
  114.     now_lst (np, &lst);
  115.     ha = hrrad(lst) - lp->l_ra;
  116.     if (sp->s_edist > 0.0) {
  117.         ehp = (2.0*6378.0/146.0e6) / sp->s_edist;
  118.         ta_par (ha, lp->l_dec, lat, height, ehp, &ha, &dec);
  119.     } else
  120.         dec = lp->l_dec;
  121.     hadec_aa (lat, ha, dec, &alt, &az);
  122.     refract (pressure, temp, alt, &alt);
  123.     sp->s_alt = alt;
  124.     sp->s_az = az;
  125.     lp->l_sky = *sp;
  126.     return (new);
  127. }
  128.  
  129. /* find local times when sun is 18 degrees below horizon.
  130.  * return 0 if just returned same stuff as previous call, else 1 if new.
  131.  */
  132. twilight_cir (np, dawn, dusk, status)
  133. Now *np;
  134. double *dawn, *dusk;
  135. int *status;
  136. {
  137.     static Now last_now = {NOMJD};
  138.     static double last_dawn, last_dusk;
  139.     static int last_status;
  140.     int new;
  141.  
  142.     if (same_cir (np, &last_now) && same_lday (np, &last_now)) {
  143.         *dawn = last_dawn;
  144.         *dusk = last_dusk;
  145.         *status = last_status;
  146.         new = 0;
  147.     } else {
  148.         double x;
  149.         (void) riset_cir (SUN,np,0,TWILIGHT,dawn,dusk,&x,&x,&x,&x,status);
  150.         last_dawn = *dawn;
  151.         last_dusk = *dusk;
  152.         last_status = *status;
  153.         last_now = *np;
  154.         new = 1;
  155.     }
  156.     return (new);
  157. }
  158.  
  159. /* find sun's circumstances now.
  160.  * as is the desired accuracy, in arc seconds; use 0.0 for best possible.
  161.  * return 0 if only alt/az changes, else 1 if all other stuff updated too.
  162.  */
  163. sun_cir (as, np, sp)
  164. double as;
  165. Now *np;
  166. Sky *sp;
  167. {
  168.     static Sky last_sky;
  169.     static Now last_now = {NOMJD};
  170.     static double last_ra, last_dec;    /* unprecessed ra/dec */
  171.     double lst, alt, az;
  172.     double ehp, ha, dec;    /* ehp: angular dia of earth from body */
  173.     int new;
  174.  
  175.     if (same_cir (np, &last_now) && about_now (np, &last_now, as*.00028)) {
  176.         *sp = last_sky;
  177.         new = 0;
  178.     } else {
  179.         double lsn, rsn;
  180.         double deps, dpsi;
  181.  
  182.         last_now = *np;
  183.         sunpos (mjd, &lsn, &rsn);        /* sun's true ecliptic long
  184.                          * and dist
  185.                          */
  186.         nutation (mjd, &deps, &dpsi);    /* correct for nutation */
  187.         lsn += dpsi;
  188.         lsn -= degrad(20.4/3600);        /* and light travel time */
  189.  
  190.         sp->s_edist = rsn;
  191.         sp->s_sdist = 0.0;
  192.         sp->s_elong = 0.0;
  193.         sp->s_size = raddeg(4.65242e-3/rsn)*3600*2;
  194.         sp->s_mag = -26.8;
  195.         sp->s_hlong = lsn-PI;    /* geo- to helio- centric */
  196.         range (&sp->s_hlong, 2*PI);
  197.         sp->s_hlat = 0.0;
  198.  
  199.         ecl_eq (mjd, 0.0, lsn, &last_ra, &last_dec);
  200.         sp->s_ra = last_ra;
  201.         sp->s_dec = last_dec;
  202.         if (epoch != EOD)
  203.         precess (mjd, epoch, &sp->s_ra, &sp->s_dec);
  204.         new = 1;
  205.     }
  206.  
  207.     now_lst (np, &lst);
  208.     ha = hrrad(lst) - last_ra;
  209.     ehp = (2.0 * 6378.0 / 146.0e6) / sp->s_edist;
  210.     ta_par (ha, last_dec, lat, height, ehp, &ha, &dec);
  211.     hadec_aa (lat, ha, dec, &alt, &az);
  212.     refract (pressure, temp, alt, &alt);
  213.     sp->s_alt = alt;
  214.     sp->s_az = az;
  215.     last_sky = *sp;
  216.     return (new);
  217. }
  218.  
  219. /* find moon's circumstances now.
  220.  * as is the desired accuracy, in arc seconds; use 0.0 for best possible.
  221.  * return 0 if only alt/az changes, else 1 if all other stuff updated too.
  222.  */
  223. moon_cir (as, np, sp)
  224. double as;
  225. Now *np;
  226. Sky *sp;
  227. {
  228.     static Sky last_sky;
  229.     static Now last_now = {NOMJD};
  230.     static double ehp;
  231.     static double last_ra, last_dec;    /* unprecessed */
  232.     double lst, alt, az;
  233.     double ha, dec;
  234.     int new;
  235.  
  236.     if (same_cir (np, &last_now) && about_now (np, &last_now, as*.000021)) {
  237.         *sp = last_sky;
  238.         new = 0;
  239.     } else {
  240.         double lam, bet;
  241.         double deps, dpsi;
  242.         double lsn, rsn;    /* sun long in rads, earth-sun dist in au */
  243.         double edistau;    /* earth-moon dist, in au */
  244.         double el;        /* elongation, rads east */
  245.  
  246.         last_now = *np;
  247.         moon (mjd, &lam, &bet, &ehp);    /* moon's true ecliptic loc */
  248.         nutation (mjd, &deps, &dpsi);    /* correct for nutation */
  249.         lam += dpsi;
  250.         range (&lam, 2*PI);
  251.  
  252.         sp->s_edist = 6378.14/sin(ehp);    /* earth-moon dist, want km */
  253.         sp->s_size = 3600*31.22512*sin(ehp);/* moon angular dia, seconds */
  254.         sp->s_hlong = lam;            /* save geo in helio fields */
  255.         sp->s_hlat = bet;
  256.  
  257.         ecl_eq (mjd, bet, lam, &last_ra, &last_dec);
  258.         sp->s_ra = last_ra;
  259.         sp->s_dec = last_dec;
  260.         if (epoch != EOD)
  261.         precess (mjd, epoch, &sp->s_ra, &sp->s_dec);
  262.  
  263.         sunpos (mjd, &lsn, &rsn);
  264.         range (&lsn, 2*PI);
  265.         elongation (lam, bet, lsn, &el);
  266.  
  267.         /* solve triangle of earth, sun, and elongation for moon-sun dist */
  268.         edistau = sp->s_edist/1.495979e8; /* km -> au */
  269.         sp->s_sdist =
  270.         sqrt (edistau*edistau + rsn*rsn - 2.0*edistau*rsn*cos(el));
  271.  
  272.         /* TODO: improve mag; this is based on a flat moon model. */
  273.         sp->s_mag = -12.7 + 2.5*(log10(PI) - log10(PI/2*(1+1.e-6-cos(el))));
  274.  
  275.         sp->s_elong = raddeg(el);    /* want degrees */
  276.         sp->s_phase = fabs(el)/PI*100.0;    /* want non-negative % */
  277.         new = 1;
  278.     }
  279.  
  280.     /* show topocentric alt/az by correcting ra/dec for parallax 
  281.      * as well as refraction.
  282.      */
  283.     now_lst (np, &lst);
  284.     ha = hrrad(lst) - last_ra;
  285.     ta_par (ha, last_dec, lat, height, ehp, &ha, &dec);
  286.     hadec_aa (lat, ha, dec, &alt, &az);
  287.     refract (pressure, temp, alt, &alt);
  288.     sp->s_alt = alt;
  289.     sp->s_az = az;
  290.     last_sky = *sp;
  291.     return (new);
  292. }
  293.  
  294. /* given geocentric ecliptic longitude and latitude, lam and bet, of some object
  295.  * and the longitude of the sun, lsn, find the elongation, el. this is the
  296.  * actual angular separation of the object from the sun, not just the difference
  297.  * in the longitude. the sign, however, IS set simply as a test on longitude
  298.  * such that el will be >0 for an evening object <0 for a morning object.
  299.  * to understand the test for el sign, draw a graph with lam going from 0-2*PI
  300.  *   down the vertical axis, lsn going from 0-2*PI across the hor axis. then
  301.  *   define the diagonal regions bounded by the lines lam=lsn+PI, lam=lsn and
  302.  *   lam=lsn-PI. the "morning" regions are any values to the lower left of the
  303.  *   first line and bounded within the second pair of lines.
  304.  * all angles in radians.
  305.  */
  306. elongation (lam, bet, lsn, el)
  307. double lam, bet, lsn;
  308. double *el;
  309. {
  310.     *el = acos(cos(bet)*cos(lam-lsn));
  311.     if (lam>lsn+PI || lam>lsn-PI && lam<lsn) *el = - *el;
  312. }
  313.  
  314. /* return whether the two Nows are for the same observing circumstances. */
  315. same_cir (n1, n2)
  316. register Now *n1, *n2;
  317. {
  318.     return (n1->n_lat == n2->n_lat
  319.         && n1->n_lng == n2->n_lng
  320.         && n1->n_temp == n2->n_temp
  321.         && n1->n_pressure == n2->n_pressure
  322.         && n1->n_height == n2->n_height
  323.         && n1->n_tz == n2->n_tz
  324.         && n1->n_epoch == n2->n_epoch);
  325. }
  326.  
  327. /* return whether the two Nows are for the same LOCAL day */
  328. same_lday (n1, n2)
  329. Now *n1, *n2;
  330. {
  331.     return (mjd_day(n1->n_mjd - n1->n_tz/24.0) ==
  332.                mjd_day(n2->n_mjd - n2->n_tz/24.0)); 
  333. }
  334.  
  335. /* return whether the mjd of the two Nows are within dt */
  336. static
  337. about_now (n1, n2, dt)
  338. Now *n1, *n2;
  339. double dt;
  340. {
  341.     return (fabs (n1->n_mjd - n2->n_mjd) <= dt/2.0);
  342. }
  343.  
  344. now_lst (np, lst)
  345. Now *np;
  346. double *lst;
  347. {
  348.     utc_gst (mjd_day(mjd), mjd_hr(mjd), lst);
  349.     *lst += radhr(lng);
  350.     range (lst, 24.0);
  351. }
  352.  
  353. /* round a time in days, *t, to the nearest second, IN PLACE. */
  354. rnd_second (t)
  355. double *t;
  356. {
  357.     *t = floor(*t*SPD+0.5)/SPD;
  358. }
  359.     
  360. double
  361. mjd_day(jd)
  362. double jd;
  363. {
  364.     return (floor(jd-0.5)+0.5);
  365. }
  366.  
  367. double
  368. mjd_hr(jd)
  369. double jd;
  370. {
  371.     return ((jd-mjd_day(jd))*24.0);
  372. }
  373.