%module Astro

%title "Astro"

%include typemaps.i

%{
#include "P_.h"
#include "astro.h"
#include "circum.h"
#include "preferences.h"
/* #include "chap95.h" - internal */
/* #include "vsop87.h" - internal */
%}

%init %{
	Tcl_PkgProvide(interp, "Astro", "3.1");
%}

/* conversions among hours (of ra), degrees and radians. */
%{
#undef degrad
#undef raddeg
#undef hrdeg
#undef deghr
#undef hrrad
#undef radhr
#undef temp
%}

%inline double degrad(double x)	{ return ((x)*PI/180.); }
	/* convert degrees into radians */
%inline double raddeg(double x)	{ return ((x)*180./PI); }
	/* convert radians into degrees */
%inline double hrdeg(double x)	{ return ((x)*15.); }
	/* convert hours into degrees */
%inline double deghr(double x)	{ return ((x)/15.); }
	/* convert degrees into hours */
%inline double hrrad(double x)	{ return degrad(hrdeg(x)); }
	/* convert hours into radians */
%inline double radhr(double x)	{ return deghr(raddeg(x)); }
	/* convert radians into hours */

#define	SIDRATE		.9972695677
	/* ratio of from synodic (solar) to sidereal (stellar) rate */

#define	SPD	(24.0*3600.0)
	/* seconds per day */

#define	MERCURY	0
	/* Mercury index for calling plans */
#define	VENUS	1
	/* Venus index for calling plans */
#define	MARS	2
	/* Mars index for calling plans */
#define	JUPITER	3
	/* Jupiter index for calling plans */
#define	SATURN	4
	/* Saturn index for calling plans */
#define	URANUS	5
	/* Uranus index for calling plans */
#define	NEPTUNE	6
	/* Neptune index for calling plans */
#define	PLUTO	7
	/* Pluto index for calling plans */
#define	SUN	8
	/* Sun index for symmetry */
#define	MOON	9
	/* Moon index for symmetry */

#define MJD0  2415020.0
	/* starting point for MJD calculations  */
#define J2000 (2451545.0 - MJD0)
	/* epoch 2000 in MJD */

/* aa_hadec.c */
%apply double *OUTPUT { double *ha };
%apply double *OUTPUT { double *dec };
extern void aa_hadec (double lat, double alt, double az, double *ha, double *dec);
/* given geographical latitude (north+, radians), lat, altitude (up+, radians),
 alt, and azimuth (angle round to the east from north+, radians),
 return hour angle (radians), ha, and declination (radians), dec.
<br>
 [returns list of hour angle and declination]
 */
%apply double *OUTPUT { double *alt };
%apply double *OUTPUT { double *az };
extern void hadec_aa (double lat, double ha, double dec, double *alt, double *az);
/* given geographical (north+, radians), lat, hour angle (radians), ha, and
 declination (radians), dec, return altitude (up+, radians), alt, and
 azimuth (angle round to the east from north+, radians),
<br>
 [returns list of altitude and azimuth]
 */

/* aberration.c */
%apply double *INOUT { double *lam };
%apply double *INOUT { double *bet };
extern void ab_ecl (double mjd, double lsn, double *lam, double *bet);
/* apply aberration correction to ecliptical coordinates lam and bet
 (in radians) for a given time mjd and handily supplied longitude of sun,
 lsn (in radians).
<br>
[returns list of ecliptic longitude and latitude]
 */
%apply double *INOUT { double *ra };
%apply double *INOUT { double *dec };
extern void ab_eq (double mjd, double lsn, double *ra, double *dec);
/* apply aberration correction to equatoreal coordinates ra and dec
 (in radians) for a given time mjd and handily supplied longitude of sun,
 lsn (in radians).
<br>
[returns list of right ascension and declination]
 */

/* anomaly.c */
%apply double *OUTPUT { double *nu };
%apply double *OUTPUT { double *ea };
extern void anomaly (double ma, double s, double *nu, double *ea);
/* given the mean anomaly, ma, and the eccentricity, s, of elliptical motion,
 find the true anomaly, *nu, and the eccentric anomaly, *ea.
 all angles in radians.
<br>
[returns list of true anomaly and eccentric anomaly]
 */

/* comet.c */
%apply double *OUTPUT { double *lpd };
%apply double *OUTPUT { double *psi };
%apply double *OUTPUT { double *rp };
%apply double *OUTPUT { double *rho };
%apply double *OUTPUT { double *lam };
%apply double *OUTPUT { double *bet };

extern void comet ( double mjd, double ep, double inc, double ap, double qp,
    double om, double *lpd, double *psi, double *rp, double *rho, double *lam,
    double *bet);
