Context Navigation

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

Visit:

Last change on this file since 2551 was 1808, checked in by cmv, 24 years ago
nouvelle interface de gestion des systemes de coordonnees et des unites cmv 3/12/01
File size: 27.5 KB

Rev	Line
[1456]	1	#include <math.h>
	2	#include <stdio.h>
	3	#include "xastropack.h"
	4
[1682]	5	// BUGS BUGS BUGS BUGS BUGS BUGS BUGS BUGS BUGS BUGS
[1791]	6	// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	7	// >>>> Corrections de divers bugs trouve dans libastro (CMV)
[1682]	8	// 1******* In the file vsop87.c line 154:
	9	// p = q/(t_abs[alpha] + alpha * t_abs[alpha-1] * 1e-4 + 1e-35);
	10	// - to avoid t_abs[-1] when alpha=0, replaced by :
	11	// if(alpha>0) p = t_abs[alpha-1]; else p=0;
	12	// p = q/(t_abs[alpha] + alpha * p * 1e-4 + 1e-35);
	13	// Mail envoye a ecdowney@ClearSkyInstitute.com
	14	// 2******* In the file eq_ecl.c line 69:
	15	// q = asin((syceps)-(cysepssx*sw));
	16	// eq_ecl.c Protection NaN dans ecleq_aux, replaced by :
	17	// q = (syceps)-(cysepssx*sw);
	18	// if(q<-1.) q = -PI/2.; else if(q>1.) q = PI/2.; else q = asin(q);
	19	// Mail envoye a ecdowney@ClearSkyInstitute.com
[1791]	20	// >>>> Corrections effectuees dans la version Xephem 3.5
	21	// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
[1682]	22
	23
[1628]	24	/*!
	25	\defgroup XAstroPack XAstroPack module
	26	This module contains simple programs to perform various
	27	astronomical computation (based on the libastro of Xephem).
	28
	29	\verbatim
[1456]	30	// TEMPS: modified Julian date (mjd) (number of days elapsed since 1900 jan 0.5)
	31	// jour [1,31] (dy)
	32	// mois [1,12] (mn)
	33	// annee (yr)
[1515]	34	// universal time [0,24[ (utc)
	35	// Greenwich mean siderial [0,24[ (gst)
	36	// Greenwich mean siderial at 0h UT [0,24[ (gst0)
[1456]	37	// EQUATORIALE: ascension droite en heures [0,24[ (ra)
	38	// declinaison en degres [-90,90] (dec)
[1682]	39	// angle horaire en heures ]-12,12] (-12=12) (ha)
[1678]	40	// temps sideral du lieu: tsid=ha+ra (ou lst)
[1456]	41	// GALACTIQUE: longitude en degres [0,360[ (glng)
	42	// latitude en degres [-90,90] (glat)
[1808]	43	// (colatitude en degres [0,180] (gcolat))
[1515]	44	// HORIZONTAL: azimuth en degres [0,360[ (az)
	45	// (angle round to the east from north+)
[1808]	46	// altitude ou elevation en degres [-90,90] (alt)
	47	// (distance zenitale en degres [0,180] (zendist))
[1515]	48	// ECLIPTIQUE: lontitude ecliptique en degres [0,360[ (eclng)
	49	// (angle round counter clockwise from the vernal equinoxe)
	50	// latitude ecliptique en degres [-90,90] (eclat)
[1808]	51	// (colatitude ecliptique en degres [0,180] (eccolat))
[1515]	52	// GEOGRAPHIE: longitude en degres ]-180,180] (geolng)
[1802]	53	// (angle <0 vers l'ouest, >0 vers l'est)
	54	// latitude en degres [-90,90] (north>0 sud<0) (geolat)
[1808]	55	// (colatitude en degres [0,180] (north=0, sud=180) (geocolat))
	56	//
	57	// --- Remarque sur la colatitude:
	58	// La latitude peut etre remplacee par la colatitude
	59	// (ou altitude/elevation par la distance zenitale):
	60	// latitude : [-90,90] avec 0=equateur, 90=pole nord, -90=pole sud
	61	// colatitude : [0,180] avec 0=pole nord, 90=equateur, 180=pole sud
	62	// colatitude = 90.-latitude , latitude = 90.-colatitude
[1628]	63	\endverbatim
	64	*/
