###################################### # nicreduce.cl ###################################### # Script to perform the standard NICMOS pipeline reduction ### # To set up for use: # task nicreduce=(path)/nicreduce.cl # set niccal=(path to calibration files)/ ### # to use: use as an Iraf task "nicreduce" # nicreduce("n46*.fits","um673H") ### # 1997.10.10 Based on B.McLeod's mmtred.cl (30 May 96 vers) # 1997.10.14 First working version, J.Lehar and B.McLeod # SCCS @(#)nicreduce.cl 1.15 02/18/99 10:30:39 ###################################### procedure nicreduce(inputimages,finalimage) ### required arguments string inputimages {"", prompt="Raw input images"} string finalimage {"", prompt="Output combined image"} ### additional files needed string darkfile {"niccal$dark.fits", prompt="dark frame (*.fits)"} string flatfile {"niccal$flath.fits", prompt="flat-field frame (*.fits)"} string maskfile {"niccal$badmask_2.fits", prompt="bad-pixel mask file (*.fits)"} string *list2 {"", prompt="list2, ignore"} ### option flags bool dodark {yes, prompt="subtract dark frame?"} bool doff1 {yes, prompt="find pixel photon-rates? (fullfit)"} bool dobfix {yes, prompt="pedestal: correct bias error estimate?"} bool donbias {yes, prompt="pedestal: correct bias levels?"} bool doff2 {yes, prompt="pedestal: adjust pixel photon-rates? (fullfit)"} bool doflat {yes, prompt="do flat-field correction?"} bool domedian {no, prompt="median filter rows&cols to remove striping?"} bool dolinfit {no, prompt="remove linear fits from each half row&col?"} bool doskydark {no, prompt="remove 'sky dark'=median of subexp?"} bool domagnif {yes, prompt="magnify images times magfactor?"} bool setmask {yes, prompt="add spots to mask file?"} bool domask {yes, prompt="construct final mask?"} bool setshift {yes, prompt="set image shifts with imexam?"} bool doshift {yes, prompt="make shifted data?"} bool setxreg {yes, prompt="set cross-corr region and weight noise box?"} bool doxreg {yes, prompt="register shifted data by cross-correlating?"} bool setcomb {yes, prompt="review and select images for combine?"} bool docomb {yes, prompt="final combine?"} ### other parameters string *list3 {"", prompt="list3, ignore"} bool dotidy {no, prompt="delete intermediate image files?"} bool dofast {no, prompt="skip pedestal effect corrections?"} bool doweight {yes, prompt="weight by image rms estimates?"} real magfactor {2., prompt="magnification factor"} string interp {"spline3", prompt="interpolation type"} real ff1sigma {3., prompt="sigma for fullfit1"} real ff2sigma {5., prompt="sigma for fullfit2"} real readnoise {10., prompt="readnoise for fullfit"} int nskyrej {1., prompt="#high pix to reject for skydark?"} real z1 {0., prompt="lowest pixel value to display"} real z2 {2., prompt="highest pixel value to display"} real spotsize {4., prompt="default spot radius in pixels"} string tempname {"temp", prompt="temp filename prefix"} ### begin ### define some variables int ipos1,ipos2 int i1, i2, j1, j2, i3, j3, i4, j4, veryear, nimage, ncomb real spotsz, u,v,w,magfactr,zz1,zz2 string images, fimage, interpol, s4,s5,s6,s7,s8,s9 images=inputimages if(docomb || domedian) fimage=finalimage magfactr=magfactor spotsz=spotsize interpol=interp ### load packages as needed images if(setshift || setmask) tv tables if(domedian) imfilter nicred set imtype=fits print("****** Making list files") path if(access("raw.list")) delete("raw.list",ver-) if(access(tempname//".list")) delete(tempname//".list",ver-) files(images,>>"raw.list") type("raw.list") i=0 list="raw.list" while(fscan(list,s1) != EOF) { i=i+1 s2=tempname//i print(s2, >> tempname//".list") } nimage=i print("List of temporary file names: ",tempname//".list") type(tempname//".list") flpr ### determine iraf version i=fscan(cl.version,s1,s1,s1,veryear) ### if(dodark) { ####################### # subtract dark frames ####################### if(access("dark.list")) delete("dark.list",ver-) files(darkfile,>>"dark.list") list="raw.list" list2="dark.list" list3=tempname//".list" while(fscan(list,s1) != EOF) { j=fscan(list2,s2) if(j==EOF) { list2="dark.list";j=fscan(list2,s2) } j=fscan(list3,s3) print("****** Darksub ",s2," from ",s1," --> ",s3) path if(access(s3//"_sub.fits")) delete(s3//"_sub.fits",verify-) darksub(s1,s2,s3//"_sub") imgets(s3//"_sub[0]","TARGNAME") hedit(s3//"_sub[0]","title",imgets.value, add+,del-,verify-,show+,update+) hedit(s3//"_sub[0]","RAWFILE",s1, add+,del-,verify-,show+,update+) hedit(s3//"_sub[0]","DARKSUB",s2, add+,del-,verify-,show+,update+) } flpr } ### if(doff1) { ####################### # first fullfit: # initial guess for count rates ####################### list=tempname//".list" while(fscan(list,s1) != EOF) { print("****** fullfit1 for ",s1) path if(access(s1//"_sub_s.fits")) delete(s1//"_sub_s.fits",verify-) if(access(s1//"_sub_i.fits")) delete(s1//"_sub_i.fits",verify-) if(access(s1//"_sub_n.fits")) delete(s1//"_sub_n.fits",verify-) fullfitbam(s1//"_sub.fits",ff1sigma, readnoise) hedit(s1//"_sub_s","FF1SIG",ff1sigma, add+,del-,verify-,show+,update+) hedit(s1//"_sub_s","FF1RN",readnoise, add+,del-,verify-,show+,update+) if(dotidy) { #delete(s1//"_sub.fits",verify-) delete(s1//"_sub_i.fits",verify-) delete(s1//"_sub_n.fits",verify-) delete(s1//"_sub.fits.log",verify-) } } flpr } ### if(dobfix && !dofast) { ####################### # biasfix: correct initial count rate guess # for the bias error ####################### list=tempname//".list" while(fscan(list,s1) != EOF) { print("****** biasfix for ",s1) path if(access(s1//"_sub_scorr.fits")) delete(s1//"_sub_scorr.fits",verify-) #biasfix(s1//"_sub_s.fits",osfn(flatfile), # s1//"_sub_scorr.fits",0.,0.) biasfix(s1//"_sub_s.fits",osfn(flatfile)//"[SCI]", s1//"_sub_scorr.fits",0.,0.) hedit(s1//"_sub_scorr","BFIXFLAT",flatfile, add+,del-,verify-,show+,update+) if(dotidy) delete(s1//"_sub_s.fits",verify-) } flpr } ### if(donbias && !dofast) { ####################### # nicbias: use corrected count rates # to correct the bias levels ####################### list=tempname//".list" while(fscan(list,s1) != EOF) { print("****** nicbias for ",s1) path if(access(s1//"_bias.fits")) delete(s1//"_bias.fits",verify-) nicbias(s1//"_sub",s1//"_sub_scorr",s1//"_bias",save-) hedit(s1//"_bias[0]","NICBIAS",s1//"_sub_scorr", add+,del-,verify-,show+,update+) if(dotidy) { delete(s1//"_sub.fits",verify-) delete(s1//"_sub_scorr.fits",verify-) } } flpr } ### if(doff2 && !dofast) { ####################### # second fullfit: # determine better count rates ####################### list=tempname//".list" while(fscan(list,s1) != EOF) { print("****** fullfit2 for ",s1) path if(access(s1//"_bias_s.fits")) delete(s1//"_bias_s.fits",verify-) fullfitbam(s1//"_bias.fits",ff2sigma,readnoise) hedit(s1//"_bias_s","FF2SIG",ff2sigma, add+,del-,verify-,show+,update+) hedit(s1//"_bias_s","FF2RN",readnoise, add+,del-,verify-,show+,update+) if(dotidy) { delete(s1//"_bias.fits",verify-) delete(s1//"_bias_i.fits",verify-) delete(s1//"_bias_n.fits",verify-) delete(s1//"_bias.fits.log",verify-) } } flpr } ### if(doflat) { ####################### # nicflatten: correct for flat field distortion. ####################### list=tempname//".list" while(fscan(list,s1) != EOF) { if(access(s1//"_flt.fits")) delete(s1//"_flt.fits",verify-) if(!dofast) { print("****** nicflatten for ",s1) path #nicflatten(s1//"_bias_s",osfn(flatfile), # s1//"_flt",0.,0.) nicflatten(s1//"_bias_s",osfn(flatfile)//"[SCI]", s1//"_flt",0.,0.) hedit(s1//"_flt","NFLTFLAT",flatfile, add+,del-,verify-,show+,update+) if(dotidy) delete(s1//"_bias_s.fits",verify-) } if(dofast) { print("****** flatten for ",s1) #imarith(s1//"_sub_s","*",osfn(flatfile),s1//"_flt") imarith(s1//"_sub_s","*",osfn(flatfile)//"[SCI]",s1//"_flt") hedit(s1//"_flt","FLATFILE",flatfile, add+,del-,verify-,show+,update+) if(dotidy) delete(s1//"_sub_s.fits",verify-) } flpr } } ### if(domedian) { ####################### # median filter each image to remove striping ####################### print("****** median filtering each image") path list=tempname//".list" if(access(tempname//"_med.fits")) delete(tempname//"_med.fits",ver-) while(fscan(list,s1) != EOF) { ### set the clipping range if(access("junk_med.dat")) delete("junk_med.dat",ver-) histpeak(s1//"_flt.fits",>>"junk_med.dat") list2="junk_med.dat" j=fscan(list2,x,y);z=y*1.5 ### median filter the columns imgets(s1//"_flt.fits", "i_naxis2") j=int(imgets.value) median(s1//"_flt.fits",tempname//"_med.fits", zloreject=x-z,zhireject=x+z, xwindow=1,ywindow=j,boundary="wrap",verb+ ) imarith(s1//"_flt.fits","-",tempname//"_med.fits",s1//"_flt.fits") delete(tempname//"_med.fits",ver-) flpr ### median filter the rows imgets(s1//"_flt.fits", "i_naxis1") j=int(imgets.value) median(s1//"_flt.fits",tempname//"_med.fits", zloreject=x-z,zhireject=x+z, xwindow=j,ywindow=1,boundary="wrap",verb+ ) imarith(s1//"_flt.fits","-",tempname//"_med.fits",s1//"_flt.fits") delete(tempname//"_med.fits",ver-) flpr } } ### if(dolinfit) { ####################### # remove linear fit to each half row/column ####################### print("****** remove linear fit to each half row/col") path list=tempname//".list" if(access(tempname//"_med.fits")) delete(tempname//"_med.fits",ver-) while(fscan(list,s1) != EOF) { ### set the clipping range if(access("junk_med.dat")) delete("junk_med.dat",ver-) histpeak(s1//"_flt.fits",>>"junk_med.dat") list2="junk_med.dat" j=fscan(list2,x,y);z=y*1.5 imgets(s1//"_flt.fits", "i_naxis2") j=int(imgets.value) ### clip out source flux if(access("junk.fits")) delete("junk.fits",ver-) copy(s1//"_flt.fits","junk.fits") imreplace("junk.fits",0.,lower=INDEF,upper=x-z) imreplace("junk.fits",0.,lower=x+z,upper=INDEF) ### clean up the low columns i1=1;i2=j;j1=1;j2=j/2 s4="["//i1//":"//i2//","//j1//":"//j2//"]" fit1d("junk.fits"//s4,tempname//"_med.fits", axis=2,type="fit",low=2,high=2, interac-,function="legendre",order=2 ) imarith(s1//"_flt.fits","-",tempname//"_med.fits", s1//"_flt.fits") imarith("junk.fits","-",tempname//"_med.fits", "junk.fits") delete(tempname//"_med.fits",ver-) flpr ### clean up the high columns i1=1;i2=j;j1=(j/2)+1;j2=j s4="["//i1//":"//i2//","//j1//":"//j2//"]" fit1d("junk.fits"//s4,tempname//"_med.fits", axis=2,type="fit",low=2,high=2, interac-,function="legendre",order=2 ) imarith(s1//"_flt.fits","-",tempname//"_med.fits", s1//"_flt.fits") imarith("junk.fits","-",tempname//"_med.fits", "junk.fits") delete(tempname//"_med.fits",ver-) flpr ### clean up the low rows i1=1;i2=j/2;j1=1;j2=j s4="["//i1//":"//i2//","//j1//":"//j2//"]" fit1d("junk.fits"//s4,tempname//"_med.fits", axis=1,type="fit",low=2,high=2, interac-,function="legendre",order=2 ) imarith(s1//"_flt.fits","-",tempname//"_med.fits", s1//"_flt.fits") imarith("junk.fits","-",tempname//"_med.fits", "junk.fits") delete(tempname//"_med.fits",ver-) flpr ### clean up the high rows i1=(j/2)+1;i2=j;j1=1;j2=j s4="["//i1//":"//i2//","//j1//":"//j2//"]" fit1d("junk.fits"//s4,tempname//"_med.fits", axis=1,type="fit",low=2,high=2, interac-,function="legendre",order=2 ) imarith(s1//"_flt.fits","-",tempname//"_med.fits", s1//"_flt.fits") imarith("junk.fits","-",tempname//"_med.fits", "junk.fits") delete(tempname//"_med.fits",ver-) delete("junk.fits",ver-) flpr } } ### if(doskydark) { ####################### # remove 'sky dark' ####################### print("****** removing 'sky dark', median of subexp") path ### make the sky dark if(access("junk.fits")) delete("junk.fits",ver-) imcombine(tempname//"?_flt.fits","junk.fits", plfile="",sigma="",logfile="STDOUT",combine="median", reject="minmax",nlow=0,nhigh=nskyrej,nkeep=1, masktype="none",scale="none",zero="none",weight="none", ) ### subtract the sky dark list=tempname//".list" while(fscan(list,s1) != EOF) { imarith(s1//"_flt.fits","-","junk.fits",s1//"_flt.fits") } ### clean up delete("junk.fits",ver-) } ### if(domagnif) { ####################### # magnify images by factor of 2 ####################### print("****** magnifying flat-fielded images") path list=tempname//".list" while(fscan(list,s1) != EOF) { ### first interpolate masked pixels print(" fixup flt for ",s1) if(access(s1//"_fix.fits")) delete(s1//"_fix.fits",verify-) nicfixpix(s1//"_flt", s1//"_fix", osfn(maskfile), 0, 0) ### then do the magnification print(" magnify flt for ",s1," by ",magfactr) if(access(s1//"_mag.fits")) delete(s1//"_mag.fits",verify-) magnify(s1//"_fix.fits",s1//"_mag.fits",magfactr,magfactr, fluxcon+,interp=interpol) hedit(s1//"_mag.fits","MAGNIFY",2., add+,del-,verify-,show+,update+) hedit(s1//"_mag.fits","BPM",s1//"_mag.pl", add+,del-,verify-,show+,update+) if(dotidy) { delete(s1//"_flt.fits",verify-) delete(s1//"_fix.fits",verify-) } ### on to next image } flpr } ### if(setmask) { ####################### # mark additional spots on mask ####################### ### get default mask if(access("junk.fits")) delete("junk.fits",ver-) chpixtype(maskfile,"junk.fits","real") magnify("junk.fits","junk.fits",magfactr,magfactr, fluxcon-,interp="linear") imreplace("junk.fits",0.,lower=INDEF,upper=0.1) imreplace("junk.fits",1.,lower=0.05,upper=INDEF) if(access("mask_mag.pl")) delete("mask_mag.pl",ver-) imcopy("junk.fits","mask_mag.pl") delete("junk.fits",ver-) print("****** marking additional mask spots") path s2="n";s5="overlay";zz1=z1;zz2=z2 list=tempname//".list" while(s2 != "q") { ### select spot location if(s2 == "n") { j=fscan(list,s1) s6=tempname//"_";i1=stridx("_",s6) s6=s1//"_";i2=stridx("_",s6) s7="spots"//substr(s1,i1,i2-1) if(j==EOF) { list=tempname//".list" j=fscan(list,s1) s6=tempname//"_";i1=stridx("_",s6) s6=s1//"_";i2=stridx("_",s6) s7="spots"//substr(s1,i1,i2-1) } } if(s2=="n" || s2=="u") { s5="overlay" ### update mask and display image if(access(s1//"_mag.pl")) delete(s1//"_mag.pl",ver-) imcopy("mask_mag.pl",s1//"_mag.pl",ver-) if(access("spots.list")) { print("Adding global spots = spots.list:") if(access(s1//"_junk.pl")) delete(s1//"_junk.pl",ver-) maskspots(s1//"_mag.pl","spots.list",s1//"_junk.pl", template="nicred$Spots/spot",verb-) delete(s1//"_mag.pl",ver-) rename(s1//"_junk.pl",s1//"_mag.pl") } if(access(s7//".list")) { print("Adding local spots = ",s7//".list:") if(access(s1//"_junk.pl")) delete(s1//"_junk.pl",ver-) maskspots(s1//"_mag.pl",s7//".list",s1//"_junk.pl", template="nicred$Spots/spot",verb-) delete(s1//"_mag.pl",ver-) rename(s1//"_junk.pl",s1//"_mag.pl") } } if(s2=="n" || s2=="m" || s2=="u" || s2=="<" || s2==">" || s2=="=") { if(s2 == "m") { s4=s5 if(s4=="none") s5="overlay" if(s4=="overlay") s5="none" } if(veryear < 1997) { display(s1//"_mag.fits",1, erase+,border-,select+,repeat-,fill-, zs-,zr-,z1=zz1,z2=zz2,ztran="log") } else { display(s1//"_mag.fits",1, bpmask="BPM",bpdispl=s5,bpcolor="red", erase+,border-,select+,repeat-,fill-, zs-,zr-,z1=zz1,z2=zz2,ztran="log") } } if(s2=="n" || s2=="?") { print("****** mark spot centers for masks") print(" image=",s1//"_mag.fits, mask=",s1//"_mag.pl") print(" --- use 'n' to see next image") print(" use 'u' to update image with new spots") print(" use 'q' to quit") print(" --- use 'g' to mark a global spot of radius=",spotsz) print(" use 's' to mark a local spot of radius=",spotsz) print(" use '+' to increase the spot size") print(" use '-' to decrease the spot size") print(" use 'd' to delete the indicated spot") print(" --- use '?' to show this help info") print(" use 'l' to show spot lists") print(" use 'z' to delete global spot list") print(" use 'x' to delete local spot list") print(" --- use 'm' to toggle mask display") print(" use '>' to double the display z scale of ",zz1,zz2) print(" use '<' to halve the display z scale of ",zz1,zz2) print(" use '=' to reset the display z scale to ",z1,z2) print(" (local=affects only this image, global=all images)" } j=fscan(imcur,x,y,k,s2) if(s2 == "+") { spotsz=spotsz+1; if(spotsz>20) spotsz=20 } if(s2 == "-") { spotsz=spotsz-1; if(spotsz<1) spotsz=1 } if(s2 == "<") { zz1=zz1/2;zz2=zz2/2;print(" display scale=",zz1,zz2) } if(s2 == ">") { zz1=zz1*2;zz2=zz2*2;print(" display scale=",zz1,zz2) } if(s2 == "=") { zz1=z1;zz2=z2;print(" display scale=",zz1,zz2) } if(s2 == "z") { delete("spots.list",ver-);print("deleted spots.list") } if(s2 == "x") { delete(s7//".list",ver-);print("deleted ",s7//".list") } if(s2 == "g") { print(" global spot x,y,size=",x,y,spotsz) print(x,y,spotsz, >>"spots.list") } if(s2 == "s") { print(" local spot x,y,size=",x,y,spotsz) print(x,y,spotsz, >>s7//".list") } if(s2 == "l") { print("Global spots.list:");type("spots.list") print("Local ",s7//".list:");type(s7//".list") } if(s2 == "d" && access("spots.list")) { ### delete from global list if(access("junk_spots.list")) delete("junk_spots.list",ver-) list2="spots.list" while(fscan(list2,u,v,w)!=EOF) { z=sqrt( (x-u)**2 + (y-v)**2 ) if(z>w) { print(u,v,w,>>"junk_spots.list") } else { print("Deleted global spot=",u,v,w) } } if(access("spots.list")) { delete("spots.list",ver-) rename("junk_spots.list","spots.list") } } if(s2 == "d" && access(s7//".list")) { ### delete from local list if(access("junk_spots.list")) delete("junk_spots.list",ver-) list2=s7//".list" while(fscan(list2,u,v,w)!=EOF) { z=sqrt( (x-u)**2 + (y-v)**2 ) if(z>w) { print(u,v,w,>>"junk_spots.list") } else { print("Deleted local spot=",u,v,w) } } if(access(s7//".list")) { delete(s7//".list",ver-) rename("junk_spots.list",s7//".list") } ### } } } flpr ### if(domask) { ####################### # construct the masks ####################### print("****** Constructing the image masks") path ### make default mask if(access("junk.fits")) delete("junk.fits",ver-) chpixtype(maskfile,"junk.fits","real") magnify("junk.fits","junk.fits",magfactr,magfactr, fluxcon-,interp="linear") imreplace("junk.fits",0.,lower=INDEF,upper=0.1) imreplace("junk.fits",1.,lower=0.05,upper=INDEF) if(access("mask_mag.pl")) delete("mask_mag.pl",ver-) imcopy("junk.fits","mask_mag.pl") delete("junk.fits",ver-) ### construct masks list=tempname//".list" while(fscan(list,s1) != EOF) { if(access(s1//"_mag.pl")) delete(s1//"_mag.pl",ver-) imcopy("mask_mag.pl",s1//"_mag.pl",ver-) ### add global spots to mask if(access("spots.list")) { print("global spots = spots.list:") if(access(s1//"_junk.pl")) delete(s1//"_junk.pl",ver-) maskspots(s1//"_mag.pl","spots.list",s1//"_junk.pl", template="nicred$Spots/spot",verb-) delete(s1//"_mag.pl",ver-) rename(s1//"_junk.pl",s1//"_mag.pl") } ### add local spots to mask s6=tempname//"_";i1=stridx("_",s6) s6=s1//"_";i2=stridx("_",s6) s7="spots"//substr(s1,i1,i2-1) if(access(s7//".list")) { print("local spots = ",s7//".list:") if(access(s1//"_junk.pl")) delete(s1//"_junk.pl",ver-) maskspots(s1//"_mag.pl",s7//".list",s1//"_junk.pl", template="nicred$Spots/spot",verb-) delete(s1//"_mag.pl",ver-) rename(s1//"_junk.pl",s1//"_mag.pl") } print(" made mask=",s1//"_mag.pl") } } ### if(setshift) { ####################### # set integer shifts to remove dithers ####################### print("****** Setting integer shifts to remove dithers") path ### get default mask if(!access("mask_mag.pl")) { if(access("junk.fits")) delete("junk.fits",ver-) chpixtype(maskfile,"junk.fits","real") magnify("junk.fits","junk.fits",magfactr,magfactr, fluxcon-,interp="linear") imreplace("junk.fits",0.,lower=INDEF,upper=0.1) imreplace("junk.fits",1.,lower=0.05,upper=INDEF) imcopy("junk.fits","mask_mag.pl") delete("junk.fits",ver-) } if(access("imexam.out")) delete("imexam.out",ver-) list=tempname//".list" while(fscan(list,s1) != EOF) { ### use default mask if none if(!access(s1//"_mag.pl")) imcopy("mask_mag.pl",s1//"_mag.pl",ver-) ### display image and get ref position print("****** mark reference object on ",s1) print(" use 'a' to mark object") print(" then use 'q' to quit") if(veryear<1997) { display(s1//"_mag.fits",1, erase+,border-,select+,repeat-,fill-, zs-,zr-,z1=z1,z2=z2,ztran="log") } else { display(s1//"_mag.fits",1, bpmask="BPM",bpdispl="overlay",bpcolor="red", erase+,border-,select+,repeat-,fill-, zs-,zr-,z1=z1,z2=z2,ztran="log") } unlearn imexamine unlearn rimexam if(veryear>=1997) rimexam.fittype="gaussian" imexam(s1//"_mag.fits",keeplog+,logfile="imexam.out") } ### compute shifts if(access("junk.list")) delete("junk.list",ver-) if(access("shifts.list")) delete("shifts.list",ver-) i=0;i3=9999;j3=9999 list2="imexam.out" while(fscan(list2,s1,s2)!=EOF) { if(s1!="#") { i=i+1 x=real(s1);y=real(s2) if(i==1) {i1=nint(x);j1=nint(y)} i2=nint(i1-x);j2=nint(j1-y) print(i2,j2,>>"junk.list") i3=min(i3,i2);j3=min(j3,j2) } } list2="junk.list" while(fscan(list2,i,j)!=EOF) { i1=i-i3;j1=j-j3 print(i1,j1,>>"shifts.list") } flpr delete("junk.list",ver-) } ### if(doshift) { ####################### # shift the images by the integer shift ####################### print("****** Shifting images with integer shifts") path ### get default mask if(!access("mask_mag.pl")) { if(access("junk.fits")) delete("junk.fits",ver-) chpixtype(maskfile,"junk.fits","real") magnify("junk.fits","junk.fits",magfactr,magfactr, fluxcon-,interp="linear") imreplace("junk.fits",0.,lower=INDEF,upper=0.1) imreplace("junk.fits",1.,lower=0.05,upper=INDEF) imcopy("junk.fits","mask_mag.pl") delete("junk.fits",ver-) } ### update old versions of shift file if(access("junk.list")) delete("junk.list",ver-) i3=9999;j3=9999 list2="shifts.list" while(fscan(list2,x,y)!=EOF) { i2=nint(x);j2=nint(y) print(i2,j2,>>"junk.list") i3=min(i3,i2);j3=min(j3,j2) } list2="junk.list" delete("shifts.list",ver-) ### current shift in i,j while(fscan(list2,i,j)!=EOF) { i1=i-i3;j1=j-j3 print(i1,j1,>>"shifts.list") } flpr delete("junk.list",ver-) ### make list of magnified images if(access("junk_mag.list")) delete("junk_mag.list",ver-) list=tempname//".list" while(fscan(list,s1) != EOF) { print(s1//"_mag.fits",>>"junk_mag.list") if(!access(s1//"_mag.pl")) imcopy("mask_mag.pl",s1//"_mag.pl",ver-) } ### k=0 # k is current image list=tempname//".list" while(fscan(list,s1) != EOF) { k=k+1 print("make shifted image:",s1//"_shift.fits") ### delete old files if(access(s1//"_shift.fits")) delete(s1//"_shift.fits",ver-) if(access(s1//"_rejec.pl")) delete(s1//"_rejec.pl",ver-) if(access(s1//"_shift.comb_log")) delete(s1//"_shift.comb_log",ver-) ### reweight data to single out this exposure if(access("junk_shiftwt.list")) delete("junk_shiftwt.list",ver-) i2=0 while(i2 < nimage) { i2=i2+1 s2=tempname//i2 ### weight other images by exposure imgets(s2//"_mag.fits", "EXPTIME") z=real(imgets.value) if(i2==k) { y=z z=z*1000. } print(z, >> "junk_shiftwt.list") } #type("junk_shiftwt.list") ### now combine imcombine( "@junk_mag.list", s1//"_shift.fits", plfile=s1//"_rejec.pl", sigma="", logfile="STDOUT", combine="average", reject="none", project=no, outtype="real", offsets="shifts.list", masktype="badvalue", maskvalue=1., blank=0., scale="none", zero="none", weight="@junk_shiftwt.list", statsec="", expname="", lthreshold=-100, hthreshold=INDEF, nlow=1, nhigh=1, nkeep=1, mclip=yes, lsigma=3., hsigma=3., rdnoise=7., gain=7., snoise=0., sigscale=0.1, pclip=-0.5, grow=0, >> s1//"_shift.comb_log" ) hedit(s1//"_shift.fits", "EXPTIME", y, add-,del-,verify-,show-,update+) type(s1//"_shift.comb_log") delete("junk_shiftwt.list",ver-) flpr ### construct mask for the shifted image print(" construct ",s1//"_shift.pl") if(access(s1//"_shift.pl")) delete(s1//"_shift.pl",ver-) imcopy(s1//"_shift.fits",s1//"_shift.pl",verb-) imreplace(s1//"_shift.pl",1.,lower=INDEF,upper=INDEF) ### unmask sampled region for final combine i2=0 list2="shifts.list" while(fscan(list2,i,j)!=EOF) { i2=i2+1 if(i2==k) { ### unmask sampled region (1pix smaller for xreg shifts) imgets(s1//"_mag.fits", "i_naxis1");i1=nint(real(imgets.value)) imgets(s1//"_mag.fits", "i_naxis2");j1=nint(real(imgets.value)) s9="["//i+1//":"//i+i1-1//","//j+1//":"//j+j1-1//"]" imreplace(s1//"_shift.pl"//s9,0.,lower=INDEF,upper=INDEF) print(" unmask ",s1//s9) flpr } } ### modify header entries hedit(s1//"_shift.fits","BPM",s1//"_shift.pl", add+,del-,verify-,show-,update+) hedit(s1//"_shift.fits","NCOMBINE",1, add+,del-,verify-,show+,update+) ### apply 1st image shift to WCS reference pixel list2="shifts.list" i2=fscan(list2,i,j) imgets(s1//"_shift.fits","CRPIX1");x=real(imgets.value)+i imgets(s1//"_shift.fits","CRPIX2");y=real(imgets.value)+j hedit(s1//"_shift.fits","CRPIX1",x, add+,del-,verify-,show-,update+) hedit(s1//"_shift.fits","CRPIX2",y, add+,del-,verify-,show-,update+) flpr } delete("junk_mag.list",ver-) } ### if(setxreg) { ####################### # set cross-correlation region and weighting noise box ####################### print("****** Setting cross-corr region and weighting noise box") path if(access("xregdat.list")) delete("xregdat.list",ver-) s3="n";s4="BEST" list=tempname//".list" i1=-1000;j1=0; i2=-1000;j2=0; i3=-1000;j3=0; i4=-10000;j4=0 while( s3 != "q" || i1==-1000 || i2==-1000 ) { ### set refimage and xcor.region and rms.noise.box if(s3 == "n") { j=fscan(list,s1) if(j==EOF) { list=tempname//".list" j=fscan(list,s1) } display(s1//"_shift.fits",1,erase+,border-,select+,repeat-,fill-, zs-,zr-,z1=z1,z2=z2,ztran="log") } if(s3=="n" || s3=="r") { print(" ... image=",s1//"_shift, refimage=",s4//"_shift") print(" use 'n' to see next image") print(" use 'r' to set this image as refimage") print(" use 'x' to reset refimage to BEST") print(" use 'q' to quit, save refimage and regions") print(" ") print(" use 'a' to mark bottom-left corner of xcor") print(" use 'b' to mark top-right corner of xcor") print(" use 'c' to mark blc of weighting noise box") print(" use 'd' to mark trc of weighting noise box") print(" ") print(" (refimage should be one of the better images)") print(" ( refimage=BEST ==> find smallest rms)") print(" (histpeak rms of whole image if no noise box)") } j=fscan(imcur,x,y,k,s3) if(s3 == "r") s4=s1 if(s3 == "x") s4="BEST" if(s3 == "a") { i1=nint(x); j1=nint(y) print(" xcor.blc=",i1,j1) } if(s3 == "b") { i2=nint(x); j2=nint(y) print(" xcor.trc=",i2,j2) } if(s3 == "c") { i3=nint(x); j3=nint(y) print(" noise.blc=",i3,j3) } if(s3 == "d") { i4=nint(x); j4=nint(y) print(" noise.trc=",i4,j4) } } ### refimage print(s4, >>"xregdat.list") ### cross-corr region s2="["//i1//":"//i2//","//j1//":"//j2//"]" print(s2, >>"xregdat.list") print("****** refimage[region]= ",s4//"_shift"//s2) ### weighting noise box if(i3>0&&j3>0) { s2="["//i3//":"//i4//","//j3//":"//j4//"]" print(s2, >>"xregdat.list") print("****** noise.box= ",s2) } } flpr ### if(doxreg) { ####################### # register images with cross-correlation ####################### print("****** Cross correlating for fine registration") path ### get refimage,xcor.region,noise.region list="xregdat.list";s4="" j=fscan(list,s2) # refimage j=fscan(list,s3) # xcor region j=fscan(list,s4) # noise box ### make list of magnified images and choose refimage if(access("junk_shift.list")) delete("junk_shift.list",ver-) if(access("junk_xreg.list")) delete("junk_xreg.list",ver-) list=tempname//".list" z=9999; ### fix to allow for filenames of form temp1_shift.fits, temp1, ### myname_forthis1, myname_forthis_temp1_shift.fits ### for compatibility with nicred1.7 xregdat.list convention ipos1=stridx("_",s2) ipos2=ipos1 while(ipos2>0) { ### find last occurrence of "_" ipos1 = ipos2 + stridx("_",substr(s2,ipos2+1,strlen(s2))) if (ipos1==ipos2) { ### found last, now truncate if necessary if(substr(s2,ipos2+1,strlen(s2))=="shift.fits") s2=substr(s2,1,ipos2-1) ipos2=0 } else { ipos2=ipos1 } } s5=s2 while(fscan(list,s1) != EOF) { ### add to lists print(s1//"_shift.fits",>>"junk_shift.list") print(s1//"_xreg.fits",>>"junk_xreg.list") ### choose refimage if "BEST" if(s2=="BEST") { if(access("imstat.out")) delete("imstat.out",ver-) if(s4!="none") { imstat(s1//"_shift.fits"//s4,fields="stddev",>>"imstat.out") list2="imstat.out" j=fscan(list2,y) j=fscan(list2,y) } if(s4=="none") { histpeak(s1//"_shift.fits",>>"imstat.out") list2="imstat.out" j=fscan(list2,x,y) } ### choose refimage with smallest rms print(" image,rms=",s1//"_shift.fits"//s4,y) if(y>"xregister.log") if (dotidy) { delete("@junk_shift.list",ver-) delete("junk_shift.list",ver-) delete("junk_xreg.list",ver-) } flpr } ### if( setcomb && nimage>1 ) { ####################### # selecting images for final combination ####################### print("****** Selecting images for final combination") ### first make up the default list if(access(fimage//".list")) delete(fimage//".list",ver-) list=tempname//".list" while(fscan(list,s1)!=EOF) { print(s1//"_xreg.fits",>> fimage//".list") } ncomb=nimage ### now inspect the images s3="n" list=fimage//".list" while( s3!="q" ) { if(s3=="n") { j=fscan(list,s1) if(j==EOF) { list=fimage//".list" j=fscan(list,s1) } display(s1,1,erase+,border-,select+,repeat-,fill-, zs-,zr-,z1=z1,z2=z2,ztran="log") } if(s3=="n" || s3=="?") { print("****** select images for final combine") print(" image=",s1) print(" --- '?' to see this help listing") print(" 'i' to inspect with imexamine") print(" --- 'l' to see the current image list") print(" 'd' to delete this image from list") print(" 'x' to re-initialize the image list") print(" --- 'n' to load next image in display 1") print(" (current image is always in display 1)") print(" '2' copy current image to display 2") print(" '3' copy current image to display 3") print(" '4' copy current image to display 4") print(" --- 'q' to quit and keep current list") print(" (options 'd' and 'x' also invoke an 'n')") } j=fscan(imcur,x,y,k,s3) if(s3 == "i") { print("---- inside imexamine") print(" 'a'==>peak params, 'r'==>radial.fit, 'q'==>quit") imexamine(s1,keeplog-) s3="?" } if(s3 == "l") { print("Current image list for combine:") type(fimage//".list") } if(s3 == "d" && ncomb==1) { print("---- ONLY 1 IMAGE IN LIST, can't delete") } if(s3 == "d" && ncomb>1) { if(access("junk.list")) delete("junk.list",ver-) list2=fimage//".list" while(fscan(list2,s4)!=EOF) if(s4!=s1) print(s4,>>"junk.list") print("revised list:") type("junk.list") delete(fimage//".list",ver-) copy("junk.list",fimage//".list",verb-) delete("junk.list",ver-) s3="n" } if(s3 == "x") { delete(fimage//".list",ver-) list2=tempname//".list" while(fscan(list2,s4)!=EOF) print(s4//"_xreg.fits",>> fimage//".list") ncomb=nimage s3="n" } if( s3=="2" || s3=="3" || s3=="4") { display(s1,int(s3),erase+,border-,select+,repeat-,fill-, zs-,zr-,z1=z1,z2=z2,ztran="log") } } print("Final combine list:") type(fimage//".list") } ### if(docomb) { ####################### # finally combine the registered images ####################### print("****** Combining the registered images") type("xregdat.list") path ### determine weighting list="xregdat.list" j=fscan(list,s2) # past refimage j=fscan(list,s2) # past xcor region j=fscan(list,s2) # get noise box if(j==EOF || s2=="") s2="none" ### construct combine list if missing if(!access(fimage//".list")){ list=tempname//".list" while(fscan(list,s1)!=EOF) print(s1//"_xreg.fits",>> fimage//".list") } ### count input images and setup weights if(access(tempname//"_combwt.list")) delete(tempname//"_combwt.list",ver-) u=0 v=0 ncomb=0 list=fimage//".list" while(fscan(list,s1) != EOF) { ncomb=ncomb+1 ### list of weights if(doweight) { ### weight images by rms estimate if(access("imstat.out")) delete("imstat.out",ver-) if(s2!="none") { imstat(s1//s2,fields="stddev",>>"imstat.out") list2="imstat.out" j=fscan(list2,z) j=fscan(list2,z);w=1.0/(z*z) } if(s2=="none") { histpeak(s1,>>"imstat.out") list2="imstat.out" j=fscan(list2,x,z);w=1.0/(z*z) } v=v+w;u=u+w*z print(w,>>tempname//"_combwt.list") } if(!doweight) { ### weight images by exposure time ncomb=ncomb+1 imgets(s1,"EXPTIME");w=real(imgets.value) print(w,>>tempname//"_combwt.list") } } if(v>0) u=u/v if(!doweight) print("****** weight by exposure") if(doweight) print("****** Weight noise box: ",s2) print("weights:") type(tempname//"_combwt.list") s3="@"//tempname//"_combwt.list" ### delete old files if(access(fimage//".fits")) delete(fimage//".fits",ver-) if(access(fimage//".pl")) delete(fimage//".pl",ver-) if(access(fimage//"_sigma.fits")) delete(fimage//"_sigma.fits",ver-) if(access(fimage//".comb_log")) delete(fimage//".comb_log",ver-) ### now combine print(" mean noise rms=",u) imcombine( "@"//fimage//".list", fimage//".fits", plfile=fimage//".pl", sigma=fimage//"_sigma.fits", logfile="STDOUT", combine="average", reject="avsigclip", project=no, outtype="real", offsets="none", masktype="badvalue", maskvalue=1., blank=0., scale="none", zero="median", weight=s3, statsec="", expname="", lthreshold=-100, hthreshold=INDEF, nlow=1, nhigh=1, nkeep=1, mclip=yes, lsigma=3., hsigma=3., rdnoise=7., gain=7., snoise=0., sigscale=0.1, pclip=-0.5, grow=0, >> fimage//".comb_log" ) ### fix the sigmas file if(v>0) imarith(fimage//"_sigma.fits","max",u,fimage//"_sigma.fits") flpr if(access("comb2.fits")) delete("comb2.fits",ver-) imexpr("a / sqrt ( b - c )","comb2.fits",1.0,ncomb,fimage//".pl") imarith(fimage//"_sigma.fits","*","comb2.fits",fimage//"_sigma.fits") delete("comb2.fits",ver-) flpr ### final stuff type(fimage//".comb_log") if(setshift || setmask || setxreg) { print("****** final combined image=",fimage//".fits") print("Type any key to continue...") path display(fimage//".fits",1,erase+,border-,select+,repeat-,fill-, zs-,zr-,z1=z1,z2=z2,ztran="log") j=fscan(imcur,x,y,k,s2) } if(dotidy) { delete(tempname//"*_xreg.fits",ver-) delete(tempname//"_combwt.list",ver-) } flpr } ### end