include include define OKVAL 0 # COADD_NTO1 -- Perform phase co-adding from N input files onto one output file # # 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_nto1 (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 int hsize pointer sp, ip # dummy pointer pointer text pointer immap() pointer imgl1r(), impl1r(), imgl1d() int imtgetim() long clktime() #============================================================================== begin # find out how long the history log may be hsize = SZ_TIME + (2 * SZ_FNAME + 30) * nfin + 50 # 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) call salloc (text, hsize, TY_CHAR) # get names of the output files nchar = imtgetim (fout, ofile, SZ_FNAME) nchar = imtgetim (ferror, errorfile, SZ_FNAME) nchar = imtgetim (foutmask, omask, SZ_FNAME) # open output files and specify their sizes 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 # 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) call cnvtime (clktime(0), Memc[text], hsize) call strcat (" created by the task COADD, from: ", Memc[text], hsize) # 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) # open input/output files ipin = immap (ifile, READ_ONLY, 0) ipphase = immap (phfile, READ_ONLY, 0) # 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) } # Do the coadding call phase_coadd (Memd[arrphase], Memr[arrdata], Memr[arrmask], maskflag, npix, nbins, Memi[npts], Memd[sum], Memd[sumwt], Memd[sumsq]) call strcat (" ", Memc[text], hsize) call strcat (ifile, Memc[text], hsize) call strcat (", and the phase file ", Memc[text], hsize) call strcat (phfile, Memc[text], hsize) call strcat (";", Memc[text], hsize) # close input files call imunmap (ipin) call imunmap (ipphase) if (maskflag) call imunmap (ipinmask) } # 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 imputh (ipout, "HISTORY", Memc[text]) call imputh (iperror, "HISTORY", Memc[text]) call imputh (ipoutmask, "HISTORY", Memc[text]) # close output files call imunmap (ipout) call imunmap (iperror) call imunmap (ipoutmask) # 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