source: Sophya/trunk/SophyaExt/XAstroPack/xastropack.cc@ 1802

Last change on this file since 1802 was 1802, checked in by cmv, 24 years ago

un bug de documentation cmv 30/11/01

File size: 18.9 KB
Line 
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// HORIZONTAL: azimuth en degres [0,360[ (az)
44// (angle round to the east from north+)
45// altitude en degres [-90,90] (alt)
46// ECLIPTIQUE: lontitude ecliptique en degres [0,360[ (eclng)
47// (angle round counter clockwise from the vernal equinoxe)
48// latitude ecliptique en degres [-90,90] (eclat)
49// GEOGRAPHIE: longitude en degres ]-180,180] (geolng)
50// (angle <0 vers l'ouest, >0 vers l'est)
51// latitude en degres [-90,90] (north>0 sud<0) (geolat)
52 \endverbatim
53*/
54
55/*! \ingroup XAstroPack
56\brief Given a coordinate type "typ", convert to standard for astropack.
57\verbatim
58 La routine convertit (in place) les coordonnees "coord1","coord2"
59 definies par le type "typ" dans les unites standard de ce systeme
60 de coordonnees.
61 "typ" code le systeme de coordonnees astronomiques et les unites utilisees
62
63 - Return : 0 = OK
64 1 = Unknown type of coordinates
65 2 = bad range for coord1
66 4 = bad range for coord2
67 6 = bad range for coord1 et coord2
68
69 Les types de coordonnees sont definies dans le enum TypAstroCoord:
70 La premiere coordonnee est de type "longitude" (alpha,longitude)
71 La deuxieme coordonnee est de type "latidude" (delta,latitude)
72 *** Definitions des unites des coordonnees et de leurs dynamiques
73 - TypCoordH0 : heure=[0,24[
74 - TypCoordH1 : heure=]-12,12]
75 - TypCoordD0 : degre=[0,360[
76 - TypCoordD1 : degre=]-180,180]
77 - TypCoordD2 : degre=[-90,90]
78 - TypCoordR0 : degre=[0,2Pi[
79 - TypCoordR1 : degre=]-Pi,Pi]
80 - TypCoordR2 : degre=[-Pi/2,Pi/2]
81 *** Definitions des combinaisons unites des coordonnees
82 - TypCoordHD : coordonnees en (heure=[0,24[,degre=[-90,90])
83 - TypCoordDD : coordonnees en (degre=[0,360[,degre=[-90,90])
84 - TypCoordRR : coordonnees en (radian=[0,2Pi[,radian=[-Pi/2,Pi/2])
85 - TypCoordH1D : coordonnees en (heure=]-12,12],degre=[-90,90])
86 - TypCoordD1D : coordonnees en (degre=]-180,180],degre=[-90,90])
87 - TypCoordR1R : coordonnees en (radian=]-Pi,Pi],radian=[-Pi/2,Pi/2])
88 *** Definitions des types de systemes de coordonnees astronomiques.
89 - TypCoordEq : Coordonnees Equatoriales alpha,delta
90 - TypCoordGal : Coordonnees Galactiques gLong, gLat
91 - TypCoordHor : Coordonnees Horizontales azimuth,altitude
92 - TypCoordEcl : Coordonnees Ecliptiques EclLong,EclLat
93 *** Definitions des unites par defaut pour les divers systemes de coordonnees.
94 - TypCoordEqStd : heure=[0,24[, degre=[-90,90]
95 - TypCoordGalStd : degre=[0,360[,degre=[-90,90]
96 - TypCoordHorStd : degre=[0,360[,degre=[-90,90]
97 - TypCoordEclStd : degre=[0,360[,degre=[-90,90]
98\endverbatim
99*/
100int CoordConvertToStd(TypAstroCoord typ,double& coord1,double& coord2)
101{
102 int rc = 0;
103
104 // ---- Equatoriales alpha,delta
105 // - standard = [0,24[ , [-90,90]
106 if(typ&TypCoordEq) {
107 if(typ&TypCoordDD) {
108 coord1 = deghr(coord1);
109 } else if(typ&TypCoordRR) {
110 coord1 = radhr(coord1);
111 coord2 = raddeg(coord2);
112 }
113 if(coord1==24.) coord1 = 0.;
114 if(coord1<0. || coord1>=24.) rc+= 2;
115 if(coord2<-90. || coord2>90. ) rc+= 4;
116
117 // ---- Galactiques gLong, gLat
118 // ---- Horizontales azimuth,altitude
119 // ---- Ecliptiques EclLong,EclLat
120 // - standard = [0,360[ , [-90,90]
121 } else if( typ&TypCoordGal || typ&TypCoordHor || typ&TypCoordEcl) {
122 if(typ&TypCoordHD) {
123 coord1 = hrdeg(coord1);
124 } else if(typ&TypCoordRR) {
125 coord1 = raddeg(coord1);
126 coord2 = raddeg(coord2);
127 }
128 if(coord1==360.) coord1 = 0.;
129 if(coord1<0. || coord1>=360.) rc+= 2;
130 if(coord2<-90. || coord2>90. ) rc+= 4;
131
132 } else { // Coordonnees non-connues
133 rc= 1;
134 }
135
136 return rc;
137}
138
139/*! \ingroup XAstroPack
140\brief Compute MJD from date
141\verbatim
142 MJD = modified Julian date (number of days elapsed since 1900 jan 0.5),
143 dy is the decimale value of the day: dy = int(dy) + utc/24.
144\endverbatim
145*/
146double MJDfrDate(double dy,int mn,int yr)
147{
148double mjd;
149cal_mjd(mn,dy,yr,&mjd);
150return mjd;
151}
152
153/*! \ingroup XAstroPack
154\brief Compute date from MJD
155*/
156void DatefrMJD(double mjd,double *dy,int *mn,int *yr)
157{
158mjd_cal(mjd,mn,dy,yr);
159}
160
161/*! \ingroup XAstroPack
162\brief Given a mjd, return the year as a double.
163*/
164double YearfrMJD(double mjd)
165{
166double yr;
167mjd_year(mjd,&yr);
168return yr;
169}
170
171/*! \ingroup XAstroPack
172\brief Given a decimal year, return mjd
173*/
174double MJDfrYear(double yr)
175{
176double mjd;
177year_mjd(yr,&mjd);
178return mjd;
179}
180
181/*! \ingroup XAstroPack
182\brief Given a mjd, return the year and number of days since 00:00 Jan 1
183\warning: if mjd = 2 January -> number of days = 1
184*/
185void YDfrMJD(double mjd,double *dy,int *yr)
186{
187mjd_dayno(mjd,yr,dy);
188}
189
190/*! \ingroup XAstroPack
191\brief Given a year,
192*/
193int IsLeapYear(int y)
194{
195return isleapyear(y);
196}
197
198/*! \ingroup XAstroPack
199\brief given an mjd, set *dow to 0..6 according to which day of the week it falls on (0=sunday).
200\return return 0 if ok else -1 if can't figure it out.
201*/
202int DayOrder(double mjd,int *dow)
203{
204return mjd_dow(mjd,dow);
205}
206
207/*! \ingroup XAstroPack
208\brief given a mjd, return the the number of days in the month.
209*/
210int DaysInMonth(double mjd)
211{
212int ndays;
213mjd_dpm(mjd,&ndays);
214return ndays;
215}
216
217/*! \ingroup XAstroPack
218\brief Given a mjd, truncate it to the beginning of the whole day
219*/
220double MJDat0hFrMJD(double mjd)
221{
222return mjd_day(mjd);
223}
224
225/*! \ingroup XAstroPack
226\brief Given a mjd, return the number of hours past midnight of the whole day
227*/
228double HfrMJD(double mjd)
229{
230return mjd_hr(mjd);
231}
232
233/*! \ingroup XAstroPack
234\brief Give GST from UTC
235\verbatim
236 Given a modified julian date, mjd, and a universally coordinated time, utc,
237 return greenwich mean siderial time, *gst.
238 N.B. mjd must be at the beginning of the day.
239\endverbatim
240*/
241double GSTfrUTC(double mjd0,double utc)
242{
243double gst;
244utc_gst(mjd0,utc,&gst);
245return gst;
246}
247
248/*! \ingroup XAstroPack
249\brief Give UTC from GST
250\verbatim
251 Given a modified julian date, mjd, and a greenwich mean siderial time, gst,
252 return universally coordinated time, *utc.
253 N.B. mjd must be at the beginning of the day.
254\endverbatim
255*/
256double UTCfrGST(double mjd0,double gst)
257{
258double utc;
259gst_utc(mjd0,gst,&utc);
260return utc;
261}
262
263/*! \ingroup XAstroPack
264\brief Given apparent altitude find airmass.
265*/
266double AirmassfrAlt(double alt)
267{
268double x;
269alt = degrad(alt);
270airmass(alt,&x);
271return x;
272}
273
274/*! \ingroup XAstroPack
275\brief given geocentric time "jd" and coords of a distant object at "ra/dec" (J2000),
276find the difference "hcp" in time between light arriving at earth vs the sun.
277\return "hcp" must be subtracted from "geocentric jd" to get "heliocentric jd".
278\warning "jd" is the TRUE Julian day (jd = mjd+MJD0).
279*/
280double HelioCorr(double jd,double ra,double dec)
281{
282double hcp;
283ra=hrrad(ra);
284dec=degrad(dec);
285heliocorr(jd,ra,dec,&hcp);
286return hcp;
287}
288
289/*! \ingroup XAstroPack
290\brief gmst0() - return Greenwich Mean Sidereal Time at 0h UT
291\param mjd0 = date at 0h UT in julian days since MJD0
292*/
293double GST0(double mjd0)
294/* Copie depuis le code de Xephem (utc_gst.c) car pas prototype*/
295{
296 double T, x;
297 T = ((int)(mjd0 - 0.5) + 0.5 - J2000)/36525.0;
298 x = 24110.54841 +
299 (8640184.812866 + (0.093104 - 6.2e-6 * T) * T) * T;
300 x /= 3600.0;
301 range(&x, 24.0);
302 return (x);
303}
304
305/*! \ingroup XAstroPack
306\brief return local sidereal time from modified julian day and longitude
307\warning nutation or obliquity correction are taken into account.
308*/
309double LSTfrMJD(double mjd,double geolng)
310{
311 double eps,lst,deps,dpsi;
312 utc_gst(mjd_day(mjd),mjd_hr(mjd),&lst);
313 lst += deghr(geolng);
314 obliquity(mjd,&eps);
315 nutation(mjd,&deps,&dpsi);
316 lst += radhr(dpsi*cos(eps+deps));
317 InRange(&lst,24.);
318 return lst;
319}
320
321/*! \ingroup XAstroPack
322\brief Give a time in h:mn:s from a decimal hour
323\verbatim
324// INPUT: hd
325// OUTPUT: h mn s (h,mn,s >=< 0)
326// REMARQUE: si hd<0 alors h<0 ET mn<0 ET s<0
327// EX: 12.51 -> h=12 mn=30 s=10 ;
328// -12.51 -> h=-12 mn=-30 s=-10 ;
329\endverbatim
330*/
331void HMSfrHdec(double hd,int *h,int *mn,double *s)
332{
333 int sgn=1;
334 if(hd<0.) {sgn=-1; hd*=-1.;}
335 *h = int(hd);
336 *mn = int((hd-(double)(*h))*60.);
337 *s = (hd - (double)(*h) - (double)(*mn)/60.)*3600.;
338 // pb precision
339 if(*s<0.) *s = 0.;
340 if(*s>60. || *s==60.) {*s-=60.; *mn+=1;} // s=double attention comparaison
341 if(*mn<0) *mn = 0;
342 if(*mn>=60) {*mn-=60; *h+=1;}
343 *h *= sgn; *mn *= sgn; *s *= (double)sgn;
344}
345
346/*! \ingroup XAstroPack
347\brief Give a decimal hour from a time in h:mn:s
348\verbatim
349// INPUT: h , mn , s (h,mn,s >=< 0)
350// RETURN: en heures decimales
351// REMARQUE: pour avoir hd=-12.51 <- h=-12 mn=-30 s=-10
352\endverbatim
353*/
354double HdecfrHMS(int h,int mn,double s)
355{
356 return ((double)h + (double)mn/60. + s/3600.);
357}
358
359/*! \ingroup XAstroPack
360\brief Give a time string from a time in h:mn:s
361\verbatim
362// INPUT: h , mn , s (h,mn,s >=< 0)
363// RETURN: string h:mn:s
364\endverbatim
365*/
366string ToStringHMS(int h,int mn,double s)
367{
368 double hd = HdecfrHMS(h,mn,s); // put in range
369 HMSfrHdec(hd,&h,&mn,&s);
370 char str[128];
371 if(hd<0.)
372 sprintf(str,"-%d:%d:%.3f",-h,-mn,-s);
373 else
374 sprintf(str,"%d:%d:%.3f",h,mn,s);
375 string dum = str;
376 return dum;
377}
378
379/*! \ingroup XAstroPack
380\brief Give a time string from a decimal hour
381*/
382string ToStringHdec(double hd)
383{
384 int h,mn; double s;
385 HMSfrHdec(hd,&h,&mn,&s);
386 return ToStringHMS(h,mn,s);
387}
388
389/*! \ingroup XAstroPack
390\brief Compute precession between 2 dates.
391*/
392void Precess(double mjd1,double mjd2,double ra1,double dec1,double *ra2,double *dec2)
393{
394 ra1 = hrrad(ra1); // radians
395 dec1 = degrad(dec1); // radians
396 precess(mjd1,mjd2,&ra1,&dec1);
397 *ra2 = radhr(ra1); InRange(ra2,24.);
398 *dec2 = raddeg(dec1);
399}
400
401/*! \ingroup XAstroPack
402\brief Convert equatorial coordinates for the given epoch into galactic coordinates
403*/
404void EqtoGal(double mjd,double ra,double dec, double *glng,double *glat)
405// Coordonnees equatoriales -> Coordonnees galactiques
406{
407 ra = hrrad(ra); // radians
408 dec = degrad(dec); // radians
409 eq_gal(mjd,ra,dec,glat,glng);
410 // Vraiment bizarre, sur Linux-g++ glng>=360 ne comprend pas glng==360 ! (CMV)
411 *glng = raddeg(*glng); InRange(glng,360.);
412 *glat = raddeg(*glat);
413}
414
415/*! \ingroup XAstroPack
416\brief Convert galactic coordinates into equatorial coordinates at the given epoch
417*/
418void GaltoEq(double mjd,double glng,double glat,double *ra,double *dec)
419// Coordonnees galactiques -> Coordonnees equatoriales
420{
421 glng = degrad(glng); // radians
422 glat = degrad(glat); // radians
423 gal_eq (mjd,glat,glng,ra,dec);
424 *ra = radhr(*ra); InRange(ra,24.);
425 *dec = raddeg(*dec);
426}
427
428/*! \ingroup XAstroPack
429\brief Convert equatorial coordinates (with hour angle instead of right ascension) into horizontal coordinates.
430*/
431void EqHtoHor(double geolat,double ha,double dec,double *az,double *alt)
432// Coordonnees equatoriales -> Coordonnees horizontales
433{
434 geolat = degrad(geolat); // radians
435 ha = hrrad(ha); // radians
436 dec = degrad(dec); // radians
437 hadec_aa (geolat,ha,dec,alt,az);
438 *alt = raddeg(*alt);
439 *az = raddeg(*az); InRange(az,360.);
440}
441
442/*! \ingroup XAstroPack
443 Convert horizontal coordinates into equatorial coordinates (with hour angle instead of right ascension).
444*/
445void HortoEqH(double geolat,double az,double alt,double *ha,double *dec)
446// Coordonnees horizontales -> Coordonnees equatoriales
447{
448 geolat = degrad(geolat); // radians
449 alt = degrad(alt); // radians
450 az = degrad(az); // radians
451 aa_hadec (geolat,alt,az,ha,dec);
452 *ha = radhr(*ha); InRange(ha,24.,12.);
453 *dec = raddeg(*dec);
454}
455
456/*! \ingroup XAstroPack
457\brief Convert equatorial coordinates into horizontal coordinates.
458*/
459void EqtoHor(double geolat,double lst,double ra,double dec,double *az,double *alt)
460// Coordonnees equatoriales -> Coordonnees horizontales
461{
462 double ha = lst - ra; InRange(&ha,24.,12.);
463 geolat = degrad(geolat); // radians
464 ha = hrrad(ha); // radians
465 dec = degrad(dec); // radians
466 hadec_aa (geolat,ha,dec,alt,az);
467 *alt = raddeg(*alt);
468 *az = raddeg(*az); InRange(az,360.);
469}
470
471/*! \ingroup XAstroPack
472 Convert horizontal coordinates into equatorial coordinates.
473*/
474void HortoEq(double geolat,double lst,double az,double alt,double *ra,double *dec)
475// Coordonnees horizontales -> Coordonnees equatoriales
476{
477 double ha;
478 geolat = degrad(geolat); // radians
479 alt = degrad(alt); // radians
480 az = degrad(az); // radians
481 aa_hadec (geolat,alt,az,&ha,dec);
482 *ra = lst - radhr(ha); InRange(ra,24.);
483 *dec = raddeg(*dec);
484}
485
486/*! \ingroup XAstroPack
487\brief Convert equatorial coordinates into geocentric ecliptic coordinates given the modified Julian date.
488\warning Correction for the effect on the angle of the obliquity due to nutation is not included.
489*/
490void EqtoEcl(double mjd,double ra,double dec,double *eclng,double *eclat)
491// Coordonnees equatoriales -> Coordonnees ecliptiques
492{
493 ra = hrrad(ra); // radians
494 dec = degrad(dec); // radians
495 eq_ecl(mjd,ra,dec,eclat,eclng);
496 *eclng = raddeg(*eclng); InRange(eclng,360.);
497 *eclat = raddeg(*eclat);
498}
499
500/*! \ingroup XAstroPack
501\brief Convert geocentric ecliptic coordinates into equatorial coordinates given the modified Julian date.
502\warning Correction for the effect on the angle of the obliquity due to nutation is not included.
503*/
504void EcltoEq(double mjd,double eclng,double eclat,double *ra,double *dec)
505// Coordonnees ecliptiques -> Coordonnees equatoriales
506{
507 eclat = degrad(eclat); // radians
508 eclng = degrad(eclng); // radians
509 ecl_eq(mjd,eclat,eclng,ra,dec);
510 *ra = radhr(*ra); InRange(ra,24.);
511 *dec = raddeg(*dec);
512}
513
514/*! \ingroup XAstroPack
515\brief Give Sun position
516\verbatim
517 given the modified JD, mjd, return the true geocentric ecliptic longitude
518 of the sun for the mean equinox of the date, *eclsn, in degres, the
519 sun-earth distance, *rsn, in AU, and the latitude *ecbsn, in degres
520 (since this is always <= 1.2 arcseconds, in can be neglected by
521 calling with ecbsn = NULL).
522 - REMARQUE:
523 * if the APPARENT ecliptic longitude is required, correct the longitude for
524 * nutation to the true equinox of date and for aberration (light travel time,
525 * approximately -9.27e7/186000/(3600*24*365)*2*pi = -9.93e-5 radians).
526\endverbatim
527*/
528void SunPos(double mjd,double *eclsn,double *ecbsn,double *rsn)
529{
530 sunpos(mjd,eclsn,rsn,ecbsn);
531 *eclsn = raddeg(*eclsn); InRange(eclsn,360.);
532 if(ecbsn!=NULL) *ecbsn = raddeg(*ecbsn);
533}
534
535/*! \ingroup XAstroPack
536\brief Give Moon position
537\verbatim
538 given the mjd, find the geocentric ecliptic longitude, lam, and latitude,
539 bet, and geocentric distance, rho in a.u. for the moon. also return
540 the sun's mean anomaly, *msp, and the moon's mean anomaly, *mdp.
541 (for the mean equinox)
542\endverbatim
543*/
544void MoonPos(double mjd,double *eclmn,double *ecbmn,double *rho)
545{
546 double msp,mdp;
547 moon(mjd,eclmn,ecbmn,rho,&msp,&mdp);
548 *eclmn = raddeg(*eclmn); InRange(eclmn,360.);
549 *ecbmn = raddeg(*ecbmn);
550}
551
552/*! \ingroup XAstroPack
553\brief Give planet position
554\verbatim
555 * given a modified Julian date, mjd, and a planet, p, find:
556 * sunecl: heliocentric longitude,
557 * sunecb: heliocentric latitude,
558 * sundist: distance from the sun to the planet,
559 * geodist: distance from the Earth to the planet,
560 * none corrected for light time, ie, they are the true values for the
561 * given instant.
562 * geoecl: geocentric ecliptic longitude,
563 * geoecb: geocentric ecliptic latitude,
564 * each corrected for light time, ie, they are the apparent values as
565 * seen from the center of the Earth for the given instant.
566 * diamang: angular diameter in arcsec at 1 AU,
567 * mag: visual magnitude when 1 AU from sun and earth at 0 phase angle.
568 * (for the mean equinox)
569 * all angles are in degres, all distances in AU.
570 *
571 * corrections for nutation and abberation must be made by the caller. The RA
572 * and DEC calculated from the fully-corrected ecliptic coordinates are then
573 * the apparent geocentric coordinates. Further corrections can be made, if
574 * required, for atmospheric refraction and geocentric parallax.
575 \endverbatim
576*/
577void PlanetPos(double mjd,int numplan,double *sunecl,double *sunecb,double *sundist
578 ,double *geodist,double *geoecl,double *geoecb
579 ,double *diamang,double *mag)
580{
581 plans(mjd,numplan,sunecl,sunecb,sundist,geodist,geoecl,geoecb,diamang,mag);
582 *geoecl = raddeg(*geoecl); InRange(geoecl,360.);
583 *geoecb = raddeg(*geoecb);
584 *sunecl = raddeg(*sunecl); InRange(sunecl,360.);
585 *sunecb = raddeg(*sunecb);
586}
Note: See TracBrowser for help on using the repository browser.