; NAME: ; darcurrent_analysis.pro (v1.7) ; ; PURPOSE: ; This procedure reads an input data file produced by ; darkcurrent_reduce.pro and uses the data therein to produce a series ; of dark current slope plots associated with images generated in a ; dark current experiment. ; CALLING SEQUENCE: ; darkcurrent_analysis,testdir, [outdir=outdir, [infile=infile ; [,jumbofile [,rd_time]]]] ; ; INPUTS: ; testdir the full path name of the directory containing an input data ; file ; ; KEYWORD PARAMETERS: ; outdir the full path name of a directory in which the output plots ; will be written. If unspecified, plots will be placed in ; the directory indicated by the testdir input variable ; ; infile the name of the input data file. If not specified, ; the default filename generated by darkcurrent_reduce.pro ; (slopevals.txt) will be used. ; ; jumbofile the name of the output data file containing values used by ; the dark jumbo plotting routine. The default value is ; 'darkjumbo.txt' ; ; rd_time keyword contains the value of the time to read the detector ; one time. The value can be read from the FITS header for ; newer Rockwell data, but is missing in Raytheon data and ; and older Rockwell data. The value is only necessary when ; processing old Fowler non-glow images. Otherwise, this can ; be ignored. ; EXAMPLE ; ; dkdir = '\\rabbit\Raytheon\SB304-008-5.0mu\cold1\dktest.15Mar03\' ; darkcurrent_analysis, dkdir ; ; ; REFERENCE: ; ; ; MODIFICATION HISTORY: ; Created by: ; Ernie Morse, IDTL, March 19, 2003 (v1.0) ; ; Modified by: ; Ernie Morse, IDTL, April 3, 2003 (v1.1) ; Fixed various bugs ; Replaced "convert PS" method of JPG writing with "z-buffer" ; method ; ; Ernie Morse, IDTL, April 10, 2003 (v1.2) ; Changed name to darkcurrent_analysis (instead of ; darkcurrent_analyze) to be consistent with naming of other ; analysis procedures. ; ; Ernie Morse, IDTL, April 30, 2003 (v1.3) ; Added code to deal with the generation of slope plots for the ; non-glow Fowler images ; Added jumbofile keyword to allow for writing the jumbo plot data ; to an arbitrary file name ; ; Ernie Morse, IDTL, October 17, 2003 (v1.4) ; Made changes to deal with new format of input files created by ; the latest version of darkcurrent_reduce.pro (v2.10) ; ; Ernie Morse, IDTL, November 5, 2003 (v1.5) ; Modified code so that the procedure stops if the rd_time value ; is not specified and not available in the header when processing ; Fowler data. Also cleaned up code to ensure that it meets IDTL ; standards. ; Ernie Morse, IDTL, July 15, 2004 (v1.6) ; Altered code to read in another column (raw signal) from ; input file. ; Altered code so that reads with raw signal values railed at the ; low end of the A/D range are not included in the fit. ; Ernie Morsie, IDTL, August 9, 2004 (v1.7) ; Removed call to month2day to translate hours since 1/1/2001 and replaced ; with call to julday to do the same thing. ; Don Figer, RIDL, May 1, 2007 ; Changed processing of infoline to allow spaces in the path name ; Don Figer, RIDL, July 21, 2007 ; Added saturation parameter ; Don Figer, RIDL, July 30, 2007 ; Added capability to accept a conversion gain (e) and pass it to rnplot pro darkcurrent_analysis, testdir, outdir=outdir, infile=infile,jumbofile=jumbofile, rd_time = rd_time, jumbovar=jumbovar,saturation=saturation,e=e ;set directory path delimiter delim=path_sep() if (not keyword_set(outdir)) then outdir = testdir if (not keyword_set(infile)) then infile = 'slopevals.txt' if (not keyword_set(jumbofile)) then jumbofile = 'darkjumbo.txt' ; add delimiter at end of directory names if necessary if (strlen(testdir)-1 ne rstrpos(testdir,delim)) then testdir=testdir+delim if (strlen(outdir)-1 ne rstrpos(outdir,delim)) then outdir=outdir+delim ; define saturation limits if (not keyword_set(saturation)) then saturation=60000.0 ; Set up for nice quality PostScript plots !p.charsize = 1.4 !x.charsize = 1.4 !y.charsize = 1.4 !p.charthick = 4.0 !x.thick = 8.0 !y.thick = 8.0 !p.thick = 8.0 !x.minor = 5 !y.minor = 5 !y.margin=[6,2] !x.margin=[14,4] ; open input data file openr, dflun, testdir+infile, /get_lun ; read first 3 lines, which contain data common to all of the images in the ; data set infoline = (ref_ver = (ref_type = (region = ''))) readf, dflun, ref_ver readf, dflun, ref_type readf, dflun, region FIRST_FILE = 1 ; loop through the rest of the file and produce plots for each set of image data read while (not eof(dflun)) do begin wait, 0.1 ; read the info line, which should contain the name of the image, ; followed by the number of read groups and the number of reads per ; group. readf, dflun, infoline ; the image name is the first part of the infoline imname = strmid(infoline,0,strpos(infoline,'.fits '))+'.fits' ; split the infoline into its three component parts infoparts = strsplit(strmid(infoline,strlen(imname),99), ' ', /extract) ; set up error handler to deal with the situation where there are ; not three elements in the infoparts array. This will happen with ; older input files where sampling information wasn't included. The ; error handler will fix the situation. ; catch, err_stat ; if (err_stat ne 0) then begin ; infoparts = strarr(3) ; infoparts[1] = '0' ; infoparts[2] = '1' ; endif ; now assign the other two parts of the infoline into variables samples = fix(infoparts[0]) reads_in_sample = fix(infoparts[1]) ; read just the header of the specified FITS image fits_read_datacube, imname, dummy, hdr, first=0,last=10 ; extract information from the header readmode = sxpar(hdr, 'readmode') numreads = sxpar(hdr, 'naxis3') detname = sxpar(hdr, 'detname') dettemp = sxpar(hdr, 'MOLYTABL') exp_time = sxpar(hdr, 'exp_time') numsamrd = sxpar(hdr, 'numsamrd') ; create arrays to hold sample times and mode values ; use dynamic memory that can be freed at the end of each loop ; iteration. Otherwise, memory use could get out of hand. ; If information on sampling isn't included in the input file, ; just assume that there is one mode value for each read of the input ; image and that the file contains no read times (only mode values). if (samples ne 0) then begin combo = ptr_new(fltarr(4, samples)) theTime = ptr_new(fltarr(samples)) theMode = ptr_new(fltarr(samples)) theMedian = ptr_new(fltarr(samples)) theCorrection = ptr_new(fltarr(samples)) rawsig = ptr_new(fltarr(samples)) readf, dflun, *combo (*thetime)[*] = (*combo)[0, *] (*themode)[*] = (*combo)[1, *] (*rawsig)[*] = (*combo)[2, *] (*themedian)[*] = (*combo)[3, *] ; (*thecorrection)[*] = (*combo)[3, *] endif else begin combo = ptr_new() theMode = ptr_new(fltarr(numreads)) readf, dflun, *theMode if (readmode eq 'Fowler') then begin if (not keyword_set(rd_time)) then begin rd_time = sxpar(hdr, 'frametime') endif if (rd_time eq 0) then begin expstr = 'Frame time not present in header and not ' + $ 'specified in keyword rd_time. Determine the ' + $ 'frame time for these data and re-try with ' + $ 'rd_time set equal to that value. ' stop, expstr endif theTime = ptr_new(findgen(numreads) * rd_time) (*theTime)[numsamrd:*] = (*theTime)[numsamrd:*] + exp_time endif else begin theTime = ptr_new(findgen(numreads) / (numreads - 1) * exp_time) endelse endelse ; IDTL FITS headers are non-standard, so the following function calls ; are necessary to extract time and date info from the header date and ; time keywords ;timeStamp = '' ;dateStamp = '' ;gettime, hdr, timeStamp ;getdate, hdr, dateStamp timestamp=sxpar(hdr,'TIME') ; define the values that will be analyzed xvals=*theTime yvals=*theMode if (detname eq 'H4RG-007-10') then yvals=*rawsig if (detname eq 'HRG-10-007') then yvals=*rawsig if (detname eq 'H4RG-10-007') then yvals=*rawsig if (detname eq 'H1RG-022-SIPIN') then yvals=*rawsig if (detname eq 'H1RG-0220SIPIN') then yvals=*rawsig if keyword_set(e) then yvals=yvals*e ; set yvals so that they start at zero yvals=yvals-yvals(0) ; Fit a line to the data. ; Consider only points where raw signal is less than saturation so that only the ; non-saturated part is involved in the fit. Skip the first point which has no information. goodPoints = where((*rawsig) lt saturation and (*rawsig) gt 0,ngoodpoints) if (ngoodpoints gt 2) then begin linePars = POLY_FIT(((xvals)[goodPoints])[1:*],((yvals)[goodPoints])[1:*], 1) endif else begin linePars = [0, 0] endelse ; set axis ranges for plots--the actual data range takes up the central ; 65% of the plot height yrng = (max(yvals) - min(yvals))/0.65 if (yrng eq 0.0) then begin yrng = 100.0 ylo = -50.0 yhi = 50.0 endif else begin ylo = min(yvals) - yrng*0.2 yhi = max(yvals) + yrng*0.15 endelse xrng = max(xvals) - min(xvals) xlo = min(xvals) - xrng*0.05 xhi = max(xvals) + xrng*0.05 xrng = xrng * 1.1 ; define filenames for PS and JPG plots ; start by extracting just the file name and discarding the path filename = strmid(imname, strpos(imname, delim, /reverse_search)+1) psfile = outdir+strmid(filename, 0, strpos(filename, '.fits'))+'.ps' jpgfile = outdir+strmid(filename, 0, strpos(filename, '.fits'))+'.jpg' ; define the readmode string to be used in plot annotation case readmode of 'Fowler' : begin rdmode_str = 'Fowler-' + strtrim(string(numsamrd),2) end 'UTR' : begin rdmode_str = 'SUTR-' + strtrim(string(numreads),2) + ',' + $ strtrim(string(numsamrd),2) end endcase ; set plot colors black = 0 white = 255 ; make the actual plots twice, saving one in PS format and one in ; JPG format ; note that the loop count has been changed to only write JPG files. change it back to 0,1 to ; get both ps and JPG for outcount=1,1 do begin if (outcount eq 0) then begin ; Open the postscript device set_plot, 'ps' device, filename = psfile, /LANDSCAPE !p.font = -1 ; use vector fonts for PS !P.MULTI = 0 ; Set up for one plot per page fsub = psfile ; set text to print at bottom of plot ; set plot margin & tick marks !Y.MINOR=0 !y.margin=[6,2] !x.margin=[14,4] ; set some size values asize = 1.12 ; charsize of annotations athick = 4.0 ; charthick of annotations lthick = 8.0 ; line thickness endif else begin ; Open the z-buffer device set_plot, 'z' ; set resolution and font device, set_resolution=[8000,6000] device, set_font='Courier' !P.MULTI = 0 ; Set up for one plot per page fsub = jpgfile ; set text to print at bottom of plot ; set plot margin & tick marks !Y.MINOR=0 !y.margin=[6,4] !x.margin=[16,4] ; set some size values asize = 8.0 ; charsize of annotations athick = 15.0 ; charthick of annotations lthick = 20.0 ; line thickness endelse ; Define a special plotting symbol to plot big, fat points _radius = asize ; Offset radius is about size of a character _X = fltarr(10) ; Place for x verticies, offset from point _Y = fltarr(10) ; Place for y verticies, offset from point _theta = 0 _delta_theta = 2 * !PI / 10 for i=0, 9 do begin _X[i] = _radius * cos(_theta) _Y[i] = _radius * sin(_theta) _theta = _theta + _delta_theta endfor ; set the symbol created above as IDL symbol #8 usersym, _X, _Y, /FILL ytitle='Signal (ADU)' slopeunits=' ADU/s' if keyword_set(e) then begin ytitle='Signal (e-)' slopeunits=' e-/s' endif ; Plot the result plot, xvals, yvals, $ TITLE='Median of Signal vs. Time: '+filename, $ XTITLE='time (seconds)', $ YTITLE=ytitle, $ psym=8, $ XRANGE=[xlo,xhi], $ XSTYLE=1, $ YRANGE=[ylo,yhi], $ YSTYLE=1, $ BACKGROUND=white, $ charsize=asize*1.25, charthick=athick, xthick=lthick, $ ythick=lthick, $ COLOR=black ; Add labels & fitted line _A = linePars[0] _B = linePars[1] slopestr = strtrim(string(linePars[1], format='(F7.3)'),2)+slopeunits xyouts, xlo + xrng * 0.025, ylo + (yrng * 0.12), $ 'Analysis Region: ' + region, charsize=asize, $ charthick=athick, color=black xyouts, xlo+xrng*0.025, ylo + (yrng*0.09), 'Refsub Version: ' + $ ref_ver, charsize=asize, charthick=athick, color=black xyouts, xlo+xrng*0.025, ylo + (yrng*0.06), 'Refsub Type: ' + $ ref_type, charsize=asize, charthick=athick, color=black xyouts, xlo+xrng*0.025, ylo + (yrng*0.03), 'Slope = '+slopestr, $ charsize=asize, charthick=athick, color=black ; xyouts, xhi-xrng*0.025, ylo + (yrng*0.03), 'Date: ' + dateStamp + $ ; ', Time: '+timeStamp + ' UT', charsize=asize, charthick=athick,$ ; color=black, alignment=1 xyouts, xhi-xrng*0.025, ylo + (yrng*0.03), 'Time: '+timeStamp + ' UT', charsize=asize, charthick=athick,$ color=black, alignment=1 xyouts, xlo+xrng*0.025, yhi - (yrng*0.135), 'Reads Per Point: ' + $ strtrim(string(reads_in_sample), 2), charsize = asize, $ charthick = athick, color = black xyouts, xlo+xrng*0.025, yhi - (yrng*0.105), 'Read Mode: ' + $ rdmode_str, charsize = asize, charthick = athick, color = black xyouts, xlo+xrng*0.025, yhi - (yrng*0.075), $ 'Detector Temperature: ' + $ strtrim(string(dettemp, format='(F7.3)'),2) + ' K', $ charsize = asize, charthick = athick, color = black xyouts, xlo+xrng*0.025, yhi - (yrng*0.045), 'Detector: ' + $ detname, charsize = asize, charthick = athick, color = black ; write the file name outside the margin of the plot xyouts, 1.0, 40, fsub, charsize = asize, charthick = athick, $ color = black, /device ; plot line with calculated slope if (goodpoints[0] ne -1) then begin plots, (xvals)[goodpoints[0]], _A + _B * (xvals)[goodpoints[0]], $ color=black, thick=lthick plots, max((xvals)[goodPoints]), $ _A + _B * max((xvals)[goodPoints]), /CONTINUE, $ color = black, thick = lthick endif ; For the JPEG plot, the plot image must be read in from memory and ; written out to a file if (outcount eq 1) then begin jpgimg = tvrd() write_jpeg, jpgfile, $ congrid(jpgimg, 1600, 1200, /center, /interp), quality=100 endif ; finish up by closing device device, /close endfor ; Write time stamp and fitted line slope to a ; results file. This is hardwired for now. if (FIRST_FILE) then begin FIRST_FILE = 0 openw, results_file_unit, outdir + 'RESULTS.TXT', /get_lun openw,djunit,outdir+jumbofile, /get_lun printf, djunit, detname printf, djunit, 'SUTR-' + strtrim(string(numreads),2) + ',' + $ strtrim(string(numsamrd),2) printf, djunit, ref_ver printf, djunit, ref_type printf, djunit, region endif printf, results_file_unit, timeStamp, ' ', _B ; Save data for darkjumbo.pro ; date_time = dateStamp + ' ' +timeStamp ; temp1 = str_sep(dateStamp,'/') ; temp2 = str_sep(timeStamp,':') ; year = fix(temp1[2]) ; month = fix(temp1[0]) ; day = fix(temp1[1]) ; hour = fix(temp2[0]) ; minu = fix(temp2[1]) ; sec = fix(temp2[2]) 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) hour=strmid(timestamp,12,2) minu=strmid(timestamp,15,2) sec=strmid(timestamp,18,2) ; calculate time in units of hours from 1/1/2001 if (strlen(timestamp) gt 10) then begin jumbo_hour = (julday(month, day, year, hour, minu, sec) - julday(1, 1, 2001, 0, 0, 0)) * 24.0d endif else begin jumbo_hour=0 endelse print,filename, dettemp, jumbo_hour, _B printf, djunit, filename, dettemp, jumbo_hour, _B,format='(A,F10.3,F15.3,F20.6)' ; read blank line separator datasep = '' readf, dflun, datasep ; free up dynamic memory ptr_free, combo, theTime, theMode, theMedian endwhile free_lun, dflun free_lun, djunit end