/* given a modified Julian date, mjd, and a set of heliocentric parabolic
 orbital elements referred to the epoch of date (mjd):
<br>   ep:   epoch of perihelion,
<br>   inc:  inclination,
<br>   ap:   argument of perihelion (equals the longitude of perihelion minus the
	   longitude of ascending node)
<br>   qp:   perihelion distance,
<br>   om:   longitude of ascending node;
<br> find:
<br>   lpd:  heliocentric longitude, 
<br>   psi:  heliocentric latitude,
<br>   rp:   distance from the sun to the planet, 
<br>   rho:  distance from the Earth to the planet,
<br>   lam:  geocentric ecliptic longitude, 
<br>   bet:  geocentric ecliptic latitude,
<br>         none are corrected for light time, ie, they are the true values for
	   the given instant.
<br>
 all angles are in radians, all distances in AU.
 mutual perturbation corrections with other solar system objects are not
 applied. corrections for nutation and abberation must be made by the caller.
 The RA and DEC calculated from the fully-corrected ecliptic coordinates are
 then the apparent geocentric coordinates. Further corrections can be made,
 if required, for atmospheric refraction and geocentric parallax.
<br>
[returns a list of lpd, psi, rp, rho, lam, and bet]
 */

/* deltat.c */
extern double deltat (double mjd);
/* calculate  DeltaT = ET - UT in seconds.  Describes the irregularities
 of the Earth rotation rate in the ET time scale.
 */

/* eq_ecl.c */
%apply double *OUTPUT { double *lat };
%apply double *OUTPUT { double *lng };
extern void eq_ecl (double mjd, double ra, double dec, double *lat, double *lng);
/* given the modified Julian date, mjd, and an equatorial ra and dec, each in
 radians, find the corresponding geocentric ecliptic latitude, *lat, and
 longititude, *lng, also each in radians.
 correction for the effect on the angle of the obliquity due to nutation is
 not included.
<br>
[returns list of lat and lng]
 */
%apply double *OUTPUT { double *ra };
%apply double *OUTPUT { double *dec };
extern void ecl_eq (double mjd, double lat, double lng, double *ra, double *dec);
/* given the modified Julian date, mjd, and a geocentric ecliptic latitude,
 *lat, and longititude, *lng, each in radians, find the corresponding
 equitorial ra and dec, also each in radians.
 correction for the effect on the angle of the obliquity due to nutation is
 not included.
<br>
[returns list of ra and dec]
 */

/* eq_gal.c */
%apply double *OUTPUT { double *lat };
%apply double *OUTPUT { double *lng };
extern void eq_gal (double mjd, double ra, double dec, double *lat, double *lng);
/* given ra and dec, each in radians, for the given epoch, find the
 corresponding galactic latitude, *lat, and longititude, *lng, also each in
 radians.
<br>
[returns list of lat and lng]
 */
%apply double *OUTPUT { double *ra };
%apply double *OUTPUT { double *dec };
extern void gal_eq (double mjd, double lat, double lng, double *ra,
    double *dec);
/* given galactic latitude, lat, and longititude, lng, each in radians, find
 the corresponding equatorial ra and dec, also each in radians, at the 
 given epoch.
<br>
[returns list of ra and dec]
 */

/* mjd.c */
%apply double *OUTPUT { double *mjd };
extern void cal_mjd (int mn, double dy, int yr, double *mjd);
/* given a date in months, mn, days, dy, years, yr,
 return the modified Julian date (number of days elapsed since 1900 jan 0.5),
 *mjd.
<br>
[returns a list of mjd]
 */
%apply int *OUTPUT { int *mn };
%apply double *OUTPUT { double *dy };
%apply int *OUTPUT { int *yr };
extern void mjd_cal (double mjd, int *mn, double *dy, int *yr);
/* given the modified Julian date (number of days elapsed since 1900 jan 0.5,),
 mjd, return the calendar date in months, *mn, days, *dy, and years, *yr.
<br>
[returns a list of mn, dy, and yr]
 */
%apply int *OUTPUT { int *dow };
extern int mjd_dow (double mjd, int *dow);
/* given an mjd, set *dow to 0..6 according to which day of the week it falls
 on (0=sunday).
 return 0 if ok else -1 if can't figure it out.
[returns a list of status and dow]
 */
%apply int *OUTPUT { int *ndays };
extern void mjd_dpm (double mjd, int *ndays);
	/* given a mjd, return the the number of days in the month.  */
%apply double *OUTPUT { double *yr };
extern void mjd_year (double mjd, double *yr);
	/* given a mjd, return the year as a double. */
%apply double *OUTPUT { double *mjd };
extern void year_mjd (double y, double *mjd);
	/* given a decimal year, return mjd */
%apply double *INOUT { double *t };
extern void rnd_second (double *t);
	/* round a time in days, *t, to the nearest second, IN PLACE. */
extern double mjd_day (double jd);
	/* given an mjd, truncate it to the beginning of the whole day */
