include include include "delaytime.h" define NLAN_EARTH 10 define NLAN_OBS 10 define INTERVAL 0.1 # in second # ODELAY_DO -- Apply correction for light delay time for STIS data # # Description: # ------------ # Compute the times of observation relative to the solar-system barycenter # of a distant object (outside the solar system). It is assumed that # the set of observations covers a short enough time that the object does # not move appreciably during the observations. Note that a one-arcsecond # error in the position of the object can result in a time-delay error # of ~2.4 millisec. It is also assumed that the set of observations covers # a short enough time interval that the relativistic correction, if required, # remains constant during the observations. # # Date Author Description # ---- ------ ----------- # 10-Nov-1984 C. D. Biemesderfer Original module # 18-Apr-1990 J.-C. Hsu rewrite in SPP # 19-Jun-1992 J.-C. Hsu read RA_TARG and DEC_TARG from header # 19-Aug-1997 J.-C. Hsu modify for STIS data # 23-Aug-2000 Phil Hodge remove out_col & add verbose arguments; # modify times in-place; update GTI table # and header keywords as well; # change INTERVAL from 1. to 0.1; # add DELAYCOR and HISTORY #------------------------------------------------------------------------------ procedure odelay_do (fin, nfin, parallax, tp_earth, tp_obs, in_col, verbose) pointer fin # i: file template pointer int nfin # i: number of files in the input template double parallax # i: parallax of the target (in arc sec) pointer tp_earth # i: pointer to the earth ephemeris table pointer tp_obs # i: pointer to the observer ephemeris table char in_col[ARB] # i: time column name bool verbose # i: print timing info? #-- pointer sp pointer history # for constructing a history record pointer orx_name # name of obs_ephem table, for history double ra, dec # RA and Dec of the target at J2000 (in deg.) double objvec[3] # unit geocentric state vector of the target double mjd1, mjd2 # MJD of the observation start/end time double epoch # MJD of an event double delta_sec # correction (sec) to be added to TIME col double t0_delay # correction at TEXPSTRT char ifile[SZ_FNAME] char ifilen[SZ_FNAME] int nchar, i, j, k pointer x_e, y_e, z_e, x_o, y_o, z_o pointer tp # pointer of the input table pointer cptr # for the TIME column in EVENTS tables pointer cp1, cp2 # for START and STOP columns in GTI table int nrows # number of rows in a table int nevents_tab # number of EVENTS tables (one less than NEXTEND) int npts_earth, npts_obs double first_earth, first_obs double last_earth, last_obs double grid_earth, grid_obs double tm # a time from the TIME column double tm_prev # the last tm for which all_delay was called char text[SZ_TIME] char delaycorr[SZ_FNAME] int frac bool obs_ephem pointer tbtopn() int tbpsta() double tbhgtd() int tbhgti() int imtgetim() long clktime() int strlen() bool streq() #============================================================================== begin call smark (sp) call salloc (history, SZ_FNAME, TY_CHAR) call salloc (orx_name, SZ_FNAME, TY_CHAR) # find out how large are the ephemerides and # allocate space for the ephemerides npts_earth = tbpsta (tp_earth, TBL_NROWS) call calloc (x_e, npts_earth, TY_DOUBLE) call calloc (y_e, npts_earth, TY_DOUBLE) call calloc (z_e, npts_earth, TY_DOUBLE) obs_ephem = (tp_obs != NULL) if (obs_ephem) { npts_obs = tbpsta (tp_obs, TBL_NROWS) call calloc (x_o, npts_obs, TY_DOUBLE) call calloc (y_o, npts_obs, TY_DOUBLE) call calloc (z_o, npts_obs, TY_DOUBLE) } # read the ephemerides, and close the tables call get_ephem (tp_earth, Memd[x_e], Memd[y_e], Memd[z_e], npts_earth, first_earth, grid_earth) if (obs_ephem) { call tbtnam (tp_obs, Memc[orx_name], SZ_FNAME) call get_ephem (tp_obs, Memd[x_o], Memd[y_o], Memd[z_o], npts_obs, first_obs, grid_obs) } # loop all input files do j = 1, nfin { # read the next file name in the template list nchar = imtgetim (fin, ifile, SZ_FNAME) if (strlen (ifile) <= 0) next if (verbose) { call printf ("odelaytime: processing %s ...\n") call pargstr (ifile) call flush(STDOUT) } call strcpy (ifile, ifilen, SZ_FNAME) call strcat ("[0]", ifilen, SZ_FNAME) # open the primary header iferr (tp = tbtopn (ifilen, READ_ONLY, 0)) { call printf ("Cannot open the input table '%s'\n") call pargstr (ifilen) call error (1, "") } # check whether this table has already been corrected iferr { call tbhgtt (tp, "DELAYCOR", delaycorr, SZ_FNAME) } then { call strcpy ("PERFORM", delaycorr, SZ_FNAME) } if (streq (delaycorr, "COMPLETE")) { call printf ("%s has already been corrected,\n") call pargstr (ifile) call printf ( " so no further processing will be applied to this table.\n") call flush (STDOUT) call tbtclo (tp) next } # get these for checking the range and for updating the GTI table mjd1 = tbhgtd (tp, "TEXPSTRT") mjd2 = tbhgtd (tp, "TEXPEND") ra = tbhgtd (tp, "RA_TARG") dec = tbhgtd (tp, "DEC_TARG") iferr { # assume (for now) that the last extension is a GTI table nevents_tab = tbhgti (tp, "NEXTEND") - 1 } then { nevents_tab = 1 } call tbtclo (tp) # check the range last_earth = first_earth + (npts_earth - 1) * grid_earth if (mjd1 < first_earth || mjd2 > last_earth ) { call printf ("epoch outside the Earth ephemeris range\n") next } if (obs_ephem) { last_obs = first_obs + (npts_obs - 1) * grid_obs if (mjd1 < first_obs || mjd2 > last_obs) { call printf ("epoch outside the obs ephemeris range\n") next } } # calculate the target's positional vector call object_pos (ra, dec, objvec) # Get the delay at the exposure start time. We'll subtract this # delay from each time that we update, if that time is relative # to EXPSTART (which must be the same as TEXPSTRT). call all_delay (mjd1, parallax, objvec, Memd[x_e], Memd[y_e], Memd[z_e], npts_earth, first_earth, grid_earth, obs_ephem, Memd[x_o], Memd[y_o], Memd[z_o], npts_obs, first_obs, grid_obs, t0_delay) # update times in GTI table call strcpy (ifile, ifilen, SZ_FNAME) call strcat ("[GTI]", ifilen, SZ_FNAME) iferr (tp = tbtopn (ifilen, READ_WRITE, 0)) { call printf ("Warning: %s has no GTI extension\n") call pargstr (ifile) # if there's no GTI, last extension is an EVENTS table nevents_tab = nevents_tab + 1 } else { call tbcfnd1 (tp, "START", cp1) call tbcfnd1 (tp, "STOP", cp2) if (cp1 == NULL || cp2 == NULL) { call tbtclo (tp) call error (1, "either START or STOP column was not found in GTI extension") } nrows = tbpsta (tp, TBL_NROWS) do i = 1, nrows { call tbegtd (tp, cp1, i, tm) epoch = mjd1 + tm / SECPERDAY call all_delay (epoch, parallax, objvec, Memd[x_e], Memd[y_e], Memd[z_e], npts_earth, first_earth, grid_earth, obs_ephem, Memd[x_o], Memd[y_o], Memd[z_o], npts_obs, first_obs, grid_obs, delta_sec) call tbeptd (tp, cp1, i, tm + (delta_sec - t0_delay)) call tbegtd (tp, cp2, i, tm) epoch = mjd1 + tm / SECPERDAY call all_delay (epoch, parallax, objvec, Memd[x_e], Memd[y_e], Memd[z_e], npts_earth, first_earth, grid_earth, obs_ephem, Memd[x_o], Memd[y_o], Memd[z_o], npts_obs, first_obs, grid_obs, delta_sec) call tbeptd (tp, cp2, i, tm + (delta_sec - t0_delay)) } if (verbose) { call printf (" GTI extension has been updated\n") call flush(STDOUT) } } # loop through all EVENTS extensions do k = 1, nevents_tab { call sprintf (ifilen, SZ_FNAME, "%s[EVENTS,%d]") call pargstr (ifile) call pargi (k) iferr (tp = tbtopn (ifilen, READ_WRITE, 0)) { call printf ("Cannot open the input table '%s'\n") call pargstr (ifilen) next } mjd1 = tbhgtd (tp, "EXPSTART") mjd2 = tbhgtd (tp, "EXPEND") # find the time (since EXPSTART) column in the table call tbcfnd1 (tp, in_col, cptr) if (cptr == NULL) { call printf ("column %s not found in %s\n") call pargstr (in_col) call pargstr (ifilen) call error (1, "") } # find how many rows are there nrows = tbpsta (tp, TBL_NROWS) tm_prev = 0. # just to have a definite starting value # print out timing info if (verbose) { call cnvtime (clktime(0), text, SZ_TIME) call printf (" start time: %s\n") call pargstr(text) call flush (STDOUT) frac = nrows / 10 call printf (" Percentage done: ") call flush(STDOUT) } # go through each row do i = 1, nrows { if (verbose) { if (mod(i,frac) == 0) { call printf (" %d") call pargi(i/frac*10) call flush(STDOUT) } } call tbegtd (tp, cptr, i, tm) # if the time is within INTERVAL of the previous time, # apply the same delaytime if (i == 1 || (tm - tm_prev) > INTERVAL) { epoch = mjd1 + tm / SECPERDAY call all_delay (epoch, parallax, objvec, Memd[x_e], Memd[y_e], Memd[z_e], npts_earth, first_earth, grid_earth, obs_ephem, Memd[x_o], Memd[y_o], Memd[z_o], npts_obs, first_obs, grid_obs, delta_sec) tm_prev = tm } # write the corrected time back to the time column call tbeptd (tp, cptr, i, tm + (delta_sec - t0_delay)) } # add delaytime to EXPSTART and EXPEND, and update header call all_delay (mjd1, parallax, objvec, Memd[x_e], Memd[y_e], Memd[z_e], npts_earth, first_earth, grid_earth, obs_ephem, Memd[x_o], Memd[y_o], Memd[z_o], npts_obs, first_obs, grid_obs, delta_sec) call tbhptd (tp, "EXPSTART", mjd1 + delta_sec/SECPERDAY) call all_delay (mjd2, parallax, objvec, Memd[x_e], Memd[y_e], Memd[z_e], npts_earth, first_earth, grid_earth, obs_ephem, Memd[x_o], Memd[y_o], Memd[z_o], npts_obs, first_obs, grid_obs, delta_sec) call tbhptd (tp, "EXPEND", mjd2 + delta_sec/SECPERDAY) # close input table call tbtclo (tp) if (verbose) { call printf ("\n") # finish the "Percentage done" call printf (" [EVENTS,%d] extension has been updated\n") call pargi (k) call cnvtime (clktime(0), text, SZ_TIME) call printf (" finish time: %s\n") call pargstr(text) call flush(STDOUT) } } # add delaytime to TEXPSTRT and TEXPEND, and update primary header call strcpy (ifile, ifilen, SZ_FNAME) call strcat ("[0]", ifilen, SZ_FNAME) tp = tbtopn (ifilen, READ_WRITE, 0) mjd1 = tbhgtd (tp, "TEXPSTRT") call all_delay (mjd1, parallax, objvec, Memd[x_e], Memd[y_e], Memd[z_e], npts_earth, first_earth, grid_earth, obs_ephem, Memd[x_o], Memd[y_o], Memd[z_o], npts_obs, first_obs, grid_obs, delta_sec) call tbhptd (tp, "TEXPSTRT", mjd1 + delta_sec/SECPERDAY) mjd2 = tbhgtd (tp, "TEXPEND") call all_delay (mjd2, parallax, objvec, Memd[x_e], Memd[y_e], Memd[z_e], npts_earth, first_earth, grid_earth, obs_ephem, Memd[x_o], Memd[y_o], Memd[z_o], npts_obs, first_obs, grid_obs, delta_sec) call tbhptd (tp, "TEXPEND", mjd2 + delta_sec/SECPERDAY) # add keyword to flag the fact that the times have been corrected call tbhadt (tp, "DELAYCOR", "COMPLETE") call tbhpcm (tp, "DELAYCOR", "delaytime has been applied") call tbhadt (tp, "HISTORY", "Times corrected to solar system barycenter;") if (obs_ephem) { call sprintf (Memc[history], SZ_FNAME, "the ORX table was %s") call pargstr (Memc[orx_name]) } else { call strcpy ("no ORX table was used.", Memc[history], SZ_FNAME) } call tbhadt (tp, "HISTORY", Memc[history]) call tbtclo (tp) if (verbose) { call printf ("... done\n") call flush(STDOUT) } } # free the memory call mfree (x_e, TY_DOUBLE) call mfree (y_e, TY_DOUBLE) call mfree (z_e, TY_DOUBLE) if (obs_ephem) { call mfree (x_o, TY_DOUBLE) call mfree (y_o, TY_DOUBLE) call mfree (z_o, TY_DOUBLE) } call sfree (sp) end # This routine computes the delaytime in seconds. procedure all_delay (epoch, parallax, objvec, x_e, y_e, z_e, npts_earth, first_earth, grid_earth, obs_ephem, x_o, y_o, z_o, npts_obs, first_obs, grid_obs, delta_sec) double epoch # i: time (MJD) of event double parallax # i: parallax of target double objvec[3] # i: unit geocentric state vector of the target double x_e[ARB], y_e[ARB], z_e[ARB] int npts_earth bool obs_ephem double first_earth, grid_earth double x_o[ARB], y_o[ARB], z_o[ARB] int npts_obs double first_obs, grid_obs double delta_sec # o: correction (sec) to be added to TIME col #-- char mess[SZ_FNAME] # for error message double geomdelt, reldelt # time corrections (sec) double xyz_obs[3] # geocentric state vector of the observer double xyz_earth[3] # baycentric state vector of the earth double telvec[3] # baycentric state vector of the observer begin # calculate delaytime due to relativistic effects call relativ (epoch, reldelt) # calculate the geometric delay time iferr (call intrp_state (epoch, x_e, y_e, z_e, first_earth, grid_earth, xyz_earth, npts_earth, NLAN_EARTH)) { call sprintf (mess, SZ_FNAME, "epoch = MJD %f, outside the earth_ephem range") call pargd (epoch) call error (1, mess) } if (obs_ephem) { iferr (call intrp_state (epoch, x_o, y_o, z_o, first_obs, grid_obs, xyz_obs, npts_obs, NLAN_OBS)) { call sprintf (mess, SZ_FNAME, "epoch = MJD %f, outside the obs_ephem range") call pargd (epoch) call error (1, mess) } } else { xyz_obs[1] = 0.d0 xyz_obs[2] = 0.d0 xyz_obs[3] = 0.d0 } call aaddd (xyz_obs, xyz_earth, telvec, 3) call geo_delay (telvec, objvec, parallax, geomdelt) delta_sec = geomdelt + reldelt end