source: Sophya/trunk/SophyaExt/XAstroPack/xastropack.cc@ 1628

#include <math.h>
#include <stdio.h>
#include "xastropack.h"

/*!
  \defgroup XAstroPack XAstroPack module
  This module contains simple programs to perform various
  astronomical computation (based on the libastro of Xephem).

  \verbatim
// TEMPS: modified Julian date (mjd) (number of days elapsed since 1900 jan 0.5)
//        jour [1,31] (dy)
//        mois [1,12] (mn)
//        annee       (yr)
//        universal time [0,24[ (utc)
//        Greenwich mean siderial [0,24[ (gst)
//        Greenwich mean siderial at 0h UT [0,24[ (gst0)
// EQUATORIALE: ascension droite en heures [0,24[ (ra)
//              declinaison en degres [-90,90]    (dec)
//              angle horaire en heures [-12,12] (-12=12) (ha)  tsid=ha+ra
// GALACTIQUE: longitude en degres [0,360[  (glng)
//             latitude en degres [-90,90]  (glat)
// HORIZONTAL: azimuth en degres [0,360[   (az)
//                 (angle round to the east from north+)
//             altitude en degres [-90,90] (alt)
// ECLIPTIQUE: lontitude ecliptique en degres [0,360[ (eclng)
//                 (angle round counter clockwise from the vernal equinoxe)
//             latitude ecliptique en degres [-90,90] (eclat)
// GEOGRAPHIE: longitude en degres ]-180,180]        (geolng)
//                 (angle + vers l'ouest, - vers l'est)
//             latitude en degres [-90,90] (north>0) (geolat)
  \endverbatim
*/

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

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

/*! \ingroup XAstroPack
\brief Compute MJD from date
*/
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
*/
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.
*/
double YearfrMJD(double mjd)
{
  double yr;
  mjd_year(mjd,&yr);
  return yr;
}

/*! \ingroup XAstroPack
\brief Given a decimal year, return mjd
*/
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
*/
void YDfrMJD(double mjd,double *dy,int *yr)
{
  mjd_dayno(mjd,yr,dy);
}

/*! \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
*/
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
*/                                                                             
double UTCfrGST(double mjd0,double gst)
{
  double utc;
  gst_utc(mjd0,gst,&utc);
  return utc;
}

/*! \ingroup XAstroPack
\brief gmst0() - return Greenwich Mean Sidereal Time at 0h UT
\param mjd = date at 0h UT in julian days since MJD0
*/                                                                             
double GST0(double mjd0)
/* Copie depuis le code de Xephem car pas prototype */
{
 double T, x;
 T = ((int)(mjd0 - 0.5) + 0.5 - J2000)/36525.0;
 x = 24110.54841 +
     (8640184.812866 + (0.093104 - 6.2e-6 * T) * T) * T;
 x /= 3600.0;
 range(&x, 24.0);
 return (x);
}

/*! \ingroup XAstroPack
\brief Compute precession between 2 dates.
*/                                                                             
void Precess(double mjd1,double mjd2,double ra1,double dec1,double *ra2,double *dec2)
{
 ra1  *= PI/12.;   // radians
 dec1 *= PI/180.;  // radians
 precess(mjd1,mjd2,&ra1,&dec1);
 *ra2 = ra1*12./PI; if(*ra2>24.) *ra2 -= 24.; if(*ra2==24.) *ra2 = 0.;
 *dec2 = dec1*180./PI;
}

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

/*! \ingroup XAstroPack
\brief Give the hour angle from sideral time and right ascencion
*/                                                                             
double HafrRaTS(double gst,double ra)
{
  double ha = gst - ra;
  // 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)
  if(ha==-12.) ha = 12.; if(ha<-12.) ha += 24.; if(ha>12.) ha -= 24.;
  return ha;
}

/*! \ingroup XAstroPack
\brief Give a time in h:mn:s from a decimal hour
\verbatim
// INPUT: hd
// OUTPUT: h mn s   (h,mn,s >=< 0)
// REMARQUE: si hd<0 alors h<0 ET mn<0 ET s<0
// EX: 12.51 -> h=12  mn=30  s=10 ; 
//    -12.51 -> h=-12 mn=-30 s=-10 ; 
\endverbatim
*/                                                                             
void HMSfrHdec(double hd,int *h,int *mn,double *s)
{
 int sgn=1;
 if(hd<0.) {sgn=-1; hd*=-1.;}
 *h  = int(hd);
 *mn = int((hd-(double)(*h))*60.);
 *s  = (hd - (double)(*h) - (double)(*mn)/60.)*3600.;
 // pb precision
 if(*s<0.) *s = 0.;
 if(*s>60. || *s==60.) {*s-=60.; *mn+=1;} // s=double attention comparaison
 if(*mn<0) *mn = 0;
 if(*mn>=60) {*mn-=60; *h+=1;}
 *h *= sgn; *mn *= sgn; *s *= (double)sgn;
}

