include include include include "delaytime.h" define SZ_HIST SZ_TIME+3*SZ_FNAME+200 # size of the history log define NLAN_EARTH 10 define NLAN_OBS 10 # DELAYTIME_DO -- Apply correction for light delay time # # 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 #------------------------------------------------------------------------------ procedure delaytime_do (fin, fout, nfin, mu_ra, mu_dec, parallax, tp_earth, tp_obs) pointer fin, fout # input: file template pointers int nfin # input: number of files in the input template double mu_ra, mu_dec # input: proper motion of the target at J2000 # (in deg./year) double parallax # input: parallax of the target (in arc sec) pointer tp_earth # input: pointer to the earth ephemeris table pointer tp_obs # input: pointer to the observer ephemeris table double ra, dec # RA and Dec of the target at J2000 (in deg.) pointer arrout # addresses of output data 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 double objvec[3] # unit geocentric state vector of the target double mjd1, new_mjd1 # MJD of the first pixel obs time (beginning of the # integration not the middle of first pixel) double epoch, geomdelt, reldelt, avetime double sample_time real cd1_1 pointer ipin, ipout char ifile[SZ_FNAME], ofile[SZ_FNAME] int npix # number of points of the input file int nchar, i, j pointer sp, ip pointer x_e, y_e, z_e, x_o, y_o, z_o int npts_earth, npts_obs double first_earth, first_obs double grid_earth, grid_obs char text[SZ_HIST] char mess[SZ_LINE] pointer immap() pointer impl1d() double imgetd() real imgetr() int tbpsta() int imtgetim() long clktime() #============================================================================== begin # find out how large are the ephemerides npts_earth = tbpsta (tp_earth, TBL_NROWS) npts_obs = tbpsta (tp_obs, TBL_NROWS) # allocate space for the ephemerides call smark (sp) call salloc (x_e, npts_earth, TY_DOUBLE) call salloc (y_e, npts_earth, TY_DOUBLE) call salloc (z_e, npts_earth, TY_DOUBLE) call salloc (x_o, npts_obs, TY_DOUBLE) call salloc (y_o, npts_obs, TY_DOUBLE) call salloc (z_o, npts_obs, TY_DOUBLE) # read the ephemerides call get_ephem (tp_earth, Memd[x_e], Memd[y_e], Memd[z_e], npts_earth, first_earth, grid_earth) 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) nchar = imtgetim (fout, ofile, SZ_FNAME) # open input/output files ipin = immap (ifile, READ_ONLY, 0) ipout = immap (ofile, NEW_FILE, 0) npix = IM_LEN(ipin, 1) IM_LEN(ipout, 1) = npix IM_PIXTYPE(ipout) = TY_DOUBLE # allocate output array space call malloc (arrout, npix, TY_DOUBLE) # input file must be 1-D data if (IM_NDIM(ipin) != 1) call error (1, "input is not 1-D data") sample_time = imgetd (ipin, "SAMPTIME") # in sec cd1_1 = imgetr (ipin, "CD1_1") # in sec ra = imgetd (ipin, "RA_TARG") dec = imgetd (ipin, "DEC_TARG") call get_mjd1 (ipin, mjd1) # calculate the average time of this observation avetime = mjd1 + npix * cd1_1 / (2.d0 * SECPERDAY) # in MJD # calculate the target's positional vector, including proper motion call object_pos (avetime, ra, dec, mu_ra, mu_dec, objvec) # calculate delaytime due to relativistic effects call relativ (avetime, reldelt) # calculate the geometric delay time, point by point do i = 1, npix { epoch = mjd1 + ((i-1) * cd1_1 + 0.5d0 * sample_time) / SECPERDAY iferr (call intrp_state (epoch, Memd[x_e], Memd[y_e], Memd[z_e], first_earth, grid_earth, xyz_earth, npts_earth, NLAN_EARTH)) { call sprintf (mess, SZ_LINE, "epoch = MJD %f, outside the earth_ephem range") call pargd (epoch) call error (1, mess) } iferr (call intrp_state (epoch, Memd[x_o], Memd[y_o], Memd[z_o], first_obs, grid_obs, xyz_obs, npts_obs, NLAN_OBS)) { call sprintf (mess, SZ_LINE, "epoch = MJD %f, outside the obs_ephem range") call pargd (epoch) call error (1, mess) } call aaddd (xyz_obs, xyz_earth, telvec, 3) call geo_delay (telvec, objvec, parallax, geomdelt) epoch = (geomdelt + reldelt) / SECPERDAY + epoch # MJD if (i == 1) new_mjd1 = epoch Memd[arrout+i-1] = (epoch - new_mjd1) * SECPERDAY # sec } # put the result to the output files and close them ip = impl1d (ipout) call amovd (Memd[arrout], Memd[ip], npix) # update header keyword(s) and history call imastr (ipout, "BUNIT", "PIXEL", SZ_BUNIT) call imaddd (ipout, "REFMJD", new_mjd1) call cnvtime (clktime(0), text, SZ_HIST) call strcat (" created by the task DELAYTIME, from the science file: ", text, SZ_HIST) call strcat (ifile, text, SZ_HIST) call imputh (ipout, "HISTORY", text) # close files call imunmap (ipin) call imunmap (ipout) # print out message of which files been created call printf ("delaytime: output file %s is created\n") call pargstr (ofile) call mfree (arrout, TY_DOUBLE) } call sfree (sp) end