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