[1456]	65
[1628]	66	/*! \ingroup XAstroPack
[1682]	67	\brief Given a coordinate type "typ", convert to standard for astropack.
	68	\verbatim
	69	La routine convertit (in place) les coordonnees "coord1","coord2"
	70	definies par le type "typ" dans les unites standard de ce systeme
	71	de coordonnees.
	72	"typ" code le systeme de coordonnees astronomiques et les unites utilisees
	73
[1808]	74	- Return : 0 = Problem
	75	TypAstroCoord du systeme de coordonnees retourne
[1682]	76
[1808]	77	- Les types de coordonnees (A,B) sont definies
	78	dans le enum TypAstroCoord (unsigned long):
	79	La premiere coordonnee "A" est de type "longitude"
	80	(alpha,longitude,azimuth)
	81	La deuxieme coordonnee "B" est de type "latidude"
	82	(delta,latitude,altitude ou elevation)
	83
	84	*** Definitions des unites des coordonnees
	85	- Coordonnee:
	86	TypCoordH en heure
	87	TypCoordD en degre
	88	TypCoordR en radian
	89	- Coordonnee "A":
	90	TypCoord1H en heure
	91	TypCoord1D en degre
	92	TypCoord1R en radian
	93	- Defaut (pas de bit leve): radians
	94	- Coordonnee "B":
	95	TypCoord2H en heure
	96	TypCoord2D en degre
	97	TypCoord2R en radian
	98	- Defaut (pas de bit leve): radians
	99	*** Definitions des types d'intervalle utilises
	100	- Coordonnee "A":
	101	TypCoord1C type intervalle
	102	[0,24[ ou [0,360[ ou [0,2Pi[
	103	TypCoord1L type intervalle centre
	104	]-12,12] ou ]-180,180] ou ]-Pi,Pi]
	105	- Defaut (pas de bit leve): TypCoord1C
	106	- Coordonnee "B":
	107	TypCoord2C type intervalle (colatitude)
	108	[0,12] ou [0,180] ou [0,Pi]
	109	TypCoord2L type intervalle centre (latitude)
	110	[-6,6] ou [-90,90] ou [-Pi/2,Pi/2]
	111	- Defaut (pas de bit leve): TypCoord2L (latitude)
	112	*** Les systemes de coordonnes astronomiques
	113	TypCoordEq coordonnees Equatoriales alpha,delta
	114	TypCoordGal coordonnees Galactiques gLong, gLat
	115	TypCoordHor coordonnees Horizontales azimuth,altitude
	116	TypCoordEcl coordonnees Ecliptiques EclLong,EclLat
	117	(Pas de defaut)
	118	*** Les systemes de coordonnes astronomiques "standard"
	119	TypCoordEqStd alpha en heure=[0,24[ delta en degre=[-90,90]
	120	TypCoordGalStd long en degre=[0,360[ lat en degre=[-90,90] (latitude)
	121	TypCoordHorStd long en degre=[0,360[ lat en degre=[-90,90] (latitude)
	122	TypCoordEclStd long en degre=[0,360[ lat en degre=[-90,90] (latitude)
[1682]	123	\endverbatim
	124	*/
