# Copyright(c) 2002 Association of Universities for Research in Astronomy, Inc. procedure gireduce(inimages) # Basic reductions of GMOS imaging data # # Version Feb 28, 2002 MT,IJ,BM v1.3 release # Aug 16, 2002 IJ debug version: "Cannot coerce string `' to real" bug (Linux bug) # Aug 26, 2002 IJ parameter encoding # Sept 20, 2002 IJ v1.4 release # Oct 31, 2002 IJ overvar[6,200] bug fixed # Feb 12, 2003 BM Use imgets to read GAIN instead of keypar to avoid problem on Linux # Mar 3, 2003 BM User imgets instead of keypar to read TRIMSEC, was the cause of segmentation fault on some Linux systems char inimages {"",prompt="Input GMOS images"} # OLDP-1 char outpref {"r",prompt="Prefix for output images"} # OLDP-1 char outimages {"",prompt="Output images"} # OLDP-1 bool fl_over {no,prompt="Subtract overscan level"} # OLDP-3 bool fl_trim {yes,prompt="Trim off the overscan section"} # OLDP-2 bool fl_bias {yes,prompt="Subtract bias image"} # OLDP-2 bool fl_dark {no,prompt="Subtract (scaled) dark image"} # OLDP-3 bool fl_flat {yes,prompt="Do flat field correction?"} # OLDP-2 bool fl_vardq {no,prompt="Create variance and data quality frames"} # OLDP-3 bool fl_addmdf {no,prompt="Add Mask Definition File? (LONGSLIT/MOS/IFU modes)"} # OLDP-3 char bias {"",prompt="Bias image name"} # OLDP-1 char dark {"",prompt="Dark image name"} # OLDP-1 char flat1 {"",prompt="Flatfield image 1"} # OLDP-1 char flat2 {"",prompt="Flatfield image 2"} # OLDP-1 char flat3 {"",prompt="Flatfield image 3"} # OLDP-1 char flat4 {"",prompt="Flatfield image 4\n"} # OLDP-1 char key_exptime {"EXPTIME",prompt="Header keyword of exposure time"} # OLDP-3 char key_biassec {"BIASSEC",prompt="Header keyword for bias section"} # OLDP-3 char key_datasec {"DATASEC",prompt="Header keyword for data section"} # OLDP-3 char rawpath {"",prompt="GPREPARE: Path for input raw images"} # OLDP-1 char gp_outpref {"g",prompt="GPREPARE: Prefix for output images"} # OLDP-1 char sci_ext {"SCI",prompt="Name of science extension"} # OLDP-3 char var_ext {"VAR",prompt="Name of variance extension"} # OLDP-3 char dq_ext {"DQ",prompt="Name of data quality extension"} # OLDP-3 char key_mdf {"MASKNAME",prompt="Header keyword for the Mask Definition File"} # OLDP-3 char mdffile {"",prompt="MDF file to use if keyword not found"} # OLDP-3 char mdfdir {"gmos$data/", prompt="MDF database directory"} # OLDP-1 char bpm {"",prompt="Bad pixel mask"} # OLDP-1 char gaindb {"gmos$data/gmosamps.dat",prompt="Database with gain data"} # OLDP-3 int sat {65000,prompt="Saturation level in raw images [ADU]"} # OLDP-3 char key_filter {"FILTER2",prompt="Header keyword of filter"} # OLDP-3 char key_ron {"RDNOISE",prompt="Header keyword for readout noise"} # OLDP-3 char key_gain {"GAIN",prompt="Header keyword for gain (e-/ADU)"} # OLDP-3 real ron {3.5,prompt="Readout noise in electrons"} # OLDP-3 real gain {2.2,prompt="Gain in e-/ADU"} # OLDP-3 bool fl_inter {no,prompt="Interactive overscan fitting?"} # OLDP-3 bool median {no,prompt="Use median instead of average in column bias?"} # OLDP-3 char function {"chebyshev",prompt="Overscan fitting function",enum="spline3|legendre|chebyshev|spline1"} # OLDP-3 int order {1,prompt="Order of overscan fitting function"} # OLDP-3 real low_reject {3.,prompt="Low sigma rejection factor in overscan fit"} # OLDP-3 real high_reject {3.,prompt="High sigma rejection factor in overscan fit"} # OLDP-3 int niterate {2,prompt="Number of rejection iterations in overscan fit"} # OLDP-3 char logfile {"",prompt="Logfile"} # OLDP-1 bool verbose {yes,prompt="Verbose?"} # OLDP-2 int status {0,prompt="Exit status (0=good)"} # OLDP-4 char *scanfile {"",prompt="For internal use only"} # OLDP-4 ################################################################################ begin # Variable declaration char l_inimages, l_outpref, l_outimages, l_sci_ext, l_var_ext, l_dq_ext char l_key_biassec, l_function, l_key_datasec, l_bias, l_dark, l_key_filter, l_key_exptime char l_flat[4], l_gp_outpref, l_bpm, l_key_ron, l_key_gain, l_rawpath char l_key_mdf, l_mdffile, l_mdfdir, l_logfile, l_gaindb bool l_fl_vardq, l_fl_over, l_median, l_fl_inter bool l_fl_trim, l_fl_bias, l_fl_dark, l_fl_flat, l_fl_addmdf, l_verbose int l_order, l_niterate, l_sat real l_low_reject, l_high_reject, l_ron, l_gain # General variables: char img, s_empty, input[200], output[200] char tmpinlist, tmpoutlist, tmplog, tmpfile char obsmodeflat[4] char l_biassec, l_datasec, tsec_img, tsec_bias[6], tsec_flat[6,4], tsec_dark[6] char thisfilter, flatfilter[4] char srms, sorder, subflat, subbias, subdark, st1, st2 # local image names and expressions char ll_bias, ll_flat, ll_dark, ll_sciexp, ll_varexp, ll_dqexp char ll_varbias, ll_varflat, ll_vardark, ll_dqbias, ll_dqflat, ll_dqdark int ll_ndq real l_darkexptime, l_darkscale char s_bias, s_flat, s_dark # for printlog only real n_gain[6], overvar[6,200], dummy, l_gainsl[6], l_gainfl[6], l_gainfh[6] real gain1, gain2, gain3, gain4, gain5, gain6 char l_readmode, l_gainmode, gaindata, tmpmdf int nbad, nfiles, nimages, ii, jj, kk int n_extflat[4], n_extbias, n_ext[200], n_extdark, n_sci int fnum, maxfiles, nx char ccd_flat[4,6], ccd_bias[6], ccd_dark[6], ccd[200,6] bool found1, found2, l_fl_raw bool canover, cantrim, canbias, candark, canflat[4], n_do[200], hasmdf[200] bool fl_tsec_img, fl_tsec_bias[6], fl_tsec_flat[6,4], badflag, fl_tsec_dark[6] struct l_struct ################################################################################ # Get parameters l_inimages=inimages; l_outpref=outpref; l_outimages=outimages l_fl_over=fl_over; l_logfile=logfile; l_verbose=verbose l_fl_bias=fl_bias; l_bias=bias; l_fl_flat=fl_flat; l_key_filter=key_filter l_flat[1]=flat1; l_flat[2]=flat2; l_flat[3]=flat3; l_flat[4]=flat4 l_fl_vardq=fl_vardq; l_sci_ext=sci_ext; l_var_ext=var_ext; l_dq_ext=dq_ext l_key_ron=key_ron; l_key_gain=key_gain ; l_key_biassec=key_biassec l_fl_dark=fl_dark; l_dark=dark ; l_fl_trim = fl_trim l_key_exptime=key_exptime # GPREPARE PARAMETERS: l_gp_outpref=gp_outpref; l_bpm=bpm; l_key_mdf=key_mdf ; l_rawpath=rawpath l_fl_addmdf=fl_addmdf; l_mdffile=mdffile; l_mdfdir=mdfdir l_key_datasec=key_datasec; l_sat=sat ; l_gaindb=gaindb l_ron=ron ; l_gain=gain # COLBIAS PARAMETERS" l_median=median; l_fl_inter=fl_inter; l_function=function l_order=order; l_niterate=niterate; l_low_reject=low_reject; l_high_reject=high_reject ################################################################################ # Initialize variables status=0 maxfiles=200 badflag=no cache("imgets","keypar","gimverify") # Make temp files tmpfile = mktemp("tmpfile") tmpinlist = mktemp("tmpinlist") tmpoutlist = mktemp("tmpoutlist") tmplog = mktemp("tmplog") tmpmdf=mktemp("tmpmdf") ################################################################################ # Define the name of the logfile s_empty=""; print(l_logfile) | scan(s_empty); l_logfile=s_empty if (l_logfile == "") { l_logfile = gmos.logfile if (l_logfile == "") { l_logfile = "gmos.log" printlog("WARNING - GIREDUCE: both gireduce.logfile and \ gmos.logfile are empty.",l_logfile,l_verbose) printlog(" Using default file gmos.log.",l_logfile,l_verbose) } } # Write to the logfile date | scan(l_struct) printlog("-----------------------------------------------------------------\ -----------",l_logfile,l_verbose) printlog("GIREDUCE -- "//l_struct,l_logfile,l_verbose) printlog("",l_logfile,l_verbose) #printlog("Debug version",l_logfile,l_verbose) ################################################################################ # Check the input images and place them in a file. nimages=0 if(substr(l_inimages,1,1)=="@") { img=substr(l_inimages,2,strlen(l_inimages)) # img is the file if(!access(img)) { printlog("ERROR - GIREDUCE: Input file "//img//" not found.",l_logfile,yes) goto crash } } files(l_inimages,sort-, > tmpfile) nbad=0; nfiles=0 scanfile = tmpfile while(fscan(scanfile,img) != EOF) { l_fl_raw=no gimverify(img); img=gimverify.outname if(gimverify.status>=1 && l_rawpath!="") { gimverify(l_rawpath//img) ; l_fl_raw=yes } if(gimverify.status==1) { printlog("ERROR - GIREDUCE: File "//img//" not found.",l_logfile,yes) nbad+=1 } else if(gimverify.status>1) { printlog("ERROR - GIREDUCE: File "//img//" not a MEF FITS image.",l_logfile,yes) nbad+=1 } else { if(l_fl_raw) imgets(l_rawpath//img//"[0]","GPREPARE", >& "dev$null") else imgets(img//"[0]","GPREPARE", >& "dev$null") if(l_fl_raw && imgets.value!="0") { printlog("ERROR - GIREDUCE: Input image "//img//" is in raw data \ directory "//l_rawpath//" and is gprepared",l_logfile,yes) goto crash } if (imgets.value == "0") { # Use already gprepared image if it exists if(access(l_gp_outpref//img//".fits")) { imgets(l_gp_outpref//img//"[0]","GPREPARE", >& "dev$null") if(imgets.value=="0") { printlog("ERROR - GIREDUCE: gprepare outpref \ parameter points to existing un-gprepared file",l_logfile,yes) goto crash } else { printlog("Using already gprepared image "//l_gp_outpref//img, l_logfile,l_verbose) } } else { gprepare(img,rawpath=l_rawpath,outimages="",outpref=l_gp_outpref, bpm=l_bpm,fl_addmdf=l_fl_addmdf,key_mdf=l_key_mdf, gaindb=l_gaindb, mdffile=l_mdffile, mdfdir=l_mdfdir, logfile=l_logfile, fl_vardq=l_fl_vardq,sci_ext=l_sci_ext,var_ext=l_var_ext, dq_ext=l_dq_ext,key_ron=l_key_ron,key_gain=l_key_gain,ron=l_ron,gain=l_gain, key_dsec=l_key_datasec,sat=l_sat,verbose=l_verbose) if(gprepare.status!=0) { printlog("ERROR - GIREDUCE: gprepare returned with error",l_logfile,yes) goto crash } } img = l_gp_outpref//img } nfiles+=1 input[nfiles] = img print(img, >> tmpinlist) } } # end while loop # Is there anything else to do? if (!l_fl_over && !l_fl_trim && !l_fl_bias && !l_fl_dark && !l_fl_flat) { goto clean } # Check for empty file list if(nfiles==0) { printlog("ERROR - GIREDUCE: No input images to process.", l_logfile,yes) nimages=0 goto crash } else if (nfiles>maxfiles) { printlog("ERROR - GIREDUCE: Maximum number of input files is "//maxfiles, l_logfile,yes) goto crash } # Exit if problems found with input files if (nbad > 0) { printlog("ERROR - GIREDUCE: "//nbad//" image(s) either do not exist or \ are not MEF files.",l_logfile, yes) goto crash } # end if (nbad > 0) printlog("",l_logfile,l_verbose) printlog("Input files:",l_logfile,l_verbose) if(l_verbose) type(tmpinlist) type(tmpinlist, >> l_logfile) printlog("",l_logfile,l_verbose) delete(tmpfile, ver-) # Input images are in tmpinlist file and an array input[ii] nimages = nfiles ################################################################################ # Check the output images. Check that they have been defined, and that they # do not exist already. Count the number of output images. If they are not # the same as the input, exit with an error (only check if they are given # as outimage). s_empty=""; print(l_outimages) | scan(s_empty); l_outimages=s_empty ii = 0 if(l_outimages!="") { if(substr(l_outimages,1,1)=="@") { img=substr(l_outimages,2,strlen(l_outimages)) if(!access(img)) { printlog("ERROR - GIREDUCE: Output file list "//img//" not found.", l_logfile,yes) goto crash } } files(l_outimages,sort-, > tmpfile) scanfile = tmpfile while(fscan(scanfile,img) != EOF) { nfiles-=1; ii+=1 gimverify(img); img=gimverify.outname if(gimverify.status!=1) { printlog("WARNING - GIREDUCE: File "//img//" already exists.", l_logfile,l_verbose) n_do[ii] = no nfiles+=1 # add 1 back to nfiles } else { output[ii] = img n_do[ii] = yes print(img, >> tmpoutlist) } } # end while loop scanfile="" delete(tmpfile, ver-) } else if (l_outimages=="") { scanfile=tmpinlist while(fscan(scanfile,img) != EOF) { nfiles-=1; ii+=1 output[ii] = l_outpref//img gimverify(output[ii]); output[ii]=gimverify.outname if(gimverify.status!=1) { printlog("WARNING - GIREDUCE: File "//output[ii]//" already exists", l_logfile,l_verbose) n_do[ii]=no nfiles+=1 # add 1 back to nfiles } else { n_do[ii]=yes print(output[ii], >> tmpoutlist) } } # end while loop scanfile="" } if(nfiles!=0) { printlog("ERROR - GIREDUCE: number of input and output images is not the same", l_logfile,yes) goto crash } printlog("Output files:",l_logfile,l_verbose) if(l_verbose) type(tmpoutlist) type(tmpoutlist, >> l_logfile) printlog("",l_logfile,l_verbose) # Output images are in array output[ii] delete (tmpinlist//","//tmpoutlist,ver-) ################################################################################ # Check the extension names # SCI extension - s_empty=""; print(l_sci_ext) | scan(s_empty); l_sci_ext=s_empty if(l_sci_ext=="") { printlog("ERROR - GIREDUCE: Science extension parameter sci_ext is missing", l_logfile,yes) goto crash } # VAR and DQ extensions - if(l_fl_vardq) { s_empty=""; print(l_var_ext) | scan(s_empty); l_var_ext=s_empty if(l_var_ext=="") { printlog("ERROR - GIREDUCE: Variance extension parameter \ var_ext is missing", l_logfile,yes) goto crash } s_empty=""; print(l_dq_ext) | scan(s_empty); l_dq_ext=s_empty if(l_dq_ext=="") { printlog("ERROR - GIREDUCE: Data quality extension parameter \ dq_ext is missing", l_logfile,yes) goto crash } } ################################################################################ # If flat-fielding: Check that filter key is defined. Check that flat field # images are ok. Get the number of science extensions. # The obsmode for the flat fields is the same as that of the input images. for(ii=1; ii<=4; ii+=1) { canflat[ii]=no } if(l_fl_flat) { # The header keyword s_empty=""; print(l_key_filter) | scan(s_empty); l_key_filter=s_empty if(l_key_filter=="") { printlog("WARNING - GIREDUCE: Cannot flat field correct.",l_logfile,l_verbose) printlog(" Header keyword parameter key_filter \ is not defined",l_logfile,l_verbose) l_fl_flat = no goto noflats } # The flat field images - for(ii=1; ii<=4; ii+=1) { s_empty=""; print(l_flat[ii]) | scan(s_empty); l_flat[ii]=s_empty if(l_flat[ii]!="") { gimverify(l_flat[ii]); l_flat[ii]=gimverify.outname if(gimverify.status==0) { canflat[ii] = yes keypar(l_flat[ii]//"[0]","GIFLAT",silent+) if(!keypar.found) { printlog("ERROR - GIREDUCE: Flat field: "//l_flat[ii]//" was \ not created with giflat task.",l_logfile, yes) badflag=yes } imgets(l_flat[ii]//"[0]","OBSMODE", >& "dev$null") obsmodeflat[ii] = imgets.value imgets(l_flat[ii]//"[0]","NSCIEXT", >& "dev$null") n_extflat[ii] = int(imgets.value) # Get the name of the filter in array flatfilter[4] keypar(l_flat[ii]//"[0]",l_key_filter,silent+) flatfilter[ii] = keypar.value if(!keypar.found) canflat[ii] = no # Place the sizes of the flats in the tsec_flat[6,4] and fl_tsec_flat[6,4] arrays. for(jj=1; jj<=n_extflat[ii]; jj+=1) { keypar(l_flat[ii]//"["//l_sci_ext//","//jj//"]","TRIMSEC",silent+) tsec_flat[jj,ii] = keypar.value; fl_tsec_flat[jj,ii]=keypar.found imgets(l_flat[ii]//"["//jj//"]","DETSEC", >& "dev$null") ccd_flat[ii,jj] = imgets.value imgets(l_flat[ii]//"["//jj//"]","CCDSUM", >& "dev$null") ccd_flat[ii,jj] = ccd_flat[ii,jj]//imgets.value } } else { printlog("ERROR - GIREDUCE: Can not find flat-field \ frame: "//l_flat[ii],l_logfile, yes) badflag=yes } } } # end of checking flat field images noflats: } # end of if(l_fl_flat) ################################################################################ # If overscan subt.: Check that bias section keywords is defined and is ok. if(l_fl_over) { s_empty=""; print(l_key_biassec) | scan(s_empty); l_key_biassec=s_empty if(l_key_biassec=="") { printlog("ERROR - GIREDUCE: bias section keyword (key_biassec) \ is an empty string.",l_logfile,yes) goto crash } } ################################################################################ # If bias subtracting: check that the bias image is defined and is ok. # If fl_vardq=yes, need to check that the bias image has var and dq planes. # Get the number of science extensions. canbias=no if(l_fl_bias) { s_empty=""; print(l_bias) | scan(s_empty); l_bias=s_empty gimverify(l_bias); l_bias=gimverify.outname if(gimverify.status==0) { canbias = yes keypar(l_bias//"[0]","GBIAS",silent+) if(!keypar.found) { printlog("ERROR - GIREDUCE: Bias frame: "//l_bias//" was not \ created with gbias task.",l_logfile, yes) badflag=yes } imgets(l_bias//"[0]","NSCIEXT", >& "dev$null") n_extbias = int(imgets.value) for(ii=1; ii<=n_extbias; ii+=1) { keypar(l_bias//"["//l_sci_ext//","//ii//"]","TRIMSEC",silent+) imgets(l_bias//"["//l_sci_ext//","//ii//"]","DETSEC", >& "dev$null") ccd_bias[ii] = imgets.value imgets(l_bias//"["//l_sci_ext//","//ii//"]","ccdsum", >& "dev$null") ccd_bias[ii] = ccd_bias[ii]//imgets.value } } else { printlog("ERROR - GIREDUCE: Can not find bias frame: "//l_bias, l_logfile,yes) badflag=yes } } ################################################################################ # If dark current correcting: check that the dark image is defined and is ok. # If fl_vardq=yes, need to check that the dark image has var and dq planes. # Get the number of science extensions. # Check that the exposure time keyword is defined candark=no;n_extdark=0 if(l_fl_dark) { s_empty=""; print(l_dark) | scan(s_empty); l_dark=s_empty gimverify(l_dark); l_dark=gimverify.outname if(gimverify.status==0) { candark = yes keypar(l_dark//"[0]","GDARK",silent+) if(!keypar.found) { printlog("ERROR - GIREDUCE: Dark frame: "//l_dark//" was not \ created with gdark task.",l_logfile, yes) badflag=yes } imgets(l_dark//"[0]","NSCIEXT", >>& "dev$null") n_extdark = int(imgets.value) for(ii=1; ii<=n_extdark; ii+=1) { keypar(l_dark//"["//l_sci_ext//","//ii//"]","TRIMSEC",silent+) imgets(l_dark//"["//l_sci_ext//","//ii//"]","DETSEC", >& "dev$null") ccd_dark[ii] = imgets.value imgets(l_dark//"["//l_sci_ext//","//ii//"]","ccdsum", >& "dev$null") ccd_dark[ii] = ccd_dark[ii]//imgets.value } imgets(l_dark//"[0]",l_key_exptime, >>& "dev$null") if(imgets.value=="0") { printlog("ERROR - GIREDUCE: Dark image "//l_dark//" contains no exposure time", l_logfile, yes) badflag=yes } else { l_darkexptime=real(imgets.value) } } else { printlog("ERROR - GIREDUCE: Can not find dark frame: "//l_dark,l_logfile, yes) badflag=yes } } ################################################################################ # If any of the calibration files is missing, exit the program. if(badflag) goto crash ################################################################################ # Redefine the flags: if(!canbias) l_fl_bias=no if(!candark) l_fl_dark=no if(!canflat[1] && !canflat[2] && !canflat[3] && !canflat[4]) l_fl_flat=no ################################################################################ # Do the actual job: # The input images name are in an array input[ii] # The output images name are in an array: output[ii] # If there is a MDF file, place it in the output MEF file. # Do the overscan of biassec with colbias. # Trim the image to datasec # Do the bias subtraction # Do the dark subtraction # Do the flat field correction # Do the job with imexpr and the expressions above. for(ii=1; ii<=nimages; ii+=1) { # Make the output file - place header, only if the image does not exist if(n_do[ii]) imcopy(input[ii]//"[0]",output[ii]//".fits",verbose-) else input[ii] = output[ii] # The number of extensions will be needed in the MDF call below and later imgets(input[ii]//"[0]","NEXTEND", >& "dev$null") n_ext[ii] = int(imgets.value) # Add the MDF to the output file if appropriate hasmdf[ii]=no imgets (input[ii]//"[0]","OBSMODE", >>& "dev$null") if (imgets.value != "IMAGE") { for (jj = 1; jj<=n_ext[ii]; jj+=1) { keypar(input[ii]//".fits["//jj//"]","EXTNAME",sile=yes,>& "dev$null") if (keypar.value == "MDF") { hasmdf[ii]=yes # fxinsert(input[ii]//".fits",output[ii]//".fits[0]",groups=jj,verbose-, >& "dev$null") n_ext[ii]=n_ext[ii]-1 # subtract 1 from number of image ext } } } # Loop over science extensions for (jj = 1; jj<=n_ext[ii]; jj+=1) { imgets(input[ii]//"["//l_sci_ext//","//jj//"]","DETSEC", >& "dev$null") ccd[ii,jj] = imgets.value imgets(input[ii]//"["//l_sci_ext//","//jj//"]","CCDSUM", >& "dev$null") ccd[ii,jj] = ccd[ii,jj]//imgets.value } } ################################################################################ # Populate the overvar[6,200] matrix, rms of overscan region for(ii=1; ii<=nimages; ii+=1) { for(jj=1; jj<=n_ext[ii]; jj+=1) overvar[jj,ii]=0. } ################################################################################ # If fl_vardq=yes, the output images have already been gprepared outside # of gireduce and they do not contain var and dq planes, fl_vardq=no if(l_fl_vardq) { for(ii=1; ii<=nimages; ii+=1) { for(jj=1; jj<=n_ext[ii]; jj+=1) { st1 = input[ii]//"["//l_var_ext//","//jj//"]" st2 = input[ii]//"["//l_dq_ext//","//jj//"]" if(!imaccess(st1) || !imaccess(st2)) { l_fl_vardq=no jj=n_ext[ii]+1 ii=nimages+1 printlog("WARNING - GIREDUCE: Some of the input gprepared images do \ not have VAR or DQ planes.", l_logfile,l_verbose) printlog(" Setting fl_vardq=no.",l_logfile,l_verbose) printlog("",l_logfile,l_verbose) } } } } ################################################################################ # Overscan correction & Trimming - they are both together which maked this a long # bit of the code, but I think this will run faster than checking the general flags # when calling each image. # First check that the overscan has not been applied yet # Second check that the image has not been trimmed yet. # If they have not been overscan corrected nor trimmed, get the # biassec and datasec keywords if(l_fl_over || l_fl_trim) { canover=l_fl_over ; cantrim=l_fl_trim for(ii=1; ii<=nimages; ii+=1) { if(canover) { keypar (input[ii]//"[0]","OVERSCAN",silent+) canover = !keypar.found if(!canover) printlog("WARNING - GIREDUCE: "//input[ii]//" already OVERSCAN corrected",l_logfile,l_verbose) } if(cantrim) { keypar (input[ii]//"[0]","TRIMMED",silent+) cantrim = !keypar.found if(!cantrim) printlog("WARNING - GIREDUCE: "//input[ii]//" already TRIMMED",l_logfile,l_verbose) } if (canover) { printlog("Image_extension Overscan_mean Overscan_rms Function Order",l_logfile,l_verbose) # Update BIASSEC header to avoid overscan contamination if (n_ext[ii]==3) { hedit(input[ii]//"["//l_sci_ext//",1]",l_key_biassec,"[1:11,*]",update+,verify-,show-) hedit(input[ii]//"["//l_sci_ext//",2]",l_key_biassec,"[1:11,*]",update+,verify-,show-) imgets(input[ii]//"["//l_sci_ext//",3]","i_naxis1", >& "dev$null") nx=int(imgets.value) hedit(input[ii]//"["//l_sci_ext//",3]",l_key_biassec,"["//(nx-10)//":"//nx//",*]",update+,verify-,show-) } # Loop through the extensions for(jj=1; jj<=n_ext[ii]; jj+=1) { keypar(input[ii]//"["//l_sci_ext//","//jj//"]",l_key_biassec,silent+) found1 = keypar.found ; l_biassec = keypar.value if(cantrim) { keypar(input[ii]//"["//l_sci_ext//","//jj//"]",l_key_datasec,silent+) found2 = keypar.found ; l_datasec = keypar.value } else { l_datasec="[]" } print("yes") | colbias(input[ii]//"["//l_sci_ext//","//jj//"]", output[ii]//"["//l_sci_ext//","//jj//",append]", bias=l_biassec,trim=l_datasec,median=l_median, interactive=l_fl_inter,function=l_function, order=l_order,low_reject=l_low_reject, high_reject=l_high_reject,niterate=l_niterate, logfile=tmplog,graphics="stdgraph",cursor="", >& "dev$null") match ("RMS",tmplog,stop-) | scan(srms,srms,srms,srms) overvar[jj,ii] = real(srms) match ("order",tmplog,stop-) | scan(sorder,sorder,sorder,sorder) l_order=int(sorder) delete(tmplog,verify-, >& "dev$null") imstat(input[ii]//"["//l_sci_ext//","//jj//"]"//l_biassec,fields="mean", format=no) | scan(dummy) printf("%-25s %6.2f %6.3f %-10s %2d\n",input[ii]//"["//l_sci_ext//","//jj//"]",dummy,overvar[jj,ii],l_function,l_order) | scan(l_struct) printlog(l_struct,l_logfile,l_verbose) overvar[jj,ii] = overvar[jj,ii]**2 # Variance gemhedit(output[ii]//"["//l_sci_ext//","//jj//"]","OVERSCAN",dummy,"Overscan mean value") gemhedit(output[ii]//"["//l_sci_ext//","//jj//"]","OVERRMS",srms, "Overscan RMS value from colbias") if(cantrim) { hedit(output[ii]//"["//l_sci_ext//","//jj//"]","TRIMSEC", l_datasec,upd+,ver-,del-,add+,show-) } } # end of loop through extensions if(cantrim) { gsetsec(output[ii],key_datsec=l_key_datasec) gemhedit(output[ii]//"[0]","TRIMMED","yes","Overscan section trimmed") printlog("Image "//output[ii]//" trimmed",l_logfile,l_verbose) } gemhedit(output[ii]//"[0]","OVERSCAN","yes","Overscan correction") } else if (!canover && cantrim) { for(jj=1; jj<=n_ext[ii]; jj+=1) { keypar(input[ii]//"["//l_sci_ext//","//jj//"]",l_key_datasec,silent+) found2 = keypar.found ; l_datasec = keypar.value imcopy(input[ii]//"["//l_sci_ext//","//jj//"]"//l_datasec, output[ii]//"["//l_sci_ext//","//jj//",append]",ver-) hedit(output[ii]//"["//l_sci_ext//","//jj//"]","TRIMSEC", l_datasec,upd+,ver-,del-,add+,show-) } } gsetsec(output[ii],key_datsec=l_key_datasec) gemhedit(output[ii]//"[0]","TRIMMED","yes","Overscan section trimmed") printlog("Image "//output[ii]//" trimmed",l_logfile,l_verbose) } # end of loop over images } # Is there anything else to do? if (!l_fl_bias && !l_fl_dark && !l_fl_flat) { goto finish } if(!l_fl_trim && !l_fl_over) { for(ii=1; ii<=nimages; ii+=1) { for(jj=1; jj<=n_ext[ii]; jj+=1) { imcopy(input[ii]//"["//l_sci_ext//","//jj//"]", output[ii]//"["//l_sci_ext//","//jj//",append]",ver+) } } } # There are 6 cases: # 1) bias subtraction only # 2) dark correction only # 3) flat field correction only # 4) bias + dark correction # 5) bias + flat field # 6) dark + flat field # Merge all 6 cases if(l_fl_bias || l_fl_flat || l_fl_dark) { canflat[1]=l_fl_flat ; canbias=l_fl_bias ; candark=l_fl_dark printlog("",l_logfile,l_verbose) if(candark) printlog("If dark subtraction: scale = exptime_input/\ exptime_dark",l_logfile,l_verbose) printf("%-18s %-22s %-22s %5s %-22s %4s %4s\n","Output image","Bias","Dark","scale","Flat","Read","Gain") | scan(l_struct) printlog(l_struct,l_logfile,l_verbose) for(ii=1; ii<=nimages; ii+=1) { canflat[1]=l_fl_flat ; canbias=l_fl_bias ; candark=l_fl_dark # BIAS: Check the image and bias frame have the same number of extensions # and detector section if(canbias) { if(n_ext[ii]!=n_extbias) { printlog("WARNING - GIREDUCE: No bias subtraction. \ "//input[ii]//" and bias image ("//l_bias//") have different \ number of extensions",l_logfile,yes) canbias=no } for(jj=1; jj<=n_ext[ii]; jj+=1) { if(ccd[ii,jj]!=ccd_bias[jj]) { printlog("WARNING - GIREDUCE: Mismatch detector section or binning of bias",l_logfile,yes) printlog("("//l_bias//") and image ("//input[ii]//")",l_logfile,yes) canbias=no } } # Check the input image has not been bias corrected before keypar(output[ii]//"[0]","BIASIM",silent+) if(keypar.found) { printlog("WARNING - GIREDUCE: Image "//input[ii]//" \ already BIAS subtracted",l_logfile,yes) canbias=no } } # end of canbias checks # DARK: Check the image and dark frame have the same number of extensions l_darkscale=0. if(candark) { if(n_ext[ii]!=n_extdark) { printlog("WARNING - GIREDUCE: No dark count correction. \ "//input[ii]//" and dark image ("//l_dark//") have different \ number of extensions",l_logfile,yes) candark=no } for(jj=1; jj<=n_ext[ii]; jj+=1) { if(ccd[ii,jj]!=ccd_dark[jj]) { printlog("WARNING - GIREDUCE: Mismatch detector section or binning of dark",l_logfile,yes) printlog("("//l_dark//") and image ("//input[ii]//").",l_logfile,yes) candark=no } } # Check the input image has not been dark corrected before keypar(output[ii]//"[0]","DARKIM",silent+) if(keypar.found) { printlog("WARNING - GIREDUCE: Image "//input[ii]//" already dark \ corrected.",l_logfile,yes) candark=no } # Get the exposure time from the image keypar(output[ii]//"[0]",l_key_exptime,silent+) if(!keypar.found) { printlog("WARNING - GIREDUCE: ("//output[ii]//") does \ not have an exposure time in header from keyword "//l_key_exptime, l_logfile,yes) candark=no } else l_darkscale=real(keypar.value)/l_darkexptime } # end of candark checks # FLAT checks if(canflat[1]) { # Check the input image has not been flat corrected before keypar(output[ii]//"[0]","FLATIM",silent+) if(keypar.found) { printlog("WARNING - GIREDUCE: Image "//input[ii]//" already \ FLAT FIELD corrected.",l_logfile,yes) canflat[1]=no } # Get the filter of the input image keypar(input[ii]//"[0]",l_key_filter,silent+) thisfilter = keypar.value if(!keypar.found) { printlog("WARNING - GIREDUCE: Filter of "//input[ii]//" \ is unknown.", l_logfile,yes) canflat[1]=no } # Pick the flat that has the same filter as the input image fnum=0 for(kk=1; kk<=4; kk+=1) { if(thisfilter==flatfilter[kk]) { fnum=kk; kk=5 } } if(fnum==0) { printlog("WARNING - GIREDUCE: Can not apply flat field correction \ to image "//input[ii]//" - no match in filters", l_logfile,l_verbose) canflat[1]=no } # Check that FLATS and inputs have the same number of extension if(n_ext[ii]!=n_extflat[fnum]) { printlog("WARNING - GIREDUCE: Can not apply flat field correction \ to image "//input[ii]//" - images have different sizes", l_logfile,yes) canflat[1]=no } for(jj=1; jj<=n_ext[ii]; jj+=1) { if(ccd[ii,jj]!=ccd_flat[fnum,jj]) { printlog("WARNING - GIREDUCE: Mismatch detector section or binning of flat",l_logfile,yes) printlog("("//flatfilter[fnum]//") and image ("//input[ii]//")",l_logfile,yes) canflat[1]=no } } # Check obsmode of image and flat keypar(input[ii]//"[0]","OBSMODE",silent+) if(obsmodeflat[fnum]!=keypar.value) { printlog("WARNING - GIREDUCE: OBSMODE parameter of "//\ input[ii]//" and flat field do not match",l_logfile,l_verbose) canflat[1]=no } } # end of canflat checks if(!canbias && !canflat[1] && !candark) goto nextimage1 if(canbias) s_bias=l_bias else s_bias="INDEF" if(candark) s_dark=l_dark else s_dark="INDEF" if(canflat[1]) s_flat=l_flat[fnum] else s_flat="INDEF" printf("%-18s %-22s %-22s %5.2f %-22s\n",output[ii],s_bias,s_dark, l_darkscale,s_flat) | scan(l_struct) printlog(l_struct,l_logfile,l_verbose) for(jj=1; jj<=n_ext[ii]; jj+=1) { # Get the gain from the header #keypar(output[ii]//"["//l_sci_ext//","//jj//"]",l_key_gain,silent+) #if(!keypar.found) { # Debug loop - not a real fix, IJ # printlog("WARNING - GIREDUCE: Cannot find "//l_key_gain//" in ext "//str(jj),l_logfile,l_verbose) # n_gain[jj]=1 #} else # n_gain[jj]=real(keypar.value) imgets(output[ii]//"["//l_sci_ext//","//jj//"]",l_key_gain, >>& "dev$null") if (imgets.value == "0") { printlog("WARNING - GIREDUCE: Cannot find "//l_key_gain//" in ext "//str(jj),l_logfile,l_verbose) n_gain[jj]=1 } else { n_gain[jj]=real(imgets.value) } # Get the trim section in the input/output image extension #keypar(output[ii]//"["//l_sci_ext//","//jj//"]","TRIMSEC",silent+) #tsec_img = keypar.value; fl_tsec_img=keypar.found imgets(output[ii]//"["//l_sci_ext//","//jj//"]","TRIMSEC", >>& "dev$null") if (imgets.value == "0") { fl_tsec_img=no } else { fl_tsec_img=yes tsec_img=imgets.value } ll_sciexp="a" if (canbias) { # Check the BIAS and inputs have the same dimensions in each extension # biasimg = l_bias//"["//l_sci_ext//","//jj//"]" subbias = "" if(fl_tsec_img && fl_tsec_bias[jj] && tsec_bias[jj]!=tsec_img) { printlog("WARNING - GIREDUCE: "//input[ii]//" \ and bias image ("//l_bias//") have been trimmed \ to different sections in the CCD.", l_logfile,l_verbose) goto nextextension1 } else if (!fl_tsec_img && fl_tsec_bias[jj]) { printlog("ERROR - GIREDUCE: Can not bias subtract.",l_logfile,yes) printlog(" Trim flag is not set but bias image \ ("//l_bias//") is trimmed to "//tsec_bias, l_logfile,yes) goto nextextension1 } else if (fl_tsec_img && !fl_tsec_bias[jj]) { printlog("WARNING - GIREDUCE: "//input[ii]//" \ is trimmed but bias image ("//l_bias//") is not.", l_logfile,l_verbose) printlog(" Using a section \ of the bias image",l_logfile,l_verbose) subbias=tsec_img } ll_bias=l_bias//"["//l_sci_ext//","//jj//"]"//subbias ll_sciexp=ll_sciexp//"-b" } else { ll_bias="INDEF" } if (candark) { # Check that DARK and inputs have the same dimensions in each extension # darkimg = l_dark//"["//l_sci_ext//","//jj//"]" subdark = "" if(fl_tsec_img && fl_tsec_dark[jj] && tsec_dark[jj]!=tsec_img) { printlog("WARNING - GIREDUCE: "//input[ii]//" \ and dark image ("//l_dark//") have been trimmed \ to different sections in the CCD.", l_logfile,l_verbose) goto nextextension4 } else if (!fl_tsec_img && fl_tsec_dark[jj]) { printlog("ERROR - GIREDUCE: Can not dark count correct.",l_logfile,yes) printlog(" Trim flag is not set but dark image \ ("//l_dark//") is trimmed to "//tsec_dark[jj], l_logfile, l_verbose) goto nextextension4 } else if (fl_tsec_img && !fl_tsec_dark[jj]) { printlog("WARNING - GIREDUCE: "//input[ii]//" \ is trimmed but dark image ("//l_dark//") is not.", l_logfile,l_verbose) printlog(" Using a section \ of the dark image",l_logfile,l_verbose) subdark=tsec_img } ll_dark=l_dark//"["//l_sci_ext//","//jj//"]"//subdark ll_sciexp=ll_sciexp//"-c*"//str(l_darkscale) } else { ll_dark="INDEF" } if(canflat[1]) { # Check the FLATS and inputs have the same dimensions in each extension #flatimg = l_flat[fnum]//"["//l_sci_ext//","//jj//"]" subflat = "" if(fl_tsec_img && fl_tsec_flat[jj,fnum] && tsec_flat[jj,fnum]!=tsec_img) { printlog("WARNING - GIREDUCE: "//input[ii]//" \ and flat image ("//l_flat[fnum]//") have been trimmed \ to different sections in the CCD.", l_logfile,l_verbose) goto nextextension1 } else if (!fl_tsec_img && fl_tsec_flat[jj,fnum]) { printlog("ERROR - GIREDUCE: Can not apply flat field correction.", l_logfile,yes) printlog(" Trim flag is set to no but flat image \ ("//l_flat[fnum]//") is trimmed to "//tsec_flat[jj,fnum], l_logfile,l_verbose) goto nextextension1 } else if (fl_tsec_img && !fl_tsec_flat[jj,fnum]) { printlog("WARNING - GIREDUCE: "//input[ii]//" \ is trimmed but flat image ("//l_flat[fnum]//") is not", l_logfile,l_verbose) printlog(" Using subsection of the\ flat image",l_logfile,l_verbose) subflat = tsec_img } ll_flat=l_flat[fnum]//"["//l_sci_ext//","//jj//"]"//subflat ll_sciexp="("//ll_sciexp//")/d*e" } else { ll_flat="INDEF" } # The bias image is: # l_bias//"["//l_sci_ext//","//jj//"]"//tsec_img # The flat field image is: # l_flat[fnum]//"["//l_sci_ext//","//jj//"]"//tsec_img imexpr(ll_sciexp,output[ii]//"["//l_sci_ext//","//jj//",overwrite]", output[ii]//"["//l_sci_ext//","//jj//"]", ll_bias,ll_dark,ll_flat,n_gain[jj], dims="auto",intype="real",outtype="real",refim="auto",bwidth=0, btype="nearest",bpixval=0.,rangecheck+,verbose-,exprdb="none") if(l_fl_vardq) { ll_dqexp="im1" ll_varexp="a+c" # c is the variance from the overscan subtraction, zero if not overscan subtracting # need to fix output to contain VAR,DQ at this point dqlist = output[ii]//"["//l_dq_ext//","//jj//"]" ll_ndq=1 if(ll_bias=="INDEF") { ll_varbias="INDEF" ; ll_dqbias="INDEF" } else { ll_varbias=l_bias//"["//l_var_ext//","//jj//"]"//subbias ll_varexp=ll_varexp//"+b" dqlist = dqlist//","//l_bias//"["//l_dq_ext//","//jj//"]"//subbias ll_ndq+=1 } if(ll_dark=="INDEF") { ll_vardark="INDEF" ; ll_dqdark="INDEF" } else { ll_vardark=l_dark//"["//l_var_ext//","//jj//"]"//subdark ll_varexp=ll_varexp//"+b+h*"//str(l_darkscale)//"*"//str(l_darkscale) dqlist = dqlist//","//l_dark//"["//l_dq_ext//","//jj//"]"//subdark ll_ndq+=1 } if(ll_flat=="INDEF") { ll_varflat="INDEF" ; ll_dqflat="INDEF" ll_varexp=ll_varexp } else { ll_varflat=l_flat[fnum]//"["//l_var_ext//","//jj//"]"//subflat ll_varexp="(("//ll_varexp//")*g*g+d*d*e*e)/(f*f)" dqlist = dqlist//","//l_flat[fnum]//"["//l_dq_ext//","//jj//"]"//subflat ll_ndq+=1 } for(kk=2;kk<=ll_ndq;kk+=1) ll_dqexp=ll_dqexp//" | im"//str(kk) print(ll_varexp) print(ll_dqexp) imexpr(ll_varexp, output[ii]//"["//l_var_ext//","//jj//",overwrite]", output[ii]//"["//l_var_ext//","//jj//"]", ll_varbias, overvar[jj,ii], output[ii]//"["//l_sci_ext//","//jj//"]", ll_varflat, ll_flat, n_gain[jj], ll_dark, dims="auto",intype="real",outtype="real",refim="auto",bwidth=0, btype="nearest",bpixval=0.,rangecheck+,verbose-,exprdb="none") addmasks(dqlist,output[ii]//"["//l_dq_ext//","//jj//",overwrite]", ll_dqexp,flags="") } if(canflat[1]) gemhedit(output[ii]//"["//l_sci_ext//","//jj//"]","GAIN",1.,"Gain") nextextension1: } if(canbias) gemhedit(output[ii]//"[0]","BIASIM",l_bias,"Bias image used by GIREDUCE") if(candark) gemhedit(output[ii]//"[0]","DARKIM",l_dark,"Dark image used by GIREDUCE") if(canflat[1]) gemhedit(output[ii]//"[0]","FLATIM",l_flat[fnum],"Flat field used by GIREDUCE") nextimage1: } # End for(ii=1; i<=nimages ... } finish: # Copy MDF for (ii=1; ii<=nimages; ii+=1) { if (hasmdf[ii]) { tcopy(input[ii]//".fits[MDF]",tmpmdf//".fits",verbose-) fxinsert(tmpmdf//".fits",output[ii]//".fits["//n_ext[ii]//"]","1",verbose-, >& "dev$null") delete(tmpmdf//".fits",verify-) } } ################################################################################ # Update the headers: GIREDUCE, GEM-TLM date | scan(l_struct) for(ii=1; ii<=nimages; ii+=1) { gemhedit(output[ii]//"[0]","GIREDUCE",l_struct,"Time stamp for GIREDUCE") gemhedit(output[ii]//"[0]","GEM-TLM",l_struct,"Last modification with GEMINI") } goto clean # Exit with error crash: status=1 for(ii=1; ii<=nimages; ii+=1) { if(n_do[ii]) imdelete(output[ii],ver-, >>& "dev$null") } # Clean up and exit clean: scanfile="" delete(tmpfile//","//tmpinlist//","//tmpoutlist//","//tmplog,ver-, >& "dev$null") printlog("",l_logfile,l_verbose) if (status == 0) { printlog("GIREDUCE exit status: good.",l_logfile,l_verbose) printlog("----------------------------------------------------------------------------",l_logfile,l_verbose) } else { printlog("GIREDUCE exit status: error.",l_logfile,l_verbose) printlog("----------------------------------------------------------------------------",l_logfile,l_verbose) } end