source: Sophya/trunk/SophyaExt/XAstroPack/xastropack.h@ 1679

#ifndef XASTROPACK_H
#define XASTROPACK_H

#include "machdefs.h"
#include <string.h>
#include <string>

#ifdef __cplusplus
extern "C" {  /* extern "C" */
#endif
#include "XAstro/P_.h"
#include "XAstro/astro.h"
#ifdef __cplusplus
}             /* extern "C" */
#endif

enum TypAstroCoord {
  TypCoordUndef  =  (unsigned long)  (0),

  // Pour indiquer que les coordonnees sont en (heure=[0,24[,degre=[-90,90])
  TypCoordHD     =  (unsigned long)  (1 << 20),
  // Pour indiquer que les coordonnees sont en (degre=[0,360[,degre=[-90,90])
  TypCoordDD     =  (unsigned long)  (1 << 21),
  // Pour indiquer que les coordonnees sont en (radian=[0,2Pi[,radian=[-Pi/2,Pi/2])
  TypCoordRR     =  (unsigned long)  (1 << 22),

  // Coordonnees Equatoriales alpha,delta
  TypCoordEq     =  (unsigned long)  (1 << 0),
  TypCoordEqStd  =  (unsigned long) ((1 << 0) | (1 << 20)),
  // Coordonnees Galactiques gLong, gLat
  TypCoordGal    =  (unsigned long)  (1 << 1),
  TypCoordGalStd =  (unsigned long) ((1 << 1) | (1 << 21)),
  // Coordonnees Horizontales azimuth,altitude
  TypCoordHor    =  (unsigned long)  (1 << 2),
  TypCoordHorStd =  (unsigned long) ((1 << 2) | (1 << 21)),
  // Coordonnees Ecliptiques EclLong,EclLat
  TypCoordEcl    =  (unsigned long)  (1 << 3),
  TypCoordEclStd =  (unsigned long) ((1 << 3) | (1 << 21))
};

// ------------------- Utilitaires -------------------
int  CoordConvertToStd(TypAstroCoord typ,double& coord1,double& coord2);

/*!
\brief Pour remettre la valeur "val" dans la dynamique [0.,range[.
       Si "vmax" different de "range", c'est la borne superieure
       qui peut etre atteinte
       (si elle est depassee, on soustrait "range").
\verbatim
  r>0 vmax>0
  r=24. vmax=24.   -> mets dans [  0,+24[ borne sup exclue
  r=24. vmax=12.   -> mets dans ]-12,+12] borne inf exclue
\endverbatim
*/
inline void InRange(double *val,double range)
             {*val-=range*floor(*val/range);}
inline void InRange(double *val,double range,double vmax)
         {InRange(val,range); if(*val>vmax) *val-=range;}

// ------------------- Gestions des temps -------------------
/*! \ingroup XAstroPack
\brief Compute true Julian day from MJD
*/
inline double TrueJDfrMJD(double mjd) {return mjd + MJD0;}

/*! \ingroup XAstroPack
\brief Compute MJD from true Julian day
*/
inline double MJDfrTrueJD(double jd) {return jd - MJD0;}

/*! \ingroup XAstroPack
\brief Compute MJD from date
\verbatim
 MJD =  modified Julian date (number of days elapsed since 1900 jan 0.5),
\endverbatim
*/
inline double MJDfrDate(double dy,int mn,int yr)
       {double mjd; cal_mjd(mn,dy,yr,&mjd); return mjd;}

/*! \ingroup XAstroPack
\brief Compute date from MJD
*/
inline void DatefrMJD(double mjd,double *dy,int *mn,int *yr)
                     {mjd_cal(mjd,mn,dy,yr);}

/*! \ingroup XAstroPack
\brief  Given a mjd, return the year as a double.
*/
inline double YearfrMJD(double mjd)
       {double yr; mjd_year(mjd,&yr); return yr;}

/*! \ingroup XAstroPack
\brief Given a decimal year, return mjd
*/
inline double MJDfrYear(double yr)
       {double mjd; year_mjd(yr,&mjd); return mjd;}

/*! \ingroup XAstroPack
\brief Given a mjd, return the year and number of days since 00:00 Jan 1
\warning: if mjd = 2 January -> number of days = 1
*/
inline void YDfrMJD(double mjd,double *dy,int *yr)
       {mjd_dayno(mjd,yr,dy);}

/*! \ingroup XAstroPack
\brief Given a year, 
*/
inline int IsLeapYear(int y)
       {return isleapyear(y);}

/*! \ingroup XAstroPack
\brief given an mjd, set *dow to 0..6 according to which day of the week it falls on (0=sunday).
\return return 0 if ok else -1 if can't figure it out.
*/
inline int DayOrder(double mjd,int *dow)
       {return mjd_dow(mjd,dow);}

/*! \ingroup XAstroPack
\brief given a mjd, return the the number of days in the month.
*/
inline int DaysInMonth(double mjd)
       {int ndays; mjd_dpm(mjd,&ndays); return ndays;}

/*! \ingroup XAstroPack
\brief Given a mjd, truncate it to the beginning of the whole day
*/
inline double MJDat0hFrMJD(double mjd) {return mjd_day(mjd);}