[1808]	125	unsigned long CoordConvertToStd(unsigned long typ,double* coord1,double* coord2)
[1682]	126	{
[1808]	127	unsigned long rc = TypCoordUndef;
[1682]	128
	129	if(typ&TypCoordEq) {
[1808]	130	// ---- Equatoriales alpha,delta
	131	if (typ&TypCoord1D) coord1 = deghr(coord1);
	132	else if(typ&TypCoord1R) coord1 = radhr(coord1);
	133	else if(!(typ&TypCoord1H)) coord1 = radhr(coord1);
[1682]	134
[1808]	135	InRange(coord1,24.);
[1682]	136
[1808]	137	if (typ&TypCoord2H) coord2 = hrdeg(coord2);
	138	else if(typ&TypCoord2R) coord2 = raddeg(coord2);
	139	else if(!(typ&TypCoord2D)) coord2 = raddeg(coord2);
	140
	141	if(typ&TypCoord2C) {
	142	InRangeCoLat(coord2,TypUniteD);
	143	ToCoLat(coord2,TypUniteD);
	144	} else InRangeLat(coord2,TypUniteD);
	145
	146	rc=TypCoordEqStd;
	147
	148	} else if(typ&TypCoordGal \|\| typ&TypCoordHor \|\| typ&TypCoordEcl) {
	149	// ---- Galactiques gLong, gLat
	150	// ---- Horizontales azimuth,altitude ou elevation
	151	// ---- Ecliptiques EclLong,EclLat
	152	if (typ&TypCoord1H) coord1 = hrdeg(coord1);
	153	else if(typ&TypCoord1R) coord1 = raddeg(coord1);
	154	else if(!(typ&TypCoord1D)) coord1 = raddeg(coord1);
	155
	156	InRange(coord1,360.);
	157
	158	if (typ&TypCoord2H) coord2 = hrdeg(coord2);
	159	else if(typ&TypCoord2R) coord2 = raddeg(coord2);
	160	else if(!(typ&TypCoord2D)) coord2 = raddeg(coord2);
	161
	162	if(typ&TypCoord2C) {
	163	InRangeCoLat(coord2,TypUniteD);
	164	ToCoLat(coord2,TypUniteD);
	165	} else InRangeLat(coord2,TypUniteD);
	166
	167	if (typ&TypCoordGal) rc=TypCoordGalStd;
	168	else if(typ&TypCoordHor) rc=TypCoordHorStd;
	169	else if(typ&TypCoordEcl) rc=TypCoordEclStd;
	170
	171	} else { // Systeme de coordonnees non-connu
	172	rc=TypCoordUndef;
[1682]	173	}
	174
	175	return rc;
	176	}
	177
	178	/*! \ingroup XAstroPack
[1808]	179	\brief Retourne te type d'unite pour la coordonnee "coordnum"
	180	pour un TypAstroCoord valant "typ"
	181	\verbatim
	182	coordnum : numero de coordonnee: 1 ou 2
	183	retourne: TypUniteH si la coordonnee est en heure
	184	TypUniteD si la coordonnee est en degre
	185	TypUniteR si la coordonnee est en radian
	186	TypUniteR par defaut
	187	TypCoordUndef si le numero de coordonnee est errone.
	188	\endverbatim
	189	*/
	190	unsigned long GetCoordUnit(int coordnum,unsigned long typ)
	191	{
	192	if(coordnum==1) {
	193	if (typ&TypCoord1H) return TypUniteH;
	194	else if(typ&TypCoord1D) return TypUniteD;
	195	else if(typ&TypCoord1R) return TypUniteR;
	196	else return TypUniteR;
	197	}
	198	if(coordnum==2) {
	199	if (typ&TypCoord2H) return TypUniteH;
	200	else if(typ&TypCoord2D) return TypUniteD;
	201	else if(typ&TypCoord2R) return TypUniteR;
	202	else return TypUniteR;
	203	}
	204	return TypCoordUndef;
	205	}
	206
	207	/*! \ingroup XAstroPack
	208	\brief Pour decoder et transcrire en TypAstroCoord une chaine
	209	donnant la structure du systeme de coordonnees.
	210	\verbatim
	211	ctype = "CAaBb"
	212	C: type de coordonnees: E Equatoriales
	213	G Galactiques
	214	H Horizontales
	215	S Ecliptiques
	216	pas de defaut
	217	A: unite de la coordonnee 1 (alpha,longitude etc...)
	218	H heure
	219	D degre
	220	R radian
	221	defaut radian
	222	a: type d'intervalle pour la coordonnee 1
	223	C intervalle [0,24[ [0,360[ [0,2*Pi[
	224	L intervalle [-12,12[ [-180,180[ [-Pi,Pi[
	225	(defaut: C)
	226	A: unite de la coordonnee 2 (delta,latitude etc...)
	227	H heure
	228	D degre
	229	R radian
	230	defaut radian
	231	a: type d'intervalle pour la coordonnee 2
	232	C intervalle [0,12] [0,180] [0,Pi]
	233	(type colatitude)
	234	L intervalle [-6,6] [-90,90][ [-Pi/2,Pi/2]
	235	(defaut: L)
	236
	237	Exemple: GDCDL : galactiques long=[0,360[ lat=[-90,90]
	238	GDxDx ou GDxD: idem
	239	Gxxxx ou G : galactiques long=[0,2*Pi[ lat=[-Pi/2,Pi/2]
	240	Exemple: EHCDL : equatoriales alpha=[0,24[ delta=[-90,90]
	241	EHxDx ou EHxD : idem
	242	Exxxx ou E : equatoriales alpha=[0,2*Pi[ delta=[-Pi/2,Pi/2]
	243
	244	- Retourne 0 si probleme dans la chaine
	245	\endverbatim
	246	*/
	247	unsigned long DecodeTypAstro(const char *ctype)
	248	{
	249	if(ctype==NULL) return TypCoordUndef;
	250	int len = strlen(ctype);
	251	if(len<1) return TypCoordUndef;
	252
	253	unsigned long typ=0;
	254	// Le type de systeme de coordonnees
	255	int i = 0;
	256	if (ctype[i]=='e' \|\| ctype[i]=='E') typ=TypCoordEq;
	257	else if(ctype[i]=='g' \|\| ctype[i]=='G') typ=TypCoordGal;
	258	else if(ctype[i]=='h' \|\| ctype[i]=='H') typ=TypCoordHor;
	259	else if(ctype[i]=='s' \|\| ctype[i]=='S') typ=TypCoordEcl;
	260	else return TypCoordUndef;
	261	// La coordonnee 1
	262	i = 1;
	263	if(i>=len)
	264	{typ \|= TypCoord1R\|TypCoord1C\|TypCoord2R\|TypCoord2L; return typ;}
	265	if (ctype[i]=='h' \|\| ctype[i]=='H') typ \|= TypCoord1H;
	266	else if(ctype[i]=='d' \|\| ctype[i]=='D') typ \|= TypCoord1D;
	267	else if(ctype[i]=='r' \|\| ctype[i]=='R') typ \|= TypCoord1R;
	268	else typ \|= TypCoord1R;
	269	i = 2;
	270	if(i>=len) {typ \|= TypCoord1C\|TypCoord2R\|TypCoord2L; return typ;}
	271	if (ctype[i]=='c' \|\| ctype[i]=='C') typ \|= TypCoord1C;
	272	else if(ctype[i]=='l' \|\| ctype[i]=='L') typ \|= TypCoord1L;
	273	else typ \|= TypCoord1C;
	274	// La coordonnee 2
	275	i = 3;
	276	if(i>=len) {typ \|= TypCoord2R\|TypCoord2L; return typ;}
	277	if (ctype[i]=='h' \|\| ctype[i]=='H') typ \|= TypCoord2H;
	278	else if(ctype[i]=='d' \|\| ctype[i]=='D') typ \|= TypCoord2D;
	279	else if(ctype[i]=='r' \|\| ctype[i]=='R') typ \|= TypCoord2R;
	280	else typ \|= TypCoord2R;
	281	i = 4;
	282	if(i>=len) {typ \|= TypCoord2L; return typ;}
	283	if (ctype[i]=='c' \|\| ctype[i]=='C') typ \|= TypCoord2C;
	284	else if(ctype[i]=='l' \|\| ctype[i]=='L') typ \|= TypCoord2L;
	285	else typ \|= TypCoord2L;
	286	// Return
	287	return typ;
	288	}
	289
	290	/*! \ingroup XAstroPack
	291	\brief Idem DecodeTypAstro(char *) mais a l'envers
	292	*/
	293	string DecodeTypAstro(unsigned long typ)
	294	{
	295	string s = "";
	296
	297	if (typ&TypCoordEq) s += "E";
	298	else if(typ&TypCoordGal) s += "G";
	299	else if(typ&TypCoordHor) s += "H";
	300	else if(typ&TypCoordEcl) s += "S";
	301	else s += "x";
	302
	303	if (typ&TypCoord1H) s += "H";
	304	else if(typ&TypCoord1D) s += "D";
	305	else if(typ&TypCoord1R) s += "R";
	306	else s += "x";
	307
	308	if (typ&TypCoord1C) s += "C";
	309	else if(typ&TypCoord1L) s += "L";
	310	else s += "x";
	311
	312	if (typ&TypCoord2H) s += "H";
	313	else if(typ&TypCoord2D) s += "D";
	314	else if(typ&TypCoord2R) s += "R";
	315	else s += "x";
	316
	317	if (typ&TypCoord2C) s += "C";
	318	else if(typ&TypCoord2L) s += "L";
	319	else s += "x";
	320
	321	return s;
	322	}
	323
	324	/*! \ingroup XAstroPack
	325	\brief Pour convertir la latitude en colatitude et vice-versa (in place)
	326	\verbatim
	327	val = valeur a convertir qui doit etre:
	328	si type "latitude" dans [-6,6] [-90,90] [-Pi/2,Pi/2]
	329	si type "colatitude" dans [0,12] [0,180] [0,Pi]
	330	typ = type d'unite: heure TypUniteH
	331	degre TypUniteD
	332	radian TypUniteR
	333	(Defaut: radian TypUniteR)
	334	\endverbatim
	335	*/
	336	void ToCoLat(double *val,unsigned long typ)
	337	{
	338	if (typ&TypUniteH) val = 6. - val;
	339	else if(typ&TypUniteD) val = 90. - val;
	340	else if(typ&TypUniteR) val = PI/2. - val;
	341	else val = PI/2. - val;
	342	}
	343
	344	/*! \ingroup XAstroPack
	345	\brief Pour remettre la valeur de la COLATITUDE "val" dans la dynamique [0.,range]
	346	\verbatim
	347	val = valeur a convertir qui doit etre mise dans [0,12] [0,180] [0,Pi]
	348	typ = type d'unite: heure TypUniteH
	349	degre TypUniteD
	350	radian TypUniteR
	351	ex en degre: 0 -> 0 , 90 -> 90 , 180 -> 180 , 270 -> 90 , 360 -> 0
	352	-90 -> 90 , -180 -> 180 , -270 -> 90 , -360 -> 0
	353	\endverbatim
	354	*/
	355	void InRangeCoLat(double *val,unsigned long typ)
	356	{
	357	double range=PI;
	358	if(typ==TypUniteH) range=12.; else if(typ==TypUniteD) range=180.;
	359	InRange(val,2.*range,range);
	360	if(val<0.) val*=-1.;
	361	}
	362
	363	/*! \ingroup XAstroPack
	364	\brief Pour remettre la valeur de la LATITUDE "val" dans la dynamique [-range/2,range/2]
	365	\verbatim
	366	val = valeur a convertir qui doit etre mise dans [-6,6] [-90,90] [-Pi/2,Pi/2]
	367	typ = type d'unite: heure TypUniteH
	368	degre TypUniteD
	369	radian TypUniteR
	370	ex en degre: 0 -> 0 , 90 -> 90 , 180 -> 0 , 270 -> -90 , 360 -> 0
	371	-90 -> -90 , -180 -> 0 , -270 -> 90 , -360 -> 0
	372	\endverbatim
	373	*/
	374	void InRangeLat(double *val,unsigned long typ)
	375	{
	376	double range = PI;
	377	if(typ==TypUniteH) range = 12.; else if(typ==TypUniteD) range = 180.;
	378	InRange(val,2.*range,range);
	379	if(val>range/2.) val = range - *val;
	380	else if(val<-range/2.) val = -(range + *val);
	381	}
	382
	383	/*! \ingroup XAstroPack
[1682]	384	\brief Compute MJD from date
	385	\verbatim
	386	MJD = modified Julian date (number of days elapsed since 1900 jan 0.5),
[1688]	387	dy is the decimale value of the day: dy = int(dy) + utc/24.
[1682]	388	\endverbatim
	389	*/
	390	double MJDfrDate(double dy,int mn,int yr)
	391	{
	392	double mjd;
	393	cal_mjd(mn,dy,yr,&mjd);
	394	return mjd;
	395	}
	396
	397	/*! \ingroup XAstroPack
	398	\brief Compute date from MJD
	399	*/
	400	void DatefrMJD(double mjd,double dy,int mn,int *yr)
	401	{
	402	mjd_cal(mjd,mn,dy,yr);
	403	}
	404
	405	/*! \ingroup XAstroPack
	406	\brief Given a mjd, return the year as a double.
	407	*/
	408	double YearfrMJD(double mjd)
	409	{
	410	double yr;
	411	mjd_year(mjd,&yr);
	412	return yr;
	413	}
	414
	415	/*! \ingroup XAstroPack
	416	\brief Given a decimal year, return mjd
	417	*/
	418	double MJDfrYear(double yr)
	419	{
	420	double mjd;
	421	year_mjd(yr,&mjd);
	422	return mjd;
	423	}
	424
	425	/*! \ingroup XAstroPack
	426	\brief Given a mjd, return the year and number of days since 00:00 Jan 1
	427	\warning: if mjd = 2 January -> number of days = 1
	428	*/
	429	void YDfrMJD(double mjd,double dy,int yr)
	430	{
	431	mjd_dayno(mjd,yr,dy);
	432	}
	433
	434	/*! \ingroup XAstroPack
	435	\brief Given a year,
	436	*/
	437	int IsLeapYear(int y)
	438	{
	439	return isleapyear(y);
	440	}
	441
	442	/*! \ingroup XAstroPack
	443	\brief given an mjd, set *dow to 0..6 according to which day of the week it falls on (0=sunday).
	444	\return return 0 if ok else -1 if can't figure it out.
	445	*/
	446	int DayOrder(double mjd,int *dow)
	447	{
	448	return mjd_dow(mjd,dow);
	449	}
	450
	451	/*! \ingroup XAstroPack
	452	\brief given a mjd, return the the number of days in the month.
	453	*/
	454	int DaysInMonth(double mjd)
	455	{
	456	int ndays;
	457	mjd_dpm(mjd,&ndays);
	458	return ndays;
	459	}
	460
	461	/*! \ingroup XAstroPack
	462	\brief Given a mjd, truncate it to the beginning of the whole day
	463	*/
	464	double MJDat0hFrMJD(double mjd)
	465	{
	466	return mjd_day(mjd);
	467	}
	468
	469	/*! \ingroup XAstroPack
	470	\brief Given a mjd, return the number of hours past midnight of the whole day
	471	*/
	472	double HfrMJD(double mjd)
	473	{
	474	return mjd_hr(mjd);
	475	}
	476
	477	/*! \ingroup XAstroPack
	478	\brief Give GST from UTC
	479	\verbatim
	480	Given a modified julian date, mjd, and a universally coordinated time, utc,
	481	return greenwich mean siderial time, *gst.
	482	N.B. mjd must be at the beginning of the day.
	483	\endverbatim
	484	*/
	485	double GSTfrUTC(double mjd0,double utc)
	486	{
	487	double gst;
	488	utc_gst(mjd0,utc,&gst);
	489	return gst;
	490	}
	491
	492	/*! \ingroup XAstroPack
	493	\brief Give UTC from GST
	494	\verbatim
	495	Given a modified julian date, mjd, and a greenwich mean siderial time, gst,
	496	return universally coordinated time, *utc.
	497	N.B. mjd must be at the beginning of the day.
	498	\endverbatim
	499	*/
	500	double UTCfrGST(double mjd0,double gst)
	501	{
	502	double utc;
	503	gst_utc(mjd0,gst,&utc);
	504	return utc;
	505	}
	506
	507	/*! \ingroup XAstroPack
	508	\brief Given apparent altitude find airmass.
	509	*/
	510	double AirmassfrAlt(double alt)
	511	{
	512	double x;
	513	alt = degrad(alt);
	514	airmass(alt,&x);
	515	return x;
	516	}
	517
	518	/*! \ingroup XAstroPack
	519	\brief given geocentric time "jd" and coords of a distant object at "ra/dec" (J2000),
	520	find the difference "hcp" in time between light arriving at earth vs the sun.
	521	\return "hcp" must be subtracted from "geocentric jd" to get "heliocentric jd".
	522	\warning "jd" is the TRUE Julian day (jd = mjd+MJD0).
	523	*/
	524	double HelioCorr(double jd,double ra,double dec)
	525	{
	526	double hcp;
	527	ra=hrrad(ra);
	528	dec=degrad(dec);
	529	heliocorr(jd,ra,dec,&hcp);
	530	return hcp;
	531	}
	532
	533	/*! \ingroup XAstroPack
[1628]	534	\brief gmst0() - return Greenwich Mean Sidereal Time at 0h UT
[1679]	535	\param mjd0 = date at 0h UT in julian days since MJD0
[1628]	536	*/
[1456]	537	double GST0(double mjd0)
[1678]	538	/* Copie depuis le code de Xephem (utc_gst.c) car pas prototype*/
[1456]	539	{
	540	double T, x;
	541	T = ((int)(mjd0 - 0.5) + 0.5 - J2000)/36525.0;
	542	x = 24110.54841 +
	543	(8640184.812866 + (0.093104 - 6.2e-6 * T) * T) * T;
	544	x /= 3600.0;
	545	range(&x, 24.0);
	546	return (x);
	547	}
	548
