1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
|
/*****************************************************************************
* Johns Hopkins University
* Center For Astrophysical Sciences
* FUSE
*****************************************************************************
*
* Synopsis: int eclipse(double pos[3], double mjdate)
*
* Description: Determines if FUSE is in the day or night portion of the orbit.
*
*
* Arguments: double pos Cartesian vector xyz coord of satellite.
* double mjdate Julian date of observation
*
* Returns: int 0 if day portion of orbit
* 1 if night portion of orbit
*
* Calls: slaDcc2s dcc2s.f
* slaDranrm dranrm.f
* slaDsep dsep.f
* slaEvp evp.f
* slaPreces preces.f
*
* History: 03/10/98 E. Murphy Begin work.
* 03/10/98 E. Murphy Initial version working
* 07/07/99 E. Murphy Replaced x,y,z in call with pos
* 07/20/99 jm Removed duplicate define statements
* 08/26/99 E. Murphy ang_sep now returned.
* 05/31/00 peb Implemented cfortran.h calls for slalib
* functions.
* 12/18/03 bjg Change calfusettag.h to calfuse.h
* 06/17/04 bjg 1.4 Corrected cfortran call to sla_preces
*
* References: This routine makes use of the Starlink set of astronomical
* subroutines (SLALIB). More information can be found at
* http://star-www.rl.ac.uk.
****************************************************************************/
#include <stdio.h>
#include <math.h>
#ifdef CFORTRAN
#include "cfortran.h"
PROTOCCALLSFSUB3(SLA_DCC2S, sla_dcc2s, DOUBLEV, PDOUBLE, PDOUBLE)
#define slaDcc2s(V, A, B) \
CCALLSFSUB3(SLA_DCC2S, sla_dcc2s, DOUBLEV, PDOUBLE, PDOUBLE, V, A, B)
PROTOCCALLSFFUN1(DOUBLE, SLA_DRANRM, sla_dranrm, DOUBLE)
#define slaDranrm(ANGLE) \
CCALLSFFUN1(SLA_DRANRM, sla_dranrm, DOUBLE, ANGLE)
PROTOCCALLSFFUN4(DOUBLE, SLA_DSEP, sla_dsep, DOUBLE, DOUBLE, DOUBLE, \
DOUBLE)
#define slaDsep(A1, B1, A2, B2) \
CCALLSFFUN4(SLA_DSEP, sla_dsep, DOUBLE, DOUBLE, DOUBLE, DOUBLE, \
A1, B1, A2, B2)
PROTOCCALLSFSUB6(SLA_EVP, sla_evp, DOUBLE, DOUBLE, DOUBLEV, DOUBLEV, \
DOUBLEV, DOUBLEV)
#define slaEvp(DATE, DEQX, DVB, DPB, DVH, DPH) \
CCALLSFSUB6(SLA_EVP, sla_evp, DOUBLE, DOUBLE, DOUBLEV, DOUBLEV, \
DOUBLEV, DOUBLEV, DATE, DEQX, DVB, DPB, DVH, DPH)
PROTOCCALLSFSUB5(SLA_PRECES, sla_preces, STRING, DOUBLE, DOUBLE, PDOUBLE, \
PDOUBLE)
#define slaPreces(SYSTEM, EP0, EP1, RA, DC) \
CCALLSFSUB5(SLA_PRECES, sla_preces, STRING, DOUBLE, DOUBLE, PDOUBLE, \
PDOUBLE, SYSTEM, EP0, EP1, RA, DC)
#else
#include "slalib.h"
#include "slamac.h"
#endif
#include "calfuse.h"
int eclipse(double *pos, double mjdate, double *ang_sep)
{
/* Define variables. */
int i, retcode;
double dvb[3], dpb[3], dvh[3], dph[3], epoch2;
double r, ra_earth, dec_earth, ra_sun, dec_sun, ang_earth, ang_sun;
char fk5_st[10]="FK5";
/* Calculate the J2000.0 apparent RA and DEC of the Earth */
r=sqrt(pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2]);
ra_earth=atan2(-pos[1],-pos[0]);
dec_earth=asin(-pos[2]/r);
/* We must precess the J2000 RA and DEC to date of observation. This
* is a very small correction, and could probably be ignored.
*/
epoch2=2000.0-((51544.0-mjdate)/365.25);
#ifdef CFORTRAN
slaPreces(fk5_st, 2000.0, epoch2, ra_earth, dec_earth);
#else
slaPreces(fk5_st, 2000.0, epoch2, &ra_earth, &dec_earth);
#endif
/* Evp returns four 3-vectors containing the barycentric velocity and
* position (dvb,dpv) and the heliocentric velocity and position
* (dvh,dph) of the Earth on the date mjdate. It requires a modified
* Julian date as input. Output has units of AU for positions and
* AU/s for velocities.
*/
slaEvp(mjdate, 2000.0, dvb, dpb, dvh, dph);
/* Convert the 3-vector position of the Sun into a J2000.0 RA and DEC */
for (i=0; i<3; i++) dph[i]*=-1.0;
#ifdef CFORTRAN
slaDcc2s(dph, ra_sun, dec_sun);
#else
slaDcc2s(dph, &ra_sun, &dec_sun);
#endif
ra_sun=slaDranrm(ra_sun);
/* We must precess the J2000 RA and DEC to date of observation. This
* is a very small correction, and could probably be ignored.
*/
epoch2=2000.0-((51544.0-mjdate)/365.25);
#ifdef CFORTRAN
slaPreces(fk5_st, 2000.0, epoch2, ra_sun, dec_sun);
#else
slaPreces(fk5_st, 2000.0, epoch2, &ra_sun, &dec_sun);
#endif
/* Compute the angular sizes of the Earth and the Sun. */
ang_earth=asin(RE/r);
ang_sun=asin(RS/(sqrt(dph[0]*dph[0]+dph[1]*dph[1]+dph[2]*dph[2])*AU));
/* Compute the angular separation of the Earth and the Sun. */
*ang_sep=slaDsep(ra_earth, dec_earth, ra_sun, dec_sun);
retcode=0;
if (*ang_sep < ang_earth+ang_sun) retcode=1;
/* Convert ang_sep into degrees. */
*ang_sep /= RADIAN;
return retcode;
}
|
