| 1 | #include <math.h>
|
|---|
| 2 | #include <stdio.h>
|
|---|
| 3 | #include "xastropack.h"
|
|---|
| 4 |
|
|---|
| 5 | // BUGS BUGS BUGS BUGS BUGS BUGS BUGS BUGS BUGS BUGS
|
|---|
| 6 | // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 7 | // >>>> Corrections de divers bugs trouve dans libastro (CMV)
|
|---|
| 8 | // 1******* In the file vsop87.c line 154:
|
|---|
| 9 | // p = q/(t_abs[alpha] + alpha * t_abs[alpha-1] * 1e-4 + 1e-35);
|
|---|
| 10 | // - to avoid t_abs[-1] when alpha=0, replaced by :
|
|---|
| 11 | // if(alpha>0) p = t_abs[alpha-1]; else p=0;
|
|---|
| 12 | // p = q/(t_abs[alpha] + alpha * p * 1e-4 + 1e-35);
|
|---|
| 13 | // Mail envoye a ecdowney@ClearSkyInstitute.com
|
|---|
| 14 | // 2******* In the file eq_ecl.c line 69:
|
|---|
| 15 | // *q = asin((sy*ceps)-(cy*seps*sx*sw));
|
|---|
| 16 | // eq_ecl.c Protection NaN dans ecleq_aux, replaced by :
|
|---|
| 17 | // *q = (sy*ceps)-(cy*seps*sx*sw);
|
|---|
| 18 | // if(*q<-1.) *q = -PI/2.; else if(*q>1.) *q = PI/2.; else *q = asin(*q);
|
|---|
| 19 | // Mail envoye a ecdowney@ClearSkyInstitute.com
|
|---|
| 20 | // >>>> Corrections effectuees dans la version Xephem 3.5
|
|---|
| 21 | // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|---|
| 22 |
|
|---|
| 23 |
|
|---|
| 24 | /*!
|
|---|
| 25 | \defgroup XAstroPack XAstroPack module
|
|---|
| 26 | This module contains simple programs to perform various
|
|---|
| 27 | astronomical computation (based on the libastro of Xephem).
|
|---|
| 28 |
|
|---|
| 29 | \verbatim
|
|---|
| 30 | // TEMPS: modified Julian date (mjd) (number of days elapsed since 1900 jan 0.5)
|
|---|
| 31 | // jour [1,31] (dy)
|
|---|
| 32 | // mois [1,12] (mn)
|
|---|
| 33 | // annee (yr)
|
|---|
| 34 | // universal time [0,24[ (utc)
|
|---|
| 35 | // Greenwich mean siderial [0,24[ (gst)
|
|---|
| 36 | // Greenwich mean siderial at 0h UT [0,24[ (gst0)
|
|---|
| 37 | // EQUATORIALE: ascension droite en heures [0,24[ (ra)
|
|---|
| 38 | // declinaison en degres [-90,90] (dec)
|
|---|
| 39 | // angle horaire en heures ]-12,12] (-12=12) (ha)
|
|---|
| 40 | // temps sideral du lieu: tsid=ha+ra (ou lst)
|
|---|
| 41 | // GALACTIQUE: longitude en degres [0,360[ (glng)
|
|---|
| 42 | // latitude en degres [-90,90] (glat)
|
|---|
| 43 | // (colatitude en degres [0,180] (gcolat))
|
|---|
| 44 | // HORIZONTAL: azimuth en degres [0,360[ (az)
|
|---|
| 45 | // (angle round to the east from north+)
|
|---|
| 46 | // altitude ou elevation en degres [-90,90] (alt)
|
|---|
| 47 | // (distance zenitale en degres [0,180] (zendist))
|
|---|
| 48 | // ECLIPTIQUE: lontitude ecliptique en degres [0,360[ (eclng)
|
|---|
| 49 | // (angle round counter clockwise from the vernal equinoxe)
|
|---|
| 50 | // latitude ecliptique en degres [-90,90] (eclat)
|
|---|
| 51 | // (colatitude ecliptique en degres [0,180] (eccolat))
|
|---|
| 52 | // GEOGRAPHIE: longitude en degres ]-180,180] (geolng)
|
|---|
| 53 | // (angle <0 vers l'ouest, >0 vers l'est)
|
|---|
| 54 | // latitude en degres [-90,90] (north>0 sud<0) (geolat)
|
|---|
| 55 | // (colatitude en degres [0,180] (north=0, sud=180) (geocolat))
|
|---|
| 56 | //
|
|---|
| 57 | // --- Remarque sur la colatitude:
|
|---|
| 58 | // La latitude peut etre remplacee par la colatitude
|
|---|
| 59 | // (ou altitude/elevation par la distance zenitale):
|
|---|
| 60 | // latitude : [-90,90] avec 0=equateur, 90=pole nord, -90=pole sud
|
|---|
| 61 | // colatitude : [0,180] avec 0=pole nord, 90=equateur, 180=pole sud
|
|---|
| 62 | // colatitude = 90.-latitude , latitude = 90.-colatitude
|
|---|
| 63 | \endverbatim
|
|---|
| 64 | */
|
|---|
| 65 |
|
|---|
| 66 | /*! \ingroup XAstroPack
|
|---|
| 67 | \brief Given a coordinate type "typ", convert to standard for astropack.
|
|---|
| 68 | \verbatim
|
|---|
| 69 | La routine convertit (in place) les coordonnees "coord1","coord2"
|
|---|
| 70 | definies par le type "typ" dans les unites standard de ce systeme
|
|---|
| 71 | de coordonnees.
|
|---|
| 72 | "typ" code le systeme de coordonnees astronomiques et les unites utilisees
|
|---|
| 73 |
|
|---|
| 74 | - Return : 0 = Problem
|
|---|
| 75 | TypAstroCoord du systeme de coordonnees retourne
|
|---|
| 76 |
|
|---|
| 77 | - Les types de coordonnees (A,B) sont definies
|
|---|
| 78 | dans le enum TypAstroCoord (unsigned long):
|
|---|
| 79 | La premiere coordonnee "A" est de type "longitude"
|
|---|
| 80 | (alpha,longitude,azimuth)
|
|---|
| 81 | La deuxieme coordonnee "B" est de type "latidude"
|
|---|
| 82 | (delta,latitude,altitude ou elevation)
|
|---|
| 83 |
|
|---|
| 84 | *** Definitions des unites des coordonnees
|
|---|
| 85 | - Coordonnee:
|
|---|
| 86 | TypCoordH en heure
|
|---|
| 87 | TypCoordD en degre
|
|---|
| 88 | TypCoordR en radian
|
|---|
| 89 | - Coordonnee "A":
|
|---|
| 90 | TypCoord1H en heure
|
|---|
| 91 | TypCoord1D en degre
|
|---|
| 92 | TypCoord1R en radian
|
|---|
| 93 | - Defaut (pas de bit leve): radians
|
|---|
| 94 | - Coordonnee "B":
|
|---|
| 95 | TypCoord2H en heure
|
|---|
| 96 | TypCoord2D en degre
|
|---|
| 97 | TypCoord2R en radian
|
|---|
| 98 | - Defaut (pas de bit leve): radians
|
|---|
| 99 | *** Definitions des types d'intervalle utilises
|
|---|
| 100 | - Coordonnee "A":
|
|---|
| 101 | TypCoord1C type intervalle
|
|---|
| 102 | [0,24[ ou [0,360[ ou [0,2Pi[
|
|---|
| 103 | TypCoord1L type intervalle centre
|
|---|
| 104 | ]-12,12] ou ]-180,180] ou ]-Pi,Pi]
|
|---|
| 105 | - Defaut (pas de bit leve): TypCoord1C
|
|---|
| 106 | - Coordonnee "B":
|
|---|
| 107 | TypCoord2C type intervalle (colatitude)
|
|---|
| 108 | [0,12] ou [0,180] ou [0,Pi]
|
|---|
| 109 | TypCoord2L type intervalle centre (latitude)
|
|---|
| 110 | [-6,6] ou [-90,90] ou [-Pi/2,Pi/2]
|
|---|
| 111 | - Defaut (pas de bit leve): TypCoord2L (latitude)
|
|---|
| 112 | *** Les systemes de coordonnes astronomiques
|
|---|
| 113 | TypCoordEq coordonnees Equatoriales alpha,delta
|
|---|
| 114 | TypCoordGal coordonnees Galactiques gLong, gLat
|
|---|
| 115 | TypCoordHor coordonnees Horizontales azimuth,altitude
|
|---|
| 116 | TypCoordEcl coordonnees Ecliptiques EclLong,EclLat
|
|---|
| 117 | (Pas de defaut)
|
|---|
| 118 | *** Les systemes de coordonnes astronomiques "standard"
|
|---|
| 119 | TypCoordEqStd alpha en heure=[0,24[ delta en degre=[-90,90]
|
|---|
| 120 | TypCoordGalStd long en degre=[0,360[ lat en degre=[-90,90] (latitude)
|
|---|
| 121 | TypCoordHorStd long en degre=[0,360[ lat en degre=[-90,90] (latitude)
|
|---|
| 122 | TypCoordEclStd long en degre=[0,360[ lat en degre=[-90,90] (latitude)
|
|---|
| 123 | \endverbatim
|
|---|
| 124 | */
|
|---|
| 125 | unsigned long CoordConvertToStd(unsigned long typ,double* coord1,double* coord2)
|
|---|
| 126 | {
|
|---|
| 127 | unsigned long rc = TypCoordUndef;
|
|---|
| 128 |
|
|---|
| 129 | if(typ&TypCoordEq) {
|
|---|
| 130 | // ---- Equatoriales alpha,delta
|
|---|
| 131 | if (typ&TypCoord1D) *coord1 = deghr(*coord1);
|
|---|
| 132 | else if(typ&TypCoord1R) *coord1 = radhr(*coord1);
|
|---|
| 133 | else if(!(typ&TypCoord1H)) *coord1 = radhr(*coord1);
|
|---|
| 134 |
|
|---|
| 135 | InRange(coord1,24.);
|
|---|
| 136 |
|
|---|
| 137 | if (typ&TypCoord2H) *coord2 = hrdeg(*coord2);
|
|---|
| 138 | else if(typ&TypCoord2R) *coord2 = raddeg(*coord2);
|
|---|
| 139 | else if(!(typ&TypCoord2D)) *coord2 = raddeg(*coord2);
|
|---|
| 140 |
|
|---|
| 141 | if(typ&TypCoord2C) {
|
|---|
| 142 | InRangeCoLat(coord2,TypUniteD);
|
|---|
| 143 | ToCoLat(coord2,TypUniteD);
|
|---|
| 144 | } else InRangeLat(coord2,TypUniteD);
|
|---|
| 145 |
|
|---|
| 146 | rc=TypCoordEqStd;
|
|---|
| 147 |
|
|---|
| 148 | } else if(typ&TypCoordGal || typ&TypCoordHor || typ&TypCoordEcl) {
|
|---|
| 149 | // ---- Galactiques gLong, gLat
|
|---|
| 150 | // ---- Horizontales azimuth,altitude ou elevation
|
|---|
| 151 | // ---- Ecliptiques EclLong,EclLat
|
|---|
| 152 | if (typ&TypCoord1H) *coord1 = hrdeg(*coord1);
|
|---|
| 153 | else if(typ&TypCoord1R) *coord1 = raddeg(*coord1);
|
|---|
| 154 | else if(!(typ&TypCoord1D)) *coord1 = raddeg(*coord1);
|
|---|
| 155 |
|
|---|
| 156 | InRange(coord1,360.);
|
|---|
| 157 |
|
|---|
| 158 | if (typ&TypCoord2H) *coord2 = hrdeg(*coord2);
|
|---|
| 159 | else if(typ&TypCoord2R) *coord2 = raddeg(*coord2);
|
|---|
| 160 | else if(!(typ&TypCoord2D)) *coord2 = raddeg(*coord2);
|
|---|
| 161 |
|
|---|
| 162 | if(typ&TypCoord2C) {
|
|---|
| 163 | InRangeCoLat(coord2,TypUniteD);
|
|---|
| 164 | ToCoLat(coord2,TypUniteD);
|
|---|
| 165 | } else InRangeLat(coord2,TypUniteD);
|
|---|
| 166 |
|
|---|
| 167 | if (typ&TypCoordGal) rc=TypCoordGalStd;
|
|---|
| 168 | else if(typ&TypCoordHor) rc=TypCoordHorStd;
|
|---|
| 169 | else if(typ&TypCoordEcl) rc=TypCoordEclStd;
|
|---|
| 170 |
|
|---|
| 171 | } else { // Systeme de coordonnees non-connu
|
|---|
| 172 | rc=TypCoordUndef;
|
|---|
| 173 | }
|
|---|
| 174 |
|
|---|
| 175 | return rc;
|
|---|
| 176 | }
|
|---|
| 177 |
|
|---|
| 178 | /*! \ingroup XAstroPack
|
|---|
| 179 | \brief Retourne te type d'unite pour la coordonnee "coordnum"
|
|---|
| 180 | pour un TypAstroCoord valant "typ"
|
|---|
| 181 | \verbatim
|
|---|
| 182 | coordnum : numero de coordonnee: 1 ou 2
|
|---|
| 183 | retourne: TypUniteH si la coordonnee est en heure
|
|---|
| 184 | TypUniteD si la coordonnee est en degre
|
|---|
| 185 | TypUniteR si la coordonnee est en radian
|
|---|
| 186 | TypUniteR par defaut
|
|---|
| 187 | TypCoordUndef si le numero de coordonnee est errone.
|
|---|
| 188 | \endverbatim
|
|---|
| 189 | */
|
|---|
| 190 | unsigned long GetCoordUnit(int coordnum,unsigned long typ)
|
|---|
| 191 | {
|
|---|
| 192 | if(coordnum==1) {
|
|---|
| 193 | if (typ&TypCoord1H) return TypUniteH;
|
|---|
| 194 | else if(typ&TypCoord1D) return TypUniteD;
|
|---|
| 195 | else if(typ&TypCoord1R) return TypUniteR;
|
|---|
| 196 | else return TypUniteR;
|
|---|
| 197 | }
|
|---|
| 198 | if(coordnum==2) {
|
|---|
| 199 | if (typ&TypCoord2H) return TypUniteH;
|
|---|
| 200 | else if(typ&TypCoord2D) return TypUniteD;
|
|---|
| 201 | else if(typ&TypCoord2R) return TypUniteR;
|
|---|
| 202 | else return TypUniteR;
|
|---|
| 203 | }
|
|---|
| 204 | return TypCoordUndef;
|
|---|
| 205 | }
|
|---|
| 206 |
|
|---|
| 207 | /*! \ingroup XAstroPack
|
|---|
| 208 | \brief Pour decoder et transcrire en TypAstroCoord une chaine
|
|---|
| 209 | donnant la structure du systeme de coordonnees.
|
|---|
| 210 | \verbatim
|
|---|
| 211 | ctype = "CAaBb"
|
|---|
| 212 | C: type de coordonnees: E Equatoriales
|
|---|
| 213 | G Galactiques
|
|---|
| 214 | H Horizontales
|
|---|
| 215 | S Ecliptiques
|
|---|
| 216 | pas de defaut
|
|---|
| 217 | A: unite de la coordonnee 1 (alpha,longitude etc...)
|
|---|
| 218 | H heure
|
|---|
| 219 | D degre
|
|---|
| 220 | R radian
|
|---|
| 221 | defaut radian
|
|---|
| 222 | a: type d'intervalle pour la coordonnee 1
|
|---|
| 223 | C intervalle [0,24[ [0,360[ [0,2*Pi[
|
|---|
| 224 | L intervalle [-12,12[ [-180,180[ [-Pi,Pi[
|
|---|
| 225 | (defaut: C)
|
|---|
| 226 | A: unite de la coordonnee 2 (delta,latitude etc...)
|
|---|
| 227 | H heure
|
|---|
| 228 | D degre
|
|---|
| 229 | R radian
|
|---|
| 230 | defaut radian
|
|---|
| 231 | a: type d'intervalle pour la coordonnee 2
|
|---|
| 232 | C intervalle [0,12] [0,180] [0,Pi]
|
|---|
| 233 | (type colatitude)
|
|---|
| 234 | L intervalle [-6,6] [-90,90][ [-Pi/2,Pi/2]
|
|---|
| 235 | (defaut: L)
|
|---|
| 236 |
|
|---|
| 237 | Exemple: GDCDL : galactiques long=[0,360[ lat=[-90,90]
|
|---|
| 238 | GDxDx ou GDxD: idem
|
|---|
| 239 | Gxxxx ou G : galactiques long=[0,2*Pi[ lat=[-Pi/2,Pi/2]
|
|---|
| 240 | Exemple: EHCDL : equatoriales alpha=[0,24[ delta=[-90,90]
|
|---|
| 241 | EHxDx ou EHxD : idem
|
|---|
| 242 | Exxxx ou E : equatoriales alpha=[0,2*Pi[ delta=[-Pi/2,Pi/2]
|
|---|
| 243 |
|
|---|
| 244 | - Retourne 0 si probleme dans la chaine
|
|---|
| 245 | \endverbatim
|
|---|
| 246 | */
|
|---|
| 247 | unsigned long DecodeTypAstro(const char *ctype)
|
|---|
| 248 | {
|
|---|
| 249 | if(ctype==NULL) return TypCoordUndef;
|
|---|
| 250 | int len = strlen(ctype);
|
|---|
| 251 | if(len<1) return TypCoordUndef;
|
|---|
| 252 |
|
|---|
| 253 | unsigned long typ=0;
|
|---|
| 254 | // Le type de systeme de coordonnees
|
|---|
| 255 | int i = 0;
|
|---|
| 256 | if (ctype[i]=='e' || ctype[i]=='E') typ=TypCoordEq;
|
|---|
| 257 | else if(ctype[i]=='g' || ctype[i]=='G') typ=TypCoordGal;
|
|---|
| 258 | else if(ctype[i]=='h' || ctype[i]=='H') typ=TypCoordHor;
|
|---|
| 259 | else if(ctype[i]=='s' || ctype[i]=='S') typ=TypCoordEcl;
|
|---|
| 260 | else return TypCoordUndef;
|
|---|
| 261 | // La coordonnee 1
|
|---|
| 262 | i = 1;
|
|---|
| 263 | if(i>=len)
|
|---|
| 264 | {typ |= TypCoord1R|TypCoord1C|TypCoord2R|TypCoord2L; return typ;}
|
|---|
| 265 | if (ctype[i]=='h' || ctype[i]=='H') typ |= TypCoord1H;
|
|---|
| 266 | else if(ctype[i]=='d' || ctype[i]=='D') typ |= TypCoord1D;
|
|---|
| 267 | else if(ctype[i]=='r' || ctype[i]=='R') typ |= TypCoord1R;
|
|---|
| 268 | else typ |= TypCoord1R;
|
|---|
| 269 | i = 2;
|
|---|
| 270 | if(i>=len) {typ |= TypCoord1C|TypCoord2R|TypCoord2L; return typ;}
|
|---|
| 271 | if (ctype[i]=='c' || ctype[i]=='C') typ |= TypCoord1C;
|
|---|
| 272 | else if(ctype[i]=='l' || ctype[i]=='L') typ |= TypCoord1L;
|
|---|
| 273 | else typ |= TypCoord1C;
|
|---|
| 274 | // La coordonnee 2
|
|---|
| 275 | i = 3;
|
|---|
| 276 | if(i>=len) {typ |= TypCoord2R|TypCoord2L; return typ;}
|
|---|
| 277 | if (ctype[i]=='h' || ctype[i]=='H') typ |= TypCoord2H;
|
|---|
| 278 | else if(ctype[i]=='d' || ctype[i]=='D') typ |= TypCoord2D;
|
|---|
| 279 | else if(ctype[i]=='r' || ctype[i]=='R') typ |= TypCoord2R;
|
|---|
| 280 | else typ |= TypCoord2R;
|
|---|
| 281 | i = 4;
|
|---|
| 282 | if(i>=len) {typ |= TypCoord2L; return typ;}
|
|---|
| 283 | if (ctype[i]=='c' || ctype[i]=='C') typ |= TypCoord2C;
|
|---|
| 284 | else if(ctype[i]=='l' || ctype[i]=='L') typ |= TypCoord2L;
|
|---|
| 285 | else typ |= TypCoord2L;
|
|---|
| 286 | // Return
|
|---|
| 287 | return typ;
|
|---|
| 288 | }
|
|---|
| 289 |
|
|---|
| 290 | /*! \ingroup XAstroPack
|
|---|
| 291 | \brief Idem DecodeTypAstro(char *) mais a l'envers
|
|---|
| 292 | */
|
|---|
| 293 | string DecodeTypAstro(unsigned long typ)
|
|---|
| 294 | {
|
|---|
| 295 | string s = "";
|
|---|
| 296 |
|
|---|
| 297 | if (typ&TypCoordEq) s += "E";
|
|---|
| 298 | else if(typ&TypCoordGal) s += "G";
|
|---|
| 299 | else if(typ&TypCoordHor) s += "H";
|
|---|
| 300 | else if(typ&TypCoordEcl) s += "S";
|
|---|
| 301 | else s += "x";
|
|---|
| 302 |
|
|---|
| 303 | if (typ&TypCoord1H) s += "H";
|
|---|
| 304 | else if(typ&TypCoord1D) s += "D";
|
|---|
| 305 | else if(typ&TypCoord1R) s += "R";
|
|---|
| 306 | else s += "x";
|
|---|
| 307 |
|
|---|
| 308 | if (typ&TypCoord1C) s += "C";
|
|---|
| 309 | else if(typ&TypCoord1L) s += "L";
|
|---|
| 310 | else s += "x";
|
|---|
| 311 |
|
|---|
| 312 | if (typ&TypCoord2H) s += "H";
|
|---|
| 313 | else if(typ&TypCoord2D) s += "D";
|
|---|
| 314 | else if(typ&TypCoord2R) s += "R";
|
|---|
| 315 | else s += "x";
|
|---|
| 316 |
|
|---|
| 317 | if (typ&TypCoord2C) s += "C";
|
|---|
| 318 | else if(typ&TypCoord2L) s += "L";
|
|---|
| 319 | else s += "x";
|
|---|
| 320 |
|
|---|
| 321 | return s;
|
|---|
| 322 | }
|
|---|
| 323 |
|
|---|
| 324 | /*! \ingroup XAstroPack
|
|---|
| 325 | \brief Pour convertir la latitude en colatitude et vice-versa (in place)
|
|---|
| 326 | \verbatim
|
|---|
| 327 | val = valeur a convertir qui doit etre:
|
|---|
| 328 | si type "latitude" dans [-6,6] [-90,90] [-Pi/2,Pi/2]
|
|---|
| 329 | si type "colatitude" dans [0,12] [0,180] [0,Pi]
|
|---|
| 330 | typ = type d'unite: heure TypUniteH
|
|---|
| 331 | degre TypUniteD
|
|---|
| 332 | radian TypUniteR
|
|---|
| 333 | (Defaut: radian TypUniteR)
|
|---|
| 334 | \endverbatim
|
|---|
| 335 | */
|
|---|
| 336 | void ToCoLat(double *val,unsigned long typ)
|
|---|
| 337 | {
|
|---|
| 338 | if (typ&TypUniteH) *val = 6. - *val;
|
|---|
| 339 | else if(typ&TypUniteD) *val = 90. - *val;
|
|---|
| 340 | else if(typ&TypUniteR) *val = PI/2. - *val;
|
|---|
| 341 | else *val = PI/2. - *val;
|
|---|
| 342 | }
|
|---|
| 343 |
|
|---|
| 344 | /*! \ingroup XAstroPack
|
|---|
| 345 | \brief Pour remettre la valeur de la COLATITUDE "val" dans la dynamique [0.,range]
|
|---|
| 346 | \verbatim
|
|---|
| 347 | val = valeur a convertir qui doit etre mise dans [0,12] [0,180] [0,Pi]
|
|---|
| 348 | typ = type d'unite: heure TypUniteH
|
|---|
| 349 | degre TypUniteD
|
|---|
| 350 | radian TypUniteR
|
|---|
| 351 | ex en degre: 0 -> 0 , 90 -> 90 , 180 -> 180 , 270 -> 90 , 360 -> 0
|
|---|
| 352 | -90 -> 90 , -180 -> 180 , -270 -> 90 , -360 -> 0
|
|---|
| 353 | \endverbatim
|
|---|
| 354 | */
|
|---|
| 355 | void InRangeCoLat(double *val,unsigned long typ)
|
|---|
| 356 | {
|
|---|
| 357 | double range=PI;
|
|---|
| 358 | if(typ==TypUniteH) range=12.; else if(typ==TypUniteD) range=180.;
|
|---|
| 359 | InRange(val,2.*range,range);
|
|---|
| 360 | if(*val<0.) *val*=-1.;
|
|---|
| 361 | }
|
|---|
| 362 |
|
|---|
| 363 | /*! \ingroup XAstroPack
|
|---|
| 364 | \brief Pour remettre la valeur de la LATITUDE "val" dans la dynamique [-range/2,range/2]
|
|---|
| 365 | \verbatim
|
|---|
| 366 | val = valeur a convertir qui doit etre mise dans [-6,6] [-90,90] [-Pi/2,Pi/2]
|
|---|
| 367 | typ = type d'unite: heure TypUniteH
|
|---|
| 368 | degre TypUniteD
|
|---|
| 369 | radian TypUniteR
|
|---|
| 370 | ex en degre: 0 -> 0 , 90 -> 90 , 180 -> 0 , 270 -> -90 , 360 -> 0
|
|---|
| 371 | -90 -> -90 , -180 -> 0 , -270 -> 90 , -360 -> 0
|
|---|
| 372 | \endverbatim
|
|---|
| 373 | */
|
|---|
| 374 | void InRangeLat(double *val,unsigned long typ)
|
|---|
| 375 | {
|
|---|
| 376 | double range = PI;
|
|---|
| 377 | if(typ==TypUniteH) range = 12.; else if(typ==TypUniteD) range = 180.;
|
|---|
| 378 | InRange(val,2.*range,range);
|
|---|
| 379 | if(*val>range/2.) *val = range - *val;
|
|---|
| 380 | else if(*val<-range/2.) *val = -(range + *val);
|
|---|
| 381 | }
|
|---|
| 382 |
|
|---|
| 383 | /*! \ingroup XAstroPack
|
|---|
| 384 | \brief Compute MJD from date
|
|---|
| 385 | \verbatim
|
|---|
| 386 | MJD = modified Julian date (number of days elapsed since 1900 jan 0.5),
|
|---|
| 387 | dy is the decimale value of the day: dy = int(dy) + utc/24.
|
|---|
| 388 | \endverbatim
|
|---|
| 389 | */
|
|---|
| 390 | double MJDfrDate(double dy,int mn,int yr)
|
|---|
| 391 | {
|
|---|
| 392 | double mjd;
|
|---|
| 393 | cal_mjd(mn,dy,yr,&mjd);
|
|---|
| 394 | return mjd;
|
|---|
| 395 | }
|
|---|
| 396 |
|
|---|
| 397 | /*! \ingroup XAstroPack
|
|---|
| 398 | \brief Compute date from MJD
|
|---|
| 399 | */
|
|---|
| 400 | void DatefrMJD(double mjd,double *dy,int *mn,int *yr)
|
|---|
| 401 | {
|
|---|
| 402 | mjd_cal(mjd,mn,dy,yr);
|
|---|
| 403 | }
|
|---|
| 404 |
|
|---|
| 405 | /*! \ingroup XAstroPack
|
|---|
| 406 | \brief Given a mjd, return the year as a double.
|
|---|
| 407 | */
|
|---|
| 408 | double YearfrMJD(double mjd)
|
|---|
| 409 | {
|
|---|
| 410 | double yr;
|
|---|
| 411 | mjd_year(mjd,&yr);
|
|---|
| 412 | return yr;
|
|---|
| 413 | }
|
|---|
| 414 |
|
|---|
| 415 | /*! \ingroup XAstroPack
|
|---|
| 416 | \brief Given a decimal year, return mjd
|
|---|
| 417 | */
|
|---|
| 418 | double MJDfrYear(double yr)
|
|---|
| 419 | {
|
|---|
| 420 | double mjd;
|
|---|
| 421 | year_mjd(yr,&mjd);
|
|---|
| 422 | return mjd;
|
|---|
| 423 | }
|
|---|
| 424 |
|
|---|
| 425 | /*! \ingroup XAstroPack
|
|---|
| 426 | \brief Given a mjd, return the year and number of days since 00:00 Jan 1
|
|---|
| 427 | \warning: if mjd = 2 January -> number of days = 1
|
|---|
| 428 | */
|
|---|
| 429 | void YDfrMJD(double mjd,double *dy,int *yr)
|
|---|
| 430 | {
|
|---|
| 431 | mjd_dayno(mjd,yr,dy);
|
|---|
| 432 | }
|
|---|
| 433 |
|
|---|
| 434 | /*! \ingroup XAstroPack
|
|---|
| 435 | \brief Given a year,
|
|---|
| 436 | */
|
|---|
| 437 | int IsLeapYear(int y)
|
|---|
| 438 | {
|
|---|
| 439 | return isleapyear(y);
|
|---|
| 440 | }
|
|---|
| 441 |
|
|---|
| 442 | /*! \ingroup XAstroPack
|
|---|
| 443 | \brief given an mjd, set *dow to 0..6 according to which day of the week it falls on (0=sunday).
|
|---|
| 444 | \return return 0 if ok else -1 if can't figure it out.
|
|---|
| 445 | */
|
|---|
| 446 | int DayOrder(double mjd,int *dow)
|
|---|
| 447 | {
|
|---|
| 448 | return mjd_dow(mjd,dow);
|
|---|
| 449 | }
|
|---|
| 450 |
|
|---|
| 451 | /*! \ingroup XAstroPack
|
|---|
| 452 | \brief given a mjd, return the the number of days in the month.
|
|---|
| 453 | */
|
|---|
| 454 | int DaysInMonth(double mjd)
|
|---|
| 455 | {
|
|---|
| 456 | int ndays;
|
|---|
| 457 | mjd_dpm(mjd,&ndays);
|
|---|
| 458 | return ndays;
|
|---|
| 459 | }
|
|---|
| 460 |
|
|---|
| 461 | /*! \ingroup XAstroPack
|
|---|
| 462 | \brief Given a mjd, truncate it to the beginning of the whole day
|
|---|
| 463 | */
|
|---|
| 464 | double MJDat0hFrMJD(double mjd)
|
|---|
| 465 | {
|
|---|
| 466 | return mjd_day(mjd);
|
|---|
| 467 | }
|
|---|
| 468 |
|
|---|
| 469 | /*! \ingroup XAstroPack
|
|---|
| 470 | \brief Given a mjd, return the number of hours past midnight of the whole day
|
|---|
| 471 | */
|
|---|
| 472 | double HfrMJD(double mjd)
|
|---|
| 473 | {
|
|---|
| 474 | return mjd_hr(mjd);
|
|---|
| 475 | }
|
|---|
| 476 |
|
|---|
| 477 | /*! \ingroup XAstroPack
|
|---|
| 478 | \brief Give GST from UTC
|
|---|
| 479 | \verbatim
|
|---|
| 480 | Given a modified julian date, mjd, and a universally coordinated time, utc,
|
|---|
| 481 | return greenwich mean siderial time, *gst.
|
|---|
| 482 | N.B. mjd must be at the beginning of the day.
|
|---|
| 483 | \endverbatim
|
|---|
| 484 | */
|
|---|
| 485 | double GSTfrUTC(double mjd0,double utc)
|
|---|
| 486 | {
|
|---|
| 487 | double gst;
|
|---|
| 488 | utc_gst(mjd0,utc,&gst);
|
|---|
| 489 | return gst;
|
|---|
| 490 | }
|
|---|
| 491 |
|
|---|
| 492 | /*! \ingroup XAstroPack
|
|---|
| 493 | \brief Give UTC from GST
|
|---|
| 494 | \verbatim
|
|---|
| 495 | Given a modified julian date, mjd, and a greenwich mean siderial time, gst,
|
|---|
| 496 | return universally coordinated time, *utc.
|
|---|
| 497 | N.B. mjd must be at the beginning of the day.
|
|---|
| 498 | \endverbatim
|
|---|
| 499 | */
|
|---|
| 500 | double UTCfrGST(double mjd0,double gst)
|
|---|
| 501 | {
|
|---|
| 502 | double utc;
|
|---|
| 503 | gst_utc(mjd0,gst,&utc);
|
|---|
| 504 | return utc;
|
|---|
| 505 | }
|
|---|
| 506 |
|
|---|
| 507 | /*! \ingroup XAstroPack
|
|---|
| 508 | \brief Given apparent altitude find airmass.
|
|---|
| 509 | */
|
|---|
| 510 | double AirmassfrAlt(double alt)
|
|---|
| 511 | {
|
|---|
| 512 | double x;
|
|---|
| 513 | alt = degrad(alt);
|
|---|
| 514 | airmass(alt,&x);
|
|---|
| 515 | return x;
|
|---|
| 516 | }
|
|---|
| 517 |
|
|---|
| 518 | /*! \ingroup XAstroPack
|
|---|
| 519 | \brief given geocentric time "jd" and coords of a distant object at "ra/dec" (J2000),
|
|---|
| 520 | find the difference "hcp" in time between light arriving at earth vs the sun.
|
|---|
| 521 | \return "hcp" must be subtracted from "geocentric jd" to get "heliocentric jd".
|
|---|
| 522 | \warning "jd" is the TRUE Julian day (jd = mjd+MJD0).
|
|---|
| 523 | */
|
|---|
| 524 | double HelioCorr(double jd,double ra,double dec)
|
|---|
| 525 | {
|
|---|
| 526 | double hcp;
|
|---|
| 527 | ra=hrrad(ra);
|
|---|
| 528 | dec=degrad(dec);
|
|---|
| 529 | heliocorr(jd,ra,dec,&hcp);
|
|---|
| 530 | return hcp;
|
|---|
| 531 | }
|
|---|
| 532 |
|
|---|
| 533 | /*! \ingroup XAstroPack
|
|---|
| 534 | \brief gmst0() - return Greenwich Mean Sidereal Time at 0h UT
|
|---|
| 535 | \param mjd0 = date at 0h UT in julian days since MJD0
|
|---|
| 536 | */
|
|---|
| 537 | double GST0(double mjd0)
|
|---|
| 538 | /* Copie depuis le code de Xephem (utc_gst.c) car pas prototype*/
|
|---|
| 539 | {
|
|---|
| 540 | double T, x;
|
|---|
| 541 | T = ((int)(mjd0 - 0.5) + 0.5 - J2000)/36525.0;
|
|---|
| 542 | x = 24110.54841 +
|
|---|
| 543 | (8640184.812866 + (0.093104 - 6.2e-6 * T) * T) * T;
|
|---|
| 544 | x /= 3600.0;
|
|---|
| 545 | range(&x, 24.0);
|
|---|
| 546 | return (x);
|
|---|
| 547 | }
|
|---|
| 548 |
|
|---|
| 549 | /*! \ingroup XAstroPack
|
|---|
| 550 | \brief return local sidereal time from modified julian day and longitude
|
|---|
| 551 | \warning nutation or obliquity correction are taken into account.
|
|---|
| 552 | */
|
|---|
| 553 | double LSTfrMJD(double mjd,double geolng)
|
|---|
| 554 | {
|
|---|
| 555 | double eps,lst,deps,dpsi;
|
|---|
| 556 | utc_gst(mjd_day(mjd),mjd_hr(mjd),&lst);
|
|---|
| 557 | lst += deghr(geolng);
|
|---|
| 558 | obliquity(mjd,&eps);
|
|---|
| 559 | nutation(mjd,&deps,&dpsi);
|
|---|
| 560 | lst += radhr(dpsi*cos(eps+deps));
|
|---|
| 561 | InRange(&lst,24.);
|
|---|
| 562 | return lst;
|
|---|
| 563 | }
|
|---|
| 564 |
|
|---|
| 565 | /*! \ingroup XAstroPack
|
|---|
| 566 | \brief Give a time in h:mn:s from a decimal hour
|
|---|
| 567 | \verbatim
|
|---|
| 568 | // INPUT: hd
|
|---|
| 569 | // OUTPUT: h mn s (h,mn,s >=< 0)
|
|---|
| 570 | // REMARQUE: si hd<0 alors h<0 ET mn<0 ET s<0
|
|---|
| 571 | // EX: 12.51 -> h=12 mn=30 s=10 ;
|
|---|
| 572 | // -12.51 -> h=-12 mn=-30 s=-10 ;
|
|---|
| 573 | \endverbatim
|
|---|
| 574 | */
|
|---|
| 575 | void HMSfrHdec(double hd,int *h,int *mn,double *s)
|
|---|
| 576 | {
|
|---|
| 577 | int sgn=1;
|
|---|
| 578 | if(hd<0.) {sgn=-1; hd*=-1.;}
|
|---|
| 579 | *h = int(hd);
|
|---|
| 580 | *mn = int((hd-(double)(*h))*60.);
|
|---|
| 581 | *s = (hd - (double)(*h) - (double)(*mn)/60.)*3600.;
|
|---|
| 582 | // pb precision
|
|---|
| 583 | if(*s<0.) *s = 0.;
|
|---|
| 584 | if(*s>60. || *s==60.) {*s-=60.; *mn+=1;} // s=double attention comparaison
|
|---|
| 585 | if(*mn<0) *mn = 0;
|
|---|
| 586 | if(*mn>=60) {*mn-=60; *h+=1;}
|
|---|
| 587 | *h *= sgn; *mn *= sgn; *s *= (double)sgn;
|
|---|
| 588 | }
|
|---|
| 589 |
|
|---|
| 590 | /*! \ingroup XAstroPack
|
|---|
| 591 | \brief Give a decimal hour from a time in h:mn:s
|
|---|
| 592 | \verbatim
|
|---|
| 593 | // INPUT: h , mn , s (h,mn,s >=< 0)
|
|---|
| 594 | // RETURN: en heures decimales
|
|---|
| 595 | // REMARQUE: pour avoir hd=-12.51 <- h=-12 mn=-30 s=-10
|
|---|
| 596 | \endverbatim
|
|---|
| 597 | */
|
|---|
| 598 | double HdecfrHMS(int h,int mn,double s)
|
|---|
| 599 | {
|
|---|
| 600 | return ((double)h + (double)mn/60. + s/3600.);
|
|---|
| 601 | }
|
|---|
| 602 |
|
|---|
| 603 | /*! \ingroup XAstroPack
|
|---|
| 604 | \brief Give a time string from a time in h:mn:s
|
|---|
| 605 | \verbatim
|
|---|
| 606 | // INPUT: h , mn , s (h,mn,s >=< 0)
|
|---|
| 607 | // RETURN: string h:mn:s
|
|---|
| 608 | \endverbatim
|
|---|
| 609 | */
|
|---|
| 610 | string ToStringHMS(int h,int mn,double s)
|
|---|
| 611 | {
|
|---|
| 612 | double hd = HdecfrHMS(h,mn,s); // put in range
|
|---|
| 613 | HMSfrHdec(hd,&h,&mn,&s);
|
|---|
| 614 | char str[128];
|
|---|
| 615 | if(hd<0.)
|
|---|
| 616 | sprintf(str,"-%d:%d:%.3f",-h,-mn,-s);
|
|---|
| 617 | else
|
|---|
| 618 | sprintf(str,"%d:%d:%.3f",h,mn,s);
|
|---|
| 619 | string dum = str;
|
|---|
| 620 | return dum;
|
|---|
| 621 | }
|
|---|
| 622 |
|
|---|
| 623 | /*! \ingroup XAstroPack
|
|---|
| 624 | \brief Give a time string from a decimal hour
|
|---|
| 625 | */
|
|---|
| 626 | string ToStringHdec(double hd)
|
|---|
| 627 | {
|
|---|
| 628 | int h,mn; double s;
|
|---|
| 629 | HMSfrHdec(hd,&h,&mn,&s);
|
|---|
| 630 | return ToStringHMS(h,mn,s);
|
|---|
| 631 | }
|
|---|
| 632 |
|
|---|
| 633 | /*! \ingroup XAstroPack
|
|---|
| 634 | \brief Compute precession between 2 dates.
|
|---|
| 635 | */
|
|---|
| 636 | void Precess(double mjd1,double mjd2,double ra1,double dec1,double *ra2,double *dec2)
|
|---|
| 637 | {
|
|---|
| 638 | ra1 = hrrad(ra1); // radians
|
|---|
| 639 | dec1 = degrad(dec1); // radians
|
|---|
| 640 | precess(mjd1,mjd2,&ra1,&dec1);
|
|---|
| 641 | *ra2 = radhr(ra1); InRange(ra2,24.);
|
|---|
| 642 | *dec2 = raddeg(dec1);
|
|---|
| 643 | }
|
|---|
| 644 |
|
|---|
| 645 | /*! \ingroup XAstroPack
|
|---|
| 646 | \brief Convert equatorial coordinates for the given epoch into galactic coordinates
|
|---|
| 647 | */
|
|---|
| 648 | void EqtoGal(double mjd,double ra,double dec, double *glng,double *glat)
|
|---|
| 649 | // Coordonnees equatoriales -> Coordonnees galactiques
|
|---|
| 650 | {
|
|---|
| 651 | ra = hrrad(ra); // radians
|
|---|
| 652 | dec = degrad(dec); // radians
|
|---|
| 653 | eq_gal(mjd,ra,dec,glat,glng);
|
|---|
| 654 | // Vraiment bizarre, sur Linux-g++ glng>=360 ne comprend pas glng==360 ! (CMV)
|
|---|
| 655 | *glng = raddeg(*glng); InRange(glng,360.);
|
|---|
| 656 | *glat = raddeg(*glat);
|
|---|
| 657 | }
|
|---|
| 658 |
|
|---|
| 659 | /*! \ingroup XAstroPack
|
|---|
| 660 | \brief Convert galactic coordinates into equatorial coordinates at the given epoch
|
|---|
| 661 | */
|
|---|
| 662 | void GaltoEq(double mjd,double glng,double glat,double *ra,double *dec)
|
|---|
| 663 | // Coordonnees galactiques -> Coordonnees equatoriales
|
|---|
| 664 | {
|
|---|
| 665 | glng = degrad(glng); // radians
|
|---|
| 666 | glat = degrad(glat); // radians
|
|---|
| 667 | gal_eq (mjd,glat,glng,ra,dec);
|
|---|
| 668 | *ra = radhr(*ra); InRange(ra,24.);
|
|---|
| 669 | *dec = raddeg(*dec);
|
|---|
| 670 | }
|
|---|
| 671 |
|
|---|
| 672 | /*! \ingroup XAstroPack
|
|---|
| 673 | \brief Convert equatorial coordinates (with hour angle instead of right ascension) into horizontal coordinates.
|
|---|
| 674 | */
|
|---|
| 675 | void EqHtoHor(double geolat,double ha,double dec,double *az,double *alt)
|
|---|
| 676 | // Coordonnees equatoriales -> Coordonnees horizontales
|
|---|
| 677 | {
|
|---|
| 678 | geolat = degrad(geolat); // radians
|
|---|
| 679 | ha = hrrad(ha); // radians
|
|---|
| 680 | dec = degrad(dec); // radians
|
|---|
| 681 | hadec_aa (geolat,ha,dec,alt,az);
|
|---|
| 682 | *alt = raddeg(*alt);
|
|---|
| 683 | *az = raddeg(*az); InRange(az,360.);
|
|---|
| 684 | }
|
|---|
| 685 |
|
|---|
| 686 | /*! \ingroup XAstroPack
|
|---|
| 687 | Convert horizontal coordinates into equatorial coordinates (with hour angle instead of right ascension).
|
|---|
| 688 | */
|
|---|
| 689 | void HortoEqH(double geolat,double az,double alt,double *ha,double *dec)
|
|---|
| 690 | // Coordonnees horizontales -> Coordonnees equatoriales
|
|---|
| 691 | {
|
|---|
| 692 | geolat = degrad(geolat); // radians
|
|---|
| 693 | alt = degrad(alt); // radians
|
|---|
| 694 | az = degrad(az); // radians
|
|---|
| 695 | aa_hadec (geolat,alt,az,ha,dec);
|
|---|
| 696 | *ha = radhr(*ha); InRange(ha,24.,12.);
|
|---|
| 697 | *dec = raddeg(*dec);
|
|---|
| 698 | }
|
|---|
| 699 |
|
|---|
| 700 | /*! \ingroup XAstroPack
|
|---|
| 701 | \brief Convert equatorial coordinates into horizontal coordinates.
|
|---|
| 702 | */
|
|---|
| 703 | void EqtoHor(double geolat,double lst,double ra,double dec,double *az,double *alt)
|
|---|
| 704 | // Coordonnees equatoriales -> Coordonnees horizontales
|
|---|
| 705 | {
|
|---|
| 706 | double ha = lst - ra; InRange(&ha,24.,12.);
|
|---|
| 707 | geolat = degrad(geolat); // radians
|
|---|
| 708 | ha = hrrad(ha); // radians
|
|---|
| 709 | dec = degrad(dec); // radians
|
|---|
| 710 | hadec_aa (geolat,ha,dec,alt,az);
|
|---|
| 711 | *alt = raddeg(*alt);
|
|---|
| 712 | *az = raddeg(*az); InRange(az,360.);
|
|---|
| 713 | }
|
|---|
| 714 |
|
|---|
| 715 | /*! \ingroup XAstroPack
|
|---|
| 716 | Convert horizontal coordinates into equatorial coordinates.
|
|---|
| 717 | */
|
|---|
| 718 | void HortoEq(double geolat,double lst,double az,double alt,double *ra,double *dec)
|
|---|
| 719 | // Coordonnees horizontales -> Coordonnees equatoriales
|
|---|
| 720 | {
|
|---|
| 721 | double ha;
|
|---|
| 722 | geolat = degrad(geolat); // radians
|
|---|
| 723 | alt = degrad(alt); // radians
|
|---|
| 724 | az = degrad(az); // radians
|
|---|
| 725 | aa_hadec (geolat,alt,az,&ha,dec);
|
|---|
| 726 | *ra = lst - radhr(ha); InRange(ra,24.);
|
|---|
| 727 | *dec = raddeg(*dec);
|
|---|
| 728 | }
|
|---|
| 729 |
|
|---|
| 730 | /*! \ingroup XAstroPack
|
|---|
| 731 | \brief Convert equatorial coordinates into geocentric ecliptic coordinates given the modified Julian date.
|
|---|
| 732 | \warning Correction for the effect on the angle of the obliquity due to nutation is not included.
|
|---|
| 733 | */
|
|---|
| 734 | void EqtoEcl(double mjd,double ra,double dec,double *eclng,double *eclat)
|
|---|
| 735 | // Coordonnees equatoriales -> Coordonnees ecliptiques
|
|---|
| 736 | {
|
|---|
| 737 | ra = hrrad(ra); // radians
|
|---|
| 738 | dec = degrad(dec); // radians
|
|---|
| 739 | eq_ecl(mjd,ra,dec,eclat,eclng);
|
|---|
| 740 | *eclng = raddeg(*eclng); InRange(eclng,360.);
|
|---|
| 741 | *eclat = raddeg(*eclat);
|
|---|
| 742 | }
|
|---|
| 743 |
|
|---|
| 744 | /*! \ingroup XAstroPack
|
|---|
| 745 | \brief Convert geocentric ecliptic coordinates into equatorial coordinates given the modified Julian date.
|
|---|
| 746 | \warning Correction for the effect on the angle of the obliquity due to nutation is not included.
|
|---|
| 747 | */
|
|---|
| 748 | void EcltoEq(double mjd,double eclng,double eclat,double *ra,double *dec)
|
|---|
| 749 | // Coordonnees ecliptiques -> Coordonnees equatoriales
|
|---|
| 750 | {
|
|---|
| 751 | eclat = degrad(eclat); // radians
|
|---|
| 752 | eclng = degrad(eclng); // radians
|
|---|
| 753 | ecl_eq(mjd,eclat,eclng,ra,dec);
|
|---|
| 754 | *ra = radhr(*ra); InRange(ra,24.);
|
|---|
| 755 | *dec = raddeg(*dec);
|
|---|
| 756 | }
|
|---|
| 757 |
|
|---|
| 758 | /*! \ingroup XAstroPack
|
|---|
| 759 | \brief Give Sun position
|
|---|
| 760 | \verbatim
|
|---|
| 761 | given the modified JD, mjd, return the true geocentric ecliptic longitude
|
|---|
| 762 | of the sun for the mean equinox of the date, *eclsn, in degres, the
|
|---|
| 763 | sun-earth distance, *rsn, in AU, and the latitude *ecbsn, in degres
|
|---|
| 764 | (since this is always <= 1.2 arcseconds, in can be neglected by
|
|---|
| 765 | calling with ecbsn = NULL).
|
|---|
| 766 | - REMARQUE:
|
|---|
| 767 | * if the APPARENT ecliptic longitude is required, correct the longitude for
|
|---|
| 768 | * nutation to the true equinox of date and for aberration (light travel time,
|
|---|
| 769 | * approximately -9.27e7/186000/(3600*24*365)*2*pi = -9.93e-5 radians).
|
|---|
| 770 | \endverbatim
|
|---|
| 771 | */
|
|---|
| 772 | void SunPos(double mjd,double *eclsn,double *ecbsn,double *rsn)
|
|---|
| 773 | {
|
|---|
| 774 | sunpos(mjd,eclsn,rsn,ecbsn);
|
|---|
| 775 | *eclsn = raddeg(*eclsn); InRange(eclsn,360.);
|
|---|
| 776 | if(ecbsn!=NULL) *ecbsn = raddeg(*ecbsn);
|
|---|
| 777 | }
|
|---|
| 778 |
|
|---|
| 779 | /*! \ingroup XAstroPack
|
|---|
| 780 | \brief Give Moon position
|
|---|
| 781 | \verbatim
|
|---|
| 782 | given the mjd, find the geocentric ecliptic longitude, lam, and latitude,
|
|---|
| 783 | bet, and geocentric distance, rho in a.u. for the moon. also return
|
|---|
| 784 | the sun's mean anomaly, *msp, and the moon's mean anomaly, *mdp.
|
|---|
| 785 | (for the mean equinox)
|
|---|
| 786 | \endverbatim
|
|---|
| 787 | */
|
|---|
| 788 | void MoonPos(double mjd,double *eclmn,double *ecbmn,double *rho)
|
|---|
| 789 | {
|
|---|
| 790 | double msp,mdp;
|
|---|
| 791 | moon(mjd,eclmn,ecbmn,rho,&msp,&mdp);
|
|---|
| 792 | *eclmn = raddeg(*eclmn); InRange(eclmn,360.);
|
|---|
| 793 | *ecbmn = raddeg(*ecbmn);
|
|---|
| 794 | }
|
|---|
| 795 |
|
|---|
| 796 | /*! \ingroup XAstroPack
|
|---|
| 797 | \brief Give planet position
|
|---|
| 798 | \verbatim
|
|---|
| 799 | * given a modified Julian date, mjd, and a planet, p, find:
|
|---|
| 800 | * sunecl: heliocentric longitude,
|
|---|
| 801 | * sunecb: heliocentric latitude,
|
|---|
| 802 | * sundist: distance from the sun to the planet,
|
|---|
| 803 | * geodist: distance from the Earth to the planet,
|
|---|
| 804 | * none corrected for light time, ie, they are the true values for the
|
|---|
| 805 | * given instant.
|
|---|
| 806 | * geoecl: geocentric ecliptic longitude,
|
|---|
| 807 | * geoecb: geocentric ecliptic latitude,
|
|---|
| 808 | * each corrected for light time, ie, they are the apparent values as
|
|---|
| 809 | * seen from the center of the Earth for the given instant.
|
|---|
| 810 | * diamang: angular diameter in arcsec at 1 AU,
|
|---|
| 811 | * mag: visual magnitude when 1 AU from sun and earth at 0 phase angle.
|
|---|
| 812 | * (for the mean equinox)
|
|---|
| 813 | * all angles are in degres, all distances in AU.
|
|---|
| 814 | *
|
|---|
| 815 | * corrections for nutation and abberation must be made by the caller. The RA
|
|---|
| 816 | * and DEC calculated from the fully-corrected ecliptic coordinates are then
|
|---|
| 817 | * the apparent geocentric coordinates. Further corrections can be made, if
|
|---|
| 818 | * required, for atmospheric refraction and geocentric parallax.
|
|---|
| 819 | \endverbatim
|
|---|
| 820 | */
|
|---|
| 821 | void PlanetPos(double mjd,int numplan,double *sunecl,double *sunecb,double *sundist
|
|---|
| 822 | ,double *geodist,double *geoecl,double *geoecb
|
|---|
| 823 | ,double *diamang,double *mag)
|
|---|
| 824 | {
|
|---|
| 825 | plans(mjd,numplan,sunecl,sunecb,sundist,geodist,geoecl,geoecb,diamang,mag);
|
|---|
| 826 | *geoecl = raddeg(*geoecl); InRange(geoecl,360.);
|
|---|
| 827 | *geoecb = raddeg(*geoecb);
|
|---|
| 828 | *sunecl = raddeg(*sunecl); InRange(sunecl,360.);
|
|---|
| 829 | *sunecb = raddeg(*sunecb);
|
|---|
| 830 | }
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