[1628]	549	/*! \ingroup XAstroPack
[1678]	550	\brief return local sidereal time from modified julian day and longitude
	551	\warning nutation or obliquity correction are taken into account.
	552	*/
	553	double LSTfrMJD(double mjd,double geolng)
	554	{
	555	double eps,lst,deps,dpsi;
	556	utc_gst(mjd_day(mjd),mjd_hr(mjd),&lst);
[1679]	557	lst += deghr(geolng);
[1678]	558	obliquity(mjd,&eps);
	559	nutation(mjd,&deps,&dpsi);
[1679]	560	lst += radhr(dpsi*cos(eps+deps));
	561	InRange(&lst,24.);
[1678]	562	return lst;
	563	}
	564
	565	/*! \ingroup XAstroPack
[1628]	566	\brief Give a time in h:mn:s from a decimal hour
	567	\verbatim
[1456]	568	// INPUT: hd
[1465]	569	// OUTPUT: h mn s (h,mn,s >=< 0)
	570	// REMARQUE: si hd<0 alors h<0 ET mn<0 ET s<0
	571	// EX: 12.51 -> h=12 mn=30 s=10 ;
	572	// -12.51 -> h=-12 mn=-30 s=-10 ;
[1628]	573	\endverbatim
	574	*/
	575	void HMSfrHdec(double hd,int h,int mn,double *s)
