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

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

InRange,degrad()... et inline cmv 10/10/2001

File size: 13.3 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 gmst0() - return Greenwich Mean Sidereal Time at 0h UT
38\param mjd0 = date at 0h UT in julian days since MJD0
39*/
40double GST0(double mjd0)
41/* Copie depuis le code de Xephem (utc_gst.c) car pas prototype*/
42{
43 double T, x;
44 T = ((int)(mjd0 - 0.5) + 0.5 - J2000)/36525.0;
45 x = 24110.54841 +
46 (8640184.812866 + (0.093104 - 6.2e-6 * T) * T) * T;
47 x /= 3600.0;
48 range(&x, 24.0);
49 return (x);
50}
51
52/*! \ingroup XAstroPack
53\brief return local sidereal time from modified julian day and longitude
54\warning nutation or obliquity correction are taken into account.
55*/
56double LSTfrMJD(double mjd,double geolng)
57{
58 double eps,lst,deps,dpsi;
59 utc_gst(mjd_day(mjd),mjd_hr(mjd),&lst);
60 lst += deghr(geolng);
61 obliquity(mjd,&eps);
62 nutation(mjd,&deps,&dpsi);
63 lst += radhr(dpsi*cos(eps+deps));
64 InRange(&lst,24.);
65 return lst;
66}
67
68/*! \ingroup XAstroPack
69\brief Give a time in h:mn:s from a decimal hour
70\verbatim
71// INPUT: hd
72// OUTPUT: h mn s (h,mn,s >=< 0)
73// REMARQUE: si hd<0 alors h<0 ET mn<0 ET s<0
74// EX: 12.51 -> h=12 mn=30 s=10 ;
75// -12.51 -> h=-12 mn=-30 s=-10 ;
76\endverbatim
77*/
78void HMSfrHdec(double hd,int *h,int *mn,double *s)
79{
80 int sgn=1;
81 if(hd<0.) {sgn=-1; hd*=-1.;}
82 *h = int(hd);
83 *mn = int((hd-(double)(*h))*60.);
84 *s = (hd - (double)(*h) - (double)(*mn)/60.)*3600.;
85 // pb precision
86 if(*s<0.) *s = 0.;
87 if(*s>60. || *s==60.) {*s-=60.; *mn+=1;} // s=double attention comparaison
88 if(*mn<0) *mn = 0;
89 if(*mn>=60) {*mn-=60; *h+=1;}
90 *h *= sgn; *mn *= sgn; *s *= (double)sgn;
91}
92
93/*! \ingroup XAstroPack
94\brief Give a decimal hour from a time in h:mn:s
95\verbatim
96// INPUT: h , mn , s (h,mn,s >=< 0)
97// RETURN: en heures decimales
98// REMARQUE: pour avoir hd=-12.51 <- h=-12 mn=-30 s=-10
99\endverbatim
100*/
101double HdecfrHMS(int h,int mn,double s)
102{
103 return ((double)h + (double)mn/60. + s/3600.);
104}
105
106/*! \ingroup XAstroPack
107\brief Give a time string from a time in h:mn:s
108\verbatim
109// INPUT: h , mn , s (h,mn,s >=< 0)
110// RETURN: string h:mn:s
111\endverbatim
112*/
113string ToStringHMS(int h,int mn,double s)
114{
115 double hd = HdecfrHMS(h,mn,s); // put in range
116 HMSfrHdec(hd,&h,&mn,&s);
117 char str[128];
118 if(hd<0.)
119 sprintf(str,"-%d:%d:%.3f",-h,-mn,-s);
120 else
121 sprintf(str,"%d:%d:%.3f",h,mn,s);
122 string dum = str;
123 return dum;
124}
125
126/*! \ingroup XAstroPack
127\brief Give a time string from a decimal hour
128*/
129string ToStringHdec(double hd)
130{
131 int h,mn; double s;
132 HMSfrHdec(hd,&h,&mn,&s);
133 return ToStringHMS(h,mn,s);
134}
135
136/*! \ingroup XAstroPack
137\brief Compute precession between 2 dates.
138*/
139void Precess(double mjd1,double mjd2,double ra1,double dec1,double *ra2,double *dec2)
140{
141 ra1 = hrrad(ra1); // radians
142 dec1 = degrad(dec1); // radians
143 precess(mjd1,mjd2,&ra1,&dec1);
144 *ra2 = radhr(ra1); InRange(ra2,24.);
145 *dec2 = raddeg(dec1);
146}
147
148/*! \ingroup XAstroPack
149\brief Convert equatorial coordinates for the given epoch into galactic coordinates
150*/
151void EqtoGal(double mjd,double ra,double dec, double *glng,double *glat)
152// Coordonnees equatoriales -> Coordonnees galactiques
153{
154 ra = hrrad(ra); // radians
155 dec = degrad(dec); // radians
156 eq_gal(mjd,ra,dec,glat,glng);
157 // Vraiment bizarre, sur Linux-g++ glng>=360 ne comprend pas glng==360 ! (CMV)
158 *glng = raddeg(*glng); InRange(glng,360.);
159 *glat = raddeg(*glat);
160}
161
162/*! \ingroup XAstroPack
163\brief Convert galactic coordinates into equatorial coordinates at the given epoch
164*/
165void GaltoEq(double mjd,double glng,double glat,double *ra,double *dec)
166// Coordonnees galactiques -> Coordonnees equatoriales
167{
168 glng = degrad(glng); // radians
169 glat = degrad(glat); // radians
170 gal_eq (mjd,glat,glng,ra,dec);
171 *ra = radhr(*ra); InRange(ra,24.);
172 *dec = raddeg(*dec);
173}
174
175/*! \ingroup XAstroPack
176\brief Convert equatorial coordinates (with hour angle instead of right ascension) into horizontal coordinates.
177*/
178void EqHtoHor(double geolat,double ha,double dec,double *az,double *alt)
179// Coordonnees equatoriales -> Coordonnees horizontales
180{
181 geolat = degrad(geolat); // radians
182 ha = hrrad(ha); // radians
183 dec = degrad(dec); // radians
184 hadec_aa (geolat,ha,dec,alt,az);
185 *alt = raddeg(*alt);
186 *az = raddeg(*az); InRange(az,360.);
187}
188
189/*! \ingroup XAstroPack
190 Convert horizontal coordinates into equatorial coordinates (with hour angle instead of right ascension).
191*/
192void HortoEqH(double geolat,double az,double alt,double *ha,double *dec)
193// Coordonnees horizontales -> Coordonnees equatoriales
194{
195 geolat = degrad(geolat); // radians
196 alt = degrad(alt); // radians
197 az = degrad(az); // radians
198 aa_hadec (geolat,alt,az,ha,dec);
199 *ha = radhr(*ha); InRange(ha,24.,12.);
200 *dec = raddeg(*dec);
201}
202
203/*! \ingroup XAstroPack
204\brief Convert equatorial coordinates into horizontal coordinates.
205*/
206void EqtoHor(double geolat,double lst,double ra,double dec,double *az,double *alt)
207// Coordonnees equatoriales -> Coordonnees horizontales
208{
209 double ha = lst - ra; InRange(&ha,24.,12.);
210 geolat = degrad(geolat); // radians
211 ha = hrrad(ha); // radians
212 dec = degrad(dec); // radians
213 hadec_aa (geolat,ha,dec,alt,az);
214 *alt = raddeg(*alt);
215 *az = raddeg(*az); InRange(az,360.);
216}
217
218/*! \ingroup XAstroPack
219 Convert horizontal coordinates into equatorial coordinates.
220*/
221void HortoEq(double geolat,double lst,double az,double alt,double *ra,double *dec)
222// Coordonnees horizontales -> Coordonnees equatoriales
223{
224 double ha;
225 geolat = degrad(geolat); // radians
226 alt = degrad(alt); // radians
227 az = degrad(az); // radians
228 aa_hadec (geolat,alt,az,&ha,dec);
229 *ra = lst - radhr(ha); InRange(ra,24.);
230 *dec = raddeg(*dec);
231}
232
233/*! \ingroup XAstroPack
234\brief Convert equatorial coordinates into geocentric ecliptic coordinates given the modified Julian date.
235\warning Correction for the effect on the angle of the obliquity due to nutation is not included.
236*/
237// Attention, j'ai modifie eq_ecl.c pour proteger NaN
238// dans ecleq_aux :
239// *q = (sy*ceps)-(cy*seps*sx*sw);
240// if(*q<-1.) *q = -PI/2.; else if(*q>1.) *q = PI/2.; else *q = asin(*q);
241void EqtoEcl(double mjd,double ra,double dec,double *eclng,double *eclat)
242// Coordonnees equatoriales -> Coordonnees ecliptiques
243{
244 ra = hrrad(ra); // radians
245 dec = degrad(dec); // radians
246 eq_ecl(mjd,ra,dec,eclat,eclng);
247 *eclng = raddeg(*eclng); InRange(eclng,360.);
248 *eclat = raddeg(*eclat);
249}
250
251/*! \ingroup XAstroPack
252\brief Convert geocentric ecliptic coordinates into equatorial coordinates given the modified Julian date.
253\warning Correction for the effect on the angle of the obliquity due to nutation is not included.
254*/
255void EcltoEq(double mjd,double eclng,double eclat,double *ra,double *dec)
256// Coordonnees ecliptiques -> Coordonnees equatoriales
257{
258 eclat = degrad(eclat); // radians
259 eclng = degrad(eclng); // radians
260 ecl_eq(mjd,eclat,eclng,ra,dec);
261 *ra = radhr(*ra); InRange(ra,24.);
262 *dec = raddeg(*dec);
263}
264
265/*! \ingroup XAstroPack
266\brief Give Sun position
267\verbatim
268 given the modified JD, mjd, return the true geocentric ecliptic longitude
269 of the sun for the mean equinox of the date, *eclsn, in degres, the
270 sun-earth distance, *rsn, in AU, and the latitude *ecbsn, in degres
271 (since this is always <= 1.2 arcseconds, in can be neglected by
272 calling with ecbsn = NULL).
273 - REMARQUE:
274 * if the APPARENT ecliptic longitude is required, correct the longitude for
275 * nutation to the true equinox of date and for aberration (light travel time,
276 * approximately -9.27e7/186000/(3600*24*365)*2*pi = -9.93e-5 radians).
277\endverbatim
278*/
279void SunPos(double mjd,double *eclsn,double *ecbsn,double *rsn)
280{
281 sunpos(mjd,eclsn,rsn,ecbsn);
282 *eclsn = raddeg(*eclsn); InRange(eclsn,360.);
283 if(ecbsn!=NULL) *ecbsn = raddeg(*ecbsn);
284}
285
286/*! \ingroup XAstroPack
287\brief Give Moon position
288\verbatim
289 given the mjd, find the geocentric ecliptic longitude, lam, and latitude,
290 bet, and geocentric distance, rho in a.u. for the moon. also return
291 the sun's mean anomaly, *msp, and the moon's mean anomaly, *mdp.
292 (for the mean equinox)
293\endverbatim
294*/
295void MoonPos(double mjd,double *eclmn,double *ecbmn,double *rho)
296{
297 double msp,mdp;
298 moon(mjd,eclmn,ecbmn,rho,&msp,&mdp);
299 *eclmn = raddeg(*eclmn); InRange(eclmn,360.);
300 *ecbmn = raddeg(*ecbmn);
301}
302
303/*! \ingroup XAstroPack
304\brief Give planet position
305\verbatim
306 * given a modified Julian date, mjd, and a planet, p, find:
307 * sunecl: heliocentric longitude,
308 * sunecb: heliocentric latitude,
309 * sundist: distance from the sun to the planet,
310 * geodist: distance from the Earth to the planet,
311 * none corrected for light time, ie, they are the true values for the
312 * given instant.
313 * geoecl: geocentric ecliptic longitude,
314 * geoecb: geocentric ecliptic latitude,
315 * each corrected for light time, ie, they are the apparent values as
316 * seen from the center of the Earth for the given instant.
317 * diamang: angular diameter in arcsec at 1 AU,
318 * mag: visual magnitude when 1 AU from sun and earth at 0 phase angle.
319 * (for the mean equinox)
320 * all angles are in degres, all distances in AU.
321 *
322 * corrections for nutation and abberation must be made by the caller. The RA
323 * and DEC calculated from the fully-corrected ecliptic coordinates are then
324 * the apparent geocentric coordinates. Further corrections can be made, if
325 * required, for atmospheric refraction and geocentric parallax.
326 \endverbatim
327*/
328void PlanetPos(double mjd,int numplan,double *sunecl,double *sunecb,double *sundist
329 ,double *geodist,double *geoecl,double *geoecb
330 ,double *diamang,double *mag)
331{
332 plans(mjd,numplan,sunecl,sunecb,sundist,geodist,geoecl,geoecb,diamang,mag);
333 *geoecl = raddeg(*geoecl); InRange(geoecl,360.);
334 *geoecb = raddeg(*geoecb);
335 *sunecl = raddeg(*sunecl); InRange(sunecl,360.);
336 *sunecb = raddeg(*sunecb);
337}
338
339/*! \ingroup XAstroPack
340\brief Given a coordinate type "typ", convert to standard for astropack
341\verbatim
342// Return : 0 = OK
343// 1 = Unknown type of coordinates
344// 2 = bad range for coord1
345// 4 = bad range for coord2
346// 6 = bad range for coord1 et coord2
347\endverbatim
348*/
349int CoordConvertToStd(TypAstroCoord typ,double& coord1,double& coord2)
350{
351 int rc = 0;
352
353 // ---- Equatoriales alpha,delta
354 // - standard = [0,24[ , [-90,90]
355 if(typ&TypCoordEq) {
356 if(typ&TypCoordDD) {
357 coord1 = deghr(coord1);
358 } else if(typ&TypCoordRR) {
359 coord1 = radhr(coord1);
360 coord2 = raddeg(coord2);
361 }
362 if(coord1==24.) coord1 = 0.;
363 if(coord1<0. || coord1>=24.) rc+= 2;
364 if(coord2<-90. || coord2>90. ) rc+= 4;
365
366 // ---- Galactiques gLong, gLat
367 // ---- Horizontales azimuth,altitude
368 // ---- Ecliptiques EclLong,EclLat
369 // - standard = [0,360[ , [-90,90]
370 } else if( typ&TypCoordGal || typ&TypCoordHor || typ&TypCoordEcl) {
371 if(typ&TypCoordHD) {
372 coord1 = hrdeg(coord1);
373 } else if(typ&TypCoordRR) {
374 coord1 = raddeg(coord1);
375 coord2 = raddeg(coord2);
376 }
377 if(coord1==360.) coord1 = 0.;
378 if(coord1<0. || coord1>=360.) rc+= 2;
379 if(coord2<-90. || coord2>90. ) rc+= 4;
380
381 } else { // Coordonnees non-connues
382 rc= 1;
383 }
384
385 return rc;
386}
Note: See TracBrowser for help on using the repository browser.