home *** CD-ROM | disk | FTP | other *** search
/ rtsi.com / 2014.01.www.rtsi.com.tar / www.rtsi.com / OS9 / OSK / EFFO / pd8.lzh / SRC / objx.c < prev    next >
Text File  |  1990-04-13  |  21KB  |  851 lines
  1. /* functions to save the user-definable objects, referred to as "x" and "y".
  2.  * this way, once defined, the objects can be quieried for position just like
  3.  * the other bodies, with obj_cir(). 
  4.  */
  5.  
  6. #include <stdio.h>
  7. #include <math.h>
  8. #include "astro.h"
  9. #include "circum.h"
  10. #include "screen.h"
  11.  
  12. extern char *strcat(), *strcpy(), *strncpy(), *getenv();
  13.  
  14. static char *dbfile;            /* !0 if set by -d option */
  15. static char dbfdef[] = "ephem.db";     /* default database file name */
  16.  
  17. /* structures to describe objects of various types.
  18.  */
  19. typedef struct {
  20.     double f_ra;    /* ra, rads, at given epoch */
  21.     double f_dec;    /* dec, rads, at given epoch */
  22.     double f_mag;    /* visual magnitude */
  23.     double f_epoch;    /* the given epoch, as an mjd */
  24. } ObjF;            /* fixed object */
  25. typedef struct {
  26.     double e_inc;    /* inclination, degrees */
  27.     double e_Om;    /* longitude of ascending node, degrees */
  28.     double e_om;    /* argument of perihelion, degress */
  29.     double e_a;        /* mean distance, aka, semi-maj axis, in AU */
  30.     double e_n;        /* daily motion, degrees/day */
  31.     double e_e;        /* eccentricity */
  32.     double e_M;        /* mean anomaly, ie, degrees from perihelion at... */
  33.     double e_cepoch;    /* epoch date (M reference), as an mjd */
  34.     double e_epoch;    /* equinox year (inc/Om/om reference), as an mjd */
  35.     double e_H, e_G;    /* magnitude model coefficients */
  36. } ObjE;            /* object in heliocentric elliptical orbit */
  37. typedef struct {
  38.     double p_ep;    /* epoch of perihelion, as an mjd */
  39.     double p_inc;    /* inclination, rads */
  40.     double p_qp;    /* perihelion distance, AU */
  41.     double p_ap;    /* argument of perihelion, rads. */
  42.     double p_om;    /* longitude of ascending node, rads */
  43.     double p_epoch;    /* reference epoch, as an mjd */
  44.     double p_g, p_k;    /* magnitude model coefficients */
  45. } ObjP;            /* object in heliocentric parabolic trajectory */
  46.  
  47. typedef struct {
  48.     int  o_type;    /* see flags, below */
  49.     int  o_on;        /* !=0 while object is active */
  50.     ObjF o_f;
  51.     ObjE o_e;
  52.     ObjP o_p;
  53. } Obj;
  54. #define    FIXED        0    /* just simple ra/dec object */
  55. #define    ELLIPTICAL    1    /* elliptical orbital elements */
  56. #define    PARABOLIC    2    /* parabolic " */
  57.  
  58. static Obj objx;
  59. static Obj objy;
  60.  
  61. /* run when Objx or y is picked from menu.
  62.  * we tell which by the planet code.
  63.  * let op define object and turn it on and off.
  64.  */
  65. obj_setup(p)
  66. int p;
  67. {
  68.     static char *pr[5] = {
  69.         "Fixed", "Elliptical", "Parabolic", "Lookup"
  70.     };
  71.     int f;
  72.     Obj *op;
  73.  
  74.     op = (p == OBJX) ? &objx : &objy;
  75.  
  76.     rechk:
  77.     /* start pointing at selection for current object type */
  78.     switch (op->o_type) {
  79.     case FIXED: f = 0; break;
  80.     case ELLIPTICAL: f = 1; break;
  81.     case PARABOLIC: f = 2; break;
  82.     }
  83.  
  84.     ask:
  85.     pr[4] = op->o_on ? "On" : "Off";
  86.     switch (f = popup (pr, f, 5)) {
  87.     case 0: obj_dfixed(op, (char**)0); goto ask;
  88.     case 1: obj_delliptical(op, (char**)0); goto ask;
  89.     case 2: obj_dparabolic(op, (char**)0); goto ask;
  90.     case 3: obj_filelookup (p, (char *)0); goto rechk;
  91.     case 4: op->o_on ^= 1; break;
  92.     }
  93. }
  94.  
  95. /* turn "on" or "off" but don't forget facts about object the object.
  96.  */
  97. obj_on (p)
  98. int p;
  99. {
  100.     if (p == OBJX)
  101.         objx.o_on = 1;
  102.     else
  103.         objy.o_on = 1;
  104. }
  105. obj_off (p)
  106. int p;
  107. {
  108.     if (p == OBJX)
  109.         objx.o_on = 0;
  110.     else
  111.         objy.o_on = 0;
  112. }
  113.  
  114. /* return true if object is now on, else 0.
  115.  */
  116. obj_ison(p)
  117. int p;
  118. {
  119.     return ((p == OBJX) ? objx.o_on : objy.o_on);
  120. }
  121.  
  122. /* set an alternate database file name.
  123.  * N.B. we assume the storage pointed to by name is permanent.
  124.  */
  125. obj_setdbfilename (name)
  126. char *name;
  127. {
  128.     dbfile = name;
  129. }
  130.  
  131. /* fill in info about object x or y.
  132.  * most arguments and conditions are the same as for plans().
  133.  * only difference is that mag is already apparent, not absolute magnitude.
  134.  * this is called by body_cir() for object x and y just like plans() is called
  135.  * for the planets.
  136.  */
  137. obj_cir (jd, p, lpd0, psi0, rp0, rho0, lam, bet, mag)
  138. double jd;    /* mjd now */
  139. int p;        /* OBJX or OBJY */
  140. double *lpd0;    /* heliocentric longitude, or NOHELIO  */
  141. double *psi0;    /* heliocentric latitude, or 0 if *lpd0 set to NOHELIO */
  142. double *rp0;    /* distance from the sun, or 0 */
  143. double *rho0;    /* true distance from the Earth, or 0 */
  144. double *lam;    /* apparent geocentric ecliptic longitude */
  145. double *bet;    /* apparent geocentric ecliptic latitude */
  146. double *mag;    /* APPARENT magnitude */
  147. {
  148.     Obj *op = (p == OBJX) ? &objx : &objy;
  149.  
  150.     switch (op->o_type) {
  151.     case FIXED: {
  152.         double xr, xd;
  153.         xr = op->o_f.f_ra;
  154.         xd = op->o_f.f_dec;
  155.         if (op->o_f.f_epoch != jd)
  156.         precess (op->o_f.f_epoch, jd, &xr, &xd);
  157.         eq_ecl (jd, xr, xd, bet, lam);
  158.  
  159.         *lpd0 = NOHELIO;
  160.         *psi0 = *rp0 = *rho0 = 0.0;
  161.         *mag = op->o_f.f_mag;
  162.         }
  163.         break;
  164.     case PARABOLIC: {
  165.         double inc, ap, om;
  166.         double lpd, psi, rp, rho;
  167.         double dt;
  168.         int pass;
  169.         /* two passes to correct lam and bet for light travel time. */
  170.         dt = 0.0;
  171.         for (pass = 0; pass < 2; pass++) {
  172.         reduce_elements (op->o_p.p_epoch, jd-dt, op->o_p.p_inc,
  173.                     op->o_p.p_ap, op->o_p.p_om, &inc, &ap, &om);
  174.         comet (jd-dt, op->o_p.p_ep, inc, ap, op->o_p.p_qp, om,
  175.                     &lpd, &psi, &rp, &rho, lam, bet);
  176.         if (pass == 0) {
  177.             *lpd0 = lpd;
  178.             *psi0 = psi;
  179.             *rp0 = rp;
  180.             *rho0 = rho;
  181.         }
  182.         dt = rho*5.775518e-3;    /* au to light-days */
  183.         }
  184.         *mag = op->o_p.p_g + 5*log10(*rho0) + 2.5*op->o_p.p_k*log10(*rp0);
  185.         }
  186.         break;
  187.     case ELLIPTICAL: {
  188.         /* this is basically the same code as pelement() and plans()
  189.          * combined and simplified for the special case of osculating
  190.          * (unperturbed) elements.
  191.          * inputs have been changed to match the Astronomical Almanac.
  192.          * we have added reduction of elements using reduce_elements().
  193.          */
  194.         double dt, lg, lsn, rsn;
  195.         double nu, ea;
  196.         double ma, rp, lo, slo, clo;
  197.         double inc, psi, spsi, cpsi;
  198.         double y, lpd, rpd, ll, rho, sll, cll;
  199.         double om;        /* arg of perihelion */
  200.         double Om;        /* long of ascending node. */
  201.         double psi_t, Psi_1, Psi_2, beta;
  202.         double e;
  203.         int pass;
  204.  
  205.         dt = 0;
  206.         sunpos (jd, &lsn, &rsn);
  207.         lg = lsn + PI;
  208.         e = op->o_e.e_e;
  209.  
  210.         for (pass = 0; pass < 2; pass++) {
  211.  
  212.         reduce_elements (op->o_e.e_epoch, jd-dt, degrad(op->o_e.e_inc),
  213.                 degrad (op->o_e.e_om), degrad (op->o_e.e_Om),
  214.                 &inc, &om, &Om);
  215.  
  216.         ma = degrad (op->o_e.e_M
  217.                 + (jd - op->o_e.e_cepoch - dt) * op->o_e.e_n);
  218.         anomaly (ma, e, &nu, &ea);
  219.         rp= op->o_e.e_a * (1-e*e) / (1+e*cos(nu));
  220.         lo = nu + om;
  221.         slo = sin(lo);
  222.         clo = cos(lo);
  223.         spsi = slo*sin(inc);
  224.         y = slo*cos(inc);
  225.         psi = asin(spsi);
  226.         lpd = atan(y/clo)+Om;
  227.         if (clo<0) lpd += PI;
  228.         range (&lpd, 2*PI);
  229.         cpsi = cos(psi);
  230.         rpd = rp*cpsi;
  231.         ll = lpd-lg;
  232.         rho = sqrt(rsn*rsn+rp*rp-2*rsn*rp*cpsi*cos(ll));
  233.         dt = rho*5.775518e-3;    /* light travel time, in days */
  234.         if (pass == 0) {
  235.             *lpd0 = lpd;
  236.             *psi0 = psi;
  237.             *rp0 = rp;
  238.             *rho0 = rho;
  239.         }
  240.         }
  241.  
  242.         sll = sin(ll);
  243.         cll = cos(ll);
  244.         *lam = atan(rsn*sll/(rpd-rsn*cll))+lpd;
  245.         /* *lam = atan(-1*rpd*sll/(rsn-rpd*cll))+lg+PI; */
  246.         range (lam, 2*PI);
  247.         *bet = atan(rpd*spsi*sin(*lam-lpd)/(cpsi*rsn*sll));
  248.  
  249.         beta = acos((rp*rp + rho*rho - rsn*rsn)/ (2*rp*rho));
  250.         psi_t = exp(log(tan(beta/2.0))*0.63);
  251.         Psi_1 = exp(-3.33*psi_t);
  252.         psi_t = exp(log(tan(beta/2.0))*1.22);
  253.         Psi_2 = exp(-1.87*psi_t);
  254.         *mag = op->o_e.e_H + 5.0*log10(rp*rho)
  255.             - 2.5*log10((1-op->o_e.e_G)*Psi_1 + op->o_e.e_G*Psi_2);
  256.  
  257.         }
  258.         break;
  259.     }
  260. }
  261.  
  262. /* define obj based on the ephem.cfg line, s.
  263.  * the "OBJX " or "OBJY " keyword is already skipped and used to set p.
  264.  * format: type,[other fields, as per corresponding ObjX typedef]
  265.  * N.B. we replace all ',' with '\0' within s IN PLACE.
  266.  * return 0 if ok, else print reason why not with f_msg() and return -1.
  267.  */
  268. obj_define (p, s)
  269. int p;    /* OBJX or OBJY */
  270. char *s;
  271. {
  272. #define    MAXARGS    20
  273.     char *av[MAXARGS];    /* point to each field for easy reference */
  274.     char c;
  275.     int ac;
  276.     Obj *op = (p == OBJX) ? &objx : &objy;
  277.  
  278.     /* parse into comma separated fields */
  279.     ac = 0;
  280.     av[0] = s;
  281.     do {
  282.         c = *s++;
  283.         if (c == ',' || c == '\0') {
  284.         s[-1] = '\0';
  285.         av[++ac] = s;
  286.         }
  287.     } while (c);
  288.  
  289.     if (ac < 1) {
  290.         f_msg ("No type for OBJX");
  291.         return (-1);
  292.     }
  293.  
  294.     switch (av[0][0]) {
  295.     case 'f':
  296.         if (ac != 5) {
  297.         f_msg ("Need ra,dec,mag,epoch for \"fixed\" OBJ");
  298.         return (-1);
  299.         }
  300.         obj_dfixed (op, av+1);
  301.         break;
  302.     case 'p':
  303.         if (ac != 9) {
  304.         f_msg ("Need ep,inc,ap,qp,om,epoch,abs,lum for \"parabolic\" OBJ");
  305.         return (-1);
  306.         }
  307.         obj_dparabolic (op, av+1);
  308.         break;
  309.     case 'e':
  310.         if (ac != 12) {
  311.         f_msg ("Need inc,lan,aop,md,dm,ecc,ma,cep,ep,h,g for \"elliptical\" OBJ");
  312.         return (-1);
  313.         }
  314.         obj_delliptical (op, av+1);
  315.         break;
  316.     default:
  317.         f_msg ("Unknown OBJ type");
  318.         return (-1);
  319.     }
  320.  
  321.     return (0);
  322. }
  323.  
  324. /* search through an ephem database file for an entry and use it to define
  325.  *   either OBJX or OBJY (as set by p).
  326.  * if a name, np, is not set then we ask for it.
  327.  * if -d was used use it; else if EPHEMDB env set use it, else use default.
  328.  */
  329. obj_filelookup (p, np)
  330. int p;            /* OBJX or OBJY */
  331. char *np;
  332. {
  333.     FILE *fp;
  334.     char *fn;
  335.     int nl;
  336.     char buf[160];
  337.     char name[64];
  338.     int found;
  339.  
  340.     /* open the database file */
  341.     if (dbfile)
  342.         fn = dbfile;
  343.     else {
  344.         fn = getenv ("EPHEMDB");
  345.         if (!fn)
  346.         fn = dbfdef;
  347.     }
  348.     fp = fopen (fn, "r");
  349. #ifdef OSK
  350.     if (!fp)
  351.         fp = fopen("/dd/sys/ephem.db","r");
  352. #endif
  353.     if (!fp) {
  354.         sprintf (buf, "Can not open database file %s", fn);
  355.         f_msg(buf);
  356.         return;
  357.     }
  358.  
  359.     /* set up object name in name with a trailing ',' */
  360.     if (np) {
  361.         (void) strncpy (name, np, sizeof(name)-2);
  362.         name[sizeof(name)-2] = '\0';    /* insure trailing '\0' */
  363.     } else  {
  364.         f_prompt ("Object name: ");
  365.         if (read_line (name, sizeof(name)-2) <= 0)
  366.         return;
  367.     }
  368.     (void) strcat (name, ",");
  369.  
  370.     /* search for line beginning with name followed by comma.
  371.      * then rest of line is the info.
  372.      */
  373.     nl = strlen (name);
  374.     found = 0;
  375.     while (fgets (buf, sizeof(buf), fp))
  376.         if (strncmp (name, buf, nl) == 0) {
  377.         found = 1;
  378.         break;
  379.         }
  380.     fclose (fp);
  381.  
  382.     if (found)
  383.         (void) obj_define (p, buf+nl);
  384.     else {
  385.         sprintf (buf, "%s not found", name);
  386.         f_msg (buf);
  387.     }
  388. }
  389.  
  390. /* define a fixed object.
  391.  * args in av, in order, are ra, dec, magnitude and reference epoch.
  392.  * if av then it is a list of strings to use for each parameter, else must
  393.  * ask for each. the av option is for cracking the ephem.cfg line.
  394.  * if asking show current settings and leave unchanged if hit RETURN.
  395.  * END aborts without making any more changes.
  396.  * o_type is set to FIXED.
  397.  * N.B. we don't error check av in any way, not even for length.
  398.  */
  399. static
  400. obj_dfixed (op, av)
  401. Obj *op;
  402. char *av[];
  403. {
  404.     char buf[NC];
  405.     char *bp;
  406.     int sts;
  407.  
  408.     op->o_type = FIXED;
  409.  
  410.     if (av) {
  411.         bp = av[0];
  412.         sts = 1;
  413.     } else {
  414.         static char p[] = "RA (h:m:s): (";
  415.         f_prompt (p);
  416.         f_ra (R_PROMPT, C_PROMPT+sizeof(p)-1, op->o_f.f_ra);
  417.         printf (") ");
  418.         sts = read_line (buf, 8+1);
  419.         if (sts < 0)
  420.         return;
  421.         bp = buf;
  422.     }
  423.     if (sts > 0) {
  424.         int h, m, s;
  425.         f_dec_sexsign (radhr(op->o_f.f_ra), &h, &m, &s);
  426.         f_sscansex (bp, &h, &m, &s);
  427.         sex_dec (h, m, s, &op->o_f.f_ra);
  428.         op->o_f.f_ra = hrrad(op->o_f.f_ra);
  429.     }
  430.  
  431.     if (av) {
  432.         bp = av[1];
  433.         sts = 1;
  434.     } else {
  435.         static char p[] = "Dec (d:m:s): (";
  436.         f_prompt (p);
  437.         f_angle (R_PROMPT, C_PROMPT+sizeof(p)-1, op->o_f.f_dec);
  438.         printf (") ");
  439.         sts = read_line (buf, 9+1);
  440.         if (sts < 0)
  441.         return;
  442.         bp = buf;
  443.     }
  444.     if (sts > 0) {
  445.         int dg, m, s;
  446.         f_dec_sexsign (raddeg(op->o_f.f_dec), &dg, &m, &s);
  447.         f_sscansex (bp, &dg, &m, &s);
  448.         sex_dec (dg, m, s, &op->o_f.f_dec);
  449.         op->o_f.f_dec = degrad(op->o_f.f_dec);
  450.     }
  451.  
  452.     if (av) {
  453.         bp = av[2];
  454.         sts = 1;
  455.     } else {
  456.         static char p[] = "Magnitude: ";
  457.         f_prompt (p);
  458.         f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", op->o_f.f_mag);
  459.         sts = read_line (buf, 8+1);
  460.         if (sts < 0)
  461.         return;
  462.         bp = buf;
  463.     }
  464.     if (sts > 0)
  465.         op->o_f.f_mag = atof (bp);
  466.  
  467.     if (av) {
  468.         bp = av[3];
  469.         sts = 1;
  470.     } else {
  471.         static char p[] = "Reference epoch (UT Date, m/d.d/y or year.d): ";
  472.         double y;
  473.         f_prompt (p);
  474.         mjd_year (op->o_f.f_epoch, &y);
  475.         printf ("(%g) ", y);
  476.         sts = read_line (buf, 20);
  477.         if (sts < 0)
  478.         return;
  479.         bp = buf;
  480.     }
  481.     if (sts > 0)
  482.         crack_year (bp, &op->o_f.f_epoch);
  483. }
  484.  
  485. /* define an object in an elliptical heliocentric orbit.
  486.  * 11 args in av, in order, are inclination, longitude of ascending node,
  487.  *   argument of perihelion, mean distance (aka semi-major axis), daily
  488.  *   motion, eccentricity, mean anomaly (ie, degrees from perihelion),
  489.  *   epoch date (ie, time of the mean anomaly value), equinox year
  490.  *   (ie, time of inc/lon/aop) absolute visual magnitude (H) and magnitude
  491.  *   slope parameter (G).
  492.  * if av then it is a list of strings to use for each parameter, else must
  493.  * ask for each. the av option is for cracking the ephem.cfg line.
  494.  * if asking show current settings and leave unchanged if hit RETURN.
  495.  * END aborts without making any more changes.
  496.  * o_type is set to ELLIPTICAL.
  497.  * N.B. we don't error check av in any way, not even for length.
  498.  */
  499. static
  500. obj_delliptical(op, av)
  501. Obj *op;
  502. char *av[];
  503. {
  504.     char buf[NC];
  505.     char *bp;
  506.     int sts;
  507.  
  508.     op->o_type = ELLIPTICAL;
  509.  
  510.     if (av) {
  511.         bp = av[0];
  512.         sts = 1;
  513.     } else {
  514.         static char p[] = "Inclination (degs):";
  515.         f_prompt (p);
  516.         f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", op->o_e.e_inc);
  517.         sts = read_line (buf, 8+1);
  518.         if (sts < 0)
  519.         return;
  520.         bp = buf;
  521.     }
  522.     if (sts > 0)
  523.         op->o_e.e_inc = atof(bp);
  524.  
  525.     if (av) {
  526.         bp = av[1];
  527.         sts = 1;
  528.     } else {
  529.         static char p[] = "Longitude of ascending node (degs):";
  530.         f_prompt (p);
  531.         f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", op->o_e.e_Om);
  532.         sts = read_line (buf, 8+1);
  533.         if (sts < 0)
  534.         return;
  535.         bp = buf;
  536.     }
  537.     if (sts > 0)
  538.         op->o_e.e_Om = atof(bp);
  539.  
  540.     if (av) {
  541.         bp = av[2];
  542.         sts = 1;
  543.     } else {
  544.         static char p[] = "Argument of Perihelion (degs):";
  545.         f_prompt (p);
  546.         f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", op->o_e.e_om);
  547.         sts = read_line (buf, 8+1);
  548.         if (sts < 0)
  549.         return;
  550.         bp = buf;
  551.     }
  552.     if (sts > 0)
  553.         op->o_e.e_om = atof(bp);
  554.  
  555.     if (av) {
  556.         bp = av[3];
  557.         sts = 1;
  558.     } else {
  559.         static char p[] = "Mean distance (AU):";
  560.         f_prompt (p);
  561.         f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", op->o_e.e_a);
  562.         sts = read_line (buf, 8+1);
  563.         if (sts < 0)
  564.         return;
  565.         bp = buf;
  566.     }
  567.     if (sts > 0)
  568.         op->o_e.e_a = atof(bp);
  569.  
  570.     if (av) {
  571.         bp = av[4];
  572.         sts = 1;
  573.     } else {
  574.         static char p[] = "Daily motion (degs/day):";
  575.         f_prompt (p);
  576.         f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", op->o_e.e_n);
  577.         sts = read_line (buf, 8+1);
  578.         if (sts < 0)
  579.         return;
  580.         bp = buf;
  581.     }
  582.     if (sts > 0)
  583.         op->o_e.e_n = atof(bp);
  584.  
  585.     if (av) {
  586.         bp = av[5];
  587.         sts = 1;
  588.     } else {
  589.         static char p[] = "Eccentricity:";
  590.         f_prompt (p);
  591.         f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", op->o_e.e_e);
  592.         sts = read_line (buf, 8+1);
  593.         if (sts < 0)
  594.         return;
  595.         bp = buf;
  596.     }
  597.     if (sts > 0)
  598.         op->o_e.e_e = atof(bp);
  599.  
  600.     if (av) {
  601.         bp = av[6];
  602.         sts = 1;
  603.     } else {
  604.         static char p[] = "Mean anomaly (degs):";
  605.         f_prompt (p);
  606.         f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", op->o_e.e_M);
  607.         sts = read_line (buf, 8+1);
  608.         if (sts < 0)
  609.         return;
  610.         bp = buf;
  611.     }
  612.     if (sts > 0)
  613.         op->o_e.e_M = atof(bp);
  614.  
  615.     if (av) {
  616.         bp = av[7];
  617.         sts = 1;
  618.     } else {
  619.         static char p[] = 
  620.         "Epoch date (UT Date, m/d.d/y or year.d): ";
  621.         int m, y;
  622.         double d;
  623.         f_prompt (p);
  624.         mjd_cal (op->o_e.e_cepoch, &m, &d, &y);
  625.         printf ("(%d/%g/%d) ", m, d, y);
  626.         sts = read_line (buf, 20);
  627.         if (sts < 0)
  628.         return;
  629.         bp = buf;
  630.     }
  631.     if (sts > 0)
  632.         crack_year (bp, &op->o_e.e_cepoch);
  633.  
  634.     if (av) {
  635.         bp = av[8];
  636.         sts = 1;
  637.     } else {
  638.         static char p[] = "Equinox year (UT Date, m/d.d/y or year.d): ";
  639.         double y;
  640.         f_prompt (p);
  641.         mjd_year (op->o_e.e_epoch, &y);
  642.         printf ("(%g) ", y);
  643.         sts = read_line (buf, 20);
  644.         if (sts < 0)
  645.         return;
  646.         bp = buf;
  647.     }
  648.     if (sts > 0)
  649.         crack_year (bp, &op->o_e.e_epoch);
  650.  
  651.  
  652.     if (av) {
  653.         bp = av[9];
  654.         sts = 1;
  655.     } else {
  656.         static char p[] = "Absolute magnitude (H):";
  657.         f_prompt (p);
  658.         f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", op->o_e.e_H);
  659.         sts = read_line (buf, 8+1);
  660.         if (sts < 0)
  661.         return;
  662.         bp = buf;
  663.     }
  664.     if (sts > 0)
  665.         op->o_e.e_H = atof(bp);
  666.  
  667.     if (av) {
  668.         bp = av[10];
  669.         sts = 1;
  670.     } else {
  671.         static char p[] = "Magnitude slope parameter (G):";
  672.         f_prompt (p);
  673.         f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", op->o_e.e_G);
  674.         sts = read_line (buf, 8+1);
  675.         if (sts < 0)
  676.         return;
  677.         bp = buf;
  678.     }
  679.     if (sts > 0)
  680.         op->o_e.e_G = atof(bp);
  681. }
  682.  
  683. /* define an object in heliocentric parabolic orbit.
  684.  * 10 args in av, in order, are epoch of perihelion, inclination, argument of
  685.  *   perihelion, perihelion distance, longitude of ascending node, reference
  686.  *   epoch, absolute magnitude and luminosity index.
  687.  * if av then it is a list of strings to use for each parameter, else must
  688.  * ask for each. the av option is for cracking the ephem.cfg line.
  689.  * if asking show current settings and leave unchanged if hit RETURN.
  690.  * END aborts without making any more changes.
  691.  * o_type is set to PARABOLIC.
  692.  * N.B. we don't error check av in any way, not even for length.
  693.  */
  694. static
  695. obj_dparabolic(op, av)
  696. Obj *op;
  697. char *av[];
  698. {
  699.     char buf[NC];
  700.     char *bp;
  701.     int sts;
  702.  
  703.     op->o_type = PARABOLIC;
  704.  
  705.     if (av) {
  706.         bp = av[0];
  707.         sts = 1;
  708.     } else {
  709.         static char p[] =
  710.         "Epoch of perihelion (UT Date, m/d.d/y or year.d): ";
  711.         int m, y;
  712.         double d;
  713.         f_prompt (p);
  714.         mjd_cal (op->o_p.p_ep, &m, &d, &y);
  715.         printf ("(%d/%g/%d) ", m, d, y);
  716.         sts = read_line (buf, 20);
  717.         if (sts < 0)
  718.         return;
  719.         bp = buf;
  720.     }
  721.     if (sts > 0)
  722.         crack_year (bp, &op->o_p.p_ep);
  723.  
  724.     if (av) {
  725.         bp = av[1];
  726.         sts = 1;
  727.     } else {
  728.         static char p[] = "Inclination (degs):";
  729.         f_prompt (p);
  730.         f_double(R_PROMPT,C_PROMPT+sizeof(p),"(%g) ",raddeg(op->o_p.p_inc));
  731.         sts = read_line (buf, 8+1);
  732.         if (sts < 0)
  733.         return;
  734.         bp = buf;
  735.     }
  736.     if (sts > 0)
  737.         op->o_p.p_inc = degrad(atof(bp));
  738.  
  739.     if (av) {
  740.         bp = av[2];
  741.         sts = 1;
  742.     } else {
  743.         static char p[] = "Argument of perihelion (degs):";
  744.         f_prompt (p);
  745.         f_double(R_PROMPT,C_PROMPT+sizeof(p),"(%g) ",raddeg(op->o_p.p_ap));
  746.         sts = read_line (buf, 8+1);
  747.         if (sts < 0)
  748.         return;
  749.         bp = buf;
  750.     }
  751.     if (sts > 0)
  752.         op->o_p.p_ap = degrad(atof(bp));
  753.  
  754.     if (av) {
  755.         bp = av[3];
  756.         sts = 1;
  757.     } else {
  758.         static char p[] = "Perihelion distance (AU):";
  759.         f_prompt (p);
  760.         f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", op->o_p.p_qp);
  761.         sts = read_line (buf, 8+1);
  762.         if (sts < 0)
  763.         return;
  764.         bp = buf;
  765.     }
  766.     if (sts > 0)
  767.         op->o_p.p_qp = atof(bp);
  768.  
  769.     if (av) {
  770.         bp = av[4];
  771.         sts = 1;
  772.     } else {
  773.         static char p[] = "Longitude of ascending node (degs):";
  774.         f_prompt (p);
  775.         f_double(R_PROMPT,C_PROMPT+sizeof(p),"(%g) ",raddeg(op->o_p.p_om));
  776.         sts = read_line (buf, 8+1);
  777.         if (sts < 0)
  778.         return;
  779.         bp = buf;
  780.     }
  781.     if (sts > 0)
  782.         op->o_p.p_om = degrad(atof(bp));
  783.  
  784.     if (av) {
  785.         bp = av[5];
  786.         sts = 1;
  787.     } else {
  788.         static char p[] = "Reference epoch (UT Date, m/d.d/y or year.d): ";
  789.         double y;
  790.         f_prompt (p);
  791.         mjd_year (op->o_p.p_epoch, &y);
  792.         printf ("(%g) ", y);
  793.         sts = read_line (buf, 20);
  794.         if (sts < 0)
  795.         return;
  796.         bp = buf;
  797.     }
  798.     if (sts > 0)
  799.         crack_year (bp, &op->o_p.p_epoch);
  800.  
  801.     if (av) {
  802.         bp = av[6];
  803.         sts = 1;
  804.     } else {
  805.         static char p[] = "Absolute magnitude (g):";
  806.         f_prompt (p);
  807.         f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", op->o_p.p_g);
  808.         sts = read_line (buf, 8+1);
  809.         if (sts < 0)
  810.         return;
  811.         bp = buf;
  812.     }
  813.     if (sts > 0)
  814.         op->o_p.p_g = atof(bp);
  815.  
  816.     if (av) {
  817.         bp = av[7];
  818.         sts = 1;
  819.     } else {
  820.         static char p[] = "Luminosity index (kappa):";
  821.         f_prompt (p);
  822.         f_double (R_PROMPT, C_PROMPT+sizeof(p), "(%g) ", op->o_p.p_k);
  823.         sts = read_line (buf, 8+1);
  824.         if (sts < 0)
  825.         return;
  826.         bp = buf;
  827.     }
  828.     if (sts > 0)
  829.         op->o_p.p_k = atof(bp);
  830. }
  831.  
  832. /* given either a decimal year (xxxx. something) or a calendar (x/x/x)
  833.  * convert it to an mjd and store it at *p;
  834.  */
  835. static
  836. crack_year (bp, p)
  837. char *bp;
  838. double *p;
  839. {
  840.     if (decimal_year(bp)) {
  841.         double y = atof (bp);
  842.         year_mjd (y, p);
  843.     } else {
  844.         int m, y;
  845.         double d;
  846.         mjd_cal (*p, &m, &d, &y);    /* init with current */
  847.         f_sscandate (bp, &m, &d, &y);
  848.         cal_mjd (m, d, y, p);
  849.     }
  850. }
  851.