[1456]	576	{
	577	int sgn=1;
	578	if(hd<0.) {sgn=-1; hd*=-1.;}
	579	*h = int(hd);
	580	mn = int((hd-(double)(h))*60.);
	581	s = (hd - (double)(h) - (double)(mn)/60.)3600.;
	582	// pb precision
	583	if(s<0.) s = 0.;
	584	if(s>60. \|\| s==60.) {s-=60.; mn+=1;} // s=double attention comparaison
	585	if(mn<0) mn = 0;
	586	if(mn>=60) {mn-=60; *h+=1;}
[1465]	587	h = sgn; mn = sgn; s = (double)sgn;
[1456]	588	}
	589
[1628]	590	/*! \ingroup XAstroPack
	591	\brief Give a decimal hour from a time in h:mn:s
	592	\verbatim
[1465]	593	// INPUT: h , mn , s (h,mn,s >=< 0)
	594	// RETURN: en heures decimales
	595	// REMARQUE: pour avoir hd=-12.51 <- h=-12 mn=-30 s=-10
[1628]	596	\endverbatim
	597	*/
	598	double HdecfrHMS(int h,int mn,double s)
[1456]	599	{
[1465]	600	return ((double)h + (double)mn/60. + s/3600.);
[1456]	601	}
	602
[1628]	603	/*! \ingroup XAstroPack
	604	\brief Give a time string from a time in h:mn:s
	605	\verbatim
[1465]	606	// INPUT: h , mn , s (h,mn,s >=< 0)
[1456]	607	// RETURN: string h:mn:s
[1628]	608	\endverbatim
	609	*/
	610	string ToStringHMS(int h,int mn,double s)
