procedure tnoise (imgFile,abbName,extn,useRimcur,frm,boxSize,skyFile, fluxThresh,ap1,ap2,ap3,ap4,SNlim,gain,outTh,outRMS,inclMeans) string imgFile {"-",prompt="List of images (full or abbreviated names)"} bool abbName {no, prompt=" Are image names abbreviated?"} string extn {prompt=" Extension appended to abbreviated image names"} bool useRimcur {yes, prompt="Use rimcursor to define sky regions?"} int frm {1, prompt=" Ximtool frame for display image"} int boxSize {70, prompt=" Side-length of region for statistics?"} string skyFile {"--", prompt=" File of sky regions if no use of rimcursor."} bool fluxThresh {yes, prompt="Calculate flux thresholds as well?"} real ap1 {4, prompt=" Aperture 1 (diameter)"} real ap2 {8, prompt=" Aperture 2 (diameter)"} real ap3 {12, prompt=" Aperture 3 (diameter)"} real ap4 {20, prompt=" Aperture 4 (diameter)"} real SNlim {1, prompt=" Photometry limit (as a S/N ratio)"} real gain {1.11, prompt=" CCD gain (e-/ADU)"} string outTh {"--", prompt=" Output file for flux threshold data"} string outRMS {"---", prompt="Output file for sky data"} bool inclMeans {yes, prompt="Include mean sky values with std devs?"} struct *flist1 struct *flist2 # list-directed structures struct *flist3 begin # Variable Declarations: string rawName # raw (abbreviated) image name string fullName # full image name int j,k,m # image counters string imgjnk # junk string to store first image name int x1, x2 int y1, y2 # dimensions of imstat region int xc, yc # central coordinates of current region string reg # string containing image stat region int Nreg # counter of the number of regions real pie # pi real area # current aperture pixel area real area1, area2, area3, area4 # input aperture areas real numer real denom # numerator and denominator in threshold expression real sig # sigma of current background real ap # current aperture diameter real Flim # current limiting flux string outStr # string of output flux threshold values string flx # current flux threshold (formatted) string padding # column spacer real rms # current rms value string c1 # first character of all filenames (for checking, since # FORTRAN will not read filenames beginning with a no.) real medSky # median sky value from all sky regions on an image bool NregOdd # = yes if Nreg is odd, = no if even real Nreal # real version of Nreg int p # median counter string jnkstr # junk string int centreIndex real centreVal # values associated with finding median pie=3.1415927 if (! defpac ("images")) { print("Package IMAGES needs to be loaded.") bye } if (! defpac ("lists")) { print("Package LISTS needs to be loaded.") bye } if (! defpac ("proto")) { print("Package PROTO needs to be loaded.") bye } # Cleaning up any files to be used: delete (files="temp_tnse*",go_ahea+,ver-) delete (files="temp_junk",go_ahea+,ver-) delete (files=outRMS,go_ahea+,ver-) if (fluxThresh==yes) { delete (files=outTh,go_ahea+,ver-) } # Checking for numeric starters in supplied file names: # (FORTRAN will not read such names if given as input!!) c1 = substr(outTh, 1, 1) if ((c1=="0")||(c1=="1")||(c1=="2")||(c1=="3")|| (c1=="4")||(c1=="5")||(c1=="6")||(c1=="7")||(c1=="8")||(c1=="9")) { print("** ERROR: first character of 'outTh' ("//outTh//") is numeric **") print(" ") print(" --> Please change incorrect file name and re-run") print(" ") print (" Done.") print (" ") print (" ") bye } # Obtaining regions: ----------------------------------------------- if (useRimcur==yes) { flist1 = imgFile imgjnk = fscan(flist1, rawName) if (abbName==yes) { fullName=rawName//extn } if (abbName==no) { fullName=rawName } display (image=fullName,frame=frm,erase+,bord-, sele+,repe-,fill-,zscal+,ztrans="linear", >> "temp_junk") print (" ") print (" *** "//fullName//" displayed in Ximtool frame "//frm//" ***") print (" ") print (" --> Press any key at each region and ctrl+ when done.") rimcursor(image=fullName,wcs="logical",wxformat="",wyformat="", cursor="",mode="ql", >> "temp_tnse1") print (" ") print (" Please wait ...") !awk '{print int($1)" "int($2)}' temp_tnse1 >temp_tnse2 flist2 = "temp_tnse2" Nreg=0 while ( fscan(flist2, xc, yc) != EOF) { Nreg=Nreg+1 x1=xc-int(boxSize/2) x2=xc+int(boxSize/2)-1 y1=yc-int(boxSize/2) y2=yc+int(boxSize/2)-1 reg = "["//x1//":"//x2//","//y1//":"//y2//"]" print (reg, >> "temp_tnse3") } } if (useRimcur==no) { print (" ") if (skyFile=="-") { print (" ** NEED TO SPECIFY AN INPUT FILE **") print (" ") print (" Please try again.") print (" ") bye } print ("Reading in regions from the file "//skyFile) flist2 = skyFile Nreg=0 while ( fscan(flist2, reg) != EOF) { Nreg=Nreg+1 print (reg, >> "temp_tnse3") } } # type(input_fi="temp_tnse2",map_cc+,device="terminal",flpar-) print (" ") type(input_fi="temp_tnse3",map_cc+,device="terminal",flpar-) print (" ") print (" There are "//Nreg//" regions; please wait ...") Nreal = 1.0*Nreg if ( (int(Nreal/2))==(Nreal/2) ) { print (" -> this is even.") NregOdd=no } else { print (" -> this is odd.") NregOdd=yes } print (" ") # Creating main image list (temp_tnse4): ---------------------------- flist1 = imgFile j=0 while ( fscan(flist1, rawName) != EOF) { j=j+1 if (abbName==yes) { fullName=rawName//extn } if (abbName==no) { fullName=rawName } print (fullName//" ", >> "temp_tnse9") # Region loop: ---------------------------- flist2 = "temp_tnse3" while ( fscan(flist2, reg) != EOF) { print (fullName//reg, >> "temp_tnse4") } imstatistics(images="@temp_tnse4",fields="image,npix,mean,stddev", lower=INDEF,upper=INDEF,binwid=0.1,format+,mode="ql", >> "temp_tnse5") type (input_fi="temp_tnse5",map_cc+,dev="terminal",flpar-) # MEAN RMS VALUES: # Obtaining a median: ----------------- delete (files="temp_median*",go_ahea+,ver-) fields(files="temp_tnse5",fields="4",lines="1-",quit_if-,print_f-, >> "temp_median1") sort(input_fi="temp_median1",column=1,ignore+,numeric+,reverse-, >> "temp_median2") if (NregOdd==yes) { centreIndex = int( (Nreal/2.0)+0.5 ) } else { centreIndex = int( (Nreal/2.0) ) } flist3 = "temp_median2" p=0 centreVal = 0 medSky = 0 while ( p <= centreIndex) { jnkstr = fscan(flist3, centreVal) p = p+1 } medSky=centreVal if (NregOdd==no) { jnkstr = fscan(flist3, centreVal) medSky=(medSky+centreVal)/2.0 } print("MEDIAN RMS = "//medSky) print (medSky, >> "temp_tnse8") # Obtaining an average: (NO LONGER USED) ------------ # fields(files="temp_tnse5",fields="4",lines="1-",quit_if-,print_f-) | average(option="new_sample", >> "temp_tnse7") # fields(files="temp_tnse7",fields=1,lines="1-",quit_if-,print_f-, >> "temp_tnse8") # MEAN SKY VALUES: ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if (inclMeans==yes) { # Obtaining a median: ----------------- # delete (files="temp_median*",go_ahea+,ver-) # fields(files="temp_tnse5",fields="3",lines="1-",quit_if-,print_f-, >> "temp_median1") # sort(input_fi="temp_median1",column=1,ignore+,numeric+,reverse-, >> "temp_median2") # if (NregOdd==yes) { # centreIndex = int( (Nreal/2.0)+0.5 ) # } else { # centreIndex = int( (Nreal/2.0) ) # } # flist3 = "temp_median2" # p=0 # centreVal = 0 # medSky = 0 # while ( p <= centreIndex) { # jnkstr = fscan(flist3, centreVal) # p = p+1 # } # medSky=centreVal # if (NregOdd==no) { # jnkstr = fscan(flist3, centreVal) # medSky=(medSky+centreVal)/2.0 # } # print("MEDIAN SKY = "//medSky) # print (medSky, >> "temp_tnse8a") # Obtaining an average: (NO LONGER USED) ------------ fields(files="temp_tnse5",fields="3",lines="1-",quit_if-,print_f-) | average(option="new_sample", >> "temp_tnse7a") fields(files="temp_tnse7a",fields=1,lines="1-",quit_if-,print_f-, >> "temp_tnse8a") } # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # if the mean sky values are included they are in files with an "a" if (inclMeans==no) { joinlines(list1="temp_tnse9,temp_tnse8",list2="",output="temp_tnse10", delim="",missing="Missing",maxchar=161,shortes+,verb-) } if (inclMeans==yes) { joinlines(list1="temp_tnse9,temp_tnse8a,temp_tnse8",list2="", output="temp_tnse10", delim="",missing="Missing",maxchar=161,shortes+,verb-) delete (files="temp_tnse7a",go_ahea+,ver-) delete (files="temp_tnse8a",go_ahea+,ver-) } delete (files="temp_tnse4",go_ahea+,ver-) delete (files="temp_tnse5",go_ahea+,ver-) delete (files="temp_tnse7",go_ahea+,ver-) delete (files="temp_tnse8",go_ahea+,ver-) delete (files="temp_tnse9",go_ahea+,ver-) } # ------------------------------------------------------------------ # Finding average of RMS values across all slices: padding = "" if (inclMeans==no) { fields(files="temp_tnse10",fields="2",lines="1-",quit_if-,print_f-) | average(option="new_sample", >> "temp_tnse11") } if (inclMeans==yes) { print("A") fields(files="temp_tnse10",fields="3",lines="1-",quit_if-,print_f-) | average(option="new_sample", >> "temp_tnse11") padding = "--- " } fields(files="temp_tnse11",fields="1",lines="1-",quit_if-,print_f-, >> "temp_tnse12") flist2 = "temp_tnse12" imgjnk = fscan(flist2, reg) print ("avg_RMS "//padding//reg, >> "temp_tnse13") !cat temp_tnse10 temp_tnse13 > temp_tnse14 delete (files="temp_tnse10",go_ahea+,ver-) delete (files="temp_tnse11",go_ahea+,ver-) delete (files="temp_tnse12",go_ahea+,ver-) delete (files="temp_tnse13",go_ahea+,ver-) rename (files="temp_tnse14",newname="temp_tnse10",field="all") print ("------------------------------") type(input_fi="temp_tnse10",map_cc+,device="terminal",flpar-) print ("------------------------------") print(" ") fields(files="temp_tnse10",fields="3",lines="1-",quit_if-,print_f-, >> "temp_tnse12") # Writing file to output (including footers): ----------------------- type(input_fi="temp_tnse10",map_cc+,device="terminal",flpar-, >> outRMS) if (inclMeans==no) { print ("# IMAGE MEAN STDDEV", >> outRMS) print ("# ("//Nreg//" REGIONS)", >> outRMS) } if (inclMeans==yes) { print ("# IMAGE MEAN SKY MEAN STDDEV", >> outRMS) print ("# ("//Nreg//" REGIONS)", >> outRMS) } # Calculating flux thresholds: ------------------------------------ if (fluxThresh==yes) { if (inclMeans==yes) { # Extracting only the image and RMS columns if mean sky was incl: fields(files="temp_tnse10",fields="1,3",lines="1-",quit_if-,print_f-, >> "temp_tnse11") delete (files="temp_tnse10",go_ahea+,ver-) rename (files="temp_tnse11",newname="temp_tnse10",field="all") } flist1 = "temp_tnse12" print (" "//SNlim//" ", > "temp_tnse13") while ( fscan(flist1, rms) != EOF) { printf ("%6.2f \n",rms, >> "temp_tnse13") } # i.e. formatting mean RMS values in temp_tnse12 into temp_tnse13 # Scan through formatted RMS values: flist2 = "temp_tnse13" outStr = "" k=0 m=j+1 while ( fscan(flist2, rms) != EOF) { k=k+1 if (k==2) { outStr = outStr//" " } if (k>m) { outStr=outStr//" " } outStr=outStr//rms//" " } print (outStr) print (outStr, >> outTh) area1 = pie*(ap1/2)*(ap1/2) area2 = pie*(ap2/2)*(ap2/2) area3 = pie*(ap3/2)*(ap3/2) area4 = pie*(ap4/2)*(ap4/2) print (ap1//" "//area1, >> "temp_tnse4") print (ap2//" "//area2, >> "temp_tnse4") print (ap3//" "//area3, >> "temp_tnse4") print (ap4//" "//area4, >> "temp_tnse4") # Scan through aperture values: ================================ flist1 = "temp_tnse4" while ( fscan(flist1, ap, area) != EOF) { # Scan through RMS sky values: flist2 = "temp_tnse10" while ( fscan(flist2, imgjnk, sig) != EOF) { numer = 1+ sqrt(1+area*(2*sig*gain/SNlim)*(2*sig*gain/SNlim)) denom = 2*gain/(SNlim*SNlim) Flim = numer/denom printf ("%6.2f \n",Flim, >> "temp_tnse5") } # Scan through formatted flux threshold values: flist2 = "temp_tnse5" outStr = ap//" " k=0 while ( fscan(flist2, flx) != EOF) { k=k+1 if (k>j) { # Adding a bit of extra space between slice 10 and mean outStr=outStr//" " } outStr=outStr//flx//" " } print (outStr) print (outStr, >> outTh) delete (files="temp_tnse5",go_ahea+,ver-) } delete (files="temp_tnse4",go_ahea+,ver-) } # =================================================================== # Cleaning up any files used: delete (files="temp_tnse*",go_ahea+,ver-) print(" ") print(" ") print(" All done.") print(" ") end