Changeset 1682 in Sophya for trunk


Ignore:
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

Location:
trunk/SophyaExt/XAstroPack
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaExt/XAstroPack/xastropack.cc

    r1679 r1682  
    22#include <stdio.h>
    33#include "xastropack.h"
     4
     5// BUGS BUGS BUGS BUGS BUGS BUGS BUGS BUGS BUGS BUGS
     6// Corrections de divers bugs trouve dans libastro (CMV)
     7// 1******* In the file vsop87.c line 154:
     8// p = q/(t_abs[alpha] + alpha * t_abs[alpha-1] * 1e-4 + 1e-35);
     9// - to avoid t_abs[-1] when alpha=0, replaced by : 
     10// if(alpha>0) p = t_abs[alpha-1]; else p=0;
     11// p = q/(t_abs[alpha] + alpha * p * 1e-4 + 1e-35);
     12// Mail envoye a ecdowney@ClearSkyInstitute.com
     13// 2******* In the file eq_ecl.c line 69:
     14// *q = asin((sy*ceps)-(cy*seps*sx*sw));
     15// eq_ecl.c Protection NaN dans ecleq_aux, replaced by :
     16// *q = (sy*ceps)-(cy*seps*sx*sw);
     17// if(*q<-1.) *q = -PI/2.; else if(*q>1.) *q = PI/2.; else *q = asin(*q);
     18// Mail envoye a ecdowney@ClearSkyInstitute.com
     19
    420
    521/*!
    … …  
    1834// EQUATORIALE: ascension droite en heures [0,24[ (ra)
    1935//              declinaison en degres [-90,90]    (dec)
    20 //              angle horaire en heures [-12,12] (-12=12) (ha)
     36//              angle horaire en heures ]-12,12] (-12=12) (ha)
    2137//              temps sideral du lieu: tsid=ha+ra (ou lst)
    2238// GALACTIQUE: longitude en degres [0,360[  (glng)
    … …  
    3551
    3652/*! \ingroup XAstroPack
     53\brief Given a coordinate type "typ", convert to standard for astropack.
     54\verbatim
     55 La routine convertit (in place) les coordonnees "coord1","coord2"
     56 definies par le type "typ" dans les unites standard de ce systeme
     57 de coordonnees.
     58 "typ" code le systeme de coordonnees astronomiques et les unites utilisees
     59
     60 - Return : 0 = OK
     61            1 = Unknown type of coordinates
     62            2 = bad range for coord1
     63            4 = bad range for coord2
     64            6 = bad range for coord1 et coord2
     65
     66 Les types de coordonnees sont definies dans le enum TypAstroCoord:
     67  La premiere coordonnee est de type "longitude" (alpha,longitude)
     68  La deuxieme coordonnee est de type "latidude" (delta,latitude)
     69 *** Definitions des unites des coordonnees et de leurs dynamiques
     70 - TypCoordH0 : heure=[0,24[
     71 - TypCoordH1 : heure=]-12,12]
     72 - TypCoordD0 : degre=[0,360[
     73 - TypCoordD1 : degre=]-180,180]
     74 - TypCoordD2 : degre=[-90,90]
     75 - TypCoordR0 : degre=[0,2Pi[
     76 - TypCoordR1 : degre=]-Pi,Pi]
     77 - TypCoordR2 : degre=[-Pi/2,Pi/2]
     78 *** Definitions des combinaisons unites des coordonnees
     79 - TypCoordHD  : coordonnees en (heure=[0,24[,degre=[-90,90])
     80 - TypCoordDD  : coordonnees en (degre=[0,360[,degre=[-90,90])
     81 - TypCoordRR  : coordonnees en (radian=[0,2Pi[,radian=[-Pi/2,Pi/2])
     82 - TypCoordH1D : coordonnees en (heure=]-12,12],degre=[-90,90])
     83 - TypCoordD1D : coordonnees en (degre=]-180,180],degre=[-90,90])
     84 - TypCoordR1R : coordonnees en (radian=]-Pi,Pi],radian=[-Pi/2,Pi/2])
     85 *** Definitions des types de systemes de coordonnees astronomiques.
     86 - TypCoordEq  : Coordonnees Equatoriales alpha,delta
     87 - TypCoordGal : Coordonnees Galactiques gLong, gLat
     88 - TypCoordHor : Coordonnees Horizontales azimuth,altitude
     89 - TypCoordEcl : Coordonnees Ecliptiques EclLong,EclLat
     90 *** Definitions des unites par defaut pour les divers systemes de coordonnees.
     91 - TypCoordEqStd  : heure=[0,24[, degre=[-90,90]
     92 - TypCoordGalStd : degre=[0,360[,degre=[-90,90]
     93 - TypCoordHorStd : degre=[0,360[,degre=[-90,90]
     94 - TypCoordEclStd : degre=[0,360[,degre=[-90,90]
     95\endverbatim
     96*/
     97int  CoordConvertToStd(TypAstroCoord typ,double& coord1,double& coord2)
     98{
     99  int rc = 0;
     100
     101  // ---- Equatoriales alpha,delta
     102  //    - standard = [0,24[ , [-90,90]
     103  if(typ&TypCoordEq) {
     104    if(typ&TypCoordDD) {
     105      coord1 = deghr(coord1);
     106    } else if(typ&TypCoordRR) {
     107      coord1 = radhr(coord1);
     108      coord2 = raddeg(coord2);
     109    }
     110    if(coord1==24.) coord1 = 0.;
     111    if(coord1<0.   || coord1>=24.) rc+= 2;
     112    if(coord2<-90. || coord2>90. ) rc+= 4;
     113
     114  // ---- Galactiques gLong, gLat
     115  // ---- Horizontales azimuth,altitude
     116  // ---- Ecliptiques EclLong,EclLat
     117  //    - standard = [0,360[ , [-90,90]
     118  } else if( typ&TypCoordGal || typ&TypCoordHor || typ&TypCoordEcl) {
     119    if(typ&TypCoordHD) {
     120      coord1 = hrdeg(coord1);
     121    } else if(typ&TypCoordRR) {
     122      coord1 = raddeg(coord1);
     123      coord2 = raddeg(coord2);
     124    }
     125    if(coord1==360.) coord1 = 0.;
     126    if(coord1<0.   || coord1>=360.) rc+= 2;
     127    if(coord2<-90. || coord2>90. )  rc+= 4;
     128
     129  } else {          // Coordonnees non-connues
     130    rc= 1;
     131  }
     132
     133  return rc;
     134}
     135
     136/*! \ingroup XAstroPack
     137\brief Compute MJD from date
     138\verbatim
     139 MJD =  modified Julian date (number of days elapsed since 1900 jan 0.5),
     140\endverbatim
     141*/
     142double MJDfrDate(double dy,int mn,int yr)
     143{
     144double mjd;
     145cal_mjd(mn,dy,yr,&mjd);
     146return mjd;
     147}
     148
     149/*! \ingroup XAstroPack
     150\brief Compute date from MJD
     151*/
     152void DatefrMJD(double mjd,double *dy,int *mn,int *yr)
     153{
     154mjd_cal(mjd,mn,dy,yr);
     155}
     156
     157/*! \ingroup XAstroPack
     158\brief  Given a mjd, return the year as a double.
     159*/
     160double YearfrMJD(double mjd)
     161{
     162double yr;
     163mjd_year(mjd,&yr);
     164return yr;
     165}
     166
     167/*! \ingroup XAstroPack
     168\brief Given a decimal year, return mjd
     169*/
     170double MJDfrYear(double yr)
     171{
     172double mjd;
     173year_mjd(yr,&mjd);
     174return mjd;
     175}
     176
     177/*! \ingroup XAstroPack
     178\brief Given a mjd, return the year and number of days since 00:00 Jan 1
     179\warning: if mjd = 2 January -> number of days = 1
     180*/
     181void YDfrMJD(double mjd,double *dy,int *yr)
     182{
     183mjd_dayno(mjd,yr,dy);
     184}
     185
     186/*! \ingroup XAstroPack
     187\brief Given a year,
     188*/
     189int IsLeapYear(int y)
     190{
     191return isleapyear(y);
     192}
     193
     194/*! \ingroup XAstroPack
     195\brief given an mjd, set *dow to 0..6 according to which day of the week it falls on (0=sunday).
     196\return return 0 if ok else -1 if can't figure it out.
     197*/
     198int DayOrder(double mjd,int *dow)
     199{
     200return mjd_dow(mjd,dow);
     201}
     202
     203/*! \ingroup XAstroPack
     204\brief given a mjd, return the the number of days in the month.
     205*/
     206int DaysInMonth(double mjd)
     207{
     208int ndays;
     209mjd_dpm(mjd,&ndays);
     210return ndays;
     211}
     212
     213/*! \ingroup XAstroPack
     214\brief Given a mjd, truncate it to the beginning of the whole day
     215*/
     216double MJDat0hFrMJD(double mjd)
     217{
     218return mjd_day(mjd);
     219}
     220
     221/*! \ingroup XAstroPack
     222\brief Given a mjd, return the number of hours past midnight of the whole day
     223*/
     224double HfrMJD(double mjd)
     225{
     226return mjd_hr(mjd);
     227}
     228
     229/*! \ingroup XAstroPack
     230\brief Give GST from UTC
     231\verbatim
     232 Given a modified julian date, mjd, and a universally coordinated time, utc,
     233 return greenwich mean siderial time, *gst.
     234 N.B. mjd must be at the beginning of the day.
     235\endverbatim
     236*/
     237double GSTfrUTC(double mjd0,double utc)
     238{
     239double gst;
     240utc_gst(mjd0,utc,&gst);
     241return gst;
     242}
     243
     244/*! \ingroup XAstroPack
     245\brief Give UTC from GST
     246\verbatim
     247 Given a modified julian date, mjd, and a greenwich mean siderial time, gst,
     248 return universally coordinated time, *utc.
     249 N.B. mjd must be at the beginning of the day.
     250\endverbatim
     251*/                                                                             
     252double UTCfrGST(double mjd0,double gst)
     253{
     254double utc;
     255gst_utc(mjd0,gst,&utc);
     256return utc;
     257}
     258
     259/*! \ingroup XAstroPack
     260\brief Given apparent altitude find airmass.
     261*/                                                                             
     262double AirmassfrAlt(double alt)
     263{
     264double x;
     265alt = degrad(alt);
     266airmass(alt,&x);
     267return x;
     268}
     269
     270/*! \ingroup XAstroPack
     271\brief given geocentric time "jd" and coords of a distant object at "ra/dec" (J2000),
     272find the difference "hcp" in time between light arriving at earth vs the sun.
     273\return "hcp" must be subtracted from "geocentric jd" to get "heliocentric jd".
     274\warning "jd" is the TRUE Julian day (jd = mjd+MJD0).
     275*/
     276double HelioCorr(double jd,double ra,double dec)
     277{
     278double hcp;
     279ra=hrrad(ra);
     280dec=degrad(dec);
     281heliocorr(jd,ra,dec,&hcp);
     282return hcp;
     283}
     284
     285/*! \ingroup XAstroPack
    37286\brief gmst0() - return Greenwich Mean Sidereal Time at 0h UT
    38287\param mjd0 = date at 0h UT in julian days since MJD0
    … …  
    235484\warning Correction for the effect on the angle of the obliquity due to nutation is not included.
    236485*/                                                                             
    237 // Attention, j'ai modifie eq_ecl.c pour proteger NaN
    238 // dans ecleq_aux :
    239 // *q = (sy*ceps)-(cy*seps*sx*sw);
    240 // if(*q<-1.) *q = -PI/2.; else if(*q>1.) *q = PI/2.; else *q = asin(*q);
    241486void EqtoEcl(double mjd,double ra,double dec,double *eclng,double *eclat)
    242487// Coordonnees equatoriales -> Coordonnees ecliptiques
    … …  
    336581 *sunecb = raddeg(*sunecb);
    337582}
    338 
    339 /*! \ingroup XAstroPack
    340 \brief Given a coordinate type "typ", convert to standard for astropack
    341 \verbatim
    342 // Return : 0 = OK
    343 //          1 = Unknown type of coordinates
    344 //          2 = bad range for coord1
    345 //          4 = bad range for coord2
    346 //          6 = bad range for coord1 et coord2
    347 \endverbatim
    348 */
    349 int  CoordConvertToStd(TypAstroCoord typ,double& coord1,double& coord2)
    350 {
    351   int rc = 0;
    352 
    353   // ---- Equatoriales alpha,delta
    354   //    - standard = [0,24[ , [-90,90]
    355   if(typ&TypCoordEq) {
    356     if(typ&TypCoordDD) {
    357       coord1 = deghr(coord1);
    358     } else if(typ&TypCoordRR) {
    359       coord1 = radhr(coord1);
    360       coord2 = raddeg(coord2);
    361     }
    362     if(coord1==24.) coord1 = 0.;
    363     if(coord1<0.   || coord1>=24.) rc+= 2;
    364     if(coord2<-90. || coord2>90. ) rc+= 4;
    365 
    366   // ---- Galactiques gLong, gLat
    367   // ---- Horizontales azimuth,altitude
    368   // ---- Ecliptiques EclLong,EclLat
    369   //    - standard = [0,360[ , [-90,90]
    370   } else if( typ&TypCoordGal || typ&TypCoordHor || typ&TypCoordEcl) {
    371     if(typ&TypCoordHD) {
    372       coord1 = hrdeg(coord1);
    373     } else if(typ&TypCoordRR) {
    374       coord1 = raddeg(coord1);
    375       coord2 = raddeg(coord2);
    376     }
    377     if(coord1==360.) coord1 = 0.;
    378     if(coord1<0.   || coord1>=360.) rc+= 2;
    379     if(coord2<-90. || coord2>90. )  rc+= 4;
    380 
    381   } else {          // Coordonnees non-connues
    382     rc= 1;
    383   }
    384 
    385   return rc;
    386 }

  • trunk/SophyaExt/XAstroPack/xastropack.h

    r1679 r1682  
    1818  TypCoordUndef  =  (unsigned long)  (0),
    1919
     20  TypCoordH0     =  (unsigned long)  (1 << 10), // heure=[0,24[
     21  TypCoordH1     =  (unsigned long)  (1 << 11), // heure=]-12,12]
     22  TypCoordD0     =  (unsigned long)  (1 << 12), // degre=[0,360[
     23  TypCoordD1     =  (unsigned long)  (1 << 13), // degre=]-180,180]
     24  TypCoordD2     =  (unsigned long)  (1 << 14), // degre=[-90,90]
     25  TypCoordR0     =  (unsigned long)  (1 << 15), // degre=[0,2Pi[
     26  TypCoordR1     =  (unsigned long)  (1 << 16), // degre=]-Pi,Pi]
     27  TypCoordR2     =  (unsigned long)  (1 << 17), // degre=[-Pi/2,Pi/2]
     28
    2029  // Pour indiquer que les coordonnees sont en (heure=[0,24[,degre=[-90,90])
    2130  TypCoordHD     =  (unsigned long)  (1 << 20),
    … …  
    2433  // Pour indiquer que les coordonnees sont en (radian=[0,2Pi[,radian=[-Pi/2,Pi/2])
    2534  TypCoordRR     =  (unsigned long)  (1 << 22),
     35  // Pour indiquer que les coordonnees sont en (heure=]-12,12],degre=[-90,90])
     36  TypCoordH1D    =  (unsigned long)  (1 << 23),
     37  // Pour indiquer que les coordonnees sont en (degre=]-180,180],degre=[-90,90])
     38  TypCoordD1D    =  (unsigned long)  (1 << 24),
     39  // Pour indiquer que les coordonnees sont en (radian=]-Pi,Pi],radian=[-Pi/2,Pi/2])
     40  TypCoordR1R    =  (unsigned long)  (1 << 25),
    2641
    2742  // Coordonnees Equatoriales alpha,delta
    … …  
    6984inline double MJDfrTrueJD(double jd) {return jd - MJD0;}
    7085
    71 /*! \ingroup XAstroPack
    72 \brief Compute MJD from date
    73 \verbatim
    74  MJD =  modified Julian date (number of days elapsed since 1900 jan 0.5),
    75 \endverbatim
    76 */
    77 inline double MJDfrDate(double dy,int mn,int yr)
    78        {double mjd; cal_mjd(mn,dy,yr,&mjd); return mjd;}
    79 
    80 /*! \ingroup XAstroPack
    81 \brief Compute date from MJD
    82 */
    83 inline void DatefrMJD(double mjd,double *dy,int *mn,int *yr)
    84                      {mjd_cal(mjd,mn,dy,yr);}
    85 
    86 /*! \ingroup XAstroPack
    87 \brief  Given a mjd, return the year as a double.
    88 */
    89 inline double YearfrMJD(double mjd)
    90        {double yr; mjd_year(mjd,&yr); return yr;}
    91 
    92 /*! \ingroup XAstroPack
    93 \brief Given a decimal year, return mjd
    94 */
    95 inline double MJDfrYear(double yr)
    96        {double mjd; year_mjd(yr,&mjd); return mjd;}
    97 
    98 /*! \ingroup XAstroPack
    99 \brief Given a mjd, return the year and number of days since 00:00 Jan 1
    100 \warning: if mjd = 2 January -> number of days = 1
    101 */
    102 inline void YDfrMJD(double mjd,double *dy,int *yr)
    103        {mjd_dayno(mjd,yr,dy);}
    104 
    105 /*! \ingroup XAstroPack
    106 \brief Given a year,
    107 */
    108 inline int IsLeapYear(int y)
    109        {return isleapyear(y);}
    110 
    111 /*! \ingroup XAstroPack
    112 \brief given an mjd, set *dow to 0..6 according to which day of the week it falls on (0=sunday).
    113 \return return 0 if ok else -1 if can't figure it out.
    114 */
    115 inline int DayOrder(double mjd,int *dow)
    116        {return mjd_dow(mjd,dow);}
    117 
    118 /*! \ingroup XAstroPack
    119 \brief given a mjd, return the the number of days in the month.
    120 */
    121 inline int DaysInMonth(double mjd)
    122        {int ndays; mjd_dpm(mjd,&ndays); return ndays;}
    123 
    124 /*! \ingroup XAstroPack
    125 \brief Given a mjd, truncate it to the beginning of the whole day
    126 */
    127 inline double MJDat0hFrMJD(double mjd) {return mjd_day(mjd);}
    128 
    129 /*! \ingroup XAstroPack
    130 \brief Given a mjd, return the number of hours past midnight of the whole day
    131 */
    132 inline double HfrMJD(double mjd) {return mjd_hr(mjd);}
    133 
    134 /*! \ingroup XAstroPack
    135 \brief Give GST from UTC
    136 \verbatim
    137  Given a modified julian date, mjd, and a universally coordinated time, utc,
    138  return greenwich mean siderial time, *gst.
    139  N.B. mjd must be at the beginning of the day.
    140 \endverbatim
    141 */
    142 inline double GSTfrUTC(double mjd0,double utc)
    143        {double gst; utc_gst(mjd0,utc,&gst) ; return gst;}
    144 
    145 /*! \ingroup XAstroPack
    146 \brief Give UTC from GST
    147 \verbatim
    148  Given a modified julian date, mjd, and a greenwich mean siderial time, gst,
    149  return universally coordinated time, *utc.
    150  N.B. mjd must be at the beginning of the day.
    151 \endverbatim
    152 */                                                                             
    153 inline double UTCfrGST(double mjd0,double gst)
    154        {double utc; gst_utc(mjd0,gst,&utc); return utc;}
     86double MJDfrDate(double dy,int mn,int yr);
     87void DatefrMJD(double mjd,double *dy,int *mn,int *yr);
     88double YearfrMJD(double mjd);
     89double MJDfrYear(double yr);
     90void YDfrMJD(double mjd,double *dy,int *yr);
     91int IsLeapYear(int y);
     92int DayOrder(double mjd,int *dow);
     93int DaysInMonth(double mjd);
     94double MJDat0hFrMJD(double mjd);
     95double HfrMJD(double mjd);
     96double GSTfrUTC(double mjd0,double utc);
     97double UTCfrGST(double mjd0,double gst);
    15598
    15699/*! \ingroup XAstroPack
    … …  
    171114// ------------------- Calculs Divers -------------------
    172115void Precess(double mjd1,double mjd2,double ra1,double dec1,double *ra2,double *dec2);
    173 
    174 /*! \ingroup XAstroPack
    175 \brief Given apparent altitude find airmass.
    176 */                                                                             
    177 inline double AirmassfrAlt(double alt)
    178        {double x; alt = degrad(alt); airmass(alt,&x); return x;}
    179 
    180 /*! \ingroup XAstroPack
    181 \brief given geocentric time "jd" and coords of a distant object at "ra/dec" (J2000),
    182 find the difference "hcp" in time between light arriving at earth vs the sun.
    183 \return "hcp" must be subtracted from "geocentric jd" to get "heliocentric jd".
    184 \warning "jd" is the TRUE Julian day (jd = mjd+MJD0).
    185 */
    186 inline double HelioCorr(double jd,double ra,double dec)
    187        {double hcp; ra=hrrad(ra); dec=degrad(dec); heliocorr(jd,ra,dec,&hcp); return hcp;}
     116double AirmassfrAlt(double alt);
     117double HelioCorr(double jd,double ra,double dec);
    188118
    189119// ------------------- Transformation de coordonnees -------------------
Note: See TracChangeset for help on using the changeset viewer.