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

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

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