Changeset 1678 in Sophya for trunk/SophyaExt/XAstroPack


Ignore:
Timestamp:
Oct 9, 2001, 11:47:57 PM (24 years ago)
Author:
cmv
Message:

EqHToHor et reverse cmv 9/10/01

Location:
trunk/SophyaExt/XAstroPack
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaExt/XAstroPack/xastropack.cc

    r1628 r1678  
    1818// EQUATORIALE: ascension droite en heures [0,24[ (ra)
    1919//              declinaison en degres [-90,90]    (dec)
    20 //              angle horaire en heures [-12,12] (-12=12) (ha)  tsid=ha+ra
     20//              angle horaire en heures [-12,12] (-12=12) (ha)
     21//              temps sideral du lieu: tsid=ha+ra (ou lst)
    2122// GALACTIQUE: longitude en degres [0,360[  (glng)
    2223//             latitude en degres [-90,90]  (glat)
    … …  
    5152/*! \ingroup XAstroPack
    5253\brief Compute MJD from date
     54\verbatim
     55 MJD =  modified Julian date (number of days elapsed since 1900 jan 0.5),
     56\endverbatim
    5357*/
    5458double MJDfrDate(double dy,int mn,int yr)
    … …  
    97101
    98102/*! \ingroup XAstroPack
     103\brief Given a year,
     104*/
     105int 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*/
     114int 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*/
     122int 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*/
     132double 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*/
     140double HfrMJD(double mjd)
     141{
     142  return mjd_hr(mjd);
     143}
     144
     145/*! \ingroup XAstroPack
    99146\brief Give GST from UTC
    100147\verbatim
    … …  
    131178*/                                                                             
    132179double GST0(double mjd0)
    133 /* Copie depuis le code de Xephem car pas prototype */
     180/* Copie depuis le code de Xephem (utc_gst.c) car pas prototype*/
    134181{
    135182 double T, x;
    … …  
    143190
    144191/*! \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*/                                                                             
     196double LSTfrGST(double gst,double geolng)
     197{
     198  double lst = gst + geolng *12./180.;
     199  InRange(&lst,24.);
     200  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*/                                                                             
     207double 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
    145219\brief Compute precession between 2 dates.
    146220*/                                                                             
    … …  
    150224 dec1 *= PI/180.;  // radians
    151225 precess(mjd1,mjd2,&ra1,&dec1);
    152  *ra2 = ra1*12./PI; if(*ra2>24.) *ra2 -= 24.; if(*ra2==24.) *ra2 = 0.;
     226 *ra2 = ra1*12./PI; InRange(ra2,24.);
    153227 *dec2 = dec1*180./PI;
    154228}
    … …  
    166240
    167241/*! \ingroup XAstroPack
    168 \brief Give the hour angle from sideral time and right ascencion
    169 */                                                                             
    170 double HafrRaTS(double gst,double ra)
    171 {
    172   double ha = gst - ra;
     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*/                                                                             
     246double HafrRaTS(double lst,double ra)
     247{
     248  double ha = lst - ra;
    173249  // Attention au probleme de la discontinuite 0h <==> 24h
    174250  // ts=1  ra=23 ; (ts-ra)=-22 <-12 --> ha = +2 = +24 + (ts-ra)
    175251  // ts=23 ra=1  ; (ts-ra)=+22 >+12 --> ha = -2 = -24 + (ts-ra)
    176   if(ha==-12.) ha = 12.; if(ha<-12.) ha += 24.; if(ha>12.) ha -= 24.;
     252  InRange(&ha,24.,12.);
    177253  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*/                                                                             
     261double 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),
     270find 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*/
     274double 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;
    178281}
    179282
    … …  
    247350
    248351/*! \ingroup XAstroPack
    249 \brief Convert equatorial coordinates into galactic coordinates
     352\brief Convert equatorial coordinates for the given epoch into galactic coordinates
    250353*/                                                                             
    251354void EqtoGal(double mjd,double ra,double dec, double *glng,double *glat)
    … …  
    256359 eq_gal(mjd,ra,dec,glat,glng);
    257360 // Vraiment bizarre, sur Linux-g++ glng>=360 ne comprend pas glng==360 ! (CMV)
    258  *glng *= 180./PI; if(*glng>360.) *glng -= 360.; if(*glng==360.) *glng = 0.;
     361 *glng *= 180./PI; InRange(glng,360.);
    259362 *glat *= 180./PI;
    260363}
    261364
    262365/*! \ingroup XAstroPack
    263 \brief Convert galactic coordinates into equatorial coordinates
     366\brief Convert galactic coordinates into equatorial coordinates at the given epoch
    264367*/                                                                             
    265368void GaltoEq(double mjd,double glng,double glat,double *ra,double *dec)
    … …  
    269372 glat *= PI/180.;  // radians
    270373 gal_eq (mjd,glat,glng,ra,dec);
    271  *ra  *= 12./PI; if(*ra>24.) *ra -= 24.; if(*ra==24.) *ra = 0.;
     374 *ra  *= 12./PI; InRange(ra,24.);
    272375 *dec *= 180./PI;
    273376}
    274377
    275378/*! \ingroup XAstroPack
    276 \brief Convert equatorial coordinates into horizontal coordinates
    277 */                                                                             
    278 void EqtoHor(double geolat,double ha,double dec,double *az,double *alt)
     379\brief Convert equatorial coordinates (with hour angle instead of right ascension) into horizontal coordinates.
     380*/                                                                             
     381void EqHtoHor(double geolat,double ha,double dec,double *az,double *alt)
    279382// Coordonnees equatoriales -> Coordonnees horizontales
    280383{
    … …  
    284387 hadec_aa (geolat,ha,dec,alt,az);
    285388 *alt *= 180./PI;
    286  *az  *= 180./PI; if(*az>360.) *az -= 360.; if(*az==360.) *az = 0.;
    287 }
    288 
    289 /*! \ingroup XAstroPack
    290  Convert horizontal coordinates into equatorial coordinates
    291 */                                                                             
    292 void HortoEq(double geolat,double az,double alt,double *ha,double *dec)
     389 *az  *= 180./PI; InRange(az,360.);
     390}
     391
     392/*! \ingroup XAstroPack
     393 Convert horizontal coordinates into equatorial coordinates (with hour angle instead of right ascension).
     394*/                                                                             
     395void HortoEqH(double geolat,double az,double alt,double *ha,double *dec)
    293396// Coordonnees horizontales -> Coordonnees equatoriales
    294397{
    … …  
    297400 az  *= PI/180.;  // radians
    298401 aa_hadec (geolat,alt,az,ha,dec);
    299  *ha  *= 12./PI;
    300   if(*ha==-12.) *ha = 12.; if(*ha<-12.) *ha += 24.; if(*ha>12.) *ha -= 24.;
     402 *ha  *= 12./PI; InRange(ha,24.,12.);
    301403 *dec *= 180./PI;
    302404}
    303405
    304406/*! \ingroup XAstroPack
    305 \brief Convert equatorial coordinates into ecliptic coordinates
     407\brief Convert equatorial coordinates into horizontal coordinates.
     408*/                                                                             
     409void EqtoHor(double geolat,double lst,double ra,double dec,double *az,double *alt)
     410// Coordonnees equatoriales -> Coordonnees horizontales
     411{
     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
     417 hadec_aa (geolat,ha,dec,alt,az);
     418 *alt *= 180./PI;
     419 *az  *= 180./PI; InRange(az,360.);
     420}
     421
     422/*! \ingroup XAstroPack
     423 Convert horizontal coordinates into equatorial coordinates.
     424*/                                                                             
     425void HortoEq(double geolat,double lst,double az,double alt,double *ra,double *dec)
     426// Coordonnees horizontales -> Coordonnees equatoriales
     427{
     428 double ha;
     429 geolat *= PI/180.;
     430 alt *= PI/180.;  // radians
     431 az  *= PI/180.;  // radians
     432 aa_hadec (geolat,alt,az,&ha,dec);
     433 ha *= 12./PI;
     434 *ra = lst - ha; InRange(ra,24.);
     435 *dec *= 180./PI;
     436}
     437
     438/*! \ingroup XAstroPack
     439\brief Convert equatorial coordinates into geocentric ecliptic coordinates given the modified Julian date.
     440\warning Correction for the effect on the angle of the obliquity due to nutation is not included.
    306441*/                                                                             
    307442// Attention, j'ai modifie eq_ecl.c pour proteger NaN
    … …  
    315450 dec *= PI/180.;  // radians
    316451 eq_ecl(mjd,ra,dec,eclat,eclng);
    317  *eclng *= 180./PI; if(*eclng>360.) *eclng -= 360.; if(*eclng==360.) *eclng = 0.;
     452 *eclng *= 180./PI; InRange(eclng,360.);
    318453 *eclat *= 180./PI;
    319454}
    320455
    321456/*! \ingroup XAstroPack
    322 \brief Convert ecliptic coordinates into equatorial coordinates
     457\brief Convert geocentric ecliptic coordinates into equatorial coordinates given the modified Julian date.
     458\warning Correction for the effect on the angle of the obliquity due to nutation is not included.
    323459*/                                                                             
    324460void EcltoEq(double mjd,double eclng,double eclat,double *ra,double *dec)
    … …  
    328464 eclng *= PI/180.;  // radians
    329465 ecl_eq(mjd,eclat,eclng,ra,dec);
    330  *ra  *= 12./PI; if(*ra>24.) *ra -= 24.; if(*ra==24.) *ra = 0.;
     466 *ra  *= 12./PI; InRange(ra,24.);
    331467 *dec *= 180./PI;
    332468}
    … …  
    344480void SunPos(double mjd,double *eclsn,double *ecbsn)
    345481{
    346   double rsn;
    347   sunpos(mjd,eclsn,&rsn,ecbsn);
    348  *eclsn *= 180./PI; if(*eclsn>360.) *eclsn -= 360.; if(*eclsn==360.) *eclsn = 0.;
     482 double rsn;
     483 sunpos(mjd,eclsn,&rsn,ecbsn);
     484 *eclsn *= 180./PI; InRange(eclsn,360.);
    349485 *ecbsn *= 180./PI;
    350486}
    … …  
    363499  double rho,msp,mdp;
    364500  moon(mjd,eclmn,ecbmn,&rho,&msp,&mdp);
    365  *eclmn *= 180./PI; if(*eclmn>360.) *eclmn -= 360.; if(*eclmn==360.) *eclmn = 0.;
     501 *eclmn *= 180./PI; InRange(eclmn,360.);
    366502 *ecbmn *= 180./PI;
    367503}
    … …  
    388524void PlanetPos(double mjd,int numplan,double *ecl,double *ecb,double *diamang)
    389525{
    390   double lpd0,psi0,rp0,rho0,mag;
    391   plans(mjd,numplan,&lpd0,&psi0,&rp0,&rho0,ecl,ecb,diamang,&mag);
    392  *ecl *= 180./PI; if(*ecl>360.) *ecl -= 360.; if(*ecl==360.) *ecl = 0.;
     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.);
    393529 *ecb *= 180./PI;
    394530}

  • trunk/SophyaExt/XAstroPack/xastropack.h

    r1515 r1678  
    4646double MJDfrYear(double yr);
    4747void YDfrMJD(double mjd,double *dy,int *yr);
     48int IsLeapYear(int y);
     49int DayOrder(double mjd,int *dow);
     50int DaysInMonth(double mjd);
     51double MJDat0hFrMJD(double mjd);
     52double HfrMJD(double mjd);
    4853double GSTfrUTC(double mjd0,double utc);
    4954double UTCfrGST(double mjd0,double gst);
    5055double GST0(double mjd0);
     56double LSTfrGST(double gst,double geolng);
     57double LSTfrMJD(double mjd,double geolng);
    5158void Precess(double mjd1,double mjd2,double ra1,double dec1,double *ra2,double *dec2);
    5259double AirmassfrAlt(double alt);
    53 double HafrRaTS(double gst,double ra);
     60double HafrRaTS(double lst,double ra);
     61double RafrHaTS(double lst,double ha);
     62double HelioCorr(double jd,double ra,double dec);
    5463void HMSfrHdec(double hd,int *h,int *mn,double *s);
    5564double HdecfrHMS(int h,int mn,double s);
    … …  
    5867void EqtoGal(double mjd,double ra,double dec,double *glng,double *glat);
    5968void GaltoEq(double mjd,double glng,double glat,double *ra,double *dec);
    60 void EqtoHor(double geolat,double ha,double dec,double *az,double *alt);
    61 void HortoEq(double geolat,double az,double alt,double *ha,double *dec);
     69void EqHtoHor(double geolat,double ha,double dec,double *az,double *alt);
     70void HortoEqH(double geolat,double az,double alt,double *ha,double *dec);
     71void EqtoHor(double geolat,double lst,double ra,double dec,double *az,double *alt);
     72void HortoEq(double geolat,double lst,double az,double alt,double *ra,double *dec);
    6273void EqtoEcl(double mjd,double ra,double dec,double *eclng,double *eclat);
    6374void EcltoEq(double mjd,double eclng,double eclat,double *ra,double *dec);
    … …  
    6980int  CoordConvertToStd(TypAstroCoord typ,double& coord1,double& coord2);
    7081
     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*/
     93inline void InRange(double *val,double range)
     94             {*val-=range*floor(*val/range);}
     95inline void InRange(double *val,double range,double vmax)
     96         {InRange(val,range); if(*val>vmax) *val-=range;}
     97
    7198#endif
Note: See TracChangeset for help on using the changeset viewer.