include "polarimetry.h" include include define SZ_HIST SZ_TIME+NAPER*SZ_FNAME*2+100 # size of the history log # POL_DO -- Perform HSP polarimetry reduction # # Description: # ------------ # # Date Author Description # ---- ------ ----------- # 06-Sep-1984 M. F. Hartman Original module # 29-Jun-1990 J.-C. Hsu rewrite in SPP # 19-Aug-1993 J.-C. Hsu change PANGAPER to PA_APER #------------------------------------------------------------------------------ procedure pol_do (fin, finmask, fout, nsets, mask, npts_poleff, npts_paoffset, aper_poleff, poleff, poleff_err, filter_pa, paoffset, paoffset_err, peff, pa0, peff_err, pa0_err, datafill) pointer fin, finmask, fout # input: file template pointers real datafill # input: fill value for bad output data pixel int nsets # input: number of input file sets bool mask # input: mask file existence flags int npts_poleff, npts_paoffset # input: sizes of the tables char aper_poleff[SZ_APER, npts_poleff] # input: tabulated aperture names real poleff[npts_poleff], poleff_err[npts_poleff] # input: tabulated polarization efficiencies and errors char filter_pa[SZ_FILTER, npts_paoffset] # input: tabulated filter names real paoffset[npts_paoffset], paoffset_err[npts_paoffset] # input: tabulated position angle offsets and errors real peff[NAPER], peff_err[NAPER] # input: polarization efficiencies and errors real pa0, pa0_err # input: position angle offset and error real paaper0 # position angle of 0 degree aperture pointer ipin[NAPER], ipinmask[NAPER], ipout[NOUT] pointer arrin[NAPER], arrmask[NAPER], arrout[NOUT] # addresses of data char oroot[SZ_FNAME] # output file root name char ifile[SZ_FNAME, NAPER], ofile[SZ_FNAME, NOUT] char imask[SZ_FNAME, NAPER] char aper[SZ_APER, NAPER] char filter[SZ_FILTER] char appen[1, NOUT] int indx[NAPER], first int nfile[NAPER] int ndim int npix # number of points of the input file real minval, maxval int nchar, i, j pointer ip, sp char text[SZ_HIST] pointer immap() pointer imgl1r() pointer impl1r() int strmatch(), stridx(), strlen() int imtgetim() real imgetr() long clktime() #============================================================================== begin # define output file appendage call strcpy ("i", appen[1, INTEN], 1) call strcpy ("p", appen[1, POL], 1) call strcpy ("a", appen[1, PA], 1) call strcpy ("q", appen[1, Q], 1) call strcpy ("u", appen[1, U], 1) call strcpy ("m", appen[1, MASK], 1) # loop all input files do j = 1, nsets { call smark (sp) # read the next set of input file names in the template list do i = 1, NAPER nfile[i] = 0 do i = 1, NAPER { nchar = imtgetim (fin, ifile[1, i], SZ_FNAME) # open input files ipin[i] = immap (ifile[1, i], READ_ONLY, 0) # read the aperture name and determine the orientation call imgstr (ipin[i], "APERTURE", aper[1, i], SZ_APER) # label each file according to the aperture name # this is HSP specific if (strmatch(aper[1,i], "P0") != 0) { indx[i] = P0 first = i paaper0 = imgetr (ipin[i], "PA_APER") nfile[P0] = nfile[P0] + 1 } else if (strmatch(aper[1,i], "P45") != 0) { indx[i] = P45 nfile[P45] = nfile[P45] + 1 } else if (strmatch(aper[1,i], "P90") != 0) { indx[i] = P90 nfile[P90] = nfile[P90] + 1 } else if (strmatch(aper[1,i], "P135") != 0) { indx[i] = P135 nfile[P135] = nfile[P135] + 1 } else { call pol_error ("unknown analyzer angle", i, j) } # get the polarization efficiency from the table if supplied if (npts_poleff != NULL) call get_poleff (aper[1,i], aper_poleff, poleff, poleff_err, npts_poleff, peff[indx[i]], peff_err[indx[i]]) # must be real data type if (IM_PIXTYPE(ipin[i]) != TY_REAL) call pol_error ("input file is not real data type", i, j) arrin[indx[i]] = imgl1r (ipin[i]) } # check there is only one file for each aperture do i = 1, NAPER { if (nfile[i] != 1) call pol_error ("incomplete set of apertures", 0, j) } # what is the filter name call strcat (aper[stridx("F", aper[1, 1]), 1], filter, SZ_FILTER) # check dimension, data type, size, and filter name against the # first file npix = IM_LEN(ipin[1], 1) ndim = IM_NDIM(ipin[1]) if (ndim != 1) call pol_error ("Not one-dimensional data", 1, j) do i = 2, NAPER { if (IM_LEN(ipin[i], 1) != npix) call pol_error ("mismatch of file size", i, j) if (IM_NDIM(ipin[i]) != ndim) call pol_error ("mismatch of dimension", i, j) if (strmatch(aper[1,i], filter) == 0) call pol_error ("mismatch of filter name", i, j) } # get the position angle offset from the table if supplied if (npts_paoffset != NULL) call get_paoffset (filter, filter_pa, paoffset, paoffset_err, npts_paoffset, pa0, pa0_err) # open input masks, if any. check dimension, data type, and size if (mask) { do i = 1, NAPER { nchar = imtgetim (finmask, imask[1, i], SZ_FNAME) ipinmask[i] = immap (imask[1, i], READ_ONLY, 0) if (IM_LEN(ipinmask[i], 1) != npix) call pol_error ("mismatch of input mask size", i, j) if (IM_NDIM(ipinmask[i]) != ndim) call pol_error ("mismatch of input mask dimension", i,j) if (IM_PIXTYPE(ipinmask[i]) != TY_REAL) call pol_error ("input mask is not real data type", i,j) arrmask[indx[i]] = imgl1r (ipinmask[i]) } } else { do i = 1, NAPER call salloc (arrmask[indx[i]], 1, TY_REAL) } # allocate output array space do i = 1, NOUT call malloc (arrout[i], npix, TY_REAL) # calculate polarization indices call calc_stokes (arrin, arrmask, arrout, npix, mask, peff, pa0, paaper0, datafill) # initialize the output file and mask nchar = imtgetim (fout, oroot, SZ_FNAME) do i= 1, NOUT { call strcpy (oroot, ofile[1, i], strlen(oroot)) call strcat (appen[1, i], ofile[1, i], SZ_FNAME) ipout[i] = immap (ofile[1, i], NEW_COPY, ipin[first]) IM_LEN(ipout[i], 1) = npix IM_PIXTYPE(ipout[i]) = TY_REAL # put the result to the output files and close them ip = impl1r (ipout[i]) call amovr (Memr[arrout[i]], Memr[ip], npix) } # update output files' header keywords call impstr (ipout[INTEN], "BUNIT", "FLUX") call impstr (ipout[POL], "BUNIT", "PERCENT") call impstr (ipout[PA], "BUNIT", "DEGREE") call impstr (ipout[Q], "BUNIT", "PERCENT") call impstr (ipout[U], "BUNIT", "PERCENT") call impstr (ipout[MASK], "BUNIT", "PIXEL") do i = 1, NOUT { # add new keywords about the calibration values used call imaddr (ipout[i], "PEFF0", peff[P0]) call imaddr (ipout[i], "PEFF45", peff[P45]) call imaddr (ipout[i], "PEFF90", peff[P90]) call imaddr (ipout[i], "PEFF135", peff[P135]) call imaddr (ipout[i], "PA0", pa0) } # update output files' history call cnvtime (clktime(0), text, SZ_HIST) call strcat (" created by the task POLARIMETRY, from the files: ", text, SZ_HIST) do i = 1, NAPER { call strcat (ifile[1, i], text, SZ_HIST) call strcat (" ", text, SZ_HIST) if (mask) { call strcat (imask[1, i], text, SZ_HIST) call strcat (" ", text, SZ_HIST) } } do i = 1, NOUT { call imputh (ipout[i], "HISTORY", text) # update minimum and maximum call alimr (Memr[arrout[i]], npix, minval, maxval) IM_MIN(ipout[i]) = minval IM_MAX(ipout[i]) = maxval IM_LIMTIME(ipout[i]) = IM_MTIME(ipout[i]) + 1 } # close files do i= 1, NAPER { call imunmap (ipin[i]) if (mask) call imunmap (ipinmask[i]) } do i = 1, NOUT { call imunmap (ipout[i]) # print out message of which files been created call printf ("polarimetry: output file %s is created\n") call pargstr (ofile[1, i]) call mfree (arrout[i], TY_REAL) } call sfree (sp) } end # POL_ERROR -- Issue errors specific to this task # # Description: # ------------ # # Date Author Description # ---- ------ ----------- # 04-Jul-1990 J.-C. Hsu coding #------------------------------------------------------------------------------ procedure pol_error (text, i, j) char text[SZ_LINE] int i, j #============================================================================== begin if (i > 0) { call eprintf ("%s at file %s, in file set no. %s\n") call pargstr (text) call pargi (i) call pargi (j) } else { call eprintf ("%s in file set no. %s\n") call pargstr (text) call pargi (j) } call error (1, "") end