Changeset 1679 in Sophya


Ignore:
Timestamp:
Oct 10, 2001, 11:55:10 PM (24 years ago)
Author:
cmv
Message:

InRange,degrad()... et inline cmv 10/10/2001

Location:
trunk/SophyaExt/XAstroPack
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaExt/XAstroPack/xastropack.cc

    r1678 r1679  
    3535
    3636/*! \ingroup XAstroPack
    37 \brief Compute true Julian day from MJD
    38 */
    39 double TrueJDfrMJD(double mjd)
    40 {
    41   return mjd + MJD0;
    42 }
    43 
    44 /*! \ingroup XAstroPack
    45 \brief Compute MJD from true Julian day
    46 */
    47 double MJDfrTrueJD(double jd)
    48 {
    49   return jd - MJD0;
    50 }
    51 
    52 /*! \ingroup XAstroPack
    53 \brief Compute MJD from date
    54 \verbatim
    55  MJD =  modified Julian date (number of days elapsed since 1900 jan 0.5),
    56 \endverbatim
    57 */
    58 double MJDfrDate(double dy,int mn,int yr)
    59 {
    60   double mjd;
    61   cal_mjd(mn,dy,yr,&mjd);
    62   return mjd;
    63 }
    64 
    65 /*! \ingroup XAstroPack
    66 \brief Compute date from MJD
    67 */
    68 void DatefrMJD(double mjd,double *dy,int *mn,int *yr)
    69 {
    70   mjd_cal(mjd,mn,dy,yr);
    71 }
    72 
    73 /*! \ingroup XAstroPack
    74 \brief  Given a mjd, return the year as a double.
    75 */
    76 double YearfrMJD(double mjd)
    77 {
    78   double yr;
    79   mjd_year(mjd,&yr);
    80   return yr;
    81 }
    82 
    83 /*! \ingroup XAstroPack
    84 \brief Given a decimal year, return mjd
    85 */
    86 double MJDfrYear(double yr)
    87 {
    88   double mjd;
    89   year_mjd(yr,&mjd);
    90   return mjd;
    91 }
    92 
    93 /*! \ingroup XAstroPack
    94 \brief Given a mjd, return the year and number of days since 00:00 Jan 1
    95 \warning: if mjd = 2 January -> number of days = 1
    96 */
    97 void YDfrMJD(double mjd,double *dy,int *yr)
    98 {
    99   mjd_dayno(mjd,yr,dy);
    100 }
    101 
    102 /*! \ingroup XAstroPack
    103 \brief Given a year,
    104 */
    105 int IsLeapYear(int y)
    106 {
    107   return isleapyear(y);
    108 }
    109 
    110 /*! \ingroup XAstroPack
    111 \brief given an mjd, set *dow to 0..6 according to which day of the week it falls on (0=sunday).
    112 \return return 0 if ok else -1 if can't figure it out.
    113 */
    114 int DayOrder(double mjd,int *dow)
    115 {
    116   return mjd_dow(mjd,dow);
    117 }
    118 
    119 /*! \ingroup XAstroPack
    120 \brief given a mjd, return the the number of days in the month.
    121 */
    122 int DaysInMonth(double mjd)
    123 {
    124   int ndays;
    125   mjd_dpm(mjd,&ndays);
    126   return ndays;
    127 }
    128 
    129 /*! \ingroup XAstroPack
    130 \brief Given a mjd, truncate it to the beginning of the whole day
    131 */
    132 double MJDat0hFrMJD(double mjd)
    133 {
    134   return mjd_day(mjd);
    135 }
    136 
    137 /*! \ingroup XAstroPack
    138 \brief Given a mjd, return the number of hours past midnight of the whole day
    139 */
    140 double HfrMJD(double mjd)
    141 {
    142   return mjd_hr(mjd);
    143 }
    144 
    145 /*! \ingroup XAstroPack
    146 \brief Give GST from UTC
    147 \verbatim
    148  Given a modified julian date, mjd, and a universally coordinated time, utc,
    149  return greenwich mean siderial time, *gst.
    150  N.B. mjd must be at the beginning of the day.
    151 \endverbatim
    152 */
    153 double GSTfrUTC(double mjd0,double utc)
    154 {
    155   double gst;
    156   utc_gst(mjd0,utc,&gst) ;
    157   return gst;
    158 }
    159 
    160 /*! \ingroup XAstroPack
    161 \brief Give UTC from GST
    162 \verbatim
    163  Given a modified julian date, mjd, and a greenwich mean siderial time, gst,
    164  return universally coordinated time, *utc.
    165  N.B. mjd must be at the beginning of the day.
    166 \endverbatim
    167 */                                                                             
    168 double UTCfrGST(double mjd0,double gst)
    169 {
    170   double utc;
    171   gst_utc(mjd0,gst,&utc);
    172   return utc;
    173 }
    174 
    175 /*! \ingroup XAstroPack
    17637\brief gmst0() - return Greenwich Mean Sidereal Time at 0h UT
    177 \param mjd = date at 0h UT in julian days since MJD0
     38\param mjd0 = date at 0h UT in julian days since MJD0
    17839*/                                                                             
    17940double GST0(double mjd0)
    … …  
    19051
    19152/*! \ingroup XAstroPack
    192 \brief return local sidereal time from greenwich mean siderial time and longitude
    193 \param precis : if not zero, then correct for obliquity and nutation
    194 \warning no nutation or obliquity correction are done.
    195 */                                                                             
    196 double LSTfrGST(double gst,double geolng)
    197 {
    198   double lst = gst + geolng *12./180.;
     53\brief return local sidereal time from modified julian day and longitude
     54\warning nutation or obliquity correction are taken into account.
     55*/                                                                             
     56double LSTfrMJD(double mjd,double geolng)
     57{
     58  double eps,lst,deps,dpsi;
     59  utc_gst(mjd_day(mjd),mjd_hr(mjd),&lst);
     60  lst += deghr(geolng);
     61  obliquity(mjd,&eps);
     62  nutation(mjd,&deps,&dpsi);
     63  lst += radhr(dpsi*cos(eps+deps));
    19964  InRange(&lst,24.);
    20065  return lst;
    201 }
    202 
    203 /*! \ingroup XAstroPack
    204 \brief return local sidereal time from modified julian day and longitude
    205 \warning nutation or obliquity correction are taken into account.
    206 */                                                                             
    207 double LSTfrMJD(double mjd,double geolng)
    208 {
    209   double eps,lst,deps,dpsi;
    210   utc_gst(mjd_day(mjd),mjd_hr(mjd),&lst);
    211   lst += geolng *12./180.;
    212   obliquity(mjd,&eps);
    213   nutation(mjd,&deps,&dpsi);
    214   lst += dpsi*cos(eps+deps) *12./M_PI;
    215   return lst;
    216 }
    217 
    218 /*! \ingroup XAstroPack
    219 \brief Compute precession between 2 dates.
    220 */                                                                             
    221 void Precess(double mjd1,double mjd2,double ra1,double dec1,double *ra2,double *dec2)
    222 {
    223  ra1  *= PI/12.;   // radians
    224  dec1 *= PI/180.;  // radians
    225  precess(mjd1,mjd2,&ra1,&dec1);
    226  *ra2 = ra1*12./PI; InRange(ra2,24.);
    227  *dec2 = dec1*180./PI;
    228 }
    229 
    230 /*! \ingroup XAstroPack
    231 \brief Given apparent altitude find airmass.
    232 */                                                                             
    233 double AirmassfrAlt(double alt)
    234 {
    235   double x;
    236   alt *= PI/180.;  // radians
    237   airmass(alt,&x);
    238   return x;
    239 }
    240 
    241 /*! \ingroup XAstroPack
    242 \brief Give the hour angle from local sideral time and right ascencion
    243 \warning right ascencion should be first precessed to date of interest
    244 \warning no nutation or obliquity correction are done.
    245 */                                                                             
    246 double HafrRaTS(double lst,double ra)
    247 {
    248   double ha = lst - ra;
    249   // Attention au probleme de la discontinuite 0h <==> 24h
    250   // ts=1  ra=23 ; (ts-ra)=-22 <-12 --> ha = +2 = +24 + (ts-ra)
    251   // ts=23 ra=1  ; (ts-ra)=+22 >+12 --> ha = -2 = -24 + (ts-ra)
    252   InRange(&ha,24.,12.);
    253   return ha;
    254 }
    255 
    256 /*! \ingroup XAstroPack
    257 \brief Give the local sideral time and the hour angle return the right ascencion
    258 \warning right ascencion is the value precessed to date of interest
    259 \warning no nutation or obliquity correction are done.
    260 */                                                                             
    261 double RafrHaTS(double lst,double ha)
    262 {
    263   double ra = lst - ha;
    264   InRange(&ra,24.);
    265   return ra;
    266 }
    267 
    268 /*! \ingroup XAstroPack
    269 \brief given geocentric time "jd" and coords of a distant object at "ra/dec" (J2000),
    270 find the difference "hcp" in time between light arriving at earth vs the sun.
    271 \return "hcp" must be subtracted from "geocentric jd" to get "heliocentric jd".
    272 \warning "jd" is the TRUE Julian day (jd = mjd+MJD0).
    273 */
    274 double HelioCorr(double jd,double ra,double dec)
    275 {
    276  double hcp;
    277  ra  *= PI/12.;   // radians
    278  dec *= PI/180.;  // radians
    279  heliocorr(jd,ra,dec,&hcp);
    280  return hcp;
    28166}
    28267
    … …  
    350135
    351136/*! \ingroup XAstroPack
     137\brief Compute precession between 2 dates.
     138*/                                                                             
     139void Precess(double mjd1,double mjd2,double ra1,double dec1,double *ra2,double *dec2)
     140{
     141 ra1  = hrrad(ra1);   // radians
     142 dec1 = degrad(dec1);  // radians
     143 precess(mjd1,mjd2,&ra1,&dec1);
     144 *ra2 = radhr(ra1); InRange(ra2,24.);
     145 *dec2 = raddeg(dec1);
     146}
     147
     148/*! \ingroup XAstroPack
    352149\brief Convert equatorial coordinates for the given epoch into galactic coordinates
    353150*/                                                                             
    … …  
    355152// Coordonnees equatoriales -> Coordonnees galactiques
    356153{
    357  ra  *= PI/12.;   // radians
    358  dec *= PI/180.;  // radians
     154 ra  = hrrad(ra);   // radians
     155 dec = degrad(dec);  // radians
    359156 eq_gal(mjd,ra,dec,glat,glng);
    360157 // Vraiment bizarre, sur Linux-g++ glng>=360 ne comprend pas glng==360 ! (CMV)
    361  *glng *= 180./PI; InRange(glng,360.);
    362  *glat *= 180./PI;
     158 *glng = raddeg(*glng); InRange(glng,360.);
     159 *glat = raddeg(*glat);
    363160}
    364161
    … …  
    369166// Coordonnees galactiques -> Coordonnees equatoriales
    370167{
    371  glng  *= PI/180.;  // radians
    372  glat *= PI/180.;  // radians
     168 glng = degrad(glng);  // radians
     169 glat = degrad(glat);  // radians
    373170 gal_eq (mjd,glat,glng,ra,dec);
    374  *ra  *= 12./PI; InRange(ra,24.);
    375  *dec *= 180./PI;
     171 *ra = radhr(*ra); InRange(ra,24.);
     172 *dec = raddeg(*dec);
    376173}
    377174
    … …  
    382179// Coordonnees equatoriales -> Coordonnees horizontales
    383180{
    384  geolat *= PI/180.;
    385  ha  *= PI/12.;   // radians
    386  dec *= PI/180.;  // radians
     181 geolat = degrad(geolat);  // radians
     182 ha  = hrrad(ha);   // radians
     183 dec = degrad(dec);  // radians
    387184 hadec_aa (geolat,ha,dec,alt,az);
    388  *alt *= 180./PI;
    389  *az  *= 180./PI; InRange(az,360.);
     185 *alt = raddeg(*alt);
     186 *az  = raddeg(*az); InRange(az,360.);
    390187}
    391188
    … …  
    396193// Coordonnees horizontales -> Coordonnees equatoriales
    397194{
    398  geolat *= PI/180.;
    399  alt *= PI/180.;  // radians
    400  az  *= PI/180.;  // radians
     195 geolat = degrad(geolat);  // radians
     196 alt = degrad(alt);  // radians
     197 az  = degrad(az);  // radians
    401198 aa_hadec (geolat,alt,az,ha,dec);
    402  *ha  *= 12./PI; InRange(ha,24.,12.);
    403  *dec *= 180./PI;
     199 *ha = radhr(*ha); InRange(ha,24.,12.);
     200 *dec = raddeg(*dec);
    404201}
    405202
    … …  
    410207// Coordonnees equatoriales -> Coordonnees horizontales
    411208{
    412  double ha = lst - ra;
    413  if(ha==-12.) ha=12.; InRange(&ha,24.,12.);
    414  geolat *= PI/180.;
    415  ha  *= PI/12.;   // radians
    416  dec *= PI/180.;  // radians
     209 double ha = lst - ra; InRange(&ha,24.,12.);
     210 geolat = degrad(geolat);  // radians
     211 ha = hrrad(ha);   // radians
     212 dec = degrad(dec);  // radians
    417213 hadec_aa (geolat,ha,dec,alt,az);
    418  *alt *= 180./PI;
    419  *az  *= 180./PI; InRange(az,360.);
     214 *alt = raddeg(*alt);
     215 *az  = raddeg(*az); InRange(az,360.);
    420216}
    421217
    … …  
    427223{
    428224 double ha;
    429  geolat *= PI/180.;
    430  alt *= PI/180.;  // radians
    431  az  *= PI/180.;  // radians
     225 geolat = degrad(geolat);  // radians
     226 alt = degrad(alt);  // radians
     227 az  = degrad(az);  // radians
    432228 aa_hadec (geolat,alt,az,&ha,dec);
    433  ha *= 12./PI;
    434  *ra = lst - ha; InRange(ra,24.);
    435  *dec *= 180./PI;
     229 *ra = lst - radhr(ha); InRange(ra,24.);
     230 *dec = raddeg(*dec);
    436231}
    437232
    … …  
    447242// Coordonnees equatoriales -> Coordonnees ecliptiques
    448243{
    449  ra  *= PI/12.;   // radians
    450  dec *= PI/180.;  // radians
     244 ra = hrrad(ra);   // radians
     245 dec = degrad(dec);  // radians
    451246 eq_ecl(mjd,ra,dec,eclat,eclng);
    452  *eclng *= 180./PI; InRange(eclng,360.);
    453  *eclat *= 180./PI;
     247 *eclng = raddeg(*eclng); InRange(eclng,360.);
     248 *eclat = raddeg(*eclat);
    454249}
    455250
    … …  
    461256// Coordonnees ecliptiques -> Coordonnees equatoriales
    462257{
    463  eclat *= PI/180.;  // radians
    464  eclng *= PI/180.;  // radians
     258 eclat = degrad(eclat);  // radians
     259 eclng = degrad(eclng);  // radians
    465260 ecl_eq(mjd,eclat,eclng,ra,dec);
    466  *ra  *= 12./PI; InRange(ra,24.);
    467  *dec *= 180./PI;
     261 *ra = radhr(*ra); InRange(ra,24.);
     262 *dec = raddeg(*dec);
    468263}
    469264
    … …  
    472267\verbatim
    473268 given the modified JD, mjd, return the true geocentric ecliptic longitude
    474  of the sun for the mean equinox of the date, *lsn, in radians, the
    475  sun-earth distance, *rsn, in AU, and the latitude *bsn, in radians
     269 of the sun for the mean equinox of the date, *eclsn, in degres, the
     270 sun-earth distance, *rsn, in AU, and the latitude *ecbsn, in degres
    476271 (since this is always <= 1.2 arcseconds, in can be neglected by
    477  calling with bsn = NULL).
    478 \endverbatim
    479 */
    480 void SunPos(double mjd,double *eclsn,double *ecbsn)
    481 {
    482  double rsn;
    483  sunpos(mjd,eclsn,&rsn,ecbsn);
    484  *eclsn *= 180./PI; InRange(eclsn,360.);
    485  *ecbsn *= 180./PI;
     272 calling with ecbsn = NULL).
     273 - REMARQUE:
     274 * if the APPARENT ecliptic longitude is required, correct the longitude for
     275 *   nutation to the true equinox of date and for aberration (light travel time,
     276 *   approximately  -9.27e7/186000/(3600*24*365)*2*pi = -9.93e-5 radians).
     277\endverbatim
     278*/
     279void SunPos(double mjd,double *eclsn,double *ecbsn,double *rsn)
     280{
     281 sunpos(mjd,eclsn,rsn,ecbsn);
     282 *eclsn = raddeg(*eclsn); InRange(eclsn,360.);
     283 if(ecbsn!=NULL) *ecbsn = raddeg(*ecbsn);
    486284}
    487285
    … …  
    495293\endverbatim
    496294*/
    497 void MoonPos(double mjd,double *eclmn,double *ecbmn)
    498 {
    499   double rho,msp,mdp;
    500   moon(mjd,eclmn,ecbmn,&rho,&msp,&mdp);
    501  *eclmn *= 180./PI; InRange(eclmn,360.);
    502  *ecbmn *= 180./PI;
     295void MoonPos(double mjd,double *eclmn,double *ecbmn,double *rho)
     296{
     297  double msp,mdp;
     298  moon(mjd,eclmn,ecbmn,rho,&msp,&mdp);
     299 *eclmn = raddeg(*eclmn); InRange(eclmn,360.);
     300 *ecbmn = raddeg(*ecbmn);
    503301}
    504302
    … …  
    507305\verbatim
    508306 * given a modified Julian date, mjd, and a planet, p, find:
    509  *   lpd0: heliocentric longitude,
    510  *   psi0: heliocentric latitude,
    511  *   rp0:  distance from the sun to the planet,
    512  *   rho0: distance from the Earth to the planet,
     307 *   sunecl: heliocentric longitude,
     308 *   sunecb: heliocentric latitude,
     309 *   sundist:  distance from the sun to the planet,
     310 *   geodist: distance from the Earth to the planet,
    513311 *         none corrected for light time, ie, they are the true values for the
    514312 *         given instant.
    515  *   lam:  geocentric ecliptic longitude,
    516  *   bet:  geocentric ecliptic latitude,
     313 *   geoecl:  geocentric ecliptic longitude,
     314 *   geoecb:  geocentric ecliptic latitude,
    517315 *         each corrected for light time, ie, they are the apparent values as
    518316 *         seen from the center of the Earth for the given instant.
    519  *   dia:  angular diameter in arcsec at 1 AU,
     317 *   diamang:  angular diameter in arcsec at 1 AU,
    520318 *   mag:  visual magnitude when 1 AU from sun and earth at 0 phase angle.
    521319 *   (for the mean equinox)
     320 * all angles are in degres, all distances in AU.
     321 *
     322 * corrections for nutation and abberation must be made by the caller. The RA
     323 *   and DEC calculated from the fully-corrected ecliptic coordinates are then
     324 *   the apparent geocentric coordinates. Further corrections can be made, if
     325 *   required, for atmospheric refraction and geocentric parallax.
    522326 \endverbatim
    523327*/
    524 void PlanetPos(double mjd,int numplan,double *ecl,double *ecb,double *diamang)
    525 {
    526  double lpd0,psi0,rp0,rho0,mag;
    527  plans(mjd,numplan,&lpd0,&psi0,&rp0,&rho0,ecl,ecb,diamang,&mag);
    528  *ecl *= 180./PI; InRange(ecl,360.);
    529  *ecb *= 180./PI;
    530 }
    531 
    532 /*! \ingroup XAstroPack
    533 \brief Give Jupiter position
    534 */
    535 void JupiterPos(double mjd,double *ecl,double *ecb,double *diamang)
    536 {
    537   PlanetPos(mjd,JUPITER,ecl,ecb,diamang);
    538 }
    539 
    540 /*! \ingroup XAstroPack
    541 \brief Give Saturn position
    542 */
    543 void SaturnPos(double mjd,double *ecl,double *ecb,double *diamang)
    544 {
    545   PlanetPos(mjd,SATURN,ecl,ecb,diamang);
     328void PlanetPos(double mjd,int numplan,double *sunecl,double *sunecb,double *sundist
     329              ,double *geodist,double *geoecl,double *geoecb
     330              ,double *diamang,double *mag)
     331{
     332 plans(mjd,numplan,sunecl,sunecb,sundist,geodist,geoecl,geoecb,diamang,mag);
     333 *geoecl = raddeg(*geoecl); InRange(geoecl,360.);
     334 *geoecb = raddeg(*geoecb);
     335 *sunecl = raddeg(*sunecl); InRange(sunecl,360.);
     336 *sunecb = raddeg(*sunecb);
    546337}
    547338
    … …  
    564355  if(typ&TypCoordEq) {
    565356    if(typ&TypCoordDD) {
    566       coord1 = coord1 / 180. * 12.;
     357      coord1 = deghr(coord1);
    567358    } else if(typ&TypCoordRR) {
    568       coord1 = coord1 / PI * 12.;
    569       coord2 = coord2 / PI * 180.;
     359      coord1 = radhr(coord1);
     360      coord2 = raddeg(coord2);
    570361    }
    571362    if(coord1==24.) coord1 = 0.;
    … …  
    579370  } else if( typ&TypCoordGal || typ&TypCoordHor || typ&TypCoordEcl) {
    580371    if(typ&TypCoordHD) {
    581       coord1 = coord1 / 12. * 180.;
     372      coord1 = hrdeg(coord1);
    582373    } else if(typ&TypCoordRR) {
    583       coord1 = coord1 / PI * 180.;
    584       coord2 = coord2 / PI * 180.;
     374      coord1 = raddeg(coord1);
     375      coord2 = raddeg(coord2);
    585376    }
    586377    if(coord1==360.) coord1 = 0.;

  • trunk/SophyaExt/XAstroPack/xastropack.h

    r1678 r1679  
    3939};
    4040
    41 double TrueJDfrMJD(double mjd);
    42 double MJDfrTrueJD(double jd);
    43 double MJDfrDate(double dy,int mn,int yr);
    44 void DatefrMJD(double mjd,double *dy,int *mn,int *yr);
    45 double YearfrMJD(double mjd);
    46 double MJDfrYear(double yr);
    47 void YDfrMJD(double mjd,double *dy,int *yr);
    48 int IsLeapYear(int y);
    49 int DayOrder(double mjd,int *dow);
    50 int DaysInMonth(double mjd);
    51 double MJDat0hFrMJD(double mjd);
    52 double HfrMJD(double mjd);
    53 double GSTfrUTC(double mjd0,double utc);
    54 double UTCfrGST(double mjd0,double gst);
     41// ------------------- Utilitaires -------------------
     42int  CoordConvertToStd(TypAstroCoord typ,double& coord1,double& coord2);
     43
     44/*!
     45\brief Pour remettre la valeur "val" dans la dynamique [0.,range[.
     46       Si "vmax" different de "range", c'est la borne superieure
     47       qui peut etre atteinte
     48       (si elle est depassee, on soustrait "range").
     49\verbatim
     50  r>0 vmax>0
     51  r=24. vmax=24.   -> mets dans [  0,+24[ borne sup exclue
     52  r=24. vmax=12.   -> mets dans ]-12,+12] borne inf exclue
     53\endverbatim
     54*/
     55inline void InRange(double *val,double range)
     56             {*val-=range*floor(*val/range);}
     57inline void InRange(double *val,double range,double vmax)
     58         {InRange(val,range); if(*val>vmax) *val-=range;}
     59
     60// ------------------- Gestions des temps -------------------
     61/*! \ingroup XAstroPack
     62\brief Compute true Julian day from MJD
     63*/
     64inline double TrueJDfrMJD(double mjd) {return mjd + MJD0;}
     65
     66/*! \ingroup XAstroPack
     67\brief Compute MJD from true Julian day
     68*/
     69inline double MJDfrTrueJD(double jd) {return jd - MJD0;}
     70
     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*/
     77inline 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*/
     83inline 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*/
     89inline 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*/
     95inline 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*/
     102inline void YDfrMJD(double mjd,double *dy,int *yr)
     103       {mjd_dayno(mjd,yr,dy);}
     104
     105/*! \ingroup XAstroPack
     106\brief Given a year,
     107*/
     108inline 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*/
     115inline 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*/
     121inline 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*/
     127inline 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*/
     132inline 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*/
     142inline 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*/                                                                             
     153inline double UTCfrGST(double mjd0,double gst)
     154       {double utc; gst_utc(mjd0,gst,&utc); return utc;}
     155
     156/*! \ingroup XAstroPack
     157\brief return local sidereal time from greenwich mean siderial time and longitude
     158\param precis : if not zero, then correct for obliquity and nutation
     159\warning no nutation or obliquity correction are done.
     160*/                                                                             
     161inline double LSTfrGST(double gst,double geolng)
     162       {double lst = gst + deghr(geolng); InRange(&lst,24.); return lst;}
     163
    55164double GST0(double mjd0);
    56 double LSTfrGST(double gst,double geolng);
    57165double LSTfrMJD(double mjd,double geolng);
    58 void Precess(double mjd1,double mjd2,double ra1,double dec1,double *ra2,double *dec2);
    59 double AirmassfrAlt(double alt);
    60 double HafrRaTS(double lst,double ra);
    61 double RafrHaTS(double lst,double ha);
    62 double HelioCorr(double jd,double ra,double dec);
    63166void HMSfrHdec(double hd,int *h,int *mn,double *s);
    64167double HdecfrHMS(int h,int mn,double s);
    65168string ToStringHMS(int h,int mn,double s);
    66169string ToStringHdec(double hd);
     170
     171// ------------------- Calculs Divers -------------------
     172void 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*/                                                                             
     177inline 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),
     182find 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*/
     186inline double HelioCorr(double jd,double ra,double dec)
     187       {double hcp; ra=hrrad(ra); dec=degrad(dec); heliocorr(jd,ra,dec,&hcp); return hcp;}
     188
     189// ------------------- Transformation de coordonnees -------------------
     190/*! \ingroup XAstroPack
     191\brief Give the hour angle from local sideral time and right ascencion
     192\warning right ascencion should be first precessed to date of interest
     193\warning no nutation or obliquity correction are done.
     194*/                                                                             
     195// Attention au probleme de la discontinuite 0h <==> 24h
     196// ts=1  ra=23 ; (ts-ra)=-22 <-12 --> ha = +2 = +24 + (ts-ra)
     197// ts=23 ra=1  ; (ts-ra)=+22 >+12 --> ha = -2 = -24 + (ts-ra)
     198inline double HafrRaTS(double lst,double ra)
     199       {double ha = lst - ra; InRange(&ha,24.,12.); return ha;}
     200
     201/*! \ingroup XAstroPack
     202\brief Give the local sideral time and the hour angle return the right ascencion
     203\warning right ascencion is the value precessed to date of interest
     204\warning no nutation or obliquity correction are done.
     205*/                                                                             
     206inline double RafrHaTS(double lst,double ha)
     207       {double ra = lst - ha; InRange(&ra,24.); return ra;}
     208
    67209void EqtoGal(double mjd,double ra,double dec,double *glng,double *glat);
    68210void GaltoEq(double mjd,double glng,double glat,double *ra,double *dec);
    … …  
    73215void EqtoEcl(double mjd,double ra,double dec,double *eclng,double *eclat);
    74216void EcltoEq(double mjd,double eclng,double eclat,double *ra,double *dec);
    75 void SunPos(double mjd,double *lsn,double *bsn);
    76 void MoonPos(double mjd,double *lmn,double *bmn);
    77 void PlanetPos(double mjd,int numplan,double *ecl,double *ecb,double *diamang);
    78 void JupiterPos(double mjd,double *ecl,double *ecb,double *diamang);
    79 void SaturnPos(double mjd,double *ecl,double *ecb,double *diamang);
    80 int  CoordConvertToStd(TypAstroCoord typ,double& coord1,double& coord2);
    81 
    82 /*!
    83 \brief Pour remettre la valeur "val" dans la dynamique [0.,range[.
    84        Si "vmax" different de "range", c'est la borne superieure
    85        qui peut etre atteinte
    86        (si elle est depassee, on soustrait "range").
    87 \verbatim
    88   r>0 vmax>0
    89   r=24. vmax=24.   -> mets dans [  0,+24[ borne sup exclue
    90   r=24. vmax=12.   -> mets dans ]-12,+12] borne inf exclue
    91 \endverbatim
    92 */
    93 inline void InRange(double *val,double range)
    94              {*val-=range*floor(*val/range);}
    95 inline void InRange(double *val,double range,double vmax)
    96          {InRange(val,range); if(*val>vmax) *val-=range;}
     217
     218// ------------------- Positions des astres -------------------
     219void SunPos(double mjd,double *eclsn,double *ecbsn,double *rsn);
     220void MoonPos(double mjd,double *lmn,double *bmn,double *rho);
     221void PlanetPos(double mjd,int numplan,double *sunecl,double *sunecb,double *sundist
     222               ,double *geodist,double *geoecl,double *geoecb
     223               ,double *diamang,double *mag);
     224/*! \ingroup XAstroPack
     225\brief Same as PlanetPos above with less arguments
     226*/
     227inline void PlanetPos(double mjd,int numplan,double *geoecl,double *geoecb
     228                      ,double *geodist,double *diamang)
     229{
     230 double sunecl,sunecb,sundist,mag;
     231 PlanetPos(mjd,numplan,&sunecl,&sunecb,&sundist,geodist,geoecl,geoecb,diamang,&mag);
     232}
     233
     234/*! \ingroup XAstroPack
     235\brief Give Jupiter position
     236*/
     237inline void JupiterPos(double mjd,double *ecl,double *ecb,double *geodist,double *diamang)
     238{
     239  PlanetPos(mjd,JUPITER,ecl,ecb,geodist,diamang);
     240}
     241
     242/*! \ingroup XAstroPack
     243\brief Give Saturn position
     244*/
     245inline void SaturnPos(double mjd,double *ecl,double *ecb,double *geodist,double *diamang)
     246{
     247  PlanetPos(mjd,SATURN,ecl,ecb,geodist,diamang);
     248}
    97249
    98250#endif
Note: See TracChangeset for help on using the changeset viewer.