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

#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),

  TypCoordH0     =  (unsigned long)  (1 << 10), // heure=[0,24[
  TypCoordH1     =  (unsigned long)  (1 << 11), // heure=]-12,12]
  TypCoordD0     =  (unsigned long)  (1 << 12), // degre=[0,360[
  TypCoordD1     =  (unsigned long)  (1 << 13), // degre=]-180,180]
  TypCoordD2     =  (unsigned long)  (1 << 14), // degre=[-90,90]
  TypCoordR0     =  (unsigned long)  (1 << 15), // degre=[0,2Pi[
  TypCoordR1     =  (unsigned long)  (1 << 16), // degre=]-Pi,Pi]
  TypCoordR2     =  (unsigned long)  (1 << 17), // degre=[-Pi/2,Pi/2]

  // 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),
  // Pour indiquer que les coordonnees sont en (heure=]-12,12],degre=[-90,90])
  TypCoordH1D    =  (unsigned long)  (1 << 23),
  // Pour indiquer que les coordonnees sont en (degre=]-180,180],degre=[-90,90])
  TypCoordD1D    =  (unsigned long)  (1 << 24),
  // Pour indiquer que les coordonnees sont en (radian=]-Pi,Pi],radian=[-Pi/2,Pi/2])
  TypCoordR1R    =  (unsigned long)  (1 << 25),

  // 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
  (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,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);
double YearfrMJD(double mjd);
double MJDfrYear(double yr);
void YDfrMJD(double mjd,double *dy,int *yr);
int IsLeapYear(int y);
int DayOrder(double mjd,int *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);
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 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);
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