/*! \ingroup XAstroPack
\brief Give a decimal hour from a time in h:mn:s
\verbatim
// INPUT: h , mn , s  (h,mn,s >=< 0)
// RETURN:  en heures decimales
// REMARQUE: pour avoir hd=-12.51 <- h=-12 mn=-30 s=-10
\endverbatim
*/                                                                             
double HdecfrHMS(int h,int mn,double s)
{
 return ((double)h + (double)mn/60. + s/3600.);
}

/*! \ingroup XAstroPack
\brief Give a time string from a time in h:mn:s
\verbatim
// INPUT: h , mn , s   (h,mn,s >=< 0)
// RETURN: string h:mn:s
\endverbatim
*/                                                                             
string ToStringHMS(int h,int mn,double s)
{
 double hd = HdecfrHMS(h,mn,s); // put in range
 HMSfrHdec(hd,&h,&mn,&s);
 char str[128];
 if(hd<0.)
   sprintf(str,"-%d:%d:%.3f",-h,-mn,-s);
 else
   sprintf(str,"%d:%d:%.3f",h,mn,s);
 string dum = str;
 return dum;
}

/*! \ingroup XAstroPack
\brief Give a time string from a decimal hour
*/                                                                             
string ToStringHdec(double hd)
{
 int h,mn; double s;
 HMSfrHdec(hd,&h,&mn,&s);
 return ToStringHMS(h,mn,s);
}

/*! \ingroup XAstroPack
\brief Convert equatorial coordinates into galactic coordinates
*/                                                                             
void EqtoGal(double mjd,double ra,double dec, double *glng,double *glat)
// Coordonnees equatoriales -> Coordonnees galactiques
{
 ra  *= PI/12.;   // radians
 dec *= PI/180.;  // radians
 eq_gal(mjd,ra,dec,glat,glng);
 // Vraiment bizarre, sur Linux-g++ glng>=360 ne comprend pas glng==360 ! (CMV)
 *glng *= 180./PI; if(*glng>360.) *glng -= 360.; if(*glng==360.) *glng = 0.;
 *glat *= 180./PI;
}

/*! \ingroup XAstroPack
\brief Convert galactic coordinates into equatorial coordinates
*/                                                                             
void GaltoEq(double mjd,double glng,double glat,double *ra,double *dec)
// Coordonnees galactiques -> Coordonnees equatoriales
{
 glng  *= PI/180.;  // radians
 glat *= PI/180.;  // radians
 gal_eq (mjd,glat,glng,ra,dec);
 *ra  *= 12./PI; if(*ra>24.) *ra -= 24.; if(*ra==24.) *ra = 0.;
 *dec *= 180./PI;
}

/*! \ingroup XAstroPack
\brief Convert equatorial coordinates into horizontal coordinates
*/                                                                             
void EqtoHor(double geolat,double ha,double dec,double *az,double *alt)
// Coordonnees equatoriales -> Coordonnees horizontales
{
 geolat *= PI/180.;
 ha  *= PI/12.;   // radians
 dec *= PI/180.;  // radians
 hadec_aa (geolat,ha,dec,alt,az);
 *alt *= 180./PI;
 *az  *= 180./PI; if(*az>360.) *az -= 360.; if(*az==360.) *az = 0.;
}

/*! \ingroup XAstroPack
 Convert horizontal coordinates into equatorial coordinates
*/                                                                             
void HortoEq(double geolat,double az,double alt,double *ha,double *dec)
// Coordonnees horizontales -> Coordonnees equatoriales
{
 geolat *= PI/180.;
 alt *= PI/180.;  // radians
 az  *= PI/180.;  // radians
 aa_hadec (geolat,alt,az,ha,dec);
 *ha  *= 12./PI; 
  if(*ha==-12.) *ha = 12.; if(*ha<-12.) *ha += 24.; if(*ha>12.) *ha -= 24.;
 *dec *= 180./PI;
}

/*! \ingroup XAstroPack
\brief Convert equatorial coordinates into ecliptic coordinates
*/                                                                             
// Attention, j'ai modifie eq_ecl.c pour proteger NaN
// dans ecleq_aux :
// *q = (sy*ceps)-(cy*seps*sx*sw);
// if(*q<-1.) *q = -PI/2.; else if(*q>1.) *q = PI/2.; else *q = asin(*q);
void EqtoEcl(double mjd,double ra,double dec,double *eclng,double *eclat)
// Coordonnees equatoriales -> Coordonnees ecliptiques
{
 ra  *= PI/12.;   // radians
 dec *= PI/180.;  // radians
 eq_ecl(mjd,ra,dec,eclat,eclng);
 *eclng *= 180./PI; if(*eclng>360.) *eclng -= 360.; if(*eclng==360.) *eclng = 0.;
 *eclat *= 180./PI;
}

