include include define SZ_PARAM 18 procedure waveoffset () # WAVEOFFSET -- Compute the offset of one spectrum from a reference spectrum # # Version 1.0 May 91 S. Hulbert Original # 1.1 Jul 91 S. Hulbert Fixed bug in writing table header char obs1[SZ_FNAME] # I: image names for first observation char wave1[SZ_FNAME] # I: image names for wavelength observation char obs2[SZ_FNAME] # I: image names for second observation char cross[SZ_FNAME] # I/O: image name for correlation char table[SZ_FNAME] # I/O: image name for output table int search # I: max distance in pixel space to search int sections # I: number of sections int i, pixin1, pixinw, pixin2, ioff, n, ic, width int fchnl1, nchnl1, xstep1, over1, pass1 int fchnl2, nchnl2, xstep2, over2, pass2 int ndelw, ndeld char det1[SZ_PARAM], det2[SZ_PARAM] char fgwa1[SZ_PARAM], aper1[SZ_PARAM], aperpos1[SZ_PARAM], polar1[SZ_PARAM] char fgwa2[SZ_PARAM], aper2[SZ_PARAM], aperpos2[SZ_PARAM], polar2[SZ_PARAM] pointer imobs1, imwave1, imobs2, bufobs1, bufobs2, bufwave1, imout, bufcc pointer diode, tp, cp_wave, cp_diod, cp_delw, cp_deld, cp_cnts pointer deltad, counts, wave, deltaw real mean_delw, mean_deld, sig_delw, sig_deld bool strne() int clgeti(), imgl1r(), aravr() pointer immap(), imps2r(), tbtopn() pointer sp begin call smark (sp) # get input from cl call clgstr ("input1", obs1, SZ_FNAME) call clgstr ("wave1", wave1, SZ_FNAME) call clgstr ("input2", obs2, SZ_FNAME) call clgstr ("table", table, SZ_FNAME) call clgstr ("cross", cross, SZ_FNAME) search = clgeti ("maxdist") sections = clgeti ("nbins") # map the images imobs1 = immap (obs1, READ_ONLY, 0) imwave1 = immap (wave1, READ_ONLY, 0) imobs2 = immap (obs2, READ_ONLY, 0) imout = immap (cross, NEW_COPY, imobs1) # set up output dimensions IM_NDIM(imout) = 2 width = search * 2 + 1 IM_LEN(imout, 1) = width IM_LEN(imout, 2) = sections # get size input pixin1 = IM_LEN (imobs1, 1) pixinw = IM_LEN (imwave1, 1) if (pixin1 != pixinw) call error (1, "Flux and wavelength images are not the same size") pixin2 = IM_LEN (imobs2, 1) if (pixin1 != pixinw) call error (1, "Flux images are not the same size") # get pattern and mode info and check for agreement call xpattern (imobs1, fchnl1, nchnl1, xstep1, over1) call obsmode (imobs1, det1, fgwa1, aper1, aperpos1, polar1, pass1) call xpattern (imobs2, fchnl2, nchnl2, xstep2, over2) call obsmode (imobs2, det2, fgwa2, aper2, aperpos2, polar2, pass2) # check to see the images are the same kind if (strne (det1, det2)) call error (1, "Detectors do not match") if (strne (fgwa1, fgwa2)) call error (1, "Gratings do not match") if (fchnl1 != fchnl2) call error (1, "First channels do not match") if (nchnl1 != nchnl2) call error (1, "Number of channels do not match") if (xstep1 != xstep2) call error (1, "X-steps do not match") # set up and fill the buffers call salloc (bufobs1, pixin1, TY_REAL) call salloc (bufobs2, pixin2, TY_REAL) call salloc (bufwave1, pixinw, TY_REAL) call amovr (Memr[imgl1r (imobs1)], Memr[bufobs1], pixin1) call amovr (Memr[imgl1r (imobs2)], Memr[bufobs2], pixin2) call amovr (Memr[imgl1r (imwave1)], Memr[bufwave1], pixinw) # do cross-correlation call salloc (bufcc, width * sections, TY_REAL) call salloc (deltad, sections, TY_REAL) call salloc (counts, sections, TY_REAL) call ycrscor (Memr[bufobs1], Memr[bufobs2], pixin1, sections, search, Memr[deltad], Memr[counts], Memr[bufcc]) # write correlation image call amovr (Memr[bufcc], Memr[imps2r (imout, 1, width, 1, sections)], width * sections) # process to put in table call salloc (wave, sections, TY_REAL) call salloc (diode, sections, TY_REAL) call salloc (deltaw, sections, TY_REAL) n = pixin1 / sections do i = 1, sections { ioff = (i - 1) * n ic = ioff + n / 2 + 1 Memr[diode + i - 1] = real (ic) / xstep1 + fchnl1 - 1 Memr[wave + i - 1] = Memr[bufwave1 + ic - 1] if (Memr[deltad + i - 1] != INDEF) { Memr[deltaw + i - 1] = Memr[deltad + i - 1] * (Memr[bufwave1 + (ic + 1) - 1] - Memr[bufwave1 + (ic - 1) - 1)) / 2.0 Memr[deltad + i - 1] = Memr[deltad + i - 1] / xstep1 } else { Memr[deltaw + i - 1] = INDEF } } # write table tp = tbtopn (table, NEW_FILE, 0) # make it a column ordered table call tbpset (tp, TBL_WHTYPE, TBL_TYPE_S_COL) # set up the columns call tbcdef (tp, cp_wave, "WAVELENGTH", "ANGSTROMS", "", TY_REAL, 1, 1) call tbcdef (tp, cp_diod, "DIODE", "", "", TY_REAL, 1, 1) call tbcdef (tp, cp_delw, "DELTAW", "ANGSTROMS", "", TY_REAL, 1, 1) call tbcdef (tp, cp_deld, "DELTAD", "", "", TY_REAL, 1, 1) call tbcdef (tp, cp_cnts, "COUNTS", "", "", TY_REAL, 1, 1) # create the table call tbtcre (tp) # write the data call tbcptr (tp, cp_wave, Memr[wave], 1, sections) call tbcptr (tp, cp_diod, Memr[diode], 1, sections) call tbcptr (tp, cp_delw, Memr[deltaw], 1, sections) call tbcptr (tp, cp_deld, Memr[deltad], 1, sections) call tbcptr (tp, cp_cnts, Memr[counts], 1, sections) # add the header keywords call tbhadt (tp, "DETECTOR", det2) call tbhadt (tp, "FGWA_ID", fgwa2) call tbhadt (tp, "APER_ID", aper2) call tbhadt (tp, "APER_POS", aperpos2) call tbhadt (tp, "POLAR_ID", polar2) call tbhadi (tp, "PASS_DIR", pass2) call tbhadt (tp, "APER_REF", aper1) call tbhadt (tp, "APOS_REF", aperpos1) call tbhadt (tp, "POLR_REF", polar1) call tbhadi (tp, "PASS_REF", pass1) call tbtclo (tp) # write to STDOUT call printf (" Wavelength Diode DeltaW DeltaD Total counts\n") call printf (" 2nd spec.\n") do i = 1, sections { call printf (" %8.3f %5.1f %8.4f %8.4f %10.1f\n") call pargr (Memr[wave + i - 1]) call pargr (Memr[diode + i - 1]) call pargr (Memr[deltaw + i - 1]) call pargr (Memr[deltad + i - 1]) call pargr (Memr[counts + i - 1]) } # find the mean values for deta wavelength and delta diode # reject at 2 sigma level ndelw = aravr (Memr[deltaw], sections, mean_delw, sig_delw, 2.0) ndeld = aravr (Memr[deltad], sections, mean_deld, sig_deld, 2.0) call printf ("\nmean DeltaW= %8.4f SD= %8.4f %2d values rejected\n") call pargr (mean_delw) call pargr (sig_delw) call pargi (sections - ndelw) call printf ("mean DeltaD= %8.4f SD= %8.4f %2d values rejected\n") call pargr (mean_deld) call pargr (sig_deld) call pargi (sections - ndeld) call sfree (sp) call imunmap (imobs1) call imunmap (imwave1) call imunmap (imobs2) call imunmap (imout) end procedure xpattern (im, fchnl, nchnls, nxsteps, overscan) # XPATTERN -- get x-pattern values from FOS header pointer im # image descriptor int fchnl # O int nchnls # O int nxsteps # O int overscan # O int imgeti() begin fchnl = imgeti (im, "FCHNL") nchnls = imgeti (im, "NCHNLS") nxsteps = imgeti (im, "NXSTEPS") overscan = imgeti (im, "OVERSCAN") end define SZ_PARAM 18 procedure obsmode (im, detector, fgwaid, aperid, aperpos, polarid, passdir) # OBSMODE -- get observing mode values from FOS header pointer im # I: image descriptor char detector[SZ_PARAM] # O char fgwaid[SZ_PARAM] # O char aperid[SZ_PARAM] # O char aperpos[SZ_PARAM] # O char polarid[SZ_PARAM] # O int passdir # O int imgeti() begin call imgstr (im, "DETECTOR", detector, SZ_PARAM) call imgstr (im, "FGWA_ID", fgwaid, SZ_PARAM) call imgstr (im, "APER_ID", aperid, SZ_PARAM) call imgstr (im, "APER_POS", aperpos, SZ_PARAM) call imgstr (im, "POLAR_ID", polarid, SZ_PARAM) passdir = imgeti (im, "PASS_DIR") end