Changeset 2551 in Sophya for trunk/SophyaExt


Ignore:
Timestamp:
Jun 15, 2004, 6:54:12 PM (21 years ago)
Author:
cmv
Message:

nouvelle version de xephem/libastro (3.6) cmv 15/6/04

Location:
trunk/SophyaExt/XephemAstroLib
Files:
14 added
5 deleted
55 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaExt/XephemAstroLib/Makefile

    r2437 r2551  
    1010        if [ ! -d $(EXTINCPATH)/XAstro ] ; then mkdir $(EXTINCPATH)/XAstro ; fi
    1111        rm -f $(EXTINCPATH)/XAstro/astro.h
     12        rm -f $(EXTINCPATH)/XAstro/bdl.h
    1213        rm -f $(EXTINCPATH)/XAstro/P_.h
    1314        rm -f $(EXTLIBPATH)/libxastro.a

  • trunk/SophyaExt/XephemAstroLib/aa_hadec.c

    r1719 r2551  
    55#include <math.h>
    66
    7 #include "P_.h"
    87#include "astro.h"
    98
    10 static void aaha_aux P_((double lat, double x, double y, double *p, double *q));
     9static void aaha_aux (double lt, double x, double y, double *p, double *q);
    1110
    12 /* given geographical latitude (n+, radians), lat, altitude (up+, radians),
     11/* given geographical latitude (n+, radians), lt, altitude (up+, radians),
    1312 * alt, and azimuth (angle round to the east from north+, radians),
    1413 * return hour angle (radians), ha, and declination (radians), dec.
    1514 */
    1615void
    17 aa_hadec (lat, alt, az, ha, dec)
    18 double lat;
    19 double alt, az;
    20 double *ha, *dec;
     16aa_hadec (
     17double lt,
     18double alt, double az,
     19double *ha, double *dec)
    2120{
    22         aaha_aux (lat, az, alt, ha, dec);
     21        aaha_aux (lt, az, alt, ha, dec);
    2322        if (*ha > PI)
    2423            *ha -= 2*PI;
    2524}
    2625
    27 /* given geographical (n+, radians), lat, hour angle (radians), ha, and
     26/* given geographical (n+, radians), lt, hour angle (radians), ha, and
    2827 * declination (radians), dec, return altitude (up+, radians), alt, and
    2928 * azimuth (angle round to the east from north+, radians),
    3029 */
    3130void
    32 hadec_aa (lat, ha, dec, alt, az)
    33 double lat;
    34 double ha, dec;
    35 double *alt, *az;
     31hadec_aa (
     32double lt,
     33double ha, double dec,
     34double *alt, double *az)
    3635{
    37         aaha_aux (lat, ha, dec, az, alt);
     36        aaha_aux (lt, ha, dec, az, alt);
    3837}
    3938
    … …  
    4342 */
    4443double
    45 geoc_lat (phi)
    46 double phi;
     44geoc_lat (
     45double phi)
    4746{
    4847#define MAXLAT  degrad(89.9999) /* avoid tan() greater than this */
    … …  
    5655 */
    5756static void
    58 aaha_aux (lat, x, y, p, q)
    59 double lat;
    60 double x, y;
    61 double *p, *q;
     57aaha_aux (
     58double lt,
     59double x, double y,
     60double *p, double *q)
    6261{
    63         static double last_lat = -3434, slat, clat;
     62        static double last_lt = -3434, slt, clt;
    6463        double cap, B;
    6564
    66         if (lat != last_lat) {
    67             slat = sin(lat);
    68             clat = cos(lat);
    69             last_lat = lat;
     65        if (lt != last_lt) {
     66            slt = sin(lt);
     67            clt = cos(lt);
     68            last_lt = lt;
    7069        }
    7170
    72         solve_sphere (-x, PI/2-y, slat, clat, &cap, &B);
     71        solve_sphere (-x, PI/2-y, slt, clt, &cap, &B);
    7372        *p = B;
    7473        *q = PI/2 - acos(cap);
    … …  
    7675
    7776/* For RCS Only -- Do Not Edit */
    78 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: aa_hadec.c,v $ $Date: 2001-10-22 12:08:26 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     77static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: aa_hadec.c,v $ $Date: 2004-06-15 16:52:37 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/aberration.c

    r1719 r2551  
    22 * based on secular unperturbed Kepler orbit
    33 *
    4  * the corrections should be applied to ra/dec and lam/beta of the
     4 * the corrections should be applied to ra/dec and lam/beta at the
    55 * epoch of date.
    66 */
    … …  
    88#include <stdio.h>
    99#include <math.h>
     10#include <stdlib.h>
    1011
    11 #if defined(__STDC__)
    12 #include <stdlib.h>
    13 #endif
    14 
    15 #include "P_.h"
    1612#include "astro.h"
    1713
    … …  
    2016#define AB_EQ_EOD       1
    2117
    22 static void ab_aux P_((double mjd, double *x, double *y, double lsn,
    23         int mode));
     18static void ab_aux (double mj, double *x, double *y, double lsn, int mode);
    2419
    2520/* apply aberration correction to ecliptical coordinates *lam and *bet
    26  * (in radians) for a given time mjd and handily supplied longitude of sun,
     21 * (in radians) for a given time m and handily supplied longitude of sun,
    2722 * lsn (in radians)
    2823 */
    2924void
    30 ab_ecl (mjd, lsn, lam, bet)
    31 double mjd, lsn, *lam, *bet;
     25ab_ecl (double mj, double lsn, double *lam, double *bet)
    3226{
    33         ab_aux(mjd, lam, bet, lsn, AB_ECL_EOD);
     27        ab_aux(mj, lam, bet, lsn, AB_ECL_EOD);
    3428}
    3529
    3630/* apply aberration correction to equatoreal coordinates *ra and *dec
    37  * (in radians) for a given time mjd and handily supplied longitude of sun,
     31 * (in radians) for a given time m and handily supplied longitude of sun,
    3832 * lsn (in radians)
    3933 */
    4034void
    41 ab_eq (mjd, lsn, ra, dec)
    42 double mjd, lsn, *ra, *dec;
     35ab_eq (double mj, double lsn, double *ra, double *dec)
    4336{
    44         ab_aux(mjd, ra, dec, lsn, AB_EQ_EOD);
     37        ab_aux(mj, ra, dec, lsn, AB_EQ_EOD);
    4538}
    4639
    … …  
    5144 */
    5245static void
    53 ab_aux (mjd, x, y, lsn, mode)
    54 double mjd, *x, *y, lsn;
    55 int mode;
     46ab_aux (double mj, double *x, double *y, double lsn, int mode)
    5647{
    57         static double lastmjd = -10000;
     48        static double lastmj = -10000;
    5849        static double eexc;     /* earth orbit excentricity */
    5950        static double leperi;   /* ... and longitude of perihelion */
    6051        static char dirty = 1;  /* flag for cached trig terms */
    6152
    62         if (mjd != lastmjd) {
     53        if (mj != lastmj) {
    6354            double T;           /* centuries since J2000 */
    6455
    65             T = (mjd - J2000)/36525.;
     56            T = (mj - J2000)/36525.;
    6657            eexc = 0.016708617 - (42.037e-6 + 0.1236e-6 * T) * T;
    6758            leperi = degrad(102.93735 + (0.71953 + 0.00046 * T) * T);
    68             lastmjd = mjd;
     59            lastmj = mj;
    6960            dirty = 1;
    7061        }
    … …  
    9990                    cp = cos(leperi);
    10091                    sp = sin(leperi);
    101                     obliquity(mjd, &eps);
     92                    obliquity(mj, &eps);
    10293                    se = sin(eps);
    10394                    ce = cos(eps);
    … …  
    127118        default:
    128119            printf ("ab_aux: bad mode: %d\n", mode);
    129             exit (1);
     120            abort();
    130121            break;
    131122
    … …  
    134125
    135126/* For RCS Only -- Do Not Edit */
    136 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: aberration.c,v $ $Date: 2001-10-22 12:08:26 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     127static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: aberration.c,v $ $Date: 2004-06-15 16:52:37 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/actan.c

    r1719 r2551  
    11#include <math.h>
    22
    3 /* @(#) $Id: actan.c,v 1.2 2001-10-22 12:08:26 cmv Exp $ */
     3/* @(#) $Id: actan.c,v 1.3 2004-06-15 16:52:37 cmv Exp $ */
    44
    55/* commonly in math.h, but not in strict ANSI C */
    … …  
    6565
    6666/* For RCS Only -- Do Not Edit */
    67 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: actan.c,v $ $Date: 2001-10-22 12:08:26 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     67static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: actan.c,v $ $Date: 2004-06-15 16:52:37 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/airmass.c

    r1719 r2551  
    11#include <math.h>
    22
    3 #include "P_.h"
    43#include "astro.h"
    54
    … …  
    109 */
    1110void
    12 airmass (aa, Xp)
    13 double aa;              /* apparent altitude, rads */
    14 double *Xp;             /* airmasses */
     11airmass (
     12double aa,              /* apparent altitude, rads */
     13double *Xp)             /* airmasses */
    1514{
    1615        double sm1;     /* secant zenith angle, minus 1 */
    … …  
    2524
    2625/* For RCS Only -- Do Not Edit */
    27 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: airmass.c,v $ $Date: 2001-10-22 12:08:26 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     26static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: airmass.c,v $ $Date: 2004-06-15 16:52:37 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/anomaly.c

    r1719 r2551  
    44#include <math.h>
    55
    6 #include "P_.h"
    76#include "astro.h"
    87
    … …  
    1615 */
    1716void
    18 anomaly (ma, s, nu, ea)
    19 double ma, s;
    20 double *nu, *ea;
     17anomaly (double ma, double s, double *nu, double *ea)
    2118{
    2219        double m, fea, corr;
    … …  
    6461
    6562/* For RCS Only -- Do Not Edit */
    66 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: anomaly.c,v $ $Date: 2001-10-22 12:08:26 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     63static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: anomaly.c,v $ $Date: 2004-06-15 16:52:37 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/ap_as.c

    r1719 r2551  
    11#include <string.h>
    22
    3 #include "P_.h"
    43#include "astro.h"
    5 #include "circum.h"
    64
    75/* convert the given apparent RA/Dec to astrometric precessed to Mjd IN PLACE.
    … …  
    119 */
    1210void
    13 ap_as (np, Mjd, rap, decp)
    14 Now *np;
    15 double Mjd;
    16 double *rap, *decp;
     11ap_as (Now *np, double Mjd, double *rap, double *decp)
    1712{
    1813        Obj o;
    … …  
    2924        *rap -= o.s_ra - *rap;
    3025        *decp -= o.s_dec - *decp;
    31         if (*decp >  PI/2) {*decp =  PI - *decp; *rap += PI; }
    32         if (*decp < -PI/2) {*decp = -PI - *decp; *rap += PI; }
    33         range (rap, 2*PI);
     26        radecrange (rap, decp);
    3427        precess (mjd, Mjd, rap, decp);
     28        radecrange (rap, decp);
    3529}
    3630
    … …  
    3933 */
    4034void
    41 as_ap (np, Mjd, rap, decp)
    42 Now *np;
    43 double Mjd;
    44 double *rap, *decp;
     35as_ap (Now *np, double Mjd, double *rap, double *decp)
    4536{
    4637        Obj o;
    … …  
    6051
    6152/* For RCS Only -- Do Not Edit */
    62 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: ap_as.c,v $ $Date: 2001-10-22 12:08:26 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     53static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: ap_as.c,v $ $Date: 2004-06-15 16:52:37 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/astro.h

    r1719 r2551  
     1#ifndef _ASTRO_H
     2#define _ASTRO_H
     3
    14#ifndef PI
    25#define PI              3.141592653589793
    … …  
    2023 * N.B. only the first 8 are valid for use with plans().
    2124 */
    22 #define MERCURY 0
    23 #define VENUS   1
    24 #define MARS    2
    25 #define JUPITER 3
    26 #define SATURN  4
    27 #define URANUS  5
    28 #define NEPTUNE 6
    29 #define PLUTO   7
    30 
    31 /* a few more handy ones */
    32 #define SUN     8
    33 #define MOON    9
    34 #define OBJX    10
    35 #define OBJY    11
    36 #define OBJZ    12
    37 #define NOBJ    13      /* total number of basic objects */
     25typedef enum {
     26    MERCURY,
     27    VENUS,
     28    MARS,
     29    JUPITER,
     30    SATURN,
     31    URANUS,
     32    NEPTUNE,
     33    PLUTO,
     34    SUN,
     35    MOON,
     36    NOBJ        /* total number of basic objects */
     37} PLCode;
     38
     39/* moon constants for pl_moon */
     40typedef enum {
     41    X_PLANET = 0,                       /* use to mean planet itself */
     42    PHOBOS = NOBJ, DEIMOS,
     43    IO, EUROPA, GANYMEDE, CALLISTO,
     44    MIMAS, ENCELADUS, TETHYS, DIONE, RHEA, TITAN, HYPERION, IAPETUS,
     45    ARIEL, UMBRIEL, TITANIA, OBERON, MIRANDA,
     46    NBUILTIN
     47} MCode;
    3848
    3949/* starting point for MJD calculations
    4050 */
    4151#define MJD0  2415020.0
    42 #define J2000 (2451545.0 - MJD0)      /* let compiler optimise */
     52#define J2000 (2451545.0 - MJD0)      /* yes, 2000 January 1 at 12h */
     53
     54/* the Now and Obj typedefs.
     55 * also, a few miscellaneous constants and declarations.
     56 */
     57
     58#define SPD     (24.0*3600.0)   /* seconds per day */
     59#define MAU     (1.4959787e11)  /* m / au */
     60#define LTAU    499.005         /* seconds light takes to travel 1 AU */
     61#define ERAD    (6.37816e6)     /* earth equitorial radius, m */
     62#define MRAD    (1.740e6)       /* moon equitorial radius, m */
     63#define SRAD    (6.95e8)        /* sun equitorial radius, m */
     64#define FTPM    3.28084         /* ft per m */
     65#define ESAT_MAG        2       /* default satellite magnitude */
     66#define FAST_SAT_RPD    0.25    /* max earth sat rev/day considered "fast" */
     67
     68#define EOD     (-9786)         /* special epoch flag: use epoch of date */
     69
     70/* info about the local observing circumstances and misc preferences */
     71typedef struct {
     72        double n_mjd;   /* modified Julian date, ie, days since
     73                         * Jan 0.5 1900 (== 12 noon, Dec 30, 1899), utc.
     74                         * enough precision to get well better than 1 second.
     75                         * N.B. if not first member, must move NOMJD inits.
     76                         */
     77        double n_lat;   /* geographic (surface-normal) lt, >0 north, rads */
     78        double n_lng;   /* longitude, >0 east, rads */
     79        double n_tz;    /* time zone, hrs behind UTC */
     80        double n_temp;  /* atmospheric temp, degrees C */
     81        double n_pressure; /* atmospheric pressure, mBar */
     82        double n_elev;  /* elevation above sea level, earth radii */
     83        double n_dip;   /* dip of sun below hzn at twilight, >0 below, rads */
     84        double n_epoch; /* desired precession display ep as an mjd, or EOD */
     85        char n_tznm[8]; /* time zone name; 7 chars or less, always 0 at end */
     86} Now;
     87
     88/* handy shorthands for fields in a Now pointer, np */
     89#define mjd     np->n_mjd
     90#define lat     np->n_lat
     91#define lng     np->n_lng
     92#define tz      np->n_tz
     93#define temp    np->n_temp
     94#define pressure np->n_pressure
     95#define elev    np->n_elev
     96#define dip     np->n_dip
     97#define epoch   np->n_epoch
     98#define tznm    np->n_tznm
     99#define mjed    mm_mjed(np)
     100
     101/* structures to describe objects of various types.
     102 */
     103
     104/* magnitude values in two different systems */
     105typedef struct {
     106    float m1, m2;       /* either g/k or H/G, depending on... */
     107    int whichm;         /* one of MAG_gk or MAG_HG */
     108} Mag;
     109
     110/* whichm */
     111#define MAG_HG          0       /* using 0 makes HG the initial default */
     112#define MAG_gk          1
     113
     114/* we actually store magnitudes times this scale factor in a short int */
     115#define MAGSCALE        100.0
     116#define set_smag(op,m)  ((op)->s_mag = (short)floor((m)*MAGSCALE + 0.5))
     117#define set_fmag(op,m)  ((op)->f_mag = (short)floor((m)*MAGSCALE + 0.5))
     118#define get_mag(op)     ((op)->s_mag / MAGSCALE)
     119
     120/* longest object name, including trailing '\0' */
     121#define MAXNM   21
     122
     123typedef unsigned char ObjType_t;
     124typedef unsigned char ObjAge_t;
     125typedef unsigned char byte;
     126
     127/* Obj is a massive union.
     128 * many fields are in common so we use macros to make things a little easier.
     129 */
     130
     131/* fields common to *all* structs in the Obj union */
     132#define OBJ_COMMON_FLDS                                                 \
     133    ObjType_t co_type;  /* current object type; see flags, below */     \
     134    byte co_flags;      /* FUSER*... used by others */                  \
     135    ObjAge_t co_age;    /* update aging code; see db.c */               \
     136    char co_name[MAXNM];/* name, including \0 */                        \
     137    float co_ra;        /* geo/topo app/mean ra, rads */                \
     138    float co_dec;       /* geo/topo app/mean dec, rads */               \
     139    float co_gaera;     /* geo apparent ra, rads */                     \
     140    float co_gaedec;    /* geo apparent dec, rads */                    \
     141    float co_az;        /* azimuth, >0 e of n, rads */                  \
     142    float co_alt;       /* altitude above topocentric horizon, rads */  \
     143    float co_elong;     /* angular sep btwen obj and sun, >0 E, degs */ \
     144    float co_size;      /* angular size, arc secs */                    \
     145    short co_mag        /* visual magnitude * MAGSCALE */
     146
     147/* fields common to all solar system objects in the Obj union */
     148#define OBJ_SOLSYS_FLDS                                                 \
     149    OBJ_COMMON_FLDS;    /* all the fixed ones plus ... */               \
     150    float so_sdist;     /* dist from object to sun, au */               \
     151    float so_edist;     /* dist from object to earth, au */             \
     152    float so_hlong;     /* heliocentric longitude, rads */              \
     153    float so_hlat;      /* heliocentric latitude, rads */               \
     154    float so_phase      /* phase, % */
     155
     156/* fields common to all fixed objects in the Obj union */
     157#define OBJ_FIXED_FLDS                                                  \
     158    char  fo_spect[2];  /* spectral codes, if appropriate */            \
     159    float fo_epoch;     /* eq of ra/dec and time when pm=0; mjd */      \
     160    float fo_ra;        /* ra, rads, in epoch frame */                  \
     161    float fo_dec;       /* dec, rads, in epoch frame */                 \
     162    float fo_pmra;      /* ra proper motion, rads/day/cos(dec) */       \
     163    float fo_pmdec;     /* dec proper motion, rads/day */               \
     164    char  fo_class      /* object class */
     165
     166/* a generic object */
     167typedef struct {
     168    OBJ_COMMON_FLDS;
     169} ObjAny;
     170
     171/* a generic sol system object */
     172typedef struct {
     173    OBJ_SOLSYS_FLDS;
     174} ObjSS;
     175
     176/* basic Fixed object info.
     177 */
     178typedef struct {
     179    OBJ_COMMON_FLDS;
     180    OBJ_FIXED_FLDS;
     181
     182    /* following are for galaxies */
     183    byte  fo_ratio;     /* minor/major diameter ratio. use s/get_ratio() */
     184    byte  fo_pa;        /* position angle, E of N, rads. use s/get_pa() */
     185} ObjF;
     186
     187/* true-orbit parameters of binary-star object type */
     188typedef struct {
     189    float bo_T;         /* epoch of periastron, years */
     190    float bo_e;         /* eccentricity */
     191    float bo_o;         /* argument of periastron, degress */
     192    float bo_O;         /* longitude of node, degrees */
     193    float bo_i;         /* inclination to plane of sky, degrees */
     194    float bo_a;         /* semi major axis, arc secs */
     195    float bo_P;         /* period, years */
     196
     197    /* companion position, computed by obj_cir() iff b_2compute */
     198    float bo_pa;        /* position angle @ ep, rads E of N */
     199    float bo_sep;       /* separation @ ep, arc secs */
     200    float bo_ra;        /* geo/topo app/mean ra, rads */
     201    float bo_dec;       /* geo/topo app/mean dec, rads */
     202} BinOrbit;
     203typedef struct {
     204    float bp_ep;        /* epoch of pa/sep, year */
     205    float bp_pa;        /* position angle @ ep, rads E of N */
     206    float bp_sep;       /* separation @ ep, arc secs */
     207
     208    /* companion position, computed by obj_cir() iff b_2compute */
     209    float bp_ra;        /* geo/topo app/mean ra, rads */
     210    float bp_dec;       /* geo/topo app/mean dec, rads */
     211} BinPos;
     212#define MAXBINPOS 2     /* max discrete epochs to store when no elements */
     213typedef struct {
     214    OBJ_COMMON_FLDS;
     215    OBJ_FIXED_FLDS;
     216
     217    byte b_2compute;    /* whether to compute secondary positions */
     218    byte b_nbp;         /* number of b_bp[] or 0 to use b_bo */
     219    short b_2mag;       /* secondary's magnitude * MAGSCALE */
     220    char b_2spect[2];   /* secondary's spectrum */
     221
     222    /* either a real orbit or a set of discrete pa/sep */
     223    union {
     224        BinOrbit b_bo;                  /* orbital elements */
     225        BinPos b_bp[MAXBINPOS];         /* table of discrete positions */
     226    } u;
     227} ObjB;
     228
     229#define fo_mag  co_mag  /* pseudonym for so_mag since it is not computed */
     230#define fo_size co_size /* pseudonym for so_size since it is not computed */
     231
     232/* macros to pack/unpack some fields */
     233#define SRSCALE         255.0           /* galaxy size ratio scale */
     234#define PASCALE         (255.0/(2*PI))  /* pos angle scale factor */
     235#define get_ratio(op)   (((int)(op)->f_ratio)/SRSCALE)
     236#define set_ratio(op,maj,min) ((op)->f_ratio = (byte)(((maj) > 0)           \
     237                                        ? ((min)*SRSCALE/(double)(maj)+0.5) \
     238                                        : 0))
     239#define get_pa(op)      ((double)(op)->f_pa/PASCALE)
     240#define set_pa(op,s)    ((op)->f_pa = (byte)((s)*PASCALE + 0.5))
     241
     242#define NCLASSES        128 /* n potential fo_classes -- allow for all ASCII */
     243
     244/* basic planet object info */
     245typedef struct {
     246    OBJ_SOLSYS_FLDS;
     247    PLCode plo_code;            /* which planet */
     248    MCode plo_moon;             /* which moon, or X_PLANET if planet */
     249    char plo_evis, plo_svis;    /* if moon: whether visible from earth, sun */
     250    double plo_x, plo_y, plo_z; /* if moon: eq dist from center, planet radii */
     251    double plo_aux1, plo_aux2;  /* various values, depending on type */
     252} ObjPl;
     253
     254/* basic info about an object in elliptical heliocentric orbit */
     255typedef struct {
     256    OBJ_SOLSYS_FLDS;
     257    float  eo_inc;      /* inclination, degrees */
     258    float  eo_Om;       /* longitude of ascending node, degrees */
     259    float  eo_om;       /* argument of perihelion, degress */
     260    float  eo_a;        /* mean distance, aka,semi-maj axis,AU */
     261    float  eo_M;        /* mean anomaly, ie, degrees from perihelion at cepoch*/
     262    float  eo_size;     /* angular size, in arc seconds at 1 AU */
     263    float  eo_startok;  /* nominal first mjd this set is ok, else 0 */
     264    float  eo_endok;    /* nominal last mjd this set is ok, else 0 */
     265    double eo_e;        /* eccentricity (double for when near 1 computing q) */
     266    double eo_cepoch;   /* epoch date (M reference), as an mjd */
     267    double eo_epoch;    /* equinox year (inc/Om/om reference), as an mjd. */
     268    Mag    eo_mag;      /* magnitude */
     269} ObjE;
     270
     271/* basic info about an object in hyperbolic heliocentric orbit */
     272typedef struct {
     273    OBJ_SOLSYS_FLDS;
     274    double ho_epoch;    /* equinox year (inc/Om/om reference), as an mjd */
     275    double ho_ep;       /* epoch of perihelion, as an mjd */
     276    float  ho_startok;  /* nominal first mjd this set is ok, else 0 */
     277    float  ho_endok;    /* nominal last mjd this set is ok, else 0 */
     278    float  ho_inc;      /* inclination, degs */
     279    float  ho_Om;       /* longitude of ascending node, degs */
     280    float  ho_om;       /* argument of perihelion, degs. */
     281    float  ho_e;        /* eccentricity */
     282    float  ho_qp;       /* perihelion distance, AU */
     283    float  ho_g, ho_k;  /* magnitude model coefficients */
     284    float  ho_size;     /* angular size, in arc seconds at 1 AU */
     285} ObjH;
     286
     287/* basic info about an object in parabolic heliocentric orbit */
     288typedef struct {
     289    OBJ_SOLSYS_FLDS;
     290    double po_epoch;    /* reference epoch, as an mjd */
     291    double po_ep;       /* epoch of perihelion, as an mjd */
     292    float  po_startok;  /* nominal first mjd this set is ok, else 0 */
     293    float  po_endok;    /* nominal last mjd this set is ok, else 0 */
     294    float  po_inc;      /* inclination, degs */
     295    float  po_qp;       /* perihelion distance, AU */
     296    float  po_om;       /* argument of perihelion, degs. */
     297    float  po_Om;       /* longitude of ascending node, degs */
     298    float  po_g, po_k;  /* magnitude model coefficients */
     299    float  po_size;     /* angular size, in arc seconds at 1 AU */
     300} ObjP;
     301
     302/* basic earth satellite object info */
     303typedef struct {
     304    OBJ_COMMON_FLDS;
     305    double eso_epoch;   /* reference epoch, as an mjd */
     306    double eso_n;       /* mean motion, rev/day
     307                         * N.B. we need double due to a sensitive differencing
     308                         * operation used to compute MeanAnomaly in
     309                         * esat_main()/satellite.c.
     310                         */
     311    float  eso_startok; /* nominal first mjd this set is ok, else 0 */
     312    float  eso_endok;   /* nominal last mjd this set is ok, else 0 */
     313    float  eso_inc;     /* inclination, degs */
     314    float  eso_raan;    /* RA of ascending node, degs */
     315    float  eso_e;       /* eccentricity */
     316    float  eso_ap;      /* argument of perigee at epoch, degs */
     317    float  eso_M;       /* mean anomaly, ie, degrees from perigee at epoch */
     318    float  eso_decay;   /* orbit decay rate, rev/day^2 */
     319    float  eso_drag;    /* object drag coefficient, (earth radii)^-1 */
     320    int    eso_orbit;   /* integer orbit number of epoch */
     321
     322    /* computed "sky" results unique to earth satellites */
     323    float  ess_elev;    /* height of satellite above sea level, m */
     324    float  ess_range;   /* line-of-site distance from observer to satellite, m*/
     325    float  ess_rangev;  /* rate-of-change of range, m/s */
     326    float  ess_sublat;  /* latitude below satellite, >0 north, rads */
     327    float  ess_sublng;  /* longitude below satellite, >0 east, rads */
     328    int    ess_eclipsed;/* 1 if satellite is in earth's shadow, else 0 */
     329} ObjES;
     330
     331typedef union {
     332    ObjAny  any;        /* these fields valid for all types */
     333    ObjSS   anyss;      /* these fields valid for all solar system types */
     334    ObjPl   pl;         /* planet */
     335    ObjF    f;          /* fixed object, plus proper motion */
     336    ObjB    b;          /* bona fide binary stars (doubles are stored in f) */
     337    ObjE    e;          /* object in heliocentric elliptical orbit */
     338    ObjH    h;          /* object in heliocentric hyperbolic trajectory */
     339    ObjP    p;          /* object in heliocentric parabolic trajectory */
     340    ObjES   es;         /* earth satellite */
     341} Obj;
     342
     343
     344/* for o_flags -- everybody must agree */
     345#define FUSER0          0x01
     346#define FUSER1          0x02
     347#define FUSER2          0x04
     348#define FUSER3          0x08
     349#define FUSER4          0x10
     350#define FUSER5          0x20
     351#define FUSER6          0x40
     352#define FUSER7          0x80
     353
     354/* mark an object as being a "field star" */
     355#define FLDSTAR         FUSER3
     356/* mark an object as circum calculation failed */
     357#define NOCIRCUM        FUSER7
     358
     359/* Obj shorthands: */
     360#define o_type          any.co_type
     361#define o_name          any.co_name
     362#define o_flags         any.co_flags
     363#define o_age           any.co_age
     364#define s_ra            any.co_ra
     365#define s_dec           any.co_dec
     366#define s_gaera         any.co_gaera
     367#define s_gaedec        any.co_gaedec
     368#define s_az            any.co_az
     369#define s_alt           any.co_alt
     370#define s_elong         any.co_elong
     371#define s_size          any.co_size
     372#define s_mag           any.co_mag
     373
     374#define s_sdist         anyss.so_sdist
     375#define s_edist         anyss.so_edist
     376#define s_hlong         anyss.so_hlong
     377#define s_hlat          anyss.so_hlat
     378#define s_phase         anyss.so_phase
     379
     380#define s_elev          es.ess_elev
     381#define s_range         es.ess_range
     382#define s_rangev        es.ess_rangev
     383#define s_sublat        es.ess_sublat
     384#define s_sublng        es.ess_sublng
     385#define s_eclipsed      es.ess_eclipsed
     386
     387#define f_class         f.fo_class
     388#define f_spect         f.fo_spect
     389#define f_ratio         f.fo_ratio
     390#define f_pa            f.fo_pa
     391#define f_epoch         f.fo_epoch
     392#define f_RA            f.fo_ra
     393#define f_pmRA          f.fo_pmra
     394#define f_dec           f.fo_dec
     395#define f_pmdec         f.fo_pmdec
     396#define f_mag           f.fo_mag
     397#define f_size          f.fo_size
     398
     399#define e_cepoch        e.eo_cepoch
     400#define e_epoch         e.eo_epoch
     401#define e_startok       e.eo_startok
     402#define e_endok         e.eo_endok
     403#define e_inc           e.eo_inc
     404#define e_Om            e.eo_Om
     405#define e_om            e.eo_om
     406#define e_a             e.eo_a
     407#define e_e             e.eo_e
     408#define e_M             e.eo_M
     409#define e_size          e.eo_size
     410#define e_mag           e.eo_mag
     411
     412#define h_epoch         h.ho_epoch
     413#define h_startok       h.ho_startok
     414#define h_endok         h.ho_endok
     415#define h_ep            h.ho_ep
     416#define h_inc           h.ho_inc
     417#define h_Om            h.ho_Om
     418#define h_om            h.ho_om
     419#define h_e             h.ho_e
     420#define h_qp            h.ho_qp
     421#define h_g             h.ho_g
     422#define h_k             h.ho_k
     423#define h_size          h.ho_size
     424
     425#define p_epoch         p.po_epoch
     426#define p_startok       p.po_startok
     427#define p_endok         p.po_endok
     428#define p_ep            p.po_ep
     429#define p_inc           p.po_inc
     430#define p_qp            p.po_qp
     431#define p_om            p.po_om
     432#define p_Om            p.po_Om
     433#define p_g             p.po_g
     434#define p_k             p.po_k
     435#define p_size          p.po_size
     436
     437#define es_epoch        es.eso_epoch
     438#define es_startok      es.eso_startok
     439#define es_endok        es.eso_endok
     440#define es_inc          es.eso_inc
     441#define es_raan         es.eso_raan
     442#define es_e            es.eso_e
     443#define es_ap           es.eso_ap
     444#define es_M            es.eso_M
     445#define es_n            es.eso_n
     446#define es_decay        es.eso_decay
     447#define es_drag         es.eso_drag
     448#define es_orbit        es.eso_orbit
     449
     450#define pl_code         pl.plo_code
     451#define pl_moon         pl.plo_moon
     452#define pl_evis         pl.plo_evis
     453#define pl_svis         pl.plo_svis
     454#define pl_x            pl.plo_x
     455#define pl_y            pl.plo_y
     456#define pl_z            pl.plo_z
     457#define pl_aux1         pl.plo_aux1
     458#define pl_aux2         pl.plo_aux2
     459
     460#define b_2compute      b.b_2compute
     461#define b_2spect        b.b_2spect
     462#define b_2mag          b.b_2mag
     463#define b_bo            b.u.b_bo
     464#define b_bp            b.u.b_bp
     465#define b_nbp           b.b_nbp
     466
     467/* insure we always refer to the fields and no monkey business */
     468#undef OBJ_COMMON_FLDS
     469#undef OBJ_SOLSYS_FLDS
     470
     471/* o_type code.
     472 * N.B. names are assigned in order in objmenu.c
     473 * N.B. if add one add switch in obj_cir().
     474 * N.B. UNDEFOBJ must be zero so new objects are undefinied by being zeroed.
     475 * N.B. maintain the bitmasks too.
     476 */
     477enum ObjType {
     478    UNDEFOBJ=0,
     479    FIXED, BINARYSTAR, ELLIPTICAL, HYPERBOLIC, PARABOLIC, EARTHSAT, PLANET,
     480    NOBJTYPES
     481};
     482
     483/* types as handy bitmasks too */
     484#define OBJTYPE2MASK(t) (1<<(t))
     485#define FIXEDM          OBJTYPE2MASK(FIXED)
     486#define BINARYSTARM     OBJTYPE2MASK(BINARYSTAR)
     487#define ELLIPTICALM     OBJTYPE2MASK(ELLIPTICAL)
     488#define HYPERBOLICM     OBJTYPE2MASK(HYPERBOLIC)
     489#define PARABOLICM      OBJTYPE2MASK(PARABOLIC)
     490#define EARTHSATM       OBJTYPE2MASK(EARTHSAT)
     491#define PLANETM         OBJTYPE2MASK(PLANET)
     492#define ALLM            (~0)
     493
     494/* rise, set and transit information.
     495 */
     496typedef struct {
     497    int rs_flags;       /* info about what has been computed and any
     498                         * special conditions; see flags, below.
     499                         */
     500    double rs_risetm;   /* mjd time of rise today */
     501    double rs_riseaz;   /* azimuth of rise, rads E of N */
     502    double rs_trantm;   /* mjd time of transit today */
     503    double rs_tranalt;  /* altitude of transit, rads up from horizon */
     504    double rs_settm;    /* mjd time of set today */
     505    double rs_setaz;    /* azimuth of set, rads E of N */
     506} RiseSet;
     507
     508/* RiseSet flags */
     509#define RS_NORISE       0x0001  /* object does not rise as such today */
     510#define RS_NOSET        0x0002  /* object does not set as such today */
     511#define RS_NOTRANS      0x0004  /* object does not transit as such today */
     512#define RS_CIRCUMPOLAR  0x0010  /* object stays up all day today */
     513#define RS_NEVERUP      0x0020  /* object never up at all today */
     514#define RS_ERROR        0x1000  /* can't figure out anything! */
     515#define RS_RISERR       (0x0100|RS_ERROR) /* error computing rise */
     516#define RS_SETERR       (0x0200|RS_ERROR) /* error computing set */
     517#define RS_TRANSERR     (0x0400|RS_ERROR) /* error computing transit */
     518
     519#define is_type(op,m)   (OBJTYPE2MASK((op)->o_type) & (m))
     520
     521/* any planet or its moons */
     522#define is_planet(op,p) (is_type(op,PLANETM) && op->pl_code == (p))
     523
     524/* any solar system object */
     525#define is_ssobj(op)    is_type(op,PLANETM|HYPERBOLICM|PARABOLICM|ELLIPTICALM)
     526
     527
     528
     529/* natural satellite support */
     530
     531typedef struct {
     532    char *full;         /* full name */
     533    char *tag;          /* Roman numeral tag */
     534    float x, y, z;      /* radii: +x:east +y:south +z:front */
     535    float ra, dec;
     536    float mag;
     537    int evis;           /* whether geometrically visible from earth */
     538    int svis;           /* whether in sun light */
     539} MoonData;
     540
     541/* separate set for each planet -- use in pl_moon */
     542
     543
     544enum _marsmoons {
     545    M_MARS = 0,                                 /* == X_PLANET */
     546    M_PHOBOS, M_DEIMOS,
     547    M_NMOONS                                    /* including planet at 0 */
     548};
     549
     550enum _jupmoons {
     551    J_JUPITER = 0,                              /* == X_PLANET */
     552    J_IO, J_EUROPA, J_GANYMEDE, J_CALLISTO,
     553    J_NMOONS                                    /* including planet */
     554};
     555
     556enum _satmoons {
     557    S_SATURN = 0,                               /* == X_PLANET */
     558    S_MIMAS, S_ENCELADUS, S_TETHYS, S_DIONE,
     559    S_RHEA, S_TITAN, S_HYPERION, S_IAPETUS,
     560    S_NMOONS                                    /* including planet */
     561};
     562
     563enum _uramoons {
     564    U_URANUS = 0,                               /* == X_PLANET */
     565    U_ARIEL, U_UMBRIEL, U_TITANIA, U_OBERON, U_MIRANDA,
     566    U_NMOONS                                    /* including planet */
     567};
     568
     569#define X_MAXNMOONS     S_NMOONS                /* N.B. chosen by hand */
    43570
    44571
    … …  
    46573/* global function declarations */
    47574
     575
    48576/* aa_hadec.c */
    49 extern void aa_hadec P_((double lat, double alt, double az, double *ha,
    50     double *dec));
    51 extern void hadec_aa P_((double lat, double ha, double dec, double *alt,
    52     double *az));
     577extern void aa_hadec (double lt, double alt, double az, double *ha,
     578    double *dec);
     579extern void hadec_aa (double lt, double ha, double dec, double *alt,
     580    double *az);
    53581
    54582/* aberration.c */
    55 extern void ab_ecl P_((double mjd, double lsn, double *lam, double *bet));
    56 extern void ab_eq P_((double mjd, double lsn, double *ra, double *dec));
     583extern void ab_ecl (double m, double lsn, double *lam, double *bet);
     584extern void ab_eq (double m, double lsn, double *ra, double *dec);
    57585
    58586/* airmass.c */
    59 extern void airmass P_((double aa, double *Xp));
     587extern void airmass (double aa, double *Xp);
    60588
    61589/* anomaly.c */
    62 extern void anomaly P_((double ma, double s, double *nu, double *ea));
     590extern void anomaly (double ma, double s, double *nu, double *ea);
     591
     592/* ap_as.c */
     593extern void ap_as ( Now *np, double Mjd, double *rap, double *decp);
     594extern void as_ap ( Now *np, double Mjd, double *rap, double *decp);
     595
     596/* atlas.c */
     597extern char *um_atlas (double ra, double dec);
     598extern char *u2k_atlas (double ra, double dec);
     599extern char *msa_atlas (double ra, double dec);
     600
     601/* aux.c */
     602extern double mm_mjed (Now *np);
    63603
    64604/* chap95.c */
    65 extern int chap95 P_((double mjd, int obj, double prec, double *ret));
     605extern int chap95 (double m, int obj, double prec, double *ret);
    66606
    67607/* chap95_data.c */
    68608
     609/* circum.c */
     610extern int obj_cir (Now *np, Obj *op);
     611
    69612/* comet.c */
    70 extern void comet P_((double mjd, double ep, double inc, double ap, double qp,
     613extern void comet (double m, double ep, double inc, double ap, double qp,
    71614    double om, double *lpd, double *psi, double *rp, double *rho, double *lam,
    72     double *bet));
     615    double *bet);
     616
     617/* constel.c */
     618#define NCNS    89
     619extern int cns_pick (double r, double d, double e);
     620extern int cns_id (char *abbrev);
     621extern char *cns_name (int id);
     622extern int cns_edges (double e, double **ra0p, double **dec0p, double **ra1p,
     623    double **dec1p);
     624extern int cns_list (double ra, double dec, double e, double rad, int ids[]);
     625extern int cns_figure (int id, double e, double ra[],double dec[],int dcodes[]);
     626extern int cns_reyfigure (int id, double e, double ra[], double dec[],
     627    int dcodes[]);
     628
     629/* dbfmt.c */
     630extern int db_crack_line (char s[], Obj *op, char nm[][MAXNM], int nnm,
     631    char whynot[]);
     632extern void db_write_line (Obj *op, char *lp);
     633extern int dbline_candidate (char line[]);
     634extern int get_fields (char *s, int delim, char *fields[]);
     635extern int db_tle (char *name, char *l1, char *l2, Obj *op);
     636extern int dateRangeOK (Now *np, Obj *op);
    73637
    74638/* deltat.c */
    75 extern double deltat P_((double mjd));
     639extern double deltat (double m);
     640
     641/* earthsat.c */
     642extern int obj_earthsat (Now *np, Obj *op);
    76643
    77644/* eq_ecl.c */
    78 extern void eq_ecl P_((double mjd, double ra, double dec, double *lat,
    79     double *lng));
    80 extern void ecl_eq P_((double mjd, double lat, double lng, double *ra,
    81     double *dec));
     645extern void eq_ecl (double m, double ra, double dec, double *lt,double *lg);
     646extern void ecl_eq (double m, double lt, double lg, double *ra,double *dec);
    82647
    83648/* eq_gal.c */
    84 extern void eq_gal P_((double mjd, double ra, double dec, double *lat,
    85     double *lng));
    86 extern void gal_eq P_((double mjd, double lat, double lng, double *ra,
    87     double *dec));
     649extern void eq_gal (double m, double ra, double dec, double *lt,double *lg);
     650extern void gal_eq (double m, double lt, double lg, double *ra,double *dec);
    88651
    89652/* formats.c */
    90 extern void fs_sexa P_((char *out, double a, int w, int fracbase));
    91 extern void fs_date P_((char out[], double jd));
    92 extern void f_scansex P_((double o, char *str, double *dp));
    93 extern void f_sscandate P_((char *bp, int pref, int *m, double *d, int *y));
    94 extern int scansex P_((char *str, double *dp));
     653extern int fs_sexa (char *out, double a, int w, int fracbase);
     654extern int fs_date (char out[], double jd);
     655extern int f_scansexa (const char *str, double *dp);
     656extern void f_sscandate (char *bp, int pref, int *m, double *d, int *y);
    95657
    96658/* helio.c */
    97 extern void heliocorr P_((double jd, double ra, double dec, double *hcp));
     659extern void heliocorr (double jd, double ra, double dec, double *hcp);
     660
     661/* jupmoon.c */
     662extern void jupiter_data (double Mjd, char dir[], Obj *eop, Obj *jop,
     663    double *jupsize, double *cmlI, double *cmlII, MoonData md[J_NMOONS]);
    98664
    99665/* libration.c */
    100 extern void llibration P_((double JD, double *llatp, double *llonp));
     666extern void llibration (double JD, double *llatp, double *llonp);
     667
     668/* magdecl.c */
     669extern int magdecl (double l, double L, double e, double y, char *dir,
     670    double *dp, char *err);
     671
     672/* marsmoon.c */
     673extern void marsm_data (double Mjd, char dir[], Obj *eop, Obj *mop,
     674    double *marssize, MoonData md[M_NMOONS]);
    101675
    102676/* misc.c */
    103 extern void zero_mem P_((void *loc, unsigned len));
    104 extern int tickmarks P_((double min, double max, int numdiv, double ticks[]));
    105 extern int lc P_((int cx, int cy, int cw, int x1, int y1, int x2, int y2,
    106     int *sx1, int *sy1, int *sx2, int *sy2));
    107 extern void hg_mag P_((double h, double g, double rp, double rho, double rsn,
    108     double *mp));
    109 extern int magdiam P_((int fmag, int magstp, double scale, double mag,
    110     double size));
    111 extern void gk_mag P_((double g, double k, double rp, double rho, double *mp));
    112 extern double atod P_((char *buf));
    113 extern void solve_sphere P_((double A, double b, double cc, double sc,
    114     double *cap, double *Bp));
    115 extern double delra P_((double dra));
     677extern void zero_mem (void *loc, unsigned len);
     678extern int tickmarks (double min, double max, int numdiv, double ticks[]);
     679extern int lc (int cx, int cy, int cw, int x1, int y1, int x2, int y2,
     680    int *sx1, int *sy1, int *sx2, int *sy2);
     681extern void hg_mag (double h, double g, double rp, double rho, double rsn,
     682    double *mp);
     683extern int magdiam (int fmag, int magstp, double scale, double mag,
     684    double size);
     685extern void gk_mag (double g, double k, double rp, double rho, double *mp);
     686extern double atod (char *buf);
     687extern void solve_sphere (double A, double b, double cc, double sc,
     688    double *cap, double *Bp);
     689extern double delra (double dra);
     690extern void now_lst (Now *np, double *lstp);
     691extern void radec2ha (Now *np, double ra, double dec, double *hap);
     692extern char *obj_description (Obj *op);
     693extern int is_deepsky (Obj *op);
    116694
    117695/* mjd.c */
    118 extern void cal_mjd P_((int mn, double dy, int yr, double *mjd));
    119 extern void mjd_cal P_((double mjd, int *mn, double *dy, int *yr));
    120 extern int mjd_dow P_((double mjd, int *dow));
    121 extern int isleapyear P_((int year));
    122 extern void mjd_dpm P_((double mjd, int *ndays));
    123 extern void mjd_year P_((double mjd, double *yr));
    124 extern void year_mjd P_((double y, double *mjd));
    125 extern void rnd_second P_((double *t));
    126 extern void mjd_dayno P_((double jd, int *yr, double *dy));
    127 extern double mjd_day P_((double jd));
    128 extern double mjd_hr P_((double jd));
    129 extern void range P_((double *v, double r));
     696extern void cal_mjd (int mn, double dy, int yr, double *m);
     697extern void mjd_cal (double m, int *mn, double *dy, int *yr);
     698extern int mjd_dow (double m, int *dow);
     699extern int isleapyear (int year);
     700extern void mjd_dpm (double m, int *ndays);
     701extern void mjd_year (double m, double *yr);
     702extern void year_mjd (double y, double *m);
     703extern void rnd_second (double *t);
     704extern void mjd_dayno (double jd, int *yr, double *dy);
     705extern double mjd_day (double jd);
     706extern double mjd_hr (double jd);
     707extern void range (double *v, double r);
     708extern void radecrange (double *ra, double *dec);
    130709
    131710/* moon.c */
    132 extern void moon P_((double mjd, double *lam, double *bet, double *rho,
    133     double *msp, double *mdp));
     711extern void moon (double m, double *lam, double *bet, double *rho,
     712    double *msp, double *mdp);
    134713
    135714/* mooncolong.c */
    136 extern void moon_colong P_((double jd, double lt, double lg, double *cp, double *kp, double *ap, double *sp));
     715extern void moon_colong (double jd, double lt, double lg, double *cp,
     716    double *kp, double *ap, double *sp);
     717
     718/* moonnf.c */
     719extern void moonnf (double mj, double *mjn, double *mjf);
    137720
    138721/* nutation.c */
    139 extern void nutation P_((double mjd, double *deps, double *dpsi));
    140 extern void nut_eq P_((double mjd, double *ra, double *dec));
     722extern void nutation (double m, double *deps, double *dpsi);
     723extern void nut_eq (double m, double *ra, double *dec);
    141724
    142725/* obliq.c */
    143 extern void obliquity P_((double mjd, double *eps));
     726extern void obliquity (double m, double *eps);
    144727
    145728/* parallax.c */
    146 extern void ta_par P_((double tha, double tdec, double phi, double ht,
    147     double *rho, double *aha, double *adec));
     729extern void ta_par (double tha, double tdec, double phi, double ht,
     730    double *rho, double *aha, double *adec);
     731
     732/* parallactic.c */
     733extern double parallacticLDA (double lt, double dec, double alt);
     734extern double parallacticLHD (double lt, double ha, double dec);
    148735
    149736/* plans.c */
    150 extern void plans P_((double mjd, int p, double *lpd0, double *psi0,
     737extern void plans (double m, PLCode p, double *lpd0, double *psi0,
    151738    double *rp0, double *rho0, double *lam, double *bet, double *dia,
    152     double *mag));
     739    double *mag);
     740
     741/* plshadow.c */
     742extern int plshadow (Now *np, Obj *op, Obj *sop, double polera,
     743double poledec, double x, double y, double z, double *sxp, double *syp);
     744
     745/* plmoon_cir.c */
     746extern int plmoon_cir (Now *np, Obj *moonop);
     747extern int getBuiltInObjs (Obj **opp);
     748extern void setMoonDir (char *dir);
    153749
    154750/* precess.c */
    155 extern void precess P_((double mjd1, double mjd2, double *ra, double *dec));
     751extern void precess (double mjd1, double mjd2, double *ra, double *dec);
    156752
    157753/* reduce.c */
    158 extern void reduce_elements P_((double mjd0, double mjd, double inc0,
    159     double ap0, double om0, double *inc, double *ap, double *om));
     754extern void reduce_elements (double mjd0, double m, double inc0,
     755    double ap0, double om0, double *inc, double *ap, double *om);
    160756
    161757/* refract.c */
    162 extern void unrefract P_((double pr, double tr, double aa, double *ta));
    163 extern void refract P_((double pr, double tr, double ta, double *aa));
     758extern void unrefract (double pr, double tr, double aa, double *ta);
     759extern void refract (double pr, double tr, double ta, double *aa);
    164760
    165761/* rings.c */
    166 extern void satrings P_((double sb, double sl, double sr, double el, double er,
    167     double JD, double *etiltp, double *stiltp));
     762extern void satrings (double sb, double sl, double sr, double el, double er,
     763    double JD, double *etiltp, double *stiltp);
    168764
    169765/* riset.c */
    170 extern void riset P_((double ra, double dec, double lat, double dis,
    171     double *lstr, double *lsts, double *azr, double *azs, int *status));
     766extern void riset (double ra, double dec, double lt, double dis,
     767    double *lstr, double *lsts, double *azr, double *azs, int *status);
     768
     769/* riset_cir.c */
     770extern void riset_cir (Now *np, Obj *op, double dis, RiseSet *rp);
     771extern void twilight_cir (Now *np, double dis, double *dawn, double *dusk,
     772    int *status);
     773
     774/* satmoon.c */
     775extern void saturn_data (double Mjd, char dir[], Obj *eop, Obj *sop,
     776    double *satsize, double *etilt, double *stlit, MoonData md[S_NMOONS]);
    172777
    173778/* sphcart.c */
    174 extern void sphcart P_((double l, double b, double r, double *x, double *y,
    175     double *z));
    176 extern void cartsph P_((double x, double y, double z, double *l, double *b,
    177     double *r));
     779extern void sphcart (double l, double b, double r, double *x, double *y,
     780    double *z);
     781extern void cartsph (double x, double y, double z, double *l, double *b,
     782    double *r);
    178783
    179784/* sun.c */
    180 extern void sunpos P_((double mjd, double *lsn, double *rsn, double *bsn));
     785extern void sunpos (double m, double *lsn, double *rsn, double *bsn);
    181786
    182787/* twobody.c */
    183 extern void vrc P_((double *v, double *r, double tp, double e, double q));
     788extern int vrc (double *v, double *r, double tp, double e, double q);
     789
     790/* umoon.c */
     791extern void uranus_data (double Mjd, char dir[], Obj *eop, Obj *uop,
     792    double *usize, MoonData md[U_NMOONS]);
    184793
    185794/* utc_gst.c */
    186 extern void utc_gst P_((double mjd, double utc, double *gst));
    187 extern void gst_utc P_((double mjd, double gst, double *utc));
     795extern void utc_gst (double m, double utc, double *gst);
     796extern void gst_utc (double m, double gst, double *utc);
    188797
    189798/* vsop87.c */
    190 extern int vsop87 P_((double mjd, int obj, double prec, double *ret));
     799extern int vsop87 (double m, int obj, double prec, double *ret);
     800
     801#endif /* _ASTRO_H */
    191802
    192803/* For RCS Only -- Do Not Edit
    193  * @(#) $RCSfile: astro.h,v $ $Date: 2001-10-22 12:08:26 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $
     804 * @(#) $RCSfile: astro.h,v $ $Date: 2004-06-15 16:52:37 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $
    194805 */

  • trunk/SophyaExt/XephemAstroLib/auxil.c

    r1719 r2551  
    33
    44#include <stdio.h>
     5#include <stdlib.h>
    56
    6 #if defined (__STDC__)
    7 #include <stdlib.h>
    8 #endif
    9 
    10 #include "P_.h"
    117#include "astro.h"
    12 #include "circum.h"
    138#include "preferences.h"
    149
     10/* default preferences */
    1511static int prefs[NPREFS] = {
    1612    PREF_TOPO, PREF_METRIC, PREF_MDY, PREF_UTCTZ, PREF_HIPREC, PREF_NOMSGBELL,
    17     PREF_PREFILL, PREF_TIPSON, PREF_CONFIRMON, PREF_WEEKSTART
     13    PREF_PREFILL, PREF_TIPSON, PREF_CONFIRMON, PREF_SUN
    1814};
    1915
    … …  
    2117 */
    2218int
    23 pref_get(pref)
    24 Preferences pref;
     19pref_get(Preferences pref)
    2520{
    2621        return (prefs[pref]);
    … …  
    3025 */
    3126int
    32 pref_set (pref, new)
    33 Preferences pref;
    34 int new;
     27pref_set (Preferences pref, int newp)
    3528{
    3629        int prior = pref_get(pref);
    37         prefs[pref] = new;
     30        prefs[pref] = newp;
    3831        return (prior);
    3932}
    … …  
    4134/* given an mjd, return it modified for terrestial dynamical time */
    4235double
    43 mm_mjed (np)
    44 Now *np;
     36mm_mjed (Now *np)
    4537{
    4638        return (mjd + deltat(mjd)/86400.0);
    … …  
    4840
    4941/* For RCS Only -- Do Not Edit */
    50 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: auxil.c,v $ $Date: 2001-10-22 12:08:26 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     42static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: auxil.c,v $ $Date: 2004-06-15 16:52:37 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/chap95.c

    r1719 r2551  
    1919
    2020#include <math.h>
    21 #include "P_.h"
     21
    2222#include "astro.h"
    2323#include "chap95.h"
    24 
    25 extern void zero_mem P_((void *loc, unsigned len));
    2624
    2725#define CHAP_MAXTPOW    2       /* NB: valid for all 5 outer planets */
    … …  
    3028 *
    3129 * input:
    32  *      mjd     modified JD; days from J1900.0 = 2415020.0
     30 *      m       modified JD; days from J1900.0 = 2415020.0
    3331 *
    3432 *      prec    precision level, in radians.
    … …  
    5048 */
    5149int
    52 chap95 (mjd, obj, prec, ret)
    53 double mjd;
    54 int obj;
    55 double prec;
    56 double *ret;
     50chap95 (double m, int obj, double prec, double *ret)
    5751{
    5852        static double a0[] = {          /* semimajor axes for precision ctrl */
    … …  
    6761
    6862        /* check parameters */
    69         if (mjd < CHAP_BEGIN || mjd > CHAP_END)
     63        if (m < CHAP_BEGIN || m > CHAP_END)
    7064                return (1);
    7165
    … …  
    7973        zero_mem ((void *)sum, sizeof(sum));
    8074
    81         T = (mjd - J2000)/36525.0;      /* centuries since J2000.0 */
     75        T = (m - J2000)/36525.0;        /* centuries since J2000.0 */
    8276
    8377        /* modify precision treshold for
    … …  
    178172
    179173/* For RCS Only -- Do Not Edit */
    180 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: chap95.c,v $ $Date: 2001-10-22 12:08:26 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     174static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: chap95.c,v $ $Date: 2004-06-15 16:52:37 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/chap95.h

    r1719 r2551  
    6161extern chap95_rec chap95_pluto[];
    6262
    63 #include "P_.h"
    64 extern int chap95 P_((double mjd, int obj, double prec, double *ret));
     63extern int chap95 (double m, int obj, double prec, double *ret);
    6564
    6665
    6766/* For RCS Only -- Do Not Edit
    68  * @(#) $RCSfile: chap95.h,v $ $Date: 2001-10-22 12:08:26 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $
     67 * @(#) $RCSfile: chap95.h,v $ $Date: 2004-06-15 16:52:37 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $
    6968 */

  • trunk/SophyaExt/XephemAstroLib/chap95_data.c

    r1719 r2551  
    781781
    782782/* For RCS Only -- Do Not Edit */
    783 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: chap95_data.c,v $ $Date: 2001-10-22 12:08:26 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     783static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: chap95_data.c,v $ $Date: 2004-06-15 16:52:37 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/circum.c

    r1719 r2551  
    99#include <stdio.h>
    1010#include <math.h>
    11 #if defined(__STDC__)
    1211#include <stdlib.h>
    13 #endif
    14 
    15 #include "P_.h"
     12
    1613#include "astro.h"
    17 #include "circum.h"
    1814#include "preferences.h"
    1915
    2016
    21 static int obj_planet P_((Now *np, Obj *op));
    22 static int obj_fixed P_((Now *np, Obj *op));
    23 static int obj_elliptical P_((Now *np, Obj *op));
    24 static int obj_hyperbolic P_((Now *np, Obj *op));
    25 static int obj_parabolic P_((Now *np, Obj *op));
    26 static int sun_cir P_((Now *np, Obj *op));
    27 static int moon_cir P_((Now *np, Obj *op));
    28 static void cir_sky P_((Now *np, double lpd, double psi, double rp, double *rho,
    29     double lam, double bet, double lsn, double rsn, Obj *op));
    30 static void cir_pos P_((Now *np, double bet, double lam, double *rho, Obj *op));
    31 static void elongation P_((double lam, double bet, double lsn, double *el));
    32 static void deflect P_((double mjd1, double lpd, double psi, double rsn,
    33     double lsn, double rho, double *ra, double *dec));
    34 static double h_albsize P_((double H));
     17static int obj_planet (Now *np, Obj *op);
     18static int obj_binary (Now *np, Obj *op);
     19static int obj_2binary (Now *np, Obj *op);
     20static int obj_fixed (Now *np, Obj *op);
     21static int obj_elliptical (Now *np, Obj *op);
     22static int obj_hyperbolic (Now *np, Obj *op);
     23static int obj_parabolic (Now *np, Obj *op);
     24static int sun_cir (Now *np, Obj *op);
     25static int moon_cir (Now *np, Obj *op);
     26static double solveKepler (double M, double e);
     27static void binaryStarOrbit (double t, double T, double e, double o, double O,
     28    double i, double a, double P, double *thetap, double *rhop);
     29static void cir_sky (Now *np, double lpd, double psi, double rp, double *rho,
     30    double lam, double bet, double lsn, double rsn, Obj *op);
     31static void cir_pos (Now *np, double bet, double lam, double *rho, Obj *op);
     32static void elongation (double lam, double bet, double lsn, double *el);
     33static void deflect (double mjd1, double lpd, double psi, double rsn,
     34    double lsn, double rho, double *ra, double *dec);
     35static double h_albsize (double H);
    3536
    3637/* given a Now and an Obj, fill in the approprirate s_* fields within Obj.
    … …  
    3839 */
    3940int
    40 obj_cir (np, op)
    41 Now *np;
    42 Obj *op;
    43 {
     41obj_cir (Now *np, Obj *op)
     42{
     43        op->o_flags &= ~NOCIRCUM;
    4444        switch (op->o_type) {
     45        case BINARYSTAR: return (obj_binary (np, op));
    4546        case FIXED:      return (obj_fixed (np, op));
    4647        case ELLIPTICAL: return (obj_elliptical (np, op));
    … …  
    5051        case PLANET:     return (obj_planet (np, op));
    5152        default:
    52             printf ("obj_cir() called with type %d\n", op->o_type);
    53             exit(1);
     53            printf ("obj_cir() called with type %d %s\n", op->o_type, op->o_name);
     54            abort();
    5455            return (-1);        /* just for lint */
    5556        }
    … …  
    5758
    5859static int
    59 obj_planet (np, op)
    60 Now *np;
    61 Obj *op;
     60obj_planet (Now *np, Obj *op)
    6261{
    6362        double lsn, rsn;        /* true geoc lng of sun; dist from sn to earth*/
    … …  
    6766        double lam, bet;        /* geocentric ecliptic long and lat */
    6867        double dia, mag;        /* angular diameter at 1 AU and magnitude */
    69         int p;
     68        PLCode p;
    7069
    7170        /* validate code and check for a few special cases */
    72         p = op->pl.pl_code;
     71        p = op->pl_code;
     72        if (p == SUN)
     73            return (sun_cir (np, op));
     74        if (p == MOON)
     75            return (moon_cir (np, op));
     76        if (op->pl_moon != X_PLANET)
     77            return (plmoon_cir (np, op));
    7378        if (p < 0 || p > MOON) {
    7479            printf ("unknown planet code: %d\n", p);
    75             exit(1);
    76         }
    77         else if (p == SUN)
    78             return (sun_cir (np, op));
    79         else if (p == MOON)
    80             return (moon_cir (np, op));
     80            abort();
     81        }
     82
     83        /* planet itself */
    8184
    8285        /* find solar ecliptical longitude and distance to sun from earth */
    … …  
    99102
    100103static int
    101 obj_fixed (np, op)
    102 Now *np;
    103 Obj *op;
     104obj_binary (Now *np, Obj *op)
     105{
     106        /* always compute circumstances of primary */
     107        if (obj_fixed (np, op) < 0)
     108            return (0);
     109
     110        /* compute secondary only if requested, and always reset request flag */
     111        if (!op->b_2compute)
     112            return (0);
     113        op->b_2compute = 0;
     114        return (obj_2binary (np, op));
     115}
     116
     117/* compute position of secondary component of a BINARYSTAR */
     118static int
     119obj_2binary (Now *np, Obj *op)
     120{
     121        if (op->b_nbp > 0) {
     122            /* we just have discrete pa/sep, project each from primary */
     123            int i;
     124            for (i = 0; i < op->b_nbp; i++) {
     125                BinPos *bp = &op->b_bp[i];
     126                bp->bp_dec = op->s_dec + bp->bp_sep*cos(bp->bp_pa);
     127                bp->bp_ra = op->s_ra + bp->bp_sep*sin(bp->bp_pa)/cos(op->s_dec);
     128            }
     129        } else {
     130            BinOrbit *bp = &op->b_bo;
     131            double t, theta, rho;
     132
     133            mjd_year (mjd, &t);
     134            binaryStarOrbit (t, bp->bo_T, bp->bo_e, bp->bo_o, bp->bo_O,
     135                                bp->bo_i, bp->bo_a, bp->bo_P, &theta, &rho);
     136            bp->bo_pa = (float)theta;
     137            bp->bo_sep = (float)rho;
     138            rho = degrad(rho/3600.);    /* arc secs to rads */
     139            bp->bo_dec = op->s_dec + rho*cos(theta);
     140            bp->bo_ra =  op->s_ra  + rho*sin(theta)/cos(op->s_dec);
     141        }
     142
     143        return (0);
     144}
     145
     146/* from W. M. Smart */
     147static void
     148binaryStarOrbit (
     149double t,               /* desired ephemeris epoch, year */
     150double T,               /* epoch of periastron, year */
     151double e,               /* eccentricity */
     152double o,               /* argument of periastron, degrees */
     153double O,               /* ascending node, degrees */
     154double i,               /* inclination, degrees */
     155double a,               /* semi major axis, arcsecs */
     156double P,               /* period, years */
     157double *thetap,         /* position angle, rads E of N */
     158double *rhop)           /* separation, arcsecs */
     159{
     160        double M, E, cosE, nu, cosnu, r, rho, theta;
     161
     162        /* find mean anomaly, insure 0..2*PI */
     163        M = 2*PI/P*(t-T);
     164        range (&M, 2*PI);
     165
     166        /* solve for eccentric anomaly */
     167        E = solveKepler (M, e);
     168        cosE = cos(E);
     169
     170        /* find true anomaly and separation */
     171        cosnu = (cosE - e)/(1.0 - e*cosE);
     172        r = a*(1.0 - e*e)/(1.0 + e*cosnu);
     173        nu = acos(cosnu);
     174        if (E > PI)
     175            nu = -nu;
     176
     177        /* project onto sky */
     178        theta = atan(tan(nu+degrad(o))*cos(degrad(i))) + degrad(O);
     179        rho = r*cos(nu+degrad(o))/cos(theta-degrad(O));
     180        if (rho < 0) {
     181            theta += PI;
     182            rho = -rho;
     183        }
     184        range (&theta, 2*PI);
     185
     186        *thetap = theta;
     187        *rhop = rho;
     188}
     189
     190/* solve kepler equation using Newton-Raphson search.
     191 * Charles and Tatum have shown it always converges starting with PI.
     192 */
     193static double
     194solveKepler (double M, double e)
     195{
     196        double E, Eprime = PI;
     197
     198        do {
     199            double cosE = cos(Eprime);
     200            E = Eprime;
     201            Eprime = (M - e*(E*cosE - sin(E)))/(1.0 - e*cosE);
     202        } while (fabs(E-Eprime) > 1e-7);
     203
     204        return (Eprime);
     205}
     206
     207static int
     208obj_fixed (Now *np, Obj *op)
    104209{
    105210        double lsn, rsn;        /* true geoc lng of sun, dist from sn to earth*/
    … …  
    108213        double el;              /* elongation */
    109214        double alt, az;         /* current alt, az */
    110         double ra, dec;         /* ra and dec at epoch of date */
     215        double ra, dec;         /* ra and dec at equinox of date */
     216        double rpm, dpm;        /* astrometric ra and dec with PM to now */
    111217        double lst;
    112218
    113         if (epoch != EOD && (float)epoch != op->f_epoch) {
    114             /* want a certain epoch -- if it's not what the database is at
    115              * we change the original to save time next time assuming the
    116              * user is likely to stick with this for a while.
    117              */
    118             double tra = op->f_RA, tdec = op->f_dec;
    119             float tepoch = (float)epoch;        /* compare w/float precision */
    120             precess (op->f_epoch, tepoch, &tra, &tdec);
    121             op->f_epoch = tepoch;
    122             op->f_RA = (float)tra;
    123             op->f_dec = (float)tdec;
    124         }
    125 
    126         /* set ra/dec to astrometric @ epoch of date */
    127         ra = op->f_RA;
    128         dec = op->f_dec;
     219        /* on the assumption that the user will stick with their chosen display
     220         * epoch for a while, we move the defining values to match and avoid
     221         * precession for every call until it is changed again.
     222         * N.B. only compare and store jd's to lowest precission (f_epoch).
     223         * N.B. maintaining J2k ref (which is arbitrary) helps avoid accum err
     224         */
     225        if (epoch != EOD && (float)epoch != (float)op->f_epoch) {
     226            double pr = op->f_RA, pd = op->f_dec, fe = (float)epoch;
     227            /* first bring back to 2k */
     228            precess (op->f_epoch, J2000, &pr, &pd);
     229            pr += op->f_pmRA*(J2000-op->f_epoch);
     230            pd += op->f_pmdec*(J2000-op->f_epoch);
     231            /* then to epoch */
     232            pr += op->f_pmRA*(fe-J2000);
     233            pd += op->f_pmdec*(fe-J2000);
     234            precess (J2000, fe, &pr, &pd);
     235            op->f_RA = (float)pr;
     236            op->f_dec = (float)pd;
     237            op->f_epoch = (float)fe;
     238        }
     239
     240        /* apply proper motion .. assume pm epoch reference equals equinox */
     241        rpm = op->f_RA + op->f_pmRA*(mjd-op->f_epoch);
     242        dpm = op->f_dec + op->f_pmdec*(mjd-op->f_epoch);
     243
     244        /* set ra/dec to astrometric @ equinox of date */
     245        ra = rpm;
     246        dec = dpm;
    129247        precess (op->f_epoch, mjed, &ra, &dec);
    130248
    … …  
    154272        } else {
    155273            /* annual parallax at time mjd is to be added here, too, but
    156              * technically in the frame of epoch (usually different from mjd)
     274             * technically in the frame of equinox (usually different from mjd)
    157275             */
    158             op->s_ra = op->f_RA;        /* already precessed */
    159             op->s_dec = op->f_dec;
     276            op->s_ra = rpm;
     277            op->s_dec = dpm;
    160278        }
    161279
    … …  
    184302 */
    185303static int
    186 obj_elliptical (np, op)
    187 Now *np;
    188 Obj *op;
     304obj_elliptical (Now *np, Obj *op)
    189305{
    190306        double lsn, rsn;        /* true geoc lng of sun; dist from sn to earth*/
    … …  
    230346
    231347            tp = mjed - dt - (op->e_cepoch - op->e_M/e_n);
    232             vrc (&nu, &rp, tp, op->e_e, op->e_a*(1-op->e_e));
     348            if (vrc (&nu, &rp, tp, op->e_e, op->e_a*(1-op->e_e)) < 0)
     349                op->o_flags |= NOCIRCUM;
    233350            nu = degrad(nu);
    234351            lo = nu + om;
    … …  
    288405 */
    289406static int
    290 obj_hyperbolic (np, op)
    291 Now *np;
    292 Obj *op;
     407obj_hyperbolic (Now *np, Obj *op)
    293408{
    294409        double lsn, rsn;        /* true geoc lng of sun; dist from sn to earth*/
    … …  
    309424        double e;               /* fast eccentricity */
    310425        double ll=0, sll, cll;  /* helio angle between object and earth */
    311         double n;               /* mean daily motion */
    312426        double mag;             /* magnitude */
    313427        double a;               /* mean distance */
    … …  
    323437        e = op->h_e;
    324438        a = op->h_qp/(e - 1.0);
    325         n = .98563/sqrt(a*a*a);
    326439
    327440        /* correct for light time by computing position at time mjd, then
    … …  
    337450
    338451            tp = mjed - dt - op->h_ep;
    339             vrc (&nu, &rp, tp, op->h_e, op->h_qp);
     452            if (vrc (&nu, &rp, tp, op->h_e, op->h_qp) < 0)
     453                op->o_flags |= NOCIRCUM;
    340454            nu = degrad(nu);
    341455            lo = nu + om;
    … …  
    382496 */
    383497static int
    384 obj_parabolic (np, op)
    385 Now *np;
    386 Obj *op;
     498obj_parabolic (Now *np, Obj *op)
    387499{
    388500        double lsn, rsn;        /* true geoc lng of sun; dist from sn to earth*/
    … …  
    422534 */
    423535static int
    424 sun_cir (np, op)
    425 Now *np;
    426 Obj *op;
     536sun_cir (Now *np, Obj *op)
    427537{
    428538        double lsn, rsn;        /* true geoc lng of sun; dist from sn to earth*/
    … …  
    452562 */
    453563static int
    454 moon_cir (np, op)
    455 Now *np;
    456 Obj *op;
     564moon_cir (Now *np, Obj *op)
    457565{
    458566        double lsn, rsn;        /* true geoc lng of sun; dist from sn to earth*/
    … …  
    480588
    481589        /* TODO: improve mag; this is based on a flat moon model. */
    482         set_smag (op, -12.7 + 2.5*(log10(PI) - log10(PI/2*(1+1.e-6-cos(el)))));
     590        i = -12.7 + 2.5*(log10(PI) - log10(PI/2*(1+1.e-6-cos(el))))
     591                                        + 5*log10(edistau/.0025) /* dist */;
     592        set_smag (op, i);
    483593
    484594        /* find phase -- allow for projection effects */
    … …  
    500610 */
    501611static void
    502 cir_sky (np, lpd, psi, rp, rho, lam, bet, lsn, rsn, op)
    503 Now *np;
    504 double lpd, psi;        /* heliocentric ecliptic long and lat */
    505 double rp;              /* dist from sun */
    506 double *rho;            /* dist from earth: in as geo, back as geo or topo */
    507 double lam, bet;        /* true geocentric ecliptic long and lat */
    508 double lsn, rsn;        /* true geoc lng of sun; dist from sn to earth*/
    509 Obj *op;
     612cir_sky (
     613Now *np,
     614double lpd,             /* heliocentric ecliptic longitude */
     615double psi,             /* heliocentric ecliptic lat */
     616double rp,              /* dist from sun */
     617double *rho,            /* dist from earth: in as geo, back as geo or topo */
     618double lam,             /* true geocentric ecliptic long */
     619double bet,             /* true geocentric ecliptic lat */
     620double lsn,             /* true geoc lng of sun */
     621double rsn,             /* dist from sn to earth*/
     622Obj *op)
    510623{
    511624        double el;              /* elongation */
    … …  
    546659 *   refract    --> alt/az      observed --> output
    547660 *
    548  * algorithm at fixed epoch:
     661 * algorithm at fixed equinox:
    549662 *   ecl_eq     --> ra/dec      geocentric mean equatoreal EOD (via mean obliq)
    550663 *   deflect    --> ra/dec        relativistic deflection [for alt/az only]
    … …  
    558671 */
    559672static void
    560 cir_pos (np, bet, lam, rho, op)
    561 Now *np;
    562 double bet, lam;/* geo lat/long (mean ecliptic of date) */
    563 double *rho;    /* in: geocentric dist in AU; out: geo- or topocentic dist */
    564 Obj *op;        /* object to set s_ra/dec as per epoch */
     673cir_pos (
     674Now *np,
     675double bet,     /* geo lat (mean ecliptic of date) */
     676double lam,     /* geo long (mean ecliptic of date) */
     677double *rho,    /* in: geocentric dist in AU; out: geo- or topocentic dist */
     678Obj *op)        /* object to set s_ra/dec as per equinox */
    565679{
    566680        double ra, dec;         /* apparent ra/dec, corrected for nut/ab */
    … …  
    649763 */
    650764static void
    651 elongation (lam, bet, lsn, el)
    652 double lam, bet, lsn;
    653 double *el;
     765elongation (double lam, double bet, double lsn, double *el)
    654766{
    655767        *el = acos(cos(bet)*cos(lam-lsn));
    … …  
    675787 */
    676788static void
    677 deflect (mjd1, lpd, psi, lsn, rsn, rho, ra, dec)
    678 double mjd1;            /* epoch */
    679 double lpd, psi;        /* heliocentric ecliptical long / lat */
    680 double rsn, lsn;        /* distance and longitude of sun */
    681 double rho;             /* geocentric distance */
    682 double *ra, *dec;       /* geocentric equatoreal */
     789deflect (
     790double mjd1,            /* equinox */
     791double lpd, double psi, /* heliocentric ecliptical long / lat */
     792double rsn, double lsn, /* distance and longitude of sun */
     793double rho,             /* geocentric distance */
     794double *ra, double *dec)/* geocentric equatoreal */
    683795{
    684796        double hra, hdec;       /* object heliocentric equatoreal */
    … …  
    742854 */
    743855static double
    744 h_albsize (H)
    745 double H;
     856h_albsize (double H)
    746857{
    747858        return (3600*raddeg(.707*1500*pow(2.51,(18-H)/2)/MAU));
    … …  
    749860
    750861/* For RCS Only -- Do Not Edit */
    751 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: circum.c,v $ $Date: 2001-10-22 12:08:26 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     862static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: circum.c,v $ $Date: 2004-06-15 16:52:37 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/comet.c

    r1719 r2551  
    11#include <math.h>
    22
    3 #include "P_.h"
    43#include "astro.h"
    54
    6 /* given a modified Julian date, mjd, and a set of heliocentric parabolic
    7  * orbital elements referred to the epoch of date (mjd):
     5/* given a modified Julian date, mj, and a set of heliocentric parabolic
     6 * orbital elements referred to the epoch of date (mj):
    87 *   ep:   epoch of perihelion,
    98 *   inc:  inclination,
    … …  
    3029 */
    3130void
    32 comet (mjd, ep, inc, ap, qp, om, lpd, psi, rp, rho, lam, bet)
    33 double mjd;
    34 double ep, inc, ap, qp, om;
    35 double *lpd, *psi, *rp, *rho, *lam, *bet;
     31comet (double mj, double ep, double inc, double ap, double qp, double om,
     32double *lpd, double *psi, double *rp, double *rho, double *lam, double *bet)
    3633{
    3734        double w, s, s2;
    … …  
    4441
    4542#define ERRLMT  0.0001
    46         w = ((mjd-ep)*3.649116e-02)/(qp*sqrt(qp));
     43        w = ((mj-ep)*3.649116e-02)/(qp*sqrt(qp));
    4744        s = w/3;
    4845        for (;;) {
    … …  
    6865        range (lpd, 2*PI);
    6966        rd = *rp * cpsi;
    70         sunpos (mjd, &lsn, &rsn, 0);
     67        sunpos (mj, &lsn, &rsn, 0);
    7168        lg = lsn+PI;
    7269        re = rsn;
    … …  
    8481
    8582/* For RCS Only -- Do Not Edit */
    86 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: comet.c,v $ $Date: 2001-10-22 12:08:26 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     83static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: comet.c,v $ $Date: 2004-06-15 16:52:38 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/compare_with_xephem.csh

    r1895 r2551  
    44  echo 'Placez vous dans le repertoire ou vous avez'
    55  echo 'mis le NOUVEAU code du libastro de Xephem.'
    6   echo '> $DPCSOURCE/XephemAstroLib/compare_with_xephem.csh'
     6  echo '> $SOPHYASOURCE/XephemAstroLib/compare_with_xephem.csh'
    77  exit -1
    88endif
    99
    10 set dir = $DPCSOURCE/XephemAstroLib/
     10set dir = $SOPHYASOURCE/XephemAstroLib/
    1111set log = `pwd`/compare.log
    1212

  • trunk/SophyaExt/XephemAstroLib/dbfmt.c

    r1719 r2551  
    44#include <ctype.h>
    55#include <math.h>
    6 
    7 #if defined(__STDC__)
    86#include <stdlib.h>
    97#include <string.h>
    10 #endif
    11 
    12 #include "P_.h"
     8
    139#include "astro.h"
    14 #include "circum.h"
    1510#include "preferences.h"
    1611
    17 extern void zero_mem P_((void *loc, unsigned len));
    18 extern double atod P_((char *buf));
    19 
    20 int get_fields P_((char *s, int delim, char *fields[]));
    21 int db_set_field P_((char bp[], int id, PrefDateFormat pref, Obj *op));
    22 int db_chk_planet P_((char name[], Obj *op));
     12
     13int get_fields (char *s, int delim, char *fields[]);
    2314
    2415#define MAXDBLINE       256     /* longest allowed db line */
    … …  
    2819#define MAXFLDS 20              /* must be more than on any expected line */
    2920
    30 #define ASIZ(a) (sizeof(a)/sizeof(a[0]))
    31 
    32 static int line_candidate P_((char *buf));
    33 static void crack_year P_((char *bp, PrefDateFormat pref, double *p));
    34 static int db_get_field P_((Obj *op, int id, char *bp));
    35 static int tle_sum P_((char *l));
    36 static double tle_fld P_((char *l, int from, int thru));
    37 static double tle_expfld P_((char *l, int start));
     21static char *enm (char *flds[MAXFLDS]);
     22static int crack_f (Obj *op, char *flds[MAXFLDS], int nf, char whynot[]);
     23static int crack_e (Obj *op, char *flds[MAXFLDS], int nf, char whynot[]);
     24static int crack_h (Obj *op, char *flds[MAXFLDS], int nf, char whynot[]);
     25static int crack_p (Obj *op, char *flds[MAXFLDS], int nf, char whynot[]);
     26static int crack_E (Obj *op, char *flds[MAXFLDS], int nf, char whynot[]);
     27static int crack_P (Obj *op, char *flds[MAXFLDS], int nf, char whynot[]);
     28static int crack_B (Obj *op, char *flds[MAXFLDS], int nf, char whynot[]);
     29static int crack_name (Obj *op, char *flds[MAXFLDS], int nf,
     30    char nm[][MAXNM], int nnm);
     31static void crack_year (char *bp, double *p);
     32static void crack_okdates (char *fld, float *startok, float *endok);
     33static int get_okdates (char *lp, float *sp, float *ep);
     34static int tle_sum (char *l);
     35static double tle_fld (char *l, int from, int thru);
     36static double tle_expfld (char *l, int start);
     37static void write_f (Obj *op, char lp[]);
     38static void write_e (Obj *op, char lp[]);
     39static void write_h (Obj *op, char lp[]);
     40static void write_p (Obj *op, char lp[]);
     41static void write_E (Obj *op, char lp[]);
     42static void write_P (Obj *op, char lp[]);
     43static void write_B (Obj *op, char lp[]);
    3844
    3945/* crack the given .edb database line into op.
    4046 * if ok
    41  *   return 0
     47 *   return number of names in nm[], or 1 if nm == NULL
    4248 * else
    4349 *   if whynot
    … …  
    4753 *       fill whynot with reason message.
    4854 *   return -1
     55 * only the first name is stored in op, all names (up to nnm) are in nm[], or
     56 *   ignored if nm == NULL.
    4957 */
    5058int
    51 db_crack_line (s, op, whynot)
    52 char s[];
    53 Obj *op;
    54 char whynot[];
    55 {
     59db_crack_line (char s[], Obj *op, char nm[][MAXNM], int nnm, char whynot[])
     60{
     61        char copy[MAXDBLINE];   /* work copy; leave s untouched */
    5662        char *flds[MAXFLDS];    /* point to each field for easy reference */
    57         char *sflds[MAXFLDS];   /* point to each sub field for easy reference */
    58         char copy[MAXDBLINE];   /* work copy; leave s untouched */
    59         int nf, nsf;            /* number of fields and subfields */
     63        int nf;
    6064        int i;
    6165
     66        /* init no response */
     67        if (whynot)
     68            whynot[0] = '\0';
     69
    6270        /* basic initial check */
    63         if (line_candidate (s) < 0) {
    64             if (whynot)
    65                 whynot[0] = '\0';
    66             return (-1);
    67         }
     71        if (dbline_candidate (s) < 0)
     72            return (-1);
    6873
    6974        /* do all the parsing on a copy */
    … …  
    7984        if (nf < 2) {
    8085            if (whynot)
    81                 sprintf (whynot, "Found only %d fields", nf);
     86                sprintf (whynot, "Bogus: %s", s);
    8287            return (-1);
    8388        }
    … …  
    8590        /* switch out on type of object - the second field */
    8691        switch (flds[1][0]) {
    87         case 'f': {
    88             static int ids[] = {F_RA, F_DEC, F_MAG};
    89             if (nf < 5 || nf > 7) {
    90                 if (whynot)
    91                     sprintf (whynot, "f needs 5-7 fields: %d", nf);
    92                 return (-1);
    93             }
    94             zero_mem ((void *)op, sizeof(ObjF));
    95             op->o_type = FIXED;
    96             nsf = get_fields(flds[1], SUBFLD, sflds);
    97             if (nsf > 1 && db_set_field (sflds[1], F_CLASS, PREF_MDY, op) < 0) {
    98                 if (whynot)
    99                     sprintf (whynot, "Bad f class: %c", sflds[1][0]);
    100                 return (-1);
    101             }
    102             if (nsf > 2)
    103                 (void) db_set_field (sflds[2], F_SPECT, PREF_MDY, op);
    104             for (i = 2; i < ASIZ(ids)+2; i++)
    105                 (void) db_set_field (flds[i], ids[i-2], PREF_MDY, op);
    106             (void) db_set_field (nf>5 && flds[5][0] ? flds[5] : "2000",
    107                                                         F_EPOCH, PREF_MDY, op);
    108             if (nf == 7)
    109                 (void) db_set_field (flds[6], F_SIZE, PREF_MDY, op);
    110             break;
    111         }
    112 
    113         case 'e': {
    114             static int ids[] = {E_INC, E_LAN, E_AOP, E_A, E_N, E_E, E_M,
    115                                                 E_CEPOCH, E_EPOCH, E_M1, E_M2
    116             };
    117             if (nf != 13 && nf != 14) {
    118                 if (whynot)
    119                     sprintf (whynot, "e needs 13 or 14 fields: %d", nf);
    120                 return (-1);
    121             }
    122             zero_mem ((void *)op, sizeof(ObjE));
    123             op->o_type = ELLIPTICAL;
    124             for (i = 2; i < ASIZ(ids)+2; i++)
    125                 (void) db_set_field (flds[i], ids[i-2], PREF_MDY, op);
    126             if (nf == 14)
    127                 (void) db_set_field (flds[13], E_SIZE, PREF_MDY, op);
    128             break;
    129         }
    130 
    131         case 'h': {
    132             static int ids[]= {H_EP,H_INC,H_LAN,H_AOP,H_E,H_QP,H_EPOCH,H_G,H_K};
    133             if (nf != 11 && nf != 12) {
    134                 if (whynot)
    135                     sprintf (whynot, "h needs 11 or 12 fields: %d", nf);
    136                 return (-1);
    137             }
    138             zero_mem ((void *)op, sizeof(ObjH));
    139             op->o_type = HYPERBOLIC;
    140             for (i = 2; i < ASIZ(ids)+2; i++)
    141                 (void) db_set_field (flds[i], ids[i-2], PREF_MDY, op);
    142             if (nf == 12)
    143                 (void) db_set_field (flds[11], H_SIZE, PREF_MDY, op);
    144             break;
    145         }
    146 
    147         case 'p': {
    148             static int ids[] = {P_EP,P_INC,P_AOP,P_QP,P_LAN,P_EPOCH,P_G,P_K};
    149             if (nf != 10 && nf != 11) {
    150                 if (whynot)
    151                     sprintf (whynot, "p needs 10 or 11 fields: %d", nf);
    152                 return (-1);
    153             }
    154             zero_mem ((void *)op, sizeof(ObjP));
    155             op->o_type = PARABOLIC;
    156             for (i = 2; i < ASIZ(ids)+2; i++)
    157                 (void) db_set_field (flds[i], ids[i-2], PREF_MDY, op);
    158             if (nf == 11)
    159                 (void) db_set_field (flds[10], P_SIZE, PREF_MDY, op);
    160             break;
    161         }
    162 
    163         case 'E': {
    164             static int ids[] = {ES_EPOCH,ES_INC,ES_RAAN,ES_E,ES_AP,ES_M,ES_N,
    165                                                             ES_DECAY,ES_ORBIT};
    166             if (nf != 11 && nf != 12) {
    167                 if (whynot)
    168                     sprintf (whynot, "E needs 11 or 12 fields: %d", nf);
    169                 return (-1);
    170             }
    171             zero_mem ((void *)op, sizeof(ObjES));
    172             op->o_type = EARTHSAT;
    173             for (i = 2; i < ASIZ(ids)+2; i++)
    174                 (void) db_set_field (flds[i], ids[i-2], PREF_MDY, op);
    175             if (nf == 12)
    176                 (void) db_set_field (flds[11], ES_DRAG, PREF_MDY, op);
    177             break;
    178         }
     92
     93        case 'f':
     94            if (crack_f (op, flds, nf, whynot) < 0)
     95                return (-1);
     96            break;
     97
     98        case 'e':
     99            if (crack_e (op, flds, nf, whynot) < 0)
     100                return (-1);
     101            break;
     102
     103        case 'h':
     104            if (crack_h (op, flds, nf, whynot) < 0)
     105                return (-1);
     106            break;
     107
     108        case 'p':
     109            if (crack_p (op, flds, nf, whynot) < 0)
     110                return (-1);
     111            break;
     112
     113        case 'B':
     114            if (crack_B (op, flds, nf, whynot) < 0)
     115                return (-1);
     116            break;
     117
     118        case 'E':
     119            if (crack_E (op, flds, nf, whynot) < 0)
     120                return (-1);
     121            break;
    179122
    180123        case 'P':
    181             /* allow, though ignore, anything after the P */
    182             i = db_chk_planet (flds[0], op);    /* does o_name too */
    183             if (i < 0) {
    184                 if (whynot)
    185                     sprintf (whynot, "Bad planet: %s", flds[0]);
    186             }
    187             return (i);
     124            if (crack_P (op, flds, nf, whynot) < 0)
     125                return (-1);
     126            break;
    188127
    189128        default:
    190129            if (whynot)
    191                 sprintf (whynot, "Unknown type: %c", flds[1][0]);
    192             return (-1);
    193         }
    194 
    195         /* load up o_name */
    196         (void) db_set_field (flds[0], O_NAME, PREF_MDY, op);
    197 
    198         return (0);
    199 }
    200 
    201 /* return 0 if TLE checksum is ok, else -1 */
    202 static int
    203 tle_sum (l)
    204 char *l;
    205 {
    206         char *lastl = l + 68;
    207         int sum;
    208 
    209         for (sum = 0; l < lastl; ) {
    210             char c = *l++;
    211             if (c == '\0')
    212                 return (-1);
    213             if (isdigit(c))
    214                 sum += c - '0';
    215             else if (c == '-')
    216                 sum++;
    217         }
    218 
    219         return (*l - '0' == (sum%10) ? 0 : -1);
    220 }
    221 
    222 /* extract the given columns and return value.
    223  * N.B. from and to are 1-based within l
     130                sprintf (whynot, "%s: Unknown type %c for %s", enm(flds),
     131                                                        flds[1][0], flds[0]);
     132            return (-1);
     133        }
     134
     135        return (crack_name (op, flds, nf, nm, nnm));
     136}
     137
     138/* write the given Obj in .edb format to lp[].
     139 * we do _not_ include a trailing '\n'.
    224140 */
    225 static double
    226 tle_fld (l, from, thru)
    227 char *l;
    228 int from, thru;
    229 {
    230         char buf[32];
    231 
    232         sprintf (buf, "%.*s", thru-from+1, l+from-1);
    233         return (atod (buf));
    234 }
    235 
    236 /* extract the exponential value starting at the given column.
    237  * N.B. start is 1-based within l
    238  */
    239 static double
    240 tle_expfld (l, start)
    241 char *l;
    242 int start;
    243 {
    244         char buf[32];
    245         double v;
    246 
    247         sprintf (buf, ".%.*s", 5, l+start);
    248         v = atod (buf) * pow (10.0, tle_fld(l, start+6, start+7));
    249         if (l[start-1] == '-')
    250             v = -v;
    251         return (v);
     141void
     142db_write_line (Obj *op, char lp[])
     143{
     144        int priorpref;
     145
     146        /* .edb format always uses MDY.
     147         * N.B. must restore old value before returning from here!
     148         */
     149        priorpref = pref_set (PREF_DATE_FORMAT, PREF_MDY);
     150
     151        switch (op->o_type) {
     152        case FIXED:
     153            write_f (op, lp);
     154            break;
     155
     156        case BINARYSTAR:
     157            write_B (op, lp);
     158            break;
     159
     160        case ELLIPTICAL:
     161            write_e (op, lp);
     162            break;
     163
     164        case HYPERBOLIC:
     165            write_h (op, lp);
     166            break;
     167
     168        case PARABOLIC:
     169            write_p (op, lp);
     170            break;
     171
     172        case EARTHSAT:
     173            write_E (op, lp);
     174            break;
     175
     176        case PLANET:
     177            write_P (op, lp);
     178            break;
     179
     180        default:
     181            printf ("Unknown type for %s: %d\n", op->o_name, op->o_type);
     182            abort();
     183        }
     184
     185        /* restore date format preference */
     186        (void) pref_set (PREF_DATE_FORMAT, priorpref);
    252187}
    253188
    … …  
    261196 */
    262197int
    263 db_tle (name, l1, l2, op)
    264 char *name, *l1, *l2;
    265 Obj *op;
     198db_tle (char *name, char *l1, char *l2, Obj *op)
    266199{
    267200        double ep;
    … …  
    306239        /* goodies from "line 1" */
    307240        op->es_drag = (float) tle_expfld (l1, 54);
     241        op->es_decay = (float) tle_fld (l1, 34, 43);
    308242        i = (int) tle_fld (l1, 19, 20);
    309243        if (i < 57)
    … …  
    325259}
    326260
    327 /* write the given Obj in .edb format to lp[].
    328  * we do _not_ include a trailing '\n'.
    329  */
    330 void
    331 db_write_line (op, lp)
    332 Obj *op;
    333 char *lp;
    334 {
    335         int priorpref;
    336         int i;
    337 
    338         /* .edb format always uses MDY.
    339          * N.B. must restore old value before returning from here!
    340          */
    341         priorpref = pref_set (PREF_DATE_FORMAT, PREF_MDY);
    342 
    343         switch (op->o_type) {
    344         case FIXED: {
    345             static int ids[] = {F_CLASS, F_SPECT, F_RA, F_DEC, F_MAG, F_EPOCH,
    346                 F_SIZE
    347             };
    348 
    349             sprintf (lp, "%s,f", op->o_name);
    350             lp += strlen(lp);
    351             for (i = 0; i < ASIZ(ids); i++)
    352                 lp += db_get_field (op, ids[i], lp);
    353             break;
    354             }
    355 
    356         case ELLIPTICAL: {
    357             static int ids[] = {E_INC, E_LAN, E_AOP, E_A, E_N, E_E, E_M,
    358                                         E_CEPOCH, E_EPOCH, E_M1, E_M2, E_SIZE
    359             };
    360 
    361             sprintf (lp, "%s,e", op->o_name);
    362             lp += strlen(lp);
    363             for (i = 0; i < ASIZ(ids); i++)
    364                 lp += db_get_field (op, ids[i], lp);
    365             break;
    366             }
    367 
    368         case HYPERBOLIC: {
    369             static int ids[]= {H_EP, H_INC, H_LAN, H_AOP, H_E, H_QP, H_EPOCH,
    370                 H_G, H_K, H_SIZE
    371             };
    372 
    373             sprintf (lp, "%s,h", op->o_name);
    374             lp += strlen(lp);
    375             for (i = 0; i < ASIZ(ids); i++)
    376                 lp += db_get_field (op, ids[i], lp);
    377             break;
    378             }
    379 
    380         case PARABOLIC: {
    381             static int ids[] = {P_EP, P_INC, P_AOP, P_QP, P_LAN, P_EPOCH, P_G,
    382                 P_K, P_SIZE
    383             };
    384 
    385             sprintf (lp, "%s,p", op->o_name);
    386             lp += strlen(lp);
    387             for (i = 0; i < ASIZ(ids); i++)
    388                 lp += db_get_field (op, ids[i], lp);
    389             break;
    390             }
    391 
    392         case EARTHSAT: {
    393             static int ids[] = {ES_EPOCH, ES_INC, ES_RAAN, ES_E, ES_AP, ES_M,
    394                 ES_N, ES_DECAY,ES_ORBIT,ES_DRAG
    395             };
    396 
    397             sprintf (lp, "%s,E", op->o_name);
    398             lp += strlen(lp);
    399             for (i = 0; i < ASIZ(ids); i++)
    400                 lp += db_get_field (op, ids[i], lp);
    401             break;
    402             }
    403 
    404         case PLANET:
    405             sprintf (lp, "%s,P", op->o_name);
    406             lp += strlen(lp);
    407             break;
    408 
    409         default:
    410             printf ("Unknown type for %s: %d\n", op->o_name, op->o_type);
    411             exit(1);
    412         }
    413 
    414         /* restore date format preference */
    415         (void) pref_set (PREF_DATE_FORMAT, priorpref);
    416 }
    417 
    418 /* given a text buffer and a field id, and a PREF_DATE_FORMAT,
    419  *   set the corresponding member in *op.
    420  * return 0 if ok, else -1.
     261/* return 0 if op has no date range information or what it does have brackets
     262 * now, else -1
    421263 */
    422264int
    423 db_set_field (bp, id, pref, op)
    424 char bp[];
    425 int id;
    426 PrefDateFormat pref;
    427 Obj *op;
    428 {
    429         double tmp;
    430 
    431         /* include all the enums and in numeric order to give us the best
    432          * possible chance the compiler will implement this as a jump table.
    433          */
    434         switch (id) {
    435         case O_TYPE:
    436             printf ("db_set_field: called with id==O_TYPE\n");
    437             exit(1);
    438             break;
    439         case O_NAME:
    440             (void) strncpy (op->o_name, bp, sizeof(op->o_name)-1);
    441             op->o_name[sizeof(op->o_name)-1] = '\0';
    442             break;
    443         case F_RA:
    444             f_scansex (radhr(op->f_RA), bp, &tmp);
    445             op->f_RA = (float) hrrad(tmp);
    446             break;
    447         case F_DEC:
    448             f_scansex (raddeg(op->f_dec), bp, &tmp);
    449             op->f_dec = (float) degrad(tmp);
    450             break;
    451         case F_EPOCH:
    452             tmp = op->f_epoch;
    453             crack_year (bp, pref, &tmp);
    454             op->f_epoch = (float) tmp;
    455             break;
    456         case F_MAG:
    457             set_fmag (op, atod(bp));
    458             break;
    459         case F_SIZE:
    460             op->f_size = (float) atod(bp);
    461             {
    462                 /* optional minor axis and position angle subfields */
    463                 char *sflds[MAXFLDS];
    464                 int nsf = get_fields(bp, SUBFLD, sflds);
    465 
    466                 if (nsf == 3) {
    467                     set_ratio(op, op->s_size, atod(sflds[1]));
    468                     set_pa(op,degrad(atod(sflds[2])));
    469                 } else {
    470                     set_ratio(op,1,1);  /* round */
    471                     set_pa(op,0.0);
    472                 }
    473             }
    474             break;
    475         case F_CLASS:
    476             switch (bp[0]) {
    477             case 'A': case 'B': case 'C': case 'D': case 'F': case 'G':
    478             case 'H': case 'K': case 'J': case 'L': case 'M': case 'N':
    479             case 'O': case 'P': case 'Q': case 'R': case 'S': case 'T':
    480             case 'U': case 'V':
    481                 op->f_class = bp[0];
    482                 break;
    483             default:
    484                 return (-1);
    485             }
    486             break;
    487         case F_SPECT: {
    488             int i, j;
    489             /* fill f_spect all the way */
    490             for (i = j = 0; i < sizeof(op->f_spect); i++)
    491                 if ((op->f_spect[i] = bp[j]) != 0)
    492                     j++;
    493             break;
    494         }
    495 
    496         case E_INC:
    497             op->e_inc = (float) atod (bp);
    498             break;
    499         case E_LAN:
    500             op->e_Om = (float) atod (bp);
    501             break;
    502         case E_AOP:
    503             op->e_om = (float) atod (bp);
    504             break;
    505         case E_A:
    506             op->e_a = (float) atod (bp);
    507             break;
    508         case E_N:
    509             /* retired */
    510             break;
    511         case E_E:
    512             op->e_e = atod (bp);
    513             break;
    514         case E_M:
    515             op->e_M = (float) atod (bp);
    516             break;
    517         case E_CEPOCH:
    518             crack_year (bp, pref, &op->e_cepoch);
    519             break;
    520         case E_EPOCH:
    521             crack_year (bp, pref, &op->e_epoch);
    522             break;
    523         case E_M1:
    524             switch (bp[0]) {
    525             case 'g':
    526                 op->e_mag.whichm = MAG_gk;
    527                 bp++;
    528                 break;
    529             case 'H':
    530                 op->e_mag.whichm = MAG_HG;
    531                 bp++;
    532                 break;
    533             default:
    534                 /* leave type unchanged if no or unrecognized prefix */
    535                 break;
    536             }
    537             op->e_mag.m1 = (float) atod(bp);
    538             break;
    539         case E_M2:
    540             switch (bp[0]) {
    541             case 'k':
    542                 op->e_mag.whichm = MAG_gk;
    543                 bp++;
    544                 break;
    545             case 'G':
    546                 op->e_mag.whichm = MAG_HG;
    547                 bp++;
    548                 break;
    549             default:
    550                 /* leave type unchanged if no or unrecognized prefix */
    551                 break;
    552             }
    553             op->e_mag.m2 = (float) atod(bp);
    554             break;
    555         case E_SIZE:
    556             op->e_size = (float) atod (bp);
    557             break;
    558 
    559         case H_EP:
    560             crack_year (bp, pref, &op->h_ep);
    561             break;
    562         case H_INC:
    563             op->h_inc = (float) atod (bp);
    564             break;
    565         case H_LAN:
    566             op->h_Om = (float) atod (bp);
    567             break;
    568         case H_AOP:
    569             op->h_om = (float) atod (bp);
    570             break;
    571         case H_E:
    572             op->h_e = (float) atod (bp);
    573             break;
    574         case H_QP:
    575             op->h_qp = (float) atod (bp);
    576             break;
    577         case H_EPOCH:
    578             crack_year (bp, pref, &op->h_epoch);
    579             break;
    580         case H_G:
    581             op->h_g = (float) atod (bp);
    582             break;
    583         case H_K:
    584             op->h_k = (float) atod (bp);
    585             break;
    586         case H_SIZE:
    587             op->h_size = (float) atod (bp);
    588             break;
    589 
    590         case P_EP:
    591             crack_year (bp, pref, &op->p_ep);
    592             break;
    593         case P_INC:
    594             op->p_inc = (float) atod (bp);
    595             break;
    596         case P_AOP:
    597             op->p_om = (float) atod (bp);
    598             break;
    599         case P_QP:
    600             op->p_qp = (float) atod (bp);
    601             break;
    602         case P_LAN:
    603             op->p_Om = (float) atod (bp);
    604             break;
    605         case P_EPOCH:
    606             crack_year (bp, pref, &op->p_epoch);
    607             break;
    608         case P_G:
    609             op->p_g = (float) atod (bp);
    610             break;
    611         case P_K:
    612             op->p_k = (float) atod (bp);
    613             break;
    614         case P_SIZE:
    615             op->p_size = (float) atod (bp);
    616             break;
    617 
    618         case ES_EPOCH:
    619             crack_year (bp, pref, &op->es_epoch);
    620             break;
    621         case ES_INC:
    622             op->es_inc = (float) atod (bp);
    623             break;
    624         case ES_RAAN:
    625             op->es_raan = (float) atod (bp);
    626             break;
    627         case ES_E:
    628             op->es_e = (float) atod (bp);
    629             break;
    630         case ES_AP:
    631             op->es_ap = (float) atod (bp);
    632             break;
    633         case ES_M:
    634             op->es_M = (float) atod (bp);
    635             break;
    636         case ES_N:
    637             op->es_n = atod (bp);
    638             break;
    639         case ES_DECAY:
    640             op->es_decay = (float) atod (bp);
    641             break;
    642         case ES_ORBIT:
    643             op->es_orbit = atoi (bp);
    644             break;
    645         case ES_DRAG:
    646             op->es_drag = (float) atod (bp);
    647             break;
    648 
     265dateRangeOK (Now *np, Obj *op)
     266{
     267        float *sp, *ep;
     268
     269        switch (op->o_type) {
     270        case ELLIPTICAL:
     271            sp = &op->e_startok;
     272            ep = &op->e_endok;
     273            break;
     274        case HYPERBOLIC:
     275            sp = &op->h_startok;
     276            ep = &op->h_endok;
     277            break;
     278        case PARABOLIC:
     279            sp = &op->p_startok;
     280            ep = &op->p_endok;
     281            break;
     282        case EARTHSAT:
     283            sp = &op->es_startok;
     284            ep = &op->es_endok;
     285            break;
    649286        default:
    650             printf ("BUG! db_set_field: bad id: %d\n", id);
    651             exit (1);
    652         }
    653 
    654         return (0);
     287            return (0);
     288        }
     289
     290        if (*sp <= mjd && (!*ep || mjd <= *ep))
     291            return (0);
     292        return (-1);
    655293}
    656294
    … …  
    662300 */
    663301int
    664 get_fields (s, delim, fields)
    665 char *s;
    666 int delim;
    667 char *fields[];
     302get_fields (char *s, int delim, char *fields[])
    668303{
    669304        int n;
    … …  
    684319}
    685320
    686 /* check name for being a planet.
    687  * if so, fill in *op and return 0, else return -1.
    688  */
    689 int
    690 db_chk_planet (name, op)
    691 char name[];
    692 Obj *op;
    693 {
    694         char namecpy[256];
    695         int i;
    696 
    697         /* these must match the order in astro.h */
    698         static char *planet_names[] = {
    699             "Mercury", "Venus", "Mars", "Jupiter", "Saturn",
    700             "Uranus", "Neptune", "Pluto", "Sun", "Moon",
    701         };
    702 
    703         /* make a copy to match our case -- strcasecmp() not entirely portable*/
    704         strcpy (namecpy, name);
    705         if (islower(namecpy[0]))
    706             namecpy[0] = toupper(namecpy[0]);
    707         for (i = 1; namecpy[i]; i++)
    708             if (isupper(namecpy[i]))
    709                 namecpy[i] = tolower(namecpy[i]);
    710 
    711         for (i = MERCURY; i <= MOON; i++) {
    712             if (strcmp (namecpy, planet_names[i]) == 0) {
    713                 zero_mem ((void *)op, sizeof(ObjPl));
    714                 op->o_type = PLANET;
    715                 (void) strcpy (op->o_name, planet_names[i]);
    716                 op->pl.pl_code = i;
    717                 return (0);
    718             }
    719         }
    720 
    721         return (-1);
    722 }
    723 
    724321/* return 0 if buf qualifies as a database line worthy of a cracking
    725322 * attempt, else -1.
    726323 */
     324int
     325dbline_candidate (char *buf)
     326{
     327        char c = buf[0];
     328
     329        return (c == '#' || c == '!' || isspace(c) ? -1 : 0);
     330}
     331
     332/* return 0 if TLE checksum is ok, else -1 */
    727333static int
    728 line_candidate (buf)
    729 char *buf;
    730 {
    731         char c = buf[0];
    732 
    733         return (c == '#' || c == '!' || isspace(c) ? -1 : 0);
    734 }
    735 
    736 /* given either a decimal year (xxxx[.xxx]) or a calendar (x/x/x)
    737  * and a DateFormat preference convert it to an mjd and store it at *p.
     334tle_sum (char *l)
     335{
     336        char *lastl = l + 68;
     337        int sum;
     338
     339        for (sum = 0; l < lastl; ) {
     340            char c = *l++;
     341            if (c == '\0')
     342                return (-1);
     343            if (isdigit(c))
     344                sum += c - '0';
     345            else if (c == '-')
     346                sum++;
     347        }
     348
     349        return (*l - '0' == (sum%10) ? 0 : -1);
     350}
     351
     352/* extract the given columns and return value.
     353 * N.B. from and to are 1-based within l
     354 */
     355static double
     356tle_fld (char *l, int from, int thru)
     357{
     358        char buf[32];
     359
     360        sprintf (buf, "%.*s", thru-from+1, l+from-1);
     361        return (atod (buf));
     362}
     363
     364/* extract the exponential value starting at the given column.
     365 * N.B. start is 1-based within l
     366 */
     367static double
     368tle_expfld (char *l, int start)
     369{
     370        char buf[32];
     371        double v;
     372
     373        sprintf (buf, ".%.*s", 5, l+start);
     374        v = atod (buf) * pow (10.0, tle_fld(l, start+6, start+7));
     375        if (l[start-1] == '-')
     376            v = -v;
     377        return (v);
     378}
     379
     380static int
     381crack_f (Obj *op, char *flds[MAXFLDS], int nf, char whynot[])
     382{
     383        char *sflds[MAXFLDS];
     384        double tmp;
     385        int nsf;
     386
     387        if (nf < 5 || nf > 7) {
     388            if (whynot)
     389                sprintf (whynot, "%s: f needs 5-7 fields, not %d",enm(flds),nf);
     390            return (-1);
     391        }
     392
     393        zero_mem ((void *)op, sizeof(ObjF));
     394        op->o_type = FIXED;
     395
     396        nsf = get_fields(flds[1], SUBFLD, sflds);
     397        if (nsf > 1) {
     398            switch (sflds[1][0]) {
     399            case 'A': case 'B': case 'C': case 'D': case 'F': case 'G':
     400            case 'H': case 'K': case 'J': case 'L': case 'M': case 'N':
     401            case 'O': case 'P': case 'Q': case 'R': case 'S': case 'T':
     402            case 'U': case 'V': case 'Y':
     403                op->f_class = sflds[1][0];
     404                if (op->f_class == 'B')
     405                    op->f_class = 'D';  /* merge B and D since BINARYSTAR */
     406                break;
     407            default:
     408                if (whynot)
     409                    sprintf (whynot, "%s: Bad f class: %c", enm(flds),
     410                                                                sflds[1][0]);
     411                return (-1);
     412            }
     413        }
     414        if (nsf > 2) {
     415            /* fill f_spect all the way */
     416            char buf[sizeof(op->f_spect)+1];
     417            memset (buf, 0, sizeof(buf));
     418            sprintf (buf, "%.*s", (int)sizeof(op->f_spect), sflds[2]);
     419            memcpy (op->f_spect, buf, (int)sizeof(op->f_spect));
     420        }
     421
     422        nsf = get_fields(flds[2], SUBFLD, sflds);
     423        f_scansexa (sflds[0], &tmp);
     424        op->f_RA = (float) hrrad(tmp);
     425        if (nsf > 1)
     426            op->f_pmRA = (float) 1.327e-11*atod(sflds[1]);/*mas/yr->rad/dy*/
     427
     428        nsf = get_fields(flds[3], SUBFLD, sflds);
     429        f_scansexa (sflds[0], &tmp);
     430        op->f_dec = (float) degrad(tmp);
     431        if (nsf > 1)
     432            op->f_pmdec = (float)1.327e-11*atod(sflds[1]);/*mas/yr->rad/dy*/
     433        if (fabs(op->f_dec) < PI/2)
     434            op->f_pmRA /= cos (op->f_dec);
     435
     436        set_fmag (op, atod(flds[4]));
     437
     438        if (nf > 5 && flds[5][0]) {
     439            tmp = op->f_epoch;
     440            crack_year (flds[5], &tmp);
     441            op->f_epoch = (float) tmp;
     442        } else
     443            op->f_epoch = J2000;        /* default */
     444
     445        if (nf > 6) {
     446            op->f_size = (float) atod(flds[6]);
     447
     448            /* optional minor axis and position angle subfields */
     449            nsf = get_fields(flds[6], SUBFLD, sflds);
     450            if (nsf == 3) {
     451                set_ratio(op, op->s_size, atod(sflds[1]));
     452                set_pa(op,degrad(atod(sflds[2])));
     453            } else {
     454                set_ratio(op,1,1);      /* round */
     455                set_pa(op,0.0);
     456            }
     457        }
     458
     459        return (0);
     460}
     461
     462static int
     463crack_e (Obj *op, char *flds[MAXFLDS], int nf, char whynot[])
     464{
     465        if (nf != 13 && nf != 14) {
     466            if (whynot)
     467                sprintf (whynot, "%s: e needs 13 or 14 fields, not %d",
     468                                                                enm(flds), nf);
     469            return (-1);
     470        }
     471
     472        zero_mem ((void *)op, sizeof(ObjE));
     473        op->o_type = ELLIPTICAL;
     474
     475        op->e_inc = (float) atod (flds[2]);
     476        op->e_Om = (float) atod (flds[3]);
     477        op->e_om = (float) atod (flds[4]);
     478        op->e_a = (float) atod (flds[5]);
     479        /* retired op->e_n = (float) atod (flds[6]); */
     480        op->e_e = atod (flds[7]);
     481        op->e_M = (float) atod (flds[8]);
     482        crack_year (flds[9], &op->e_cepoch);
     483        crack_okdates (flds[9], &op->e_startok, &op->e_endok);
     484        crack_year (flds[10], &op->e_epoch);
     485
     486        /* magnitude model gk or HG(default). allow prefixes in either field */
     487        op->e_mag.whichm = flds[11][0] == 'g' ? MAG_gk : MAG_HG;
     488        if (isdigit(flds[11][0]))
     489            op->e_mag.m1 = (float) atod(&flds[11][0]);
     490        else
     491            op->e_mag.m1 = (float) atod(&flds[11][1]);
     492        if (isdigit(flds[12][0]))
     493            op->e_mag.m2 = (float) atod(&flds[12][0]);
     494        else
     495            op->e_mag.m2 = (float) atod(&flds[12][1]);
     496
     497        if (nf == 14)
     498            op->e_size = (float) atod (flds[13]);
     499
     500        return (0);
     501}
     502
     503static int
     504crack_h (Obj *op, char *flds[MAXFLDS], int nf, char whynot[])
     505{
     506        if (nf != 11 && nf != 12) {
     507            if (whynot)
     508                sprintf (whynot, "%s: h needs 11 or 12 fields, not %d",
     509                                                                enm(flds), nf);
     510            return (-1);
     511        }
     512
     513        zero_mem ((void *)op, sizeof(ObjH));
     514        op->o_type = HYPERBOLIC;
     515
     516        crack_year (flds[2], &op->h_ep);
     517        crack_okdates (flds[2], &op->h_startok, &op->h_endok);
     518        op->h_inc = (float) atod (flds[3]);
     519        op->h_Om = (float) atod (flds[4]);
     520        op->h_om = (float) atod (flds[5]);
     521        op->h_e = (float) atod (flds[6]);
     522        op->h_qp = (float) atod (flds[7]);
     523        crack_year (flds[8], &op->h_epoch);
     524        op->h_g = (float) atod (flds[9]);
     525        op->h_k = (float) atod (flds[10]);
     526
     527        if (nf == 12)
     528            op->h_size = (float) atod (flds[11]);
     529
     530        return (0);
     531}
     532
     533static int
     534crack_p (Obj *op, char *flds[MAXFLDS], int nf, char whynot[])
     535{
     536        if (nf != 10 && nf != 11) {
     537            if (whynot)
     538                sprintf (whynot, "%s: p needs 10 or 11 fields, not %d",
     539                                                                enm(flds), nf);
     540            return (-1);
     541        }
     542
     543        zero_mem ((void *)op, sizeof(ObjP));
     544        op->o_type = PARABOLIC;
     545
     546        crack_year (flds[2], &op->p_ep);
     547        crack_okdates (flds[2], &op->p_startok, &op->p_endok);
     548        op->p_inc = (float) atod (flds[3]);
     549        op->p_om = (float) atod (flds[4]);
     550        op->p_qp = (float) atod (flds[5]);
     551        op->p_Om = (float) atod (flds[6]);
     552        crack_year (flds[7], &op->p_epoch);
     553        op->p_g = (float) atod (flds[8]);
     554        op->p_k = (float) atod (flds[9]);
     555
     556        if (nf == 11)
     557            op->p_size = (float) atod (flds[10]);
     558
     559        return (0);
     560}
     561
     562static int
     563crack_E (Obj *op, char *flds[MAXFLDS], int nf, char whynot[])
     564{
     565        if (nf != 11 && nf != 12) {
     566            if (whynot)
     567                sprintf (whynot, "%s: E needs 11 or 12 fields, not %d",
     568                                                            enm(flds), nf);
     569            return (-1);
     570        }
     571
     572        zero_mem ((void *)op, sizeof(ObjES));
     573        op->o_type = EARTHSAT;
     574        crack_year (flds[2], &op->es_epoch);
     575        crack_okdates (flds[2], &op->es_startok, &op->es_endok);
     576        op->es_inc = (float) atod (flds[3]);
     577        op->es_raan = (float) atod (flds[4]);
     578        op->es_e = (float) atod (flds[5]);
     579        op->es_ap = (float) atod (flds[6]);
     580        op->es_M = (float) atod (flds[7]);
     581        op->es_n = atod (flds[8]);
     582        op->es_decay = (float) atod (flds[9]);
     583        op->es_orbit = atoi (flds[10]);
     584        if (nf == 12)
     585            op->es_drag = (float) atod (flds[11]);
     586
     587        return (0);
     588}
     589
     590static int
     591crack_P (Obj *op, char *flds[MAXFLDS], int nf, char whynot[])
     592{
     593        Obj *bi;
     594        int nbi;
     595        int i;
     596
     597        nbi = getBuiltInObjs (&bi);
     598
     599        for (i = 0; i < nbi; i++) {
     600            Obj *bop = bi + i;
     601            if (is_type(bop,PLANETM) && !strcmp (flds[0], bop->o_name)) {
     602                memcpy ((void *)op, bop, sizeof(ObjPl));
     603                return (0);
     604            }
     605        }
     606
     607        if (whynot)
     608            sprintf (whynot, "%s: Unknown planet or moon", enm(flds));
     609        return (-1);
     610}
     611
     612static int
     613crack_B (Obj *op, char *flds[MAXFLDS], int nf, char whynot[])
     614{
     615        char *sflds[MAXFLDS];
     616        double tmp;
     617        int nsf;
     618
     619        if (nf != 7) {
     620            if (whynot)
     621                sprintf (whynot, "%s: B need 7 fields, not %d", enm(flds), nf);
     622            return (-1);
     623        }
     624
     625        zero_mem ((void *)op, sizeof(ObjB));
     626        op->o_type = BINARYSTAR;
     627
     628        nsf = get_fields(flds[1], SUBFLD, sflds);
     629        if (nsf > 1) {
     630            switch (sflds[1][0]) {
     631            case 'a': case 'c': case 'e': case 'x': case 'y': case 'o':
     632            case 's': case 't': case 'u': case 'v': case 'b': case 'd':
     633            case 'q': case 'r': case 'p': case 'U': case 'V': case 'Y':
     634                op->f_class = sflds[1][0];
     635                break;
     636            default:
     637                if (whynot)
     638                    sprintf (whynot, "%s: Bad B class: %c", enm(flds),
     639                                                                sflds[1][0]);
     640                return (-1);
     641            }
     642        }
     643        if (nsf > 2) {
     644            /* fill f_spect all the way */
     645            char buf[sizeof(op->f_spect)+1];
     646            memset (buf, 0, sizeof(buf));
     647            sprintf (buf, "%.*s", (int)sizeof(op->f_spect), sflds[2]);
     648            memcpy (op->f_spect, buf, (int)sizeof(op->f_spect));
     649        }
     650        if (nsf > 3) {
     651            /* fill b_2spect all the way */
     652            char buf[sizeof(op->b_2spect)+1];
     653            memset (buf, 0, sizeof(buf));
     654            sprintf (buf, "%.*s", (int)sizeof(op->b_2spect), sflds[3]);
     655            memcpy (op->b_2spect, buf, (int)sizeof(op->b_2spect));
     656        }
     657
     658        nsf = get_fields(flds[2], SUBFLD, sflds);
     659        f_scansexa (sflds[0], &tmp);
     660        op->f_RA = (float) hrrad(tmp);
     661        if (nsf > 1)
     662            op->f_pmRA = (float) 1.327e-11*atod(sflds[1]);/*mas/yr->rad/dy*/
     663
     664        nsf = get_fields(flds[3], SUBFLD, sflds);
     665        f_scansexa (sflds[0], &tmp);
     666        op->f_dec = (float) degrad(tmp);
     667        if (nsf > 1)
     668            op->f_pmdec = (float)1.327e-11*atod(sflds[1]);/*mas/yr->rad/dy*/
     669        if (fabs(op->f_dec) < PI/2)
     670            op->f_pmRA /= cos (op->f_dec);
     671
     672        nsf = get_fields(flds[4], SUBFLD, sflds);
     673        if (nsf > 0)
     674            set_fmag (op, atod(sflds[0]));
     675        if (nsf > 1)
     676            op->b_2mag = (short)floor((atod(sflds[1]))*MAGSCALE + 0.5);
     677
     678        if (flds[5][0]) {
     679            tmp = op->f_epoch;
     680            crack_year (flds[5], &tmp);
     681            op->f_epoch = (float) tmp;
     682        } else
     683            op->f_epoch = J2000;        /* default */
     684
     685        nsf = get_fields(flds[6], SUBFLD, sflds);
     686        if (nsf == 7) {
     687            int l;
     688            char c;
     689
     690            op->b_bo.bo_a = atod(sflds[0]);
     691            op->b_bo.bo_i = atod(sflds[1]);
     692            op->b_bo.bo_O = atod(sflds[2]);
     693            op->b_bo.bo_e = atod(sflds[3]);
     694            op->b_bo.bo_T = atod(sflds[4]);
     695            op->b_bo.bo_o = atod(sflds[5]);
     696            op->b_bo.bo_P = atod(sflds[6]);
     697
     698            /* reject some weird entries actually seen in real lists */
     699            if (op->b_bo.bo_a <= 0) {
     700                if (whynot)
     701                    sprintf (whynot, "%s: Bogus B semi major axis: %g",
     702                                                    enm(flds), op->b_bo.bo_a);
     703                return (-1);
     704            }
     705            if (op->b_bo.bo_P <= 0) {
     706                if (whynot)
     707                    sprintf (whynot, "%s: Bogus B period: %g", enm(flds),
     708                                                                op->b_bo.bo_P);
     709                return (-1);
     710            }
     711
     712            /* scale period */
     713            l = strlen (sflds[6]);
     714            c = sflds[6][l-1];
     715            switch (c) {
     716            case 'y': case 'Y':
     717                break;
     718            case 'h': case 'H':
     719                op->b_bo.bo_P /= (24.0*365.25);
     720                break;
     721            case 'd': case 'D':
     722                op->b_bo.bo_P /= 365.25;
     723                break;
     724            default:
     725                if (c != ' ' && !isdigit(c)) {
     726                    if (whynot)
     727                        sprintf (whynot,"%s: B period suffix not Y, D or H: %c",
     728                                                                enm(flds), c);
     729                    return (-1);
     730                }
     731            }
     732
     733        } else if (nsf==3 || nsf==6 || nsf==9) {
     734            double yr;
     735            int i;
     736
     737            op->b_nbp = nsf/3;
     738            for (i = 0; i < nsf; i += 3) {
     739                tmp = 0;
     740                crack_year (sflds[i+0], &tmp);
     741                mjd_year (tmp, &yr);
     742                op->b_bp[i/3].bp_ep = (float)yr;
     743                op->b_bp[i/3].bp_sep = atod(sflds[i+1]);
     744                op->b_bp[i/3].bp_pa = degrad(atod(sflds[i+2]));
     745            }
     746        } else {
     747            if (whynot)
     748                sprintf (whynot,
     749                        "%s: B needs 3,6 or 7 subfields in field 7, not %d",
     750                                                                enm(flds), nsf);
     751            return (-1);
     752        }
     753
     754        return (0);
     755}
     756
     757/* put all names in nm but load only the first into o_name */
     758static int
     759crack_name (Obj *op, char *flds[MAXFLDS], int nf, char nm[][MAXNM], int nnm)
     760{
     761        char *sflds[MAXFLDS];
     762        int nsf;
     763        int i;
     764       
     765        nsf = get_fields (flds[0], SUBFLD, sflds);
     766        for (i = 0; nm && i < nsf && i < nnm; i++) {
     767            strncpy (nm[i], sflds[i], MAXNM);
     768            nm[i][MAXNM-1] = '\0';
     769        }
     770        strncpy (op->o_name, sflds[0], MAXNM-1);
     771        return (nsf);
     772}
     773
     774/* simple name cracker just for error messages */
     775static char *
     776enm (char *flds[MAXFLDS])
     777{
     778        char *sflds[MAXFLDS];
     779        int nsf = get_fields (flds[0], SUBFLD, sflds);
     780        return (nsf > 0 ? sflds[0] : "Unknown");
     781}
     782
     783/* given either a decimal year (xxxx[.xxx]) or a calendar (x/x/x) date
     784 * convert it to an mjd and store it at *p.
    738785 */
    739786static void
    740 crack_year (bp, pref, p)
    741 char *bp;
    742 PrefDateFormat pref;
    743 double *p;
     787crack_year (char *bp, double *p)
    744788{
    745789        int m, y;
    … …  
    747791
    748792        mjd_cal (*p, &m, &d, &y);       /* init with current */
    749         f_sscandate (bp, pref, &m, &d, &y);
     793        f_sscandate (bp, PREF_MDY, &m, &d, &y);
    750794        cal_mjd (m, d, y, p);
    751795}
    752796
    753 /* given an *op and a field id, add it to the given text buffer bp.
    754  * return the number of characters added to bp.
     797/* crack the startok and endok date fields found in several Obj types.
     798 * set to 0 if blank or any problems.
     799 */
     800static void
     801crack_okdates (char *fld, float *startok, float *endok)
     802{
     803        char *sflds[MAXFLDS];
     804        double tmp;
     805        int m, y;
     806        double d;
     807        int nsf;
     808
     809        *startok = *endok = 0;
     810        nsf = get_fields(fld, SUBFLD, sflds);
     811        if (nsf > 1) {
     812            d = m = y = 0;
     813            f_sscandate (sflds[1], PREF_MDY, &m, &d, &y);
     814            cal_mjd (m, d, y, &tmp);
     815            *startok = (float)tmp;
     816            if (nsf > 2) {
     817                d = m = y = 0;
     818                f_sscandate (sflds[2], PREF_MDY, &m, &d, &y);
     819                cal_mjd (m, d, y, &tmp);
     820                *endok = (float)tmp;
     821            }
     822        }
     823}
     824
     825/* add startok and endok to string at lp if non-zero.
     826 * return number of characters added.
    755827 */
    756828static int
    757 db_get_field (op, id, bp)
    758 Obj *op;
    759 int id;
    760 char *bp;
    761 {
    762         char *bpsave = bp;
     829get_okdates (char *lp, float *sp, float *ep)
     830{
     831        char *lp0 = lp;
     832
     833        if (*sp || *ep) {
     834            *lp++ = '|';
     835            if (*sp)
     836                lp += fs_date (lp, *sp);
     837            if (*ep) {
     838                *lp++ = '|';
     839                lp += fs_date (lp, *ep);
     840            }
     841        }
     842
     843        return (lp - lp0);
     844}
     845
     846static void
     847write_f (Obj *op, char lp[])
     848{
    763849        double tmp;
    764850
    765         /* include all the enums and in numeric order to give us the best
    766          * possible chance the compiler will implement this as a jump table.
    767          */
    768         switch (id) {
    769         case F_RA:
    770             sprintf (bp, ",");
    771             bp += strlen(bp);
    772             fs_sexa (bp, radhr(op->f_RA), 2, 36000);
    773             bp += strlen(bp);
    774             break;
    775         case F_DEC:
    776             sprintf (bp, ",");
    777             bp += strlen(bp);
    778             fs_sexa (bp, raddeg(op->f_dec), 3, 3600);
    779             bp += strlen(bp);
    780             break;
    781         case F_EPOCH:
    782             mjd_year (op->f_epoch, &tmp);
    783             sprintf (bp, ",%.6g", tmp); /* %.7g gives 2000.001 */
    784             bp += strlen(bp);
    785             break;
    786         case F_MAG:
    787             sprintf (bp, ",%6.2f", get_mag(op));
    788             bp += strlen(bp);
    789             break;
    790         case F_SIZE:
    791             sprintf (bp, ",%.7g", op->f_size);
    792             bp += strlen(bp);
    793             if (op->f_ratio || op->f_pa) {
    794                 sprintf (bp,"|%g|%g", op->f_size*get_ratio(op),
     851        lp += sprintf (lp, "%s,f", op->o_name);
     852        if (op->f_class)
     853            lp += sprintf (lp, "|%c", op->f_class);
     854        if (op->f_spect[0])
     855            lp += sprintf (lp, "|%.*s", (int)sizeof(op->f_spect), op->f_spect);
     856        *lp++ = ',';
     857        lp += fs_sexa (lp, radhr(op->f_RA), 2, 36000);
     858        if (op->f_pmRA)
     859            lp += sprintf (lp, "|%.6g",cos(op->f_dec)*op->f_pmRA/1.327e-11);
     860        *lp++ = ',';
     861        lp += fs_sexa (lp, raddeg(op->f_dec), 3, 3600);
     862        if (op->f_pmdec)
     863            lp += sprintf (lp, "|%.6g", op->f_pmdec/1.327e-11);
     864        lp += sprintf (lp, ",%.2f", get_mag(op));
     865        mjd_year (op->f_epoch, &tmp);
     866        lp += sprintf (lp, ",%.6g", tmp); /* %.7g gives 2000.001 */
     867        lp += sprintf (lp, ",%.7g", op->f_size);
     868        if (op->f_size && (op->f_ratio || op->f_pa))
     869            lp += sprintf (lp,"|%g|%g", op->f_size*get_ratio(op),
    795870                                                            raddeg(get_pa(op)));
    796                 bp += strlen(bp);
    797             }
    798             break;
    799         case F_CLASS:
    800             if (op->f_class) {
    801                 sprintf (bp, "|%c", op->f_class);
    802                 bp += strlen(bp);
    803             }
    804             break;
    805         case F_SPECT:
    806             if (op->f_spect[0] != '\0') {
    807                 sprintf (bp, "|%c", op->f_spect[0]);
    808                 bp += strlen(bp);
    809                 if (op->f_spect[1] != '\0') {
    810                     sprintf (bp, "%c", op->f_spect[1]);
    811                     bp += strlen(bp);
    812                 }
    813             }
    814             break;
    815 
    816         case E_INC:
    817             sprintf (bp, ",%.7g", op->e_inc);
    818             bp += strlen(bp);
    819             break;
    820         case E_LAN:
    821             sprintf (bp, ",%.7g", op->e_Om);
    822             bp += strlen(bp);
    823             break;
    824         case E_AOP:
    825             sprintf (bp, ",%.7g", op->e_om);
    826             bp += strlen(bp);
    827             break;
    828         case E_A:
    829             sprintf (bp, ",%.7g", op->e_a);
    830             bp += strlen(bp);
    831             break;
    832         case E_N:
    833             /* retired */
    834             sprintf (bp, ",");
    835             bp += strlen(bp);
    836             break;
    837         case E_E:
    838             sprintf (bp, ",%.7g", op->e_e);
    839             bp += strlen(bp);
    840             break;
    841         case E_M:
    842             sprintf (bp, ",%.7g", op->e_M);
    843             bp += strlen(bp);
    844             break;
    845         case E_CEPOCH:
    846             sprintf (bp, ",");
    847             bp += strlen(bp);
    848             fs_date (bp, op->e_cepoch);
    849             bp += strlen(bp);
    850             break;
    851         case E_EPOCH:
    852             sprintf (bp, ",");
    853             bp += strlen(bp);
    854             fs_date (bp, op->e_epoch);
    855             bp += strlen(bp);
    856             break;
    857         case E_M1:
    858             if (op->e_mag.whichm == MAG_gk) {
    859                 sprintf (bp, ",g%.7g", op->e_mag.m1);
    860                 bp += strlen(bp);
    861             } else if (op->e_mag.whichm == MAG_HG) {
    862                 sprintf (bp, ",H%.7g", op->e_mag.m1);
    863                 bp += strlen(bp);
    864             } else {
    865                 sprintf (bp, ",%.7g", op->e_mag.m1);
    866                 bp += strlen(bp);
    867             }
    868             break;
    869         case E_M2:
    870             sprintf (bp, ",%.7g", op->e_mag.m2);
    871             bp += strlen(bp);
    872             break;
    873         case E_SIZE:
    874             sprintf (bp, ",%.7g", op->e_size);
    875             bp += strlen(bp);
    876             break;
    877 
    878         case H_EP:
    879             sprintf (bp, ",");
    880             bp += strlen(bp);
    881             fs_date (bp, op->h_ep);
    882             bp += strlen(bp);
    883             break;
    884         case H_INC:
    885             sprintf (bp, ",%.7g", op->h_inc);
    886             bp += strlen(bp);
    887             break;
    888         case H_LAN:
    889             sprintf (bp, ",%.7g", op->h_Om);
    890             bp += strlen(bp);
    891             break;
    892         case H_AOP:
    893             sprintf (bp, ",%.7g", op->h_om);
    894             bp += strlen(bp);
    895             break;
    896         case H_E:
    897             sprintf (bp, ",%.7g", op->h_e);
    898             bp += strlen(bp);
    899             break;
    900         case H_QP:
    901             sprintf (bp, ",%.7g", op->h_qp);
    902             bp += strlen(bp);
    903             break;
    904         case H_EPOCH:
    905             sprintf (bp, ",");
    906             bp += strlen(bp);
    907             fs_date (bp, op->h_epoch);
    908             bp += strlen(bp);
    909             break;
    910         case H_G:
    911             sprintf (bp, ",%.7g", op->h_g);
    912             bp += strlen(bp);
    913             break;
    914         case H_K:
    915             sprintf (bp, ",%.7g", op->h_k);
    916             bp += strlen(bp);
    917             break;
    918         case H_SIZE:
    919             sprintf (bp, ",%.7g", op->h_size);
    920             bp += strlen(bp);
    921             break;
    922 
    923         case P_EP:
    924             sprintf (bp, ",");
    925             bp += strlen(bp);
    926             fs_date (bp, op->p_ep);
    927             bp += strlen(bp);
    928             break;
    929         case P_INC:
    930             sprintf (bp, ",%.7g", op->p_inc);
    931             bp += strlen(bp);
    932             break;
    933         case P_AOP:
    934             sprintf (bp, ",%.7g", op->p_om);
    935             bp += strlen(bp);
    936             break;
    937         case P_QP:
    938             sprintf (bp, ",%.7g", op->p_qp);
    939             bp += strlen(bp);
    940             break;
    941         case P_LAN:
    942             sprintf (bp, ",%.7g", op->p_Om);
    943             bp += strlen(bp);
    944             break;
    945         case P_EPOCH:
    946             sprintf (bp, ",");
    947             bp += strlen(bp);
    948             fs_date (bp, op->p_epoch);
    949             bp += strlen(bp);
    950             break;
    951         case P_G:
    952             sprintf (bp, ",%.7g", op->p_g);
    953             bp += strlen(bp);
    954             break;
    955         case P_K:
    956             sprintf (bp, ",%.7g", op->p_k);
    957             bp += strlen(bp);
    958             break;
    959         case P_SIZE:
    960             sprintf (bp, ",%.7g", op->p_size);
    961             bp += strlen(bp);
    962             break;
    963 
    964         case ES_EPOCH:
    965             sprintf (bp, ",");
    966             bp += strlen(bp);
    967             fs_date (bp, op->es_epoch);
    968             bp += strlen(bp);
    969             break;
    970         case ES_INC:
    971             sprintf (bp, ",%.7g", op->es_inc);
    972             bp += strlen(bp);
    973             break;
    974         case ES_RAAN:
    975             sprintf (bp, ",%.7g", op->es_raan);
    976             bp += strlen(bp);
    977             break;
    978         case ES_E:
    979             sprintf (bp, ",%.7g", op->es_e);
    980             bp += strlen(bp);
    981             break;
    982         case ES_AP:
    983             sprintf (bp, ",%.7g", op->es_ap);
    984             bp += strlen(bp);
    985             break;
    986         case ES_M:
    987             sprintf (bp, ",%.7g", op->es_M);
    988             bp += strlen(bp);
    989             break;
    990         case ES_N:
    991             sprintf (bp, ",%.7g", op->es_n);
    992             bp += strlen(bp);
    993             break;
    994         case ES_DECAY:
    995             sprintf (bp, ",%.7g", op->es_decay);
    996             bp += strlen(bp);
    997             break;
    998         case ES_ORBIT:
    999             sprintf (bp, ",%d", op->es_orbit);
    1000             bp += strlen(bp);
    1001             break;
    1002         case ES_DRAG:
    1003             sprintf (bp, ",%.7g", op->es_drag);
    1004             bp += strlen(bp);
    1005             break;
    1006 
    1007         default:
    1008             printf ("BUG! db_set_field: bad id: %d\n", id);
    1009             exit (1);
    1010         }
    1011 
    1012         return (bp - bpsave);
     871}
     872
     873static void
     874write_e (Obj *op, char lp[])
     875{
     876        lp += sprintf (lp, "%s,e", op->o_name);
     877        lp += sprintf (lp, ",%.7g", op->e_inc);
     878        lp += sprintf (lp, ",%.7g", op->e_Om);
     879        lp += sprintf (lp, ",%.7g", op->e_om);
     880        lp += sprintf (lp, ",%.7g", op->e_a);
     881        lp += sprintf (lp, ",%.7g", 0.0);               /* retired op->e_n */
     882        lp += sprintf (lp, ",%.7g", op->e_e);
     883        lp += sprintf (lp, ",%.7g", op->e_M);
     884        *lp++ = ',';
     885        lp += fs_date (lp, op->e_cepoch);
     886        lp += get_okdates (lp, &op->e_startok, &op->e_endok);
     887        *lp++ = ',';
     888        lp += fs_date (lp, op->e_epoch);
     889        if (op->e_mag.whichm == MAG_gk)
     890            lp += sprintf (lp, ",g%.7g", op->e_mag.m1);
     891        else if (op->e_mag.whichm == MAG_HG)
     892            lp += sprintf (lp, ",H%.7g", op->e_mag.m1);
     893        else
     894            lp += sprintf (lp, ",%.7g", op->e_mag.m1);
     895        lp += sprintf (lp, ",%.7g", op->e_mag.m2);
     896        lp += sprintf (lp, ",%.7g", op->e_size);
     897}
     898
     899static void
     900write_h (Obj *op, char lp[])
     901{
     902        lp += sprintf (lp, "%s,h", op->o_name);
     903        *lp++ = ',';
     904        lp += fs_date (lp, op->h_ep);
     905        lp += get_okdates (lp, &op->h_startok, &op->h_endok);
     906        lp += sprintf (lp, ",%.7g", op->h_inc);
     907        lp += sprintf (lp, ",%.7g", op->h_Om);
     908        lp += sprintf (lp, ",%.7g", op->h_om);
     909        lp += sprintf (lp, ",%.7g", op->h_e);
     910        lp += sprintf (lp, ",%.7g", op->h_qp);
     911        *lp++ = ',';
     912        lp += fs_date (lp, op->h_epoch);
     913        lp += sprintf (lp, ",%.7g", op->h_g);
     914        lp += sprintf (lp, ",%.7g", op->h_k);
     915        lp += sprintf (lp, ",%.7g", op->h_size);
     916}
     917
     918static void
     919write_p (Obj *op, char lp[])
     920{
     921        lp += sprintf (lp, "%s,p", op->o_name);
     922        *lp++ = ',';
     923        lp += fs_date (lp, op->p_ep);
     924        lp += get_okdates (lp, &op->p_startok, &op->p_endok);
     925        lp += sprintf (lp, ",%.7g", op->p_inc);
     926        lp += sprintf (lp, ",%.7g", op->p_om);
     927        lp += sprintf (lp, ",%.7g", op->p_qp);
     928        lp += sprintf (lp, ",%.7g", op->p_Om);
     929        *lp++ = ',';
     930        lp += fs_date (lp, op->p_epoch);
     931        lp += sprintf (lp, ",%.7g", op->p_g);
     932        lp += sprintf (lp, ",%.7g", op->p_k);
     933        lp += sprintf (lp, ",%.7g", op->p_size);
     934}
     935
     936static void
     937write_E (Obj *op, char lp[])
     938{
     939        lp += sprintf (lp, "%s,E", op->o_name);
     940        *lp++ = ',';
     941        lp += fs_date (lp, op->es_epoch);
     942        lp += get_okdates (lp, &op->es_startok, &op->es_endok);
     943        lp += sprintf (lp, ",%.7g", op->es_inc);
     944        lp += sprintf (lp, ",%.7g", op->es_raan);
     945        lp += sprintf (lp, ",%.7g", op->es_e);
     946        lp += sprintf (lp, ",%.7g", op->es_ap);
     947        lp += sprintf (lp, ",%.7g", op->es_M);
     948        lp += sprintf (lp, ",%.7g", op->es_n);
     949        lp += sprintf (lp, ",%.7g", op->es_decay);
     950        lp += sprintf (lp, ",%d", op->es_orbit);
     951        lp += sprintf (lp, ",%.7g", op->es_drag);
     952}
     953
     954static void
     955write_B (Obj *op, char lp[])
     956{
     957        double tmp;
     958
     959        lp += sprintf (lp, "%s,B", op->o_name);
     960        if (op->f_class)
     961            lp += sprintf (lp, "|%c", op->f_class);
     962        if (op->f_spect[0])
     963            lp += sprintf (lp, "|%.*s", (int)sizeof(op->f_spect), op->f_spect);
     964        if (op->b_2spect[0])
     965            lp += sprintf (lp, "|%.*s", (int)sizeof(op->b_2spect),op->b_2spect);
     966        *lp++ = ',';
     967        lp += fs_sexa (lp, radhr(op->f_RA), 2, 36000);
     968        if (op->f_pmRA)
     969            lp += sprintf (lp, "|%.6g",cos(op->f_dec)*op->f_pmRA/1.327e-11);
     970        *lp++ = ',';
     971        lp += fs_sexa (lp, raddeg(op->f_dec), 3, 3600);
     972        if (op->f_pmdec)
     973            lp += sprintf (lp, "|%.6g", op->f_pmdec/1.327e-11);
     974        lp += sprintf (lp, ",%.2f", get_mag(op));
     975        lp += sprintf (lp, "|%.2f", op->b_2mag/MAGSCALE);
     976        mjd_year (op->f_epoch, &tmp);
     977        lp += sprintf (lp, ",%.6g", tmp); /* %.7g gives 2000.001 */
     978        if (op->b_nbp == 0) {
     979            lp += sprintf (lp, ",%.6g",  op->b_bo.bo_a);
     980            lp += sprintf (lp, "|%.6g",  op->b_bo.bo_i);
     981            lp += sprintf (lp, "|%.6g",  op->b_bo.bo_O);
     982            lp += sprintf (lp, "|%.6g",  op->b_bo.bo_e);
     983            lp += sprintf (lp, "|%.6g",  op->b_bo.bo_T);
     984            lp += sprintf (lp, "|%.6g",  op->b_bo.bo_o);
     985            lp += sprintf (lp, "|%.6gy", op->b_bo.bo_P);
     986        } else {
     987            int i;
     988
     989            for (i = 0; i < op->b_nbp; i++) {
     990                BinPos *bp = &op->b_bp[i];
     991                lp += sprintf (lp, "%c%.6g", i==0?',':'|', bp->bp_ep);
     992                lp += sprintf (lp, "|%.6g", bp->bp_sep);
     993                lp += sprintf (lp, "|%.6g", raddeg(bp->bp_pa));
     994            }
     995        }
     996}
     997
     998static void
     999write_P (Obj *op, char lp[])
     1000{
     1001
     1002        lp += sprintf (lp, "%s,P", op->o_name);
    10131003}
    10141004
    10151005/* For RCS Only -- Do Not Edit */
    1016 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: dbfmt.c,v $ $Date: 2001-10-22 12:08:26 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     1006static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: dbfmt.c,v $ $Date: 2004-06-15 16:52:38 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/deep.c

    r1719 r2551  
    788788
    789789/* For RCS Only -- Do Not Edit */
    790 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: deep.c,v $ $Date: 2001-10-22 12:08:26 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     790static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: deep.c,v $ $Date: 2004-06-15 16:52:38 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/deepconst.h

    r1719 r2551  
    22#define _CONST_H
    33
    4 /* $Id: deepconst.h,v 1.2 2001-10-22 12:08:26 cmv Exp $ */
     4/* $Id: deepconst.h,v 1.3 2004-06-15 16:52:38 cmv Exp $ */
    55
    66
    … …  
    3131
    3232/* For RCS Only -- Do Not Edit
    33  * @(#) $RCSfile: deepconst.h,v $ $Date: 2001-10-22 12:08:26 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $
     33 * @(#) $RCSfile: deepconst.h,v $ $Date: 2004-06-15 16:52:38 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $
    3434 */

  • trunk/SophyaExt/XephemAstroLib/deltat.c

    r1719 r2551  
    2929 * in 2130.
    3030 *
    31  * Input is mjd (modified julian date from MJD0 on). [stern]
    32  *  Note that xephem uses a different epoch for this "mjd" than the
     31 * Input is mj (modified julian date from MJD0 on). [stern]
     32 *  Note that xephem uses a different epoch for this "mj" than the
    3333 *  normal value of JD=240000.5.
    3434 * See AA page B4.
    3535 *
    36  * Output double deltat(mjd) is ET-UT1 in seconds.
     36 * Output double deltat(mj) is ET-UT1 in seconds.
    3737 *
    3838 *
    … …  
    6767 *   - adopted #include's for xephem
    6868 *   - made dt[] static
    69  *   - made mjd the time argument [was: year Y].
     69 *   - made mj the time argument [was: year Y].
    7070 *   - updated observed and extrapolated data from tables at
    7171 *      ftp://maia.usno.navy.mil/ser7/ -- data deviated by up to 0.8 s
    … …  
    7474 *   - replaced treatment after TABEND by linear extrapolation instead
    7575 *      of second order version
    76  *   - installed lastmjd cache (made ans static)
     76 *   - installed lastmj cache (made ans static)
    7777 *
    7878 *   - no changes to table interpolation scheme and past extrapolations */
    7979
    80 #include "P_.h"
     80#include <math.h>
     81
    8182#include "astro.h"
    8283
    … …  
    161162 * of the Earth rotation rate in the ET time scale.
    162163 */
    163 double deltat(mjd)
    164 double mjd;
     164double deltat(double mj)
    165165{
    166166        double Y;
    … …  
    168168        int d[6];
    169169        int i, iy, k;
    170         double floor();
    171170        static double ans;
    172         static double lastmjd = -10000;
    173 
    174         if (mjd == lastmjd) {
     171        static double lastmj = -10000;
     172
     173        if (mj == lastmj) {
    175174            return(ans);
    176175        }
    177         lastmjd = mjd;
    178 
    179         Y = 2000.0 + (mjd - J2000)/365.25;
     176        lastmj = mj;
     177
     178        Y = 2000.0 + (mj - J2000)/365.25;
    180179
    181180        if( Y > TABEND  &&  Y < 2130.0 ) {
    … …  
    305304
    306305/* For RCS Only -- Do Not Edit */
    307 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: deltat.c,v $ $Date: 2001-10-22 12:08:26 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     306static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: deltat.c,v $ $Date: 2004-06-15 16:52:38 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/earthsat.c

    r1719 r2551  
    2323#include <math.h>
    2424#include <string.h>
    25 
    26 #if defined(__STDC__)
    2725#include <stdlib.h>
    28 #endif
    29 
    30 #include "P_.h"
     26
    3127#include "astro.h"
    32 #include "circum.h"
    3328#include "preferences.h"
    3429
    … …  
    4237typedef double MAT3x3[3][3];
    4338
    44 static void esat_prop P_((Now *np, Obj *op, double *SatX, double *SatY, double
    45     *SatZ, double *SatVX, double *SatVY, double *SatVZ));
    46 static void GetSatelliteParams P_((Obj *op));
    47 static void GetSiteParams P_((Now *np));
    48 static double Kepler P_((double MeanAnomaly, double Eccentricity));
    49 static void GetSubSatPoint P_((double SatX, double SatY, double SatZ,
    50     double T, double *Latitude, double *Longitude, double *Height));
    51 static void GetSatPosition P_((double EpochTime, double EpochRAAN,
     39static int crazyOp (Now *np, Obj *op);
     40static void esat_prop (Now *np, Obj *op, double *SatX, double *SatY, double
     41    *SatZ, double *SatVX, double *SatVY, double *SatVZ);
     42static void GetSatelliteParams (Obj *op);
     43static void GetSiteParams (Now *np);
     44static double Kepler (double MeanAnomaly, double Eccentricity);
     45static void GetSubSatPoint (double SatX, double SatY, double SatZ,
     46    double T, double *Latitude, double *Longitude, double *Height);
     47static void GetSatPosition (double EpochTime, double EpochRAAN,
    5248    double EpochArgPerigee, double SemiMajorAxis, double Inclination,
    5349    double Eccentricity, double RAANPrecession, double PerigeePrecession,
    5450    double T, double TrueAnomaly, double *X, double *Y, double *Z,
    55     double *Radius, double *VX, double *VY, double *VZ));
    56 static void GetSitPosition P_((double SiteLat, double SiteLong,
     51    double *Radius, double *VX, double *VY, double *VZ);
     52static void GetSitPosition (double SiteLat, double SiteLong,
    5753    double SiteElevation, double CrntTime, double *SiteX, double *SiteY,
    58     double *SiteZ, double *SiteVX, double *SiteVY, MAT3x3 SiteMatrix));
    59 static void GetRange P_((double SiteX, double SiteY, double SiteZ,
     54    double *SiteZ, double *SiteVX, double *SiteVY, MAT3x3 SiteMatrix);
     55static void GetRange (double SiteX, double SiteY, double SiteZ,
    6056    double SiteVX, double SiteVY, double SatX, double SatY, double SatZ,
    6157    double SatVX, double SatVY, double SatVZ, double *Range,
    62     double *RangeRate));
    63 static void GetTopocentric P_((double SatX, double SatY, double SatZ,
     58    double *RangeRate);
     59static void GetTopocentric (double SatX, double SatY, double SatZ,
    6460    double SiteX, double SiteY, double SiteZ, MAT3x3 SiteMatrix, double *X,
    65     double *Y, double *Z));
    66 static void GetBearings P_((double SatX, double SatY, double SatZ,
     61    double *Y, double *Z);
     62static void GetBearings (double SatX, double SatY, double SatZ,
    6763    double SiteX, double SiteY, double SiteZ, MAT3x3 SiteMatrix,
    68     double *Azimuth, double *Elevation));
    69 static int Eclipsed P_((double SatX, double SatY, double SatZ,
    70     double SatRadius, double CrntTime));
    71 static void InitOrbitRoutines P_((double EpochDay, int AtEod));
     64    double *Azimuth, double *Elevation);
     65static int Eclipsed (double SatX, double SatY, double SatZ,
     66    double SatRadius, double CrntTime);
     67static void InitOrbitRoutines (double EpochDay, int AtEod);
    7268
    7369#ifdef USE_ORBIT_PROPAGATOR
    74 static void GetPrecession P_((double SemiMajorAxis, double Eccentricity,
    75     double Inclination, double *RAANPrecession, double *PerigeePrecession));
     70static void GetPrecession (double SemiMajorAxis, double Eccentricity,
     71    double Inclination, double *RAANPrecession, double *PerigeePrecession);
    7672#endif /* USE_ORBIT_PROPAGATOR */
    7773
    … …  
    135131 */
    136132int
    137 obj_earthsat (np, op)
    138 Now *np;
    139 Obj *op;
     133obj_earthsat (Now *np, Obj *op)
    140134{
    141135        double Radius;              /* From geocenter                  */
    … …  
    267261 */
    268262static void
    269 esat_prop (np, op, SatX, SatY, SatZ, SatVX, SatVY, SatVZ)
    270 Now *np;
    271 Obj *op;
    272 double *SatX,*SatY,*SatZ;
    273 double *SatVX,*SatVY,*SatVZ;
     263esat_prop (Now *np, Obj *op, double *SatX, double *SatY, double *SatZ,
     264double *SatVX, double *SatVY, double *SatVZ)
    274265{
    275266#ifdef USE_ORBIT_PROPAGATOR
    … …  
    284275        double Radius;
    285276        double CrntTime;
     277
     278        if (crazyOp (np, op)) {
     279            *SatX = *SatY = *SatZ = *SatVX = *SatVY = *SatVZ = 0;
     280            return;
     281        }
    286282
    287283        SemiMajorAxis = 331.25 * exp(2*log(MinutesPerDay/epochMeanMotion)/3);
    … …  
    331327        double dt;
    332328        int yr;
     329
     330        if (crazyOp (np, op)) {
     331            *SatX = *SatY = *SatZ = *SatVX = *SatVY = *SatVZ = 0;
     332            return;
     333        }
    333334
    334335        /* init */
    … …  
    392393}
    393394
     395/* return 1 if op is crazy @ np */
     396static int
     397crazyOp (Now *np, Obj *op)
     398{
     399        /* figure its daft if long enough to decay 1 rev/day */
     400        return (fabs(op->es_epoch - mjd) > 1/op->es_decay);
     401}
    394402
    395403/* grab the xephem stuff from op and copy into orbit's globals.
    396404 */
    397405static void
    398 GetSatelliteParams(op)
    399 Obj *op;
     406GetSatelliteParams(Obj *op)
    400407{
    401408        /* the following are for the orbit functions */
    … …  
    429436 
    430437static void
    431 GetSiteParams(np)
    432 Now *np;
     438GetSiteParams(Now *np)
    433439{
    434440        SiteLat = lat;
    … …  
    465471 
    466472static
    467 double Kepler(MeanAnomaly,Eccentricity)
    468 register double MeanAnomaly,Eccentricity;
    469  
     473double Kepler(double MeanAnomaly, double Eccentricity)
    470474{
    471475register double E;              /* Eccentric Anomaly                    */
    … …  
    494498 
    495499static void
    496 GetSubSatPoint(SatX,SatY,SatZ,T,Latitude,Longitude,Height)
    497 double SatX,SatY,SatZ,T;
    498 double *Latitude,*Longitude,*Height;
     500GetSubSatPoint(double SatX, double SatY, double SatZ, double T,
     501double *Latitude, double *Longitude, double *Height)
    499502{
    500503    double r;
    … …  
    524527#ifdef USE_ORBIT_PROPAGATOR
    525528static void
    526 GetPrecession(SemiMajorAxis,Eccentricity,Inclination,
    527         RAANPrecession,PerigeePrecession)
    528 double SemiMajorAxis,Eccentricity,Inclination;
    529 double *RAANPrecession,*PerigeePrecession;
     529GetPrecession(double SemiMajorAxis, double Eccentricity, double Inclination,
     530double *RAANPrecession, double *PerigeePrecession)
    530531{
    531532  *RAANPrecession = 9.95*pow(EarthRadius/SemiMajorAxis,3.5) * cos(Inclination)
    … …  
    544545
    545546static void
    546 GetSatPosition(EpochTime,EpochRAAN,EpochArgPerigee,SemiMajorAxis,
    547         Inclination,Eccentricity,RAANPrecession,PerigeePrecession,
    548         T,TrueAnomaly,X,Y,Z,Radius,VX,VY,VZ)
    549  
    550 double EpochTime,EpochRAAN,EpochArgPerigee;
    551 double SemiMajorAxis,Inclination,Eccentricity;
    552 double RAANPrecession,PerigeePrecession,T, TrueAnomaly;
    553 double *X,*Y,*Z,*Radius,*VX,*VY,*VZ;
     547GetSatPosition(double EpochTime, double EpochRAAN, double EpochArgPerigee,
     548double SemiMajorAxis, double Inclination, double Eccentricity,
     549double RAANPrecession, double PerigeePrecession, double T,
     550double TrueAnomaly, double *X, double *Y, double *Z, double *Radius,
     551double *VX, double *VY, double *VZ)
    554552
    555553{
    … …  
    607605
    608606static void
    609 GetSitPosition(SiteLat,SiteLong,SiteElevation,CrntTime,
    610              SiteX,SiteY,SiteZ,SiteVX,SiteVY,SiteMatrix)
    611 
    612 double SiteLat,SiteLong,SiteElevation,CrntTime;
    613 double *SiteX,*SiteY,*SiteZ,*SiteVX,*SiteVY;
    614 MAT3x3 SiteMatrix;
    615 
     607GetSitPosition(double SiteLat, double SiteLong, double SiteElevation,
     608double CrntTime, double *SiteX, double *SiteY, double *SiteZ, double *SiteVX,
     609double *SiteVY, MAT3x3 SiteMatrix)
    616610{
    617611    static double G1,G2; /* Used to correct for flattening of the Earth */
    … …  
    663657
    664658static void
    665 GetRange(SiteX,SiteY,SiteZ,SiteVX,SiteVY,
    666         SatX,SatY,SatZ,SatVX,SatVY,SatVZ,Range,RangeRate)
    667 
    668 double SiteX,SiteY,SiteZ,SiteVX,SiteVY;
    669 double SatX,SatY,SatZ,SatVX,SatVY,SatVZ;
    670 double *Range,*RangeRate;
     659GetRange(double SiteX, double SiteY, double SiteZ, double SiteVX,
     660double SiteVY, double SatX, double SatY, double SatZ, double SatVX,
     661double SatVY, double SatVZ, double *Range, double *RangeRate)
    671662{
    672663    double DX,DY,DZ;
    … …  
    684675
    685676static void
    686 GetTopocentric(SatX,SatY,SatZ,SiteX,SiteY,SiteZ,SiteMatrix,X,Y,Z)
    687 double SatX,SatY,SatZ,SiteX,SiteY,SiteZ;
    688 double *X,*Y,*Z;
    689 MAT3x3 SiteMatrix;
     677GetTopocentric(double SatX, double SatY, double SatZ, double SiteX,
     678double SiteY, double SiteZ, MAT3x3 SiteMatrix, double *X, double *Y,
     679double *Z)
    690680{
    691681    SatX -= SiteX;
    … …  
    702692
    703693static void
    704 GetBearings(SatX,SatY,SatZ,SiteX,SiteY,SiteZ,SiteMatrix,Azimuth,Elevation)
    705 double SatX,SatY,SatZ,SiteX,SiteY,SiteZ;
    706 MAT3x3 SiteMatrix;
    707 double *Azimuth,*Elevation;
     694GetBearings(double SatX, double SatY, double SatZ, double SiteX,
     695double SiteY, double SiteZ, MAT3x3 SiteMatrix, double *Azimuth,
     696double *Elevation)
    708697{
    709698    double x,y,z;
    … …  
    720709
    721710static int
    722 Eclipsed(SatX,SatY,SatZ,SatRadius,CrntTime)
    723 double SatX,SatY,SatZ,SatRadius,CrntTime;
     711Eclipsed(double SatX, double SatY, double SatZ, double SatRadius,
     712double CrntTime)
    724713{
    725714    double MeanAnomaly,TrueAnomaly;
    … …  
    751740
    752741static void
    753 InitOrbitRoutines(EpochDay, AtEod)
    754 double EpochDay;
    755 int AtEod;
     742InitOrbitRoutines(double EpochDay, int AtEod)
    756743{
    757744    double T,T2,T3,Omega;
    … …  
    797784
    798785/* For RCS Only -- Do Not Edit */
    799 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: earthsat.c,v $ $Date: 2001-10-22 12:08:27 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     786static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: earthsat.c,v $ $Date: 2004-06-15 16:52:38 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/eq_ecl.c

    r1719 r2551  
    22#include <math.h>
    33
    4 #include "P_.h"
    54#include "astro.h"
    65
    7 static void ecleq_aux P_((int sw, double mjd, double x, double y,
    8     double *p, double *q));
     6static void ecleq_aux (int sw, double mj, double x, double y,
     7    double *p, double *q);
    98
    109#define EQtoECL 1
    … …  
    1211
    1312
    14 /* given the modified Julian date, mjd, and an equitorial ra and dec, each in
    15  * radians, find the corresponding geocentric ecliptic latitude, *lat, and
    16  * longititude, *lng, also each in radians.
     13/* given the modified Julian date, mj, and an equitorial ra and dec, each in
     14 * radians, find the corresponding geocentric ecliptic latitude, *lt, and
     15 * longititude, *lg, also each in radians.
    1716 * correction for the effect on the angle of the obliquity due to nutation is
    1817 * not included.
    1918 */
    2019void
    21 eq_ecl (mjd, ra, dec, lat, lng)
    22 double mjd;
    23 double ra, dec;
    24 double *lat, *lng;
     20eq_ecl (double mj, double ra, double dec, double *lt, double *lg)
    2521{
    26         ecleq_aux (EQtoECL, mjd, ra, dec, lng, lat);
     22        ecleq_aux (EQtoECL, mj, ra, dec, lg, lt);
    2723}
    2824
    29 /* given the modified Julian date, mjd, and a geocentric ecliptic latitude,
    30  * *lat, and longititude, *lng, each in radians, find the corresponding
     25/* given the modified Julian date, mj, and a geocentric ecliptic latitude,
     26 * *lt, and longititude, *lg, each in radians, find the corresponding
    3127 * equitorial ra and dec, also each in radians.
    3228 * correction for the effect on the angle of the obliquity due to nutation is
    … …  
    3430 */
    3531void
    36 ecl_eq (mjd, lat, lng, ra, dec)
    37 double mjd;
    38 double lat, lng;
    39 double *ra, *dec;
     32ecl_eq (double mj, double lt, double lg, double *ra, double *dec)
    4033{
    41         ecleq_aux (ECLtoEQ, mjd, lng, lat, ra, dec);
     34        ecleq_aux (ECLtoEQ, mj, lg, lt, ra, dec);
    4235}
    4336
    4437static void
    45 ecleq_aux (sw, mjd, x, y, p, q)
    46 int sw;                 /* +1 for eq to ecliptic, -1 for vv. */
    47 double mjd;
    48 double x, y;            /* sw==1: x==ra, y==dec.  sw==-1: x==lng, y==lat. */
    49 double *p, *q;          /* sw==1: p==lng, q==lat. sw==-1: p==ra, q==dec. */
     38ecleq_aux (
     39int sw,                 /* +1 for eq to ecliptic, -1 for vv. */
     40double mj,
     41double x, double y,     /* sw==1: x==ra, y==dec.  sw==-1: x==lg, y==lt. */
     42double *p, double *q)   /* sw==1: p==lg, q==lt. sw==-1: p==ra, q==dec. */
    5043{
    51         static double lastmjd = -10000; /* last mjd calculated */
     44        static double lastmj = -10000;  /* last mj calculated */
    5245        static double seps, ceps;       /* sin and cos of mean obliquity */
    5346        double sx, cx, sy, cy, ty, sq;
    5447
    55         if (mjd != lastmjd) {
     48        if (mj != lastmj) {
    5649            double eps;
    57             obliquity (mjd, &eps);              /* mean obliquity for date */
     50            obliquity (mj, &eps);               /* mean obliquity for date */
    5851            seps = sin(eps);
    5952            ceps = cos(eps);
    60             lastmjd = mjd;
     53            lastmj = mj;
    6154        }
    6255
    … …  
    7770
    7871/* For RCS Only -- Do Not Edit */
    79 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: eq_ecl.c,v $ $Date: 2001-10-22 12:08:27 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     72static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: eq_ecl.c,v $ $Date: 2004-06-15 16:52:38 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/eq_gal.c

    r1719 r2551  
    44#include <math.h>
    55
    6 #include "P_.h"
    76#include "astro.h"
    87
    9 static void galeq_aux P_((int sw, double x, double y, double *p, double *q));
    10 static void galeq_init P_((void));
     8static void galeq_aux (int sw, double x, double y, double *p, double *q);
     9static void galeq_init (void);
    1110
    1211#define EQtoGAL 1
    … …  
    1817static double gpd = degrad(27.12825);   /* Dec of  " */
    1918static double cgpd, sgpd;               /* cos() and sin() of gpd */
    20 static double mjd2000;                  /* mjd of 2000 */
     19static double mj2000;                   /* mj of 2000 */
    2120static int before;                      /* whether these have been set yet */
    2221
    2322/* given ra and dec, each in radians, for the given epoch, find the
    24  * corresponding galactic latitude, *lat, and longititude, *lng, also each in
     23 * corresponding galactic latitude, *lt, and longititude, *lg, also each in
    2524 * radians.
    2625 */
    2726void
    28 eq_gal (mjd, ra, dec, lat, lng)
    29 double mjd, ra, dec;
    30 double *lat, *lng;
     27eq_gal (double mj, double ra, double dec, double *lt, double *lg)
    3128{
    3229        galeq_init();
    33         precess (mjd, mjd2000, &ra, &dec);
    34         galeq_aux (EQtoGAL, ra, dec, lng, lat);
     30        precess (mj, mj2000, &ra, &dec);
     31        galeq_aux (EQtoGAL, ra, dec, lg, lt);
    3532}
    3633
    37 /* given galactic latitude, lat, and longititude, lng, each in radians, find
     34/* given galactic latitude, lt, and longititude, lg, each in radians, find
    3835 * the corresponding equitorial ra and dec, also each in radians, at the
    3936 * given epoch.
    4037 */
    4138void
    42 gal_eq (mjd, lat, lng, ra, dec)
    43 double mjd, lat, lng;
    44 double *ra, *dec;
     39gal_eq (double mj, double lt, double lg, double *ra, double *dec)
    4540{
    4641        galeq_init();
    47         galeq_aux (GALtoEQ, lng, lat, ra, dec);
    48         precess (mjd2000, mjd, ra, dec);
     42        galeq_aux (GALtoEQ, lg, lt, ra, dec);
     43        precess (mj2000, mj, ra, dec);
    4944}
    5045
    5146static void
    52 galeq_aux (sw, x, y, p, q)
    53 int sw;                 /* +1 for eq to gal, -1 for vv. */
    54 double x, y;            /* sw==1: x==ra, y==dec.  sw==-1: x==lng, y==lat. */
    55 double *p, *q;          /* sw==1: p==lng, q==lat. sw==-1: p==ra, q==dec. */
     47galeq_aux (
     48int sw,                 /* +1 for eq to gal, -1 for vv. */
     49double x, double y,     /* sw==1: x==ra, y==dec.  sw==-1: x==lg, y==lt. */
     50double *p, double *q)   /* sw==1: p==lg, q==lt. sw==-1: p==ra, q==dec. */
    5651{
    5752        double sy, cy, a, ca, sa, b, sq, c, d;
    … …  
    9691            cgpd = cos (gpd);
    9792            sgpd = sin (gpd);
    98             mjd2000 = J2000;
     93            mj2000 = J2000;
    9994            before = 1;
    10095        }
    … …  
    10297
    10398/* For RCS Only -- Do Not Edit */
    104 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: eq_gal.c,v $ $Date: 2001-10-22 12:08:27 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     99static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: eq_gal.c,v $ $Date: 2004-06-15 16:52:38 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/formats.c

    r1719 r2551  
    55#include <string.h>
    66
    7 #include "P_.h"
    87#include "astro.h"
    98#include "preferences.h"
    … …  
    1716 *      600:    <w>:mm.m
    1817 *      60:     <w>:mm
    19  */
    20 void
    21 fs_sexa (out, a, w, fracbase)
    22 char *out;
    23 double a;
    24 int w;
    25 int fracbase;
    26 {
     18 * return number of characters written to out, not counting final '\0'.
     19 */
     20int
     21fs_sexa (char *out, double a, int w, int fracbase)
     22{
     23        char *out0 = out;
    2724        unsigned long n;
    2825        int d;
    … …  
    4441        /* form the whole part; "negative 0" is a special case */
    4542        if (isneg && d == 0)
    46             (void) sprintf (out, "%*s-0", w-2, "");
     43            out += sprintf (out, "%*s-0", w-2, "");
    4744        else
    48             (void) sprintf (out, "%*d", w, isneg ? -d : d);
    49         out += w;
     45            out += sprintf (out, "%*d", w, isneg ? -d : d);
    5046
    5147        /* do the rest */
    … …  
    5349        case 60:        /* dd:mm */
    5450            m = f/(fracbase/60);
    55             (void) sprintf (out, ":%02d", m);
     51            out += sprintf (out, ":%02d", m);
    5652            break;
    5753        case 600:       /* dd:mm.m */
    58             (void) sprintf (out, ":%02d.%1d", f/10, f%10);
     54            out += sprintf (out, ":%02d.%1d", f/10, f%10);
    5955            break;
    6056        case 3600:      /* dd:mm:ss */
    6157            m = f/(fracbase/60);
    6258            s = f%(fracbase/60);
    63             (void) sprintf (out, ":%02d:%02d", m, s);
     59            out += sprintf (out, ":%02d:%02d", m, s);
    6460            break;
    6561        case 36000:     /* dd:mm:ss.s*/
    6662            m = f/(fracbase/60);
    6763            s = f%(fracbase/60);
    68             (void) sprintf (out, ":%02d:%02d.%1d", m, s/10, s%10);
     64            out += sprintf (out, ":%02d:%02d.%1d", m, s/10, s%10);
    6965            break;
    7066        case 360000:    /* dd:mm:ss.ss */
    7167            m = f/(fracbase/60);
    7268            s = f%(fracbase/60);
    73             (void) sprintf (out, ":%02d:%02d.%02d", m, s/100, s%100);
     69            out += sprintf (out, ":%02d:%02d.%02d", m, s/100, s%100);
    7470            break;
    7571        default:
    7672            printf ("fs_sexa: unknown fracbase: %d\n", fracbase);
    77             exit(1);
     73            abort();
    7874        }
     75
     76        return (out - out0);
    7977}
    8078
    8179/* put the given modified Julian date, jd, in out[] according to preference
    8280 * format.
    83  */
    84 void
    85 fs_date (out, jd)
    86 char out[];
    87 double jd;
     81 * return number of characters written to out, not counting final '\0'.
     82 */
     83int
     84fs_date (char out[], double jd)
    8885{
    8986        int p = pref_get (PREF_DATE_FORMAT);
     87        char *out0 = out;
    9088        int m, y;
    9189        double d;
    … …  
    10098        switch (p) {
    10199        case PREF_YMD:
    102             (void) sprintf (out, "%4d/%02d/%02.6g", y, m, d);
     100            out += sprintf (out, "%4d/%02d/%02.6g", y, m, d);
    103101            break;
    104102        case PREF_DMY:
    105             (void) sprintf (out, "%2.6g/%02d/%-4d", d, m, y);
     103            out += sprintf (out, "%2.6g/%02d/%-4d", d, m, y);
    106104            break;
    107105        case PREF_MDY:
    108             (void) sprintf (out, "%2d/%02.6g/%-4d", m, d, y);
     106            out += sprintf (out, "%2d/%02.6g/%-4d", m, d, y);
    109107            break;
    110108        default:
    111109            printf ("fs_date: bad date pref: %d\n", p);
    112             exit (1);
     110            abort();
    113111        }
    114 }
    115 
    116 /* scan a sexagesimal string and update a double. the string is of the form
    117  *   H:M:S. a negative value may be indicated by a '-' char before any
    118  *   component. All components may be integral or real. In addition to ':' the
    119  *    separator may also be '/' or ';' or ',' or '-'.
    120  * any components not specified in bp[] are copied from old's.
    121  *   eg:  ::10  only changes S
    122  *        10    only changes H
    123  *        10:0  changes H and M
    124  */
    125 void
    126 f_scansex (o, bp, np)
    127 double o;       /* old value */
    128 char bp[];      /* input string */
    129 double *np;     /* new value */
    130 {
    131         double oh, om, os;
    132         double nh, nm, ns;
    133         int nneg;
    134         double tmp;
    135         char c;
    136 
    137         /* fast macro to scan each segment of bp */
    138 #define SCANSEX(new,old)                                                \
    139         if (*bp == '-') {nneg = 1; bp++;}                               \
    140         if (sscanf (bp, "%lf", &new) != 1) new = old;                   \
    141         while ((c=(*bp)) != '\0') {                                     \
    142             bp++;                                                       \
    143             if (c==':' || c=='/' || c==';' || c==',' || c=='-')         \
    144                 break;                                                  \
    145         }
    146 
    147         /* get h:m:s components of o in case bp[] defers.
    148          * and constrain to positive for now.
    149          */
    150         if (o < 0.0)
    151             o = -o;
    152         oh = floor(o);
    153         tmp = (o - oh)*60.0;
    154         om = floor(tmp);
    155         os = (tmp - om)*60.0;
    156        
    157         /* round to small portion of a second to reduce differencing errors. */
    158         if (os > 59.99999) {
    159             os = 0.0;
    160             if ((om += 1.0) >= 60.0) {
    161                 om = 0.0;
    162                 oh += 1.0;
    163             }
    164         }
    165 
    166         /* scan each component of the buffer */
    167         while (*bp == ' ') bp++;
    168         nneg = 0;
    169         SCANSEX (nh, oh)
    170         SCANSEX (nm, om)
    171         SCANSEX (ns, os)
    172 
    173         /* back to hours */
    174         tmp = ns/3600.0 + nm/60.0 + nh;
    175         if (nneg)
    176             tmp = -tmp;
    177 
    178         *np = tmp;
     112
     113        return (out - out0);
     114}
     115
     116
     117/* convert sexagesimal string str AxBxC to double.
     118 *   x can be anything non-numeric. Any missing A, B or C will be assumed 0.
     119 *   optional - and + can be anywhere.
     120 * return 0 if ok, -1 if can't find a thing.
     121 */
     122int
     123f_scansexa (
     124const char *str0,       /* input string */
     125double *dp)             /* cracked value, if return 0 */
     126{
     127        double a = 0, b = 0, c = 0;
     128        char str[128];
     129        char *neg;
     130        int r;
     131
     132        /* copy str0 so we can play with it */
     133        strncpy (str, str0, sizeof(str)-1);
     134        str[sizeof(str)-1] = '\0';
     135
     136        neg = strchr(str, '-');
     137        if (neg)
     138            *neg = ' ';
     139        r = sscanf (str, "%lf%*[^0-9]%lf%*[^0-9]%lf", &a, &b, &c);
     140        if (r < 1)
     141            return (-1);
     142        *dp = a + b/60 + c/3600;
     143        if (neg)
     144            *dp *= -1;
     145        return (0);
    179146}
    180147
    … …  
    182149 *   PREF_DATE_FORMAT preference into its components. allow the day to be a
    183150 *   floating point number,
    184  * a lone component is always a year if it contains a decimal point or pref
    185  *   is MDY or DMY and it is not a reasonable month or day value, respectively.
    186  * leave any unspecified component unchanged.
    187  * actually, the slashes may be anything but digits or a decimal point.
     151 * actually the slashes may be anything but digits.
     152 * a lone component with a decimal point is considered a year.
    188153 */
    189154void
    190 f_sscandate (bp, pref, m, d, y)
    191 char *bp;
    192 int pref;       /* one of PREF_X for PREF_DATE_FORMAT */
    193 int *m, *y;
    194 double *d;
    195 {
    196         double comp[3]; /* the three components */
    197         int set[3];     /* set[i] is 1 if component i is present */
    198         int in;         /* scan state: 1 while we are in a component */
    199         int ncomp;      /* incremented as each component is discovered */
    200         int ndp;        /* number of decimal points in last component */
    201         char *bp0 = NULL, c;
    202 
    203         set[0] = set[1] = set[2] = 0;
    204         ncomp = 0;
    205         in = 0;
    206         ndp = 0;
    207 
    208         /* scan the input string.
    209          * '\0' counts as a component terminator too.
    210          */
    211         do {
    212             /* here, *bp is the next character to be investigated */
    213             c = *bp;
    214             if (c == '.' || isdigit(c) || (c == '-' && !in)) {
    215                 /* found what passes for a floating point number */
    216                 if (in == 0) {
    217                     /* save the beginning of it in bp0 and init ndp */
    218                     bp0 = bp;
    219                     in = 1;
    220                     ndp = 0;
    221                 }
    222                 if (c == '.')
    223                     ndp++;
    224             } else if (c != ' ') {
    225                 /* not in a number now ... */
    226                 if (in) {
    227                     /* ... but we *were* in a number, so it counts */
    228                     if (ncomp >= 3)
    229                         return; /* too many components.. just bail */
    230                     comp[ncomp] = atod (bp0);
    231                     set[ncomp] = 1;
    232                     in = 0;
    233                 }
    234 
    235                 /* regardless, a non-blank non-float means another component */
    236                 ncomp++;
     155f_sscandate (
     156char *bp,
     157int pref,       /* one of PREF_X for PREF_DATE_FORMAT */
     158int *m,
     159double *d,
     160int *y)
     161{
     162        double c[3]; /* the three components */
     163        int n;
     164
     165        memset (c, 0, sizeof(c));
     166        n = sscanf (bp, "%lf%*[^0-9]%lf%*[^0-9]%lf", &c[0], &c[1], &c[2]);
     167        if (n == 1 && (strchr(bp, '.')
     168                        || (pref == PREF_MDY && (c[0] < 1 || c[0] > 12))
     169                        || (pref == PREF_DMY && (c[0] < 1 || c[0] > 31)))) {
     170            double Mjd;
     171            year_mjd (c[0], &Mjd);
     172            mjd_cal (Mjd, m, d, y);
     173        } else {
     174            switch (pref) {
     175            case PREF_MDY:
     176                if (n > 0 && c[0] != 0)
     177                    *m = (int)c[0];
     178                if (n > 1 && c[1] != 0)
     179                    *d = c[1];
     180                if (n > 2 && c[2] != 0)
     181                    *y = (int)c[2];
     182                break;
     183            case PREF_YMD:
     184                if (n > 0 && c[0] != 0)
     185                    *y = (int)c[0];
     186                if (n > 1 && c[1] != 0)
     187                    *m = (int)c[1];
     188                if (n > 2 && c[2] != 0)
     189                    *d = c[2];
     190                break;
     191            case PREF_DMY:
     192                if (n > 1 && c[0] != 0)
     193                    *d = c[0];
     194                if (n > 2 && c[1] != 0)
     195                    *m = (int)c[1];
     196                if (n > 3 && c[2] != 0)
     197                    *y = (int)c[2];
     198                break;
    237199            }
    238             bp++;
    239         } while (c);
    240 
    241         /* it's a decimal year if there is exactly one component and
    242          *   it contains a decimal point
    243          *   or we are in MDY format and it's not in the range 1..12
    244          *   or we are in DMY format and it's not int the rangle 1..31
    245          */
    246         if (ncomp == 1 &&
    247                 (ndp > 0
    248                  || (pref == PREF_MDY && !(comp[0] >= 1 && comp[0] <= 12))
    249                  || (pref == PREF_DMY && !(comp[0] >= 1 && comp[0] <= 31)))) {
    250             double Mjd;
    251             year_mjd (comp[0], &Mjd);
    252             mjd_cal (Mjd, m, d, y);
    253             return;
    254200        }
    255 
    256         switch (pref) {
    257         case PREF_MDY:
    258             if (set[0]) *m = (int)comp[0];
    259             if (set[1]) *d = comp[1];
    260             if (set[2]) *y = (int)comp[2];
    261             break;
    262         case PREF_YMD:
    263             if (set[0]) *y = (int)comp[0];
    264             if (set[1]) *m = (int)comp[1];
    265             if (set[2]) *d = comp[2];
    266             break;
    267         case PREF_DMY:
    268             if (set[0]) *d = comp[0];
    269             if (set[1]) *m = (int)comp[1];
    270             if (set[2]) *y = (int)comp[2];
    271             break;
    272         }
    273 }
    274 
    275 /* given a string of the form xx:mm[:ss] or xx:mm.dd, convert it to a double at
    276  * *dp. actually, ':' may also be ';', '/' or ',' too. all components may be
    277  * floats.
    278  * return -1 if trouble, else 0
    279  */
    280 int
    281 scansex (str, dp)
    282 char *str;
    283 double *dp;
    284 {
    285         double x, m = 0.0, s = 0.0;
    286         int isneg;
    287         int nf;
    288 
    289         while (isspace(*str))
    290             str++;
    291         if (*str == '-') {
    292             isneg = 1;
    293             str++;
    294         } else
    295             isneg = 0;
    296 
    297         nf = sscanf (str, "%lf%*[:;/,]%lf%*[:;/,]%lf", &x, &m, &s);
    298         if (nf < 1)
    299             return (-1);
    300         *dp = x + m/60.0 + s/3600.0;
    301         if (isneg)
    302             *dp = - *dp;
    303         return (0);
    304201}
    305202
    306203/* For RCS Only -- Do Not Edit */
    307 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: formats.c,v $ $Date: 2001-10-22 12:08:27 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     204static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: formats.c,v $ $Date: 2004-06-15 16:52:38 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/helio.c

    r1719 r2551  
    33#include <math.h>
    44
    5 #include "P_.h"
    65#include "astro.h"
    76
    … …  
    1211 */
    1312void
    14 heliocorr (jd, ra, dec, hcp)
    15 double jd, ra, dec;
    16 double *hcp;
     13heliocorr (double jd, double ra, double dec, double *hcp)
    1714{
    1815        double e;       /* obliquity of ecliptic */
    … …  
    5249
    5350/* For RCS Only -- Do Not Edit */
    54 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: helio.c,v $ $Date: 2001-10-22 12:08:27 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     51static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: helio.c,v $ $Date: 2004-06-15 16:52:38 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/libration.c

    r1719 r2551  
    3535#include <math.h>
    3636
    37 #include "P_.h"
    3837#include "astro.h"
    3938
    … …  
    7170   answer in arc seconds.  */
    7271static double
    73 mods3600 (x)
    74      double x;
     72mods3600 (double x)
    7573{
    7674  double y;
    … …  
    8583   n is the highest harmonic to compute.  */
    8684static int
    87 sscc (k, arg, n)
    88      int k;
    89      double arg;
    90      int n;
     85sscc (int k, double arg, int n)
    9186{
    9287  double cu, su, cv, sv, s;
    … …  
    123118
    124119static int
    125 dargs (J, plan)
    126      double J;
    127      struct plantbl *plan;
     120dargs (double J, struct plantbl *plan)
    128121{
    129122  double T2, w;
    … …  
    229222
    230223static double
    231 gplan (J, plan)
    232      double J;
    233      struct plantbl *plan;
     224gplan (double J, struct plantbl *plan)
    234225{
    235226  double su, cu, sv, cv;
    … …  
    22132204 */
    22142205void
    2215 llibration (JD, llatp, llonp)
    2216 double JD;
    2217 double *llatp;
    2218 double *llonp;
     2206llibration (double JD, double *llatp, double *llonp)
    22192207{
    2220         double lon, lat;        /* arc seconds */
     2208        double lg, lt;  /* arc seconds */
    22212209       
    2222         lon = gplan (JD, &liblon);
    2223         lat = gplan (JD, &liblat);
    2224 
    2225         *llonp = degrad (lon/3600.0);
    2226         *llatp = degrad (lat/3600.0);
     2210        lg = gplan (JD, &liblon);
     2211        lt = gplan (JD, &liblat);
     2212
     2213        *llonp = degrad (lg/3600.0);
     2214        *llatp = degrad (lt/3600.0);
    22272215}
    22282216
    22292217/* For RCS Only -- Do Not Edit */
    2230 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: libration.c,v $ $Date: 2001-10-22 12:08:27 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     2218static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: libration.c,v $ $Date: 2004-06-15 16:52:39 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/misc.c

    r1719 r2551  
    66#include <stdio.h>
    77#include <math.h>
    8 
    9 #if defined(__STDC__)
    108#include <stdlib.h>
    119#include <string.h>
    12 #else
    13 extern double atof();
    14 #endif
    15 
    16 #include "P_.h"
     10
    1711#include "astro.h"
    18 #include "circum.h"
    1912
    2013/* zero from loc for len bytes */
    2114void
    22 zero_mem (loc, len)
    23 void *loc;
    24 unsigned len;
     15zero_mem (void *loc, unsigned len)
    2516{
    2617        (void) memset (loc, 0, len);
    … …  
    3324 */
    3425int
    35 tickmarks (min, max, numdiv, ticks)
    36 double min, max;
    37 int numdiv;
    38 double ticks[];
     26tickmarks (double min, double max, int numdiv, double ticks[])
    3927{
    4028        static int factor[] = { 1, 2, 5 };
    … …  
    4735        minscale = fabs(max - min);
    4836        delta = minscale/numdiv;
    49         for (n=0; n < sizeof(factor)/sizeof(factor[0]); n++) {
     37        for (n=0; n < (int)(sizeof(factor)/sizeof(factor[0])); n++) {
    5038            double scale;
    5139            double x = delta/factor[n];
    … …  
    6755 */
    6856char *
    69 obj_description (op)
    70 Obj *op;
    71 {
    72         static struct {
    73             char class;
     57obj_description (Obj *op)
     58{
     59        typedef struct {
     60            char classcode;
    7461            char *desc;
    75         } fixed_class_map[] = {
     62        } CC;
     63
     64#define NFCM    ((int)(sizeof(fixed_class_map)/sizeof(fixed_class_map[0])))
     65        static CC fixed_class_map[] = {
    7666            {'A', "Cluster of Galaxies"},
    77             {'B', "Binary Star"},
     67            {'B', "Binary System"},
    7868            {'C', "Globular Cluster"},
    7969            {'D', "Double Star"},
    … …  
    9484            {'U', "Cluster, with nebulosity"},
    9585            {'V', "Variable Star"},
     86            {'Y', "Supernova"},
    9687        };
    97 #define NFCM    (sizeof(fixed_class_map)/sizeof(fixed_class_map[0]))
     88
     89#define NBCM    ((int)(sizeof(binary_class_map)/sizeof(binary_class_map[0])))
     90        static CC binary_class_map[] = {
     91            {'a', "Astrometric binary"},
     92            {'c', "Cataclysmic variable"},
     93            {'e', "Eclipsing binary"},
     94            {'x', "High-mass X-ray binary"},
     95            {'y', "Low-mass X-ray binary"},
     96            {'o', "Occultation binary"},
     97            {'s', "Spectroscopic binary"},
     98            {'t', "1-line spectral binary"},
     99            {'u', "2-line spectral binary"},
     100            {'v', "Spectrum binary"},
     101            {'b', "Visual binary"},
     102            {'d', "Visual binary, apparent"},
     103            {'q', "Visual binary, optical"},
     104            {'r', "Visual binary, physical"},
     105            {'p', "Exoplanet"},
     106        };
    98107
    99108        switch (op->o_type) {
    … …  
    102111                int i;
    103112                for (i = 0; i < NFCM; i++)
    104                     if (fixed_class_map[i].class == op->f_class)
     113                    if (fixed_class_map[i].classcode == op->f_class)
    105114                        return (fixed_class_map[i].desc);
    106115            }
    107             return (op->o_type == FIXED ? "Fixed" : "Field Star");
     116            return ("Fixed");
    108117        case PARABOLIC:
    109118            return ("Solar - Parabolic");
    … …  
    112121        case ELLIPTICAL:
    113122            return ("Solar - Elliptical");
    114         case PLANET:
    115             return ("Planet");
     123        case BINARYSTAR:
     124            if (op->f_class) {
     125                int i;
     126                for (i = 0; i < NFCM; i++)
     127                    if (binary_class_map[i].classcode == op->f_class)
     128                        return (binary_class_map[i].desc);
     129            }
     130            return ("Binary system");
     131        case PLANET: {
     132            static char nsstr[16];
     133            static Obj *biop;
     134
     135            if (op->pl_code == SUN)
     136                return ("Star");
     137            if (op->pl_code == MOON)
     138                return ("Moon of Earth");
     139            if (op->pl_moon == X_PLANET)
     140                return ("Planet");
     141            if (!biop)
     142                getBuiltInObjs (&biop);
     143            sprintf (nsstr, "Moon of %s", biop[op->pl_code].o_name);
     144            return (nsstr);
     145            }
    116146        case EARTHSAT:
    117147            return ("Earth Sat");
    118148        default:
    119149            printf ("obj_description: unknown type: 0x%x\n", op->o_type);
    120             exit (1);
     150            abort();
    121151            return (NULL);      /* for lint */
    122152        }
    … …  
    126156 */
    127157void
    128 now_lst (np, lstp)
    129 Now *np;
    130 double *lstp;
     158now_lst (Now *np, double *lstp)
    131159{
    132160        static double last_mjd = -23243, last_lng = 121212, last_lst;
    … …  
    156184 */
    157185void
    158 radec2ha (np, ra, dec, hap)
    159 Now *np;
    160 double ra, dec;
    161 double *hap;
     186radec2ha (Now *np, double ra, double dec, double *hap)
    162187{
    163188        double ha, lst;
    … …  
    182207 */
    183208int
    184 lc (cx, cy, cw, x1, y1, x2, y2, sx1, sy1, sx2, sy2)
    185 int cx, cy, cw;                 /* circle bounding box corner and width */
    186 int x1, y1, x2, y2;             /* line segment endpoints */
    187 int *sx1, *sy1, *sx2, *sy2;     /* segment inside the circle */
     209lc (
     210int cx, int cy, int cw,                 /* circle bbox corner and width */
     211int x1, int y1, int x2, int y2,         /* line segment endpoints */
     212int *sx1, int *sy1, int *sx2, int *sy2) /* segment inside the circle */
    188213{
    189214        int dx = x2 - x1;
    … …  
    240265 */
    241266void
    242 hg_mag (h, g, rp, rho, rsn, mp)
    243 double h, g;
    244 double rp;      /* sun-obj dist, AU */
    245 double rho;     /* earth-obj dist, AU */
    246 double rsn;     /* sun-earth dist, AU */
    247 double *mp;
     267hg_mag (
     268double h, double g,
     269double rp,      /* sun-obj dist, AU */
     270double rho,     /* earth-obj dist, AU */
     271double rsn,     /* sun-earth dist, AU */
     272double *mp)
    248273{
    249274        double psi_t, Psi_1, Psi_2, beta;
    … …  
    265290        psi_t = pow (tb2, 1.22);
    266291        Psi_2 = exp(-1.87*psi_t);
    267         *mp = h + 5.0*log10(rp*rho) - 2.5*log10((1-g)*Psi_1 + g*Psi_2);
     292        *mp = h + 5.0*log10(rp*rho);
     293        if (Psi_1 || Psi_2) *mp -= 2.5*log10((1-g)*Psi_1 + g*Psi_2);
    268294}
    269295
    … …  
    273299 */
    274300int
    275 magdiam (fmag, magstp, scale, mag, size)
    276 int fmag;       /* faintest mag */
    277 int magstp;     /* mag range per dot size */
    278 double scale;   /* rads per pixel */
    279 double mag;     /* magnitude */
    280 double size;    /* rads, or 0 */
     301magdiam (
     302int fmag,       /* faintest mag */
     303int magstp,     /* mag range per dot size */
     304double scale,   /* rads per pixel */
     305double mag,     /* magnitude */
     306double size)    /* rads, or 0 */
    281307{
    282308        int diam, sized;
    … …  
    295321 */
    296322void
    297 gk_mag (g, k, rp, rho, mp)
    298 double g, k;
    299 double rp;      /* sun-obj dist, AU */
    300 double rho;     /* earth-obj dist, AU */
    301 double *mp;
     323gk_mag (
     324double g, double k,
     325double rp,      /* sun-obj dist, AU */
     326double rho,     /* earth-obj dist, AU */
     327double *mp)
    302328{
    303329        *mp = g + 5.0*log10(rho) + 2.5*k*log10(rp);
    … …  
    309335 */
    310336double
    311 atod (buf)
    312 char *buf;
    313 {
    314         return (atof (buf));
     337atod (char *buf)
     338{
     339        return (strtod (buf, NULL));
    315340}
    316341
    … …  
    331356 */
    332357void
    333 solve_sphere (A, b, cc, sc, cap, Bp)
    334 double A, b;
    335 double cc, sc;
    336 double *cap, *Bp;
     358solve_sphere (double A, double b, double cc, double sc, double *cap, double *Bp)
    337359{
    338360        double cb = cos(b), sb = sin(b);
    339         double cA = cos(A);
     361        double sA, cA = cos(A);
     362        double x, y;
    340363        double ca;
    341364        double B;
    … …  
    350373            return;
    351374
    352         if (cc > .99999) {
    353             /* as c approaches 0, B approaches pi - A */
    354             B = PI - A;
    355         } else if (cc < -.99999) {
    356             /* as c approaches PI, B approaches A */
    357             B = A;
    358         } else {
    359             /* compute cB and sB and remove common factor of sa from quotient.
    360              * be careful where B causes atan to blow.
    361              */
    362             double sA = sin(A);
    363             double x, y;
    364 
    365             y = sA*sb*sc;
    366             x = cb - ca*cc;
    367        
    368             if (fabs(x) < 1e-5)
    369                 B = y < 0 ? 3*PI/2 : PI/2;
    370             else
    371                 B = atan2 (y, x);
    372         }
     375        sA = sin(A);
     376        y = sA*sb*sc;
     377        x = cb - ca*cc;
     378        B = y ? (x ? atan2(y,x) : (y>0 ? PI/2 : -PI/2)) : (x>=0 ? 0 : PI);
    373379
    374380        *Bp = B;
    … …  
    435441 */
    436442double
    437 delra (dra)
    438 double dra;
     443delra (double dra)
    439444{
    440445        double fdra = fmod(fabs(dra), 2*PI);
    … …  
    449454 */
    450455int
    451 is_deepsky (op)
    452 Obj *op;
     456is_deepsky (Obj *op)
    453457{
    454458        int deepsky = 0;
    … …  
    473477
    474478/* For RCS Only -- Do Not Edit */
    475 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: misc.c,v $ $Date: 2001-10-22 12:08:27 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     479static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: misc.c,v $ $Date: 2004-06-15 16:52:39 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/mjd.c

    r1719 r2551  
    44#include <math.h>
    55
    6 #include "P_.h"
    76#include "astro.h"
    87
    … …  
    1211 */
    1312void
    14 cal_mjd (mn, dy, yr, mjd)
    15 int mn, yr;
    16 double dy;
    17 double *mjd;
     13cal_mjd (int mn, double dy, int yr, double *mjp)
    1814{
    1915        static double last_mjd, last_dy;
    … …  
    2319
    2420        if (mn == last_mn && yr == last_yr && dy == last_dy) {
    25             *mjd = last_mjd;
     21            *mjp = last_mjd;
    2622            return;
    2723        }
    … …  
    4945        d = (int)(30.6001*(m+1));
    5046
    51         *mjd = b + c + d + dy - 0.5;
     47        *mjp = b + c + d + dy - 0.5;
    5248
    5349        last_mn = mn;
    5450        last_dy = dy;
    5551        last_yr = yr;
    56         last_mjd = *mjd;
     52        last_mjd = *mjp;
    5753}
    5854
    5955/* given the modified Julian date (number of days elapsed since 1900 jan 0.5,),
    60  * mjd, return the calendar date in months, *mn, days, *dy, and years, *yr.
    61  */
    62 void
    63 mjd_cal (mjd, mn, dy, yr)
    64 double mjd;
    65 int *mn, *yr;
    66 double *dy;
    67 {
    68         static double last_mjd, last_dy;
     56 * mj, return the calendar date in months, *mn, days, *dy, and years, *yr.
     57 */
     58void
     59mjd_cal (double mj, int *mn, double *dy, int *yr)
     60{
     61        static double last_mj, last_dy;
    6962        static int last_mn, last_yr;
    7063        double d, f;
    … …  
    7467         * 0 is noon the last day of 1899.
    7568         */
    76         if (mjd == 0.0) {
     69        if (mj == 0.0) {
    7770            *mn = 12;
    7871            *dy = 31.5;
    … …  
    8174        }
    8275
    83         if (mjd == last_mjd) {
     76        if (mj == last_mj) {
    8477            *mn = last_mn;
    8578            *yr = last_yr;
    … …  
    8881        }
    8982
    90         d = mjd + 0.5;
     83        d = mj + 0.5;
    9184        i = floor(d);
    9285        f = d-i;
    … …  
    118111        last_dy = *dy;
    119112        last_yr = *yr;
    120         last_mjd = mjd;
     113        last_mj = mj;
    121114}
    122115
    … …  
    126119 */
    127120int
    128 mjd_dow (mjd, dow)
    129 double mjd;
    130 int *dow;
     121mjd_dow (double mj, int *dow)
    131122{
    132123        /* cal_mjd() uses Gregorian dates on or after Oct 15, 1582.
    … …  
    137128         * can not easily be accounted for from the cal_mjd() number...
    138129         */
    139         if (mjd < -53798.5) {
     130        if (mj < -53798.5) {
    140131            /* pre sept 14, 1752 too hard to correct |:-S */
    141132            return (-1);
    142133        }
    143         *dow = ((long)floor(mjd-.5) + 1) % 7;/* 1/1/1900 (mjd 0.5) is a Monday*/
     134        *dow = ((long)floor(mj-.5) + 1) % 7;/* 1/1/1900 (mj 0.5) is a Monday*/
    144135        if (*dow < 0)
    145136            *dow += 7;
    … …  
    149140/* given a year, return whether it is a leap year */
    150141int
    151 isleapyear (y)
    152 int y;
     142isleapyear (int y)
    153143{
    154144        return ((y%4==0 && y%100!=0) || y%400==0);
    … …  
    157147/* given a mjd, return the the number of days in the month.  */
    158148void
    159 mjd_dpm (mjd, ndays)
    160 double mjd;
    161 int *ndays;
     149mjd_dpm (double mj, int *ndays)
    162150{
    163151        static short dpm[] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
    … …  
    165153        double d;
    166154
    167         mjd_cal (mjd, &m, &d, &y);
     155        mjd_cal (mj, &m, &d, &y);
    168156        *ndays = (m==2 && isleapyear(y)) ? 29 : dpm[m-1];
    169157}
    … …  
    171159/* given a mjd, return the year and number of days since 00:00 Jan 1 */
    172160void
    173 mjd_dayno (mjd, yr, dy)
    174 double mjd;
    175 int *yr;
    176 double *dy;
     161mjd_dayno (double mj, int *yr, double *dy)
    177162{
    178163        double yrd;
    … …  
    180165        int dpy;
    181166
    182         mjd_year (mjd, &yrd);
     167        mjd_year (mj, &yrd);
    183168        *yr = yri = (int)yrd;
    184169        dpy = isleapyear(yri) ? 366 : 365;
    … …  
    188173/* given a mjd, return the year as a double. */
    189174void
    190 mjd_year (mjd, yr)
    191 double mjd;
    192 double *yr;
    193 {
    194         static double last_mjd, last_yr;
     175mjd_year (double mj, double *yr)
     176{
     177        static double last_mj, last_yr;
    195178        int m, y;
    196179        double d;
    197180        double e0, e1;  /* mjd of start of this year, start of next year */
    198181
    199         if (mjd == last_mjd) {
     182        if (mj == last_mj) {
    200183            *yr = last_yr;
    201184            return;
    202185        }
    203186
    204         mjd_cal (mjd, &m, &d, &y);
     187        mjd_cal (mj, &m, &d, &y);
    205188        if (y == -1) y = -2;
    206189        cal_mjd (1, 1.0, y, &e0);
    207190        cal_mjd (1, 1.0, y+1, &e1);
    208         *yr = y + (mjd - e0)/(e1 - e0);
    209 
    210         last_mjd = mjd;
     191        *yr = y + (mj - e0)/(e1 - e0);
     192
     193        last_mj = mj;
    211194        last_yr = *yr;
    212195}
    … …  
    214197/* given a decimal year, return mjd */
    215198void
    216 year_mjd (y, mjd)
    217 double y;
    218 double *mjd;
     199year_mjd (double y, double *mjp)
    219200{
    220201        double e0, e1;  /* mjd of start of this year, start of next year */
    … …  
    224205        cal_mjd (1, 1.0, yf, &e0);
    225206        cal_mjd (1, 1.0, yf+1, &e1);
    226         *mjd = e0 + (y - yf)*(e1-e0);
     207        *mjp = e0 + (y - yf)*(e1-e0);
    227208}
    228209
    229210/* round a time in days, *t, to the nearest second, IN PLACE. */
    230211void
    231 rnd_second (t)
    232 double *t;
     212rnd_second (double *t)
    233213{
    234214        *t = floor(*t*SPD+0.5)/SPD;
    … …  
    237217/* given an mjd, truncate it to the beginning of the whole day */
    238218double
    239 mjd_day(jd)
    240 double jd;
    241 {
    242         return (floor(jd-0.5)+0.5);
     219mjd_day(double mj)
     220{
     221        return (floor(mj-0.5)+0.5);
    243222}
    244223
    245224/* given an mjd, return the number of hours past midnight of the whole day */
    246225double
    247 mjd_hr(jd)
    248 double jd;
    249 {
    250         return ((jd-mjd_day(jd))*24.0);
     226mjd_hr(double mj)
     227{
     228        return ((mj-mjd_day(mj))*24.0);
    251229}
    252230
    … …  
    254232 */
    255233void
    256 range (v, r)
    257 double *v, r;
     234range (double *v, double r)
    258235{
    259236        *v -= r*floor(*v/r);
    260237}
    261238
     239/* insure 0 <= ra < 2PI and -PI/2 <= dec <= PI/2. if dec needs
     240 * complimenting, reflect ra too
     241 */
     242void
     243radecrange (double *ra, double *dec)
     244{
     245        if (*dec < -PI/2) {
     246            *dec = -PI - *dec;
     247            *ra += PI;
     248        } else if (*dec > PI/2) {
     249            *dec = PI - *dec;
     250            *ra += PI;
     251        }
     252        range (ra, 2*PI);
     253}
     254
    262255/* For RCS Only -- Do Not Edit */
    263 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: mjd.c,v $ $Date: 2001-10-22 12:08:27 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     256static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: mjd.c,v $ $Date: 2004-06-15 16:52:39 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/moon.c

    r1719 r2551  
    6262 */
    6363
     64#include <stdlib.h>
     65#include <stdio.h>
    6466#include <math.h>
    65 #if defined(__STDC__)
    66 #include <stdlib.h>
    67 #endif
    68 
    69 #include "P_.h"
     67
    7068#include "astro.h"
    7169
    … …  
    8684};
    8785
    88 static double mods3600 P_((double x));
    89 static void mean_elements P_((double JED));
    90 static int sscc P_((int k, double arg, int n));
    91 static int g2plan P_((double J, struct plantbl *plan, double *pobj, int flag));
    92 static double g1plan P_((double J, struct plantbl *plan));
    93 static int gecmoon P_((double J, struct plantbl *lrtab,
    94         struct plantbl *lattab, double *pobj));
     86static double mods3600 (double x);
     87static void mean_elements (double JED);
     88static int sscc (int k, double arg, int n);
     89static int g2plan (double J, struct plantbl *plan, double *pobj, int flag);
     90static double g1plan (double J, struct plantbl *plan);
     91static int gecmoon (double J, struct plantbl *lrtab,
     92        struct plantbl *lattab, double *pobj);
    9593
    9694/* time points */
    … …  
    103101static double Args[NARGS];
    104102static double LP_equinox;
    105 static double NF_arcsec;
    106103static double Ea_arcsec;
    107104static double pA_precession;
    … …  
    28322829   answer in arc seconds  */
    28332830static double
    2834 mods3600(x)
    2835      double x;
     2831mods3600(double x)
    28362832{
    28372833  double y;
    … …  
    28482844
    28492845static void
    2850 mean_elements (JED)
    2851      double JED;
     2846mean_elements (double JED)
    28522847{
    28532848  double x, T, T2;
    … …  
    29322927              - 7.5311878482337989e-04) * T /* F, t^3 */
    29332928             - 1.3117809789650071e+01) * T2; /* F, t^2 */
    2934   NF_arcsec = x;
    29352929  Args[10] = x;
    29362930
    … …  
    30083002 */
    30093003static int
    3010 sscc (k, arg, n)
    3011      int k;
    3012      double arg;
    3013      int n;
     3004sscc (int k, double arg, int n)
    30143005{
    30153006  double cu, su, cv, sv, s;
    … …  
    30403031   of the same list of arguments.  */
    30413032static int
    3042 g2plan (J, plan, pobj, flag)
    3043      double J;
    3044      struct plantbl *plan;
    3045      double pobj[];
    3046      int flag;
     3033g2plan (double J, struct plantbl *plan, double pobj[], int flag)
    30473034{
    30483035  int i, j, k, m, k1, ip, np, nt;
    … …  
    31683155
    31693156static double
    3170 g1plan (J, plan)
    3171      double J;
    3172      struct plantbl *plan;
     3157g1plan (double J, struct plantbl *plan)
    31733158{
    31743159  int i, j, k, m, k1, ip, np, nt;
    … …  
    32713256 */
    32723257static int
    3273 gecmoon (J, lrtab, lattab, pobj)
    3274      double J;
    3275      struct plantbl *lrtab, *lattab;
    3276      double pobj[];
     3258gecmoon (double J, struct plantbl *lrtab, struct plantbl *lattab, double pobj[])
    32773259{
    32783260  double x;
    … …  
    32943276/*********** end stephen moshier's moon code ****************/
    32953277
    3296 static void moon_fast P_((double mjd, double *lam, double *bet,
    3297         double *hp, double *msp, double *mdp));
     3278static void moon_fast (double mj, double *lam, double *bet,
     3279        double *hp, double *msp, double *mdp);
    32983280
    32993281/* previous version (elwood):
    … …  
    33093291 */
    33103292static void
    3311 moon_fast (mjd, lam, bet, hp, msp, mdp)
    3312 double mjd;
    3313 double *lam, *bet, *hp;
    3314 double *msp, *mdp;
     3293moon_fast (double mj, double *lam, double *bet, double *hp, double *msp,
     3294double *mdp)
    33153295{
    33163296        double t, t2;
    … …  
    33243304        double m1, m2, m3, m4, m5, m6;
    33253305
    3326         t = mjd/36525.;
     3306        t = mj/36525.;
    33273307        t2 = t*t;
    33283308
    3329         m1 = mjd/27.32158213;
     3309        m1 = mj/27.32158213;
    33303310        m1 = 360.0*(m1-(long)m1);
    3331         m2 = mjd/365.2596407;
     3311        m2 = mj/365.2596407;
    33323312        m2 = 360.0*(m2-(long)m2);
    3333         m3 = mjd/27.55455094;
     3313        m3 = mj/27.55455094;
    33343314        m3 = 360.0*(m3-(long)m3);
    3335         m4 = mjd/29.53058868;
     3315        m4 = mj/29.53058868;
    33363316        m4 = 360.0*(m4-(long)m4);
    3337         m5 = mjd/27.21222039;
     3317        m5 = mj/27.21222039;
    33383318        m5 = 360.0*(m5-(long)m5);
    3339         m6 = mjd/6798.363307;
     3319        m6 = mj/6798.363307;
    33403320        m6 = 360.0*(m6-(long)m6);
    33413321
    … …  
    34563436 */
    34573437void
    3458 moon (mjd, lam, bet, rho, msp, mdp)
    3459 double mjd;
    3460 double *lam, *bet, *rho;
    3461 double *msp, *mdp;
     3438moon (double mj, double *lam, double *bet, double *rho, double *msp,
     3439double *mdp)
    34623440{
    34633441        double pobj[3], dt;
    34643442        double hp;
    34653443
    3466         if (mjd >= MOSHIER_BEGIN && mjd <= MOSHIER_END) {
     3444        if (mj >= MOSHIER_BEGIN && mj <= MOSHIER_END) {
    34673445                /* retard for light time */
    3468                 moon_fast (mjd, lam, bet, &hp, msp, mdp);
     3446                moon_fast (mj, lam, bet, &hp, msp, mdp);
    34693447                *rho = EarthRadius/AUKM/sin(hp);
    34703448                dt = *rho * 5.7755183e-3;       /* speed of light in a.u/day */
    3471                 gecmoon(mjd + MJD0 - dt, &moonlr, &moonlat, pobj);
     3449                gecmoon(mj + MJD0 - dt, &moonlr, &moonlat, pobj);
    34723450
    34733451                *lam = pobj[0];
    … …  
    34783456                *mdp = STR * Args[12];
    34793457        } else {
    3480                 moon_fast (mjd, lam, bet, &hp, msp, mdp);
     3458                moon_fast (mj, lam, bet, &hp, msp, mdp);
    34813459                *rho = EarthRadius/AUKM/sin(hp);
    34823460
    … …  
    34853463
    34863464/* For RCS Only -- Do Not Edit */
    3487 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: moon.c,v $ $Date: 2001-10-22 12:08:27 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     3465static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: moon.c,v $ $Date: 2004-06-15 16:52:39 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/mooncolong.c

    r1719 r2551  
    33
    44#include <stdio.h>
     5#include <stdlib.h>
    56#include <math.h>
    6 #if defined(__STDC__)
    7 #include <stdlib.h>
    8 #endif
    9 
    10 #include "P_.h"
     7
    118#include "astro.h"
    129
    13 static void Librations P_((double RAD, double LAMH, double BH, double OM,
    14     double F, double L, double L1, double *L0, double *B0));
    15 static void Moon P_((double RAD, double T, double T2, double LAM0, double R,
     10static void Librations (double RAD, double LAMH, double BH, double OM,
     11    double F, double L, double L1, double *L0, double *B0);
     12static void Moon (double RAD, double T, double T2, double LAM0, double R,
    1613    double M, double *F, double *L1, double *OM, double *LAM, double *B,
    17     double *DR, double *LAMH, double *BH));
    18 static void Sun P_((double RAD, double T, double T2, double *L, double *M,
    19     double *R, double *LAM0));
     14    double *DR, double *LAMH, double *BH);
     15static void Sun (double RAD, double T, double T2, double *L, double *M,
     16    double *R, double *LAM0);
    2017
    2118/* given a Julian date and a lunar location, find selenographic colongitude of
    … …  
    2623 */
    2724void
    28 moon_colong (jd, lt, lg, cp, kp, ap, sp)
    29 double jd;      /* jd */
    30 double lt, lg;  /* lat/long of location on moon, rads +N +E */
    31 double *cp;     /* selenographic colongitude (-lng of rising sun), rads */
    32 double *kp;     /* illuminated fraction of surface from Earth */
    33 double *ap;     /* sun altitude at location, rads */
    34 double *sp;     /* lunar latitude of subsolar point, rads */
     25moon_colong (
     26double jd,      /* jd */
     27double lt,      /* lat of location on moon, rads +N +E */
     28double lg,      /* long of location on moon, rads +N +E */
     29double *cp,     /* selenographic colongitude (-lng of rising sun), rads */
     30double *kp,     /* illuminated fraction of surface from Earth */
     31double *ap,     /* sun altitude at location, rads */
     32double *sp)     /* lunar latitude of subsolar point, rads */
    3533{
    3634        double RAD = .0174533;
    … …  
    8987
    9088static void
    91 Librations (RAD, LAMH, BH, OM, F, L, L1, L0, B0)
    92 double RAD;
    93 double LAMH;
    94 double BH;
    95 double OM;
    96 double F;
    97 double L;
    98 double L1;
    99 double *L0;
    100 double *B0;
     89Librations (double RAD, double LAMH, double BH, double OM, double F,
     90double L, double L1, double *L0, double *B0)
    10191{
    10292        double I, PSI, W, NUM, DEN, A, TEMP;
    … …  
    123113
    124114static void
    125 Moon (RAD, T, T2, LAM0, R, M, F, L1, OM, LAM, B, DR, LAMH, BH)
    126 double RAD;
    127 double T;
    128 double T2;
    129 double LAM0;
    130 double R;
    131 double M;
    132 double *F;
    133 double *L1;
    134 double *OM;
    135 double *LAM;
    136 double *B;
    137 double *DR;
    138 double *LAMH;
    139 double *BH;
     115Moon (double RAD, double T, double T2, double LAM0, double R, double M,
     116double *F, double *L1, double *OM, double *LAM, double *B, double *DR,
     117double *LAMH, double *BH)
    140118{
    141119        double T3, M1, D2, SUMR, SUML, DIST;
    … …  
    175153
    176154static void
    177 Sun (RAD, T, T2, L, M, R, LAM0)
    178 double RAD;
    179 double T;
    180 double T2;
    181 double *L;
    182 double *M;
    183 double *R;
    184 double *LAM0;
     155Sun (double RAD, double T, double T2, double *L, double *M, double *R,
     156double *LAM0)
    185157{
    186158        double T3, C, V, E, THETA, OM;
    … …  
    238210        if (ac != 2) {
    239211            fprintf (stderr, "%s: JD\n", av[0]);
    240             exit (1);
     212            abort();
    241213        }
    242214
    … …  
    262234
    263235/* For RCS Only -- Do Not Edit */
    264 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: mooncolong.c,v $ $Date: 2001-10-22 12:08:27 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     236static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: mooncolong.c,v $ $Date: 2004-06-15 16:52:39 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/nutation.c

    r1719 r2551  
    33 * (given to 0.001") EXACTLY over 750 days (1995 and 1996)
    44 */
     5#include <stdio.h>
    56#include <math.h>
    67
    7 #include "P_.h"
    88#include "astro.h"
    99
    … …  
    268268};
    269269
    270 /* given the modified JD, mjd, find the nutation in obliquity, *deps, and
     270/* given the modified JD, mj, find the nutation in obliquity, *deps, and
    271271 * the nutation in longitude, *dpsi, each in radians.
    272272 */
    273273void
    274 nutation (mjd, deps, dpsi)
    275 double mjd;
    276 double *deps;   /* on input:  precision parameter in arc seconds */
    277 double *dpsi;
     274nutation (
     275double mj,
     276double *deps,   /* on input:  precision parameter in arc seconds */
     277double *dpsi)
    278278{
    279         static double lastmjd = -10000, lastdeps, lastdpsi;
     279        static double lastmj = -10000, lastdeps, lastdpsi;
    280280        double T, T2, T3, T10;                  /* jul cent since J2000 */
    281281        double prec;                            /* series precis in arc sec */
    … …  
    287287                         */
    288288
    289         if (mjd == lastmjd) {
     289        if (mj == lastmj) {
    290290            *deps = lastdeps;
    291291            *dpsi = lastdpsi;
    … …  
    304304        prec *= NUT_SCALE/10;
    305305
    306         T = (mjd - J2000)/36525.;
     306        T = (mj - J2000)/36525.;
    307307        T2 = T * T;
    308308        T3 = T2 * T;
    … …  
    362362        lastdeps = degrad(lastdeps/3600./NUT_SCALE);
    363363
    364         lastmjd = mjd;
     364        lastmj = mj;
    365365        *deps = lastdeps;
    366366        *dpsi = lastdpsi;
    367367}
    368368
    369 /* given the modified JD, mjd, correct, IN PLACE, the right ascension *ra
     369/* given the modified JD, mj, correct, IN PLACE, the right ascension *ra
    370370 * and declination *dec (both in radians) for nutation.
    371371 */
    372372void
    373 nut_eq (mjd, ra, dec)
    374 double mjd, *ra, *dec;
     373nut_eq (double mj, double *ra, double *dec)
    375374{
    376         static double lastmjd = -10000;
     375        static double lastmj = -10000;
    377376        static double a[3][3];          /* rotation matrix */
    378377        double xold, yold, zold, x, y, z;
    379378
    380         if (mjd != lastmjd) {
     379        if (mj != lastmj) {
    381380            double epsilon, dpsi, deps;
    382381            double se, ce, sp, cp, sede, cede;
    383382
    384             obliquity(mjd, &epsilon);
    385             nutation(mjd, &deps, &dpsi);
     383            obliquity(mj, &epsilon);
     384            nutation(mj, &deps, &dpsi);
    386385
    387386            /* the rotation matrix a applies the nutation correction to
    … …  
    428427            a[2][2] = sede*cp*se+cede*ce;
    429428
    430             lastmjd = mjd;
     429            lastmj = mj;
    431430        }
    432431
    … …  
    440439
    441440/* For RCS Only -- Do Not Edit */
    442 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: nutation.c,v $ $Date: 2001-10-22 12:08:27 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     441static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: nutation.c,v $ $Date: 2004-06-15 16:52:39 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/obliq.c

    r1719 r2551  
    11#include <stdio.h>
    22
    3 #include "P_.h"
    43#include "astro.h"
    54
    6 /* given the modified Julian date, mjd, find the mean obliquity of the
     5/* given the modified Julian date, mj, find the mean obliquity of the
    76 * ecliptic, *eps, in radians.
    87 *
    … …  
    109 */
    1110void
    12 obliquity (mjd, eps)
    13 double mjd;
    14 double *eps;
     11obliquity (double mj, double *eps)
    1512{
    16         static double lastmjd = -16347, lasteps;
     13        static double lastmj = -16347, lasteps;
    1714
    18         if (mjd != lastmjd) {
    19             double t = (mjd - J2000)/36525.;    /* centuries from J2000 */
     15        if (mj != lastmj) {
     16            double t = (mj - J2000)/36525.;     /* centuries from J2000 */
    2017            lasteps = degrad(23.4392911 +       /* 23^ 26' 21".448 */
    2118                            t * (-46.8150 +
    2219                            t * ( -0.00059 +
    2320                            t * (  0.001813 )))/3600.0);
    24             lastmjd = mjd;
     21            lastmj = mj;
    2522        }
    2623        *eps = lasteps;
    … …  
    2825
    2926/* For RCS Only -- Do Not Edit */
    30 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: obliq.c,v $ $Date: 2001-10-22 12:08:27 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     27static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: obliq.c,v $ $Date: 2004-06-15 16:52:39 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/parallax.c

    r1719 r2551  
    22#include <math.h>
    33
    4 #include "P_.h"
    54#include "astro.h"
    65
    … …  
    1413 */
    1514void
    16 ta_par (tha, tdec, phi, ht, rho, aha, adec)
    17 double tha, tdec, phi, ht, *rho;
    18 double *aha, *adec;
     15ta_par (double tha, double tdec, double phi, double ht, double *rho,
     16double *aha, double *adec)
    1917{
    2018        static double last_phi = 1000.0, last_ht = -1000.0, xobs, zobs;
    … …  
    4240
    4341/* For RCS Only -- Do Not Edit */
    44 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: parallax.c,v $ $Date: 2001-10-22 12:08:27 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     42static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: parallax.c,v $ $Date: 2004-06-15 16:52:39 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/plans.c

    r1719 r2551  
    55#include <math.h>
    66
    7 #include "P_.h"
    87#include "astro.h"
    98#include "vsop87.h"
    109#include "chap95.h"
    1110
    12 static void pluto_ell P_((double mjd, double *ret));
    13 static void chap_trans P_((double mjd, double *ret));
    14 static void planpos P_((double mjd, int obj, double prec, double *ret));
     11static void pluto_ell (double mj, double *ret);
     12static void chap_trans (double mj, double *ret);
     13static void planpos (double mj, int obj, double prec, double *ret);
    1514
    1615/* coordinate transformation
    … …  
    2120 */
    2221static void
    23 chap_trans (mjd, ret)
    24 double mjd;     /* destination epoch */
    25 double *ret;    /* vector to be transformed _IN PLACE_ */
     22chap_trans (
     23double mj,      /* destination epoch */
     24double *ret)    /* vector to be transformed _IN PLACE_ */
    2625{
    2726        double ra, dec, r, eps;
    … …  
    2928
    3029        cartsph(ret[0], ret[1], ret[2], &ra, &dec, &r);
    31         precess(J2000, mjd, &ra, &dec);
    32         obliquity(mjd, &eps);
     30        precess(J2000, mj, &ra, &dec);
     31        obliquity(mj, &eps);
    3332        sr = sin(ra); cr = cos(ra);
    3433        sd = sin(dec); cd = cos(dec);
    … …  
    4342 */
    4443static void
    45 pluto_ell (mjd, ret)
    46 double mjd;     /* epoch */
    47 double *ret;    /* ecliptic coordinates {l,b,r} at equinox of date */
     44pluto_ell (
     45double mj,      /* epoch */
     46double *ret)    /* ecliptic coordinates {l,b,r} at equinox of date */
    4847{
    4948        /* mean orbital elements of Pluto.
    … …  
    5554                Om0 = 110.307,                  /* long asc node, deg */
    5655                omeg0 = 113.768,                /* arg of perihel, deg */
    57                 mjdp = 2448045.539 - MJD0,      /* epoch of perihel */
    58                 mjdeq = J2000,                  /* equinox of elements */
     56                mjp = 2448045.539 - MJD0,       /* epoch of perihel */
     57                mjeq = J2000,                   /* equinox of elements */
    5958                n = 144.9600/36525.;            /* daily motion, deg */
    6059
    … …  
    6362        double lo, slo, clo;    /* longitude in orbit from asc node */
    6463
    65         reduce_elements(mjdeq, mjd, degrad(inc0), degrad(omeg0), degrad(Om0),
     64        reduce_elements(mjeq, mj, degrad(inc0), degrad(omeg0), degrad(Om0),
    6665                                &inc, &omeg, &Om);
    67         ma = degrad((mjd - mjdp) * n);
     66        ma = degrad((mj - mjp) * n);
    6867        anomaly(ma, e, &nu, &ea);
    6968        ret[2] = a * (1.0 - e*cos(ea));                 /* r */
    … …  
    8180 */
    8281static void
    83 planpos (mjd, obj, prec, ret)
    84 double mjd;
    85 int obj;
    86 double prec;
    87 double *ret;
    88 {
    89         if (mjd >= CHAP_BEGIN && mjd <= CHAP_END) {
     82planpos (double mj, int obj, double prec, double *ret)
     83{
     84        if (mj >= CHAP_BEGIN && mj <= CHAP_END) {
    9085            if (obj >= JUPITER) {               /* prefer Chapront */
    91                 chap95(mjd, obj, prec, ret);
    92                 chap_trans (mjd, ret);
     86                chap95(mj, obj, prec, ret);
     87                chap_trans (mj, ret);
    9388            } else {                            /* VSOP for inner planets */
    94                 vsop87(mjd, obj, prec, ret);
     89                vsop87(mj, obj, prec, ret);
    9590            }
    9691        } else {                                /* outside Chapront time: */
    9792            if (obj != PLUTO) {                 /* VSOP for all but Pluto */
    98                 vsop87(mjd, obj, prec, ret);
     93                vsop87(mj, obj, prec, ret);
    9994            } else {                            /* Pluto mean elliptic orbit */
    100                 pluto_ell(mjd, ret);
     95                pluto_ell(mj, ret);
    10196            }
    10297        }
    … …  
    126121};
    127122
    128 /* given a modified Julian date, mjd, and a planet, p, find:
     123/* given a modified Julian date, mj, and a planet, p, find:
    129124 *   lpd0: heliocentric longitude,
    130125 *   psi0: heliocentric latitude,
    … …  
    148143 */
    149144void
    150 plans (mjd, p, lpd0, psi0, rp0, rho0, lam, bet, dia, mag)
    151 double mjd;
    152 int p;
    153 double *lpd0, *psi0, *rp0, *rho0, *lam, *bet, *dia, *mag;
    154 {
    155         static double lastmjd = -10000;
     145plans (double mj, PLCode p, double *lpd0, double *psi0, double *rp0,
     146double *rho0, double *lam, double *bet, double *dia, double *mag)
     147{
     148        static double lastmj = -10000;
    156149        static double lsn, bsn, rsn;    /* geocentric coords of sun */
    157150        static double xsn, ysn, zsn;    /* cartesian " */
    … …  
    163156        int pass;
    164157
    165         /* get sun cartesian; needed only once at mjd */
    166         if (mjd != lastmjd) {
    167             sunpos (mjd, &lsn, &rsn, &bsn);
     158        /* get sun cartesian; needed only once at mj */
     159        if (mj != lastmj) {
     160            sunpos (mj, &lsn, &rsn, &bsn);
    168161            sphcart (lsn, bsn, rsn, &xsn, &ysn, &zsn);
    169             lastmjd = mjd;
     162            lastmj = mj;
    170163        }
    171164
    172         /* first find the true position of the planet at mjd.
     165        /* first find the true position of the planet at mj.
    173166         * then repeat a second time for a slightly different time based
    174167         * on the position found in the first pass to account for light-travel
    … …  
    183176             * alternative option:  vsop allows calculating rates.
    184177             */
    185             planpos(mjd - dt, p, 0.0, ret);
     178            planpos(mj - dt, p, 0.0, ret);
    186179
    187180            lp = ret[0];
    … …  
    225218        if (p == SATURN) {
    226219            double et, st, set;
    227             satrings (bp, lp, rp, lsn+PI, rsn, mjd+MJD0, &et, &st);
     220            satrings (bp, lp, rp, lsn+PI, rsn, mj+MJD0, &et, &st);
    228221            set = sin(fabs(et));
    229222            *mag += (-2.60 + 1.25*set)*set;
    … …  
    232225
    233226/* For RCS Only -- Do Not Edit */
    234 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: plans.c,v $ $Date: 2001-10-22 12:08:27 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     227static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: plans.c,v $ $Date: 2004-06-15 16:52:39 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/precess.c

    r1719 r2551  
    22#include <math.h>
    33
    4 #include "P_.h"
    54#include "astro.h"
    65
    7 static void precess_hiprec P_((double mjd1, double mjd2, double *ra,
    8     double *dec));
     6static void precess_hiprec (double mjd1, double mjd2, double *ra, double *dec);
    97
    108
    … …  
    1917 */
    2018void
    21 precess (mjd1, mjd2, ra, dec)
    22 double mjd1, mjd2;      /* initial and final epoch modified JDs */
    23 double *ra, *dec;       /* ra/dec for mjd1 in, for mjd2 out */
     19precess (
     20double mjd1, double mjd2,       /* initial and final epoch modified JDs */
     21double *ra, double *dec)        /* ra/dec for mjd1 in, for mjd2 out */
    2422{
    2523        precess_hiprec (mjd1, mjd2, ra, dec);
    … …  
    4240 */
    4341static void
    44 precess_hiprec (mjd1, mjd2, ra, dec)
    45 double mjd1, mjd2;      /* initial and final epoch modified JDs */
    46 double *ra, *dec;       /* ra/dec for mjd1 in, for mjd2 out */
     42precess_hiprec (
     43double mjd1, double mjd2,       /* initial and final epoch modified JDs */
     44double *ra, double *dec)        /* ra/dec for mjd1 in, for mjd2 out */
    4745{
    4846        static double last_mjd1 = -213.432, last_from;
    … …  
    130128#if 0
    131129static void
    132 precess_fast (mjd1, mjd2, ra, dec)
    133 double mjd1, mjd2;      /* initial and final epoch modified JDs */
    134 double *ra, *dec;       /* ra/dec for mjd1 in, for mjd2 out */
     130precess_fast (
     131double mjd1, double mjd2,       /* initial and final epoch modified JDs */
     132double *ra, double *dec)        /* ra/dec for mjd1 in, for mjd2 out */
    135133{
    136134#define N       degrad (20.0468/3600.0)
    … …  
    146144
    147145/* For RCS Only -- Do Not Edit */
    148 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: precess.c,v $ $Date: 2001-10-22 12:08:27 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     146static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: precess.c,v $ $Date: 2004-06-15 16:52:39 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/preferences.h

    r1719 r2551  
    2323typedef enum {PREF_SAT, PREF_SUN, PREF_MON} PrefWeekStart;
    2424
    25 extern int pref_get P_((Preferences p));
    26 extern int pref_set P_((Preferences p, int new));
     25extern int pref_get (Preferences p);
     26extern int pref_set (Preferences p, int newp);
    2727
    2828#endif /* _PREFERENCES_H */
    2929
    3030/* For RCS Only -- Do Not Edit
    31  * @(#) $RCSfile: preferences.h,v $ $Date: 2001-10-22 12:08:27 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $
     31 * @(#) $RCSfile: preferences.h,v $ $Date: 2004-06-15 16:52:39 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $
    3232 */

  • trunk/SophyaExt/XephemAstroLib/reduce.c

    r1719 r2551  
    11#include <math.h>
    22
    3 #include "P_.h"
     3#include <stdio.h>
     4
    45#include "astro.h"
    56
    6 /* convert those orbital elements that change from epoch mjd0 to epoch mjd.
     7/* convert those orbital elements that change from epoch mj0 to epoch mj.
    78 */
    89void
    9 reduce_elements (mjd0, mjd, inc0, ap0, om0, inc, ap, om)
    10 double mjd0;    /* initial epoch */
    11 double mjd;     /* desired epoch */
    12 double inc0;    /* initial inclination, rads */
    13 double ap0;     /* initial argument of perihelion, as an mjd */
    14 double om0;     /* initial long of ascending node, rads */
    15 double *inc;    /* desired inclination, rads */
    16 double *ap;     /* desired epoch of perihelion, as an mjd */
    17 double *om;     /* desired long of ascending node, rads */
     10reduce_elements (
     11double mj0,     /* initial epoch */
     12double mj,      /* desired epoch */
     13double inc0,    /* initial inclination, rads */
     14double ap0,     /* initial argument of perihelion, as an mj */
     15double om0,     /* initial long of ascending node, rads */
     16double *inc,    /* resultant inclination, rads */
     17double *ap,     /* resultant arg of perihelion, as an mj */
     18double *om)     /* resultant long of ascending node, rads */
    1819{
    1920        double t0, t1;
    … …  
    2627        double seta, ceta;
    2728
    28         if (fabs(mjd - mjd0) < 1e-5) {
     29        if (fabs(mj - mj0) < 1e-5) {
    2930            /* sin(eta) blows for inc < 10 degrees -- anyway, no need */
    3031            *inc = inc0;
    … …  
    3435        }
    3536
    36         t0 = mjd0/365250.0;
    37         t1 = mjd/365250.0;
     37        t0 = mj0/365250.0;
     38        t1 = mj/365250.0;
    3839
    3940        tt = t1-t0;
    … …  
    7576
    7677/* For RCS Only -- Do Not Edit */
    77 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: reduce.c,v $ $Date: 2001-10-22 12:08:27 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     78static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: reduce.c,v $ $Date: 2004-06-15 16:52:40 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/refract.c

    r1895 r2551  
    22#include <math.h>
    33
    4 #include "P_.h"
    54#include "astro.h"
    65
    7 static void unrefractLT15 P_((double pr, double tr, double aa, double *ta));
    8 static void unrefractGE15 P_((double pr, double tr, double aa, double *ta));
     6static void unrefractLT15 (double pr, double tr, double aa, double *ta);
     7static void unrefractGE15 (double pr, double tr, double aa, double *ta);
    98
    109void
    11 unrefract (pr, tr, aa, ta)
    12 double pr, tr;
    13 double aa;
    14 double *ta;
     10unrefract (double pr, double tr, double aa, double *ta)
    1511{
    1612#define LTLIM   14.5
    … …  
    3531
    3632static void
    37 unrefractGE15 (pr, tr, aa, ta)
    38 double pr, tr;
    39 double aa;
    40 double *ta;
     33unrefractGE15 (double pr, double tr, double aa, double *ta)
    4134{
    4235        double r;
    … …  
    4740
    4841static void
    49 unrefractLT15 (pr, tr, aa, ta)
    50 double pr, tr;
    51 double aa;
    52 double *ta;
     42unrefractLT15 (double pr, double tr, double aa, double *ta)
    5343{
    5444        double aadeg = raddeg(aa);
    … …  
    6757 */
    6858void
    69 refract (pr, tr, ta, aa)
    70 double pr, tr;
    71 double ta;
    72 double *aa;
     59refract (double pr, double tr, double ta, double *aa)
    7360{
    7461#define MAXRERR degrad(0.1/3600.)       /* desired accuracy, rads */
    … …  
    10289
    10390/* For RCS Only -- Do Not Edit */
    104 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: refract.c,v $ $Date: 2002-02-07 09:27:06 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};
     91static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: refract.c,v $ $Date: 2004-06-15 16:52:40 $ $Revision: 1.4 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/rings.c

    r1719 r2551  
    11#include <stdio.h>
    22#include <math.h>
    3 #if defined(__STDC__)
    43#include <stdlib.h>
    5 #endif
    64
    7 #include "P_.h"
    85#include "astro.h"
    9 #include "circum.h"
    106
    117/*  RINGS OF SATURN by Olson, et al, BASIC Code from Sky & Telescope, May 1995.
    … …  
    1410 */
    1511void
    16 satrings (sb, sl, sr, el, er, JD, etiltp, stiltp)
    17 double sb, sl, sr;              /* Saturn hlat, hlong, sun dist */
    18 double el, er;                  /* Earth hlong, sun dist */
    19 double JD;                      /* Julian date */
    20 double *etiltp, *stiltp;        /* tilt from earth and sun, rads southward */
     12satrings (
     13double sb, double sl, double sr,        /* Saturn hlat, hlong, sun dist */
     14double el, double er,                   /* Earth hlong, sun dist */
     15double JD,                              /* Julian date */
     16double *etiltp, double *stiltp)         /* tilt from earth and sun, rads south*/
    2117{
    2218        double t, i, om;
    … …  
    4541        *stiltp = bp;
    4642}
     43
     44/* For RCS Only -- Do Not Edit */
     45static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: rings.c,v $ $Date: 2004-06-15 16:52:40 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/riset.c

    r1719 r2551  
    22#include <math.h>
    33
    4 #include "P_.h"
    54#include "astro.h"
    65
    76/* given the true geocentric ra and dec of an object, the observer's latitude,
    8  *   lat, and a horizon displacement correction, dis, all in radians, find the
     7 *   lt, and a horizon displacement correction, dis, all in radians, find the
    98 *   local sidereal times and azimuths of rising and setting, lstr/s
    109 *   and azr/s, also all in radians, respectively.
    … …  
    2524 */
    2625void
    27 riset (ra, dec, lat, dis, lstr, lsts, azr, azs, status)
    28 double ra, dec;
    29 double lat, dis;
    30 double *lstr, *lsts;
    31 double *azr, *azs;
    32 int *status;
     26riset (double ra, double dec, double lt, double dis, double *lstr,
     27double *lsts, double *azr, double *azs, int *status)
    3328{
    3429#define EPS     (1e-9)  /* math rounding fudge - always the way, eh? */
    … …  
    4035        int shemi;              /* flag for southern hemisphere reflection */
    4136
    42         /* reflect lat and dec if in southern hemisphere, then az back later */
    43         if ((shemi= (lat < 0.)) != 0) {
    44             lat = -lat;
     37        /* reflect lt and dec if in southern hemisphere, then az back later */
     38        if ((shemi= (lt < 0.)) != 0) {
     39            lt = -lt;
    4540            dec = -dec;
    4641        }
    … …  
    4843        /* establish zenith angle, and its extrema */
    4944        z = (PI/2.) + dis;
    50         zmin = fabs (dec - lat);
    51         zmax = PI - fabs(dec + lat);
     45        zmin = fabs (dec - lt);
     46        zmax = PI - fabs(dec + lt);
    5247
    5348        /* first consider special cases.
    … …  
    6459
    6560        /* compute rising hour angle -- beware found off */
    66         cos_h = (cos(z)-sin(lat)*sin(dec))/(cos(lat)*cos(dec));
     61        cos_h = (cos(z)-sin(lt)*sin(dec))/(cos(lt)*cos(dec));
    6762        if (cos_h >= 1.)
    6863            h =  0.;
    … …  
    7368
    7469        /* compute setting azimuth -- beware found off */
    75         xaz = sin(dec)*cos(lat)-cos(dec)*cos(h)*sin(lat);
     70        xaz = sin(dec)*cos(lt)-cos(dec)*cos(h)*sin(lt);
    7671        yaz = -1.*cos(dec)*sin(h);
    7772        if (xaz == 0.) {
    … …  
    10398
    10499/* For RCS Only -- Do Not Edit */
    105 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: riset.c,v $ $Date: 2001-10-22 12:08:28 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     100static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: riset.c,v $ $Date: 2004-06-15 16:52:40 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/riset_cir.c

    r1719 r2551  
    33#include <stdio.h>
    44#include <math.h>
    5 #if defined(__STDC__)
    65#include <stdlib.h>
    76#include <string.h>
    8 #endif
    9 
    10 #include "P_.h"
     7
    118#include "astro.h"
    12 #include "circum.h"
    139
    1410#define TMACC   (10./3600./24.0)        /* convergence accuracy, days */
    1511
    16 static void e_riset_cir P_((Now *np, Obj *op, double dis, RiseSet *rp));
    17 static int find_0alt P_((double dt, double dis, Now *np, Obj *op));
    18 static int find_transit P_((double dt, Now *np, Obj *op));
    19 static int find_max P_((Now *np, Obj *op, double tr, double ts, double *tp,
    20     double *alp));
     12static void e_riset_cir (Now *np, Obj *op, double dis, RiseSet *rp);
     13static int find_0alt (double dt, double dis, Now *np, Obj *op);
     14static int find_transit (double dt, Now *np, Obj *op);
     15static int find_max (Now *np, Obj *op, double tr, double ts, double *tp,
     16    double *alp);
    2117
    2218/* find where and when an object, op, will rise and set and
    … …  
    2622 */
    2723void
    28 riset_cir (np, op, dis, rp)
    29 Now *np;
    30 Obj *op;
    31 double dis;
    32 RiseSet *rp;
     24riset_cir (Now *np, Obj *op, double dis, RiseSet *rp)
    3325{
    3426        double mjdn;    /* mjd of local noon */
    … …  
    6456        /* first approximation is to find rise/set times of a fixed object
    6557         * at the current epoch in its position at local noon.
    66          * N.B. add typical refraction for initial go/no-go test. if it
    67          *   passes, real code does refraction rigorously.
     58         * N.B. add typical refraction if dis is above horizon for initial
     59         *   go/no-go test. if it passes, real code does refraction rigorously.
    6860         */
    6961        n.n_mjd = mjdn;
    … …  
    7365        }
    7466        ran = o.s_gaera;
    75         riset (o.s_gaera, o.s_gaedec, lat, dis+.01, &lr, &ls, &ar, &as, &rss);
     67        riset (o.s_gaera, o.s_gaedec, lat, dis+(dis>.01 ? 0 : .01), &lr, &ls,
     68                                                                &ar, &as, &rss);
    7669        switch (rss) {
    7770        case  0:  break;
    … …  
    133126 */
    134127void
    135 twilight_cir (np, dis, dawn, dusk, status)
    136 Now *np;
    137 double dis;
    138 double *dawn, *dusk;
    139 int *status;
     128twilight_cir (Now *np, double dis, double *dawn, double *dusk, int *status)
    140129{
    141130        RiseSet rs;
    142131        Obj o;
    143132
     133        memset (&o, 0, sizeof(o));
    144134        o.o_type = PLANET;
    145         o.pl.pl_code = SUN;
     135        o.pl_code = SUN;
    146136        (void) strcpy (o.o_name, "Sun");
    147137        riset_cir (np, &o, dis, &rs);
    … …  
    160150 */
    161151static void
    162 e_riset_cir (np, op, dis, rp)
    163 Now *np;
    164 Obj *op;
    165 double dis;
    166 RiseSet *rp;
     152e_riset_cir (Now *np, Obj *op, double dis, RiseSet *rp)
    167153{
    168154#define DEGSTEP 5               /* time step is about this many degrees */
    … …  
    272258 */
    273259static int
    274 find_0alt (dt, dis, np, op)
    275 double dt;      /* hours from noon to first guess at event */
    276 double dis;     /* horizon displacement, rads */
    277 Now *np;        /* working Now -- starts with mjd is noon, returns as answer */
    278 Obj *op;        /* working object -- returns as answer */
     260find_0alt (
     261double dt,      /* hours from noon to first guess at event */
     262double dis,     /* horizon displacement, rads */
     263Now *np,        /* working Now -- starts with mjd is noon, returns as answer */
     264Obj *op)        /* working object -- returns as answer */
    279265{
    280266#define MAXPASSES       20              /* max iterations to try */
    … …  
    286272
    287273        /* insure initial guess is today -- if not, move by 24 hours */
    288         if (dt < -12.0)
    289             dt += 24.0;
    290         if (dt > 12.0)
    291             dt -= 24.0;
     274        if (dt < -12.0 && !find_0alt (dt+24, dis, np, op))
     275            return (0);
     276        mjd = mjdn;
     277        if (dt > 12.0 && !find_0alt (dt-24, dis, np, op))
     278            return (0);
     279        mjd = mjdn;
    292280       
    293281        /* convert dt to days for remainder of algorithm */
    … …  
    325313 */
    326314static int
    327 find_transit (dt, np, op)
    328 double dt;
    329 Now *np;
    330 Obj *op;
     315find_transit (double dt, Now *np, Obj *op)
    331316{
    332317#define MAXLOOPS        10
    … …  
    370355 */
    371356static int
    372 find_max (np, op, tr, ts, tp, alp)
    373 Now *np;
    374 Obj *op;
    375 double tr, ts;          /* times of rise and set */
    376 double *tp, *alp;       /* time of max altitude, and that altitude */
     357find_max (
     358Now *np,
     359Obj *op,
     360double tr, double ts,           /* times of rise and set */
     361double *tp, double *alp)        /* time of max altitude, and that altitude */
    377362{
    378363        mjd = (ts + tr)/2;
    … …  
    385370
    386371/* For RCS Only -- Do Not Edit */
    387 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: riset_cir.c,v $ $Date: 2001-10-22 12:08:28 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     372static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: riset_cir.c,v $ $Date: 2004-06-15 16:52:40 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/satlib.h

    r1719 r2551  
    22#define __SATLIB_H
    33
    4 /* $Id: satlib.h,v 1.2 2001-10-22 12:08:28 cmv Exp $ */
     4/* $Id: satlib.h,v 1.3 2004-06-15 16:52:40 cmv Exp $ */
    55
    66typedef struct _SatElem {
    … …  
    203203
    204204/* For RCS Only -- Do Not Edit
    205  * @(#) $RCSfile: satlib.h,v $ $Date: 2001-10-22 12:08:28 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $
     205 * @(#) $RCSfile: satlib.h,v $ $Date: 2004-06-15 16:52:40 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $
    206206 */

  • trunk/SophyaExt/XephemAstroLib/satspec.h

    r1719 r2551  
    22#define __SATSPEC_H
    33
    4 /* $Id: satspec.h,v 1.2 2001-10-22 12:08:28 cmv Exp $ */
     4/* $Id: satspec.h,v 1.3 2004-06-15 16:52:40 cmv Exp $ */
    55
    66#include "sattypes.h"
    … …  
    4040
    4141/* For RCS Only -- Do Not Edit
    42  * @(#) $RCSfile: satspec.h,v $ $Date: 2001-10-22 12:08:28 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $
     42 * @(#) $RCSfile: satspec.h,v $ $Date: 2004-06-15 16:52:40 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $
    4343 */

  • trunk/SophyaExt/XephemAstroLib/sattypes.h

    r1719 r2551  
    22#define __SATTYPES_H
    33
    4 /* $Id: sattypes.h,v 1.2 2001-10-22 12:08:28 cmv Exp $ */
     4/* $Id: sattypes.h,v 1.3 2004-06-15 16:52:40 cmv Exp $ */
    55
    66typedef struct _Vec3 {
    … …  
    2525
    2626/* For RCS Only -- Do Not Edit
    27  * @(#) $RCSfile: sattypes.h,v $ $Date: 2001-10-22 12:08:28 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $
     27 * @(#) $RCSfile: sattypes.h,v $ $Date: 2004-06-15 16:52:40 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $
    2828 */

  • trunk/SophyaExt/XephemAstroLib/sdp4.c

    r1719 r2551  
    9696
    9797void
    98 sdp4 (sat, pos, dpos, TSINCE)
    99 SatData *sat;
    100 Vec3 *pos;
    101 Vec3 *dpos;
    102 double TSINCE;
     98sdp4 (SatData *sat, Vec3 *pos, Vec3 *dpos, double TSINCE)
    10399{
    104100    int i;
    … …  
    109105
    110106    /* private temporary variables */
    111     double A,AXN,AYN,AYNL,BETA,BETAL,CAPU,COS2U,COSEPW,
     107    double A,AXN,AYN,AYNL,BETA,BETAL,CAPU,COS2U,COSEPW=0,
    112108        COSIK,COSNOK,COSU,COSUK,E,ECOSE,ELSQ,EM=0,EPW,ESINE,OMGADF,PL,
    113         R,RDOT,RDOTK,RFDOT,RFDOTK,RK,SIN2U,SINEPW,SINIK,SINNOK,
    114         SINU,SINUK,TEMP,TEMP1,TEMP2,TEMP3,TEMP4,TEMP5,TEMP6,TEMPA,
     109        R,RDOT,RDOTK,RFDOT,RFDOTK,RK,SIN2U,SINEPW=0,SINIK,SINNOK,
     110        SINU,SINUK,TEMP,TEMP1,TEMP2,TEMP3=0,TEMP4=0,TEMP5=0,TEMP6=0,TEMPA,
    115111        TEMPE,TEMPL,TSQ,U,UK,UX,UY,UZ,VX,VY,VZ,XINC=0,XINCK,XL,XLL,XLT,
    116112        XMAM,XMDF,XMX,XMY,XN,XNODDF,XNODE,XNODEK;
    … …  
    432428
    433429/* For RCS Only -- Do Not Edit */
    434 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: sdp4.c,v $ $Date: 2001-10-22 12:08:28 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     430static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: sdp4.c,v $ $Date: 2004-06-15 16:52:40 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/sgp4.c

    r1719 r2551  
    7979 */
    8080void
    81 sgp4(sat, pos, dpos, TSINCE)
    82 SatData *sat;
    83 Vec3 *pos;
    84 Vec3 *dpos;
    85 double TSINCE;
     81sgp4(SatData *sat, Vec3 *pos, Vec3 *dpos, double TSINCE)
    8682{
    8783    int i;
    … …  
    8985    double A1, A3OVK2, AO, BETAO, BETAO2, C1SQ, C2, C3, COEF, COEF1,
    9086        DEL1, DELO, EETA, EOSQ, ETASQ, PERIGE, PINVSQ, PSISQ, QOMS24,
    91         S4, TEMP, TEMP1, TEMP2, TEMP3, THETA2, THETA4, TSI, X1M5TH,
     87        S4, TEMP, TEMP1, TEMP2, TEMP3=0, THETA2, THETA4, TSI, X1M5TH,
    9288        XHDOT1;
    9389
    94     double A, AXN, AYN, AYNL, BETA, BETAL, CAPU, COS2U, COSEPW, COSIK,
     90    double A, AXN, AYN, AYNL, BETA, BETAL, CAPU, COS2U, COSEPW=0, COSIK,
    9591        COSNOK, COSU, COSUK, DELM, DELOMG, E, ECOSE, ELSQ, EPW, ESINE,
    9692        OMEGA, OMGADF, PL, R, RDOT, RDOTK, RFDOT, RFDOTK, RK, SIN2U,
    97         SINEPW, SINIK, SINNOK, SINU, SINUK, TCUBE, TEMP4, TEMP5, TEMP6,
     93        SINEPW=0, SINIK, SINNOK, SINU, SINUK, TCUBE, TEMP4=0, TEMP5=0, TEMP6=0,
    9894        TEMPA, TEMPE, TEMPL, TFOUR, TSQ, U, UK, UX, UY, UZ, VX, VY, VZ,
    9995        XINCK, XL, XLL, XLT, XMDF, XMP, XMX, XMY, XN, XNODDF, XNODE,
    … …  
    403399
    404400/* For RCS Only -- Do Not Edit */
    405 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: sgp4.c,v $ $Date: 2001-10-22 12:08:28 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     401static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: sgp4.c,v $ $Date: 2004-06-15 16:52:40 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/sphcart.c

    r1719 r2551  
    11#include <math.h>
    2 #include "P_.h"
     2#include <stdio.h>
     3
    34#include "astro.h"
    45
    56/* transformation from spherical to cartesian coordinates */
    67void
    7 sphcart (l, b, r, x, y, z)
    8 double l, b, r;         /* source: spherical coordinates */
    9 double *x, *y, *z;      /* result: rectangular coordinates */
     8sphcart (
     9double l, double b, double r,           /* source: spherical coordinates */
     10double *x, double *y, double *z)        /* result: rectangular coordinates */
    1011{
    1112        double rcb = r * cos(b);
    … …  
    1819/* transformation from cartesian to spherical coordinates */
    1920void
    20 cartsph (x, y, z, l, b, r)
    21 double x, y, z;         /* source: rectangular coordinates */
    22 double *l, *b, *r;      /* result: spherical coordinates */
     21cartsph (
     22double x, double y, double z,           /* source: rectangular coordinates */
     23double *l, double *b, double *r)        /* result: spherical coordinates */
    2324{
    2425        double rho = x*x + y*y;
    … …  
    4041
    4142/* For RCS Only -- Do Not Edit */
    42 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: sphcart.c,v $ $Date: 2001-10-22 12:08:28 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     43static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: sphcart.c,v $ $Date: 2004-06-15 16:52:40 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/sun.c

    r1719 r2551  
    22#include <math.h>
    33
    4 #include "P_.h"
    54#include "astro.h"
    65#include "vsop87.h"
    76
    8 /* given the modified JD, mjd, return the true geocentric ecliptic longitude
     7/* given the modified JD, mj, return the true geocentric ecliptic longitude
    98 *   of the sun for the mean equinox of the date, *lsn, in radians, the
    109 *   sun-earth distance, *rsn, in AU, and the latitude *bsn, in radians
    … …  
    1716 */
    1817void
    19 sunpos (mjd, lsn, rsn, bsn)
    20 double mjd;
    21 double *lsn, *rsn, *bsn;
     18sunpos (double mj, double *lsn, double *rsn, double *bsn)
    2219{
    23         static double last_mjd = -3691, last_lsn, last_rsn, last_bsn;
     20        static double last_mj = -3691, last_lsn, last_rsn, last_bsn;
    2421        double ret[6];
    2522
    26         if (mjd == last_mjd) {
     23        if (mj == last_mj) {
    2724            *lsn = last_lsn;
    2825            *rsn = last_rsn;
    … …  
    3128        }
    3229
    33         vsop87(mjd, SUN, 0.0, ret);     /* full precision earth pos */
     30        vsop87(mj, SUN, 0.0, ret);      /* full precision earth pos */
    3431
    3532        *lsn = ret[0] - PI;             /* revert to sun pos */
    … …  
    3936        last_rsn = *rsn = ret[2];
    4037        last_bsn = -ret[1];
    41         last_mjd = mjd;
     38        last_mj = mj;
    4239
    4340        if (bsn) *bsn = last_bsn;       /* assign only if non-NULL pointer */
    … …  
    4542
    4643/* For RCS Only -- Do Not Edit */
    47 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: sun.c,v $ $Date: 2001-10-22 12:08:28 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     44static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: sun.c,v $ $Date: 2004-06-15 16:52:40 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/testxephem.c

    r2395 r2551  
    44#include <math.h>
    55#include <unistd.h>
    6 #include "P_.h"
     6/* include "P_.h" */
    77#include "astro.h"
    88
    … …  
    1111int main(int narg,char** arg)
    1212{
    13   int year=1989,month=7,day=21; /* UTC*/
     13 int year=1989,month=7,day=21; /* UTC*/
    1414 double hour=12.5;
    1515 int npnum=8, pnum[8]={MERCURY,VENUS,MARS,JUPITER,SATURN,URANUS,NEPTUNE,PLUTO};
    1616 int i;
    17  double mjd,sunecl,sunecb,sundist,geodist,geoecl,geoecb,diamang,mag,msp,mdp;
     17 double mymjd,sunecl,sunecb,sundist,geodist,geoecl,geoecb,diamang,mag,msp,mdp;
    1818
    19  cal_mjd(month,day+hour/24.,year,&mjd);
    20  printf("date %d/%d/%d %gh  -> mjd=%f\n",day,month,year,hour,mjd);
     19 cal_mjd(month,day+hour/24.,year,&mymjd);
     20 printf("date %d/%d/%d %gh  -> mjd=%f\n",day,month,year,hour,mymjd);
    2121
    2222 /* Soleil */
    23  sunpos(mjd,&geoecl,&geodist,&geoecb);
     23 sunpos(mymjd,&geoecl,&geodist,&geoecb);
    2424 printf("Sun: geo ecl=%g ecb=%g dist=%g\n",geoecl*R2D,geoecb*R2D,geodist);
    2525
    2626 /* Lune */
    27  moon(mjd,&geoecl,&geoecb,&geodist,&msp,&mdp);
     27 moon(mymjd,&geoecl,&geoecb,&geodist,&msp,&mdp);
    2828 printf("Moon: geo ecl=%g ecb=%g dist=%g\n",geoecl*R2D,geoecb*R2D,geodist);
    2929
    … …  
    3131 printf("--- Planets\n");
    3232 for(i=0;i<npnum;i++) {
    33    plans(mjd,pnum[i],&sunecl,&sunecb,&sundist,&geodist,&geoecl,&geoecb,&diamang,&mag);
     33   plans(mymjd,pnum[i],&sunecl,&sunecb,&sundist,&geodist,&geoecl,&geoecb,&diamang,&mag);
    3434   printf("pnum=%d: sun ecl=%g ecb=%g dist=%g,",pnum[i],sunecl*R2D,sunecb*R2D,sundist);
    3535   printf(" geo ecl=%g ecb=%g dist=%g,",geoecl*R2D,geoecb*R2D,geodist);

  • trunk/SophyaExt/XephemAstroLib/thetag.c

    r1719 r2551  
    33#include "deepconst.h"
    44
    5 /* @(#) $Id: thetag.c,v 1.2 2001-10-22 12:08:28 cmv Exp $ */
     5/* @(#) $Id: thetag.c,v 1.3 2004-06-15 16:52:40 cmv Exp $ */
    66
    77
    … …  
    8888
    8989/* For RCS Only -- Do Not Edit */
    90 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: thetag.c,v $ $Date: 2001-10-22 12:08:28 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     90static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: thetag.c,v $ $Date: 2004-06-15 16:52:40 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/twobody.c

    r1719 r2551  
    2222 
    2323/* Constants used when solving Kepler's equation */
     24#undef  EPSILON
    2425#define EPSILON   3E-8
     26#undef  INFINITY
    2527#define INFINITY  1E+10
    2628 
    2729/* Math constants */
     30#undef  PI
    2831#define PI      3.14159265358979323846
    2932#define RADEG   ( 180.0 / PI )
    … …  
    125128 
    126129 
    127 void vrc( double *v, double *r, double tp, double e, double q )
     130/* return 0 if ok, else -1 */
     131int vrc( double *v, double *r, double tp, double e, double q )
    128132/*
    129133 * Elliptic, hyperbolic and near-parabolic orbits:
    … …  
    145149        *v = 0.0;
    146150        *r = q;
    147         return;
     151        return 0;
    148152    }
    149153 
    … …  
    176180                printf( "\nNear-parabolic orbit: inaccurate result."
    177181                        "\n  e = %f, lambda = %f, w = %f", e, lambda, w );
    178                 exit(1);
     182                return -1;
    179183            }
    180184            else
    … …  
    202206        *v = 2.0 * atand(w);
    203207        *r = q * (1+w2) / ( 1.0 + w2*lambda );
    204         return;  /* Near-parabolic orbit */
     208        return 0;  /* Near-parabolic orbit */
    205209    }
    206210 
    … …  
    232236        *r = q * (1.0+e) / ( 1.0 + e*cosd(*v) );
    233237    }
     238    return 0;
    234239 
    235240} /* vrc */
    236241
    237242/* For RCS Only -- Do Not Edit */
    238 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: twobody.c,v $ $Date: 2001-10-22 12:08:28 $ $Revision: 1.1 $ $Name: not supported by cvs2svn $"};
     243static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: twobody.c,v $ $Date: 2004-06-15 16:52:41 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/utc_gst.c

    r1719 r2551  
    1 #include "P_.h"
    21#include "astro.h"
    32
    4 static double gmst0 P_((double mjd));
     3static double gmst0 (double mj);
    54
    6 /* given a modified julian date, mjd, and a universally coordinated time, utc,
     5/* given a modified julian date, mj, and a universally coordinated time, utc,
    76 * return greenwich mean siderial time, *gst.
    8  * N.B. mjd must be at the beginning of the day.
     7 * N.B. mj must be at the beginning of the day.
    98 */
    109void
    11 utc_gst (mjd, utc, gst)
    12 double mjd;
    13 double utc;
    14 double *gst;
     10utc_gst (double mj, double utc, double *gst)
    1511{
    16         static double lastmjd = -18981;
     12        static double lastmj = -18981;
    1713        static double t0;
    1814
    19         if (mjd != lastmjd) {
    20             t0 = gmst0(mjd);
    21             lastmjd = mjd;
     15        if (mj != lastmj) {
     16            t0 = gmst0(mj);
     17            lastmj = mj;
    2218        }
    2319        *gst = (1.0/SIDRATE)*utc + t0;
    … …  
    2521}
    2622
    27 /* given a modified julian date, mjd, and a greenwich mean siderial time, gst,
     23/* given a modified julian date, mj, and a greenwich mean siderial time, gst,
    2824 * return universally coordinated time, *utc.
    29  * N.B. mjd must be at the beginning of the day.
     25 * N.B. mj must be at the beginning of the day.
    3026 */
    3127void
    32 gst_utc (mjd, gst, utc)
    33 double mjd;
    34 double gst;
    35 double *utc;
     28gst_utc (double mj, double gst, double *utc)
    3629{
    37         static double lastmjd = -10000;
     30        static double lastmj = -10000;
    3831        static double t0;
    3932
    40         if (mjd != lastmjd) {
    41             t0 = gmst0 (mjd);
    42             lastmjd = mjd;
     33        if (mj != lastmj) {
     34            t0 = gmst0 (mj);
     35            lastmj = mj;
    4336        }
    4437        *utc = gst - t0;
    … …  
    5043 */
    5144static double
    52 gmst0 (mjd)
    53 double mjd;     /* date at 0h UT in julian days since MJD0 */
     45gmst0 (
     46double mj)      /* date at 0h UT in julian days since MJD0 */
    5447{
    5548        double T, x;
    5649
    57         T = ((int)(mjd - 0.5) + 0.5 - J2000)/36525.0;
     50        T = ((int)(mj - 0.5) + 0.5 - J2000)/36525.0;
    5851        x = 24110.54841 +
    5952                (8640184.812866 + (0.093104 - 6.2e-6 * T) * T) * T;
    … …  
    6760/* original routine by elwood; has a secular drift of 0.08s/cty */
    6861static double
    69 tnaught (mjd)
    70 double mjd;     /* julian days since 1900 jan 0.5 */
     62tnaught (mj)
     63double mj;      /* julian days since 1900 jan 0.5 */
    7164{
    72         double dmjd;
     65        double dmj;
    7366        int m, y;
    7467        double d;
    7568        double t, t0;
    7669
    77         mjd_cal (mjd, &m, &d, &y);
    78         cal_mjd (1, 0., y, &dmjd);
    79         t = dmjd/36525;
    80         t0 = 6.57098e-2 * (mjd - dmjd) -
     70        mjd_cal (mj, &m, &d, &y);
     71        cal_mjd (1, 0., y, &dmj);
     72        t = dmj/36525;
     73        t0 = 6.57098e-2 * (mj - dmj) -
    8174             (24 - (6.6460656 + (5.1262e-2 + (t * 2.581e-5))*t) -
    8275                   (2400 * (t - (((double)y - 1900)/100))));
    … …  
    9083  char *argv[];
    9184{
    92         double mjd, gst;
    93         while (scanf("%lf", &mjd) == 1) {
    94                 mjd -= MJD0;
    95                 gst = tnaught(mjd);
    96                 printf("%17.9f %10.7f %10.7f\n", mjd + MJD0, gst, gmst0(mjd));
     85        double mj, gst;
     86        while (scanf("%lf", &mj) == 1) {
     87                mj -= MJD0;
     88                gst = tnaught(mj);
     89                printf("%17.9f %10.7f %10.7f\n", mj + MJD0, gst, gmst0(mj));
    9790        }
    9891}
    … …  
    10093
    10194/* For RCS Only -- Do Not Edit */
    102 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: utc_gst.c,v $ $Date: 2001-10-22 12:08:28 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     95static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: utc_gst.c,v $ $Date: 2004-06-15 16:52:41 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/vector.h

    r1719 r2551  
    22#define __SATVECTOR_H
    33
    4 /* $Id: vector.h,v 1.2 2001-10-22 12:08:28 cmv Exp $ */
     4/* $Id: vector.h,v 1.3 2004-06-15 16:52:41 cmv Exp $ */
    55
    66#define dotp(A,B) ((A).x*(B).x+(A).y*(B).y+(A).z*(B).z)
    … …  
    1616
    1717/* For RCS Only -- Do Not Edit
    18  * @(#) $RCSfile: vector.h,v $ $Date: 2001-10-22 12:08:28 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $
     18 * @(#) $RCSfile: vector.h,v $ $Date: 2004-06-15 16:52:41 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $
    1919 */

  • trunk/SophyaExt/XephemAstroLib/vsop87.c

    r1719 r2551  
    3535
    3636#include <math.h>
    37 #include "P_.h"
     37
    3838#include "astro.h"
    3939#include "vsop87.h"
    … …  
    5151 *    Input :
    5252 *
    53  *    mjd      modified julian date, counted from J1900.0
     53 *    mj       modified julian date, counted from J1900.0
    5454 *             time scale : dynamical time TDB.
    5555 *
    … …  
    9898 ******************************************************************/
    9999int
    100 vsop87 (mjd, obj, prec, ret)
    101 double mjd;
    102 int obj;
    103 double prec;
    104 double *ret;
     100vsop87 (double mj, int obj, double prec, double *ret)
    105101{
    106102    static double (*vx_map[])[3] = {            /* data tables */
    … …  
    134130    /* time and its powers */
    135131    t[0] = 1.0;
    136     t[1] = (mjd - J2000)/VSOP_A1000;
     132    t[1] = (mj - J2000)/VSOP_A1000;
    137133    for (i = 2; i <= VSOP_MAXALPHA; ++i) t[i] = t[i-1] * t[1];
    138134    t_abs[0] = 1.0;
    … …  
    211207
    212208/* For RCS Only -- Do Not Edit */
    213 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: vsop87.c,v $ $Date: 2001-10-22 12:08:28 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};
     209static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: vsop87.c,v $ $Date: 2004-06-15 16:52:41 $ $Revision: 1.4 $ $Name: not supported by cvs2svn $"};

  • trunk/SophyaExt/XephemAstroLib/vsop87.h

    r1719 r2551  
    8686extern int vn_neptune[][3];
    8787
    88 extern int vsop87 P_((double mjd, int obj, double prec, double *ret));
     88extern int vsop87 (double mj, int obj, double prec, double *ret);
    8989
    9090
    9191/* For RCS Only -- Do Not Edit
    92  * @(#) $RCSfile: vsop87.h,v $ $Date: 2001-10-22 12:08:28 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $
     92 * @(#) $RCSfile: vsop87.h,v $ $Date: 2004-06-15 16:52:41 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $
    9393 */

  • trunk/SophyaExt/XephemAstroLib/vsop87_data.c

    r1719 r2551  
    69866986
    69876987/* For RCS Only -- Do Not Edit */
    6988 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: vsop87_data.c,v $ $Date: 2001-10-22 12:08:28 $ $Revision: 1.2 $ $Name: not supported by cvs2svn $"};
     6988static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: vsop87_data.c,v $ $Date: 2004-06-15 16:52:41 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};
Note: See TracChangeset for help on using the changeset viewer.