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

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

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