include include define SZ_HIST SZ_TIME+3*SZ_FNAME+200 # size of the history log define OKVAL 0 # COADD_NTON -- Perform phase co-adding from N input files onto N output files # # Description: # ------------ # This is mostly a file housekeeping routine for the phase coadding task. # It takes care of input/output files open/close, check for dimensions, read/ # write data from/to files, allocate memory spaces etc. # # Date Author Description # ---- ------ ----------- # 20-Feb-1990 J.-C. Hsu rewrite in SPP #------------------------------------------------------------------------------ procedure coadd_nton (fin, fphase, finmask, fout, ferror, foutmask, nfin, nbins, maskflag, fillval) pointer fin, fphase, finmask, fout, ferror, foutmask # input: file template pointers int nfin # input: number of files in the input template int nbins # input: number of bins for coaddition bool maskflag # input: is there input mask file(s)? real fillval # input: fill value for invalid pixels in output pointer arrdata, arrphase, arrmask # addresses of data pointer arrout, arrerror, arromask # addresses of output data pointer npts, sum, sumwt, sumsq # addresses of statistics buffers pointer ipin, ipphase, ipinmask, ipout, iperror, ipoutmask char ifile[SZ_FNAME], phfile[SZ_FNAME], imask[SZ_FNAME], ofile[SZ_FNAME], errorfile[SZ_FNAME], omask[SZ_FNAME] int npix # number of points of the input file int nchar, i, j pointer sp, ip # dummy pointer char text[SZ_HIST] pointer immap() pointer imgl1r(), imgl1d(), impl1r() int imtgetim() long clktime() #============================================================================== begin # allocate array space call smark (sp) call salloc (arrout, nbins, TY_REAL) call salloc (arrerror, nbins, TY_REAL) call salloc (arromask, nbins, TY_REAL) call salloc (npts, nbins, TY_INT) call salloc (sum, nbins, TY_DOUBLE) call salloc (sumwt, nbins, TY_DOUBLE) call salloc (sumsq, nbins, TY_DOUBLE) # loop all input files do i = 1, nfin { # read the next file name in the template list nchar = imtgetim (fin, ifile, SZ_FNAME) nchar = imtgetim (fphase, phfile, SZ_FNAME) nchar = imtgetim (fout, ofile, SZ_FNAME) nchar = imtgetim (ferror, errorfile, SZ_FNAME) nchar = imtgetim (foutmask, omask, SZ_FNAME) # open input/output files ipin = immap (ifile, READ_ONLY, 0) ipphase = immap (phfile, READ_ONLY, 0) ipout = immap (ofile, NEW_FILE, 0) iperror = immap (errorfile, NEW_FILE, 0) ipoutmask = immap (omask, NEW_FILE, 0) IM_LEN(ipout,1) = nbins IM_LEN(iperror,1) = nbins IM_LEN(ipoutmask,1) = nbins # check the file size npix = IM_LEN(ipin, 1) if (IM_LEN(ipphase,1) != npix) call error (1, "mismatch of file size between input data and phase file") # 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") # input phase file must be 1-D data and double precision data type if (IM_NDIM(ipphase) != 1) call error (1, "input phase file is not 1-D data") if (IM_PIXTYPE(ipphase) != TY_DOUBLE) call error (1, "input phase file is not double data type") # read data from the input data and phase files arrdata = imgl1r (ipin) arrphase = imgl1d (ipphase) # 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") arrmask = imgl1r (ipinmask) } # reset summation buffers to zeros call amovki (0, Memi[npts], nbins) call amovkd (0., Memd[sum], nbins) call amovkd (0., Memd[sumwt], nbins) call amovkd (0., Memd[sumsq], nbins) # Do the coadding call phase_coadd (Memd[arrphase], Memr[arrdata], Memr[arrmask], maskflag, npix, nbins, Memi[npts], Memd[sum], Memd[sumwt], Memd[sumsq]) # calculate the mean and mean error do j = 0, nbins-1 { if (Memi[npts+j] > 1) { Memr[arromask+j] = OKVAL Memr[arrout+j] = real(Memd[sum+j]/Memd[sumwt+j]) Memr[arrerror+j] = sqrt((real(Memd[sumsq+j]/Memd[sumwt+j]) - Memr[arrout+j]**2)/real(Memi[npts+j]-1)) } else if (Memi[npts+j] == 1) { Memr[arromask+j] = OKVAL Memr[arrout+j] = real(Memd[sum+j]/Memd[sumwt+j]) Memr[arrerror+j] = 0. } else { Memr[arromask+j] = fillval Memr[arrout+j] = 0. Memr[arrerror+j] = 0. } } # put the result to the output files and close them ip = impl1r (ipout) call amovr (Memr[arrout], Memr[ip], nbins) ip = impl1r (iperror) call amovr (Memr[arrerror], Memr[ip], nbins) ip = impl1r (ipoutmask) call amovr (Memr[arromask], Memr[ip], nbins) # update header keyword(s) and history call imastr (ipout, "BUNIT", "PHASE", SZ_BUNIT) call cnvtime (clktime(0), text, SZ_HIST) call strcat (" created by the task COADD, from the science file: ", text, SZ_HIST) call strcat (ifile, text, SZ_HIST) call strcat (", and the phase file: ", text, SZ_HIST) call strcat (phfile, text, SZ_HIST) call imputh (ipout, "HISTORY", text) call imputh (iperror, "HISTORY", text) call imputh (ipoutmask, "HISTORY", text) # close files call imunmap (ipin) call imunmap (ipphase) call imunmap (ipout) call imunmap (iperror) call imunmap (ipoutmask) if (maskflag) call imunmap (ipinmask) # print out message of which files been created call printf ("coadd: output file %s is created\n") call pargstr (ofile) call printf ("coadd: output error file %s is created\n") call pargstr (errorfile) call printf ("coadd: output mask %s is created\n") call pargstr (omask) } call sfree (sp) end