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

#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"
#undef mjd
#undef lat
#undef lng
#undef tz
#undef temp
#undef pressure
#undef elev
#undef dip
#undef epoch
#undef tznm
#undef mjed
#ifdef __cplusplus
}             /* extern "C" */
#endif

enum TypAstroCoord {
  //------ Coordonnees de type (A,B) avec:
  //      A,B <==> alpha,delta
  //      A,B <==> longitude,latitude
  //      A,B <==> azimuth,altitude ou elevation

  //------ Les unites des coordonnees
  TypUniteH       =  (unsigned long)  (1 << 0), // coord en heure
  TypUniteD       =  (unsigned long)  (1 << 1), // coord en degre
  TypUniteR       =  (unsigned long)  (1 << 2), // coord en radian

  //------ Les unites des coordonnees A et B
  TypCoord1H      =  (unsigned long)  (1 << 3), // coord A en heure
  TypCoord1D      =  (unsigned long)  (1 << 4), // coord A en degre
  TypCoord1R      =  (unsigned long)  (1 << 5), // coord A en radian
  TypCoord2H      =  (unsigned long)  (1 << 6), // coord B en heure
  TypCoord2D      =  (unsigned long)  (1 << 7), // coord B en degre
  TypCoord2R      =  (unsigned long)  (1 << 8), // coord B en radian

  //------ Le type d'intervalle des coordonnees A et B
  // Coord A intervalle... ex: [0,24[ ou [0,360[ ou [0,2Pi[
  TypCoord1C      =  (unsigned long)  (1 << 10),
  // Coord A intervalle centre... ex: ]-12,12] ou ]-180,180] ou ]-Pi,Pi]
  TypCoord1L      =  (unsigned long)  (1 << 11),
  // Coord B intervalle... ex:  [0,12] ou [0,180] ou [0,Pi] (colatitude)
  TypCoord2C      =  (unsigned long)  (1 << 12),
  // Coord B intervalle centre... ex: [-6,6] ou [-90,90] ou [-Pi/2,Pi/2] (latitude)
  TypCoord2L      =  (unsigned long)  (1 << 13),

  //------ Les systemes de coordonnees astronomiques
  // Coordonnees Equatoriales alpha,delta
  TypCoordEq     =  (unsigned long)  (1 << 15),
  // Coordonnees Galactiques gLong, gLat
  TypCoordGal    =  (unsigned long)  (1 << 16),
  // Coordonnees Horizontales azimuth,altitude
  TypCoordHor    =  (unsigned long)  (1 << 17),
  // Coordonnees Ecliptiques EclLong,EclLat
  TypCoordEcl    =  (unsigned long)  (1 << 18),

  //------ Les systemes de coordonnees astronomiques "standard"
  // Coordonnees Equatoriales alpha,delta
  TypCoordEqStd     =  (unsigned long)  (TypCoordEq |TypCoord1C|TypCoord1H|TypCoord2L|TypCoord2D),
  // Coordonnees Galactiques gLong, gLat
  TypCoordGalStd    =  (unsigned long)  (TypCoordGal|TypCoord1C|TypCoord1D|TypCoord2L|TypCoord2D),
  // Coordonnees Horizontales azimuth,altitude
  TypCoordHorStd    =  (unsigned long)  (TypCoordHor|TypCoord1C|TypCoord1D|TypCoord2L|TypCoord2D),
  // Coordonnees Ecliptiques EclLong,EclLat
  TypCoordEclStd    =  (unsigned long)  (TypCoordEcl|TypCoord1C|TypCoord1D|TypCoord2L|TypCoord2D),

  //------ undef pour la forme
  TypCoordUndef  =  (unsigned long)  (0)
};

// ------------------- Utilitaires -------------------
unsigned long CoordConvertToStd(unsigned long typ,double* coord1,double* coord2);
inline unsigned long CoordConvertToStd(unsigned long typ,double &coord1,double &coord2)
       {return CoordConvertToStd(typ,&coord1,&coord2);}

unsigned long GetCoordUnit(int coordnum,unsigned long typ);

unsigned long DecodeTypAstro(const char *ctype);
inline unsigned long DecodeTypAstro(const string stype)
       {return DecodeTypAstro(stype.c_str());}
string        DecodeTypAstro(unsigned long typ);

void ToCoLat(double* val,unsigned long typ);
inline void ToCoLat(double &val,unsigned long typ)
       {ToCoLat(&val,typ);}

void InRangeLat(double* val,unsigned long typ);
inline void InRangeLat(double &val,unsigned long typ)
       {InRangeLat(&val,typ);}

void InRangeCoLat(double* val,unsigned long typ);
inline void InRangeCoLat(double &val,unsigned long typ)
       {InRangeCoLat(&val,typ);}

/*!
\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 "val" dans [  0,+24[ borne sup exclue
  r=24. vmax=12.   -> mets "val" dans ]-12,+12] borne inf exclue
  (ATTENTION: ca ne marche pas pour la latitude [-90,90]
              car -90 ne peut etre inclus dans l'intervalle)
\endverbatim
*/
inline void InRange(double *val,double range)
             {*val-=range*floor(*val/range);}
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;}
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;}

double MJDfrDate(double dy,int mn,int yr);

