include include define SZ_FITTYPE 8 define SZ_HIST SZ_TIME+3*SZ_FNAME+100 # SUBSKY_DO -- Perform backgrouund subtraction # # Description: # ------------ # # Date Author Description # ---- ------ ----------- # 23-Mar-1990 J.-C. Hsu rewrite in SPP #------------------------------------------------------------------------------ procedure subsky_do (fin, fsky, fout, finmask, fskymask, foutmask, nfin, fittype, order, nptsfit, nskyfiles, maskflag, bmaskflag) pointer fin, fsky, fout, finmask, fskymask, foutmask # input: file template pointers int nfin # input: number of files in the input template char fittype[SZ_FITTYPE] # input: fitting method int order # input: order of the fitting polynomial int nptsfit # input: number of fitting points int nskyfiles # input: number of sky files per data file bool maskflag # input: is there input mask file(s)? bool bmaskflag # input: is there sky mask file(s)? pointer indata, inmask, sky, skymask, outdata, outmask # addresses of data pointer intime, skytime pointer ipin, ipout, ipinmask, ipoutmask, ipsky1, ipsky2, ipsmask1, ipsmask2 char ifile[SZ_FNAME], ofile[SZ_FNAME], imask[SZ_FNAME], omask[SZ_FNAME] char skyfile1[SZ_FNAME], skyfile2[SZ_FNAME], smask1[SZ_FNAME] char smask2[SZ_FNAME] int npix # number of points in the time series int npsky1, npsky2 # number of points in sky file(s) int nchar, i pointer sp, ip # dummy pointer char text[SZ_HIST] pointer immap() pointer imgl1r(), impl1r() int imtgetim() bool streq() long clktime() #============================================================================== begin # loop all input files do i = 1, nfin { npsky2 = 0 # read the next file name in the template list nchar = imtgetim (fin, ifile, SZ_FNAME) nchar = imtgetim (fout, ofile, SZ_FNAME) nchar = imtgetim (foutmask, omask, SZ_FNAME) nchar = imtgetim (fsky, skyfile1, SZ_FNAME) # open input, output, and sky files ipin = immap (ifile, READ_ONLY, 0) ipout = immap (ofile, NEW_COPY, ipin) ipoutmask = immap (omask, NEW_COPY, ipin) ipsky1 = immap (skyfile1, READ_ONLY, 0) # if there are two sky files per science file, open the second # sky file if (nskyfiles == 2) { nchar = imtgetim (fsky, skyfile2, SZ_FNAME) ipsky2 = immap (skyfile2, READ_ONLY, 0) npsky2 = IM_LEN(ipsky2, 1) } npix = IM_LEN(ipin, 1) npsky1 = IM_LEN(ipsky1, 1) # input file must be 1-D data and real data type if (IM_NDIM(ipin) != 1) call error (1, "input is not 1-D data") if (IM_PIXTYPE(ipin) != TY_REAL) call error (1, "input file is not real data type") # allocate array space call smark (sp) call salloc (outdata, npix, TY_REAL) call salloc (outmask, npix, TY_REAL) call salloc (intime, npix, TY_DOUBLE) call salloc (sky, npsky1+npsky2, TY_REAL) call salloc (skytime, npsky1+npsky2, TY_DOUBLE) # read data from the input file and sky file(s) indata = imgl1r (ipin) ip = imgl1r (ipsky1) call amovr (Memr[ip], Memr[sky], npsky1) if (nskyfiles == 2) { ip = imgl1r (ipsky2) call amovr (Memr[ip], Memr[sky+npsky1], npsky2) } # if there is input mask, read its data if (maskflag) { nchar = imtgetim (finmask, imask, SZ_FNAME) ipinmask = immap (imask, READ_ONLY, 0) if (IM_LEN(ipinmask,1) != npix) call error (1, "mismatch of file size between input and input mask") inmask = imgl1r (ipinmask) } # if there is sky mask, read its data if (bmaskflag) { call salloc (skymask, npsky1+npsky2, TY_REAL) nchar = imtgetim (fskymask, smask1, SZ_FNAME) ipsmask1 = immap (smask1, READ_ONLY, 0) if (IM_LEN(ipsmask1,1) != npsky1) call error (1, "mismatch of file size between sky and sky mask") if (nskyfiles == 2) { nchar = imtgetim (fskymask, smask2, SZ_FNAME) ipsmask2 = immap (smask2, READ_ONLY, 0) if (IM_LEN(ipsmask2,1) != npsky2) call error (1, "mismatch of file size between sky and sky mask") } ip = imgl1r (ipsmask1) call amovr (Memr[ip], Memr[skymask], npsky1) if (nskyfiles == 2) { ip = imgl1r (ipsmask2) call amovr (Memr[ip], Memr[skymask+npsky1], npsky2) } } # get the time array call get_times (ipin, ipsky1, ipsky2, nskyfiles, Memd[intime], Memd[skytime], npix, npsky1, npsky2) # overall fitting case if (streq(fittype, "overall")) call subglobal(Memr[indata], Memd[intime], Memr[sky], Memd[skytime], Memr[inmask], Memr[skymask], Memr[outdata], Memr[outmask], order, maskflag, bmaskflag, npix, npsky1+npsky2) # other cases else call fitnsub (Memr[indata], Memd[intime], Memr[sky], Memd[skytime], Memr[inmask], Memr[skymask], Memr[outdata], Memr[outmask], fittype, nptsfit, order, maskflag, bmaskflag, npix, npsky1+npsky2) # put the result to the output files and close them ip = impl1r (ipout) call amovr (Memr[outdata], Memr[ip], npix) ip = impl1r (ipoutmask) call amovr (Memr[outmask], Memr[ip], npix) # update header keyword(s) and history call cnvtime (clktime(0), text, SZ_HIST) call strcat (" created by the task SUBSKY, ", text, SZ_HIST) call strcat ("from the science file: ", text, SZ_HIST) call strcat (ifile, text, SZ_HIST) call strcat (", and the sky file(s): ", text, SZ_HIST) call strcat (skyfile1, text, SZ_HIST) if (nskyfiles == 2) { call strcat (", ", text, SZ_HIST) call strcat (skyfile2, text, SZ_HIST) } call imputh (ipout, "HISTORY", text) call imputh (ipoutmask, "HISTORY", text) # close files call imunmap (ipin) call imunmap (ipout) call imunmap (ipoutmask) call imunmap (ipsky1) if (nskyfiles == 2) call imunmap (ipsky2) if (maskflag) call imunmap (ipinmask) if (bmaskflag) { call imunmap (ipsmask1) if (nskyfiles == 2) call imunmap (ipsmask2) } # print out message of which files been created call printf ("subsky: output file %s is created\n") call pargstr (ofile) call printf ("subsky: output mask %s is created\n") call pargstr (omask) call sfree (sp) } end