procedure asncombine(association) file association {"", prompt="list of images", mode="q"} file mask {"nictools$mask.pl", prompt="bad pixel mask"} bool doerrors {yes, prompt="do error computation?"} file background{"", prompt= "background to subtract from each image before shift"} string keyback {"", prompt= "background keyword"} bool fheader {no, prompt= "use background scaling from header if available?"} real fmin {0.5, prompt= "mimimum factor to multiply background"} real fmax {1.5, prompt= "maximum factor to multiply background"} int nf {11, prompt= "number of factors to try"} bool bckdisplay { yes, prompt= "display different background subtracted images?"} bool regtable {yes, prompt= "use regions from table if available?"} string region1 {"[1:128,*]", prompt="first region on image to statistics on"} string region2 {"[129:256,*]", prompt="second region on image to statistics on"} bool submedian {yes, prompt="subtract the mode from each column?"} bool modbefore {no, prompt= "use previously deteremined modes if avaiable?"} bool modreg {yes, prompt="use rows in region 1/2 for column medians?"} string bckregion {"*", prompt="rows (format n1:n2) to compute mode for each column"} string modcol {"gauss", prompt= "what fit to modes should be used?", \ enum="imsurfit|subtracted|gauss|median"} bool moddisplay { yes, prompt= "display modes subtracted from the images?"} bool subquad {yes, prompt="try to balance quadrants?"} bool qheader {no, prompt= "use quadrants stats from header if available?"} string regquad1 {"[1:128,129:256]", prompt="upper left quadrant for balancing"} string regquad2 {"[129:256,129:256]", prompt="upper right quadrant for balancing"} string regquad3 {"[1:128,1:128]", prompt="lower left quadrant for balancing"} string regquad4 {"[129:256,1:128]", prompt="lower right quadrant for balancing"} bool quaddisplay { yes, prompt= "display quadrants subtracted from the images?"} int resolution {4, prompt= "resolution in fraction of pixel"} int rebin {4, prompt= "rebinning after processing (1=keep full resolution)"} int wcsref {1, prompt= "reference image for wcs (< 1 if none)"} file logfile {"asncombine.log", prompt= "log file with shifts and factors"} int nsigma {5, prompt= "nsigma rejection for stats ( <0 to do no rejection)"} int nleft {1, prompt= "order of background fit for x < 129 "} int nright {3, prompt= "order of background fit for x > 128 "} #string xreg1 {"129-200", prompt= "first region for linear fits"} #string xreg2 {"205-225", prompt= "first region for linear fits"} #string xreg3 {"230-256", prompt= "first region for linear fits"} real smlength {5, prompt= "smoothing length for medians"} bool split {yes, prompt="combine 2 subset of images?"} bool hupdate {yes, prompt="update header with results?"} string weightregion {"[10:246,10:246]", prompt="region on subsampled image to compute weight"} file subdir {"subtracted", prompt= "optional directory to save subtracted images"} int inssize {-1, prompt="size of image to imbed data (<1 for none)"} bool overwrite {no, prompt= "overwrite existing output images?"} bool stopplot {no, prompt="stop after each plot?"} bool restart {no, prompt="restart of a previous run?"} bool resdisp {yes, prompt="display final result?"} bool continue {yes, prompt="continue?", mode = "q"} bool verbose {no, prompt="verbose output for debugging?"} struct *imagelist begin file comb, names, masks, names_tmp, names_tmp2, masks_tmp, outname file comb1, comb2, tmp_table, tmp_table1, tmp_table2, tmp_table3 file names_tmperr2, names_tmperr, tabledir file comb_err, backused, masked_bck file backtry, backlist, nbest, member, igi_data, igi_script, line_tmp file backtry2, listpix, masklist file names_mask, names_mask2, insimage, insmask string backregion1, backregion2, insregion int xoff, yoff, l, len, kf, i, nmembers, ins1, ins2 real factor, fbest, sigmabest, sigma, med, sigma1, sigma2, npix1, npix2 real sigmin, sigmax, facmin, facmax, newval, max_size real mod1bck, mod2bck, mod3bck, mod4bck real fbest1, mode1, mode2, low, up, mod1, mod2, mod3, mod4 int ix1, ix2, iy1, iy2 file bck_tmp, bck_fit, names_tmp3, bck_fit1, bck_fit2, memnam, names_tmp4 file fit_image, fit_image2, fit_image3, fit_image4, fit_image5 bool nomod, oldtable cache ("tinfo") cache ("tabpar") cache ("keypar") max_size= 1000*resolution # maximum size of shift (any big integer number) fbest= 0 outname= association len= strlen(association) if (substr(association, len-4, len) == "_asn" ) { outname= substr(association, 1, len-4) } if (substr(association, len-4, len) == "_asc" ) { outname= substr(association, 1, len-4) } if (substr(association, len-8, len) == "_asc.fits" ) { outname= substr(association, 1, len-9) } if (substr(association, len-8, len) == "_asn.fits" ) { outname= substr(association, 1, len-9) } if (access(outname//".fits")) { if (overwrite) del(outname//".fits", veri-) else { print("output file ", outname, " already exists.") print(" use asncombine overwrite=yes") bye } } print( "outname=", outname) set tmp = "tmp/" if (!access("tmp")) mkdir("tmp") listpix= mktemp("tmp$listpix") insimage= mktemp("tmp$insimg")//".fits" insmask= mktemp("tmp$insmsk")//".pl" names_mask= mktemp("tmp$nm_msk")//".pl" names_mask2= mktemp("tmp$nm_msk2_")//".pl" masklist= mktemp("tmp$masklist") names= mktemp("tmp$rep") masks= mktemp("tmp$msk") masked_bck= mktemp("tmp$mskbck")//".fits" names_tmp= mktemp("tmp$rep_tmp")//".fits" names_tmperr= mktemp("tmp$tmp_err")//".fits" names_tmperr2= mktemp("tmp$tmp_err2_")//".fits" names_tmp2= mktemp("tmp$rep_tm2")//".fits" names_tmp3= mktemp("tmp$rep_tm3")//".fits" names_tmp4= mktemp("tmp$rep_tm4")//".fits" bck_tmp= mktemp("tmp$bck")//".fits" bck_fit= mktemp("tmp$bckfit")//".fits" bck_fit1= mktemp("tmp$bckfit1_")//".fits" bck_fit2= mktemp("tmp$bckfit2_")//".fits" masks_tmp= mktemp("tmp$msk_tmp")//".pl" comb= mktemp("tmp$com")//".list" comb_err= mktemp("tmp$com_err")//".list" comb1= mktemp("tmp$com")//".list" comb2= mktemp("tmp$com")//".list" tmp_table= mktemp("tmp$table")//".fits" tmp_table1= mktemp("tmp$table1_")//".fits" tmp_table2= mktemp("tmp$table2_")//".fits" tmp_table3= mktemp("tmp$table3_")//".fits" fit_image= mktemp("tmp$fitim")//".fits" fit_image2= mktemp("tmp$fitim2_")//".fits" fit_image3= mktemp("tmp$fitim3_")//".fits" fit_image4= mktemp("tmp$fitim4_")//".fits" fit_image5= mktemp("tmp$fitim5_")//".fits" igi_data= mktemp("tmp$igid") igi_script= mktemp("tmp$script") if (subdir != "" && !access(subdir)) mkdir(subdir) # if (inssize > 0) { mkpattern(insimage, ndim=2, ncol=inssize, nline=inssize,\ v1=1, v2=1, pattern="constant") ins1= inssize/2-128 ins2= ins1+256 insregion= "["//ins1//":"//ins2//","//ins1//":"//ins2//"]" if (verbose) print("region to insert:", insregion) imcopy(mask, insimage//insregion, verb= verbose) imcopy(insimage, insmask, verb= verbose) mkpattern(insimage, ndim=2, ncol=inssize, nline=inssize,\ v1=0, v2=0, pattern="constant", option="replace") } # backregion1= bckregion backregion2= bckregion # # BIG loop: for each member in asn file: # tinfo(outname//"_asn.fits", > "dev$null") nmembers= tinfo.nrows-1 print("number of images to combine: ", nmembers) for (l=1; l< tinfo.nrows; l+=1) { tabpar(outname//"_asn.fits", "MEMNAME", l ) member= tabpar.value//"_cal" memnam= tabpar.value # # restart after a crash # if (restart) print("looking for ", subdir//"/"//memnam//"_cal.fits") if (restart && access(subdir//"/"//memnam//"_cal.fits")) { print("reusing ", subdir//"/"//memnam//"_cal.fits") copy(subdir//"/"//memnam//"_cal.pl", ".") if (l == 1) { print(subdir//"/"//memnam//"_cal.fits", > comb) if (doerrors) print(subdir//"/"//memnam//"_err.fits", > comb_err) } else { print(subdir//"/"//memnam//"_cal.fits", >> comb) if (doerrors) print(subdir//"/"//memnam//"_err.fits", >> comb_err) } if (nmembers > 1 && split && (l <= tinfo.nrows/2.) ) \ print(subdir//"/"//memnam//"_cal.fits", >> comb1) if (nmembers > 1 && split && (l > tinfo.nrows/2.) ) \ print(subdir//"/"//memnam//"_cal.fits", >> comb2) next } # # # print("processing ", member) if ( regtable || modreg ) { keypar(association//"[1]", "BCKREG1", silent=verbose) if (!keypar.found) keypar(association//"[0]", "BCKREG1", silent=verbose) if (keypar.found) { if (verbose) print(keypar.value) print(keypar.value) | scan(ix1, ix2, iy1, iy2) tabpar(outname//"_asn.fits", "XOFFSET", l ) ix1= int(ix1+real(tabpar.value)) ix2= int(ix2+real(tabpar.value)) tabpar(outname//"_asn.fits", "YOFFSET", l ) iy1= int(iy1+real(tabpar.value)) iy2= int(iy2+real(tabpar.value)) if (iy1 < 1) iy1=1; if (iy1 > 256) iy1= 256 if (ix1 < 1) ix1=1; if (ix1 > 256) ix1= 256 if (iy2 < 1) iy2=1; if (iy2 > 256) iy2= 256 if (ix2 < 1) ix2=1; if (ix2 > 256) ix2= 256 region1= "["//ix1//":"//ix2//","//iy1//":"//iy2//"]" if (verbose) print(region1) backregion1= str(int(min(iy1, iy2)))//":"//int(max(iy1, iy2)) } keypar(association//"[1]", "BCKREG2", silent=verbose) if (!keypar.found) keypar(association//"[0]", "BCKREG2", silent=verbose) if (keypar.found) { if (verbose) print(keypar.value) print(keypar.value) | scan(ix1, ix2, iy1, iy2) tabpar(outname//"_asn.fits", "XOFFSET", l ) ix1= int(ix1+real(tabpar.value)) ix2= int(ix2+real(tabpar.value)) tabpar(outname//"_asn.fits", "YOFFSET", l ) iy1= int(iy1+real(tabpar.value)) iy2= int(iy2+real(tabpar.value)) if (iy1 < 1) iy1=1; if (iy1 > 256) iy1= 256 if (ix1 < 1) ix1=1; if (ix1 > 256) ix1= 256 if (iy2 < 1) iy2=1; if (iy2 > 256) iy2= 256 if (ix2 < 1) ix2=1; if (ix2 > 256) ix2= 256 region2= "["//ix1//":"//ix2//","//iy1//":"//iy2//"]" backregion2= str(int(min(iy1, iy2)))//":"//int(max(iy1, iy2)) if (verbose) print(region2) } } if (background == "" && keyback == "") imcopy(member//".fits[sci]", names_tmp3, verb=verbose) else { if (keyback == "") backused= background else { imgets(member//".fits[0]", "bckfile") backused= imgets.value } if (!access(masked_bck) ) { imcopy(backused, masked_bck, verbose=verbose) if (verbose) print("masking ", masked_bck, " with ", mask) fixpix(masked_bck, mask) } sigmabest=9999999 sigmin=9999999 sigmax=-9999999 facmin=9999999 facmax=-9999999 backlist= mktemp("tmp$bcklist") print("...subtracting background ", backused) imgets.value="9999" if (fheader) { imgets( member//".fits[0]", "factor", > "dev$null") } if (imgets.value != "0" && fheader) fbest= real(imgets.value) else { # # # now let's try different background scalings: for (kf=1; kf <= nf; kf+=1) { backtry= mktemp("tmp$bcktry")//".fits" backtry2= mktemp("tmp$bcktry")//".fits" print(backtry, >> backlist) if (nf > 1) factor= fmin + (kf-1)* (fmax-fmin)/(nf-1) else factor= (fmin + fmax) /2. #imcalc(tabpar.value//"_cal.fits[sci],"//backused, \ # backtry, "im1-"//factor//"*im2", verb-) imari(backused, "*", factor, backtry2) imari(memnam//"_cal.fits[sci]", "-", backtry2, backtry) del(backtry2, veri-) if (verbose) print("masking ", backtry, " with ", mask) fixpix(backtry, mask) imstat(backtry//region1, form-, field="stddev") | \ scan(sigma1) imstat(backtry//region2, form-, field="stddev") | \ scan(sigma2) imstat(backtry//region1, form-, field="npix") | \ scan(npix1) imstat(backtry//region2, form-, field="npix") | \ scan(npix2) if (nsigma > 0) { imstat(backtry//region1, form-, field="mode") | \ scan(mode1) imstat(backtry//region2, form-, field="mode") | \ scan(mode2) up= (mode1+mode2+nsigma*(sigma1+sigma2))/2 low= (mode1+mode2-nsigma*(sigma1+sigma2))/2 imstat(backtry//region1, form-, field="stddev", \ low=low, up=up ) | scan(sigma1) imstat(backtry//region2, form-, field="stddev", \ low=low, up=up ) | scan(sigma2) } sigma= sqrt( sigma1**2*npix1 + sigma2**2*npix2) / \ sqrt(npix1 + npix2) print(factor, " ", sigma) print(factor, " ", sigma, " ", sigma1, " ", sigma2, >> igi_data) if (sigma < sigmabest ) { sigmabest= sigma fbest= factor nbest= backtry } sigmin= min(sigmin, sigma, sigma1, sigma2) sigmax= max(sigmax, sigma, sigma1, sigma2) facmin= min(facmin, factor) facmax= max(facmax, factor) } # now let's try different background scalings: fbest1= fbest if (nf > 1) { for (kf=0; kf <= 6; kf+=1) { if (kf == 3) next backtry= mktemp("tmp$bcktry")//".fits" backtry2= mktemp("tmp$bcktry2")//".fits" print(backtry, >> backlist) factor= fbest1 + ((kf-1)/4.-0.5)* (fmax-fmin)/(nf-1) #imcalc(tabpar.value//"_cal.fits[sci],"//backused, \ # backtry2, "im1-"//factor//"*im2", verb-) imari(backused, "*", factor, backtry2) imari(memnam//"_cal.fits[sci]", "-", backtry2, backtry) del(backtry2, veri-) if (verbose) print("masking ", backtry, " with ", mask) fixpix(backtry, mask) imstat(backtry//region1, form-, field="stddev") | \ scan(sigma1) imstat(backtry//region2, form-, field="stddev") | \ scan(sigma2) imstat(backtry//region1, form-, field="npix") | \ scan(npix1) imstat(backtry//region2, form-, field="npix") | \ scan(npix2) if (nsigma > 0) { imstat(backtry//region1, form-, field="mode") | \ scan(mode1) imstat(backtry//region2, form-, field="mode") | \ scan(mode2) up= (mode1+mode2+nsigma*(sigma1+sigma2))/2 low= (mode1+mode2-nsigma*(sigma1+sigma2))/2 imstat(backtry//region1, form-, field="stddev", \ low=low, up=up ) | scan(sigma1) imstat(backtry//region2, form-, field="stddev", \ low=low, up=up ) | scan(sigma2) } sigma= sqrt( sigma1**2*npix1 + sigma2**2*npix2) / \ sqrt(npix1 + npix2) print(factor, " ", sigma) print(factor, " ", sigma, " ", sigma1, " ", sigma2, >> igi_data) if (sigma < sigmabest ) { sigmabest= sigma fbest= factor nbest= backtry } } } } # save best factor and regions in image header if (hupdate) { hedit(member//".fits[0]", "factor", fbest, add+, veri-, upd+, show-) hedit(member//".fits[0]", "region1", region1, add+, veri-, upd+, show-) hedit(member//".fits[0]", "region2", region2, add+, veri-, upd+, show-) } # display the results if desired: if (bckdisplay && !fheader) { print( "data ", igi_data, >> igi_script) print( "xc 1", >> igi_script) print( "yc 2", >> igi_script) sigmin= sigmin- (sigmax-sigmin)/10. sigmax= sigmax+ (sigmax-sigmin)/10. print( "lim ", facmin, " ", facmax, " ", \ sigmin, " ", sigmax, >> igi_script) print( "points", >> igi_script) print( "ptype 1 1", >> igi_script) print("color 2", >> igi_script) print( "yc 3 ", >> igi_script) print( "points", >> igi_script) print("color 3", >> igi_script) print( "yc 4 ", >> igi_script) print( "points", >> igi_script) print("color 1", >> igi_script) print( "box", >> igi_script) print( "xlabel factor", >> igi_script) print( "ylabel sigma", >> igi_script) print( "rel ", fbest, " ", 1.05*sigmabest, >> igi_script) print( "draw ", fbest, " ", 1.51*sigmabest, >> igi_script) print( "rel ", fbest, " ", 1.05*sigmabest, >> igi_script) print( "draw ", (fbest-0.02), " ", 1.07*sigmabest, >> igi_script) print( "rel ", fbest, " ", 1.05*sigmabest, >> igi_script) print( "draw ", (fbest+0.02), " ", 1.07*sigmabest, >> igi_script) print( "end ", >> igi_script) igi( device="stdgraph", metacode="", < igi_script ) gflush del(igi_script, veri-) mosaic(memnam//"_cal.fits[sci],@"//backlist//","//nbest) } if (!fheader) { if (verbose) { print("deleting:", backlist, " ", igi_data) type(backlist) } del("@"//backlist, veri-) del(backlist, veri-) del(igi_data, veri-) } print("using factor= ", fbest) #imcalc(tabpar.value//"_cal.fits[sci],"//backused, \ # names_tmp3, "im1-"//fbest//"*im2", verb-) imari(backused, "*", fbest, names_tmp4) #imari(tabpar.value//"_cal.fits[sci]", "-", names_tmp3, backtry) imari(memnam//"_cal.fits[sci]", "-", names_tmp4, names_tmp3) del(names_tmp4, veri-) } # balance quadrants if desired: if (subquad) { #imcopy(names_tmp3, bck_fit, > "dev$null") imcopy(names_tmp3, bck_fit, verbose=verbose) imrename(names_tmp3, names_tmp2, > "dev$null") imgets( member//".fits[0]", "quad1", > "dev$null") if (qheader) { imgets( member//".fits[0]", "quad1", > "dev$null") mod1= real(imgets.value) imgets( member//".fits[0]", "quad2", > "dev$null") mod2= real(imgets.value) imgets( member//".fits[0]", "quad3", > "dev$null") mod3= real(imgets.value) imgets( member//".fits[0]", "quad4", > "dev$null") mod4= real(imgets.value) mod1bck= 1 mod2bck= 1 mod3bck= 1 mod4bck= 1 } else { if (verbose) print("masking ", names_tmp2, " with ", mask) fixpix(names_tmp2, mask) imstat(names_tmp2//regquad1, form-, field="midpt") \ | scan(mod1) imstat(names_tmp2//regquad2, form-, field="midpt") \ | scan(mod2) imstat(names_tmp2//regquad3, form-, field="midpt") \ | scan(mod3) imstat(names_tmp2//regquad4, form-, field="midpt") \ | scan(mod4) if (access(masked_bck) ) { imstat(masked_bck//regquad1, form-, field="midpt") \ | scan(mod1bck) imstat(masked_bck//regquad2, form-, field="midpt") \ | scan(mod2bck) imstat(masked_bck//regquad3, form-, field="midpt") \ | scan(mod3bck) imstat(masked_bck//regquad4, form-, field="midpt") \ | scan(mod4bck) del(masked_bck, veri-) } else { mod1bck= 1 mod2bck= 1 mod3bck= 1 mod4bck= 1 } } print("quadrants: ", mod1, mod2, mod3, mod4) print(" / bck : ", (mod1bck-mod1)/mod1bck, (mod2bck-mod2)/mod2bck, \ (mod3bck-mod3)/mod3bck, (mod4bck-mod4)/mod4bck ) if (verbose) print(" ... replacing q1") imreplace(bck_fit//"[1:128,1:128]", (mod3-mod1)) if (verbose) print(" ... replacing q2") imreplace(bck_fit//"[129:256,1:128]", (mod4-mod2)) if (verbose) print(" ... replacing q3") imreplace(bck_fit//"[*,129:256]", 0) imari(names_tmp2, "-", bck_fit, names_tmp3) if (stopplot) { continue= yes print("\n\nimage display shows the ", l, "th ", \ "image in association ", member, \ " with different ", \ "scalings for the subtracted background.",\ "The last displayed image is the one ",\ "which will be used.\n\n") if (!continue) bye } if (quaddisplay) { mosaic(names_tmp2//","//bck_fit//","//names_tmp3) if (stopplot) { continue= yes print("\n\nimage display shows the ", l, "th ", \ "image in association before and after ", \ "balancing of the quadrants.\n",\ "upper left: before balancing\n", "upper right: values subtracted\n", "lower left: after balancing\n\n") if (!continue) bye } } del(bck_fit, veri-, > "dev$null") del(names_tmp2, veri-, > "dev$null") if (hupdate) { hedit(member//".fits[0]", "quad1", mod1, add+, veri-, upd+, show-) hedit(member//".fits[0]", "quad2", mod2, add+, veri-, upd+, show-) hedit(member//".fits[0]", "quad3", mod3, add+, veri-, upd+, show-) hedit(member//".fits[0]", "quad4", mod4, add+, veri-, upd+, show-) hedit(member//".fits[0]", "regquad1", regquad1, add+, veri-, \ upd+, show-) hedit(member//".fits[0]", "regquad2", regquad2, add+, veri-, \ upd+, show-) hedit(member//".fits[0]", "regquad3", regquad3, add+, veri-, \ upd+, show-) hedit(member//".fits[0]", "regquad4", regquad4, add+, veri-, \ upd+, show-) } } # subtract a fit to background if (submedian) { if (verbose) print("subtracting median") if (modbefore && access(subdir//"/"//memnam//"_mod.fits")) { print("using previous file") tabim(subdir//"/"//memnam//"_mod.fits", bck_fit1, \ "modes", 1, 256) imstack(bck_fit1, bck_fit2) if (verbose) imhead(bck_fit2) blkrep(bck_fit2, bck_tmp, 1, 5) del(bck_fit1, veri-, > "dev$null") del(bck_fit2, veri-, > "dev$null") } else { # ok, we actually have to determined the modes print("taking mode column by column, please be patient") imcopy( names_tmp3, bck_tmp, verbose=verbose) if (verbose) print("masking ", names_tmp3, " with ", mask) fixpix(names_tmp3, mask) for (i=1; i <= 256; i+=1) { if (i == 128) print(" half way through") med= 0 if ( (modreg==yes) && (regtable)) { #if (verbose) print("regions: ", backregion1, " ", backregion2) listpix(names_tmp3//"["//i//","//backregion1//"]", > listpix) listpix(names_tmp3//"["//i//","//backregion2//"]", >> listpix) } else listpix(names_tmp3//"["//i//","//bckregion//"]", > listpix) tstat(listpix, "c2", > "dev$null") med= tstat.median #if (verbose) print("median 1:", med) up= med+ nsigma*tstat.stddev low= med- nsigma*tstat.stddev tstat(listpix, "c2", > "dev$null", low=low, high=up) med= tstat.median up= med+ nsigma*tstat.stddev low= med- nsigma*tstat.stddev tstat(listpix, "c2", > "dev$null", low=low, high=up) #if (verbose) print("median 2 and 3:", med, " ", tstat.median) del(listpix, veri-) imreplace(bck_tmp//"["//i//",1:5]", tstat.median) } } # imsurfitting modes print("imsurfitting modes...") print(" left...") imsurfit(bck_tmp//"[1:128,1:5]", bck_fit1, nleft, 1, cross-, type_out="fit",\ regions="all") print(" right...") imsurfit(bck_tmp//"[129:256,1:5]", bck_fit2, nright, 1, cross-,type_out="fit",\ regions="all") if (verbose) { imhead(bck_tmp) imstat(bck_tmp) #imari(bck_tmp//"[*,1:5]", "*", 0., bck_fit) mkpattern(bck_fit, pattern="constant", v1=0, v2=0,\ ndim=2, ncols=256, nlines=5) imhead(bck_fit) imstat(bck_fit) imcopy(bck_fit1, bck_fit//"[1:128,1:5]") imcopy(bck_fit2, bck_fit//"[129:256,1:5]") imhead(bck_fit) imstat(bck_fit) } else { mkpattern(bck_fit, pattern="constant", v1=0, v2=0,\ ndim=2, ncols=256, nlines=5) flpr imcopy(bck_fit1, bck_fit//"[1:128,1:5]", verb=verbose) imcopy(bck_fit2, bck_fit//"[129:256,1:5]", verb=verbose) } del(bck_fit1, veri-, > "dev$null") del(bck_fit2, veri-, > "dev$null") # save the fits if subdir is selected, otherwise use tmp direc: if (subdir == "") tabledir= "tmp$" else tabledir= subdir # create output table if (modbefore && access(tabledir//"/"//memnam//"_mod.fits")) { tproject(tabledir//"/"//memnam//"_mod.fits", \ tmp_table, "modes,col1", uniq-) del(tabledir//"/"//memnam//"_mod.fits", veri-) oldtable= yes } else { if (access(tabledir//"/"//memnam//"_mod.fits")) \ del(tabledir//"/"//memnam//"_mod.fits", veri-) imtab(bck_tmp//"[*,1]", \ tmp_table, "modes", format="", pname="col", wcs="logical") oldtable= no } imtab(bck_fit//"[*,1]", \ tmp_table, \ "imsurfit", format="") if (!oldtable && modcol == "subtracted") { print("no column previously subtracted, using imsurfit") tcalc(tmp_table, "subtracted", "imsurfit") } del(bck_fit, veri-, > "dev$null") del(bck_tmp, veri-, > "dev$null") # second method is linear fits of regions: rename(tmp_table, tabledir//"/"//memnam//"_mod.fits") # third method is smoothing: if (verbose) print("gauss") tabim( tabledir//"/"//memnam//"_mod.fits", \ fit_image, "modes", 1, 256 ) imcopy(fit_image//"[1:128]", fit_image4, verbose=verbose) gauss(fit_image4, fit_image2, smlength, ratio=1, theta=0) del(fit_image4, veri-, > "dev$null") imcopy(fit_image//"[129:256]", fit_image4, verb=verbose) gauss(fit_image4, fit_image3, smlength, ratio=1, theta=0) del(fit_image4, veri-, > "dev$null") flpr copy(fit_image, fit_image4, verb=verbose) if (verbose) imhead(fit_image) if (verbose) imhead(fit_image4) imcopy(fit_image2, fit_image4//"[1:128]", verb=verbose) imcopy(fit_image3, fit_image4//"[129:256]", verb=verbose) imtab(fit_image4, tmp_table2, "gauss", \ format="", pname="col", wcs="logical") tmerge(tabledir//"/"//memnam//"_mod.fits,"//tmp_table2, \ tmp_table3, "merge", allcol+, allrows=256) delete(tabledir//"/"//memnam//"_mod.fits", veri-, > "dev$null") rename(tmp_table3, tabledir//"/"//memnam//"_mod.fits", > "dev$null") del(tmp_table2, veri-, > "dev$null") del(fit_image, veri-, > "dev$null") del(fit_image2, veri-, > "dev$null") del(fit_image3, veri-, > "dev$null") del(fit_image4, veri-, > "dev$null") # fourth method is median: if (verbose) print("median") tabim( tabledir//"/"//memnam//"_mod.fits", \ fit_image, "modes", 1, 256 ) imcopy(fit_image//"[1:128]", fit_image4, verb=verbose) median(fit_image4, fit_image2, smlength,1, bound="reflect", > "dev$null") del(fit_image4, veri-, > "dev$null") imcopy(fit_image//"[129:256]", fit_image4, verb=verbose) median(fit_image4, fit_image3, smlength,1, bound="reflect", > "dev$null") del(fit_image4, veri-, > "dev$null") copy(fit_image, fit_image5, verb=verbose) if (verbose) imhead(fit_image) if (verbose) imhead(fit_image5) imcopy(fit_image2, fit_image5//"[1:128]", verb=verbose) imcopy(fit_image3, fit_image5//"[129:256]", verb=verbose) imtab(fit_image5, tmp_table2, "median", \ format="", pname="col", wcs="logical") tmerge(tabledir//"/"//memnam//"_mod.fits,"//tmp_table2, \ tmp_table3, "merge", allcol+, allrows=256) delete(tabledir//"/"//memnam//"_mod.fits", veri-, > "dev$null") rename(tmp_table3, tabledir//"/"//memnam//"_mod.fits", > "dev$null") del(tmp_table2, veri-, > "dev$null") del(fit_image, veri-, > "dev$null") del(fit_image2, veri-, > "dev$null") del(fit_image3, veri-, > "dev$null") del(fit_image5, veri-, > "dev$null") # let's put in table the column we are going to subtract if (modcol == "subtracted" ) tcalc( tabledir//"/"//memnam//"_mod.fits", "subtracted", \ "imsurfit") else tcalc( tabledir//"/"//memnam//"_mod.fits", "subtracted", \ modcol) # create background if (verbose) print("creating background") tabim(subdir//"/"//memnam//"_mod.fits", bck_fit1, \ modcol, 1, 256) imcopy("nictools$line.fits", bck_fit2, verb=verbose) imcopy( bck_fit1, bck_fit2//"[1:256,1]", verb=verbose) if (verbose) imhead(bck_fit2) blkrep(bck_fit2, bck_fit, 1, 256) del(bck_fit1, veri-, > "dev$null") del(bck_fit2, veri-, > "dev$null") # do subtraction: Start by creating subtraction image: imari(names_tmp3, "-", bck_fit, names_tmp2) # first method is the imari fit if (moddisplay) { if (modcol!="") mosaic(bck_fit//","//names_tmp3// \ ","//names_tmp2) else mosaic(bck_tmp//","//bck_fit//","//names_tmp3// \ ","//names_tmp2) } # display fits if desired: if (moddisplay) { print("data ",tabledir//"/"//memnam//"_mod.fits", > igi_script) print("xc col1", >> igi_script) print("yc modes", >> igi_script) print("lim", >> igi_script) print("points", >> igi_script) print("yc imsurfit", >> igi_script) print("con", >> igi_script) print("box", >> igi_script) print("xlabel column", >> igi_script) print("ylabel mode", >> igi_script) print("color 6", >> igi_script) print("yc gauss", >> igi_script) print("con", >> igi_script) print("color 7", >> igi_script) print("yc median", >> igi_script) print("con", >> igi_script) # finally the one we used: print("color 5", >> igi_script) print("yc subtracted", >> igi_script) print("con", >> igi_script) print("end", >> igi_script) igi( < igi_script ) del(igi_script, veri-) } del(bck_fit, veri-) del(names_tmp3, veri-) } else { # no median to subtract if (verbose) print(" no modes to be subtracted.") names_tmp2= names_tmp3 } # if no backgrounds have been subtracted: if (!access(names_tmp2) && background =="" && \ keyback =="" && !submedian && !subquad) imcopy(member//".fits[sci]", names_tmp2) # resampling: if (verbose) print(" resampling ", resolution) if (inssize > 0) insertim(names_tmp2, insimage, insregion) blkrep(names_tmp2, names_tmp, resolution, resolution) if (verbose) imhead(names_tmp, long-) del(names_tmp2, veri-) # resample errors: if (doerrors) { if (verbose) print(" resampling errors ", resolution) if (access(names_tmp//"_er1.fits")) del(names_tmp//"_er1.fits", veri-) imcopy(member//".fits[err]", names_tmp2, verb=verbose) if (verbose) print("masking ", names_tmp2, " with ", mask) fixpix(names_tmp2, mask, verb-) flpr if (verbose) imhead(names_tmp2) if (inssize > 0) insertim(names_tmp2, insimage, insregion) blkrep(names_tmp2, names_tmp//"_er1", resolution, resolution) imdel(names_tmp2, veri-) flpr #imcalc(names_tmp//"_er1.fits", names_tmp//"_err.fits", "im1**2", verb-) imreplace(names_tmp//"_er1.fits", 9999., low=9999., up=INDEF) imari(names_tmp//"_er1.fits", "*", names_tmp//"_er1.fits", \ names_tmp//"_err.fits") if (access(names_tmp//"_er1.fits")) del(names_tmp//"_er1.fits", veri-) } if (inssize > 0) blkrep(insmask, masks_tmp, resolution, resolution) else blkrep(mask, masks_tmp, resolution, resolution) tabpar(outname//"_asn.fits", "XOFFSET", l ) xoff= int(real(tabpar.value)*resolution+0.5+max_size)-max_size tabpar(outname//"_asn.fits", "YOFFSET", l ) yoff= int(real(tabpar.value)*resolution+0.5+max_size)-max_size print("shifting by ", xoff, yoff) print(member, " ", (real(xoff)/resolution), (real(yoff)/resolution), fbest, \ >> logfile) dpar("asncombine") | head("STDIN", nline=15 , >> logfile) print(" ", >> logfile) imshift(names_tmp, names//l, xoff, yoff, bound="constant", \ interp_type="nearest", constant=0) if (doerrors) imshift(names_tmp//"_err", names//l//"_err", xoff, yoff, bound="constant", \ interp_type="nearest", constant=0) if (verbose) print("computing weights in ", names//l, weightregion) minmax(names//l//weightregion) if (minmax.minval < minmax.maxval) imstat(names//l//weightregion, form-, field="stddev") | \ scan(sigma1) else { print("WARNING: cannot compute weight in ", \ names//l//weightregion, " constant=", minmax.minval) sigma1=1 } hedit(names//l, "weight", (1./sigma1**2), show-, add+, veri-, upd+) imshift(masks_tmp, masks//l//".pl", xoff, yoff, bound="constant",\ interp_type="nearest", constant=1) print( masks//l//".pl", >> masklist) if (subdir != "") { if (!access(subdir)) mkdir(subdir) if (access(subdir//"/"//member//".fits")) \ del(subdir//"/"//member//".fits", veri-) if (access(subdir//"/"//member//"_err.fits")) \ del(subdir//"/"//member//"_err.fits", veri-) if (access(subdir//"/"//member//".pl")) \ del(subdir//"/"//member//".pl", veri-) blkavg(names//l, subdir//"/"//member//".fits", rebin, rebin) if (doerrors) \ blkavg(names//l//"_err", subdir//"/"//member//"_err.fits", rebin, rebin) blkavg(masks//l, subdir//"/"//member//".pl", rebin, rebin) hedit(subdir//"/"//member//".fits", "bpm", \ member//".pl", show-, add+, veri-, upd+) } del(names_tmp, veri-) del(masks_tmp, veri-) if (doerrors) del(names_tmp//"_err.fits", veri-) hedit(names//l, "bpm", masks//l//".pl", show-, add+, veri-, upd+) if (doerrors) hedit(names//l//"_err.fits", "bpm", masks//l//".pl", show-, add+, veri-, upd+) if (l == 1) { print(names//l//".fits", > comb) if (doerrors) print(names//l//"_err.fits", > comb_err) } else { print(names//l//".fits", >> comb) if (doerrors) print(names//l//"_err.fits", >> comb_err) } if (nmembers > 1 && split && (l <= tinfo.nrows/2.) ) print(names//l, >> comb1) if (nmembers > 1 && split && l > tinfo.nrows/2.) print(names//l, >> comb2) if (stopplot) { continue= yes print("\n\nimage display shows the ", l, "th ", \ "image in association ",member," before and after ", \ "subtraction of the mode fits.\n",\ "upper left: subtracted background \n", \ "upper right: image before mode subtraction\n",\ "lower left: image after mode subtraction\n",\ "plot shows mode as function of lines with ", \ "various fits, cyan is the one which will be used.") if (!continue) bye } } imcombine("@"//comb, names_tmp, masktyp="badvalue", maskval=1) flpr blkavg(names_tmp, outname//".fits", rebin, rebin) del(names_tmp, veri-) # rewrite wcs if desired: if (wcsref > 0 && wcsref < tinfo.nrows-1) { if (verbose) print("...updating wcs") tabpar(outname//"_asn.fits", "xoffset", wcsref ) xoff= ( int(real(tabpar.value)*resolution+0.5+max_size)- \ max_size ) tabpar(outname//"_asn.fits", "yoffset", wcsref ) yoff= ( int(real(tabpar.value)*resolution+0.5+max_size)- \ max_size ) tabpar(outname//"_asn.fits", "MEMNAME", wcsref ) member= tabpar.value//"_cal.fits[sci]" imgets(member, "CD1_1") newval= real(imgets.value)/resolution*rebin hedit(outname, "CD1_1", newval, veri-, upd+, show=verbose) imgets(member, "CD1_2") newval= real(imgets.value)/resolution*rebin hedit(outname, "CD1_2", newval, veri-, upd+, show=verbose) imgets(member, "CD2_1") newval= real(imgets.value)/resolution*rebin hedit(outname, "CD2_1", newval, veri-, upd+, show=verbose) imgets(member, "CD2_2") newval= real(imgets.value)/resolution*rebin hedit(outname, "CD2_2", newval, veri-, upd+, show=verbose) imgets(member, "CRPIX1") if (verbose) print("xoffset:", xoff, " crpix1:", imgets.value) newval= (real(imgets.value)*resolution+xoff+ \ 1./resolution-1)/rebin - 1/rebin+1 hedit(outname, "CRPIX1", newval, veri-, upd+, show=verbose) if (verbose) print("yoffset:", yoff, " crpix2:", imgets.value) imgets(member, "CRPIX2") newval= (real(imgets.value)*resolution+yoff+ \ 1./resolution-1)/rebin - 1/rebin+1 hedit(outname, "CRPIX2", newval, veri-, upd+, show=verbose) } if (doerrors) { if (access(outname//"_err.fits")) del(outname//"_err.fits", veri-) if (verbose) { print("adding files in", comb_err, ", \ output=", names_tmp//"_err") type(comb_err) } imcombine("@"//comb_err, names_tmperr, combine="average", \ weight="none", scale="none", reject="none") imcombine("@"//masklist, names_mask, combine="average", \ weight="none", scale="none", reject="none") del(masklist, veri-) imari(names_mask, "*", nmembers, names_mask2) imdel(names_mask, veri-) imari(nmembers, "-", names_mask2, names_mask) imdel(names_mask2, veri-) if (verbose) imstat(names_mask) imari(names_tmperr, "/", names_mask, names_tmperr2, divzero=100) imdel(names_tmperr, veri-) imdel(names_mask, veri-) imfunc(names_tmperr2, names_tmperr, "sqrt", verbose=verbose) blkavg(names_tmperr, outname//"_err.fits", rebin, rebin) del(names_tmperr, veri-) del(names_tmperr2, veri-) del("@"//comb_err, veri-) del(comb_err, veri-) # make mos image if (access(outname//"_mos.fits")) del(outname//"_mos.fits", veri-) imcopy(outname//".fits", outname//"_mos.fits[sci,inherit]", \ verb=verbose) imcopy(outname//"_err.fits", \ outname//"_mos.fits[err,inherit,append]", verb=verbose) #imcopy(names_tmp, outname//"_mos.fits[dq,append]") } if (split && nmembers > 1 ) { if (access(subdir//"/"//outname//"_asn.fits")) \ del (subdir//"/"//outname//"_asn.fits", veri-) copy(outname//"_asn.fits", subdir//"/"//outname//"_asn.fits") if (access(outname//"_1.fits")) del(outname//"_1.fits", veri-) if (access(outname//"_2.fits")) del(outname//"_2.fits", veri-) if (access(outname//"_diff.fits")) del(outname//"_diff.fits", veri-) imcombine("@"//comb1, names_tmp, masktyp="badvalue", maskval=1) flpr blkavg(names_tmp, outname//"_1.fits", rebin, rebin) if (access(names_tmp//".fits")) del(names_tmp//".fits", veri-) else del(names_tmp, veri-) imcombine("@"//comb2, names_tmp, masktyp="badvalue", maskval=1) flpr blkavg(names_tmp, outname//"_2.fits", rebin, rebin) imari(outname//"_1.fits", "-", outname//"_2.fits", \ outname//"_diff.fits") del(comb1, veri-) del(comb2, veri-) del(names_tmp, veri-) } del("@"//comb, veri-) del(comb, veri-) if (verbose) print("deleting :", masks//"*"//".pl") del(masks//"*"//".pl", veri-) hedit(outname//".fits","bpm", " ", delete+, veri-, upd+, show-) if (split && nmembers > 1) { hedit(outname//"_1.fits","bpm", " ", delete+, veri-, upd+, show-) hedit(outname//"_2.fits","bpm", " ", delete+, veri-,upd+, show-) hedit(outname//"_diff.fits","bpm", " ", delete+, veri-, \ upd+, show-) if (resdisp) mosaic(outname//".fits,"//outname//"_diff.fits,"// \ outname//"_1.fits,"//outname//"_2.fits") print("\n\nimage display shows the comb image in ", \ "association ", outname, ":\n", \ "upper left: combinened image\n",\ "upper right: difference between lower two\n",\ "lower left: first half of images combinened\n",\ "lower right: second half of images combinened\n") } else { if (resdisp) display(outname//".fits", 1, fill+) } end