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
149
150
151
152
153
154
155
156
157
158
159
|
/*****************************************************************************
* Johns Hopkins University
* Center For Astrophysical Sciences
* FUSE
*****************************************************************************
*
* Synopsis: double state_limb(double pos[3], double mjdate,
* double ra, double dec, double *zdis,tint day_limb)
*
* Description: Computes the earth limb angle from J2000 RA and DEC given
* the 3-vector position of FUSE.
*
* Arguments: double pos (km) 3-vector position of satellite
* double mjdate (days) Modified Julian date
* double ra,dec (deg) J2000.0 RA and DEC
* double zdist (deg) zenith angle
* integer day_limb output : 1 if bright limb, 0 if dark limb
*
* Returns: double (deg) Earth limb angle
* integer day_limb will have the value 1 or 0 when
* returning in the main program.
*
* Calls: slaPreces preces.f
* slaDcc2s
* slaDranrm
* slaEvp
*
* History: 03/10/98 E. Murphy Begin work.
* 03/10/98 E. Murphy Initial version working
* 07/07/99 E. Murphy Changed call to use pos instead of x,y,z
* 08/26/99 E. Murphy Added zenith distance.
* 05/31/00 peb Implemented cfortran.h calls for slalib
* functions.
* Ake, T. 1998 in The Scientific Impact of the Goddard
* High Resolution Spectrograph, ed. J. C. Brandt et al.,
* ASP Conference Series, in preparation.
*
* 09/23/00 v1.5 ma Now this function determine bright or
* dark limb
* To do so, the procedure has a new
* input argument: day_limb.
* day_limb = 1 if bright, 0 if dark.
* 09/24/00 v1.6 ma Fixed a bug in the input argument
* 09/27/00 v1.7 ma Fixed bug in limbvec calc.
* 10/02/00 v1.8 jwk Add fortran wrappers for slaDcc2s,
* slaEvp, and slaDranrm
* 12/18/03 bjg 1.3 Change calfusettag.h to calfuse.h
* 06/17/04 bjg Corrected cfortran call to sla_preces
* 07/22/04 bjg 1.4 Remove unused variables
****************************************************************************/
#include <math.h>
#ifdef CFORTRAN
#include "cfortran.h"
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)
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)
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)
#else
#include "slalib.h"
#include "slamac.h"
#endif
#include "calfuse.h"
double state_limb(double pos[3], double mjdate,
double ra, double dec, double *zdist,int *day_limb)
{
double dvb[3], dpb[3], dvh[3], dph[3], limbvec[3], epoch2;
double ra_earth, dec_earth, r, ndist;
double ra_sun, dec_sun;
int i;
char fk5_st[10]="FK5";
ra*=RADIAN;
dec*=RADIAN;
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
/* Now we need the position of the sun (same code as in eclipse.c) */
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
/* ndist is the angular distance (RADIAN) from the nadir to the target */
ndist = (acos(sin(dec)*sin(dec_earth)+
cos(dec)*cos(dec_earth)*cos(ra-ra_earth)));
/* Now we want to calculate the limb vector (originating at the center
of Earth and crossing the target-FUSE line perpendicularily)*/
limbvec[2]=pos[2]+r*cos(ndist)*sin(dec);
limbvec[1]=pos[1]+r*cos(ndist)*cos(dec)*sin(ra);
limbvec[0]=pos[0]+r*cos(ndist)*cos(dec)*cos(ra);
/* Now we do the scalar product with the sun vector; then
if the result is <0 the limb is in the day zone (less
than 90 degrees between sun vector and limb vector)*/
if ((limbvec[0]*cos(dec_sun)*cos(ra_sun) +
limbvec[1]*cos(dec_sun)*sin(ra_sun)+limbvec[2]*sin(dec_sun)) >0) {
*day_limb=1;
} else {
*day_limb=0;
}
*zdist = 180.0 - ndist/RADIAN;
return ndist/RADIAN-(asin(RE/r)/RADIAN);
}
|
