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