home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 October / usenetsourcesnewsgroupsinfomagicoctober1994disk2.iso / misc / volume9 / ephem2 / part02 < prev    next >
Text File  |  1989-11-27  |  56KB  |  2,265 lines
  1. Newsgroups: comp.sources.misc
  2. From: allbery@uunet.UU.NET (Brandon S. Allbery - comp.sources.misc)
  3. Subject: v09i039: ephem, v4.8, 2 of 5
  4. Reply-To: ecd@umn-cs.cs.umn.edu@ncs-med.UUCP (Elwood C. Downey)
  5.  
  6. Posting-number: Volume 9, Issue 39
  7. Submitted-by: ecd@umn-cs.cs.umn.edu@ncs-med.UUCP (Elwood C. Downey)
  8. Archive-name: ephem2/part02
  9.  
  10. # This is a "shell archive" file; run it with sh.
  11. # This is file 2.
  12. echo x screen.h
  13. cat > screen.h << 'xXx'
  14. /* screen layout details
  15.  *
  16.  * it looks better if the fields are drawn in some nice order so it you
  17.  * rearrange the fields, check the menu printing functions.
  18.  */
  19.  
  20. /* size of screen */
  21. #define    NR    24
  22. #define    NC    80
  23.  
  24. #define    ASPECT    (4./3.)    /* screen width to height dimensions ratio */
  25.  
  26. #define    GAP    6    /* gap between field name and value */
  27.  
  28. #define    COL1        1
  29. #define    COL2        27
  30. #define    COL3        44
  31. #define    COL4        61    /* calendar */
  32.  
  33. #define    R_PROMPT    1    /* prompt row */
  34. #define    C_PROMPT    COL1
  35.  
  36. #define    R_NEWCIR    2
  37. #define    C_NEWCIR    ((NC-17)/2) /* 17 is length of the message */
  38.  
  39. #define    R_TOP        3    /* first row of top menu items */
  40.  
  41. #define    R_TZN    (R_TOP+0)
  42. #define    C_TZN    COL1
  43. #define    R_LT    R_TZN
  44. #define    C_LT    (C_TZN+GAP-2)
  45. #define    R_LD    R_TZN
  46. #define    C_LD    (C_TZN+13)
  47.  
  48. #define    R_UT    (R_TOP+1)
  49. #define    C_UT    COL1
  50. #define    C_UTV    (C_UT+GAP-2)
  51. #define    R_UD    R_UT
  52. #define    C_UD    (C_UT+13)
  53.  
  54. #define    R_JD    (R_TOP+2)
  55. #define    C_JD    COL1
  56. #define    C_JDV    (C_JD+GAP+3)
  57.  
  58. #define    R_LST    (R_TOP)
  59. #define    C_LST    COL2
  60. #define    C_LSTV    (C_LST+GAP)
  61.  
  62. #define    R_LAT    (R_TOP+0)
  63. #define    C_LAT    COL3
  64. #define    C_LATV    (C_LAT+4)
  65.  
  66. #define    R_DAWN    (R_TOP+2)
  67. #define    C_DAWN    COL2
  68. #define    C_DAWNV    (C_DAWN+GAP+3)
  69.  
  70. #define    R_STPSZ    (R_TOP+7)
  71. #define    C_STPSZ    COL2
  72. #define    C_STPSZV (C_STPSZ+GAP-1)
  73.  
  74. #define    R_HEIGHT (R_TOP+2)
  75. #define    C_HEIGHT COL3
  76. #define    C_HEIGHTV (C_HEIGHT+GAP)
  77.  
  78. #define    R_PRES    (R_TOP+4)
  79. #define    C_PRES    COL3
  80. #define    C_PRESV    (C_PRES+GAP)
  81.  
  82. #define    R_WATCH    (R_TOP+4)
  83. #define    C_WATCH    COL1
  84.  
  85. #define    R_SRCH    (R_TOP+5)
  86. #define    C_SRCH    COL1
  87. #define    C_SRCHV    (C_SRCH+16)
  88.  
  89. #define    R_PLOT    (R_TOP+6)
  90. #define    C_PLOT    (COL1)
  91. #define    C_PLOTV (C_PLOT+20)
  92.  
  93. #define    R_ALTM    (R_TOP+7)
  94. #define    C_ALTM    COL1
  95. #define    C_ALTMV    (C_ALTM+10)
  96.  
  97. #define    R_TZONE    (R_TOP+5)
  98. #define    C_TZONE    COL3
  99. #define    C_TZONEV (C_TZONE+GAP-1)
  100.  
  101. #define    R_LONG    (R_TOP+1)
  102. #define    C_LONG    COL3
  103. #define    C_LONGV    (C_LONG+4)
  104.  
  105. #define    R_DUSK    (R_TOP+3)
  106. #define    C_DUSK    COL2
  107. #define    C_DUSKV    (C_DUSK+GAP+3)
  108.  
  109. #define    R_NSTEP (R_TOP+6)
  110. #define    C_NSTEP    COL2
  111. #define    C_NSTEPV (C_NSTEP+GAP)
  112.  
  113. #define    R_TEMP    (R_TOP+3)
  114. #define    C_TEMP    COL3
  115. #define    C_TEMPV    (C_TEMP+GAP)
  116.  
  117. #define    R_EPOCH        (R_TOP+6)
  118. #define    C_EPOCH        COL3
  119. #define    C_EPOCHV    (C_EPOCH+GAP)
  120.  
  121. #define    R_MNUDEP    (R_TOP+6)
  122. #define    C_MNUDEP    COL3
  123. #define    C_MNUDEPV    (C_EPOCH+GAP)
  124.  
  125. #define    R_LON    (R_TOP+4)
  126. #define    C_LON    COL2
  127. #define    C_LONV    (C_LON+GAP+3)
  128.  
  129. #define    R_CAL    R_TOP
  130. #define    C_CAL   COL4
  131.  
  132. /* menu 1 info table */
  133. #define    R_PLANTAB    (R_TOP+9)
  134. #define    R_SUN        (R_PLANTAB+2)
  135. #define    R_MOON        (R_PLANTAB+3)
  136. #define    R_MERCURY    (R_PLANTAB+4)
  137. #define    R_VENUS        (R_PLANTAB+5)
  138. #define    R_MARS        (R_PLANTAB+6)
  139. #define    R_JUPITER    (R_PLANTAB+7)
  140. #define    R_SATURN    (R_PLANTAB+8)
  141. #define    R_URANUS    (R_PLANTAB+9)
  142. #define    R_NEPTUNE    (R_PLANTAB+10)
  143. #define    R_PLUTO        (R_PLANTAB+11)
  144. #define    R_OBJX        (R_PLANTAB+12)
  145. #define    C_OBJ        1
  146. #define    C_RA        4
  147. #define    C_DEC        12
  148. #define    C_AZ        19
  149. #define    C_ALT        26
  150. #define    C_HLONG        33
  151. #define    C_HLAT        40
  152. #define    C_EDIST        47
  153. #define C_SDIST     54
  154. #define    C_ELONG        61
  155. #define    C_SIZE        68
  156. #define    C_MAG        73
  157. #define    C_PHASE        78
  158.  
  159. /* menu 2 screen items */
  160. #define    C_RISETM    10
  161. #define    C_RISEAZ    18
  162. #define    C_TRANSTM    31
  163. #define    C_TRANSALT    39
  164. #define    C_SETTM        52
  165. #define    C_SETAZ        60
  166. #define    C_TUP        72
  167.  
  168. /* menu 3 items */
  169. #define    C_SUN        4
  170. #define    C_MOON        11
  171. #define    C_MERCURY    18
  172. #define    C_VENUS        25
  173. #define    C_MARS        32
  174. #define    C_JUPITER    39
  175. #define    C_SATURN    46
  176. #define    C_URANUS    53
  177. #define    C_NEPTUNE    60
  178. #define    C_PLUTO        67
  179. #define    C_OBJX        74
  180.  
  181. #define    MAXOBJXNM    2    /* object x's name. does not include trail 0 */
  182.  
  183. #define    PW    (NC-C_PROMPT+1)    /* total prompt line width */
  184.  
  185. /* macros to pack a row/col and menu selection flags all into an int.
  186.  * (use this rather than a structure because we can compare them so easily.
  187.  * could use bit fields and a union, but then can't init them or use switch.)
  188.  * bit field defs: [15..14]=menu  [13..12]=flags  [11..7]=row  [6..0]=column.
  189.  * see sel_fld.c.
  190.  * F_MNUX also used in main to manage which bottom menu is up.
  191.  */
  192. #define    F_CHG        (1<<12)    /* field may be picked for changing */
  193. #define    F_PLT        (1<<13)    /* field may be picked for plotting */
  194. #define    F_MMNU        (0<<14)    /* field is on main menu */
  195. #define    F_MNU1        (1<<14)    /* field is on menu 1 */
  196. #define    F_MNU2         (2<<14)    /* field is on menu 2 */
  197. #define    F_MNU3        (3<<14)    /* field is on menu 3 */
  198. #define    rcfpack(r,c,f)    ((f) | ((r) << 7) | (c))
  199. #define    unpackr(p)    (((p) >> 7) & 0x1f)
  200. #define    unpackc(p)    ((p) & 0x7f)
  201. #define    unpackrc(p)    ((p) & 0xfff)
  202. #define    tstpackf(p,f)    (((p) & ((f)&0x3000)) && \
  203.             (((p)&0xc000) == ((f)&0xc000) || ((p)&0xc000) == 0))
  204.  
  205. /* additions to the planet defines from astro.h.
  206.  * must not conflict, and must fit in range 0..15.
  207.  */
  208. #define    SUN    (PLUTO+1)
  209. #define    MOON    (PLUTO+2)
  210. #define    OBJX    (PLUTO+3)    /* the user-defined object */
  211. #define    NOBJ    (OBJX+1)    /* total number of objects */
  212.  
  213. #define    cntrl(x)    ((x) & 037)
  214. #define    QUIT        cntrl('d')    /* char to exit program */
  215. #define    HELP        '?'        /* char to give help message */
  216. #define    REDRAW        cntrl('l')    /* char to redraw (like vi) */
  217. #define    VERSION        cntrl('v')    /* char to display version number */
  218. #define    END        'q'        /* char to quit current mode */
  219. xXx
  220. echo x aa_hadec.c
  221. cat > aa_hadec.c << 'xXx'
  222. #include <stdio.h>
  223. #include <math.h>
  224. #include "astro.h"
  225.  
  226. /* given latitude (n+, radians), lat, altitude (up+, radians), alt, and
  227.  * azimuth (angle round to the east from north+, radians),
  228.  * return hour angle (radians), ha, and declination (radians), dec.
  229.  */
  230. aa_hadec (lat, alt, az, ha, dec)
  231. double lat;
  232. double alt, az;
  233. double *ha, *dec;
  234. {
  235.     aaha_aux (lat, az, alt, ha, dec);
  236. }
  237.  
  238. /* given latitude (n+, radians), lat, hour angle (radians), ha, and declination
  239.  * (radians), dec,
  240.  * return altitude (up+, radians), alt, and
  241.  * azimuth (angle round to the east from north+, radians),
  242.  */
  243. hadec_aa (lat, ha, dec, alt, az)
  244. double lat;
  245. double ha, dec;
  246. double *alt, *az;
  247. {
  248.     aaha_aux (lat, ha, dec, az, alt);
  249. }
  250.  
  251. /* the actual formula is the same for both transformation directions so
  252.  * do it here once for each way.
  253.  * N.B. all arguments are in radians.
  254.  */
  255. static
  256. aaha_aux (lat, x, y, p, q)
  257. double lat;
  258. double x, y;
  259. double *p, *q;
  260. {
  261.     static double lastlat = -1000.;
  262.     static double sinlastlat, coslastlat;
  263.     double sy, cy;
  264.     double sx, cx;
  265.     double sq, cq;
  266.     double a;
  267.     double cp;
  268.  
  269.     /* latitude doesn't change much, so try to reuse the sin and cos evals.
  270.      */
  271.     if (lat != lastlat) {
  272.         sinlastlat = sin (lat);
  273.         coslastlat = cos (lat);
  274.         lastlat = lat;
  275.     }
  276.  
  277.     sy = sin (y);
  278.     cy = cos (y);
  279.     sx = sin (x);
  280.     cx = cos (x);
  281.  
  282. /* define GOODATAN2 if atan2 returns full range -PI through +PI.
  283.  */
  284. #ifdef GOODATAN2
  285.     *q = asin ((sy*sinlastlat) + (cy*coslastlat*cx));
  286.     *p = atan2 (-cy*sx, -cy*cx*sinlastlat + sy*coslastlat);
  287. #else
  288. #define    EPS    (1e-20)
  289.     sq = (sy*sinlastlat) + (cy*coslastlat*cx);
  290.     *q = asin (sq);
  291.     cq = cos (*q);
  292.     a = coslastlat*cq;
  293.     if (a > -EPS && a < EPS)
  294.         a = a < 0 ? -EPS : EPS; /* avoid / 0 */
  295.     cp = (sy - (sinlastlat*sq))/a;
  296.     if (cp >= 1.0)    /* the /a can be slightly > 1 */
  297.         *p = 0.0;
  298.     else if (cp <= -1.0)
  299.         *p = PI;
  300.     else
  301.         *p = acos ((sy - (sinlastlat*sq))/a);
  302.     if (sx>0) *p = 2.0*PI - *p;
  303. #endif
  304. }
  305. xXx
  306. echo x altmenus.c
  307. cat > altmenus.c << 'xXx'
  308. /* printing routines for the three alternative bottom half menus,
  309.  * "menu1", "menu2" and "menu3".
  310.  */
  311.  
  312. #include <stdio.h>
  313. #include <math.h>
  314. #include "astro.h"
  315. #include "circum.h"
  316. #include "screen.h"
  317.  
  318. static int altmenu = F_MNU1;    /* which alternate menu is up; one of F_MNUi */
  319. static int alt2_stdhzn;    /* whether to use STDHZN (aot ADPHZN) horizon algthm  */
  320. static int alt3_geoc;    /* whether to use geocentric (aot topocentric) vantage*/
  321.  
  322. /* table of screen rows given a body #define from astro/h or screen.h */
  323. static short bodyrow[NOBJ] = {
  324.     R_MERCURY, R_VENUS, R_MARS, R_JUPITER, R_SATURN,
  325.     R_URANUS, R_NEPTUNE, R_PLUTO, R_SUN, R_MOON, R_OBJX
  326. };
  327. /* table of screen cols for third menu format, given body #define ... */
  328. static short bodycol[NOBJ] = {
  329.     C_MERCURY, C_VENUS, C_MARS, C_JUPITER, C_SATURN,
  330.     C_URANUS, C_NEPTUNE, C_PLUTO, C_SUN, C_MOON, C_OBJX
  331. };
  332.  
  333. /* let op decide which alternate menu should be up,
  334.  * including any menu-specific setup they might require.
  335.  * return 0 if things changed to require updating the alt menu; else -1.
  336.  */
  337. altmenu_setup()
  338. {
  339.     static char *flds[5] = {
  340.         "Data menu", "Rise/Set menu", "Separations menu"
  341.     };
  342.     int newmenu = altmenu, newhzn = alt2_stdhzn, newgeoc = alt3_geoc;
  343.     int new;
  344.     int fn = altmenu == F_MNU1 ? 0 : altmenu == F_MNU2 ? 1 : 2;
  345.  
  346.     ask:
  347.     flds[3]= newhzn ? "Standard hzn" : "Adaptive hzn";
  348.     flds[4]= newgeoc? "Geocentric" : "Topocentric";
  349.  
  350.     switch (popup (flds, fn, 5)) {
  351.     case 0: newmenu = F_MNU1; break;
  352.     case 1: newmenu = F_MNU2; break;
  353.     case 2: newmenu = F_MNU3; break;
  354.     case 3: newhzn ^= 1; fn = 3; goto ask;
  355.     case 4: newgeoc ^= 1; fn = 4; goto ask;
  356.     default: return (-1);
  357.     }
  358.  
  359.     new = 0;
  360.     if (newmenu != altmenu) {
  361.         altmenu = newmenu;
  362.         new++;
  363.     }
  364.     if (newhzn != alt2_stdhzn) {
  365.         alt2_stdhzn = newhzn;
  366.         if (newmenu == F_MNU2)
  367.         new++;
  368.     }
  369.     if (newgeoc != alt3_geoc) {
  370.         alt3_geoc = newgeoc;
  371.         if (newmenu == F_MNU3)
  372.         new++;
  373.     }
  374.     return (new ? 0 : -1);
  375. }
  376.  
  377. /* erase the info for the given planet */
  378. alt_nobody (p)
  379. int p;
  380. {
  381.     c_pos (bodyrow[p], C_RA);
  382.     c_eol();
  383. }
  384.  
  385. alt_body (b, force, np)
  386. int b;        /* which body, ala astro.h and screen.h defines */
  387. int force;    /* if !0 then draw for sure, else just if changed since last */
  388. Now *np;
  389. {
  390.     switch (altmenu) {
  391.     case F_MNU1: alt1_body (b, force, np); break;
  392.     case F_MNU2: alt2_body (b, force, np); break;
  393.     case F_MNU3: alt3_body (b, force, np); break;
  394.     }
  395. }
  396.  
  397. /* draw the labels for the current alternate menu format */
  398. alt_labels ()
  399. {
  400.     switch (altmenu) {
  401.     case F_MNU1: alt1_labels (); break;
  402.     case F_MNU2: alt2_labels (); break;
  403.     case F_MNU3: alt3_labels (); break;
  404.     }
  405. }
  406.  
  407. /* erase the labels for the current alternate menu format */
  408. alt_nolabels ()
  409. {
  410.     int i;
  411.  
  412.     f_string (R_ALTM, C_ALTMV, "             ");
  413.     for (i = R_PLANTAB; i < R_SUN; i++) {
  414.         c_pos (i, C_RA);
  415.         c_eol();
  416.     }
  417. }
  418.  
  419. alt_menumask()
  420. {
  421.     return (altmenu);
  422. }
  423.  
  424. /* handy function to return the next planet in the order in which they are
  425.  * displayed in the lower half of the screen.
  426.  * input is a given planet, return is the next planet.
  427.  * if input is not legal, then first planet is returned; when input is the
  428.  * last planet, then -1 is returned.
  429.  * typical usage is something like:
  430.  *   for (p = nxtbody(-1); p != -1; p = nxtbody(p))
  431.  */
  432. nxtbody(c)
  433. int c;
  434. {
  435.     static short nxtpl[NOBJ] = {
  436.         VENUS, MARS, JUPITER, SATURN, URANUS,
  437.         NEPTUNE, PLUTO, OBJX, MOON, MERCURY, -1
  438.     };
  439.  
  440.     if (c < MERCURY || c >= NOBJ)
  441.         return (SUN);
  442.     else
  443.         return (nxtpl[c]);
  444. }
  445.  
  446. static
  447. alt1_labels()
  448. {
  449.     f_string (R_ALTM, C_ALTMV, "  Planet Data");
  450.  
  451.     f_string (R_PLANTAB,    C_RA,    "R.A.");
  452.     f_string (R_PLANTAB+1,    C_RA,    "Hr:Mn.d");
  453.     f_string (R_PLANTAB,    C_DEC,    "Dec");
  454.     f_string (R_PLANTAB+1,    C_DEC,    "Deg:Mn");
  455.     f_string (R_PLANTAB,    C_AZ,    "Az");
  456.     f_string (R_PLANTAB+1,    C_AZ,    "Deg E");
  457.     f_string (R_PLANTAB,    C_ALT,    "Alt");
  458.     f_string (R_PLANTAB+1,    C_ALT,    "Deg Up");
  459.     f_string (R_PLANTAB,    C_HLONG,"Helio");
  460.     f_string (R_PLANTAB+1,    C_HLONG,"Long");
  461.     f_string (R_PLANTAB,    C_HLAT,    "Helio");
  462.     f_string (R_PLANTAB+1,    C_HLAT,    "Lat");
  463.     f_string (R_PLANTAB,    C_EDIST,"Ea Dst");
  464.     f_string (R_PLANTAB+1,    C_EDIST,"AU(mi)");
  465.     f_string (R_PLANTAB,    C_SDIST,"Sn Dst");
  466.     f_string (R_PLANTAB+1,    C_SDIST,"AU");
  467.     f_string (R_PLANTAB,    C_ELONG,"Elong");
  468.     f_string (R_PLANTAB+1,    C_ELONG,"Deg E");
  469.     f_string (R_PLANTAB,    C_SIZE,    "Size");
  470.     f_string (R_PLANTAB+1,    C_SIZE,    "ArcS");
  471.     f_string (R_PLANTAB,    C_MAG,    "VMag");
  472.     f_string (R_PLANTAB,    C_PHASE,"Phs");
  473.     f_char (R_PLANTAB+1,    C_PHASE,'%');
  474. }
  475.  
  476. static
  477. alt2_labels()
  478. {
  479.     f_string (R_ALTM, C_ALTMV, "Rise/Set Info");
  480.  
  481.     f_string (R_PLANTAB,    C_RISETM+5,    "Rise");
  482.     f_string (R_PLANTAB+1,    C_RISETM+1,    "Time");
  483.     f_string (R_PLANTAB+1,    C_RISEAZ+2,    "Az");
  484.     f_string (R_PLANTAB,    C_TRANSTM+3,    "Transit");
  485.     f_string (R_PLANTAB+1,    C_TRANSTM+1,    "Time");
  486.     f_string (R_PLANTAB+1,    C_TRANSALT+2,    "Alt");
  487.     f_string (R_PLANTAB,    C_SETTM+5,    "Set");
  488.     f_string (R_PLANTAB+1,    C_SETTM+1,    "Time");
  489.     f_string (R_PLANTAB+1,    C_SETAZ+2,    "Az");
  490.     f_string (R_PLANTAB,    C_TUP,        "Hrs Up");
  491. }
  492.  
  493. static
  494. alt3_labels()
  495. {
  496.     f_string (R_ALTM, C_ALTMV, "  Separations");
  497.  
  498.     f_string (R_PLANTAB,    C_SUN+2,    "Sun");
  499.     f_string (R_PLANTAB,    C_MOON+1,    "Moon");
  500.     f_string (R_PLANTAB,    C_MERCURY+1,    "Merc");
  501.     f_string (R_PLANTAB,    C_VENUS+1,    "Venus");
  502.     f_string (R_PLANTAB,    C_MARS+1,    "Mars");
  503.     f_string (R_PLANTAB,    C_JUPITER+2,    "Jup");
  504.     f_string (R_PLANTAB,    C_SATURN,    "Saturn");
  505.     f_string (R_PLANTAB,    C_URANUS,    "Uranus");
  506.     f_string (R_PLANTAB,    C_NEPTUNE+2,    "Nep");
  507.     f_string (R_PLANTAB,    C_PLUTO+1,    "Pluto");
  508.  
  509.     if (objx_ison()) {
  510.         char xnm[MAXOBJXNM+1];
  511.         char buf[10];
  512.         objx_get ((double *)0, (double *)0, (double *)0, xnm);
  513.         sprintf (buf, "%-2.2s", xnm); /* set all spaces */
  514.         f_string (R_PLANTAB, C_OBJX, buf);
  515.     }
  516. }
  517.  
  518. /* print body info in first menu format */
  519. static
  520. alt1_body (p, force, np)
  521. int p;        /* which body, as in astro.h/screen.h defines */
  522. int force;    /* whether to print for sure or only if things have changed */
  523. Now *np;
  524. {
  525.     Sky sky;
  526.     double as = plot_ison() || srch_ison() ? 0.0 : 60.0;
  527.     int row = bodyrow[p];
  528.  
  529.     if (body_cir (p, as, np, &sky) || force) {
  530.         if (p == OBJX)
  531.         pr_objxname();
  532.         f_ra (row, C_RA, sky.s_ra);
  533.         f_angle (row, C_DEC, sky.s_dec);
  534.         if (p != MOON && p != OBJX) {
  535.         f_angle (row, C_HLONG, sky.s_hlong);
  536.         f_angle (row, C_HLAT, sky.s_hlat);
  537.         }
  538.  
  539.         if (p == MOON) {
  540.         /* distance is on km, show in miles */
  541.         f_double (R_MOON, C_EDIST, "%6.0f", sky.s_edist/1.609344);
  542.         } else if (p != OBJX) {
  543.         /* show distance in au */
  544.         f_double (row, C_EDIST,(sky.s_edist>=10.0)?"%6.3f":"%6.4f",
  545.                                 sky.s_edist);
  546.         }
  547.         if (p != OBJX)
  548.         f_double (row, C_SDIST, (sky.s_sdist>=10.0)?"%6.3f":"%6.4f",
  549.                                 sky.s_sdist);
  550.         f_double (row, C_ELONG, "%6.1f", sky.s_elong);
  551.         if (p != OBJX) {
  552.         f_double (row, C_SIZE, (p==MOON||p==SUN)?"%4.0f":"%4.1f",
  553.                                 sky.s_size);
  554.         f_double (row, C_MAG, (p==MOON||p==SUN)?"%4.0f":"%4.1f",
  555.                                 sky.s_mag);
  556.         if (p != SUN) {
  557.             char buf[10];
  558.             /* would just do this if Turbo-C 2.0 "%?.0f" worked:
  559.              * f_double (row, C_PHASE, "%3.0f", sky.s_phase);
  560.              */
  561.             sprintf (buf, "%3d", (int)(sky.s_phase+0.5));
  562.             f_string (row, C_PHASE, buf);
  563.         }
  564.         }
  565.     }
  566.  
  567.     f_angle (row, C_AZ, sky.s_az);
  568.     f_angle (row, C_ALT, sky.s_alt);
  569. }
  570.  
  571. /* print body info in the second menu format */
  572. static
  573. alt2_body (p, force, np)
  574. int p;        /* which body, as in astro.h/screen.h defines */
  575. int force;    /* whether to print for sure or only if things have changed */
  576. Now *np;
  577. {
  578.     double ltr, lts, ltt, azr, azs, altt;
  579.     int row = bodyrow[p];
  580.     int status;
  581.     double tmp;
  582.     int today_tup = 0;
  583.  
  584.     if (!riset_cir (p, np, alt2_stdhzn?STDHZN:ADPHZN, <r, <s, <t,
  585.                     &azr, &azs, &altt, &status) && !force)
  586.         return;
  587.  
  588.     alt_nobody (p);
  589.     if (p == OBJX)
  590.         pr_objxname();
  591.  
  592.     if (status & RS_ERROR) {
  593.         /* can not find where body is! */
  594.         f_string (row, C_RISETM, "?Error?");
  595.         return;
  596.     }
  597.     if (status & RS_CIRCUMPOLAR) {
  598.         /* body is up all day */
  599.         f_string (row, C_RISETM, "Circumpolar");
  600.         f_mtime (row, C_TRANSTM, ltt);
  601.         if (status & RS_2TRANS)
  602.         f_char (row, C_TRANSTM+5, '+');
  603.         f_angle (row, C_TRANSALT, altt);
  604.         f_string (row, C_TUP, "24:00"); /*f_mtime() changes to 0:00 */
  605.         return;
  606.     }
  607.     if (status & RS_NEVERUP) {
  608.         /* body never up at all today */
  609.         f_string (row, C_RISETM, "Never up");
  610.         f_mtime (row, C_TUP, 0.0);
  611.         return;
  612.     }
  613.  
  614.     if (status & RS_NORISE) {
  615.         /* object does not rise as such today */
  616.         f_string (row, C_RISETM, "Never rises");
  617.         ltr = 0.0; /* for TUP */
  618.         today_tup = 1;
  619.     } else {
  620.         f_mtime (row, C_RISETM, ltr);
  621.         if (status & RS_2RISES) {
  622.         /* object rises more than once today */
  623.         f_char (row, C_RISETM+5, '+');
  624.         }
  625.         f_angle (row, C_RISEAZ, azr);
  626.     }
  627.  
  628.     if (status & RS_NOTRANS)
  629.         f_string (row, C_TRANSTM, "No transit");
  630.     else {
  631.         f_mtime (row, C_TRANSTM, ltt);
  632.         if (status & RS_2TRANS)
  633.         f_char (row, C_TRANSTM+5, '+');
  634.         f_angle (row, C_TRANSALT, altt);
  635.     }
  636.  
  637.     if (status & RS_NOSET) {
  638.         /* object does not set as such today */
  639.         f_string (row, C_SETTM, "Never sets");
  640.         lts = 24.0;    /* for TUP */
  641.         today_tup = 1;
  642.     } else {
  643.         f_mtime (row, C_SETTM, lts);
  644.         if (status & RS_2SETS)
  645.         f_char (row, C_SETTM+5, '+');
  646.         f_angle (row, C_SETAZ, azs);
  647.     }
  648.  
  649.     tmp = lts - ltr;
  650.     if (tmp < 0)
  651.         tmp = 24.0 + tmp;
  652.     f_mtime (row, C_TUP, tmp);
  653.     if (today_tup)
  654.         f_char (row, C_TUP+5, '+');
  655. }
  656.  
  657. /* print body info in third menu format. this may be either the geocentric
  658.  *   or topocentric angular separation between object p and each of the others.
  659.  *   the latter, of course, includes effects of refraction and so can change
  660.  *   quite rapidly near the time of each planets rise or set.
  661.  * for now, we don't save old values so we always redo everything and ignore
  662.  *  the "force" argument. this isn't that bad since body_cir() has memory and
  663.  *   will avoid most computations as we hit them again in the lower triangle.
  664.  */
  665. /*ARGSUSED*/
  666. static
  667. alt3_body (p, force, np)
  668. int p;        /* which body, as in astro.h/screen.h defines */
  669. int force;    /* whether to print for sure or only if things have changed */
  670. Now *np;
  671. {
  672.     int row = bodyrow[p];
  673.     Sky skyp, skyq;
  674.     double spy, cpy, px, *qx, *qy;
  675.     int wantx = objx_ison();
  676.     double as = plot_ison() || srch_ison() ? 0.0 : 60.0;
  677.     int q;
  678.  
  679.     (void) body_cir (p, as, np, &skyp);
  680.     if (alt3_geoc) {
  681.         /* use ra for "x", dec for "y". */
  682.         spy = sin (skyp.s_dec);
  683.         cpy = cos (skyp.s_dec);
  684.         px = skyp.s_ra;
  685.         qx = &skyq.s_ra;
  686.         qy = &skyq.s_dec;
  687.     } else {
  688.         /* use azimuth for "x", altitude for "y". */
  689.         spy = sin (skyp.s_alt);
  690.         cpy = cos (skyp.s_alt);
  691.         px = skyp.s_az;
  692.         qx = &skyq.s_az;
  693.         qy = &skyq.s_alt;
  694.     }
  695.     if (p == OBJX)
  696.         pr_objxname();
  697.     for (q = nxtbody(-1); q != -1; q = nxtbody(q))
  698.         if (q != p && (q != OBJX || wantx)) {
  699.         double csep;
  700.         (void) body_cir (q, as, np, &skyq);
  701.         csep = spy*sin(*qy) + cpy*cos(*qy)*cos(px-*qx);
  702.         f_angle (row, bodycol[q], acos(csep));
  703.         }
  704. }
  705.  
  706. static
  707. pr_objxname()
  708. {
  709.     char n[MAXOBJXNM+1];
  710.     char buf[10];
  711.     objx_get ((double *)0, (double *)0, (double *)0, n);
  712.     sprintf (buf, "%-2.2s", n); /* set all spaces */
  713.     f_string (R_OBJX, C_OBJ, buf);
  714. }
  715. xXx
  716. echo x anomaly.c
  717. cat > anomaly.c << 'xXx'
  718. #include <stdio.h>
  719. #include <math.h>
  720. #include "astro.h"
  721.  
  722. #define    TWOPI    (2*PI)
  723.  
  724. /* given the mean anomaly, ma, and the eccentricity, s, of elliptical motion,
  725.  * find the true anomaly, *nu, and the eccentric anomaly, *ea.
  726.  * all angles in radians.
  727.  */
  728. anomaly (ma, s, nu, ea)
  729. double ma, s;
  730. double *nu, *ea;
  731. {
  732.     double m, dla, fea;
  733.  
  734.     m = ma-TWOPI*(int)(ma/TWOPI);
  735.     fea = m;
  736.     while (1) {
  737.         dla = fea-(s*sin(fea))-m;
  738.         if (fabs(dla)<1e-6)
  739.         break;
  740.         dla /= 1-(s*cos(fea));
  741.         fea -= dla;
  742.     }
  743.  
  744.     *nu = 2*atan(sqrt((1+s)/(1-s))*tan(fea/2));
  745.     *ea = fea;
  746. }
  747. xXx
  748. echo x cal_mjd.c
  749. cat > cal_mjd.c << 'xXx'
  750. #include <stdio.h>
  751. #include <math.h>
  752. #include "astro.h"
  753.  
  754. /* given a date in months, mn, days, dy, years, yr,
  755.  * return the modified Julian date (number of days elapsed since 1900 jan 0.5),
  756.  * *mjd.
  757.  */
  758. cal_mjd (mn, dy, yr, mjd)
  759. int mn, yr;
  760. double dy;
  761. double *mjd;
  762. {
  763.     int b, d, m, y;
  764.     long c;
  765.  
  766.     m = mn;
  767.     y = (yr < 0) ? yr + 1 : yr;
  768.     if (mn < 3) {
  769.         m += 12;
  770.         y -= 1;
  771.     }
  772.  
  773.     if (yr < 1582 || yr == 1582 && (mn < 10 || mn == 10 && dy < 15)) 
  774.         b = 0;
  775.     else {
  776.         int a;
  777.         a = y/100;
  778.         b = 2 - a + a/4;
  779.     }
  780.  
  781.     if (y < 0)
  782.         c = (long)((365.25*y) - 0.75) - 694025L;
  783.     else
  784.         c = (long)(365.25*y) - 694025L;
  785.  
  786.     d = 30.6001*(m+1);
  787.  
  788.     *mjd = b + c + d + dy - 0.5;
  789. }
  790.  
  791. /* given the modified Julian date (number of days elapsed since 1900 jan 0.5,),
  792.  * mjd, return the calendar date in months, *mn, days, *dy, and years, *yr.
  793.  */
  794. mjd_cal (mjd, mn, dy, yr)
  795. double mjd;
  796. int *mn, *yr;
  797. double *dy;
  798. {
  799.     double d, f;
  800.     double i, a, b, ce, g;
  801.  
  802.     d = mjd + 0.5;
  803.     i = floor(d);
  804.     f = d-i;
  805.     if (f == 1) {
  806.         f = 0;
  807.         i += 1;
  808.     }
  809.  
  810.     if (i > -115860.0) {
  811.         a = floor((i/36524.25)+.9983573)+14;
  812.         i += 1 + a - floor(a/4.0);
  813.     }
  814.  
  815.     b = floor((i/365.25)+.802601);
  816.     ce = i - floor((365.25*b)+.750001)+416;
  817.     g = floor(ce/30.6001);
  818.     *mn = g - 1;
  819.     *dy = ce - floor(30.6001*g)+f;
  820.     *yr = b + 1899;
  821.  
  822.     if (g > 13.5)
  823.         *mn = g - 13;
  824.     if (*mn < 2.5)
  825.         *yr = b + 1900;
  826.     if (*yr < 1)
  827.         *yr -= 1;
  828. }
  829.  
  830. /* given an mjd, set *dow to 0..6 according to which dayof the week it falls
  831.  * on (0=sunday) or set it to -1 if can't figure it out.
  832.  */
  833. mjd_dow (mjd, dow)
  834. double mjd;
  835. int *dow;
  836. {
  837.     /* cal_mjd() uses Gregorian dates on or after Oct 15, 1582.
  838.      * (Pope Gregory XIII dropped 10 days, Oct 5..14, and improved the leap-
  839.      * year algorithm). however, Great Britian and the colonies did not
  840.      * adopt it until Sept 14, 1752 (they dropped 11 days, Sept 3-13,
  841.      * due to additional accumulated error). leap years before 1752 thus
  842.      * can not easily be accounted for from the cal_mjd() number...
  843.      */
  844.     if (mjd < -53798.5) {
  845.         /* pre sept 14, 1752 too hard to correct */
  846.         *dow = -1;
  847.         return;
  848.     }
  849.     *dow = ((int)floor(mjd-.5) + 1) % 7; /* 1/1/1900 (mjd 0.5) is a Monday*/
  850.     if (*dow < 0)
  851.         *dow += 7;
  852. }
  853.  
  854. /* given a mjd, return the the number of days in the month.  */
  855. mjd_dpm (mjd, ndays)
  856. double mjd;
  857. int *ndays;
  858. {
  859.     static short dpm[] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
  860.     int m, y;
  861.     double d;
  862.  
  863.     mjd_cal (mjd, &m, &d, &y);
  864.     *ndays = (m==2 && ((y%4==0 && y%100!=0)||y%400==0)) ? 29 : dpm[m-1];
  865. }
  866.  
  867.  
  868. /* given a mjd, return the year as a double. */
  869. mjd_year (mjd, yr)
  870. double mjd;
  871. double *yr;
  872. {
  873.     int m, y;
  874.     double d;
  875.     double e0, e1;    /* mjd of start of this year, start of next year */
  876.  
  877.     mjd_cal (mjd, &m, &d, &y);
  878.     cal_mjd (1, 1.0, y, &e0);
  879.     cal_mjd (12, 31.0, y, &e1); e1 += 1.0;
  880.     *yr = y + (mjd - e0)/(e1 - e0);
  881. }
  882. xXx
  883. echo x circum.c
  884. cat > circum.c << 'xXx'
  885. /* fill in a Sky struct with all we know about each object.
  886.  *(object-x is in objx.c)
  887.  */
  888.  
  889. #include <stdio.h>
  890. #include <math.h>
  891. #include "astro.h"
  892. #include "circum.h"
  893. #include "screen.h"    /* just for SUN and MOON */
  894.  
  895. /* find body p's circumstances now.
  896.  * to save some time the caller may specify a desired accuracy, in arc seconds.
  897.  * if, based on its mean motion, it would not have moved this much since the
  898.  * last time we were called we only recompute altitude and azimuth and avoid
  899.  * recomputing the planet's heliocentric position. use 0.0 for best possible.
  900.  * return 0 if only alt/az changes, else 1 if all other stuff updated too.
  901.  * TODO: beware of fast retrograde motions.
  902.  */
  903. body_cir (p, as, np, sp)
  904. int p;
  905. double as;
  906. Now *np;
  907. Sky *sp;
  908. {
  909.     typedef struct {
  910.         double l_dpas;    /* mean days per arc second */
  911.         Now l_now;        /* when l_sky was found */
  912.         Sky l_sky;
  913.     } Last;
  914.     /* must be in same order as the astro.h #define's */
  915.     static Last last[8] =
  916.         {{.000068},{.00017},{.00053},{.0034},{.0092},{.027},{.046},{.069}};
  917.     double lst, alt, az;
  918.     Last *lp;
  919.     int new;
  920.  
  921.     switch (p) {
  922.     case SUN:
  923.         return (sun_cir (as, np, sp));
  924.     case MOON:
  925.         return (moon_cir (as, np, sp));
  926.     case OBJX:
  927.         return (objx_cir (as, np, sp));
  928.     }
  929.  
  930.     /* if less than l_every days from last time for this planet
  931.      * just redo alt/az.
  932.      */
  933.     lp = last + p;
  934.     if (same_cir(np,&lp->l_now) && about_now(np,&lp->l_now,as*lp->l_dpas)) {
  935.         *sp = lp->l_sky;
  936.         new = 0;
  937.     } else {
  938.         double lpd0, psi0; /* heliocentric ecliptic longitude and latitude */
  939.         double rp0;        /* dist from sun */
  940.         double rho0;    /* dist from earth */
  941.         double lam, bet;    /* geocentric ecliptic long and lat */
  942.         double dia, mag;    /* angular diameter at 1 AU and magnitude */
  943.         double lsn, rsn;    /* true geoc lng of sun, dist from sn to earth*/
  944.         double deps, dpsi;
  945.         double a, ca, sa;
  946.         double el;    /* elongation */
  947.         double f;    /* phase from earth */
  948.  
  949.         lp->l_now = *np;
  950.         plans (mjd, p, &lpd0, &psi0, &rp0, &rho0, &lam, &bet, &dia, &mag);
  951.         nutation (mjd, &deps, &dpsi); /* correct for nutation */
  952.         lam += dpsi;
  953.         sunpos (mjd, &lsn, &rsn);
  954.         /* correct for 20.4" aberration */
  955.         a = lsn-lam;
  956.         ca = cos(a);
  957.         sa = sin(a);
  958.         lam -= degrad(20.4/3600)*ca/cos(bet);
  959.         bet -= degrad(20.4/3600)*sa*sin(bet);
  960.  
  961.         ecl_eq (mjd, bet, lam, &sp->s_ra, &sp->s_dec);
  962.         if (epoch != EOD)
  963.         precess (mjd, epoch, &sp->s_ra, &sp->s_dec);
  964.         sp->s_edist = rho0;
  965.         sp->s_sdist = rp0;
  966.         elongation (lam, bet, lsn, &el);
  967.         el = raddeg(el);
  968.         sp->s_elong = el;
  969.         sp->s_size = dia/rho0;
  970.         f = 0.5*(1+cos(lam-lpd0));
  971.         sp->s_phase = f*100.0; /* percent */
  972.         sp->s_mag = 5.0*log(rp0*rho0/sqrt(f))/log(10.0) + mag;
  973.         sp->s_hlong = lpd0;
  974.         sp->s_hlat = psi0;
  975.         new = 1;
  976.     }
  977.  
  978.     /* alt, az; correct for refraction, in place */
  979.     now_lst (np, &lst);
  980.     hadec_aa (lat, hrrad(lst) - sp->s_ra, sp->s_dec, &alt, &az);
  981.     refract (pressure, temp, alt, &alt);
  982.     sp->s_alt = alt;
  983.     sp->s_az = az;
  984.     lp->l_sky = *sp;
  985.     return (new);
  986. }
  987.  
  988. /* find local times when sun is 18 degrees below horizon.
  989.  * return 0 if just returned same stuff as previous call, else 1 if new.
  990.  */
  991. twilight_cir (np, dawn, dusk, status)
  992. Now *np;
  993. double *dawn, *dusk;
  994. int *status;
  995. {
  996.     static Now last_now;
  997.     static double last_dawn, last_dusk;
  998.     static int last_status;
  999.     int new;
  1000.  
  1001.     if (same_cir (np, &last_now) && same_lday (np, &last_now)) {
  1002.         *dawn = last_dawn;
  1003.         *dusk = last_dusk;
  1004.         *status = last_status;
  1005.         new = 0;
  1006.     } else {
  1007.         double x;
  1008.         (void) riset_cir (SUN,np,TWILIGHT,dawn,dusk,&x,&x,&x,&x,status);
  1009.         last_dawn = *dawn;
  1010.         last_dusk = *dusk;
  1011.         last_status = *status;
  1012.         last_now = *np;
  1013.         new = 1;
  1014.     }
  1015.     return (new);
  1016. }
  1017.  
  1018. /* find sun's circumstances now.
  1019.  * as is the desired accuracy, in arc seconds; use 0.0 for best possible.
  1020.  * return 0 if only alt/az changes, else 1 if all other stuff updated too.
  1021.  */
  1022. sun_cir (as, np, sp)
  1023. double as;
  1024. Now *np;
  1025. Sky *sp;
  1026. {
  1027.     static Sky last_sky;
  1028.     static Now last_now;
  1029.     double lst, alt, az;
  1030.     int new;
  1031.  
  1032.     if (same_cir (np, &last_now) && about_now (np, &last_now, as*.00028)) {
  1033.         *sp = last_sky;
  1034.         new = 0;
  1035.     } else {
  1036.         double lsn, rsn;
  1037.         double deps, dpsi;
  1038.  
  1039.         last_now = *np;
  1040.         sunpos (mjd, &lsn, &rsn);        /* sun's true ecliptic long
  1041.                          * and dist
  1042.                          */
  1043.         nutation (mjd, &deps, &dpsi);    /* correct for nutation */
  1044.         lsn += dpsi-degrad(20.4/3600);    /* and 20.4" aberration */
  1045.  
  1046.         sp->s_edist = rsn;
  1047.         sp->s_sdist = 0.0;
  1048.         sp->s_elong = 0.0;
  1049.         sp->s_size = raddeg(4.65242e-3/rsn)*3600*2;
  1050.         sp->s_mag = -26.8;
  1051.         sp->s_hlong = lsn-PI;    /* geo- to helio- centric */
  1052.         range (&sp->s_hlong, 2*PI);
  1053.         sp->s_hlat = 0.0;
  1054.  
  1055.         ecl_eq (mjd, 0.0, lsn, &sp->s_ra, &sp->s_dec);
  1056.         if (epoch != EOD)
  1057.         precess (mjd, epoch, &sp->s_ra, &sp->s_dec);
  1058.         new = 1;
  1059.     }
  1060.  
  1061.     now_lst (np, &lst);
  1062.     hadec_aa (lat, hrrad(lst) - sp->s_ra, sp->s_dec, &alt, &az);
  1063.     refract (pressure, temp, alt, &alt);
  1064.     sp->s_alt = alt;
  1065.     sp->s_az = az;
  1066.     last_sky = *sp;
  1067.     return (new);
  1068. }
  1069.  
  1070. /* find moon's circumstances now.
  1071.  * as is the desired accuracy, in arc seconds; use 0.0 for best possible.
  1072.  * return 0 if only alt/az changes, else 1 if all other stuff updated too.
  1073.  */
  1074. moon_cir (as, np, sp)
  1075. double as;
  1076. Now *np;
  1077. Sky *sp;
  1078. {
  1079.     static Sky last_sky;
  1080.     static Now last_now;
  1081.     static double ehp;
  1082.     double lst, alt, az;
  1083.     double ha, dec;
  1084.     int new;
  1085.  
  1086.     if (same_cir (np, &last_now) && about_now (np, &last_now, as*.000021)) {
  1087.         *sp = last_sky;
  1088.         new = 0;
  1089.     } else {
  1090.         double lam, bet;
  1091.         double deps, dpsi;
  1092.         double lsn, rsn;    /* sun long in rads, earth-sun dist in au */
  1093.         double edistau;    /* earth-moon dist, in au */
  1094.         double el;        /* elongation, rads east */
  1095.  
  1096.         last_now = *np;
  1097.         moon (mjd, &lam, &bet, &ehp);    /* moon's true ecliptic loc */
  1098.         nutation (mjd, &deps, &dpsi);    /* correct for nutation */
  1099.         lam += dpsi;
  1100.         range (&lam, 2*PI);
  1101.  
  1102.         sp->s_edist = 6378.14/sin(ehp);    /* earth-moon dist, want km */
  1103.         sp->s_size = 3600*31.22512*sin(ehp);/* moon angular dia, seconds */
  1104.  
  1105.         ecl_eq (mjd, bet, lam, &sp->s_ra, &sp->s_dec);
  1106.         if (epoch != EOD)
  1107.         precess (mjd, epoch, &sp->s_ra, &sp->s_dec);
  1108.  
  1109.         sunpos (mjd, &lsn, &rsn);
  1110.         range (&lsn, 2*PI);
  1111.         elongation (lam, bet, lsn, &el);
  1112.  
  1113.         /* solve triangle of earth, sun, and elongation for moon-sun dist */
  1114.         edistau = sp->s_edist/1.495979e8; /* km -> au */
  1115.         sp->s_sdist =
  1116.         sqrt (edistau*edistau + rsn*rsn - 2.0*edistau*rsn*cos(el));
  1117.  
  1118.         /* TODO: improve mag; this is based on a flat moon model. */
  1119.         sp->s_mag = -12 + 2.5*(log10(PI) - log10(PI/2*(1+1.e-6-cos(el))));
  1120.  
  1121.         sp->s_elong = raddeg(el);    /* want degrees */
  1122.         sp->s_phase = fabs(el)/PI*100.0;    /* want non-negative % */
  1123.         sp->s_hlong = sp->s_hlat = 0.0;
  1124.         new = 1;
  1125.     }
  1126.  
  1127.     /* show topocentric alt/az by correcting ra/dec for parallax 
  1128.      * as well as refraction.
  1129.      */
  1130.     now_lst (np, &lst);
  1131.     ha = hrrad(lst) - sp->s_ra;
  1132.     ta_par (ha, sp->s_dec, lat, height, ehp, &ha, &dec);
  1133.     hadec_aa (lat, ha, dec, &alt, &az);
  1134.     refract (pressure, temp, alt, &alt);
  1135.     sp->s_alt = alt;
  1136.     sp->s_az = az;
  1137.     last_sky = *sp;
  1138.     return (new);
  1139. }
  1140.  
  1141. /* given geocentric ecliptic longitude and latitude, lam and bet, of some object
  1142.  * and the longitude of the sun, lsn, find the elongation, el. this is the
  1143.  * actual angular separation of the object from the sun, not just the difference
  1144.  * in the longitude. the sign, however, IS set simply as a test on longitude
  1145.  * such that el will be >0 for an evening object <0 for a morning object.
  1146.  * to understand the test for el sign, draw a graph with lam going from 0-2*PI
  1147.  *   down the vertical axis, lsn going from 0-2*PI across the hor axis. then
  1148.  *   define the diagonal regions bounded by the lines lam=lsn+PI, lam=lsn and
  1149.  *   lam=lsn-PI. the "morning" regions are any values to the lower left of the
  1150.  *   first line and bounded within the second pair of lines.
  1151.  * all angles in radians.
  1152.  */
  1153. elongation (lam, bet, lsn, el)
  1154. double lam, bet, lsn;
  1155. double *el;
  1156. {
  1157.     *el = acos(cos(bet)*cos(lam-lsn));
  1158.     if (lam>lsn+PI || lam>lsn-PI && lam<lsn) *el = - *el;
  1159. }
  1160.  
  1161. /* return whether the two Nows are for the same observing circumstances. */
  1162. same_cir (n1, n2)
  1163. register Now *n1, *n2;
  1164. {
  1165.     return (n1->n_lat == n2->n_lat
  1166.         && n1->n_lng == n2->n_lng
  1167.         && n1->n_temp == n2->n_temp
  1168.         && n1->n_pressure == n2->n_pressure
  1169.         && n1->n_height == n2->n_height
  1170.         && n1->n_tz == n2->n_tz
  1171.         && n1->n_epoch == n2->n_epoch);
  1172. }
  1173.  
  1174. /* return whether the two Nows are for the same LOCAL day */
  1175. same_lday (n1, n2)
  1176. Now *n1, *n2;
  1177. {
  1178.     return (mjd_day(n1->n_mjd - n1->n_tz/24.0) ==
  1179.         mjd_day(n2->n_mjd - n2->n_tz/24.0)); 
  1180. }
  1181.  
  1182. /* return whether the mjd of the two Nows are within dt */
  1183. static
  1184. about_now (n1, n2, dt)
  1185. Now *n1, *n2;
  1186. double dt;
  1187. {
  1188.     return (fabs (n1->n_mjd - n2->n_mjd) <= dt/2.0);
  1189. }
  1190.  
  1191. now_lst (np, lst)
  1192. Now *np;
  1193. double *lst;
  1194. {
  1195.     utc_gst (mjd_day(mjd), mjd_hr(mjd), lst);
  1196.     *lst += radhr(lng);
  1197.     range (lst, 24.0);
  1198. }
  1199.  
  1200. /* round a time in days, *t, to the nearest second, IN PLACE. */
  1201. rnd_second (t)
  1202. double *t;
  1203. {
  1204.     *t = floor(*t*SPD+0.5)/SPD;
  1205. }
  1206.     
  1207. double mjd_day(jd)
  1208. double jd;
  1209. {
  1210.     return (floor(jd-0.5)+0.5);
  1211. }
  1212.  
  1213. double mjd_hr(jd)
  1214. double jd;
  1215. {
  1216.     return ((jd-mjd_day(jd))*24.0);
  1217. }
  1218. xXx
  1219. echo x compiler.c
  1220. cat > compiler.c << 'xXx'
  1221. /* module to compile and execute a c-style arithmetic expression.
  1222.  * public entry points are compile_expr() and execute_expr().
  1223.  *
  1224.  * one reason this is so nice and tight is that all opcodes are the same size
  1225.  * (an int) and the tokens the parser returns are directly usable as opcodes,
  1226.  * for the most part. constants and variables are compiled as an opcode
  1227.  * with an offset into the auxiliary opcode tape, opx.
  1228.  */
  1229.  
  1230. #include <math.h>
  1231. #include "screen.h"
  1232.  
  1233. /* parser tokens and opcodes, as necessary */
  1234. #define    HALT    0    /* good value for HALT since program is inited to 0 */
  1235. /* binary operators (precedences in table, below) */
  1236. #define    ADD    1
  1237. #define    SUB    2
  1238. #define    MULT    3
  1239. #define    DIV    4
  1240. #define    AND    5
  1241. #define    OR    6
  1242. #define    GT    7
  1243. #define    GE    8
  1244. #define    EQ    9
  1245. #define    NE    10
  1246. #define    LT    11
  1247. #define    LE    12
  1248. /* unary op, precedence in NEG_PREC #define, below */
  1249. #define    NEG    13
  1250. /* symantically operands, ie, constants, variables and all functions */
  1251. #define    CONST    14    
  1252. #define    VAR    15
  1253. #define    ABS    16    /* add functions if desired just like this is done */
  1254. /* purely tokens - never get compiled as such */
  1255. #define    LPAREN    255
  1256. #define    RPAREN    254
  1257. #define    ERR    (-1)
  1258.  
  1259. /* precedence of each of the binary operators.
  1260.  * in case of a tie, compiler associates left-to-right.
  1261.  * N.B. each entry's index must correspond to its #define!
  1262.  */
  1263. static int precedence[] = {0,5,5,6,6,2,1,4,4,3,3,4,4};
  1264. #define    NEG_PREC    7    /* negation is highest */
  1265.  
  1266. /* execute-time operand stack */
  1267. #define    MAX_STACK    16
  1268. static double stack[MAX_STACK], *sp;
  1269.  
  1270. /* space for compiled opcodes - the "program".
  1271.  * opcodes go in lower 8 bits.
  1272.  * when an opcode has an operand (as CONST and VAR) it is really in opx[] and
  1273.  *   the index is in the remaining upper bits.
  1274.  */
  1275. #define    MAX_PROG 32
  1276. static int program[MAX_PROG], *pc;
  1277. #define    OP_SHIFT    8
  1278. #define    OP_MASK        0xff
  1279.  
  1280. /* auxiliary operand info.
  1281.  * the operands (all but lower 8 bits) of CONST and VAR are really indeces
  1282.  * into this array. thus, no point in making this any longer than you have
  1283.  * bits more than 8 in your machine's int to index into it, ie, make
  1284.  *    MAX_OPX <= 1 << ((sizeof(int)-1)*8)
  1285.  * also, the fld's must refer to ones being flog'd, so not point in more
  1286.  * of these then that might be used for plotting and srching combined.
  1287.  */
  1288. #define    MAX_OPX    16
  1289. typedef union {
  1290.     double opu_f;        /* value when opcode is CONST */
  1291.     int opu_fld;        /* rcfpack() of field when opcode is VAR */
  1292. } OpX;
  1293. static OpX opx[MAX_OPX];
  1294. static int opxidx;
  1295.  
  1296. /* these are global just for easy/rapid access */
  1297. static int parens_nest;    /* to check that parens end up nested */
  1298. static char *err_msg;    /* caller provides storage; we point at it with this */
  1299. static char *cexpr, *lcexpr; /* pointers that move along caller's expression */
  1300. static int good_prog;    /* != 0 when program appears to be good */
  1301.  
  1302. /* compile the given c-style expression.
  1303.  * return 0 and set good_prog if ok,
  1304.  * else return -1 and a reason message in errbuf.
  1305.  */
  1306. compile_expr (ex, errbuf)
  1307. char *ex;
  1308. char *errbuf;
  1309. {
  1310.     int instr;
  1311.  
  1312.     /* init the globals.
  1313.      * also delete any flogs used in the previous program.
  1314.      */
  1315.     cexpr = ex;
  1316.     err_msg = errbuf;
  1317.     pc = program;
  1318.     opxidx = 0;
  1319.     parens_nest = 0;
  1320.     do {
  1321.         instr = *pc++;
  1322.         if ((instr & OP_MASK) == VAR)
  1323.         flog_delete (opx[instr >> OP_SHIFT].opu_fld);
  1324.     } while (instr != HALT);
  1325.  
  1326.     pc = program;
  1327.     if (compile(0) == ERR) {
  1328.         sprintf (err_msg + strlen(err_msg), " at \"%.10s\"", lcexpr);
  1329.         good_prog = 0;
  1330.         return (-1);
  1331.     }
  1332.     *pc++ = HALT;
  1333.     good_prog = 1;
  1334.     return (0);
  1335. }
  1336.  
  1337. /* execute the expression previously compiled with compile_expr().
  1338.  * return 0 with *vp set to the answer if ok, else return -1 with a reason
  1339.  * why not message in errbuf.
  1340.  */
  1341. execute_expr (vp, errbuf)
  1342. double *vp;
  1343. char *errbuf;
  1344. {
  1345.     int s;
  1346.  
  1347.     err_msg = errbuf;
  1348.     sp = stack + MAX_STACK;    /* grows towards lower addresses */
  1349.     pc = program;
  1350.     s = execute(vp);
  1351.     if (s < 0)
  1352.         good_prog = 0;
  1353.     return (s);
  1354. }
  1355.  
  1356. /* this is a way for the outside world to ask whether there is currently a
  1357.  * reasonable program compiled and able to execute.
  1358.  */
  1359. prog_isgood()
  1360. {
  1361.     return (good_prog);
  1362. }
  1363.  
  1364. /* get and return the opcode corresponding to the next token.
  1365.  * leave with lcexpr pointing at the new token, cexpr just after it.
  1366.  * also watch for mismatches parens and proper operator/operand alternation.
  1367.  */
  1368. static
  1369. next_token ()
  1370. {
  1371.     static char toomt[] = "More than %d terms";
  1372.     static char badop[] = "Illegal operator";
  1373.     int tok = ERR;    /* just something illegal */
  1374.     char c;
  1375.  
  1376.     while ((c = *cexpr) == ' ')
  1377.         cexpr++;
  1378.     lcexpr = cexpr++;
  1379.  
  1380.     /* mainly check for a binary operator */
  1381.     switch (c) {
  1382.     case '\0': --cexpr; tok = HALT; break; /* keep returning HALT */
  1383.     case '+': tok = ADD; break; /* compiler knows when it's really unary */
  1384.     case '-': tok = SUB; break; /* compiler knows when it's really negate */
  1385.     case '*': tok = MULT; break;
  1386.     case '/': tok = DIV; break;
  1387.     case '(': parens_nest++; tok = LPAREN; break;
  1388.     case ')':
  1389.         if (--parens_nest < 0) {
  1390.             sprintf (err_msg, "Too many right parens");
  1391.         return (ERR);
  1392.         } else
  1393.         tok = RPAREN;
  1394.         break;
  1395.     case '|':
  1396.         if (*cexpr == '|') { cexpr++; tok = OR; }
  1397.         else { sprintf (err_msg, badop); return (ERR); }
  1398.         break;
  1399.     case '&':
  1400.         if (*cexpr == '&') { cexpr++; tok = AND; }
  1401.         else { sprintf (err_msg, badop); return (ERR); }
  1402.         break;
  1403.     case '=':
  1404.         if (*cexpr == '=') { cexpr++; tok = EQ; }
  1405.         else { sprintf (err_msg, badop); return (ERR); }
  1406.         break;
  1407.     case '!':
  1408.         if (*cexpr == '=') { cexpr++; tok = NE; }
  1409.         else { sprintf (err_msg, badop); return (ERR); }
  1410.         break;
  1411.     case '<':
  1412.         if (*cexpr == '=') { cexpr++; tok = LE; }
  1413.         else tok = LT;
  1414.         break;
  1415.     case '>':
  1416.         if (*cexpr == '=') { cexpr++; tok = GE; }
  1417.         else tok = GT;
  1418.         break;
  1419.     }
  1420.  
  1421.     if (tok != ERR)
  1422.         return (tok);
  1423.  
  1424.     /* not op so check for a constant, variable or function */
  1425.     if (isdigit(c) || c == '.') {
  1426.         if (opxidx > MAX_OPX) {
  1427.         sprintf (err_msg, toomt, MAX_OPX);
  1428.         return (ERR);
  1429.         }
  1430.         opx[opxidx].opu_f = atof (lcexpr);
  1431.         tok = CONST | (opxidx++ << OP_SHIFT);
  1432.         skip_double();
  1433.     } else if (isalpha(c)) {
  1434.         /* check list of functions */
  1435.         if (strncmp (lcexpr, "abs", 3) == 0) {
  1436.         cexpr += 2;
  1437.         tok = ABS;
  1438.         } else {
  1439.         /* not a function, so assume it's a variable */
  1440.         int fld;
  1441.         if (opxidx > MAX_OPX) {
  1442.             sprintf (err_msg, toomt, MAX_OPX);
  1443.             return (ERR);
  1444.         }
  1445.         fld = parse_fieldname ();
  1446.         if (fld < 0) {
  1447.             sprintf (err_msg, "Unknown field");
  1448.             return (ERR);
  1449.         } else {
  1450.             if (flog_add (fld) < 0) { /* register with field logger */
  1451.             sprintf (err_msg, "Sorry; too many fields");
  1452.             return (ERR);
  1453.             }
  1454.             opx[opxidx].opu_fld = fld;
  1455.             tok = VAR | (opxidx++ << OP_SHIFT);
  1456.         }
  1457.         }
  1458.     }
  1459.  
  1460.     return (tok);
  1461. }
  1462.  
  1463. /* move cexpr on past a double.
  1464.  * allow sci notation.
  1465.  * no need to worry about a leading '-' or '+' but allow them after an 'e'.
  1466.  * TODO: this handles all the desired cases, but also admits a bit too much
  1467.  *   such as things like 1eee2...3. geeze; to skip a double right you almost
  1468.  *   have to go ahead and crack it!
  1469.  */
  1470. static
  1471. skip_double()
  1472. {
  1473.     int sawe = 0;    /* so we can allow '-' or '+' right after an 'e' */
  1474.  
  1475.     while (1) {
  1476.         char c = *cexpr;
  1477.         if (isdigit(c) || c=='.' || (sawe && (c=='-' || c=='+'))) {
  1478.         sawe = 0;
  1479.         cexpr++;
  1480.         } else if (c == 'e') {
  1481.         sawe = 1;
  1482.         cexpr++;
  1483.         } else
  1484.         break;
  1485.     }
  1486. }
  1487.  
  1488. /* call this whenever you want to dig out the next (sub)expression.
  1489.  * keep compiling instructions as long as the operators are higher precedence
  1490.  * than prec, then return that "look-ahead" token that wasn't (higher prec).
  1491.  * if error, fill in a message in err_msg[] and return ERR.
  1492.  */
  1493. static
  1494. compile (prec)
  1495. int prec;
  1496. {
  1497.     int expect_binop = 0;    /* set after we have seen any operand.
  1498.                  * used by SUB so it can tell if it really 
  1499.                  * should be taken to be a NEG instead.
  1500.                  */
  1501.     int tok = next_token ();
  1502.  
  1503.         while (1) {
  1504.         int p;
  1505.         if (tok == ERR)
  1506.         return (ERR);
  1507.         if (pc - program >= MAX_PROG) {
  1508.         sprintf (err_msg, "Program is too long");
  1509.         return (ERR);
  1510.         }
  1511.  
  1512.         /* check for special things like functions, constants and parens */
  1513.             switch (tok & OP_MASK) {
  1514.             case HALT: return (tok);
  1515.         case ADD:
  1516.         if (expect_binop)
  1517.             break;    /* procede with binary addition */
  1518.         /* just skip a unary positive(?) */
  1519.         tok = next_token();
  1520.         continue;
  1521.         case SUB:
  1522.         if (expect_binop)
  1523.             break;    /* procede with binary subtract */
  1524.         tok = compile (NEG_PREC);
  1525.         *pc++ = NEG;
  1526.         expect_binop = 1;
  1527.         continue;
  1528.             case ABS: /* other funcs would be handled the same too ... */
  1529.         /* eat up the function parenthesized argument */
  1530.         if (next_token() != LPAREN || compile (0) != RPAREN) {
  1531.             sprintf (err_msg, "Function arglist error");
  1532.             return (ERR);
  1533.         }
  1534.         /* then handled same as ... */
  1535.             case CONST: /* handled same as... */
  1536.         case VAR:
  1537.         *pc++ = tok;
  1538.         tok = next_token();
  1539.         expect_binop = 1;
  1540.         continue;
  1541.             case LPAREN:
  1542.         if (compile (0) != RPAREN) {
  1543.             sprintf (err_msg, "Unmatched left paren");
  1544.             return (ERR);
  1545.         }
  1546.         tok = next_token();
  1547.         expect_binop = 1;
  1548.         continue;
  1549.             case RPAREN:
  1550.         return (RPAREN);
  1551.             }
  1552.  
  1553.         /* everything else is a binary operator */
  1554.         p = precedence[tok];
  1555.             if (p > prec) {
  1556.                 int newtok = compile (p);
  1557.         if (newtok == ERR)
  1558.             return (ERR);
  1559.                 *pc++ = tok;
  1560.         expect_binop = 1;
  1561.                 tok = newtok;
  1562.             } else
  1563.                 return (tok);
  1564.         }
  1565. }
  1566.  
  1567. /* "run" the program[] compiled with compile().
  1568.  * if ok, return 0 and the final result,
  1569.  * else return -1 with a reason why not message in err_msg.
  1570.  */
  1571. static
  1572. execute(result)
  1573. double *result;
  1574. {
  1575.     int instr; 
  1576.  
  1577.     do {
  1578.         instr = *pc++;
  1579.         switch (instr & OP_MASK) {
  1580.         /* put these in numberic order so hopefully even the dumbest
  1581.          * compiler will choose to use a jump table, not a cascade of ifs.
  1582.          */
  1583.         case HALT: break;    /* outer loop will stop us */
  1584.         case ADD:  sp[1] = sp[1] +  sp[0]; sp++; break;
  1585.         case SUB:  sp[1] = sp[1] -  sp[0]; sp++; break;
  1586.         case MULT: sp[1] = sp[1] *  sp[0]; sp++; break;
  1587.         case DIV:  sp[1] = sp[1] /  sp[0]; sp++; break;
  1588.         case AND:  sp[1] = sp[1] && sp[0] ? 1 : 0; sp++; break;
  1589.         case OR:   sp[1] = sp[1] || sp[0] ? 1 : 0; sp++; break;
  1590.         case GT:   sp[1] = sp[1] >  sp[0] ? 1 : 0; sp++; break;
  1591.         case GE:   sp[1] = sp[1] >= sp[0] ? 1 : 0; sp++; break;
  1592.         case EQ:   sp[1] = sp[1] == sp[0] ? 1 : 0; sp++; break;
  1593.         case NE:   sp[1] = sp[1] != sp[0] ? 1 : 0; sp++; break;
  1594.         case LT:   sp[1] = sp[1] <  sp[0] ? 1 : 0; sp++; break;
  1595.         case LE:   sp[1] = sp[1] <= sp[0] ? 1 : 0; sp++; break;
  1596.         case NEG:  *sp = -*sp; break;
  1597.         case CONST: *--sp = opx[instr >> OP_SHIFT].opu_f; break;
  1598.         case VAR:
  1599.         if (flog_get (opx[instr >> OP_SHIFT].opu_fld, --sp) < 0) {
  1600.             sprintf (err_msg, "Bug! VAR field not logged");
  1601.             return (-1);
  1602.         }
  1603.         break;
  1604.         case ABS:  *sp = fabs (*sp); break;
  1605.         default:
  1606.         sprintf (err_msg, "Bug! bad opcode: 0x%x", instr);
  1607.         return (-1);
  1608.         }
  1609.         if (sp < stack) {
  1610.         sprintf (err_msg, "Runtime stack overflow");
  1611.         return (-1);
  1612.         } else if (sp - stack > MAX_STACK) {
  1613.         sprintf (err_msg, "Bug! runtime stack underflow");
  1614.         return (-1);
  1615.         }
  1616.     } while (instr != HALT);
  1617.  
  1618.     /* result should now be on top of stack */
  1619.     if (sp != &stack[MAX_STACK - 1]) {
  1620.         sprintf (err_msg, "Bug! stack has %d items",MAX_STACK-(sp-stack));
  1621.         return (-1);
  1622.     }
  1623.     *result = *sp;
  1624.     return (0);
  1625. }
  1626.  
  1627. static
  1628. isdigit(c)
  1629. char c;
  1630. {
  1631.     return (c >= '0' && c <= '9');
  1632. }
  1633.  
  1634. static
  1635. isalpha (c)
  1636. char c;
  1637. {
  1638.     return ((c >= 'a' && c <= 'z') || (c >=  'A' && c <= 'Z'));
  1639. }
  1640.  
  1641. /* starting with lcexpr pointing at a string expected to be a field name,
  1642.  * return an rcfpack(r,c,0) of the field else -1 if bad.
  1643.  * when return, leave lcexpr alone but move cexpr to just after the name.
  1644.  */
  1645. static
  1646. parse_fieldname ()
  1647. {
  1648.     int r = -1, c = -1;     /* anything illegal */
  1649.     char *fn = lcexpr;    /* likely faster than using the global */
  1650.     char f0, f1;
  1651.     char *dp;
  1652.  
  1653.     /* search for first thing not an alpha char.
  1654.      * leave it in f0 and leave dp pointing to it.
  1655.      */
  1656.     dp = fn;
  1657.     while (isalpha(f0 = *dp))
  1658.         dp++;
  1659.  
  1660.     /* crack the new field name.
  1661.      * when done trying, leave dp pointing at first char just after it.
  1662.      * set r and c if we recognized it.
  1663.      */
  1664.     if (f0 == '.') {
  1665.         /* planet.column pair.
  1666.          * first crack the planet portion (pointed to by fn): set r.
  1667.          * then the second portion (pointed to by dp+1): set c.
  1668.          */
  1669.         char xn[MAXOBJXNM+1];
  1670.  
  1671.         f0 = fn[0];
  1672.         f1 = fn[1];
  1673.  
  1674.         /* if there is an object-x we insist on a perfect match */
  1675.         objx_get ((double *)0, (double *)0, (double *)0, xn);
  1676.         if (xn[0] && f0 == xn[0] && (!xn[1] || f1 == xn[1]))
  1677.         r = R_OBJX;
  1678.         else
  1679.         switch (f0) {
  1680.         case 'j':
  1681.                     r = R_JUPITER;
  1682.             break;
  1683.         case 'm':
  1684.             if (f1 == 'a')      r = R_MARS;
  1685.             else if (f1 == 'e') r = R_MERCURY;
  1686.             else if (f1 == 'o') r = R_MOON;
  1687.             break;
  1688.         case 'n':
  1689.                     r = R_NEPTUNE;
  1690.             break;
  1691.         case 'p':
  1692.                     r = R_PLUTO;
  1693.             break;
  1694.         case 's':
  1695.             if (f1 == 'a')      r = R_SATURN;
  1696.             else if (f1 == 'u') r = R_SUN;
  1697.             break;
  1698.         case 'u':
  1699.                     r = R_URANUS;
  1700.             break;
  1701.         case 'v':
  1702.                     r = R_VENUS;
  1703.             break;
  1704.         }
  1705.  
  1706.         /* now crack the column (stuff after the dp) */
  1707.         dp++;    /* point at good stuff just after the decimal pt */
  1708.         f0 = dp[0];
  1709.         f1 = dp[1];
  1710.         switch (f0) {
  1711.         case 'a':
  1712.         if (f1 == 'l')        c = C_ALT;
  1713.         else if (f1 == 'z')   c = C_AZ;
  1714.         break;
  1715.         case 'd':
  1716.                       c = C_DEC;
  1717.         break;
  1718.         case 'e':
  1719.         if (f1 == 'd')        c = C_EDIST;
  1720.         else if (f1 == 'l')   c = C_ELONG;
  1721.         break;
  1722.         case 'h':
  1723.         if (f1 == 'l') {
  1724.             if (dp[2] == 'a')              c = C_HLAT;
  1725.             else if (dp[2] == 'o')         c = C_HLONG;
  1726.         } else if (f1 == 'r' || f1 == 'u') c = C_TUP;
  1727.         break;
  1728.         case 'j':
  1729.                       c = C_JUPITER;
  1730.         break;
  1731.         case 'm':
  1732.         if (f1 == 'a')        c = C_MARS;
  1733.         else if (f1 == 'e')   c = C_MERCURY;
  1734.         else if (f1 == 'o')   c = C_MOON;
  1735.         break;
  1736.         case 'n':
  1737.                       c = C_NEPTUNE;
  1738.         break;
  1739.         case 'p':
  1740.         if (f1 == 'h')        c = C_PHASE;
  1741.         else if (f1 == 'l')   c = C_PLUTO;
  1742.         break;
  1743.         case 'r':
  1744.         if (f1 == 'a') {
  1745.             if (dp[2] == 'z') c = C_RISEAZ;
  1746.             else           c = C_RA;
  1747.         } else if (f1 == 't') c = C_RISETM;
  1748.         break;
  1749.         case 's':
  1750.         if (f1 == 'a') {
  1751.             if (dp[2] == 'z') c = C_SETAZ;
  1752.             else          c = C_SATURN;
  1753.         } else if (f1 == 'd') c = C_SDIST;
  1754.         else if (f1 == 'i')   c = C_SIZE;
  1755.         else if (f1 == 't')   c = C_SETTM;
  1756.         else if (f1 == 'u')   c = C_SUN;
  1757.         break;
  1758.         case 't':
  1759.         if (f1 == 'a')        c = C_TRANSALT;
  1760.         else if (f1 == 't')   c = C_TRANSTM;
  1761.         break;
  1762.         case 'u':
  1763.                       c = C_URANUS;
  1764.         break;
  1765.         case 'v':
  1766.         if (f1 == 'e')        c = C_VENUS;
  1767.         else if (f1 == 'm')   c = C_MAG;
  1768.         break;
  1769.         }
  1770.  
  1771.         /* now skip dp on past the column stuff */
  1772.         while (isalpha(*dp))
  1773.         dp++;
  1774.     } else {
  1775.         /* no decimal point; some field in the top of the screen */
  1776.         f0 = fn[0];
  1777.         f1 = fn[1];
  1778.         switch (f0) {
  1779.         case 'd':
  1780.         if (f1 == 'a')      r = R_DAWN, c = C_DAWN;
  1781.         else if (f1 == 'u') r = R_DUSK, c = C_DUSK;
  1782.         break;
  1783.         case 'n':
  1784.         r = R_LON, c = C_LON;
  1785.         break;
  1786.         }
  1787.     }
  1788.  
  1789.     cexpr = dp;
  1790.     if (r <= 0 || c <= 0) return (-1);
  1791.     return (rcfpack (r, c, 0));
  1792. }
  1793. xXx
  1794. echo x eq_ecl.c
  1795. cat > eq_ecl.c << 'xXx'
  1796. #include <stdio.h>
  1797. #include <math.h>
  1798. #include "astro.h"
  1799.  
  1800. #define    EQtoECL    1
  1801. #define    ECLtoEQ    (-1)
  1802.  
  1803. /* given the modified Julian date, mjd, and an equitorial ra and dec, each in
  1804.  * radians, find the corresponding geocentric ecliptic latitude, *lat, and
  1805.  * longititude, *lng, also each in radians.
  1806.  * the effects of nutation and the angle of the obliquity are included.
  1807.  */
  1808. eq_ecl (mjd, ra, dec, lat, lng)
  1809. double mjd, ra, dec;
  1810. double *lat, *lng;
  1811. {
  1812.     ecleq_aux (EQtoECL, mjd, ra, dec, lng, lat);
  1813. }
  1814.  
  1815. /* given the modified Julian date, mjd, and a geocentric ecliptic latitude,
  1816.  * *lat, and longititude, *lng, each in radians, find the corresponding
  1817.  * equitorial ra and dec, also each in radians.
  1818.  * the effects of nutation and the angle of the obliquity are included.
  1819.  */
  1820. ecl_eq (mjd, lat, lng, ra, dec)
  1821. double mjd, lat, lng;
  1822. double *ra, *dec;
  1823. {
  1824.     ecleq_aux (ECLtoEQ, mjd, lng, lat, ra, dec);
  1825. }
  1826.  
  1827. static
  1828. ecleq_aux (sw, mjd, x, y, p, q)
  1829. int sw;            /* +1 for eq to ecliptic, -1 for vv. */
  1830. double mjd, x, y;    /* sw==1: x==ra, y==dec.  sw==-1: x==lng, y==lat. */
  1831. double *p, *q;        /* sw==1: p==lng, q==lat. sw==-1: p==ra, q==dec. */
  1832. {
  1833.     static double lastmjd;        /* last mjd calculated */
  1834.     static double seps, ceps;    /* sin and cos of mean obliquity */
  1835.     double sx, cx, sy, cy, ty;
  1836.  
  1837.     if (mjd != lastmjd) {
  1838.         double deps, dpsi, eps;
  1839.         obliquity (mjd, &eps);        /* mean obliquity for date */
  1840. #ifdef NONUTATION
  1841. #define NONUTATION
  1842.         deps = 0;                /* checkout handler did not
  1843.                          * include nutation correction.
  1844.                          */
  1845. #else
  1846.         nutation (mjd, &deps, &dpsi);    /* correctin due to nutation */
  1847. #endif
  1848.         eps += deps;            /* true obliquity for date */
  1849.             seps = sin(eps);
  1850.         ceps = cos(eps);
  1851.         lastmjd = mjd;
  1852.     }
  1853.  
  1854.     sy = sin(y);
  1855.     cy = cos(y);                /* always non-negative */
  1856.         if (fabs(cy)<1e-20) cy = 1e-20;        /* insure > 0 */
  1857.         ty = sy/cy;
  1858.     cx = cos(x);
  1859.     sx = sin(x);
  1860.         *q = asin((sy*ceps)-(cy*seps*sx*sw));
  1861.         *p = atan(((sx*ceps)+(ty*seps*sw))/cx);
  1862.         if (cx<0) *p += PI;        /* account for atan quad ambiguity */
  1863.     range (p, 2*PI);
  1864. }
  1865. xXx
  1866. echo x flog.c
  1867. cat > flog.c << 'xXx'
  1868. /* this is a simple little package to manage the saving and retrieving of
  1869.  * field values, which we call field logging or "flogs". a flog consists of a
  1870.  * field location, ala rcfpack(), and its value as a double. you can reset the
  1871.  * list of flogs, add to and remove from the list of registered fields and log
  1872.  * a field if it has been registered.
  1873.  *
  1874.  * this is used by the plotting and searching facilities of ephem to maintain
  1875.  * the values of the fields that are being plotted or used in search
  1876.  * expressions.
  1877.  *
  1878.  * a field can be in use for more than one
  1879.  * thing at a time (eg, all the X plot values may the same time field, or
  1880.  * searching and plotting might be on at one time using the same field) so
  1881.  * we consider the field to be in use as long a usage count is > 0.
  1882.  */
  1883.  
  1884. #include "screen.h"
  1885.  
  1886. #define    NFLOGS    32
  1887.  
  1888. typedef struct {
  1889.     int fl_usagecnt;    /* number of "users" logging to this field */
  1890.     int fl_fld;        /* an rcfpack(r,c,0) */
  1891.     double fl_val;
  1892. } FLog;
  1893.  
  1894. static FLog flog[NFLOGS];
  1895.  
  1896. /* add fld to the list. if already there, just increment usage count.
  1897.  * return 0 if ok, else -1 if no more room.
  1898.  */
  1899. flog_add (fld)
  1900. int fld;
  1901. {
  1902.     FLog *flp, *unusedflp = 0;
  1903.  
  1904.     /* scan for fld already in list, or find an unused one along the way */
  1905.     for (flp = &flog[NFLOGS]; --flp >= flog; ) {
  1906.         if (flp->fl_usagecnt > 0) {
  1907.         if (flp->fl_fld == fld) {
  1908.             flp->fl_usagecnt++;
  1909.             return (0);
  1910.         }
  1911.         } else
  1912.         unusedflp = flp;
  1913.     }
  1914.     if (unusedflp) {
  1915.         unusedflp->fl_fld = fld;
  1916.         unusedflp->fl_usagecnt = 1;
  1917.         return (0);
  1918.     }
  1919.     return (-1);
  1920. }
  1921.  
  1922. /* decrement usage count for flog for fld. if goes to 0 take it out of list.
  1923.  * ok if not in list i guess...
  1924.  */
  1925. flog_delete (fld)
  1926. int fld;
  1927. {
  1928.     FLog *flp;
  1929.  
  1930.     for (flp = &flog[NFLOGS]; --flp >= flog; )
  1931.         if (flp->fl_fld == fld && flp->fl_usagecnt > 0) {
  1932.         if (--flp->fl_usagecnt <= 0) {
  1933.             flp->fl_usagecnt = 0;
  1934.         }
  1935.         break;
  1936.         }
  1937. }
  1938.  
  1939. /* if plotting or searching is active then
  1940.  * if rcfpack(r,c,0) is in the fld list, set its value to val.
  1941.  * return 0 if ok, else -1 if not in list.
  1942.  */
  1943. flog_log (r, c, val)
  1944. int r, c;
  1945. double val;
  1946. {
  1947.     if (plot_ison() || srch_ison()) {
  1948.         FLog *flp;
  1949.         int fld = rcfpack (r, c, 0);
  1950.         for (flp = &flog[NFLOGS]; --flp >= flog; )
  1951.         if (flp->fl_fld == fld && flp->fl_usagecnt > 0) {
  1952.             flp->fl_val = val;
  1953.             return(0);
  1954.         }
  1955.         return (-1);
  1956.     } else
  1957.         return (0);
  1958. }
  1959.  
  1960. /* search for fld in list. if find it return its value.
  1961.  * return 0 if found it, else -1 if not in list.
  1962.  */
  1963. flog_get (fld, vp)
  1964. int fld;
  1965. double *vp;
  1966. {
  1967.     FLog *flp;
  1968.  
  1969.     for (flp = &flog[NFLOGS]; --flp >= flog; )
  1970.         if (flp->fl_fld == fld && flp->fl_usagecnt > 0) {
  1971.         *vp = flp->fl_val;
  1972.         return (0);
  1973.         }
  1974.     return (-1);
  1975. }
  1976. xXx
  1977. echo x formats.c
  1978. cat > formats.c << 'xXx'
  1979. /* basic formating routines.
  1980.  * all the screen oriented printing should go through here.
  1981.  */
  1982.  
  1983. #include <stdio.h>
  1984. #include <math.h>
  1985. #include "astro.h"
  1986. #include "screen.h"
  1987.  
  1988. /* suppress screen io if this is true, but always flog stuff.
  1989.  */
  1990. static int f_scrnoff;
  1991. f_on ()
  1992. {
  1993.     f_scrnoff = 0;
  1994. }
  1995. f_off ()
  1996. {
  1997.     f_scrnoff = 1;
  1998. }
  1999.  
  2000. /* draw n blanks at the given cursor position.  */
  2001. f_blanks (r, c, n)
  2002. int r, c, n;
  2003. {
  2004.     if (f_scrnoff)
  2005.         return;
  2006.     c_pos (r, c);
  2007.     while (--n >= 0)
  2008.         putchar (' ');
  2009. }
  2010.  
  2011. /* print the given value, v, in "sexadecimal" format at [r,c]
  2012.  * ie, in the form A:m.P, where A is a digits wide, P is p digits.
  2013.  * if p == 0, then no decimal point either.
  2014.  */
  2015. f_sexad (r, c, a, p, mod, v)
  2016. int r, c;
  2017. int a, p;    /* left space, min precision */
  2018. int mod;    /* don't let whole portion get this big */
  2019. double v;
  2020. {
  2021.     char astr[32], str[32];
  2022.     int dec;
  2023.     double frac;
  2024.     int visneg;
  2025.  
  2026.     (void) flog_log (r, c, v);
  2027.  
  2028.     if (f_scrnoff)
  2029.         return;
  2030.  
  2031.     if (v >= 0.0)
  2032.         visneg = 0;
  2033.     else {
  2034.         visneg = 1;
  2035.         v = -v;
  2036.     }
  2037.  
  2038.     dec = v;
  2039.     frac = (v - dec)*60.0;
  2040.     sprintf (str, "59.%.*s5", p, "999999999");
  2041.     if (frac >= atof (str)) {
  2042.         dec += 1;
  2043.         frac = 0.0;
  2044.     }
  2045.     dec %= mod;
  2046.     if (dec == 0 && visneg)
  2047.         strcpy (str, "-0");
  2048.     else
  2049.         sprintf (str, "%d", visneg ? -dec : dec);
  2050.  
  2051.     /* would just do this if Turbo-C 2.0 %?.0f" worked:
  2052.      * sprintf (astr, "%*s:%0*.*f", a, str, p == 0 ? 2 : p+3, p, frac);
  2053.      */
  2054.     if (p == 0)
  2055.         sprintf (astr, "%*s:%02d", a, str, (int)(frac+0.5));
  2056.     else
  2057.         sprintf (astr, "%*s:%0*.*f", a, str, p+3, p, frac);
  2058.     f_string (r, c, astr);
  2059. }
  2060.  
  2061. /* print the given value, t, in sexagesimal format at [r,c]
  2062.  * ie, in the form T:mm:ss, where T is nd digits wide.
  2063.  * N.B. we assume nd >= 2.
  2064.  */
  2065. f_sexag (r, c, nd, t)
  2066. int r, c, nd;
  2067. double t;
  2068. {
  2069.     char tstr[32];
  2070.     int h, m, s;
  2071.     int tisneg;
  2072.     
  2073.     (void) flog_log (r, c, t);
  2074.     if (f_scrnoff)
  2075.         return;
  2076.     dec_sex (t, &h, &m, &s, &tisneg);
  2077.     if (h == 0 && tisneg)
  2078.         sprintf (tstr, "%*s-0:%02d:%02d", nd-2, "", m, s);
  2079.     else
  2080.         sprintf (tstr, "%*d:%02d:%02d", nd, tisneg ? -h : h, m, s);
  2081.     f_string (r, c, tstr);
  2082. }
  2083.  
  2084. /* print angle ra, in radians, in ra hours as hh:mm.m at [r,c]
  2085.  * N.B. we assume ra is >= 0.
  2086.  */
  2087. f_ra (r, c, ra)
  2088. int r, c;
  2089. double ra;
  2090. {
  2091.     f_sexad (r, c, 2, 1, 24, radhr(ra));
  2092. }
  2093.  
  2094. /* print time, t, as hh:mm:ss */
  2095. f_time (r, c, t)
  2096. int r, c;
  2097. double t;
  2098. {
  2099.     f_sexag (r, c, 2, t);
  2100. }
  2101.  
  2102. /* print time, t, as +/-hh:mm:ss (don't show leading +) */
  2103. f_signtime (r, c, t)
  2104. int r, c;
  2105. double t;
  2106. {
  2107.     f_sexag (r, c, 3, t);
  2108. }
  2109.  
  2110. /* print time, t, as hh:mm */
  2111. f_mtime (r, c, t)
  2112. int r, c;
  2113. double t;
  2114. {
  2115.     f_sexad (r, c, 2, 0, 24, t);
  2116. }
  2117.  
  2118. /* print angle, a, in rads, as degress at [r,c] in form ddd:mm */
  2119. f_angle(r, c, a)
  2120. int r, c;
  2121. double a;
  2122. {
  2123.     f_sexad (r, c, 3, 0, 360, raddeg(a));
  2124. }
  2125.  
  2126. /* print angle, a, in rads, as degress at [r,c] in form dddd:mm:ss */
  2127. f_gangle(r, c, a)
  2128. int r, c;
  2129. double a;
  2130. {
  2131.     f_sexag (r, c, 4, raddeg(a));
  2132. }
  2133.  
  2134. /* print the given modified Julian date, jd, as the starting date at [r,c]
  2135.  * in the form mm/dd/yyyy.
  2136.  */
  2137. f_date (r, c, jd)
  2138. int r, c;
  2139. double jd;
  2140. {
  2141.     char dstr[32];
  2142.     int m, y;
  2143.     double d, tmp;
  2144.  
  2145.     /* shadow to the plot subsystem as years. */
  2146.     mjd_year (jd, &tmp);
  2147.     (void) flog_log (r, c, tmp);
  2148.     if (f_scrnoff)
  2149.         return;
  2150.  
  2151.     mjd_cal (jd, &m, &d, &y);
  2152.  
  2153.     sprintf (dstr, "%2d/%02d/%04d", m, (int)(d), y);
  2154.     f_string (r, c, dstr);
  2155. }
  2156.  
  2157. f_char (row, col, c)
  2158. int row, col;
  2159. char c;
  2160. {
  2161.     if (f_scrnoff)
  2162.         return;
  2163.     c_pos (row, col);
  2164.     putchar (c);
  2165. }
  2166.  
  2167. f_string (r, c, s)
  2168. int r, c;
  2169. char *s;
  2170. {
  2171.     if (f_scrnoff)
  2172.         return;
  2173.     c_pos (r, c);
  2174.     fputs (s, stdout);
  2175. }
  2176.  
  2177. f_double (r, c, fmt, f)
  2178. int r, c;
  2179. char *fmt;
  2180. double f;
  2181. {
  2182.     char str[80];
  2183.     (void) flog_log (r, c, f);
  2184.     sprintf (str, fmt, f);
  2185.     f_string (r, c, str);
  2186. }
  2187.  
  2188. /* print prompt line */
  2189. f_prompt (p)
  2190. char *p;
  2191. {
  2192.     c_pos (R_PROMPT, C_PROMPT);
  2193.     c_eol ();
  2194.     c_pos (R_PROMPT, C_PROMPT);
  2195.     fputs (p, stdout);
  2196. }
  2197.  
  2198. /* print a message and wait for op to hit any key */
  2199. f_msg (m)
  2200. char *m;
  2201. {
  2202.     f_prompt (m);
  2203.     (void) read_char();
  2204. }
  2205.  
  2206. /* crack a line of the form X?X?X into its components,
  2207.  *   where X is an integer and ? can be any character except '0-9' or '-',
  2208.  *   such as ':' or '/'.
  2209.  * only change those fields that are specified:
  2210.  *   eg:  ::10    only changes *s
  2211.  *        10    only changes *d
  2212.  *        10:0  changes *d and *m
  2213.  * if see '-' anywhere, first non-zero component will be made negative.
  2214.  */
  2215. f_sscansex (bp, d, m, s)
  2216. char *bp;
  2217. int *d, *m, *s;
  2218. {
  2219.     char c;
  2220.     int *p = d;
  2221.     int *nonzp = 0;
  2222.     int sawneg = 0;
  2223.     int innum = 0;
  2224.  
  2225.     while (c = *bp++)
  2226.         if (c >= '0' && c <= '9') {
  2227.         if (!innum) {
  2228.             *p = 0;
  2229.             innum = 1;
  2230.         }
  2231.         *p = *p*10 + (c - '0');
  2232.         if (*p && !nonzp)
  2233.             nonzp = p;
  2234.         } else if (c == '-') {
  2235.         sawneg = 1;
  2236.         } else if (c != ' ') {
  2237.         /* advance to next component */
  2238.         p = (p == d) ? m : s;
  2239.         innum = 0;
  2240.         }
  2241.  
  2242.     if (sawneg && nonzp)
  2243.         *nonzp = -*nonzp;
  2244. }
  2245.  
  2246. /* just like dec_sex() but makes the first non-zero element negative if
  2247.  * x is negative (instead of returning a sign flag).
  2248.  */
  2249. f_dec_sexsign (x, h, m, s)
  2250. double x;
  2251. int *h, *m, *s;
  2252. {
  2253.     int n;
  2254.     dec_sex (x, h, m, s, &n);
  2255.     if (n) {
  2256.         if (*h)
  2257.         *h = -*h;
  2258.         else if (*m)
  2259.         *m = -*m;
  2260.         else
  2261.         *s = -*s;
  2262.     }
  2263. }
  2264. xXx
  2265.