#include <math.h>
#include <stdio.h>
#include "xastropack.h"
// TEMPS: modified Julian date (mjd) (number of days elapsed since 1900 jan 0.5)
// jour [1,31] (dy)
// mois [1,12] (mn)
// annee (yr)
// universal time [0,24[ (utc)
// Greenwich mean siderial [0,24[ (gst)
// Greenwich mean siderial at 0h UT [0,24[ (gst0)
// EQUATORIALE: ascension droite en heures [0,24[ (ra)
// declinaison en degres [-90,90] (dec)
// angle horaire en heures [-12,12] (-12=12) (ha) tsid=ha+ra
// GALACTIQUE: longitude en degres [0,360[ (glng)
// latitude en degres [-90,90] (glat)
// HORIZONTAL: azimuth en degres [0,360[ (az)
// (angle round to the east from north+)
// altitude en degres [-90,90] (alt)
// ECLIPTIQUE: lontitude ecliptique en degres [0,360[ (eclng)
// (angle round counter clockwise from the vernal equinoxe)
// latitude ecliptique en degres [-90,90] (eclat)
// GEOGRAPHIE: longitude en degres ]-180,180] (geolng)
// (angle + vers l'ouest, - vers l'est)
// latitude en degres [-90,90] (north>0) (geolat)
double TrueJDfrMJD(double mjd)
{
return mjd + MJD0;
}
double MJDfrTrueJD(double jd)
{
return jd - MJD0;
}
double MJDfrDate(double dy,int mn,int yr)
{
double mjd;
cal_mjd(mn,dy,yr,&mjd);
return mjd;
}
void DatefrMJD(double mjd,double *dy,int *mn,int *yr)
{
mjd_cal(mjd,mn,dy,yr);
}
/* given a mjd, return the year as a double. */
double YearfrMJD(double mjd)
{
double yr;
mjd_year(mjd,&yr);
return yr;
}
/* given a decimal year, return mjd */
double MJDfrYear(double yr)
{
double mjd;
year_mjd(yr,&mjd);
return mjd;
}
/* given a mjd, return the year and number of days since 00:00 Jan 1 */
/* Attention: si mjd = 2 Janvier -> number of days = 1 */
void YDfrMJD(double mjd,double *dy,int *yr)
{
mjd_dayno(mjd,yr,dy);
}
/* given a modified julian date, mjd, and a universally coordinated time, utc,
* return greenwich mean siderial time, *gst.
* N.B. mjd must be at the beginning of the day.
*/
double GSTfrUTC(double mjd0,double utc)
{
double gst;
utc_gst(mjd0,utc,&gst) ;
return gst;
}
/* given a modified julian date, mjd, and a greenwich mean siderial time, gst,
* return universally coordinated time, *utc.
* N.B. mjd must be at the beginning of the day.
*/
double UTCfrGST(double mjd0,double gst)
{
double utc;
gst_utc(mjd0,gst,&utc);
return utc;
}
/* gmst0() - return Greenwich Mean Sidereal Time at 0h UT */
/* mjd = date at 0h UT in julian days since MJD0 */
double GST0(double mjd0)
/* Copie depuis le code de Xephem car pas prototype */
{
double T, x;
T = ((int)(mjd0 - 0.5) + 0.5 - J2000)/36525.0;
x = 24110.54841 +
(8640184.812866 + (0.093104 - 6.2e-6 * T) * T) * T;
x /= 3600.0;
range(&x, 24.0);
return (x);
}
void Precess(double mjd1,double mjd2,double ra1,double dec1,double *ra2,double *dec2)
{
ra1 *= PI/12.; // radians
dec1 *= PI/180.; // radians
precess(mjd1,mjd2,&ra1,&dec1);
*ra2 = ra1*12./PI; if(*ra2>24.) *ra2 -= 24.; if(*ra2==24.) *ra2 = 0.;
*dec2 = dec1*180./PI;
}
/* given apparent altitude find airmass. */
double AirmassfrAlt(double alt)
{
double x;
alt *= PI/180.; // radians
airmass(alt,&x);
return x;
}
/* donne l'angle horaire a partir du temps sideral et de l'ascension droite */
double HafrRaTS(double gst,double ra)
{
double ha = gst - ra;
// Attention au probleme de la discontinuite 0h <==> 24h
// ts=1 ra=23 ; (ts-ra)=-22 <-12 --> ha = +2 = +24 + (ts-ra)
// ts=23 ra=1 ; (ts-ra)=+22 >+12 --> ha = -2 = -24 + (ts-ra)
if(ha==-12.) ha = 12.; if(ha<-12.) ha += 24.; if(ha>12.) ha -= 24.;
return ha;
}
void HMSfrHdec(double hd,int *h,int *mn,double *s)
// INPUT: hd
// OUTPUT: h mn s (h,mn,s >=< 0)
// REMARQUE: si hd<0 alors h<0 ET mn<0 ET s<0
// EX: 12.51 -> h=12 mn=30 s=10 ;
// -12.51 -> h=-12 mn=-30 s=-10 ;
{
int sgn=1;
if(hd<0.) {sgn=-1; hd*=-1.;}
*h = int(hd);
*mn = int((hd-(double)(*h))*60.);
*s = (hd - (double)(*h) - (double)(*mn)/60.)*3600.;
// pb precision
if(*s<0.) *s = 0.;
if(*s>60. || *s==60.) {*s-=60.; *mn+=1;} // s=double attention comparaison
if(*mn<0) *mn = 0;
if(*mn>=60) {*mn-=60; *h+=1;}
*h *= sgn; *mn *= sgn; *s *= (double)sgn;
}
double HdecfrHMS(int h,int mn,double s)
// INPUT: h , mn , s (h,mn,s >=< 0)
// RETURN: en heures decimales
// REMARQUE: pour avoir hd=-12.51 <- h=-12 mn=-30 s=-10
{
return ((double)h + (double)mn/60. + s/3600.);
}
string ToStringHMS(int h,int mn,double s)
// INPUT: h , mn , s (h,mn,s >=< 0)
// RETURN: string h:mn:s
{
double hd = HdecfrHMS(h,mn,s); // put in range
HMSfrHdec(hd,&h,&mn,&s);
char str[128];
if(hd<0.)
sprintf(str,"-%d:%d:%.3f",-h,-mn,-s);
else
sprintf(str,"%d:%d:%.3f",h,mn,s);
string dum = str;
return dum;
}
string ToStringHdec(double hd)
{
int h,mn; double s;
HMSfrHdec(hd,&h,&mn,&s);
return ToStringHMS(h,mn,s);
}
void EqtoGal(double mjd,double ra,double dec, double *glng,double *glat)
// Coordonnees equatoriales -> Coordonnees galactiques
{
ra *= PI/12.; // radians
dec *= PI/180.; // radians
eq_gal(mjd,ra,dec,glat,glng);
// Vraiment bizarre, sur Linux-g++ glng>=360 ne comprend pas glng==360 ! (CMV)
*glng *= 180./PI; if(*glng>360.) *glng -= 360.; if(*glng==360.) *glng = 0.;
*glat *= 180./PI;
}
void GaltoEq(double mjd,double glng,double glat,double *ra,double *dec)
// Coordonnees galactiques -> Coordonnees equatoriales
{
glng *= PI/180.; // radians
glat *= PI/180.; // radians
gal_eq (mjd,glat,glng,ra,dec);
*ra *= 12./PI; if(*ra>24.) *ra -= 24.; if(*ra==24.) *ra = 0.;
*dec *= 180./PI;
}
void EqtoHor(double geolat,double ha,double dec,double *az,double *alt)
// Coordonnees equatoriales -> Coordonnees horizontales
{
geolat *= PI/180.;
ha *= PI/12.; // radians
dec *= PI/180.; // radians
hadec_aa (geolat,ha,dec,alt,az);
*alt *= 180./PI;
*az *= 180./PI; if(*az>360.) *az -= 360.; if(*az==360.) *az = 0.;
}
void HortoEq(double geolat,double az,double alt,double *ha,double *dec)
// Coordonnees horizontales -> Coordonnees equatoriales
{
geolat *= PI/180.;
alt *= PI/180.; // radians
az *= PI/180.; // radians
aa_hadec (geolat,alt,az,ha,dec);
*ha *= 12./PI;
if(*ha==-12.) *ha = 12.; if(*ha<-12.) *ha += 24.; if(*ha>12.) *ha -= 24.;
*dec *= 180./PI;
}
// Attention, j'ai modifie eq_ecl.c pour proteger NaN
// dans ecleq_aux :
// *q = (sy*ceps)-(cy*seps*sx*sw);
// if(*q<-1.) *q = -PI/2.; else if(*q>1.) *q = PI/2.; else *q = asin(*q);
void EqtoEcl(double mjd,double ra,double dec,double *eclng,double *eclat)
// Coordonnees equatoriales -> Coordonnees ecliptiques
{
ra *= PI/12.; // radians
dec *= PI/180.; // radians
eq_ecl(mjd,ra,dec,eclat,eclng);
*eclng *= 180./PI; if(*eclng>360.) *eclng -= 360.; if(*eclng==360.) *eclng = 0.;
*eclat *= 180./PI;
}
void EcltoEq(double mjd,double eclng,double eclat,double *ra,double *dec)
// Coordonnees ecliptiques -> Coordonnees equatoriales
{
eclat *= PI/180.; // radians
eclng *= PI/180.; // radians
ecl_eq(mjd,eclat,eclng,ra,dec);
*ra *= 12./PI; if(*ra>24.) *ra -= 24.; if(*ra==24.) *ra = 0.;
*dec *= 180./PI;
}
/* given the modified JD, mjd, return the true geocentric ecliptic longitude
* of the sun for the mean equinox of the date, *lsn, in radians, the
* sun-earth distance, *rsn, in AU, and the latitude *bsn, in radians
* (since this is always <= 1.2 arcseconds, in can be neglected by
* calling with bsn = NULL). */
void SunPos(double mjd,double *eclsn,double *ecbsn)
{
double rsn;
sunpos(mjd,eclsn,&rsn,ecbsn);
*eclsn *= 180./PI; if(*eclsn>360.) *eclsn -= 360.; if(*eclsn==360.) *eclsn = 0.;
*ecbsn *= 180./PI;
}
/* given the mjd, find the geocentric ecliptic longitude, lam, and latitude,
* bet, and geocentric distance, rho in a.u. for the moon. also return
* the sun's mean anomaly, *msp, and the moon's mean anomaly, *mdp.
* (for the mean equinox) */
void MoonPos(double mjd,double *eclmn,double *ecbmn)
{
double rho,msp,mdp;
moon(mjd,eclmn,ecbmn,&rho,&msp,&mdp);
*eclmn *= 180./PI; if(*eclmn>360.) *eclmn -= 360.; if(*eclmn==360.) *eclmn = 0.;
*ecbmn *= 180./PI;
}
void PlanetPos(double mjd,int numplan,double *ecl,double *ecb,double *diamang)
/* given a modified Julian date, mjd, and a planet, p, find:
* lpd0: heliocentric longitude,
* psi0: heliocentric latitude,
* rp0: distance from the sun to the planet,
* rho0: distance from the Earth to the planet,
* none corrected for light time, ie, they are the true values for the
* given instant.
* lam: geocentric ecliptic longitude,
* bet: geocentric ecliptic latitude,
* each corrected for light time, ie, they are the apparent values as
* seen from the center of the Earth for the given instant.
* dia: angular diameter in arcsec at 1 AU,
* mag: visual magnitude when 1 AU from sun and earth at 0 phase angle.
* (for the mean equinox) */
{
double lpd0,psi0,rp0,rho0,mag;
plans(mjd,numplan,&lpd0,&psi0,&rp0,&rho0,ecl,ecb,diamang,&mag);
*ecl *= 180./PI; if(*ecl>360.) *ecl -= 360.; if(*ecl==360.) *ecl = 0.;
*ecb *= 180./PI;
}
void JupiterPos(double mjd,double *ecl,double *ecb,double *diamang)
{
PlanetPos(mjd,JUPITER,ecl,ecb,diamang);
}
void SaturnPos(double mjd,double *ecl,double *ecb,double *diamang)
{
PlanetPos(mjd,SATURN,ecl,ecb,diamang);
}
/* Given a coordinate type "typ", convert to standard for astropack */
int CoordConvertToStd(TypAstroCoord typ,double& coord1,double& coord2)
// Return : 0 = OK
// 1 = Type de coordonnees non connu
// 2 = Mauvais range pour coord1
// 4 = Mauvais range pour coord2
// 6 = Mauvais range pour coord1 et coord2
{
int rc = 0;
// ---- Equatoriales alpha,delta
// - standard = [0,24[ , [-90,90]
if(typ&TypCoordEq) {
if(typ&TypCoordDD) {
coord1 = coord1 / 180. * 12.;
} else if(typ&TypCoordRR) {
coord1 = coord1 / PI * 12.;
coord2 = coord2 / PI * 180.;
}
if(coord1==24.) coord1 = 0.;
if(coord1<0. || coord1>=24.) rc+= 2;
if(coord2<-90. || coord2>90. ) rc+= 4;
// ---- Galactiques gLong, gLat
// ---- Horizontales azimuth,altitude
// ---- Ecliptiques EclLong,EclLat
// - standard = [0,360[ , [-90,90]
} else if( typ&TypCoordGal || typ&TypCoordHor || typ&TypCoordEcl) {
if(typ&TypCoordHD) {
coord1 = coord1 / 12. * 180.;
} else if(typ&TypCoordRR) {
coord1 = coord1 / PI * 180.;
coord2 = coord2 / PI * 180.;
}
if(coord1==360.) coord1 = 0.;
if(coord1<0. || coord1>=360.) rc+= 2;
if(coord2<-90. || coord2>90. ) rc+= 4;
} else { // Coordonnees non-connues
rc= 1;
}
return rc;
}