Changeset 1682 in Sophya for trunk/SophyaExt/XAstroPack/xastropack.cc

Timestamp:

Oct 11, 2001, 2:42:00 PM (24 years ago)

Author:

cmv

Message:

pb avec inline contenants des routines libastro cmv 11/10/01

File:

• trunk/SophyaExt/XAstroPack/xastropack.cc (modified) (5 diffs)

trunk/SophyaExt/XAstroPack/xastropack.cc

@@ -2 +2 @@
 #include <stdio.h>
 #include "xastropack.h"
+
+// BUGS BUGS BUGS BUGS BUGS BUGS BUGS BUGS BUGS BUGS
+// Corrections de divers bugs trouve dans libastro (CMV)
+// 1******* In the file vsop87.c line 154:
+// p = q/(t_abs[alpha] + alpha * t_abs[alpha-1] * 1e-4 + 1e-35);
+// - to avoid t_abs[-1] when alpha=0, replaced by :  
+// if(alpha>0) p = t_abs[alpha-1]; else p=0;
+// p = q/(t_abs[alpha] + alpha * p * 1e-4 + 1e-35);
+// Mail envoye a ecdowney@ClearSkyInstitute.com
+// 2******* In the file eq_ecl.c line 69:
+// *q = asin((sy*ceps)-(cy*seps*sx*sw));
+// eq_ecl.c Protection NaN dans ecleq_aux, replaced by : 
+// *q = (sy*ceps)-(cy*seps*sx*sw);
+// if(*q<-1.) *q = -PI/2.; else if(*q>1.) *q = PI/2.; else *q = asin(*q);
+// Mail envoye a ecdowney@ClearSkyInstitute.com
+
 
 /*!
@@ -18 +34 @@
 // EQUATORIALE: ascension droite en heures [0,24[ (ra)
 //              declinaison en degres [-90,90]    (dec)
-//              angle horaire en heures [-12,12] (-12=12) (ha) 
+//              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)
@@ -35 +51 @@
 
 /*! \ingroup XAstroPack
+\brief Given a coordinate type "typ", convert to standard for astropack.
+\verbatim
+ La routine convertit (in place) les coordonnees "coord1","coord2"
+ definies par le type "typ" dans les unites standard de ce systeme
+ de coordonnees.
+ "typ" code le systeme de coordonnees astronomiques et les unites utilisees
+
+ - 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
+
+ Les types de coordonnees sont definies dans le enum TypAstroCoord:
+  La premiere coordonnee est de type "longitude" (alpha,longitude)
+  La deuxieme coordonnee est de type "latidude" (delta,latitude)
+ *** Definitions des unites des coordonnees et de leurs dynamiques
+ - TypCoordH0 : heure=[0,24[
+ - TypCoordH1 : heure=]-12,12]
+ - TypCoordD0 : degre=[0,360[
+ - TypCoordD1 : degre=]-180,180]
+ - TypCoordD2 : degre=[-90,90]
+ - TypCoordR0 : degre=[0,2Pi[
+ - TypCoordR1 : degre=]-Pi,Pi]
+ - TypCoordR2 : degre=[-Pi/2,Pi/2]
+ *** Definitions des combinaisons unites des coordonnees
+ - TypCoordHD  : coordonnees en (heure=[0,24[,degre=[-90,90])
+ - TypCoordDD  : coordonnees en (degre=[0,360[,degre=[-90,90])
+ - TypCoordRR  : coordonnees en (radian=[0,2Pi[,radian=[-Pi/2,Pi/2])
+ - TypCoordH1D : coordonnees en (heure=]-12,12],degre=[-90,90])
+ - TypCoordD1D : coordonnees en (degre=]-180,180],degre=[-90,90])
+ - TypCoordR1R : coordonnees en (radian=]-Pi,Pi],radian=[-Pi/2,Pi/2])
+ *** Definitions des types de systemes de coordonnees astronomiques.
+ - TypCoordEq  : Coordonnees Equatoriales alpha,delta
+ - TypCoordGal : Coordonnees Galactiques gLong, gLat
+ - TypCoordHor : Coordonnees Horizontales azimuth,altitude
+ - TypCoordEcl : Coordonnees Ecliptiques EclLong,EclLat
+ *** Definitions des unites par defaut pour les divers systemes de coordonnees.
+ - TypCoordEqStd  : heure=[0,24[, degre=[-90,90]
+ - TypCoordGalStd : degre=[0,360[,degre=[-90,90]
+ - TypCoordHorStd : degre=[0,360[,degre=[-90,90]
+ - TypCoordEclStd : degre=[0,360[,degre=[-90,90]
+\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 = deghr(coord1);
+    } else if(typ&TypCoordRR) {
+      coord1 = radhr(coord1);
+      coord2 = raddeg(coord2);
+    }
+    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 = hrdeg(coord1);
+    } else if(typ&TypCoordRR) {
+      coord1 = raddeg(coord1);
+      coord2 = raddeg(coord2);
+    }
+    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;
+}
+
+/*! \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 Given apparent altitude find airmass.
+*/                                                                             
+double AirmassfrAlt(double alt)
+{
+double x;
+alt = degrad(alt);
+airmass(alt,&x);
+return x;
+}
+
+/*! \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=hrrad(ra);
+dec=degrad(dec);
+heliocorr(jd,ra,dec,&hcp);
+return hcp;
+}
+
+/*! \ingroup XAstroPack
 \brief gmst0() - return Greenwich Mean Sidereal Time at 0h UT
 \param mjd0 = date at 0h UT in julian days since MJD0
@@ -235 +484 @@
 \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
-// 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
@@ -336 +581 @@
  *sunecb = raddeg(*sunecb);
 }
-
-/*! \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 = deghr(coord1);
-    } else if(typ&TypCoordRR) {
-      coord1 = radhr(coord1);
-      coord2 = raddeg(coord2);
-    }
-    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 = hrdeg(coord1);
-    } else if(typ&TypCoordRR) {
-      coord1 = raddeg(coord1);
-      coord2 = raddeg(coord2);
-    }
-    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;
-}