Index: /trunk/SophyaExt/XAstroPack/Makefile
===================================================================
--- /trunk/SophyaExt/XAstroPack/Makefile	(revision 1456)
+++ /trunk/SophyaExt/XAstroPack/Makefile	(revision 1456)
@@ -0,0 +1,15 @@
+include ../Mgr/Makefile.h
+ 
+all: $(LIB)libXAstroPack.a
+clean:
+	rm -f $(OBJ)xastropack.o
+	rm -f $(LIB)libXAstroPack.a
+$(LIB)libXAstroPack.a : $(OBJ)xastropack.o
+	$(ARCXX) $(ARCXXFLAGS) $@ $?
+ifeq ($(CXX),cxx)
+	$(CXX) $? $(CPPFLAGS) $(CXXFLAGS) -o $(OBJ)xx.x -Hf
+endif 
+ 
+$(OBJ)xastropack.o: xastropack.cc xastropack.h \
+ $(EXTINCPATH)XAstro/P_.h \
+ $(EXTINCPATH)XAstro/astro.h
Index: /trunk/SophyaExt/XAstroPack/extlib_list
===================================================================
--- /trunk/SophyaExt/XAstroPack/extlib_list	(revision 1456)
+++ /trunk/SophyaExt/XAstroPack/extlib_list	(revision 1456)
@@ -0,0 +1,1 @@
+xastro
Index: /trunk/SophyaExt/XAstroPack/extslb_list
===================================================================
--- /trunk/SophyaExt/XAstroPack/extslb_list	(revision 1456)
+++ /trunk/SophyaExt/XAstroPack/extslb_list	(revision 1456)
@@ -0,0 +1,1 @@
+xastro
Index: /trunk/SophyaExt/XAstroPack/objlist.list
===================================================================
--- /trunk/SophyaExt/XAstroPack/objlist.list	(revision 1456)
+++ /trunk/SophyaExt/XAstroPack/objlist.list	(revision 1456)
@@ -0,0 +1,1 @@
+xastropack.o
Index: /trunk/SophyaExt/XAstroPack/xastropack.cc
===================================================================
--- /trunk/SophyaExt/XAstroPack/xastropack.cc	(revision 1456)
+++ /trunk/SophyaExt/XAstroPack/xastropack.cc	(revision 1456)
@@ -0,0 +1,313 @@
+#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,23[ (utc)
+//        Greenwich mean siderial [0,23[ (gst)
+//        Greenwich mean siderial at 0h UT [0,23[ (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: altitude en degres [-90,90] (alt)
+//             azimuth en degres [0,360[   (az)
+//               (angle round to the east from north+)
+// ECLIPTIQUE: latitude ecliptique en degres [-90,90] (eclat)
+//             lontitude ecliptique en degres [0,360[ (eclng)
+//               (angle round counter clockwise from the vernal equinoxe)
+// GEOGRAPHIE: latitude en degres [-90,90] (north>0) (geolat)
+//             longitude en degres ]-180,180]        (geolng)
+//               (angle + vers l'ouest, - vers l'est)
+
+double TrueMJDfrMJD(double mjd)
+{
+  return mjd + 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 HdectoHMS(double hd,int *h,int *mn,double *s)
+// INPUT: hd
+// OUTPUT: h:mn:s
+// REMARQUE: si hd<0 alors h<0 mais toujours mn,s>=0
+{
+ 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;
+}
+
+double HMStoHdec(int h,int mn,double s)
+// INPUT: h , mn , s
+// RETURN: h:|mn|:|s| en heures decimales
+// REMARQUE: si h<0  return -h:mn:s
+// ERROR: mn<0 ou s<0 n'est pas correct, le programme les remet>0
+{
+ if(mn<0) {
+   cout<<"HMStoHdec: mn out of range <0 : "<<mn<<" changed to abs()"<<endl;
+   mn *= -1;
+ }
+ if(s<0.) {
+   cout<<"HMStoHdec: s out of range <0 : "<<s<<" changed to abs()"<<endl;
+   s *= -1.;
+ }
+ int sgn=1;  if(h<0) {sgn=-1; h*=-1;}
+ return ((double)h + (double)mn/60. + s/3600.)*(double)sgn;
+}
+
+string ToStringHMS(int h,int mn,double s)
+// INPUT: h , mn>=0 , s >=0
+// RETURN: string h:mn:s
+// REMARQUE: si h<0  return -h:mn:s
+// ERROR: mn<0 ou s<0 est une erreur !
+//        on prend la valeur absolue mn->|mn| , s->|s|
+{
+ double hd = HMStoHdec(h,mn,s); // put in range
+ HdectoHMS(hd,&h,&mn,&s);
+ char str[128];
+ sprintf(str,"%d:%d:%.3f",h,mn,s);
+ string dum = str;
+ return dum;
+}
+
+string ToStringHdec(double hd)
+{
+ int h,mn; double s;
+ HdectoHMS(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);
+}
Index: /trunk/SophyaExt/XAstroPack/xastropack.h
===================================================================
--- /trunk/SophyaExt/XAstroPack/xastropack.h	(revision 1456)
+++ /trunk/SophyaExt/XAstroPack/xastropack.h	(revision 1456)
@@ -0,0 +1,44 @@
+#ifndef XASTROPACK_H
+#define XASTROPACK_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
+
+double TrueMJDfrMJD(double mjd);
+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);
+double GSTfrUTC(double mjd0,double utc);
+double UTCfrGST(double mjd0,double gst);
+double GST0(double mjd0);
+void Precess(double mjd1,double mjd2,double ra1,double dec1,double *ra2,double *dec2);
+double AirmassfrAlt(double alt);
+double HafrRaTS(double gst,double ra);
+void HdectoHMS(double hd,int *h,int *mn,double *s);
+double HMStoHdec(int h,int mn,double s);
+string ToStringHMS(int h,int mn,double s);
+string ToStringHdec(double hd);
+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 EqtoEcl(double mjd,double ra,double dec,double *eclng,double *eclat);
+void EcltoEq(double mjd,double eclng,double eclat,double *ra,double *dec);
+void SunPos(double mjd,double *lsn,double *bsn);
+void MoonPos(double mjd,double *lmn,double *bmn);
+void PlanetPos(double mjd,int numplan,double *ecl,double *ecb,double *diamang);
+void JupiterPos(double mjd,double *ecl,double *ecb,double *diamang);
+void SaturnPos(double mjd,double *ecl,double *ecb,double *diamang);
+
+#endif