[1456]	611	{
[1465]	612	double hd = HdecfrHMS(h,mn,s); // put in range
	613	HMSfrHdec(hd,&h,&mn,&s);
[1456]	614	char str[128];
[1465]	615	if(hd<0.)
	616	sprintf(str,"-%d:%d:%.3f",-h,-mn,-s);
	617	else
	618	sprintf(str,"%d:%d:%.3f",h,mn,s);
[1456]	619	string dum = str;
	620	return dum;
	621	}
	622
[1628]	623	/*! \ingroup XAstroPack
	624	\brief Give a time string from a decimal hour
	625	*/
[1456]	626	string ToStringHdec(double hd)
	627	{
	628	int h,mn; double s;
[1465]	629	HMSfrHdec(hd,&h,&mn,&s);
[1456]	630	return ToStringHMS(h,mn,s);
	631	}
	632
[1628]	633	/*! \ingroup XAstroPack
[1679]	634	\brief Compute precession between 2 dates.
	635	*/
	636	void Precess(double mjd1,double mjd2,double ra1,double dec1,double ra2,double dec2)
	637	{
	638	ra1 = hrrad(ra1); // radians
	639	dec1 = degrad(dec1); // radians
	640	precess(mjd1,mjd2,&ra1,&dec1);
	641	*ra2 = radhr(ra1); InRange(ra2,24.);
	642	*dec2 = raddeg(dec1);
	643	}
	644
	645	/*! \ingroup XAstroPack
[1678]	646	\brief Convert equatorial coordinates for the given epoch into galactic coordinates
[1628]	647	*/
[1456]	648	void EqtoGal(double mjd,double ra,double dec, double glng,double glat)
	649	// Coordonnees equatoriales -> Coordonnees galactiques
	650	{
[1679]	651	ra = hrrad(ra); // radians
	652	dec = degrad(dec); // radians
[1456]	653	eq_gal(mjd,ra,dec,glat,glng);
	654	// Vraiment bizarre, sur Linux-g++ glng>=360 ne comprend pas glng==360 ! (CMV)
[1679]	655	glng = raddeg(glng); InRange(glng,360.);
	656	glat = raddeg(glat);
[1456]	657	}
	658
