diff options
Diffstat (limited to 'src/libcf/eclipse.c')
| -rw-r--r-- | src/libcf/eclipse.c | 148 |
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; +} |