# Copyright(c) 2001-2004 Association of Universities for Research in Astronomy, Inc. procedure gdark(inimages,outdark) # # Takes GMOS dark images, RAW or gprepared, 3 or 6 amps and combines # to obtain a final dark (overscan and bias corrected), with # optional trimming for image mode reduction, and VAR/DQ output # # Calls gprepare (if necessary); colbias for overscan correction of # the individual images (may be changed to correct final dark if # overscan found to be stable); subtracts zero level image # # Updates GAIN and RDNOISE in the final (combined) image (SCI exts) # # Assumes that DARK frames do not contain an MDF # # # Date Version Author Comments # # 2000/12/18 1.0 cwinge created # 2001/04/23 1.1 cwinge Revised and included NSCIEXT # 2001/07/08 1.2 cwinge some more revision # 2001/07/09 1.3 ij getsec, syntax errors # 2001/07/11 1.4 cwinge # 2001/07/12 1.5 ij,cwinge merged version # 2001/07/13 1.6 cwinge gimverify, use of sed for consistency # with gprepare # July 13, 2001 IJ refixed the bug for one input image # tmpdq* and dqlist* deleted, fl_vardq+ and one input image fixed # cleaned ERROR/WARNING msg, cleaned extra {} etc. # TRIMSEC checked correctly for VAR & DQ and one image # 2001/07/15 1.6 cwinge fl_over- will cause imarith error. Corrected # July 19,2001 IJ # fxinsert with groups="" for VAR/DQ to make it work on new lab data # 2001/07/24 1.7 cwinge removed tsec from parameters, hardwired # for the moment, the task exits if fl_trim=yes # added fl_keepgprep and fl_ovstrim, reset fl_vardq # "mask size" and "nsci_img undefined" errors corrected # clean tmp files if exiting with error # gimverify and gemhedit # 2001/07/25 1.7 cwinge typo fixed for fl_vardq+ and no VAR/DQ in gprepd files # July 30, 2001 IJ fixed gprepare fl_vardq interaction l.405 # (warning msg is now in correct, but the code is doing the right thing to raw images) # # WARNING: Under STSDAS v2.2, imcalc operating over a MEF will create a # MEF, so further operations even on temporary images must # specify explicitly the extension (even if there is only one). # To avoid conflicts, this version uses imutil.imexpr # # Aug 25, 2003 KL IRAF2.12 - new/modified parameters # hedit: addonly # imcombine: headers,bpmasks,expmasks,outlimits # rejmask->rejmasks, plfile->nrejmasks string inimages {"",prompt="Input GMOS images or list"} string outdark {"",prompt="Output dark image"} string bias {"",prompt="Bias (zero level) calibration image"} string logfile {"",prompt="Logfile"} bool fl_over {yes,prompt="Aplly overscan strip correction?"} bool fl_bias {yes,prompt="Aplly zero level (bias) correction?"} string kw_osec {"BIASSEC",prompt="Header keyword for overscan strip image section."} string kw_dsec {"DATASEC",prompt="Header keyword for data section (excludes the overscan)"} bool fl_trim {no,prompt="Trim non-illuminated frame borders (Image mode only)?"} bool fl_ovstrim {no,prompt="Trim overscan section"} bool fl_keepgprep {no,prompt="Keep gprepared images"} #string tsec_1 {"",prompt="Optional data section for amplifier 1"} #string tsec_2 {"",prompt="Optional data section for amplifier 2"} #string tsec_3 {"",prompt="Optional data section for amplifier 3"} #string tsec_4 {"",prompt="Optional data section for amplifier 4"} #string tsec_5 {"",prompt="Optional data section for amplifier 5"} #string tsec_6 {"",prompt="Optional data section for amplifier 6"} string kw_ron {"RDNOISE",prompt="Header keyword for readout noise"} string kw_gain {"GAIN",prompt="Header keyword for gain (e-/ADU"} real ron {0.,prompt="Readout noise value to use if keyword not found"} real gain {1.,prompt="Gain value to use if keyword not found\n"} # Gprepare bool fl_vardq {no,prompt="Create variance and data quality frames?"} string sci_ext {"SCI",prompt="Name of science extension"} string var_ext {"VAR",prompt="Name of variance extension"} string dq_ext {"DQ",prompt="Name of data quality extension"} string gp_bpm {"",prompt="Bad Pixel Mask filename"} int gp_sat {65000,prompt="Saturation level in raw images\n"} # Colbias bool ovs_fliter {no,prompt="Interactive overscan fitting?"} bool ovs_med {no,prompt="Use median instead of average in column bias?"} string ovs_func {"chebyshev",prompt="Overscan fitting function."} int ovs_order {10,prompt="Order of overscan fitting function."} real ovs_lowr {3.,prompt="Low sigma rejection factor."} real ovs_highr {3.,prompt="High sigma rejection factor."} int ovs_niter {0,prompt="Number of rejection iterations.\n"} # Imcombine string im_combine {"average",prompt="Type of combination operation"} string im_reject {"minmax",prompt="Type of rejection algorithm"} real im_lowthresh {INDEF,prompt="Lower threshold"} real im_highthresh {INDEF,prompt="Upper threshold"} string im_scale {"exposure",prompt="Image scaling"} string im_zero {"none",prompt="Image zero point offset"} string im_weight {"none",prompt="Image weights"} string im_statsec {"[*,*]",prompt="Image region for computing statistics"} string im_expname {"EXPTIME",prompt="Header keyword for exposure time"} int im_nlow {0,min=0,prompt="minmax: Number of low pixels to reject"} int im_nhigh {1,min=1,prompt="minmax: Number of high pixels to reject"} int im_nkeep {1,min=0,prompt="Minimum to keep or maximum to reject\n"} bool im_mclip {yes,prompt="Use median in sigma clipping algorithms?"} real im_lsigma {3.,prompt="Lower sigma clipping factor"} real im_hsigma {3.,prompt="Upper sigma clipping factor"} string im_snoise {"0.0",prompt="ccdclip: Sensitivity noise (electrons)"} real im_sigscale {0.1,prompt="Tolerance for sigma clipping scaling correction"} real im_pclip {-0.5,prompt="pclip: Percentile clipping parameter"} real im_grow {0.0,prompt="Radius (pixels) for neighbor rejection\n"} bool verbose {no,prompt="Verbose output?"} int status {0,prompt="Exit status (0=good)"} struct* flist {"",prompt="Internal use only"} begin # Local variables: string l_inimages, l_outdark, l_bias, l_logfile, l_kw_osec, l_kw_dsec string l_tsec[6], l_sci_ext string l_var_ext, l_dq_ext, l_gp_bpm, l_kw_ron, l_kw_gain, l_ovs_func string l_im_combine string l_im_reject, l_im_scale, l_im_zero, l_im_weight, l_im_statsec string l_im_expname, l_im_snoise bool l_fl_over, l_fl_bias, l_fl_trim, l_fl_vardq, l_ovs_fliter, l_ovs_med bool l_im_mclip bool l_verbose, l_fl_keepgprep, l_fl_ovstrim int l_gp_sat, l_ovs_order, l_ovs_niter, l_im_nlow, l_im_nhigh, l_im_nkeep real l_ron, l_gain, l_ovs_lowr, l_ovs_highr, l_im_lowthresh,l_im_highthresh real l_im_lsigma, l_im_hsigma, l_im_sigscale, l_im_pclip, l_im_grow struct l_struct # Runtime variables: string tmplist, inlist, suf, img, inimg[200], test, stsec, dtsec string tmpsci[200], tmpdq string biaslog, dqlist, tmplog, sciout[50], new_reject,dqout[50] string varout[50], tsec string fxsexp, fxvexp, fxdexp, oversec, tmpbias, tmpdqb string tmpdq2, dqexp string tmpdq0, outpref bool fl_gprep[200], fl_fitover[200], fl_biassub[200] int maxfiles, maxscext, ninp, ngp, len, i, j, k, n, nsci_old, nsci_img int next_out, varpos, dqpos, n_var, n_dq, temp,nsci_bias, xval, yval real rms, var, gaineff, roneff, xptime, temp2 struct strout # Temporaries tmplist = mktemp("inlist") dqlist = mktemp("tmpdqlist") biaslog = mktemp("tmpbiaslog") tmpdq = mktemp("tmpdq") tmpdq2 = mktemp("tmpdq2") tmpdq0 = mktemp("tmpdq0") tmplog = mktemp("tmplog") # Set values: l_inimages=inimages; l_outdark=outdark; l_bias=bias; l_logfile=logfile l_kw_osec=kw_osec l_kw_dsec=kw_dsec; l_sci_ext=sci_ext; l_var_ext=var_ext; l_dq_ext=dq_ext l_gp_bpm=gp_bpm; l_kw_ron=kw_ron; l_kw_gain=kw_gain; l_ovs_func=ovs_func l_im_combine=im_combine; l_im_reject=im_reject; l_im_scale=im_scale l_im_zero=im_zero l_im_weight=im_weight; l_im_expname=im_expname; l_im_snoise=im_snoise l_im_statsec=im_statsec l_fl_over=fl_over; l_fl_bias=fl_bias; l_fl_trim=fl_trim l_fl_vardq=fl_vardq l_ovs_fliter=ovs_fliter; l_ovs_med=ovs_med; l_im_mclip=im_mclip l_verbose=verbose l_fl_keepgprep = fl_keepgprep; l_fl_ovstrim=fl_ovstrim l_gp_sat=gp_sat; l_ovs_order=ovs_order; l_ovs_niter=ovs_niter l_im_nlow=im_nlow l_im_nhigh=im_nhigh; l_im_nkeep=im_nkeep l_ron=ron; l_gain=gain; l_ovs_lowr=ovs_lowr; l_ovs_highr=ovs_highr l_im_lowthresh=im_lowthresh; l_im_highthresh=im_highthresh l_im_lsigma=im_lsigma; l_im_hsigma=im_hsigma l_im_sigscale=im_sigscale; l_im_pclip=im_pclip; l_im_grow=im_grow status = 0 maxfiles = 200 maxscext = 50 ngp = 0 nsci_img = 0 ninp = 1 fl_gprep[1] = no # Define the tsec's, for now keep the whole image l_tsec[1]="[*,*]"; l_tsec[2]="[*,*]"; l_tsec[3]="[*,*]" l_tsec[4]="[*,*]"; l_tsec[5]="[*,*]"; l_tsec[6]="[*,*]" cache("imgets") # Tests the logfile: if (l_logfile == "" || stridx(" ",l_logfile)>0) { l_logfile = gmos.logfile if (l_logfile == "" || stridx(" ",l_logfile)>0) { l_logfile = "gmos.log" printlog("WARNING - GDARK: both gdark.logfile and gmos.logfile are empty.", logfile=l_logfile, verbose=l_verbose) printlog(" Using default file gmos.log.", logfile=l_logfile, verbose=yes) } } # Now start logging: (there are 76 dashes) date | scan(l_struct) printlog("----------------------------------------------------------------------------", logfile=l_logfile, verbose=l_verbose) printlog("GDARK -- "//l_struct, logfile=l_logfile, verbose=l_verbose) printlog("", logfile=l_logfile, verbose=l_verbose) printlog("Darklist = "//l_inimages, logfile=l_logfile, verbose=l_verbose) printlog("Output Dark = "//l_outdark,logfile=l_logfile, verbose=l_verbose) printlog("Bias image = "//l_bias,logfile=l_logfile, verbose=l_verbose) printlog(" ",logfile=l_logfile, verbose=l_verbose) # Tests the SCI extension name: if (l_sci_ext == "" || stridx(" ",l_sci_ext)>0) { printlog("ERROR - GDARK: Science extension name parameter sci_ext is missing", logfile=l_logfile, verbose=l_verbose) goto crash } # Tests the VAR and DQ extension names, if required: if (l_fl_vardq) { if (l_var_ext == "" || stridx(" ",l_var_ext)>0 || l_dq_ext == "" || stridx(" ",l_dq_ext)>0) { printlog("ERROR - GDARK: Variance and/or data quality extension name (var_ext/dq_ext) missing", logfile=l_logfile, verbose=l_verbose) goto crash } } # Be sure RDNOISE and GAIN keywords are defined (we need them) if (l_kw_gain == "" || stridx(" ",l_kw_gain)>0) { printlog("ERROR - GDARK: The image gain keyword kw_gain is missing", logfile=l_logfile, verbose=l_verbose) goto crash } if (l_kw_ron == "" || stridx(" ",l_kw_ron)>0) { printlog("ERROR - GDARK: The image readout noise keyword kw_ron is missing", logfile=l_logfile, verbose=l_verbose) goto crash } # Test the data sections - keyword can be null if image trimming selected: # Trim secs now hardwired so we don't need the test # But for the moment, exits if fl_trim=yes if (l_fl_trim) { printlog("ERROR - GDARK: No image mode trimming allowed for the moment. ", logfile=l_logfile, verbose=l_verbose) goto crash } #if (l_fl_trim) { # if(l_tsec[1] == "" || stridx(" ",l_tsec[1])>0) { # printlog("ERROR - GDARK: The image mode trimming is selected but ", logfile=l_logfile, verbose=l_verbose) # printlog(" the data sections are not defined", logfile=l_logfile, verbose=yes) # goto crash # } #} else if (l_kw_dsec =="" || stridx(" ",l_kw_dsec)>0) { printlog("ERROR - GDARK: Data section keyword is missing and", logfile=l_logfile, verbose=l_verbose) printlog(" image mode trimming is not selected.", logfile=l_logfile, verbose=yes) goto crash } # If overscan subtracting, test the overscan section if (l_fl_over) { if (l_kw_osec == "" || stridx(" ",l_kw_osec) > 0) { printlog("ERROR - GDARK: The overscan subtraction is selected but ", logfile=l_logfile, verbose=l_verbose) printlog(" the overscan section keyword is missing", logfile=l_logfile, verbose=yes) goto crash } } # If bias subtracting, test the bias image if (l_fl_bias) { if (l_bias == "" || stridx(" ",l_bias)>0) { printlog("ERROR - GDARK: Zero level image is an empty string.", logfile=l_logfile,verbose=l_verbose) goto crash } if (!imaccess(l_bias) ) { printlog("ERROR - GDARK: Zero level image "//l_bias//" not found.", logfile=l_logfile,verbose=l_verbose) goto crash } # Test the Bias image: gimverify (l_bias) if (gimverify.status == 0) { # Must be .fits (gimverify strips the .fits if it was there, so we # have to put it back) l_bias = gimverify.outname//".fits" } else if (gimverify.status == 1) { printlog("ERROR - GDARK: Bias image "//l_bias//" does not exist", logfile=l_logfile,verbose=yes) goto crash } else { printlog("ERROR - GDARK: Bias image "//l_bias//" is not MEF", logfile=l_logfile,verbose=yes) goto crash } if (l_fl_vardq) { if (!imaccess(l_bias//"["//l_var_ext//"]") || !imaccess(l_bias//"["//l_dq_ext//"]") ) { printlog("WARNING - GDARK: Could not access VAR or DQ extensions for bias image." , logfile=l_logfile, verbose=l_verbose) printlog(" Setting fl_vardq=no", logfile=l_logfile,verbose=yes) l_fl_vardq = no } } # And get the number of bias extensions imgets (l_bias//"[0]", "NSCIEXT", >>& "dev$null") if (imgets.value == "0") { printlog ("ERROR - GDARK: Could not access NSCIEXT keyword for bias image.", logfile=l_logfile, verbose=l_verbose) goto crash } nsci_bias = int(imgets.value) } # Load up input name list inlist=l_inimages # Empty string is an error if (inlist == "" || stridx(" ",inlist)>0) { printlog("ERROR - GDARK: Input file not specified.", logfile=l_logfile, verbose=l_verbose) goto crash } # Test for @. Changed the rest as done for gprepare (avoids the use of sed) if (substr(inlist,1,1) == "@" && !access(substr(inlist,2,strlen(inlist))) ) { printlog("ERROR - GDARK: The input list "//substr(inlist,2,len)//" does not exist", logfile=l_logfile, verbose=l_verbose) goto crash } files(inlist,sort-, > tmplist) flist=tmplist # If keeping the gprepared images, default outpref to "g"; otherwise set # to "tmp" if (l_fl_keepgprep) outpref="g" else outpref="tmp" # Now we test the images: ninp = 0 while(fscan(flist, img) != EOF) { len = strlen(img) suf = substr(img,len-3,len) # Image must exist, but cannot be GEIS or OIF # changed this to use gimverify, per consistency with gprepare gimverify (img) if (gimverify.status == 0) { # Must be .fits (gimverify strips the .fits if it was there, so # we have to put it back) img = gimverify.outname//".fits" ninp = ninp+1 if (ninp > maxfiles) { printlog("ERROR - GDARK: Maximum number of input images exceeded", logfile=l_logfile, verbose=l_verbose) goto crash } inimg[ninp] = img fl_gprep[ninp] = no } else if (gimverify.status == 1) { printlog("ERROR - GDARK: Input image "//img//" does not exist", logfile=l_logfile,verbose=l_verbose) goto crash } else { printlog("ERROR - GDARK: Input image "//img//" is not MEF", logfile=l_logfile, verbose=l_verbose) goto crash } # The first check relates to the current test data do not having # EXTNAME and EXTVER correctly set. # Gprepare takes care of it and adds the NSCIEXT, which we query here. # Check for VAR and DQ extensions in case image already gprepared. # If not found, issues a warning and disables VAR/DQ creation imgets (img//"[0]", "NSCIEXT", >>& "dev$null") # Changed this, the warning msg is not always correct, but at least # it is doing the right thing w/ raw images, IJ # if (!l_fl_vardq && imgets.value == "0") { if (imgets.value == "0") { printlog ("WARNING - GDARK: Image "//img//" has not been gprepared.", logfile=l_logfile, verbose=l_verbose) printlog (" Flagging for GPREPARE.", logfile=l_logfile, verbose=yes) fl_gprep[ninp] = yes ngp = ngp + 1 } if (l_fl_vardq && !fl_gprep[ninp]) { if (!imaccess(img//"["//l_var_ext//"]") || !imaccess(img//"["//l_dq_ext//"]")) { printlog("WARNING - GDARK: Could not access Variance and/or Data Quality extensions for "//img, logfile=l_logfile, verbose=l_verbose) printlog(" Setting fl_vardq = no", logfile=l_logfile, verbose=yes) l_fl_vardq = no } } } # end of while(img) loop flist="" delete(tmplist,ver-, >>& "dev$null") # Call gprepare if needed: if (ngp >0) { i = 1 while (i <= ninp) { if (fl_gprep[i]) { print(inimg[i], >> tmplist) inimg[i] = outpref//inimg[i] } i = i+1 } flist = tmplist printlog("GDARK: calling GPREPARE.",logfile=l_logfile,verbose=l_verbose) while (fscan(flist, img) != EOF) { gprepare (img,outimages="",outpref=outpref,bpm=l_gp_bpm, fl_addmdf=no,kw_mdf="MDFFILE",mdffile="",logfile=l_logfile, fl_vardq=l_fl_vardq,sci_ext=l_sci_ext,var_ext=l_var_ext, dq_ext=l_dq_ext,kw_ron=l_kw_ron,kw_gain=l_kw_gain,ron=l_ron, gain=l_gain,sat=l_gp_sat,verbose=l_verbose,status=0) if (gprepare.status > 0) { if (!imaccess(outpref//img) ) { printlog("ERROR - GDARK: "//outpref//img//" does not exist.", logfile=l_logfile, verbose=l_verbose) goto crash } else printlog("WARNING - GDARK: Using a pre-existing GPREPAREd image", logfile=l_logfile, verbose=l_verbose) } } } flist="" delete(tmplist,ver-, >>& "dev$null") # Tests the output file: if (l_outdark == "" || stridx(" ",l_outdark)>0) { printlog("ERROR - GDARK: The output filename outdark is not defined", logfile=l_logfile,verbose=l_verbose) goto crash } # Tests if it already exists: if (imaccess(l_outdark)) { printlog("ERROR - GDARK: Output file "//l_outdark//" already exists.", logfile=l_logfile,verbose=l_verbose) goto crash } # Must be .fits len = strlen(l_outdark) suf = substr(l_outdark,len-4,len) if (suf !=".fits" ) l_outdark = l_outdark//".fits" # Check the number of extensions and the data sections # Check the need for overscan correction # Check the need for bias subtraction i = 1 nsci_old = 0 while (i <= ninp) { img = inimg[i] nsci_img = 0 # Gets the number of SCI extensions from the header imgets(img//"[0]","NSCIEXT", >>& "dev$null") nsci_img = int(imgets.value) if (nsci_img > maxscext) { printlog("ERROR - GDARK: Maximum number of "//l_sci_ext//" extensions exceeded.", logfile=l_logfile, verbose=l_verbose) goto crash } if (nsci_old!= 0 && nsci_img != nsci_old) { printlog("ERROR - GDARK: The input images do not have all the same number of Science extensions.", logfile=l_logfile, verbose=l_verbose) goto crash } nsci_old = nsci_img j = 1 if (!l_fl_trim) { while (j<=nsci_img) { imgets(img//"["//l_sci_ext//","//j//"]",l_kw_dsec, >>& "dev$null") if (imgets.value == "0") { printlog("ERROR - GDARK: The data section keyword for image "//img//"["//l_sci_ext//","//j//"] is missing.", logfile=l_logfile, verbose=l_verbose) goto crash } j = j+1 } } # Check if the bias and images have the same number of SCI extensions if (l_fl_bias) { if (nsci_img !=nsci_bias) { printlog("ERROR - GDARK: The zero and image frames do not have the same number of extensions", logfile=l_logfile, verbose=l_verbose) goto crash } } # Overscan subtraction flag: imgets (img//"[0]","OVERSCAN", >>& "dev$null") if (!l_fl_over || imgets.value != "0") { printlog("WARNING - GDARK: No overscan correction selected or the image has", logfile=l_logfile, verbose=l_verbose) printlog(" already been corrected for it.", logfile=l_logfile, verbose=yes) fl_fitover[i]=no } else fl_fitover[i]=yes # Bias subtraction flag: imgets (img//"[0]","BIASSUB", >>& "dev$null") if (!l_fl_bias || imgets.value != "0") { printlog("WARNING - GDARK: No bias subtraction selected or the image has", logfile=l_logfile, verbose=l_verbose) printlog(" already been corrected for it.", logfile=l_logfile, verbose=yes) fl_biassub[i]=no } else fl_biassub[i]=yes i = i + 1 } # End of while(i) loop # Now for the main loop: j = 1 n = 1 dqexp = "((a==b) ?1 : 0)" while (j <= nsci_img) { stsec = "["//l_sci_ext//","//j//"]" dtsec = "["//l_dq_ext//","//j//"]" i = 1 gaineff = 0.0 roneff = 0.0 while (i <= ninp) { # Create tmp FITS file names used within this loop tmpbias = mktemp("tmpbias") tmpdqb = mktemp("tmpdqb") img = inimg[i] # If not trimming the overscan section, tsec is the full image # and gsetsec is disabled if (!l_fl_ovstrim) { imgets (img//stsec,"i_naxis1", >>& "dev$null") xval = int(imgets.value) imgets (img//stsec,"i_naxis2", >>& "dev$null") yval = int(imgets.value) tsec = "[1:"//str(xval)//",1:"//str(yval)//"]" # Otherwise, if not trimming the image borders, use DATASEC } else if (!l_fl_trim) { imgets(img//stsec, l_kw_dsec, >>& "dev$null") tsec = imgets.value } else { tsec=l_tsec[j] } # Get the exposure time from the PHU imgets (img//"[0]",l_im_expname, >>& "dev$null") xptime = real(imgets.value) if (xptime == 0. ) { printlog("ERROR - GDARK: image "//img//" exposure time is zero or the", logfile=l_logfile, verbose=l_verbose) printlog(" "//l_expname//" keyword is missing in the image header.", logfile=l_logfile, verbose=yes) goto crash } tmpsci[i] = mktemp("tmpsci") # Overscan fit and subtraction: if (fl_fitover[i]) { imgets(img//stsec, l_kw_osec, >>& "dev$null") oversec = imgets.value colbias (img//stsec,tmpsci[i],bias=oversec,trim=tsec, median=l_ovs_med,interactive=l_ovs_fliter, function=l_ovs_func,order=l_ovs_order,low_reject=l_ovs_lowr, high_reject=l_ovs_highr,niterate=l_ovs_niter, logfiles=biaslog,graphics="stdgraph",cursor="") if (l_verbose) type (biaslog, >> l_logfile) if (i == 1) { type (biaslog) | match ("RMS","STDIN",stop-) | scan (strout) len = strlen(strout) rms = real(substr(strout,len-6,len)) var = rms*rms date | scan(l_struct) hedit(tmpsci[i],"OVERFIT", "Overscan section "//oversec// " fit rms "//rms,add+, addonly-,delete-,version-,show-,update+) } delete (biaslog, ver-, >>& "dev$null") } else imcopy (img//stsec//tsec, tmpsci[i], ver-) hedit (tmpsci[i],"TRIMSEC",tsec,add+,addonly-,delete-,verify-, show-,update+) hedit (tmpsci[i],l_im_expname,xptime,add+,addonly-,delete-,verify-, show-,update+) # Create the BPM for each extension/image if (l_fl_vardq) { imgets(img//dtsec,"TRIMSEC", >>& "dev$null") if (imgets.value == "0") { imcopy(img//dtsec//tsec,tmpdq//"_"//n//".pl", verbose=no) hedit(tmpdq//"_"//n//".pl","TRIMSEC",tsec,add+,addonly-, delete-,verify-, show-,update+) } else { test = imgets.value if (test != tsec) { printlog("ERROR - GDARK: The DQ plane for extension "//j//" has been trimmed", logfile=l_logfile, verbose=l_verbose) printlog(" to a different section than the science extension.", logfile=l_logfile, verbose=yes) goto crash } else { printlog("WARNING - GDARK: The DQ plane for extension "//j//" has already", logfile=l_logfile, verbose=l_verbose) printlog(" been trimmed. Trim section not applied.", logfile=l_logfile,verbose=yes) imcopy(img//dtsec,tmpdq//"_"//n//".pl", verbose=no) } } print (tmpdq//"_"//n//".pl" , >> dqlist) hedit(tmpsci[i],"BPM",tmpdq//"_"//n//".pl",add+,addonly-, delete-,verify-,show-,update+) } # Subtract the bias if (fl_biassub[i]) { imgets(l_bias//stsec,"TRIMSEC", >>& "dev$null") if (imgets.value == "0") { imcopy(l_bias//stsec//tsec,tmpbias,verbose=no) if (l_fl_vardq) imcopy(l_bias//dtsec//tsec,tmpdqb, verbose=no) } else { test = imgets.value if (test != tsec) { printlog("ERROR - GDARK: The bias image for extension "//j//" has been trimmed", logfile=l_logfile, verbose=l_verbose) printlog(" to a different section than the dark image.", logfile=l_logfile, verbose=yes) goto crash } else { printlog("WARNING - GDARK: The bias image for extension "//j//" has already", logfile=l_logfile,verbose=l_verbose) printlog(" been trimmed. Trim section not applied.", logfile=l_logfile,verbose=yes) imcopy(l_bias//stsec,tmpbias, verbose=no) if (l_fl_vardq) imcopy(l_bias//dtsec,tmpdqb, verbose=no) } } imarith (tmpsci[i], "-", tmpbias, tmpsci[i], ver-) if (i == 1) hedit(tmpsci[i],"BIASIMG","Bias image "//l_bias,add+, addonly-,delete-,verify-,show-,update+) imdelete (tmpbias,ver-, >>& "dev$null") # Combine the dark and bias DQ planes if VAR/DQ is enabled if (l_fl_vardq) { imrename (tmpdq//"_"//n//".pl",tmpdq0, ver-) addmasks(tmpdq0//","//tmpdqb, tmpdq//"_"//n//".pl","im1 || im2", flags="") imdelete (tmpdqb//","//tmpdq0,ver-, >>& "dev$null") } } print (tmpsci[i], >> tmplist) # Calculates the rdnoise and gain: imgets(img//stsec,l_kw_gain, >>& "dev$null") if (imgets.value == "0") { printlog("WARNING - GDARK: The keyword for image gain "//l_kw_gain//" was not found in "//img//dtsec, logfile = l_logfile, verbose=l_verbose) printlog(" Using gain of "//str(l_gain), logfile=l_logfile, verbose=yes) hedit(img//dtsec,l_kw_gain,l_gain,add+,addonly-,delete-,verify-, show-,update+) gaineff = gaineff + l_gain } else gaineff = gaineff + real(imgets.value) imgets(img//stsec, l_kw_ron, >>& "dev$null") if (imgets.value == "0") { printlog("WARNING - GDARK: The keyword for image readout noise "//l_kw_ron//" was not found", logfile=l_logfile, verbose=l_verbose) printlog(" in "//img//dtsec//" Using ron of "//str(l_ron), logfile=l_logfile,verbose=yes) hedit(img//dtsec,l_kw_ron,l_ron, add+,addonly-,delete-,verify-, show-,update+) roneff = sqrt(roneff**2 + l_ron**2) } else roneff = sqrt(roneff**2 + real(imgets.value)**2) i = i+1 n = n+1 # delete tmp files if (imaccess(tmpbias)) imdelete(tmpbias,ver-, >>& "dev$null") if (imaccess(tmpdqb)) imdelete(tmpdqb,ver-, >>& "dev$null") } # End of while(i) loop # Now combine sciout[j] = mktemp("tmpsciout") dqout[j] = mktemp("tmpdqout") varout[j] = mktemp("tmpvarout") new_reject = l_im_reject if (ninp == 1) { printlog("WARNING - GDARK: Only one input image", logfile=l_logfile,verbose=l_verbose) imcopy("@"//tmplist,sciout[j],verbose-) if (l_fl_vardq) { imcopy("@"//dqlist, dqout[j], verbose-) imcopy(img//"["//l_var_ext//","//j//"]"//tsec, varout[j], verbose-) if (fl_fitover[1]) imarith(varout[j],"+",var, varout[j], ver-) if (fl_biassub[1]) { # Need to check sci_ext, TRIMSEC not in other extensions, IJ imgets(l_bias//"["//l_sci_ext//","//j//"]","TRIMSEC", >>& "dev$null") if (imgets.value == "0") test = tsec else test="" imarith (varout[j], "+", l_bias//"["//l_var_ext//","//j//"]"//test, varout[j], ver-) } } } else { if (ninp <=5) { if (ninp ==2) { printlog("WARNING - GDARK: Only two images to combine, turning off rejection.", logfile=l_logfile, verbose=l_verbose) new_reject = "none" } else printlog("WARNING - GDARK: Five or less images to combine", logfile=l_logfile,verbose=l_verbose) } if (l_fl_vardq) { # Combine with scale=exposure then normalize to the EXPTIME in # the output header, which is the weighted average of the # input values imcombine("@"//tmplist,sciout[j],headers="",bpmasks="", rejmasks="",nrejmasks=tmpdq2,expmasks="",sigmas=varout[j], logfile=tmplog,combine=l_im_combine,reject=new_reject, project=no,outtype="real",outlimits="",offsets="none", masktype="goodvalue",maskvalue=0.,blank=0,scale=l_im_scale, zero=l_im_zero,weight=l_im_weight,statsec=l_im_statsec, expname=l_im_expname,lthreshold=l_im_lowthresh, hthreshold=l_im_highthresh,nlow=l_im_nlow,nhigh=l_im_nhigh, nkeep=l_im_nkeep,mclip=l_im_mclip,lsigma=l_im_lsigma, hsigma=l_im_hsigma,rdnoise=l_ron,gain=l_gain, snoise=l_im_snoise,sigscale=l_im_sigscale,pclip=l_im_pclip, grow=l_im_grow) type (tmplog, >> l_logfile) imgets(sciout[j], l_im_expname, >>& "dev$null") xptime = real(imgets.value) imarith (sciout[j], "/", xptime, sciout[j], ver-) hedit (sciout[j], "EXPTIME", "1",add+,addonly-,delete-,verify-, show-,update+) delete (tmplog, ver-, >>& "dev$null") # Calculate the VAR plane from the sigma image imarith(varout[j],"*",varout[j],varout[j], ver-) # Calculate the DQ plane imexpr (dqexp, dqout[j], tmpdq2, ninp, dims="auto", intype="auto", outtype="int",refim="auto",bwidth=0, btype="nearest",bpixval=0.,rangecheck+,verbose-, exprdb="none",lastout="") imdelete (tmpdq2, ver-, >>& "dev$null") } else { imcombine("@"//tmplist,sciout[j],headers="",bpmasks="", rejmasks="",nrejmasks="",expmasks="",sigmas="", logfile=tmplog,combine=l_im_combine,reject=new_reject, project=no,outtype="real",outlimits="",offsets="none", masktype="goodvalue",maskvalue=0.,blank=0,scale=l_im_scale, zero=l_im_zero,weight=l_im_weight,statsec=l_im_statsec, expname=l_im_expname,lthreshold=l_im_lowthresh, hthreshold=l_im_highthresh,nlow=l_im_nlow,nhigh=l_im_nhigh, nkeep=l_im_nkeep,mclip=l_im_mclip,lsigma=l_im_lsigma, hsigma=l_im_hsigma,rdnoise=l_ron,gain=l_gain, snoise=l_im_snoise,sigscale=l_im_sigscale,pclip=l_im_pclip, grow=l_im_grow) type (tmplog, >> l_logfile) imgets (sciout[j],"EXPTIME", >>& "dev$null") xptime = real(imgets.value) imarith (sciout[j], "/", xptime, sciout[j], ver-) hedit (sciout[j],"EXPTIME","1",add+,addonly-,delete-,verify-, show-,update+) delete (tmplog, ver-, >>& "dev$null") } # end of "else" for fl_vardq # Update ReadNoise and Gain if more than one image imgets(sciout[j],"GAINORIG", >>& "dev$null") if (imgets.value == "0") { imgets(sciout[j], l_kw_gain, >>& "dev$null") hedit(sciout[j],"GAINORIG",real(imgets.value),add+,addonly-, delete-,verify-,show-,update+) } imgets(sciout[j],"RONORIG", >>& "dev$null") if (imgets.value == "0") { imgets(sciout[j], l_kw_ron, >>& "dev$null") hedit(sciout[j],"RONORIG",real(imgets.value),add+,addonly-, delete-,verify-,show-,update+) } hedit(sciout[j],l_kw_gain,gaineff,add+,addonly-,delete-,verify-, show-,update+) hedit(sciout[j],l_kw_ron,roneff,add+,addonly-,delete-,verify-, show-,update+) } # end of combining imdelete("@"//tmplist, ver-, >>& "dev$null") if (access(dqlist)) imdelete("@"//dqlist, ver-, >>& "dev$null") delete(tmplist//","//dqlist, ver-, >>& "dev$null") j = j+1 } #End of while(j) loop # Build the expressions for fxinsert next_out = nsci_img fxsexp = sciout[1]//".fits" if (l_fl_vardq) { next_out = 3*nsci_img varpos = nsci_img dqpos = 2*nsci_img fxvexp = varout[1]//".fits" fxdexp = dqout[1]//".fits" } j = 2 while (j<=nsci_img) { fxsexp = fxsexp//","//sciout[j]//".fits" if (l_fl_vardq) { fxvexp = fxvexp//","//varout[j]//".fits" fxdexp = fxdexp//","//dqout[j]//".fits" } j = j+1 } # Creates the image, using PHU from first image on input list fxcopy(inimg[1],l_outdark,groups="0",new_file=yes, verbose=no, >>& "dev$null") # And insert fxinsert(fxsexp,l_outdark//"[0]","",verbose=no, >>& "dev$null") if (l_fl_vardq) { fxinsert(fxvexp,l_outdark//"["//varpos//"]","",verbose=no, >>& "dev$null") fxinsert(fxdexp,l_outdark//"["//dqpos//"]","",verbose=no, >>& "dev$null") } # Update the headers and clean up temporaries j = 1 while(j<=nsci_img) { #hedit(l_outdark//"["//j//"]","NEXTEND",next_out,add-,addonly-,del-,ver-,show-,upd+) hedit(l_outdark//"["//j//"]","EXTNAME",l_sci_ext,add+,addonly-,delete-, verify-,show-,update+) hedit(l_outdark//"["//j//"]","EXTVER",j,add+,addonly-,delete-,verify-, show-,update+) imdelete(sciout[j],ver-, >>& "dev$null") if (l_fl_vardq) { n_var = varpos+j n_dq = dqpos+j #hedit(l_outdark//"["//n_var//"]","NEXTEND",next_out,add-,addonly-,del-,ver-,show-,upd+) hedit(l_outdark//"["//n_var//"]","EXTNAME",l_var_ext,add+,addonly-, delete-,verify-,show-,update+) hedit(l_outdark//"["//n_var//"]","EXTVER",j,add+,addonly-,delete-, verify-,show-,update+) #hedit(l_outdark//"["//n_dq//"]","NEXTEND",next_out,add-,addonly-,del-,ver-,show-,upd+) hedit(l_outdark//"["//n_dq//"]","EXTNAME",l_dq_ext,add+,addonly-, delete-,verify-,show-,update+) hedit(l_outdark//"["//n_dq//"]","EXTVER",j,add+,addonly-,delete-, verify-,show-,update+) imdelete(varout[j]//","//dqout[j],ver-, >>& "dev$null") } j = j+1 } # Set the data section of the output image, IJ # But skip if not trimming the overscan if (l_fl_ovstrim) gsetsec(l_outdark,key_datsec="DATASEC") # Update the PHU date | scan(l_struct) gemhedit(l_outdark//"[0]","NCOMBINE",ninp,"Number of extensions") gemhedit(l_outdark//"[0]","NEXTEND",next_out,"Number of science extensions") gemhedit(l_outdark//"[0]","GDARK",l_struct,"Time stamp for GDARK") gemhedit(l_outdark//"[0]","GEM-TLM",l_struct, "Last modification with GEMINI") imgets(l_outdark//"[0]","OVERSCAN", >>& "dev$null") if (l_fl_over && (imgets.value == "0") ) gemhedit(l_outdark//"[0]","OVERSCAN",l_struct, "Time stamp for overscan subtraction") imgets(l_outdark//"[0]","BIASSUB", >>& "dev$null") if (l_fl_bias && (imgets.value == "0") ) gemhedit(l_outdark//"[0]","BIASSUB",l_struct//" with "//l_bias, "Time stamp and image for bias subtraction") # Clean up: goto clean crash: # Exit with error subroutine status=1 goto clean clean: # wrap up i = 1 while (i <= ninp) { if (!l_fl_keepgprep && fl_gprep[i]) imdelete (inimg[i], ver-, >>& "dev$null") imdelete (tmpsci[i], ver-, >>& "dev$null") i = i+1 } imdelete ("tmpdq_*.pl", ver-, >>& "dev$null") delete (tmplist//","//biaslog//","//dqlist, ver-, >>& "dev$null") flist="" # close log file if (status == 0) test="no" else test="one or more" printlog(" ", logfile=l_logfile, verbose=l_verbose) printlog("GDARK finished with "//test//" errors.", logfile=l_logfile, verbose=l_verbose) #if (ngp > 0) { # printlog("One or more images were processed with GPREPARE and ",logfile=l_logfile, verbose=l_verbose) # printlog("left in the current directory as tmp.", logfile=l_logfile, verbose=l_verbose) #} printlog("----------------------------------------------------------------------------", logfile=l_logfile, verbose=l_verbose) end