/*! \ingroup XAstroPack
\brief Given a mjd, return the number of hours past midnight of the whole day
*/
inline double HfrMJD(double mjd) {return mjd_hr(mjd);}

/*! \ingroup XAstroPack
\brief Give GST from UTC
\verbatim
 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.
\endverbatim
*/
inline double GSTfrUTC(double mjd0,double utc)
       {double gst; utc_gst(mjd0,utc,&gst) ; return gst;}

/*! \ingroup XAstroPack
\brief Give UTC from GST
\verbatim
 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.
\endverbatim
*/                                                                             
inline double UTCfrGST(double mjd0,double gst)
       {double utc; gst_utc(mjd0,gst,&utc); return utc;}

/*! \ingroup XAstroPack
\brief return local sidereal time from greenwich mean siderial time and longitude
\param precis : if not zero, then correct for obliquity and nutation
\warning no nutation or obliquity correction are done.
*/                                                                             
inline double LSTfrGST(double gst,double geolng)
       {double lst = gst + deghr(geolng); InRange(&lst,24.); return lst;}

double GST0(double mjd0);
double LSTfrMJD(double mjd,double geolng);
void HMSfrHdec(double hd,int *h,int *mn,double *s);
double HdecfrHMS(int h,int mn,double s);
string ToStringHMS(int h,int mn,double s);
string ToStringHdec(double hd);

// ------------------- Calculs Divers -------------------
void Precess(double mjd1,double mjd2,double ra1,double dec1,double *ra2,double *dec2);

/*! \ingroup XAstroPack
\brief Given apparent altitude find airmass.
*/                                                                             
inline double AirmassfrAlt(double alt)
       {double x; alt = degrad(alt); airmass(alt,&x); return x;}

/*! \ingroup XAstroPack
\brief given geocentric time "jd" and coords of a distant object at "ra/dec" (J2000),
find the difference "hcp" in time between light arriving at earth vs the sun.
\return "hcp" must be subtracted from "geocentric jd" to get "heliocentric jd".
\warning "jd" is the TRUE Julian day (jd = mjd+MJD0).
*/
inline double HelioCorr(double jd,double ra,double dec)
       {double hcp; ra=hrrad(ra); dec=degrad(dec); heliocorr(jd,ra,dec,&hcp); return hcp;}

// ------------------- Transformation de coordonnees -------------------
/*! \ingroup XAstroPack
\brief Give the hour angle from local sideral time and right ascencion
\warning right ascencion should be first precessed to date of interest
\warning no nutation or obliquity correction are done.
*/                                                                             
// Attention au probleme de la discontinuite 0h <==> 24h
// ts=1  ra=23 ; (ts-ra)=-22 <-12 --> ha = +2 = +24 + (ts-ra)
// ts=23 ra=1  ; (ts-ra)=+22 >+12 --> ha = -2 = -24 + (ts-ra)
inline double HafrRaTS(double lst,double ra)
       {double ha = lst - ra; InRange(&ha,24.,12.); return ha;}

/*! \ingroup XAstroPack
\brief Give the local sideral time and the hour angle return the right ascencion
\warning right ascencion is the value precessed to date of interest
\warning no nutation or obliquity correction are done.
*/                                                                             
inline double RafrHaTS(double lst,double ha)
       {double ra = lst - ha; InRange(&ra,24.); return ra;}

void EqtoGal(double mjd,double ra,double dec,double *glng,double *glat);
void GaltoEq(double mjd,double glng,double glat,double *ra,double *dec);
void EqHtoHor(double geolat,double ha,double dec,double *az,double *alt);
void HortoEqH(double geolat,double az,double alt,double *ha,double *dec);
void EqtoHor(double geolat,double lst,double ra,double dec,double *az,double *alt);
void HortoEq(double geolat,double lst,double az,double alt,double *ra,double *dec);
void EqtoEcl(double mjd,double ra,double dec,double *eclng,double *eclat);
void EcltoEq(double mjd,double eclng,double eclat,double *ra,double *dec);

// ------------------- Positions des astres -------------------
void SunPos(double mjd,double *eclsn,double *ecbsn,double *rsn);
void MoonPos(double mjd,double *lmn,double *bmn,double *rho);
void PlanetPos(double mjd,int numplan,double *sunecl,double *sunecb,double *sundist
               ,double *geodist,double *geoecl,double *geoecb
               ,double *diamang,double *mag);
/*! \ingroup XAstroPack
\brief Same as PlanetPos above with less arguments
*/
inline void PlanetPos(double mjd,int numplan,double *geoecl,double *geoecb
                      ,double *geodist,double *diamang)
{
 double sunecl,sunecb,sundist,mag;
 PlanetPos(mjd,numplan,&sunecl,&sunecb,&sundist,geodist,geoecl,geoecb,diamang,&mag);
}

/*! \ingroup XAstroPack
\brief Give Jupiter position
*/
inline void JupiterPos(double mjd,double *ecl,double *ecb,double *geodist,double *diamang)
{
  PlanetPos(mjd,JUPITER,ecl,ecb,geodist,diamang);
}

/*! \ingroup XAstroPack
\brief Give Saturn position
*/
inline void SaturnPos(double mjd,double *ecl,double *ecb,double *geodist,double *diamang)
{
  PlanetPos(mjd,SATURN,ecl,ecb,geodist,diamang);
}

#endif