Index: trunk/SophyaExt/XAstroPack/xastropack.cc
===================================================================
--- trunk/SophyaExt/XAstroPack/xastropack.cc	(revision 1678)
+++ trunk/SophyaExt/XAstroPack/xastropack.cc	(revision 1679)
@@ -35,145 +35,6 @@
 
 /*! \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
-\verbatim
- MJD =  modified Julian date (number of days elapsed since 1900 jan 0.5),
-\endverbatim
-*/
-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 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
- 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
+\param mjd0 = date at 0h UT in julian days since MJD0
 */                                                                             
 double GST0(double mjd0)
@@ -190,93 +51,17 @@
 
 /*! \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.;
+\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 += deghr(geolng);
+  obliquity(mjd,&eps);
+  nutation(mjd,&deps,&dpsi);
+  lst += radhr(dpsi*cos(eps+deps));
   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.
-*/                                                                             
-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; InRange(ra2,24.);
- *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 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)
-  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;
 }
 
@@ -350,4 +135,16 @@
 
 /*! \ingroup XAstroPack
+\brief Compute precession between 2 dates.
+*/                                                                             
+void Precess(double mjd1,double mjd2,double ra1,double dec1,double *ra2,double *dec2)
+{
+ ra1  = hrrad(ra1);   // radians
+ dec1 = degrad(dec1);  // radians
+ precess(mjd1,mjd2,&ra1,&dec1);
+ *ra2 = radhr(ra1); InRange(ra2,24.);
+ *dec2 = raddeg(dec1);
+}
+
+/*! \ingroup XAstroPack
 \brief Convert equatorial coordinates for the given epoch into galactic coordinates
 */                                                                             
@@ -355,10 +152,10 @@
 // Coordonnees equatoriales -> Coordonnees galactiques
 {
- ra  *= PI/12.;   // radians
- dec *= PI/180.;  // radians
+ ra  = hrrad(ra);   // radians
+ dec = degrad(dec);  // radians
  eq_gal(mjd,ra,dec,glat,glng);
  // Vraiment bizarre, sur Linux-g++ glng>=360 ne comprend pas glng==360 ! (CMV)
- *glng *= 180./PI; InRange(glng,360.);
- *glat *= 180./PI;
+ *glng = raddeg(*glng); InRange(glng,360.);
+ *glat = raddeg(*glat);
 }
 
@@ -369,9 +166,9 @@
 // Coordonnees galactiques -> Coordonnees equatoriales
 {
- glng  *= PI/180.;  // radians
- glat *= PI/180.;  // radians
+ glng = degrad(glng);  // radians
+ glat = degrad(glat);  // radians
  gal_eq (mjd,glat,glng,ra,dec);
- *ra  *= 12./PI; InRange(ra,24.);
- *dec *= 180./PI;
+ *ra = radhr(*ra); InRange(ra,24.);
+ *dec = raddeg(*dec);
 }
 
@@ -382,10 +179,10 @@
 // Coordonnees equatoriales -> Coordonnees horizontales
 {
- geolat *= PI/180.;
- ha  *= PI/12.;   // radians
- dec *= PI/180.;  // radians
+ geolat = degrad(geolat);  // radians
+ ha  = hrrad(ha);   // radians
+ dec = degrad(dec);  // radians
  hadec_aa (geolat,ha,dec,alt,az);
- *alt *= 180./PI;
- *az  *= 180./PI; InRange(az,360.);
+ *alt = raddeg(*alt);
+ *az  = raddeg(*az); InRange(az,360.);
 }
 
@@ -396,10 +193,10 @@
 // Coordonnees horizontales -> Coordonnees equatoriales
 {
- geolat *= PI/180.;
- alt *= PI/180.;  // radians
- az  *= PI/180.;  // radians
+ geolat = degrad(geolat);  // radians
+ alt = degrad(alt);  // radians
+ az  = degrad(az);  // radians
  aa_hadec (geolat,alt,az,ha,dec);
- *ha  *= 12./PI; InRange(ha,24.,12.);
- *dec *= 180./PI;
+ *ha = radhr(*ha); InRange(ha,24.,12.);
+ *dec = raddeg(*dec);
 }
 
@@ -410,12 +207,11 @@
 // 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
+ double ha = lst - ra; InRange(&ha,24.,12.);
+ geolat = degrad(geolat);  // radians
+ ha = hrrad(ha);   // radians
+ dec = degrad(dec);  // radians
  hadec_aa (geolat,ha,dec,alt,az);
- *alt *= 180./PI;
- *az  *= 180./PI; InRange(az,360.);
+ *alt = raddeg(*alt);
+ *az  = raddeg(*az); InRange(az,360.);
 }
 
@@ -427,11 +223,10 @@
 {
  double ha;
- geolat *= PI/180.;
- alt *= PI/180.;  // radians
- az  *= PI/180.;  // radians
+ geolat = degrad(geolat);  // radians
+ alt = degrad(alt);  // radians
+ az  = degrad(az);  // radians
  aa_hadec (geolat,alt,az,&ha,dec);
- ha *= 12./PI;
- *ra = lst - ha; InRange(ra,24.);
- *dec *= 180./PI;
+ *ra = lst - radhr(ha); InRange(ra,24.);
+ *dec = raddeg(*dec);
 }
 
