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