[1628]	659	/*! \ingroup XAstroPack
[1678]	660	\brief Convert galactic coordinates into equatorial coordinates at the given epoch
[1628]	661	*/
[1456]	662	void GaltoEq(double mjd,double glng,double glat,double ra,double dec)
	663	// Coordonnees galactiques -> Coordonnees equatoriales
	664	{
[1679]	665	glng = degrad(glng); // radians
	666	glat = degrad(glat); // radians
[1456]	667	gal_eq (mjd,glat,glng,ra,dec);
[1679]	668	ra = radhr(ra); InRange(ra,24.);
	669	dec = raddeg(dec);
[1456]	670	}
	671
[1628]	672	/*! \ingroup XAstroPack
[1678]	673	\brief Convert equatorial coordinates (with hour angle instead of right ascension) into horizontal coordinates.
[1628]	674	*/
[1678]	675	void EqHtoHor(double geolat,double ha,double dec,double az,double alt)
[1456]	676	// Coordonnees equatoriales -> Coordonnees horizontales
	677	{
[1679]	678	geolat = degrad(geolat); // radians
	679	ha = hrrad(ha); // radians
	680	dec = degrad(dec); // radians
[1456]	681	hadec_aa (geolat,ha,dec,alt,az);
[1679]	682	alt = raddeg(alt);
	683	az = raddeg(az); InRange(az,360.);
[1456]	684	}
	685
[1628]	686	/*! \ingroup XAstroPack
[1678]	687	Convert horizontal coordinates into equatorial coordinates (with hour angle instead of right ascension).
[1628]	688	*/
[1678]	689	void HortoEqH(double geolat,double az,double alt,double ha,double dec)
[1456]	690	// Coordonnees horizontales -> Coordonnees equatoriales
	691	{
[1679]	692	geolat = degrad(geolat); // radians
	693	alt = degrad(alt); // radians
	694	az = degrad(az); // radians
[1456]	695	aa_hadec (geolat,alt,az,ha,dec);
[1679]	696	ha = radhr(ha); InRange(ha,24.,12.);
	697	dec = raddeg(dec);
[1456]	698	}
	699
[1628]	700	/*! \ingroup XAstroPack
[1678]	701	\brief Convert equatorial coordinates into horizontal coordinates.
[1628]	702	*/
[1678]	703	void EqtoHor(double geolat,double lst,double ra,double dec,double az,double alt)
	704	// Coordonnees equatoriales -> Coordonnees horizontales
	705	{
[1679]	706	double ha = lst - ra; InRange(&ha,24.,12.);
	707	geolat = degrad(geolat); // radians
	708	ha = hrrad(ha); // radians
	709	dec = degrad(dec); // radians
[1678]	710	hadec_aa (geolat,ha,dec,alt,az);
[1679]	711	alt = raddeg(alt);
	712	az = raddeg(az); InRange(az,360.);
[1678]	713	}
	714
	715	/*! \ingroup XAstroPack
	716	Convert horizontal coordinates into equatorial coordinates.
	717	*/
	718	void HortoEq(double geolat,double lst,double az,double alt,double ra,double dec)
	719	// Coordonnees horizontales -> Coordonnees equatoriales
	720	{
	721	double ha;
[1679]	722	geolat = degrad(geolat); // radians
	723	alt = degrad(alt); // radians
	724	az = degrad(az); // radians
[1678]	725	aa_hadec (geolat,alt,az,&ha,dec);
[1679]	726	*ra = lst - radhr(ha); InRange(ra,24.);
	727	dec = raddeg(dec);
[1678]	728	}
	729
	730	/*! \ingroup XAstroPack
	731	\brief Convert equatorial coordinates into geocentric ecliptic coordinates given the modified Julian date.
	732	\warning Correction for the effect on the angle of the obliquity due to nutation is not included.
	733	*/
[1456]	734	void EqtoEcl(double mjd,double ra,double dec,double eclng,double eclat)
	735	// Coordonnees equatoriales -> Coordonnees ecliptiques
	736	{
[1679]	737	ra = hrrad(ra); // radians
	738	dec = degrad(dec); // radians
[1456]	739	eq_ecl(mjd,ra,dec,eclat,eclng);
[1679]	740	eclng = raddeg(eclng); InRange(eclng,360.);
	741	eclat = raddeg(eclat);
[1456]	742	}
	743