void DatefrMJD(double mjd,double *dy,int *mn,int *yr);
inline void DatefrMJD(double mjd,double &dy,int &mn,int &yr)
       {DatefrMJD(mjd,&dy,&mn,&yr);}

double YearfrMJD(double mjd);

double MJDfrYear(double yr);

void YDfrMJD(double mjd,double *dy,int *yr);
inline void YDfrMJD(double mjd,double &dy,int &yr)
       {YDfrMJD(mjd,&dy,&yr);}

int IsLeapYear(int y);

int DayOrder(double mjd,int *dow);
inline int DayOrder(double mjd,int &dow)
       {return DayOrder(mjd,&dow);}

int DaysInMonth(double mjd);
double MJDat0hFrMJD(double mjd);
double HfrMJD(double mjd);
double GSTfrUTC(double mjd0,double utc);
double UTCfrGST(double mjd0,double gst);

/*! \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);
inline void HMSfrHdec(double hd,int &h,int &mn,double &s)
       {HMSfrHdec(hd,&h,&mn,&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);
inline void Precess(double mjd1,double mjd2,double ra1,double dec1,double &ra2,double &dec2)
       {Precess(mjd1,mjd2,ra1,dec1,&ra2,&dec2);}

inline void PrecessInPlace(double mjd1,double mjd2,double* ra,double* dec) {
  double ra1=*ra, dec1=*dec;
  Precess(mjd1,mjd2,ra1,dec1,ra,dec);
}
inline void PrecessInPlace(double mjd1,double mjd2,double &ra,double &dec) {
  double ra1=ra, dec1=dec;
  Precess(mjd1,mjd2,ra1,dec1,ra,dec);
}

double AirmassfrAlt(double alt);
double HelioCorr(double jd,double ra,double dec);

// ------------------- Transformation de coordonnees -------------------
/*! \ingroup XAstroPack
\brief Give the hour angle from local sideral time and right ascension
\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);
inline void EqtoGal(double mjd,double ra,double dec,double &glng,double &glat)
       {EqtoGal(mjd,ra,dec,&glng,&glat);}


void GaltoEq(double mjd,double glng,double glat,double *ra,double *dec);
inline void GaltoEq(double mjd,double glng,double glat,double &ra,double &dec)
       {GaltoEq( mjd, glng, glat, &ra, &dec);}

void EqHtoHor(double geolat,double ha,double dec,double *az,double *alt);
inline void EqHtoHor(double geolat,double ha,double dec,double &az,double &alt)
       {EqHtoHor( geolat, ha, dec, &az, &alt);}

void HortoEqH(double geolat,double az,double alt,double *ha,double *dec);
inline void HortoEqH(double geolat,double az,double alt,double &ha,double &dec)
       {HortoEqH( geolat, az, alt, &ha, &dec);}

void EqtoHor(double geolat,double lst,double ra,double dec,double *az,double *alt);
inline void EqtoHor(double geolat,double lst,double ra,double dec,double &az,double &alt)
       {EqtoHor( geolat, lst, ra, dec, &az, &alt);}

void HortoEq(double geolat,double lst,double az,double alt,double *ra,double *dec);
inline void HortoEq(double geolat,double lst,double az,double alt,double &ra,double &dec)
       {HortoEq( geolat, lst, az, alt, &ra, &dec);}

void EqtoEcl(double mjd,double ra,double dec,double *eclng,double *eclat);
inline void EqtoEcl(double mjd,double ra,double dec,double &eclng,double &eclat)
       {EqtoEcl( mjd, ra, dec, &eclng, &eclat);}

void EcltoEq(double mjd,double eclng,double eclat,double *ra,double *dec);
inline void EcltoEq(double mjd,double eclng,double eclat,double &ra,double &dec)
       {EcltoEq( mjd, eclng, eclat, &ra, &dec);}

// ------------------- Positions des astres -------------------
void SunPos(double mjd,double *eclsn,double *ecbsn,double *rsn);
inline void SunPos(double mjd,double &eclsn,double &ecbsn,double &rsn)
       {SunPos( mjd, &eclsn, &ecbsn, &rsn);}

void MoonPos(double mjd,double *lmn,double *bmn,double *rho);
inline void MoonPos(double mjd,double &lmn,double &bmn,double &rho)
       {MoonPos( mjd, &lmn, &bmn, &rho);}

void PlanetPos(double mjd,PLCode numplan,double *sunecl,double *sunecb,double *sundist
               ,double *geodist,double *geoecl,double *geoecb
               ,double *diamang,double *mag);
inline void PlanetPos(double mjd,PLCode numplan,double &sunecl,double &sunecb,double &sundist
               ,double &geodist,double &geoecl,double &geoecb
               ,double &diamang,double &mag)
       {PlanetPos(mjd,numplan,&sunecl,&sunecb,&sundist,&geodist,&geoecl,&geoecb,&diamang,&mag);}

/*! \ingroup XAstroPack
\brief Same as PlanetPos above with less arguments
*/
inline void PlanetPos(double mjd,PLCode 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);
}
inline void PlanetPos(double mjd,PLCode 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);
}
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);
}
inline void SaturnPos(double mjd,double &ecl,double &ecb,double &geodist,double &diamang)
{
  PlanetPos(mjd,SATURN,ecl,ecb,geodist,diamang);
}

#endif