diff options
Diffstat (limited to 'src/libcf/cf_velang.c')
| -rw-r--r-- | src/libcf/cf_velang.c | 147 |
1 files changed, 147 insertions, 0 deletions
diff --git a/src/libcf/cf_velang.c b/src/libcf/cf_velang.c new file mode 100644 index 0000000..8395738 --- /dev/null +++ b/src/libcf/cf_velang.c @@ -0,0 +1,147 @@ +/***************************************************************************** + * Johns Hopkins University + * Center For Astrophysical Sciences + * FUSE + ***************************************************************************** + * + * Synopsis: cf_velang(fitsfile *outfits, double mjd) + * + * Description: Calculates the sun angle, moon angle, and conversions to + * heliocentric and LSR velocities at given mjd, which should be + * the middle of the exposure. For FUV images, the needed + * data is in the primary header, so we have to go back + * up to the top. For FES images, the data should be in + * each HDU, so you should call cf_velang_calc instead. + * + * Arguments: fitsfile *outfits Pointer to FITS file containing the + * input data and orbital elements + * double mjd The Modified Julian Date of the + * photon arrival time. + * + * Calls: slaDsep dsep.f + * slaRdplan rdplan.f + * + * History: 07/13/98 emm Begin work + * 06/10/99 emm Added keyword header update, added + * calculation of MOONANG. + * 06/11/99 emm Fixed some bugs that Jayant found. + * 06/21/99 jm Corrected name + * 07/16/99 emm Changed GEO_LON to GEO_LONG + * 07/21/99 peb RA and Dec now use RA_TARG and DEC_TARG + * 07/22/99 emm Added V_GEOCEN update. + * 07/23/99 emm Added cf_velang_calc for FES images. + * 05/31/00 peb Implemented cfortran.h calls for slalib + * functions. + * 04/01/03 v1.3 wvd Replaced printf with cf_verbose + * 04/04/03 v1.4 wvd Modifed calls to cf_verbose. + * 12/18/03 v1.5 bjg Change calfusettag.h to calfuse.h + * 07/06/06 v1.6 wvd Add call to pole_ang and populate + * header keyword POLEANGL. + * 03/07/07 v1.7 wvd Remove pos from call to space_vel. + * + ****************************************************************************/ + +#include <stdio.h> + +#ifdef CFORTRAN +#include "cfortran.h" +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) +PROTOCCALLSFSUB7(SLA_RDPLAN, sla_rdplan, DOUBLE, INT, DOUBLE, DOUBLE, \ + PDOUBLE, PDOUBLE, PDOUBLE) +#define slaRdplan(DATE, NP, ELONG, PHI, RA, DEC, DIAM) \ + CCALLSFSUB7(SLA_RDPLAN, sla_rdplan, DOUBLE, INT, DOUBLE, DOUBLE, \ + PDOUBLE, PDOUBLE, PDOUBLE, \ + DATE, NP, ELONG, PHI, RA, DEC, DIAM) +#else +#include "slalib.h" +#include "slamac.h" +#endif + +#include "calfuse.h" +#include "sgp4.h" + +SGP4 set_orbit_parms(fitsfile *); + +static void +cf_velang_calc(fitsfile *outfits, double mjd) +{ + int status=0; + char comment[FLEN_CARD]; + double ra, dec, geo_lon, geo_lat, mag_lat, pos[3], vel[3], gmst, helvel; + double ra_moon, dec_moon, moon_ang, dia_moon, vlsr1, vlsr2, sun_ang; + double geo_vel, pole_angle; + SGP4 sgp4; + + /* get the state vector at time mjd */ + sgp4 = set_orbit_parms(outfits); + SGP4_getStateVector(sgp4, mjd, pos, vel); + SGP4_precess(pos, mjd, MJD2000); SGP4_precess(vel, mjd, MJD2000); + + cf_verbose(2,"State vector at MJD=%14.6f days", mjd); + cf_verbose(2,"x=%14.6f y=%14.6f z=%14.6f km", pos[0], pos[1], pos[2]); + cf_verbose(2,"vx=%14.6f vy=%14.6f vz=%14.6f km/s",vel[0], vel[1], vel[2]); + + state_geod(pos, mjd, &geo_lon, &geo_lat, &gmst); + + mag_lat = geod_mag(geo_lon, geo_lat); + /* + * Note: this function has been changed to use RA_TARG and DEC_TARG + * keywords, because the RA_APER and DEC_APER keywords are currently + * not being set. + */ + FITS_read_key(outfits, TDOUBLE, "RA_TARG", &ra, comment, &status); + FITS_read_key(outfits, TDOUBLE, "DEC_TARG", &dec, comment, &status); + + geo_vel= space_vel(vel, ra, dec); + helvel = helio_vel(mjd, ra, dec); + vlsr1 = lsrk_vel(ra, dec); + vlsr2 = lsrd_vel(ra, dec); + + /* Calculate sun angle */ + sun_ang = solar_ang(mjd, ra, dec); + + /* Calculate moon angle */ +#ifdef CFORTRAN + slaRdplan(mjd, 3, geo_lon, geo_lat, ra_moon, dec_moon, dia_moon); +#else + slaRdplan(mjd, 3, geo_lon, geo_lat, &ra_moon, &dec_moon, &dia_moon); +#endif + moon_ang = slaDsep(ra*RADIAN, dec*RADIAN, ra_moon, dec_moon)/RADIAN; + + /* Calculate sun angle */ + pole_angle = pole_ang(pos, vel, ra, dec); + + FITS_update_key(outfits, TDOUBLE, "SUNANGLE", &sun_ang, NULL, &status); + FITS_update_key(outfits, TDOUBLE, "MOONANGL", &moon_ang, NULL, &status); + FITS_update_key(outfits, TDOUBLE, "POLEANGL", &pole_angle, NULL, &status); + FITS_update_key(outfits, TDOUBLE, "GEO_LONG", &geo_lon, NULL, &status); + FITS_update_key(outfits, TDOUBLE, "GEO_LAT", &geo_lat, NULL, &status); + FITS_update_key(outfits, TDOUBLE, "MAG_LAT", &mag_lat, NULL, &status); + FITS_update_key(outfits, TDOUBLE, "V_GEOCEN", &geo_vel, NULL, &status); + FITS_update_key(outfits, TDOUBLE, "V_HELIO", &helvel, NULL, &status); + FITS_update_key(outfits, TDOUBLE, "V_LSRDYN", &vlsr2, NULL, &status); + FITS_update_key(outfits, TDOUBLE, "V_LSRSTD", &vlsr1, NULL, &status); +} + + +void cf_velang(fitsfile *outfits, double mjd) +{ + int status=0, hdunum=0, hdutype; + /* + * The outfits file pointer may be down in one of the extensions + * whereas the orbital info is in the primary hdu. Therefore, + * record the current position, move to the primary header, + * read the orbital data, and return to the previous hdu. + */ + FITS_get_hdu_num(outfits, &hdunum); + FITS_movabs_hdu(outfits, 1, &hdutype, &status); + + cf_velang_calc(outfits, mjd); + + FITS_movabs_hdu(outfits, hdunum, &hdutype, &status); + +} |