Index: /trunk/SophyaExt/XAstroPack/xastropack.cc
===================================================================
--- /trunk/SophyaExt/XAstroPack/xastropack.cc	(revision 1677)
+++ /trunk/SophyaExt/XAstroPack/xastropack.cc	(revision 1678)
@@ -18,5 +18,6 @@
 // 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
+//              angle horaire en heures [-12,12] (-12=12) (ha) 
+//              temps sideral du lieu: tsid=ha+ra (ou lst)
 // GALACTIQUE: longitude en degres [0,360[  (glng)
 //             latitude en degres [-90,90]  (glat)
@@ -51,4 +52,7 @@
 /*! \ingroup XAstroPack
 \brief Compute MJD from date
+\verbatim
+ MJD =  modified Julian date (number of days elapsed since 1900 jan 0.5),
+\endverbatim
 */
 double MJDfrDate(double dy,int mn,int yr)
@@ -97,4 +101,47 @@
 
 /*! \ingroup XAstroPack
+\brief Given a year, 
+*/
+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.
+*/
+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.
+*/
+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
+*/
+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
+*/
+double HfrMJD(double mjd)
+{
+  return mjd_hr(mjd);
+}
+
+/*! \ingroup XAstroPack
 \brief Give GST from UTC
 \verbatim
@@ -131,5 +178,5 @@
 */                                                                             
 double GST0(double mjd0)
-/* Copie depuis le code de Xephem car pas prototype */
+/* Copie depuis le code de Xephem (utc_gst.c) car pas prototype*/
 {
  double T, x;
@@ -143,4 +190,31 @@
 
 /*! \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.
+*/                                                                             
+double LSTfrGST(double gst,double geolng)
+{
+  double lst = gst + geolng *12./180.;
+  InRange(&lst,24.);
+  return lst;
+}
+
+/*! \ingroup XAstroPack
+\brief return local sidereal time from modified julian day and longitude
+\warning nutation or obliquity correction are taken into account.
+*/                                                                             
+double LSTfrMJD(double mjd,double geolng)
+{
+  double eps,lst,deps,dpsi;
+  utc_gst(mjd_day(mjd),mjd_hr(mjd),&lst);
+  lst += geolng *12./180.;
+  obliquity(mjd,&eps);
+  nutation(mjd,&deps,&dpsi);
+  lst += dpsi*cos(eps+deps) *12./M_PI;
+  return lst;
+}
+
+/*! \ingroup XAstroPack
 \brief Compute precession between 2 dates.
 */                                                                             
@@ -150,5 +224,5 @@
  dec1 *= PI/180.;  // radians
  precess(mjd1,mjd2,&ra1,&dec1);
- *ra2 = ra1*12./PI; if(*ra2>24.) *ra2 -= 24.; if(*ra2==24.) *ra2 = 0.;
+ *ra2 = ra1*12./PI; InRange(ra2,24.);
  *dec2 = dec1*180./PI;
 }
@@ -166,14 +240,43 @@
 
 /*! \ingroup XAstroPack
-\brief Give the hour angle from sideral time and right ascencion
-*/                                                                             
-double HafrRaTS(double gst,double ra)
-{
-  double ha = gst - ra;
+\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.
+*/                                                                             
+double HafrRaTS(double lst,double ra)
+{
+  double ha = lst - 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.;
+  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.
+*/                                                                             
+double RafrHaTS(double lst,double ha)
+{
+  double ra = lst - ha;
+  InRange(&ra,24.);
+  return ra;
+}
+
+/*! \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).
+*/
+double HelioCorr(double jd,double ra,double dec)
+{
+ double hcp;
+ ra  *= PI/12.;   // radians
+ dec *= PI/180.;  // radians
+ heliocorr(jd,ra,dec,&hcp);
+ return hcp;
 }
 
@@ -247,5 +350,5 @@
 
 /*! \ingroup XAstroPack
-\brief Convert equatorial coordinates into galactic coordinates
+\brief Convert equatorial coordinates for the given epoch into galactic coordinates
 */                                                                             
 void EqtoGal(double mjd,double ra,double dec, double *glng,double *glat)
@@ -256,10 +359,10 @@
  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.;
+ *glng *= 180./PI; InRange(glng,360.);
  *glat *= 180./PI;
 }
 
 /*! \ingroup XAstroPack
-\brief Convert galactic coordinates into equatorial coordinates
+\brief Convert galactic coordinates into equatorial coordinates at the given epoch
 */                                                                             
 void GaltoEq(double mjd,double glng,double glat,double *ra,double *dec)
@@ -269,12 +372,12 @@
  glat *= PI/180.;  // radians
  gal_eq (mjd,glat,glng,ra,dec);
- *ra  *= 12./PI; if(*ra>24.) *ra -= 24.; if(*ra==24.) *ra = 0.;
+ *ra  *= 12./PI; InRange(ra,24.);
  *dec *= 180./PI;
 }
 
 /*! \ingroup XAstroPack
-\brief Convert equatorial coordinates into horizontal coordinates
-*/                                                                             
-void EqtoHor(double geolat,double ha,double dec,double *az,double *alt)
+\brief Convert equatorial coordinates (with hour angle instead of right ascension) into horizontal coordinates.
+*/                                                                             
+void EqHtoHor(double geolat,double ha,double dec,double *az,double *alt)
 // Coordonnees equatoriales -> Coordonnees horizontales
 {
@@ -284,11 +387,11 @@
  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)