extern double mjd_hr (double jd);
	/* given an mjd, return the number of hours past midnight of the whole day */
%apply double *INOUT { double *v };
extern void range (double *v, double r);
	/* insure 0 <= *v < r. */

/* moon.c */

%apply double *OUTPUT { double *lam };
%apply double *OUTPUT { double *bet };
%apply double *OUTPUT { double *rho };
%apply double *OUTPUT { double *msp };
%apply double *OUTPUT { double *mdp };
extern void moon (double mjd, double *lam, double *bet, double *rho, double *msp, double *mdp);
/*
 given the mjd, find the geocentric ecliptic longitude, lam, and latitude,
   bet, and geocentric distance, rho in a.u. for the moon.  also return
   the sun's mean anomaly, *msp, and the moon's mean anomaly, *mdp.
<br>
[returns a list of lam, bet, rho, msp, and mdp]
*/

/* mooncolong.c */
%apply double *OUTPUT { double *cp };
%apply double *OUTPUT { double *kp };
%apply double *OUTPUT { double *ap };
%apply double *OUTPUT { double *sp };
extern void moon_colong (double mjd, double lt, double lg, double *cp, double *kp, double *ap, double *sp);
/* given a Julian date and a lunar location, find selenographic colongitude of
 rising sun, lunar latitude of subsolar point, illuminated fraction, and alt
 of sun at the given location. Any pointer may be 0 if not interested.
 From Bruning and Talcott, October 1995 _Astronomy_, page 76.
 N.B. lunar coordinates use +E, but selenograhic colongs are +W.
<br>
[returns a list of colongitude, latitude, fraction and altitude]
 */

/* nutation.c */
%apply double *OUTPUT { double *deps };
%apply double *OUTPUT { double *dpsi };
extern void nutation (double mjd, double *deps, double *dpsi);
/* given the modified JD, mjd, find the nutation in obliquity, *deps, and
 the nutation in longitude, *dpsi, each in radians.
<br>
[returns a list of deps and dpsi]
 */
%apply double *INOUT { double *ra };
%apply double *INOUT { double *dec };
extern void nut_eq(double mjd, double *ra, double *dec);
/* given the modified JD, mjd, correct, IN PLACE, the right ascension *ra
 and declination *dec (both in radians) for nutation.
<br>
[returns a list of ra and dec]
 */

/* obliq.c */
%apply double *OUTPUT { double *eps };
extern void obliquity(double mjd, double *eps);
/* given the modified Julian date, mjd, find the mean obliquity of the
 ecliptic, *eps, in radians.
<br>
[returns a list of eps]
 */

/* parallax.c */
%apply double *IN { double *rho };
%apply double *OUTPUT { double *aha };
%apply double *OUTPUT { double *adec };
extern void ta_par(double tha, double tdec, double phi, double ht, double *rho, double *aha, double *adec);
/* given true ha and dec, tha and tdec, the geographical latitude, phi, the
 height above sea-level (as a fraction of the earths radius, 6378.16km),
 ht, and the geocentric distance rho in Earth radii(!), find the apparent
 ha and dec, aha and adec allowing for parallax.
 all angles in radians. ehp is the angle subtended at the body by the
 earth's equator.
<br>
[returns a list of aha and adec]
 */

/* plans.c */
%apply double *OUTPUT { double *lpd0 };
%apply double *OUTPUT { double *psi0 };
%apply double *OUTPUT { double *rp0 };
%apply double *OUTPUT { double *rho0 };
%apply double *OUTPUT { double *lam };
%apply double *OUTPUT { double *bet };
%apply double *OUTPUT { double *dia };
%apply double *OUTPUT { double *mag };
extern void plans(double mjd, int p, double *lpd0, double *psi0,
    double *rp0, double *rho0, double *lam, double *bet, double *dia,
    double *mag);
/* given a modified Julian date, mjd, and a planet, p, find:
<br>   lpd0: heliocentric longitude, 
<br>   psi0: heliocentric latitude,
<br>   rp0:  distance from the sun to the planet, 
<br>   rho0: distance from the Earth to the planet,
<br>         none corrected for light time, ie, they are the true values for the
         given instant.
<br>   lam:  geocentric ecliptic longitude, 
<br>   bet:  geocentric ecliptic latitude,
<br>         each corrected for light time, ie, they are the apparent values as
	   seen from the center of the Earth for the given instant.
<br>   dia:  angular diameter in arcsec at 1 AU, 
<br>   mag:  visual magnitude when 1 AU from sun and earth at 0 phase angle.
<br>
 all angles are in radians, all distances in AU.
<br>
 corrections for nutation and abberation must be made by the caller. The RA 
   and DEC calculated from the fully-corrected ecliptic coordinates are then
   the apparent geocentric coordinates. Further corrections can be made, if
   required, for atmospheric refraction and geocentric parallax.
<br>
[returns a list of lpd0, psi0, rp0, rho0, lam, bet, dia, and mag]
*/

