diff options
Diffstat (limited to 'src/libcf/read_tle.c')
| -rw-r--r-- | src/libcf/read_tle.c | 303 |
1 files changed, 303 insertions, 0 deletions
diff --git a/src/libcf/read_tle.c b/src/libcf/read_tle.c new file mode 100644 index 0000000..38db069 --- /dev/null +++ b/src/libcf/read_tle.c @@ -0,0 +1,303 @@ +/****************************************************************************** + * Johns Hopkins University + * Center For Astrophysical Sciences + * FUSE + ****************************************************************************** + * + * Synopsis: read_tle(fitsfile *fptr) + * + * Description: read_tle will read the standard tle format file FUSE.TLE + * and place the orbital elements into the header of the + * fitsfile. The FUSE.TLE file will contain multiple orbital + * element sets, so the set closest in time to the observation + * will be used. + * + * Arguments: fitsfile *fptr Pointer to input file + * + * History: 11/11/98 emurphy Begin work. + * 03/18/98 emurphy Added To do list. + * 04/08/99 peb Tidied code and added necessary + * include files + * 06/07/99 peb Added reporting of version control. + * 06/22/99 peb Added FITS_ wrappers. + * 08/05/99 emm Added statements to print date + * of obs and TLE if dtime > 5. + * Also modified program to work with + * FES keywords TEXPSTRT and TEXPEND + * 05/31/00 peb Implemented cfortran.h calls for slalib + * functions. + * 02/13/03 v1.4 wvd Change 0 to NULL in FITS_update_key + * 03/01/03 v1.5 wvd Correct use of pointer in FITS_read_key + * 04/01/03 v1.6 wvd Replace cf_errmsg with cf_if_warning, + * printf with cf_verbose + * 12/18/03 v1.7 bjg Change calfusettag.h to calfuse.h + * 04/07/07 v1.8 wvd Initialize min_mjd to zero to silence + * compiler warnings. + * + *****************************************************************************/ + +#include <math.h> +#include <stdio.h> +#include <stdlib.h> +#include <string.h> + +#ifdef CFORTRAN +#include "cfortran.h" +PROTOCCALLSFSUB5(SLA_CLDJ, sla_cldj, INT, INT, INT, PDOUBLE, PINT) +#define slaCldj(IY, IM, ID, DJM, J) \ + CCALLSFSUB5(SLA_CLDJ, sla_cldj, INT, INT, INT, PDOUBLE, PINT, \ + IY, IM, ID, DJM, J) +#else +#include "slalib.h" +#include "slamac.h" +#endif + +#include "calfuse.h" + +#define MAXCHARS 120 +static char CF_PRGM_ID[] = "read_tle"; +static char CF_VER_NUM[] = "1.8"; + +void read_tle(fitsfile *fptr) +{ + char line1[MAXCHARS], line2[MAXCHARS], line3[MAXCHARS]; + char sat_num[10], security_class[10], international_num[10]; + char inchar[15][15], sp[15][2], sgn[4][2]; + char comment[FLEN_CARD], eccen_str[20]; + char n6_mant_str[20], drag_mant_str[20], instrument[FLEN_CARD]; + int status=0, hdutype=0, card_num, epoch_year, n6_exponent; + int drag_exponent, ephemeris, elset_num, checksum1, checksum2; + int rev_num, e_month, e_day; + double n6_mantissa, drag_mantissa, epoch_day, n2, inclin, raan; + double aop, mean_anom, mean_mot, n6, drag, semimaj, revspersec; + double expstart, expend, expmiddle, eccentr, dtime=1.0e9; + double i_day, f_day, mjd, mjd_f, mjd_d, min_mjd=0; + FILE *ftle; + + /* Enter a timestamp into the log. */ + cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Begin reading TLE file"); + + /* Open the file with the standard two-line elements. */ + ftle = NULL; + ftle = fopen(cf_cal_file("FUSE.TLE"),"r"); + if (ftle == NULL) + cf_if_error("Cannot find file FUSE.TLE, which contains the" + " orbital elements."); + + /* Rewind the fits file and get the start time of the exposure */ + FITS_movabs_hdu(fptr, 1, &hdutype, &status); + FITS_read_key(fptr, TSTRING, "INSTRUME", instrument, comment, &status); + + /* In the primary header, the FES files use keywords with a different + name. This is because the FES may have many exposures in a file. + The TEXPSTRT is the start of the first exposure and TEXPEND is the + end of the last exposure. */ + if (!strncmp(instrument,"FES",3)) { + /* Read in FES keywords */ + FITS_read_key(fptr, TDOUBLE, "TEXPSTRT", &expstart, comment, &status); + FITS_read_key(fptr, TDOUBLE, "TEXPEND", &expend, comment, &status); + } else { + /* Read in FUV keywords */ + FITS_read_key(fptr, TDOUBLE, "EXPSTART", &expstart, comment, &status); + FITS_read_key(fptr, TDOUBLE, "EXPEND", &expend, comment, &status); + } + + /* expmiddle is in mjd */ + expmiddle = (expstart+expend)/2.0; + + /* Read in FUSE.TLE file by getting one line at a time */ + while (fgets(line1, MAXCHARS,ftle) != NULL) { + fgets(line2, MAXCHARS, ftle); + fgets(line3, MAXCHARS, ftle); + +#ifdef DEBUG + printf("\n\n%-80.80s\n",line1); + printf("%-80.80s\n",line2); + printf("%-80.80s\n",line3); +#endif + /* + * This may look like a stupid way to do this, but it has to be + * done this way. You see, the fields can contain leading blanks + * which should be interpreted as zeros. However, c skips over + * leading blanks when reading %f and %d. + */ + sscanf(line2,"%1c%1c%5c%1c%1c%8c%1c%2c%12c%1c%1c%9c%1c%1c%5c%2c%1c" + "%1c%5c%2c%1c%1c%1c%4c%1c", + inchar[1],sp[1],inchar[2],inchar[3],sp[2],inchar[4],sp[3], + inchar[5],inchar[6],sp[4],sgn[1],inchar[7],sp[5],sgn[2], + inchar[8],inchar[9],sp[6],sgn[3],inchar[10],inchar[11], + sp[7],inchar[12], sp[8],inchar[13],inchar[14]); + + inchar[1][1]='\0'; + card_num=atoi(inchar[1]); + + inchar[2][5]='\0'; + strncpy(sat_num,inchar[2],5); + + inchar[3][1]='\0'; + strncpy(security_class,inchar[3],1); + + inchar[4][8]='\0'; + strncpy(international_num,inchar[4],8); + + inchar[5][2]='\0'; + epoch_year=atoi(inchar[5]); + if (epoch_year < 50) + epoch_year+=2000; + else + epoch_year+=1900; + + inchar[6][12]='\0'; + epoch_day=atof(inchar[6]); + + inchar[7][9]='\0'; + n2=atof(inchar[7]); + if (sgn[1][0] == '-') + n2*=-1.0; + + inchar[8][5]='\0'; + strcpy(n6_mant_str,"0."); + strncat(n6_mant_str,inchar[8],7); +#ifdef DEBUG + printf("n6=%-10.10s\n",n6_mant_str); +#endif + n6_mantissa=atof(n6_mant_str); + if (sgn[2][0] == '-') + n6_mantissa*=-1.0; + + inchar[9][2]='\0'; + n6_exponent=atoi(inchar[9]); + + inchar[10][5]='\0'; + strcpy(drag_mant_str,"0."); + strncat(drag_mant_str,inchar[10],7); +#ifdef DEBUG + printf("Drag=%-10.10s\n",drag_mant_str); +#endif + drag_mantissa=atof(drag_mant_str); + if (sgn[3][0] == '-') + drag_mantissa*=-1.0; + + inchar[11][2]='\0'; + drag_exponent=atoi(inchar[11]); + + inchar[12][1]='\0'; + ephemeris=atoi(inchar[12]); + + inchar[13][4]='\0'; + elset_num=atoi(inchar[13]); + + inchar[14][1]='\0'; + checksum1=atoi(inchar[14]); + + n6=n6_mantissa*pow(10.0,(double) n6_exponent); + drag=drag_mantissa*pow(10.0,(double) drag_exponent); + +#ifdef DEBUG + printf("n6=>%10.5f %5d %15.6E\n", n6_mantissa, n6_exponent, n6); + printf("drag=>%10.5f %5d %15.6E\n", drag_mantissa, + drag_exponent, drag); + for(i=1; i<=14; i++) + printf("%2d %-15.15s\n",i,inchar[i]); + for(i=1; i<=3; i++) + printf("%2d %-5.5s\n",i,sgn[i]); +#endif + sscanf(line3,"%1c%1c%5c%1c%8c%1c%8c%1c%7c%1c%8c%1c%8c%1c%11c%5c%1c", + inchar[1],sp[1],inchar[2],sp[2],inchar[3],sp[3],inchar[4], + sp[4],inchar[5],sp[5],inchar[6],sp[6],inchar[7],sp[7], + inchar[8],inchar[9],inchar[10]); + inchar[1][1]='\0'; + card_num=atoi(inchar[1]); + + inchar[2][5]='\0'; + strncpy(sat_num,inchar[2],5); + + inchar[3][8]='\0'; + inclin=atof(inchar[3]); + + inchar[4][8]='\0'; + raan=atof(inchar[4]); + + inchar[5][7]='\0'; + /* + * strcpy(eccen_str,"0."); + * strncat(eccen_str,inchar[5],8); + */ + sprintf(eccen_str,"0.%-8s",inchar[5]); +#ifdef DEBUG + printf("Ecstr=%-20.20s\n",eccen_str); +#endif + eccentr=atof(eccen_str); + + inchar[6][8]='\0'; + aop=atof(inchar[6]); + + inchar[7][8]='\0'; + mean_anom=atof(inchar[7]); + + inchar[8][11]='\0'; + mean_mot=atof(inchar[8]); + + inchar[9][5]='\0'; + rev_num=atof(inchar[9]); + + inchar[10][1]='\0'; + checksum2=atoi(inchar[10]); + +#ifdef DEBUG + for(i=1; i<=10; i++) + printf("%2d %-15.15s\n",i,inchar[i]); + printf("%4d\n%12.8f\n%10.8f\n%10.5E\n%10.5E%5d%5d%5d\n", + epoch_year, epoch_day, n2, n6, drag, ephemeris, + elset_num, checksum1); + printf("\n%8.4f\n%8.4f\n%10.8f\n%8.4f\n%8.4f\n%11.8f\n%5d\n%1d\n", + inclin, raan, (float) eccentr, aop, mean_anom, + mean_mot, rev_num, checksum2); +#endif + /* Convert the year and day of year into a year month day */ + + f_day = modf(epoch_day, &i_day); + + month_day(epoch_year, (int) i_day, &e_month, &e_day); + + /* Convert year month day into a Julian date. */ +#ifdef CFORTRAN + slaCldj(epoch_year, e_month, e_day, mjd, status); +#else + slaCldj(epoch_year, e_month, e_day, &mjd, &status); +#endif + + mjd += f_day; + + if (fabs(expmiddle-mjd) < dtime) { + dtime = fabs(expmiddle-mjd); + min_mjd=mjd; + revspersec = mean_mot/(24.0*60.0*60.0); + semimaj = pow(MU/(4.0*PI*PI*revspersec*revspersec),1.0/3.0); + mjd_f=modf(mjd, &mjd_d); + FITS_update_key(fptr, TDOUBLE, "EPCHTIMD", &mjd_d, NULL, &status); + FITS_update_key(fptr, TDOUBLE, "EPCHTIMF", &mjd_f, NULL, &status); + FITS_update_key(fptr, TDOUBLE, "INCLINAT", &inclin, NULL, &status); + FITS_update_key(fptr, TDOUBLE, "ECCENTRY", &eccentr, NULL, &status); + FITS_update_key(fptr, TDOUBLE, "MEANANOM", &mean_anom, NULL,&status); + FITS_update_key(fptr, TDOUBLE, "ARGPERIG", &aop, NULL, &status); + FITS_update_key(fptr, TDOUBLE, "RASCASCN", &raan, NULL, &status); + FITS_update_key(fptr, TDOUBLE, "SEMIMAJR", &semimaj, NULL, &status); + FITS_update_key(fptr, TDOUBLE, "MEANMOTN", &mean_mot, NULL, &status); + FITS_update_key(fptr, TDOUBLE, "FDM2COEF", &n2, NULL, &status); + FITS_update_key(fptr, TDOUBLE, "SDM6COEF", &n6, NULL, &status); + FITS_update_key(fptr, TDOUBLE, "DRAGCOEF", &drag, NULL, &status); + FITS_update_key(fptr, TSTRING, "PROPMODL", "SGP4", NULL, &status); + } + } /* endwhile */ + + if (dtime > 5.0) { + cf_verbose(2, "MJD OBS=%12.5f MJD TLE=%12.5f DT=%12.5f\n", + expmiddle, min_mjd, dtime); + cf_if_warning("Orbital elements are more than 5 days old."); + } + + /* Enter a timestamp into the log. */ + cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Done reading TLE file"); +} + |