[1628]	744	/*! \ingroup XAstroPack
[1678]	745	\brief Convert geocentric ecliptic coordinates into equatorial coordinates given the modified Julian date.
	746	\warning Correction for the effect on the angle of the obliquity due to nutation is not included.
[1628]	747	*/
[1456]	748	void EcltoEq(double mjd,double eclng,double eclat,double ra,double dec)
	749	// Coordonnees ecliptiques -> Coordonnees equatoriales
	750	{
[1679]	751	eclat = degrad(eclat); // radians
	752	eclng = degrad(eclng); // radians
[1456]	753	ecl_eq(mjd,eclat,eclng,ra,dec);
[1679]	754	ra = radhr(ra); InRange(ra,24.);
	755	dec = raddeg(dec);
[1456]	756	}
	757
[1628]	758	/*! \ingroup XAstroPack
	759	\brief Give Sun position
	760	\verbatim
	761	given the modified JD, mjd, return the true geocentric ecliptic longitude
[1679]	762	of the sun for the mean equinox of the date, *eclsn, in degres, the
	763	sun-earth distance, rsn, in AU, and the latitude ecbsn, in degres
[1628]	764	(since this is always <= 1.2 arcseconds, in can be neglected by
[1679]	765	calling with ecbsn = NULL).
	766	- REMARQUE:
	767	* if the APPARENT ecliptic longitude is required, correct the longitude for
	768	* nutation to the true equinox of date and for aberration (light travel time,
	769	* approximately -9.27e7/186000/(360024365)2pi = -9.93e-5 radians).
[1628]	770	\endverbatim
	771	*/
[1679]	772	void SunPos(double mjd,double eclsn,double ecbsn,double *rsn)
[1456]	773	{
[1679]	774	sunpos(mjd,eclsn,rsn,ecbsn);
	775	eclsn = raddeg(eclsn); InRange(eclsn,360.);
	776	if(ecbsn!=NULL) ecbsn = raddeg(ecbsn);
[1456]	777	}
	778
[1628]	779	/*! \ingroup XAstroPack
	780	\brief Give Moon position
	781	\verbatim
	782	given the mjd, find the geocentric ecliptic longitude, lam, and latitude,
	783	bet, and geocentric distance, rho in a.u. for the moon. also return
	784	the sun's mean anomaly, msp, and the moon's mean anomaly, mdp.
	785	(for the mean equinox)
	786	\endverbatim
	787	*/
[1679]	788	void MoonPos(double mjd,double eclmn,double ecbmn,double *rho)
[1456]	789	{
[1679]	790	double msp,mdp;
	791	moon(mjd,eclmn,ecbmn,rho,&msp,&mdp);
	792	eclmn = raddeg(eclmn); InRange(eclmn,360.);
	793	ecbmn = raddeg(ecbmn);
[1456]	794	}
	795
[1628]	796	/*! \ingroup XAstroPack
	797	\brief Give planet position
	798	\verbatim
	799	* given a modified Julian date, mjd, and a planet, p, find:
[1679]	800	* sunecl: heliocentric longitude,
	801	* sunecb: heliocentric latitude,
	802	* sundist: distance from the sun to the planet,
	803	* geodist: distance from the Earth to the planet,
[1456]	804	* none corrected for light time, ie, they are the true values for the
	805	* given instant.
[1679]	806	* geoecl: geocentric ecliptic longitude,
	807	* geoecb: geocentric ecliptic latitude,
[1456]	808	* each corrected for light time, ie, they are the apparent values as
	809	* seen from the center of the Earth for the given instant.
[1679]	810	* diamang: angular diameter in arcsec at 1 AU,
[1456]	811	* mag: visual magnitude when 1 AU from sun and earth at 0 phase angle.
[1628]	812	* (for the mean equinox)
[1679]	813	* all angles are in degres, all distances in AU.
	814	*
	815	* corrections for nutation and abberation must be made by the caller. The RA
	816	* and DEC calculated from the fully-corrected ecliptic coordinates are then
	817	* the apparent geocentric coordinates. Further corrections can be made, if
	818	* required, for atmospheric refraction and geocentric parallax.
[1628]	819	\endverbatim
	820	*/
[1679]	821	void PlanetPos(double mjd,int numplan,double sunecl,double sunecb,double *sundist
	822	,double geodist,double geoecl,double *geoecb
	823	,double diamang,double mag)
[1456]	824	{
[1679]	825	plans(mjd,numplan,sunecl,sunecb,sundist,geodist,geoecl,geoecb,diamang,mag);
	826	geoecl = raddeg(geoecl); InRange(geoecl,360.);
	827	geoecb = raddeg(geoecb);
	828	sunecl = raddeg(sunecl); InRange(sunecl,360.);
	829	sunecb = raddeg(sunecb);
[1456]	830	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: