Context Navigation

xastropack.cc@ 3309

Visit:

Last change on this file since 3309 was 2615, checked in by cmv, 21 years ago
using namespace sophya enleve de machdefs.h, nouveau sopnamsp.h cmv 10/09/2004
File size: 27.5 KB

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

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

Context Navigation

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

Download in other formats: