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

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

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