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

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

File:
1 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.;
Note: See TracChangeset for help on using the changeset viewer.