#include #include #include "sopnamsp.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 // >>>> Corrections effectuees dans la version Xephem 3.5 // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> /*! \defgroup XAstroPack XAstroPack module This module contains simple programs to perform various astronomical computation (based on the libastro of Xephem). \verbatim // TEMPS: modified Julian date (mjd) (number of days elapsed since 1900 jan 0.5) // jour [1,31] (dy) // mois [1,12] (mn) // annee (yr) // universal time [0,24[ (utc) // Greenwich mean siderial [0,24[ (gst) // Greenwich mean siderial at 0h UT [0,24[ (gst0) // EQUATORIALE: ascension droite en heures [0,24[ (ra) // declinaison en degres [-90,90] (dec) // 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) // latitude en degres [-90,90] (glat) // (colatitude en degres [0,180] (gcolat)) // HORIZONTAL: azimuth en degres [0,360[ (az) // (angle round to the east from north+) // altitude ou elevation en degres [-90,90] (alt) // (distance zenitale en degres [0,180] (zendist)) // ECLIPTIQUE: lontitude ecliptique en degres [0,360[ (eclng) // (angle round counter clockwise from the vernal equinoxe) // latitude ecliptique en degres [-90,90] (eclat) // (colatitude ecliptique en degres [0,180] (eccolat)) // GEOGRAPHIE: longitude en degres ]-180,180] (geolng) // (angle <0 vers l'ouest, >0 vers l'est) // latitude en degres [-90,90] (north>0 sud<0) (geolat) // (colatitude en degres [0,180] (north=0, sud=180) (geocolat)) // // --- Remarque sur la colatitude: // La latitude peut etre remplacee par la colatitude // (ou altitude/elevation par la distance zenitale): // latitude : [-90,90] avec 0=equateur, 90=pole nord, -90=pole sud // colatitude : [0,180] avec 0=pole nord, 90=equateur, 180=pole sud // colatitude = 90.-latitude , latitude = 90.-colatitude \endverbatim */ /*! \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 = Problem TypAstroCoord du systeme de coordonnees retourne - Les types de coordonnees (A,B) sont definies dans le enum TypAstroCoord (unsigned long): La premiere coordonnee "A" est de type "longitude" (alpha,longitude,azimuth) La deuxieme coordonnee "B" est de type "latidude" (delta,latitude,altitude ou elevation) *** Definitions des unites des coordonnees - Coordonnee: TypCoordH en heure TypCoordD en degre TypCoordR en radian - Coordonnee "A": TypCoord1H en heure TypCoord1D en degre TypCoord1R en radian - Defaut (pas de bit leve): radians - Coordonnee "B": TypCoord2H en heure TypCoord2D en degre TypCoord2R en radian - Defaut (pas de bit leve): radians *** Definitions des types d'intervalle utilises - Coordonnee "A": TypCoord1C type intervalle [0,24[ ou [0,360[ ou [0,2Pi[ TypCoord1L type intervalle centre ]-12,12] ou ]-180,180] ou ]-Pi,Pi] - Defaut (pas de bit leve): TypCoord1C - Coordonnee "B": TypCoord2C type intervalle (colatitude) [0,12] ou [0,180] ou [0,Pi] TypCoord2L type intervalle centre (latitude) [-6,6] ou [-90,90] ou [-Pi/2,Pi/2] - Defaut (pas de bit leve): TypCoord2L (latitude) *** Les systemes de coordonnes astronomiques TypCoordEq coordonnees Equatoriales alpha,delta TypCoordGal coordonnees Galactiques gLong, gLat TypCoordHor coordonnees Horizontales azimuth,altitude TypCoordEcl coordonnees Ecliptiques EclLong,EclLat (Pas de defaut) *** Les systemes de coordonnes astronomiques "standard" TypCoordEqStd alpha en heure=[0,24[ delta en degre=[-90,90] TypCoordGalStd long en degre=[0,360[ lat en degre=[-90,90] (latitude) TypCoordHorStd long en degre=[0,360[ lat en degre=[-90,90] (latitude) TypCoordEclStd long en degre=[0,360[ lat en degre=[-90,90] (latitude) \endverbatim */ unsigned long CoordConvertToStd(unsigned long typ,double* coord1,double* coord2) { unsigned long rc = TypCoordUndef; if(typ&TypCoordEq) { // ---- Equatoriales alpha,delta if (typ&TypCoord1D) *coord1 = deghr(*coord1); else if(typ&TypCoord1R) *coord1 = radhr(*coord1); else if(!(typ&TypCoord1H)) *coord1 = radhr(*coord1); InRange(coord1,24.); if (typ&TypCoord2H) *coord2 = hrdeg(*coord2); else if(typ&TypCoord2R) *coord2 = raddeg(*coord2); else if(!(typ&TypCoord2D)) *coord2 = raddeg(*coord2); if(typ&TypCoord2C) { InRangeCoLat(coord2,TypUniteD); ToCoLat(coord2,TypUniteD); } else InRangeLat(coord2,TypUniteD); rc=TypCoordEqStd; } else if(typ&TypCoordGal || typ&TypCoordHor || typ&TypCoordEcl) { // ---- Galactiques gLong, gLat // ---- Horizontales azimuth,altitude ou elevation // ---- Ecliptiques EclLong,EclLat if (typ&TypCoord1H) *coord1 = hrdeg(*coord1); else if(typ&TypCoord1R) *coord1 = raddeg(*coord1); else if(!(typ&TypCoord1D)) *coord1 = raddeg(*coord1); InRange(coord1,360.); if (typ&TypCoord2H) *coord2 = hrdeg(*coord2); else if(typ&TypCoord2R) *coord2 = raddeg(*coord2); else if(!(typ&TypCoord2D)) *coord2 = raddeg(*coord2); if(typ&TypCoord2C) { InRangeCoLat(coord2,TypUniteD); ToCoLat(coord2,TypUniteD); } else InRangeLat(coord2,TypUniteD); if (typ&TypCoordGal) rc=TypCoordGalStd; else if(typ&TypCoordHor) rc=TypCoordHorStd; else if(typ&TypCoordEcl) rc=TypCoordEclStd; } else { // Systeme de coordonnees non-connu rc=TypCoordUndef; } return rc; } /*! \ingroup XAstroPack \brief Retourne te type d'unite pour la coordonnee "coordnum" pour un TypAstroCoord valant "typ" \verbatim coordnum : numero de coordonnee: 1 ou 2 retourne: TypUniteH si la coordonnee est en heure TypUniteD si la coordonnee est en degre TypUniteR si la coordonnee est en radian TypUniteR par defaut TypCoordUndef si le numero de coordonnee est errone. \endverbatim */ unsigned long GetCoordUnit(int coordnum,unsigned long typ) { if(coordnum==1) { if (typ&TypCoord1H) return TypUniteH; else if(typ&TypCoord1D) return TypUniteD; else if(typ&TypCoord1R) return TypUniteR; else return TypUniteR; } if(coordnum==2) { if (typ&TypCoord2H) return TypUniteH; else if(typ&TypCoord2D) return TypUniteD; else if(typ&TypCoord2R) return TypUniteR; else return TypUniteR; } return TypCoordUndef; } /*! \ingroup XAstroPack \brief Pour decoder et transcrire en TypAstroCoord une chaine donnant la structure du systeme de coordonnees. \verbatim ctype = "CAaBb" C: type de coordonnees: E Equatoriales G Galactiques H Horizontales S Ecliptiques pas de defaut A: unite de la coordonnee 1 (alpha,longitude etc...) H heure D degre R radian defaut radian a: type d'intervalle pour la coordonnee 1 C intervalle [0,24[ [0,360[ [0,2*Pi[ L intervalle [-12,12[ [-180,180[ [-Pi,Pi[ (defaut: C) A: unite de la coordonnee 2 (delta,latitude etc...) H heure D degre R radian defaut radian a: type d'intervalle pour la coordonnee 2 C intervalle [0,12] [0,180] [0,Pi] (type colatitude) L intervalle [-6,6] [-90,90][ [-Pi/2,Pi/2] (defaut: L) Exemple: GDCDL : galactiques long=[0,360[ lat=[-90,90] GDxDx ou GDxD: idem Gxxxx ou G : galactiques long=[0,2*Pi[ lat=[-Pi/2,Pi/2] Exemple: EHCDL : equatoriales alpha=[0,24[ delta=[-90,90] EHxDx ou EHxD : idem Exxxx ou E : equatoriales alpha=[0,2*Pi[ delta=[-Pi/2,Pi/2] - Retourne 0 si probleme dans la chaine \endverbatim */ unsigned long DecodeTypAstro(const char *ctype) { if(ctype==NULL) return TypCoordUndef; int len = strlen(ctype); if(len<1) return TypCoordUndef; unsigned long typ=0; // Le type de systeme de coordonnees int i = 0; if (ctype[i]=='e' || ctype[i]=='E') typ=TypCoordEq; else if(ctype[i]=='g' || ctype[i]=='G') typ=TypCoordGal; else if(ctype[i]=='h' || ctype[i]=='H') typ=TypCoordHor; else if(ctype[i]=='s' || ctype[i]=='S') typ=TypCoordEcl; else return TypCoordUndef; // La coordonnee 1 i = 1; if(i>=len) {typ |= TypCoord1R|TypCoord1C|TypCoord2R|TypCoord2L; return typ;} if (ctype[i]=='h' || ctype[i]=='H') typ |= TypCoord1H; else if(ctype[i]=='d' || ctype[i]=='D') typ |= TypCoord1D; else if(ctype[i]=='r' || ctype[i]=='R') typ |= TypCoord1R; else typ |= TypCoord1R; i = 2; if(i>=len) {typ |= TypCoord1C|TypCoord2R|TypCoord2L; return typ;} if (ctype[i]=='c' || ctype[i]=='C') typ |= TypCoord1C; else if(ctype[i]=='l' || ctype[i]=='L') typ |= TypCoord1L; else typ |= TypCoord1C; // La coordonnee 2 i = 3; if(i>=len) {typ |= TypCoord2R|TypCoord2L; return typ;} if (ctype[i]=='h' || ctype[i]=='H') typ |= TypCoord2H; else if(ctype[i]=='d' || ctype[i]=='D') typ |= TypCoord2D; else if(ctype[i]=='r' || ctype[i]=='R') typ |= TypCoord2R; else typ |= TypCoord2R; i = 4; if(i>=len) {typ |= TypCoord2L; return typ;} if (ctype[i]=='c' || ctype[i]=='C') typ |= TypCoord2C; else if(ctype[i]=='l' || ctype[i]=='L') typ |= TypCoord2L; else typ |= TypCoord2L; // Return return typ; } /*! \ingroup XAstroPack \brief Idem DecodeTypAstro(char *) mais a l'envers */ string DecodeTypAstro(unsigned long typ) { string s = ""; if (typ&TypCoordEq) s += "E"; else if(typ&TypCoordGal) s += "G"; else if(typ&TypCoordHor) s += "H"; else if(typ&TypCoordEcl) s += "S"; else s += "x"; if (typ&TypCoord1H) s += "H"; else if(typ&TypCoord1D) s += "D"; else if(typ&TypCoord1R) s += "R"; else s += "x"; if (typ&TypCoord1C) s += "C"; else if(typ&TypCoord1L) s += "L"; else s += "x"; if (typ&TypCoord2H) s += "H"; else if(typ&TypCoord2D) s += "D"; else if(typ&TypCoord2R) s += "R"; else s += "x"; if (typ&TypCoord2C) s += "C"; else if(typ&TypCoord2L) s += "L"; else s += "x"; return s; } /*! \ingroup XAstroPack \brief Pour convertir la latitude en colatitude et vice-versa (in place) \verbatim val = valeur a convertir qui doit etre: si type "latitude" dans [-6,6] [-90,90] [-Pi/2,Pi/2] si type "colatitude" dans [0,12] [0,180] [0,Pi] typ = type d'unite: heure TypUniteH degre TypUniteD radian TypUniteR (Defaut: radian TypUniteR) \endverbatim */ void ToCoLat(double *val,unsigned long typ) { if (typ&TypUniteH) *val = 6. - *val; else if(typ&TypUniteD) *val = 90. - *val; else if(typ&TypUniteR) *val = PI/2. - *val; else *val = PI/2. - *val; } /*! \ingroup XAstroPack \brief Pour remettre la valeur de la COLATITUDE "val" dans la dynamique [0.,range] \verbatim val = valeur a convertir qui doit etre mise dans [0,12] [0,180] [0,Pi] typ = type d'unite: heure TypUniteH degre TypUniteD radian TypUniteR ex en degre: 0 -> 0 , 90 -> 90 , 180 -> 180 , 270 -> 90 , 360 -> 0 -90 -> 90 , -180 -> 180 , -270 -> 90 , -360 -> 0 \endverbatim */ void InRangeCoLat(double *val,unsigned long typ) { double range=PI; if(typ==TypUniteH) range=12.; else if(typ==TypUniteD) range=180.; InRange(val,2.*range,range); if(*val<0.) *val*=-1.; } /*! \ingroup XAstroPack \brief Pour remettre la valeur de la LATITUDE "val" dans la dynamique [-range/2,range/2] \verbatim val = valeur a convertir qui doit etre mise dans [-6,6] [-90,90] [-Pi/2,Pi/2] typ = type d'unite: heure TypUniteH degre TypUniteD radian TypUniteR ex en degre: 0 -> 0 , 90 -> 90 , 180 -> 0 , 270 -> -90 , 360 -> 0 -90 -> -90 , -180 -> 0 , -270 -> 90 , -360 -> 0 \endverbatim */ void InRangeLat(double *val,unsigned long typ) { double range = PI; if(typ==TypUniteH) range = 12.; else if(typ==TypUniteD) range = 180.; InRange(val,2.*range,range); if(*val>range/2.) *val = range - *val; else if(*val<-range/2.) *val = -(range + *val); } /*! \ingroup XAstroPack \brief Compute MJD from date \verbatim MJD = modified Julian date (number of days elapsed since 1900 jan 0.5), dy is the decimale value of the day: dy = int(dy) + utc/24. \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 */ double GST0(double mjd0) /* Copie depuis le code de Xephem (utc_gst.c) car pas prototype*/ { double T, x; T = ((int)(mjd0 - 0.5) + 0.5 - J2000)/36525.0; x = 24110.54841 + (8640184.812866 + (0.093104 - 6.2e-6 * T) * T) * T; x /= 3600.0; range(&x, 24.0); return (x); } /*! \ingroup XAstroPack \brief return local sidereal time from modified julian day and longitude \warning nutation or obliquity correction are taken into account. */ double LSTfrMJD(double mjd,double geolng) { double eps,lst,deps,dpsi; utc_gst(mjd_day(mjd),mjd_hr(mjd),&lst); lst += deghr(geolng); obliquity(mjd,&eps); nutation(mjd,&deps,&dpsi); lst += radhr(dpsi*cos(eps+deps)); InRange(&lst,24.); return lst; } /*! \ingroup XAstroPack \brief Give a time in h:mn:s from a decimal hour \verbatim // INPUT: hd // OUTPUT: h mn s (h,mn,s >=< 0) // REMARQUE: si hd<0 alors h<0 ET mn<0 ET s<0 // EX: 12.51 -> h=12 mn=30 s=10 ; // -12.51 -> h=-12 mn=-30 s=-10 ; \endverbatim */ void HMSfrHdec(double hd,int *h,int *mn,double *s) { int sgn=1; if(hd<0.) {sgn=-1; hd*=-1.;} *h = int(hd); *mn = int((hd-(double)(*h))*60.); *s = (hd - (double)(*h) - (double)(*mn)/60.)*3600.; // pb precision if(*s<0.) *s = 0.; if(*s>60. || *s==60.) {*s-=60.; *mn+=1;} // s=double attention comparaison if(*mn<0) *mn = 0; if(*mn>=60) {*mn-=60; *h+=1;} *h *= sgn; *mn *= sgn; *s *= (double)sgn; } /*! \ingroup XAstroPack \brief Give a decimal hour from a time in h:mn:s \verbatim // INPUT: h , mn , s (h,mn,s >=< 0) // RETURN: en heures decimales // REMARQUE: pour avoir hd=-12.51 <- h=-12 mn=-30 s=-10 \endverbatim */ double HdecfrHMS(int h,int mn,double s) { return ((double)h + (double)mn/60. + s/3600.); } /*! \ingroup XAstroPack \brief Give a time string from a time in h:mn:s \verbatim // INPUT: h , mn , s (h,mn,s >=< 0) // RETURN: string h:mn:s \endverbatim */ string ToStringHMS(int h,int mn,double s) { double hd = HdecfrHMS(h,mn,s); // put in range HMSfrHdec(hd,&h,&mn,&s); char str[128]; if(hd<0.) sprintf(str,"-%d:%d:%.3f",-h,-mn,-s); else sprintf(str,"%d:%d:%.3f",h,mn,s); string dum = str; return dum; } /*! \ingroup XAstroPack \brief Give a time string from a decimal hour */ string ToStringHdec(double hd) { int h,mn; double s; HMSfrHdec(hd,&h,&mn,&s); return ToStringHMS(h,mn,s); } /*! \ingroup XAstroPack \brief Compute precession between 2 dates. */ void Precess(double mjd1,double mjd2,double ra1,double dec1,double *ra2,double *dec2) { ra1 = hrrad(ra1); // radians dec1 = degrad(dec1); // radians precess(mjd1,mjd2,&ra1,&dec1); *ra2 = radhr(ra1); InRange(ra2,24.); *dec2 = raddeg(dec1); } /*! \ingroup XAstroPack \brief Convert equatorial coordinates for the given epoch into galactic coordinates */ void EqtoGal(double mjd,double ra,double dec, double *glng,double *glat) // Coordonnees equatoriales -> Coordonnees galactiques { ra = hrrad(ra); // radians dec = degrad(dec); // radians eq_gal(mjd,ra,dec,glat,glng); // Vraiment bizarre, sur Linux-g++ glng>=360 ne comprend pas glng==360 ! (CMV) *glng = raddeg(*glng); InRange(glng,360.); *glat = raddeg(*glat); } /*! \ingroup XAstroPack \brief Convert galactic coordinates into equatorial coordinates at the given epoch */ void GaltoEq(double mjd,double glng,double glat,double *ra,double *dec) // Coordonnees galactiques -> Coordonnees equatoriales { glng = degrad(glng); // radians glat = degrad(glat); // radians gal_eq (mjd,glat,glng,ra,dec); *ra = radhr(*ra); InRange(ra,24.); *dec = raddeg(*dec); } /*! \ingroup XAstroPack \brief Convert equatorial coordinates (with hour angle instead of right ascension) into horizontal coordinates. */ void EqHtoHor(double geolat,double ha,double dec,double *az,double *alt) // Coordonnees equatoriales -> Coordonnees horizontales { geolat = degrad(geolat); // radians ha = hrrad(ha); // radians dec = degrad(dec); // radians hadec_aa (geolat,ha,dec,alt,az); *alt = raddeg(*alt); *az = raddeg(*az); InRange(az,360.); } /*! \ingroup XAstroPack Convert horizontal coordinates into equatorial coordinates (with hour angle instead of right ascension). */ void HortoEqH(double geolat,double az,double alt,double *ha,double *dec) // Coordonnees horizontales -> Coordonnees equatoriales { geolat = degrad(geolat); // radians alt = degrad(alt); // radians az = degrad(az); // radians aa_hadec (geolat,alt,az,ha,dec); *ha = radhr(*ha); InRange(ha,24.,12.); *dec = raddeg(*dec); } /*! \ingroup XAstroPack \brief Convert equatorial coordinates into horizontal coordinates. */ void EqtoHor(double geolat,double lst,double ra,double dec,double *az,double *alt) // Coordonnees equatoriales -> Coordonnees horizontales { double ha = lst - ra; InRange(&ha,24.,12.); geolat = degrad(geolat); // radians ha = hrrad(ha); // radians dec = degrad(dec); // radians hadec_aa (geolat,ha,dec,alt,az); *alt = raddeg(*alt); *az = raddeg(*az); InRange(az,360.); } /*! \ingroup XAstroPack Convert horizontal coordinates into equatorial coordinates. */ void HortoEq(double geolat,double lst,double az,double alt,double *ra,double *dec) // Coordonnees horizontales -> Coordonnees equatoriales { double ha; geolat = degrad(geolat); // radians alt = degrad(alt); // radians az = degrad(az); // radians aa_hadec (geolat,alt,az,&ha,dec); *ra = lst - radhr(ha); InRange(ra,24.); *dec = raddeg(*dec); } /*! \ingroup XAstroPack \brief Convert equatorial coordinates into geocentric ecliptic coordinates given the modified Julian date. \warning Correction for the effect on the angle of the obliquity due to nutation is not included. */ void EqtoEcl(double mjd,double ra,double dec,double *eclng,double *eclat) // Coordonnees equatoriales -> Coordonnees ecliptiques { ra = hrrad(ra); // radians dec = degrad(dec); // radians eq_ecl(mjd,ra,dec,eclat,eclng); *eclng = raddeg(*eclng); InRange(eclng,360.); *eclat = raddeg(*eclat); } /*! \ingroup XAstroPack \brief Convert geocentric ecliptic coordinates into equatorial coordinates given the modified Julian date. \warning Correction for the effect on the angle of the obliquity due to nutation is not included. */ void EcltoEq(double mjd,double eclng,double eclat,double *ra,double *dec) // Coordonnees ecliptiques -> Coordonnees equatoriales { eclat = degrad(eclat); // radians eclng = degrad(eclng); // radians ecl_eq(mjd,eclat,eclng,ra,dec); *ra = radhr(*ra); InRange(ra,24.); *dec = raddeg(*dec); } /*! \ingroup XAstroPack \brief Give Sun position \verbatim given the modified JD, mjd, return the true geocentric ecliptic longitude of the sun for the mean equinox of the date, *eclsn, in degres, the sun-earth distance, *rsn, in AU, and the latitude *ecbsn, in degres (since this is always <= 1.2 arcseconds, in can be neglected by calling with ecbsn = NULL). - REMARQUE: * if the APPARENT ecliptic longitude is required, correct the longitude for * nutation to the true equinox of date and for aberration (light travel time, * approximately -9.27e7/186000/(3600*24*365)*2*pi = -9.93e-5 radians). \endverbatim */ void SunPos(double mjd,double *eclsn,double *ecbsn,double *rsn) { sunpos(mjd,eclsn,rsn,ecbsn); *eclsn = raddeg(*eclsn); InRange(eclsn,360.); if(ecbsn!=NULL) *ecbsn = raddeg(*ecbsn); } /*! \ingroup XAstroPack \brief Give Moon position \verbatim given the mjd, find the geocentric ecliptic longitude, lam, and latitude, bet, and geocentric distance, rho in a.u. for the moon. also return the sun's mean anomaly, *msp, and the moon's mean anomaly, *mdp. (for the mean equinox) \endverbatim */ void MoonPos(double mjd,double *eclmn,double *ecbmn,double *rho) { double msp,mdp; moon(mjd,eclmn,ecbmn,rho,&msp,&mdp); *eclmn = raddeg(*eclmn); InRange(eclmn,360.); *ecbmn = raddeg(*ecbmn); } /*! \ingroup XAstroPack \brief Give planet position \verbatim * given a modified Julian date, mjd, and a planet, p, find: * sunecl: heliocentric longitude, * sunecb: heliocentric latitude, * sundist: distance from the sun to the planet, * geodist: distance from the Earth to the planet, * none corrected for light time, ie, they are the true values for the * given instant. * geoecl: geocentric ecliptic longitude, * geoecb: geocentric ecliptic latitude, * each corrected for light time, ie, they are the apparent values as * seen from the center of the Earth for the given instant. * diamang: angular diameter in arcsec at 1 AU, * mag: visual magnitude when 1 AU from sun and earth at 0 phase angle. * (for the mean equinox) * all angles are in degres, all distances in AU. * * corrections for nutation and abberation must be made by the caller. The RA * and DEC calculated from the fully-corrected ecliptic coordinates are then * the apparent geocentric coordinates. Further corrections can be made, if * required, for atmospheric refraction and geocentric parallax. \endverbatim */ void PlanetPos(double mjd,PLCode numplan,double *sunecl,double *sunecb,double *sundist ,double *geodist,double *geoecl,double *geoecb ,double *diamang,double *mag) { plans(mjd,numplan,sunecl,sunecb,sundist,geodist,geoecl,geoecb,diamang,mag); *geoecl = raddeg(*geoecl); InRange(geoecl,360.); *geoecb = raddeg(*geoecb); *sunecl = raddeg(*sunecl); InRange(sunecl,360.); *sunecb = raddeg(*sunecb); }