+ *az  *= 180./PI; InRange(az,360.);
+}
+
+/*! \ingroup XAstroPack
+ Convert horizontal coordinates into equatorial coordinates (with hour angle instead of right ascension).
+*/                                                                             
+void HortoEqH(double geolat,double az,double alt,double *ha,double *dec)
 // Coordonnees horizontales -> Coordonnees equatoriales
 {
@@ -297,11 +400,43 @@
  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.;
+ *ha  *= 12./PI; InRange(ha,24.,12.);
  *dec *= 180./PI;
 }
 
 /*! \ingroup XAstroPack
-\brief Convert equatorial coordinates into ecliptic coordinates
+\brief Convert equatorial coordinates into horizontal coordinates.
+*/                                                                             
+void EqtoHor(double geolat,double lst,double ra,double dec,double *az,double *alt)
+// Coordonnees equatoriales -> Coordonnees horizontales
+{
+ double ha = lst - ra;
+ if(ha==-12.) ha=12.; InRange(&ha,24.,12.);
+ geolat *= PI/180.;
+ ha  *= PI/12.;   // radians
+ dec *= PI/180.;  // radians
+ hadec_aa (geolat,ha,dec,alt,az);
+ *alt *= 180./PI;
+ *az  *= 180./PI; InRange(az,360.);
+}
+
+/*! \ingroup XAstroPack
+ Convert horizontal coordinates into equatorial coordinates.
+*/                                                                             
+void HortoEq(double geolat,double lst,double az,double alt,double *ra,double *dec)
+// Coordonnees horizontales -> Coordonnees equatoriales
+{
+ double ha;
+ geolat *= PI/180.;
+ alt *= PI/180.;  // radians
+ az  *= PI/180.;  // radians
+ aa_hadec (geolat,alt,az,&ha,dec);
+ ha *= 12./PI;
+ *ra = lst - ha; InRange(ra,24.);
+ *dec *= 180./PI;
+}
+
+/*! \ingroup XAstroPack
+\brief Convert equatorial coordinates into geocentric ecliptic coordinates given the modified Julian date.
+\warning Correction for the effect on the angle of the obliquity due to nutation is not included.
 */                                                                             
 // Attention, j'ai modifie eq_ecl.c pour proteger NaN
@@ -315,10 +450,11 @@
  dec *= PI/180.;  // radians
  eq_ecl(mjd,ra,dec,eclat,eclng);
- *eclng *= 180./PI; if(*eclng>360.) *eclng -= 360.; if(*eclng==360.) *eclng = 0.;
+ *eclng *= 180./PI; InRange(eclng,360.);
  *eclat *= 180./PI;
 }
 
 /*! \ingroup XAstroPack
-\brief Convert ecliptic coordinates into equatorial coordinates
+\brief Convert geocentric ecliptic coordinates into equatorial coordinates given the modified Julian date.
+\warning Correction for the effect on the angle of the obliquity due to nutation is not included.
 */                                                                             
 void EcltoEq(double mjd,double eclng,double eclat,double *ra,double *dec)
@@ -328,5 +464,5 @@
  eclng *= PI/180.;  // radians
  ecl_eq(mjd,eclat,eclng,ra,dec); 
- *ra  *= 12./PI; if(*ra>24.) *ra -= 24.; if(*ra==24.) *ra = 0.;
+ *ra  *= 12./PI; InRange(ra,24.);
  *dec *= 180./PI;
 }
@@ -344,7 +480,7 @@
 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.;
+ double rsn;
+ sunpos(mjd,eclsn,&rsn,ecbsn);
+ *eclsn *= 180./PI; InRange(eclsn,360.);
  *ecbsn *= 180./PI;
 }
@@ -363,5 +499,5 @@
   double rho,msp,mdp;
   moon(mjd,eclmn,ecbmn,&rho,&msp,&mdp);
- *eclmn *= 180./PI; if(*eclmn>360.) *eclmn -= 360.; if(*eclmn==360.) *eclmn = 0.;
+ *eclmn *= 180./PI; InRange(eclmn,360.);
  *ecbmn *= 180./PI;
 }
@@ -388,7 +524,7 @@
 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.;
+ double lpd0,psi0,rp0,rho0,mag;
+ plans(mjd,numplan,&lpd0,&psi0,&rp0,&rho0,ecl,ecb,diamang,&mag);
+ *ecl *= 180./PI; InRange(ecl,360.);
  *ecb *= 180./PI;
 }
Index: /trunk/SophyaExt/XAstroPack/xastropack.h
===================================================================
--- /trunk/SophyaExt/XAstroPack/xastropack.h	(revision 1677)
+++ /trunk/SophyaExt/XAstroPack/xastropack.h	(revision 1678)
@@ -46,10 +46,19 @@
 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);
 double GST0(double mjd0);
+double LSTfrGST(double gst,double geolng);
+double LSTfrMJD(double mjd,double geolng);
 void Precess(double mjd1,double mjd2,double ra1,double dec1,double *ra2,double *dec2);
 double AirmassfrAlt(double alt);
-double HafrRaTS(double gst,double ra);
+double HafrRaTS(double lst,double ra);
+double RafrHaTS(double lst,double ha);
+double HelioCorr(double jd,double ra,double dec);
 void HMSfrHdec(double hd,int *h,int *mn,double *s);
 double HdecfrHMS(int h,int mn,double s);
@@ -58,6 +67,8 @@
 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 EqtoHor(double geolat,double ha,double dec,double *az,double *alt);
-void HortoEq(double geolat,double az,double alt,double *ha,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);
@@ -69,3 +80,19 @@
 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;}
+
 #endif