@@ -447,9 +242,9 @@
 // Coordonnees equatoriales -> Coordonnees ecliptiques
 {
- ra  *= PI/12.;   // radians
- dec *= PI/180.;  // radians
+ ra = hrrad(ra);   // radians
+ dec = degrad(dec);  // radians
  eq_ecl(mjd,ra,dec,eclat,eclng);
- *eclng *= 180./PI; InRange(eclng,360.);
- *eclat *= 180./PI;
+ *eclng = raddeg(*eclng); InRange(eclng,360.);
+ *eclat = raddeg(*eclat);
 }
 
@@ -461,9 +256,9 @@
 // Coordonnees ecliptiques -> Coordonnees equatoriales
 {
- eclat *= PI/180.;  // radians
- eclng *= PI/180.;  // radians
+ eclat = degrad(eclat);  // radians
+ eclng = degrad(eclng);  // radians
  ecl_eq(mjd,eclat,eclng,ra,dec); 
- *ra  *= 12./PI; InRange(ra,24.);
- *dec *= 180./PI;
+ *ra = radhr(*ra); InRange(ra,24.);
+ *dec = raddeg(*dec);
 }
 
@@ -472,16 +267,19 @@
 \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
+ of the sun for the mean equinox of the date, *eclsn, in degres, the
+ sun-earth distance, *rsn, in AU, and the latitude *ecbsn, in degres
  (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; InRange(eclsn,360.);
- *ecbsn *= 180./PI;
+ calling with ecbsn = NULL).
+ - REMARQUE:
+ * if the APPARENT ecliptic longitude is required, correct the longitude for
+ *   nutation to the true equinox of date and for aberration (light travel time,
+ *   approximately  -9.27e7/186000/(3600*24*365)*2*pi = -9.93e-5 radians).
+\endverbatim
+*/
+void SunPos(double mjd,double *eclsn,double *ecbsn,double *rsn)
+{
+ sunpos(mjd,eclsn,rsn,ecbsn);
+ *eclsn = raddeg(*eclsn); InRange(eclsn,360.);
+ if(ecbsn!=NULL) *ecbsn = raddeg(*ecbsn);
 }
 
@@ -495,10 +293,10 @@
 \endverbatim
 */
-void MoonPos(double mjd,double *eclmn,double *ecbmn)
-{
-  double rho,msp,mdp;
-  moon(mjd,eclmn,ecbmn,&rho,&msp,&mdp);
- *eclmn *= 180./PI; InRange(eclmn,360.);
- *ecbmn *= 180./PI;
+void MoonPos(double mjd,double *eclmn,double *ecbmn,double *rho)
+{
+  double msp,mdp;
+  moon(mjd,eclmn,ecbmn,rho,&msp,&mdp);
+ *eclmn = raddeg(*eclmn); InRange(eclmn,360.);
+ *ecbmn = raddeg(*ecbmn);
 }
 
@@ -507,41 +305,34 @@
 \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,
+ *   sunecl: heliocentric longitude, 
+ *   sunecb: heliocentric latitude,
+ *   sundist:  distance from the sun to the planet, 
+ *   geodist: 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,
+ *   geoecl:  geocentric ecliptic longitude, 
+ *   geoecb:  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, 
+ *   diamang:  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)
+ * all angles are in degres, all distances in AU.
+ *
+ * corrections for nutation and abberation must be made by the caller. The RA 
+ *   and DEC calculated from the fully-corrected ecliptic coordinates are then
+ *   the apparent geocentric coordinates. Further corrections can be made, if
+ *   required, for atmospheric refraction and geocentric parallax.
  \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; InRange(ecl,360.);
- *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);
+void PlanetPos(double mjd,int numplan,double *sunecl,double *sunecb,double *sundist
+              ,double *geodist,double *geoecl,double *geoecb
+              ,double *diamang,double *mag)
+{
+ plans(mjd,numplan,sunecl,sunecb,sundist,geodist,geoecl,geoecb,diamang,mag);
+ *geoecl = raddeg(*geoecl); InRange(geoecl,360.);
+ *geoecb = raddeg(*geoecb);
+ *sunecl = raddeg(*sunecl); InRange(sunecl,360.);
+ *sunecb = raddeg(*sunecb);
 }
 
@@ -564,8 +355,8 @@
   if(typ&TypCoordEq) {
     if(typ&TypCoordDD) {
-      coord1 = coord1 / 180. * 12.;
+      coord1 = deghr(coord1);
     } else if(typ&TypCoordRR) {
-      coord1 = coord1 / PI * 12.;
-      coord2 = coord2 / PI * 180.;
+      coord1 = radhr(coord1);
+      coord2 = raddeg(coord2);
     }
     if(coord1==24.) coord1 = 0.;
@@ -579,8 +370,8 @@
   } else if( typ&TypCoordGal || typ&TypCoordHor || typ&TypCoordEcl) {
     if(typ&TypCoordHD) {
-      coord1 = coord1 / 12. * 180.;
+      coord1 = hrdeg(coord1);
     } else if(typ&TypCoordRR) {
-      coord1 = coord1 / PI * 180.;
-      coord2 = coord2 / PI * 180.;
+      coord1 = raddeg(coord1);
+      coord2 = raddeg(coord2);
     }
     if(coord1==360.) coord1 = 0.;
