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

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

// 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)

double TrueJDfrMJD(double mjd)
{
  return mjd + MJD0;
}

double MJDfrTrueJD(double jd)
{
  return jd - MJD0;
}

double MJDfrDate(double dy,int mn,int yr)
{
  double mjd;
  cal_mjd(mn,dy,yr,&mjd);
  return mjd;
}

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

/* given a mjd, return the year as a double. */
double YearfrMJD(double mjd)
{
  double yr;
  mjd_year(mjd,&yr);
  return yr;
}

/* given a decimal year, return mjd */ 
double MJDfrYear(double yr)
{
  double mjd;
  year_mjd(yr,&mjd);
  return mjd;
}

/* given a mjd, return the year and number of days since 00:00 Jan 1 */
/* Attention: si mjd = 2 Janvier -> number of days = 1               */
void YDfrMJD(double mjd,double *dy,int *yr)
{
  mjd_dayno(mjd,yr,dy);
}

/* 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.
 */
double GSTfrUTC(double mjd0,double utc)
{
  double gst;
  utc_gst(mjd0,utc,&gst) ;
  return gst;
}

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

/* gmst0() - return Greenwich Mean Sidereal Time at 0h UT */
/* 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);
}

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;
}

/* given apparent altitude find airmass. */
double AirmassfrAlt(double alt)
{
  double x;
  alt *= PI/180.;  // radians
  airmass(alt,&x);
  return x;
}

/* donne l'angle horaire a partir du temps sideral et de l'ascension droite */
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;
}

void HMSfrHdec(double hd,int *h,int *mn,double *s)
// 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 ; 
{
 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;
}

double HdecfrHMS(int h,int mn,double s)
// INPUT: h , mn , s  (h,mn,s >=< 0)
// RETURN:  en heures decimales
// REMARQUE: pour avoir hd=-12.51 <- h=-12 mn=-30 s=-10
{
 return ((double)h + (double)mn/60. + s/3600.);
}

string ToStringHMS(int h,int mn,double s)
// INPUT: h , mn , s   (h,mn,s >=< 0)
// RETURN: string h:mn: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;
}

string ToStringHdec(double hd)
{
 int h,mn; double s;
 HMSfrHdec(hd,&h,&mn,&s);
 return ToStringHMS(h,mn,s);
}

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;
}

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;
}

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.;
}

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;
}

// 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;
}

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;
}

/* 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). */
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;
}

/* 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) */
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;
}

void PlanetPos(double mjd,int numplan,double *ecl,double *ecb,double *diamang)
/* 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) */
{
  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;
}

void JupiterPos(double mjd,double *ecl,double *ecb,double *diamang)
{
  PlanetPos(mjd,JUPITER,ecl,ecb,diamang);
}

void SaturnPos(double mjd,double *ecl,double *ecb,double *diamang)
{
  PlanetPos(mjd,SATURN,ecl,ecb,diamang);
}

/* Given a coordinate type "typ", convert to standard for astropack */
int  CoordConvertToStd(TypAstroCoord typ,double& coord1,double& coord2)
// Return : 0 = OK
//          1 = Type de coordonnees non connu
//          2 = Mauvais range pour coord1
//          4 = Mauvais range pour coord2
//          6 = Mauvais range pour coord1 et 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;
}