Initial commit

1 files changed, 148 insertions, 0 deletions

diff --git a/src/libcf/eclipse.c b/src/libcf/eclipse.c
new file mode 100644
index 0000000..a1e0b44
--- /dev/null
+++ b/src/libcf/eclipse.c
@@ -0,0 +1,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;
+}