aboutsummaryrefslogtreecommitdiff
path: root/src/libcf/eclipse.c
blob: a1e0b440bb801e96a66da1af931d1f1be8a7dc81 (plain) (blame)
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;
}