; NAME: ; readnoise_reduce (v2.12) ; ; PURPOSE: ; This program reads input image filenames from a file, calculates readnoise ; on pairs of images, and writes the results to a file. ; ; CALLING SEQUENCE: ; IDL> readnoise_reduce, testdir, [outdir, infile, refcor, plottitle, region, ; minexp, novar, varonly, tnint] ; ; INPUTS: ; testdir - string containing name of directory holding list of input ; files and the files themselves ; ; KEYWORD PARAMETERS: ; outdir - string containing name of directory to which output should be ; written (default value is testdir + Readnoise_Results) ; infile - input file name. File must be in directory specified by ; testdir. Default name is 'filelist.txt' ; noref - set to disable reference pixel correction ; ma_refcor - set for maximal averaging reference pixel correction ; plottitle - string containing title for plot. This string is appended to ; 'Readnoise test for ' by the rnplot procedure. Default value ; is the detector name read from the FITS header. ; ; region - a 4-vector of integers specifying the region of a 2D image ; plane over which readnoise is calculated. For an array ; section A[x1:x2, y1:y2], region should be specified as ; [x1, x2, y1, y2]. Default value is assigned based on detector ; type. ; minexp - set if read noise in continuous reads is desired (default is ; total noise in 1000-seconds) ; novar - set if variance vs. time should not be calculated (can be ; used with default or minexp set) ; varonly - set if only variance vs. time should be calculated (value of ; minexp ignored if this keyword is set) ; tnint - exposure time to use for total noise option--default is ; calculated and is the maximum possible value. ; rowinc - row increment for region ; colinc - column increment for region ; ; NOTE: The default behavior is to apply spatial averaging reference pixel ; correction. If noref is set, the status of ma_refcor is ignored. ; ; EXAMPLE ; ; IDL> rndir = '\\rabbit\rockwell\SB-304-003\warm1\readnoise.15Nov02' ; IDL> resultsdir = rndir + '\Results' ; IDL> readnoise_reduce, rndir, outdir = resultsdir, /refcor, ; region = [4, 2051, 0, 2047] ; ; REFERENCE: ; ; MODIFICATION HISTORY: ; Written by: ; Mike Telewicz, IDTL ; Modified by: ; Ernie Morse, IDTL, November 15, 2002 ; Added header ; Added call to procedure to create PS and GIF plots ; Updated reference pixel correction to call refsub by B.J. Rauscher ; Eliminated sigma clipping and good pixel mask from input (no longer ; used by readnoise.pro) ; Ernie Morse, IDTL, February 4, 2003 ; Changed to reading input FITS images one plane at a time in order ; to allow for handling of very large files. Commented out call to ; mkdiffimg and instead process CDS and Fowler-sampled images within ; the body of this procedure. ; Ernie Morse, IDTL, March 27, 2003 (v2.0) ; Changed to allow processing on long dark ramps ; Ernie Morse & Don Figer, IDTL, April 10, 2003 (v2.1) ; Changed to use readcol to read all images at once instead of ; reading one-by-one in loop ; Changed name to readnoise_reduce.pro ; Added flush command to force buffered output to file immediately ; Ernie Morse, IDTL, April 14, 2003 (v2.2) ; Changed calls to refsub--must now pass sca name instead of sca type ; Ernie Morse, IDTL, April 16, 2003 (v2.3) ; Fixed problem of railed pixels by marking them as NaN before feeding ; to readnoise.pro ; Bernie Rauscher, IDTL, April 23, 2003 (v2.4) ; Added minexp keyword. ; Ernie Morse, IDTL, April 21 - 27, 2003 (v2.5) ; Added a whole bunch of stuff to support calculating variance vs. ; time ; New keywords novar and varonly added ; Ernie Morse, IDTL, May 1, 2003 (v2.6) ; Corrected some bugs introduced in the change to v2.5 ; Ernie Morse, IDTL, May 5, 2003 (kept version at v2.6) ; Added if-then statement to prevent railedpix from being applied as ; an array index if no railed pixels were present (this came up when ; attempting to do a shorted resistor read noise test.) ; Ernie Morse, IDTL, May 6, 2003 (v2.7) ; Eliminated sa_refcor keyword--spatial average reference pixel ; correction is now default ; Added noref keyword--can be set to explicitly disable reference ; pixel correction ; Made changes in code to reflect new default case (spatial averaging). ; Changed default output directory to testdir + 'ReadNoise_Results' ; Ernie Morse, IDTL, February 5, 2004 (v2.8) ; Edited to ensure compliance with IDTL coding standards ; Sito Balleza, IDTL, March 20, 2004 ; Added SIPIN in detector names. ; Ernie Morse, IDTL, May 18, 2004 (v2.9) ; Added total noise interval keyword (tnint) to define the exposure time ; to use for the total noise test. ; Ernie Morse, IDTL, June 21, 2004 (v2.10) ; Changed code so that tnint is calculated if not provided. ; Ernie Morse, IDTL, June 29, 2004 (v2.11) ; Added printing of final newline character to variance output file-- ; variance_analysis crashes without it. ; Ernie Morse, IDTL, August 9, 2004 (v2.12) ; Added rowinc and colinc keywords to allow for non-contiguous ; regions ; Don Figer, RIDL, July 30, 2007 ; Added capability to accept a conversion gain (e) and pass it to rnplot pro readnoise_reduce, testdir, outdir = outdir, infile = infile, $ ma_refcor = ma_refcor, noref = noref, plottitle = plottitle, $ region = region, minexp = minexp, novar = novar, varonly = varonly, $ tnint = tnint, colinc = colinc, rowinc = rowinc, e=e ; specify file path delimiter based pathdelim = path_sep() ; check to see if input directory name ends in path delimiter--if not, ; add it to the end if (strmid(testdir, strlen(testdir) - 1, 1) ne pathdelim) then testdir = testdir + pathdelim ; set outdir to testdir if not set by user if (keyword_set(outdir) eq 0) then outdir = testdir + 'ReadNoise_Results' + pathdelim if (not keyword_set(infile)) then infile = 'filelist.txt' if (not keyword_set(rowinc)) then rowinc = 1 if (not keyword_set(colinc)) then colinc = 1 ; check to see if outdirectory name ends in path delimiter--if not, \ ; add it to the end if (strmid(outdir, strlen(outdir) - 1, 1) ne pathdelim) then begin outdir = outdir + pathdelim endif ; open files for output if (not keyword_set(varonly)) then begin if (keyword_set(minexp)) then begin datfile = outdir + 'rdnoisevals.txt' plotfile = outdir + 'rdnoiseplot' endif else begin datfile = outdir + 'noisevals.txt' plotfile = outdir + 'noiseplot' endelse endif if (not(keyword_set(novar))) then begin varfile = outdir + 'variancevals.txt' varplot = outdir + 'varianceplot' endif ; if reference pixel correction is to be done, query the refsub function for ; the version number, otherwise, just set a flag indicating it was not done if (not keyword_set(noref)) then begin ; create non-empty string needed as input to refsub version parameter ; version # returned into ref_ver by refsub ref_ver = ' ' refsub, 0, sca = ' ', version = ref_ver if (not keyword_set(ma_refcor)) then begin ref_type = 'Spatial averaging' endif else begin ref_type = 'Maximal averaging' endelse endif else begin ref_ver = 'Not Used' ref_type = 'N/A' endelse ; open input file and read all filenames into an array readcol, testdir + infile, imgfiles, format = '(A)' ; loop through the images in the file list, performing read noise on difference images ; paircount keeps track of # of pairs of "like" images processed paircount = 0 ; readsum holds a running sum of the calculated readnoise-- it is divided by ; paircount later to determine average read noise over all the pairs of ; "like" images processed readsum = 0 ; ; for imgcount = 0, n_elements(imgfiles) - 2 do begin filename1 = imgfiles[imgcount] filename2 = imgfiles[imgcount + 1] ; Determine the readmode of both our images (read header only) and get ; some other needed values from the image headers. ; first get the info for the current image fits_read_datacube, testdir + filename1, image1, header1, first=1, last=10 ;readmode1 = sxpar(header1, "readmode") ;sample1 = sxpar(header1, "numsamrd") ;numutr1 = sxpar(header1, 'numutr') ;temp1 = sxpar(header1, 'dettemp') naxis1 = sxpar(header1, 'naxis1') naxis2 = sxpar(header1, 'naxis2') naxis3 = sxpar(header1, 'naxis3') exptime1 = sxpar(header1, 'exp_time') readmode1='UTR' numutr1 = fix(sxpar(header1, 'NGROUPS')) sample1 = fix(sxpar(header1, 'NREADS')) if (sample1 gt numutr1) then begin sample1 = fix(sxpar(header1, 'NGROUPS')) numutr1 = fix(sxpar(header1, 'NREADS')) endif ndrops = fix(sxpar(header1, 'NDROPS')) temp1 = float(sxpar(header1, 'MOLYTABL')) exptime1 = 4096.*4096./32.*10.e-6*(ndrops*numutr1+sample1*numutr1/2.) ; get info for the next image in the list fits_read_datacube, testdir + filename2, image2, header2, first=1, last=10 ;readmode2 = sxpar(header2, "readmode") readmode2='UTR' ;sample2 = sxpar(header2, "numsamrd") ;numutr2 = sxpar(header2, 'numutr') ;temp2 = sxpar(header2, 'dettemp') numutr2 = fix(sxpar(header2, 'NGROUPS')) sample2 = fix(sxpar(header2, 'NREADS')) if (sample2 gt numutr2) then begin sample2 = fix(sxpar(header2, 'NGROUPS')) numutr2 = fix(sxpar(header2, 'NREADS')) endif temp2 = float(sxpar(header2, 'MOLYTABL')) ; if temp1 equals 0, that means the dettemp header key was not present ; and the experiment was performed without temperature control. ; Assume a warm test (295 K) in this situation if (temp1 eq 0) then begin temp1 = 295 endif if (temp2 eq 0) then begin temp2 = 295 endif ; if this is the first image overall, get some information common to all ; the images from the header if (imgcount eq 0) then begin ;read detector name from header scaName = sxpar(header1, 'detname') ;parse sca type from detector name scaType = (strsplit(scaName, '-', /extract))[0] ; set default region based on scatype if (not keyword_set(region)) then begin case scaType OF 'HRG' : region = [132, 4091, 4, 4091] 'H4RG' : region = [132, 4091, 4, 4091] 'H1RG' : region = [4, 1019, 4, 1019] 'H2RG' : region = [4, 2043, 4, 2043] 'SIPIN': region = [4, 2043, 4, 2043] 'SB304': region = [4, 2051, 0, 2047] ELSE : begin print, 'Detector type ' + scatype + ' not supported.' $ + ' Exiting procedure.' return end endcase endif ; get parameters for plot routine if (keyword_set(plottitle) eq 0) then plottitle = scaName ; sxpar doesn't get the date value from IDTL FITS header, so the ; next two lines find the header element that contains the date and ; extract the text after the ' = ' character ; plotdate = (header1[where(strpos(header1, 'DATE =') ne -1)])[0] ; plotdate = getwrd(plotdate, 2) timestamp=sxpar(header1,'TIME') timestamp=' '+timestamp MONTHS = ['JAN','FEB','MAR','APR','MAY','JUN','JUL','AUG', 'SEP','OCT','NOV','DEC'] MONTH = long(where(strupcase(strmid(timestamp,5,3)) eq MONTHS))+1 day=strmid(timestamp,9,2) year=strmid(timestamp,21,4) plotdate=string(month)+'/'+string(day)+'/'+string(year) ; create string representation of region, to be used later as a ; plotting label region_str = [strtrim(string(region[0:1]), 2), $ strtrim(string(colinc), 2), $ strtrim(string(region[2:3]), 2), $ strtrim(string(rowinc), 2)] ; write some information to the header of the output file if (not keyword_set(varonly)) then begin openw, lun2, datfile, /get_lun printf, lun2, plottitle printf, lun2, ref_ver printf, lun2, ref_type printf, lun2, plotdate printf, lun2, region_str, $ format = '("[", A, ":", A, ":", A, ", "' + ', A, ":", A, ":", A, "]")' endif ; open the file for the variance data, unless the keyword was set ; to disable it ; print header information if (not keyword_set(novar)) then begin openw, varlun, varfile, /get_lun printf, varlun, plottitle printf, varlun, ref_ver printf, varlun, ref_type printf, varlun, plotdate printf, varlun, region_str, $ format = '("[", A, ":", A, ":", A, ", "' + ', A, ":", A, ":", A, "]")' endif endif ; now begin processing all the data ; check each image against the next one in the queue--if they have the ; same readmode, # of reads and detector temperature, ; find read noise on their difference image. Otherwise, go on to the next image. if (readmode2 eq readmode1 AND sample1 eq sample2 AND numutr1 eq numutr2 AND round(temp1) eq round(temp2)) then begin ; if we are just starting with images of a certain readmode, number ; of reads and temperature, initialize some values ; readsum is an array that holds a running sum of read noise values ; numreps determines the number of Fowler pairs that should be ; processed from the pair of images ; backstartval is the # of the read representing the 1st read of the ; 2nd group for Fowler sampling ; pairs is an array that holds the Fowler samples that will be ; processed on a pair of images (i.e, [1, 2, 4, 8] to do Fowler-1, ; Fowler-2, Fowler-4 and Fowler-8) if (paircount eq 0) then begin ; for Fowler exposures if (readmode1 ne 'UTR') then begin numreps = 1 readsum = fltarr(1) variancesum = fltarr(2, numutr1) ; find the exposure time for each read, save in first column ; of variance array variancesum[0, *] = (findgen(numutr1) / (numutr1 - 1)) * exptime1 backstartval = sample1 pairs = [sample1] utrflag = 0 ; for UTR exposures, simulated Fowler sampling is used endif else begin variancesum = fltarr(2, numutr1) ; find the exposure time for each read, save in first column ; of variance array variancesum[0, *] = (findgen(numutr1) / (numutr1-1)) * exptime1 ; max number of fowler samples is one less than highest power of 2 ; less than the number of reads, with a maximum of 32 numpairs = (2 ^ (floor(alog10(numutr1) * 3.32193) - 1)) < 32 if (not keyword_set(minexp)) then begin ; targetread will be the read that is the starting point from ; which the 2nd group of Fowlers will be read for ; tnint-second exposure time processing ; if tnint is set, the target read is the read which ; occurred closest to that time if (keyword_set(tnint)) then begin mindiff = min(abs(variancesum[0, *] - tnint), targetread) endif else begin ; if tnint is not set, the target read is the one numpairs ; prior to the final read targetread = numutr1 - numpairs tnint = variancesum[0, targetread] endelse endif ; set up the pairs array pairs = indgen(round(alog10(numpairs) * 3.32193) + 1) pairs = 2 ^ pairs ;the pairs array should now contain something like ; [1, 2, 4, 8, 16, 32], which is the Fowlers that we want to ; consider numreps = n_elements(pairs) ; we will process all the Fowlers on a pair of images before ; moving on to the next pair.\ ; readsum will hold a running total of the read noise values ; for each Fowler-N under consideration. readsum = fltarr(numreps) endelse endif ; if we're doing variance, set all exposure times to negative to ; indicate that all reads need to be processed if (not keyword_set(novar)) then variancesum[0, *] = -1 * ((findgen(numutr1) / (numutr1-1)) * exptime1) ; define arrays to hold first and last groups of Fowler data for ; like images if (not keyword_set(varonly)) then begin ; 1st Fowler group of 1st image group1_1 = lonarr(naxis1, naxis2) ; 2nd Fowler group of 1st image group1_2 = lonarr(naxis1, naxis2) ; 1st Fowler group of 2nd image group2_1 = lonarr(naxis1, naxis2) ; 2nd Fowler group of 2nd image group2_2 = lonarr(naxis1, naxis2) endif ; stepval determines how many pixels should be read in order to ; read one frame stepval = long(naxis1) * naxis2 - 1 ; this variable will keep track of how many frames have been read ; into each of the group arrays. this is used to prevent ; re-reading the same data multiple times frames_read = 0 ; loop through the various Fowler modes for a pair of images for fowlercount = 0, numreps-1 do begin ; if minexp is set, we need to set a different value for ; backstartval every time through the fowlercount loop ; otherwise, we always use the same value (the read represented ; by targetread) if (keyword_set(minexp)) then begin backstartval = pairs[fowlercount] endif else begin backstartval = targetread endelse ; read planes for 1st group and store their running total in ; arrays ; set beginning point to beginning of frame after the last read ; frame startwd = (stepval + 1) * frames_read ; if minexp is set, the first Fowler group should be read only ; for the Fowler-1 case. ; Otherwise, we can just move add group 2 to group 1 to get the ; new group 1. if (keyword_set(minexp)) then begin if (fowlercount eq 0) then begin group1end = 1 endif else begin group1end = 0 endelse endif else begin group1end = pairs[fowlercount] - frames_read endelse for i = 1, group1end do begin endwd = startwd + stepval ; read from the current read from both images print,'Reading 1st read of ',testdir + filename1 fits_read_datacube, testdir + filename1, image1, header1, first = startwd, last = endwd print,'Reading 1st read of ',testdir + filename2 fits_read_datacube, testdir + filename2, image2, header2, first = startwd, last = endwd ; determine the absolute read # readnum = (i + frames_read) - 1 print, 'Processing pair ' + strtrim(string(paircount), 2) + ', read # ' + strtrim(string(readnum), 2) + ' in 1st Fowler group' ; find the railed pixels in the first read of the 1st image if (frames_read eq 0) then begin railedpix = where(image1 gt 65535 OR image1 lt 0) firstread1 = image1 firstread2 = image2 endif else begin if (not(keyword_set(novar))) then begin if (readnum ne 0) then begin cds1 = image1 - firstread1 cds2 = image2 - firstread2 cds1 = reform(cds1, naxis1, naxis2, /overwrite) cds2 = reform(cds2, naxis1, naxis2, /overwrite) ; apply reference pixel correction if requested ; by the user if (not keyword_set(noref)) then begin if (not keyword_set(ma_refcor)) then begin refsub, cds1, sca = scaName refsub, cds2, sca = scaName endif else begin refsub, cds1, sca = scaName, /simple refsub, cds2, sca = scaName, /simple endelse endif if (railedpix[0] ne -1) then begin cds1[railedpix] = !values.f_nan cds2[railedpix] = !values.f_nan endif varianceval = readnoise(cds1, cds2, region, colinc, rowinc) variancesum[1, readnum] = variancesum[1, readnum] + varianceval ; make the time positive to indicate that this ; read has been processed variancesum[0, readnum] = -1 * variancesum[0, readnum] endif endif endelse if (not keyword_set(varonly)) then begin ; add the current reads to the group arrays group1_1 = group1_1 + image1 group2_1 = group2_1 + image2 endif startwd = endwd + 1 endfor ;find start read position for 2nd group startwd = backstartval * (stepval + 1) + (stepval + 1) * frames_read ; if minexp is set, we should read the whole N reads into group ; 2--if not, we can just read the ones that aren't already there ; this will be accomplished by setting frames_read to 0 if ; minexp is set. this is done later in the code ; read planes for 2nd group and store their running total in arrays for j = 1, pairs[fowlercount] - frames_read do begin endwd = startwd + stepval print,'Reading 2nd read of ',testdir + filename1 fits_read_datacube, testdir + filename1, image1, header1, first = startwd, last = endwd print,'Reading 2nd read of ',testdir + filename2 fits_read_datacube, testdir + filename2, image2, header1, first = startwd, last = endwd ; determine the absolute read # readnum = (j + frames_read + backstartval) - 1 print, 'Processing pair ' + strtrim(string(paircount), 2) + ', read # ' + strtrim(string(readnum), 2) wait,0.2 if (not(keyword_set(novar))) then begin cds1 = float(image1) - float(firstread1) cds2 = float(image2) - float(firstread2) cds1 = reform(cds1, naxis1, naxis2, /overwrite) cds2 = reform(cds2, naxis1, naxis2, /overwrite) ; apply reference pixel correction if requested by the user if (not keyword_set(noref)) then begin if (not keyword_set(ma_refcor)) then begin refsub, cds1, sca = scaName refsub, cds2, sca = scaName endif else begin refsub, cds1, sca = scaName, /simple refsub, cds2, sca = scaName, /simple endelse endif if (railedpix[0] ne -1) then begin cds1[railedpix] = !values.f_nan cds2[railedpix] = !values.f_nan endif varianceval = readnoise(cds1, cds2, region, colinc, rowinc) variancesum[1, readnum] = variancesum[1, readnum] + varianceval ; make the time positive to indicate that this read has been processed variancesum[0, readnum] = -1 * variancesum[0, readnum] endif if (not keyword_set(varonly)) then begin group1_2 = group1_2 + image1 group2_2 = group2_2 + image2 endif startwd = endwd + 1 endfor if (not keyword_set(varonly)) then begin ; construct Fowler-subtracted image for 1st image image1a = (float(group1_2) - float(group1_1)) / (pairs[fowlercount]) ; construct Fowler-subtracted image for 2nd image image2a = (float(group2_2) - float(group2_1)) / (pairs[fowlercount]) ; apply reference pixel correction if requested by the user if (not keyword_set(noref)) then begin if (not keyword_set(ma_refcor)) then begin refsub, image1a, sca = scaName refsub, image2a, sca = scaName endif else begin refsub, image1a, sca = scaName, /simple refsub, image2a, sca = scaName, /simple endelse endif ; calculate the read noise and print the result ; set the railed pixels to NaN before feeding to readnoise if (railedpix[0] ne -1) then begin image1a[railedpix] = !values.f_nan image2a[railedpix] = !values.f_nan endif readnoiseval = readnoise(image1a, image2a, region, colinc, rowinc) readsum[fowlercount] = readnoiseval + readsum[fowlercount] endif ; update the variable frames_read--it should now be equal to the ; Fowler-N value if minexp is not set ; if minexp is set, it should be 0 and the 2nd group should be ; added to the first if (keyword_set(minexp)) then begin frames_read = 0 group1_1 = group1_1 + group1_2 group2_1 = group2_1 + group2_2 group1_2[*] = 0.0 group2_2[*] = 0.0 endif else begin frames_read = pairs[fowlercount] endelse endfor ; do the variance calculation on all reads that didn't get done ; during the Fowler processing if (not keyword_set(novar)) then begin reads_notdone = where(variancesum[0, *] lt 0.0) for varcount = 0, n_elements(reads_notdone) - 1 do begin framenum = reads_notdone[varcount] print, 'Processing pair ' + strtrim(string(paircount), 2) $ + ', read # ' + strtrim(string(framenum), 2) $ + ' post-Fowler' fits_read_datacube, testdir + filename1, image1, header1, first = framenum * (stepval + 1), $ last = (framenum + 1)*(stepval + 1) - 1 fits_read_datacube, testdir + filename2, image2, header1, first = framenum * (stepval + 1), $ last = (framenum + 1) * (stepval + 1) - 1 cds1 = float(image1) - float(firstread1) cds2 = float(image2) - float(firstread2) cds1 = reform(cds1, naxis1, naxis2, /overwrite) cds2 = reform(cds2, naxis1, naxis2, /overwrite) ; apply reference pixel correction if requested by the user if (not keyword_set(noref)) then begin if (not keyword_set(ma_refcor)) then begin refsub, cds1, sca = scaName refsub, cds2, sca = scaName endif else begin refsub, cds1, sca = scaName, /simple refsub, cds2, sca = scaName, /simple endelse endif if (railedpix[0] ne -1) then begin cds1[railedpix] = !values.f_nan cds2[railedpix] = !values.f_nan endif varianceval = readnoise(cds1, cds2, region, colinc, rowinc) variancesum[1, framenum] = variancesum[1, framenum] $ + varianceval ; set time back to positive variancesum[0, framenum] = -1 * variancesum[0, framenum] endfor endif paircount = paircount + 1 ; if the two images aren't the same, then continue to the next after ; printing results endif else begin ; print results to file ; paircount is set to 0 after results are printed, so if it is ; now 0, there are no new results for this temperature if (paircount ne 0) then begin if (not keyword_set(varonly)) then begin for rescount = 0, n_elements(pairs)-1 do begin ; print the average read noise values for all pairs ; processed at this detector temperature; ; this is accomplished by dividing the summed readnoise ; by the number of pairs processed. print, pairs[rescount], round(temp1), $ (readsum[rescount] / paircount) printf, lun2, pairs[rescount], round(temp1), $ (readsum[rescount] / paircount) endfor ;force buffered output to file flush, lun2 endif if (not keyword_set(novar)) then begin ; the value we want to print is the square of the average ; noise from all the pairs processed at this detector ; temperature variancesum[1, *] = (variancesum[1, *] / paircount) ^ 2 printf, varlun, round(temp1), numutr1 printf, varlun, variancesum printf, varlun flush, varlun endif endif paircount = 0 endelse endfor ; if at end of file, print the final result to output file if (paircount ne 0) then begin if (not keyword_set(varonly)) then begin for rescount = 0, n_elements(pairs)-1 do begin ; print the average read noise values for all pairs processed at ; this detector temperature; ; this is accomplished by dividing the summed readnoise by the ; number of pairs processed. printf, lun2, pairs[rescount], round(temp1), $ (readsum[rescount] / paircount) print, pairs[rescount], round(temp1), $ (readsum[rescount] / paircount) endfor endif if (not keyword_set(novar)) then begin ; the value we want to print is the square of the average noise from ; all the pairs processed at this detector temperature variancesum[1, *] = (variancesum[1, *] / paircount) ^ 2 printf, varlun, round(temp1), numutr1 printf, varlun, variancesum printf, varlun endif endif if (not keyword_set(varonly)) then free_lun, lun2 if (not keyword_set(novar)) then free_lun, varlun if (not keyword_set(varonly)) then begin ; call plot routine to make plots if (keyword_set(minexp)) then begin ; If minexp keyword is set, make read noise plot rnplot, datfile, plotfile, minexp = minexp, e=e endif else begin ; Else make total noise plot (the default for RNPLOT.PRO) rnplot, datfile, plotfile, tnint = tnint, e=e endelse endif ; make the variance plots if (not keyword_set(novar)) then begin variance_analysis, outdir endif end