/*! \ingroup XAstroPack
\brief Convert ecliptic coordinates into equatorial coordinates
*/                                                                             
void EcltoEq(double mjd,double eclng,double eclat,double *ra,double *dec)
// Coordonnees ecliptiques -> Coordonnees equatoriales
{
 eclat *= PI/180.;  // radians
 eclng *= PI/180.;  // radians
 ecl_eq(mjd,eclat,eclng,ra,dec); 
 *ra  *= 12./PI; if(*ra>24.) *ra -= 24.; if(*ra==24.) *ra = 0.;
 *dec *= 180./PI;
}

/*! \ingroup XAstroPack
\brief Give Sun position
\verbatim
 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).
\endverbatim
*/
void SunPos(double mjd,double *eclsn,double *ecbsn)
{
  double rsn;
  sunpos(mjd,eclsn,&rsn,ecbsn);
 *eclsn *= 180./PI; if(*eclsn>360.) *eclsn -= 360.; if(*eclsn==360.) *eclsn = 0.;
 *ecbsn *= 180./PI;
}

/*! \ingroup XAstroPack
\brief Give Moon position
\verbatim
 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. 
 (for the mean equinox)
\endverbatim
*/
void MoonPos(double mjd,double *eclmn,double *ecbmn)
{
  double rho,msp,mdp;
  moon(mjd,eclmn,ecbmn,&rho,&msp,&mdp);
 *eclmn *= 180./PI; if(*eclmn>360.) *eclmn -= 360.; if(*eclmn==360.) *eclmn = 0.;
 *ecbmn *= 180./PI;
}

/*! \ingroup XAstroPack
\brief Give planet position
\verbatim
 * given a modified Julian date, mjd, and a planet, p, find:
 *   lpd0: heliocentric longitude, 
 *   psi0: heliocentric latitude,
 *   rp0:  distance from the sun to the planet, 
 *   rho0: distance from the Earth to the planet,
 *         none corrected for light time, ie, they are the true values for the
 *         given instant.
 *   lam:  geocentric ecliptic longitude, 
 *   bet:  geocentric ecliptic latitude,
 *         each corrected for light time, ie, they are the apparent values as
 *         seen from the center of the Earth for the given instant.
 *   dia:  angular diameter in arcsec at 1 AU, 
 *   mag:  visual magnitude when 1 AU from sun and earth at 0 phase angle.
 *   (for the mean equinox)
 \endverbatim
*/
void PlanetPos(double mjd,int numplan,double *ecl,double *ecb,double *diamang)
{
  double lpd0,psi0,rp0,rho0,mag;
  plans(mjd,numplan,&lpd0,&psi0,&rp0,&rho0,ecl,ecb,diamang,&mag);
 *ecl *= 180./PI; if(*ecl>360.) *ecl -= 360.; if(*ecl==360.) *ecl = 0.;
 *ecb *= 180./PI;
}

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

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

/*! \ingroup XAstroPack
\brief Given a coordinate type "typ", convert to standard for astropack
\verbatim
// Return : 0 = OK
//          1 = Unknown type of coordinates
//          2 = bad range for coord1
//          4 = bad range for coord2
//          6 = bad range for coord1 et coord2
\endverbatim
*/
int  CoordConvertToStd(TypAstroCoord typ,double& coord1,double& coord2)
{
  int rc = 0;

  // ---- Equatoriales alpha,delta
  //    - standard = [0,24[ , [-90,90]
  if(typ&TypCoordEq) {
    if(typ&TypCoordDD) {
      coord1 = coord1 / 180. * 12.;
    } else if(typ&TypCoordRR) {
      coord1 = coord1 / PI * 12.;
      coord2 = coord2 / PI * 180.;
    }
    if(coord1==24.) coord1 = 0.;
    if(coord1<0.   || coord1>=24.) rc+= 2;
    if(coord2<-90. || coord2>90. ) rc+= 4;

  // ---- Galactiques gLong, gLat
  // ---- Horizontales azimuth,altitude
  // ---- Ecliptiques EclLong,EclLat
  //    - standard = [0,360[ , [-90,90]
  } else if( typ&TypCoordGal || typ&TypCoordHor || typ&TypCoordEcl) {
    if(typ&TypCoordHD) {
      coord1 = coord1 / 12. * 180.;
    } else if(typ&TypCoordRR) {
      coord1 = coord1 / PI * 180.;
      coord2 = coord2 / PI * 180.;
    }
    if(coord1==360.) coord1 = 0.;
    if(coord1<0.   || coord1>=360.) rc+= 2;
    if(coord2<-90. || coord2>90. )  rc+= 4;

  } else {          // Coordonnees non-connues
    rc= 1;
  }

  return rc;
}