/* precess.c */
%apply double *INOUT { double *ra };
%apply double *INOUT { double *dec };
extern void precess(double mjd1, double mjd2, double *ra, double *dec);
/* corrects ra and dec, both in radians, for precession from epoch 1 to epoch 2.
 the epochs are given by their modified JDs, mjd1 and mjd2, respectively.
 N.B. ra and dec are modifed IN PLACE.
<br>
[returns a list of ra and dec]
 */

/* reduce.c */
%apply double *OUTPUT { double *inc }
%apply double *OUTPUT { double *ap }
%apply double *OUTPUT { double *om }
extern void reduce_elements(double mjd0, double mjd, double inc0,
    double ap0, double om0, double *inc, double *ap, double *om);
/* convert those orbital elements that change from epoch mjd0 to epoch mjd.
 */

/* refract.c */
%apply double *OUTPUT { double *aa }
extern void refract(double pr, double tr, double ta, double *aa);
/* correct the true altitude, ta, for refraction to the apparent altitude, aa,
 each in radians, given the local atmospheric pressure, pr, in mbars, and
 the temperature, tr, in degrees C.
 */
%apply double *OUTPUT { double *ta }
extern void unrefract(double pr, double tr, double aa, double *ta);
/* correct the apparent altitude, aa, for refraction to the true altitude, ta,
 each in radians, given the local atmospheric pressure, pr, in mbars, and
 the temperature, tr, in degrees C.
 */

/* riset.c */
%apply double *OUTPUT { double *lstr }
%apply double *OUTPUT { double *lsts }
%apply double *OUTPUT { double *azr }
%apply double *OUTPUT { double *azs }
%apply int *OUTPUT { int *status }
extern void riset(double ra, double dec, double lat, double dis,
    double *lstr, double *lsts, double *azr, double *azs, int *status);
/* given the true geocentric ra and dec of an object, the observer's latitude,
   lat, and a horizon displacement correction, dis, all in radians, find the
   local sidereal times and azimuths of rising and setting, lstr/s
   and azr/s, also all in radians, respectively.
 dis is the vertical displacement from the true position of the horizon. it
   is positive if the apparent position is higher than the true position.
   said another way, it is positive if the shift causes the object to spend
   longer above the horizon. for example, atmospheric refraction is typically
   assumed to produce a vertical shift of 34 arc minutes at the horizon; dis
   would then take on the value +9.89e-3 (radians). On the other hand, if
   your horizon has hills such that your apparent horizon is, say, 1 degree
   above sea level, you would allow for this by setting dis to -1.75e-2
   (radians).

 This version contributed by Konrad Bernloehr, Nov. 1996

 status: 0=normal; 1=never rises; -1=circumpolar.
 In case of non-zero status, all other returned variables are undefined.
 */

/* sphcart.c */
%apply double *OUTPUT { double *x }
%apply double *OUTPUT { double *y }
%apply double *OUTPUT { double *z }
extern void sphcart(double l, double b, double r, double *x, double *y, double *z);
/* transformation from spherical to cartesian coordinates */
%apply double *OUTPUT { double *l }
%apply double *OUTPUT { double *b }
%apply double *OUTPUT { double *r }
extern void cartsph(double x, double y, double z, double *l, double *b, double *r);
/* transformation from cartesian to spherical coordinates */

/* sun.c */
%apply double *OUTPUT { double *lsn }
%apply double *OUTPUT { double *rsn }
%apply double *OUTPUT { double *bsn }
extern void sunpos(double mjd, double *lsn, double *rsn, double *bsn);
/* given the modified JD, mjd, return the true geocentric ecliptic longitude
   of the sun for the mean equinox of the date, *lsn, in radians, the
   sun-earth distance, *rsn, in AU, and the latitude *bsn, in radians
   (since this is always <= 1.2 arcseconds, in can be neglected by
   calling with bsn = NULL).

 if the APPARENT ecliptic longitude is required, correct the longitude for
   nutation to the true equinox of date and for aberration (light travel time,
   approximately  -9.27e7/186000/(3600*24*365)*2*pi = -9.93e-5 radians).
 */

/* utc_gst.c */
%apply double *OUTPUT { double *gst }
extern void utc_gst(double mjd, double utc, double *gst);
/* given a modified julian date, mjd, and a universally coordinated time, utc,
 return greenwich mean siderial time, *gst.
 N.B. mjd must be at the beginning of the day.
 */
%apply double *OUTPUT { double *utc }
extern void gst_utc(double mjd, double gst, double *utc);
/* given a modified julian date, mjd, and a greenwich mean siderial time, gst,
 return universally coordinated time, *utc.
 N.B. mjd must be at the beginning of the day.
 */