; NAME: ; darcurrent_analysis.pro (v1.5) ; ; 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. FUNCTION month2days,month,year ;function written by Sungsoo Kim, taken from darkplot2.pro ; leap month calculation valid until 2099 leap = 0 if((year mod 4) eq 0) then leap = 1 CASE month OF 1: days = 0 2: days = 31 3: days = 59+leap 4: days = 90+leap 5: days = 120+leap 6: days = 151+leap 7: days = 181+leap 8: days = 212+leap 9: days = 243+leap 10: days = 273+leap 11: days = 304+leap 12: days = 334+leap ELSE: begin print,'Problem with month2days: ',month stop end ENDCASE return, days END ;----------------------------------------------------------------------------------- pro darkcurrent_analysis, testdir, outdir=outdir, infile=infile, $ jumbofile=jumbofile, rd_time = rd_time ;set directory path delimiter if (!version.os_family eq 'unix') then begin delim='/' endif else begin delim='\' endelse 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 ; 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 ; 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 ; split the infoline into its three component parts infoparts = strsplit(infoline, ' ', /extract) ; the image name is the first part of the infoline imname = infoparts[0] ; 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[1]) reads_in_sample = fix(infoparts[2]) ; read just the header of the specified FITS image fits_read, imname, 0, hdr, /header_only ; extract information from the header readmode = sxpar(hdr, 'readmode') numreads = sxpar(hdr, 'naxis3') detname = sxpar(hdr, 'detname') dettemp = sxpar(hdr, 'dettemp') 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(2, samples)) theTime = ptr_new(fltarr(samples)) theMode = ptr_new(fltarr(samples)) readf, dflun, *combo (*thetime)[*] = (*combo)[0, *] (*themode)[*] = (*combo)[1, *] 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 ; Fit a line to the data. ; Consider only points where signal is less than 30000 so that only the ; non-saturated part is involved in the fit goodPoints = where((*themode) lt 30000.0) linePars = POLY_FIT(((*thetime)[goodPoints])[1:*], $ ((*themode)[goodPoints])[1:*], 1) ; set axis ranges for plots--the actual data range takes up the central ; 65% of the plot height yrng = (max(*themode) - min(*themode))/0.65 ylo = min(*themode) - yrng*0.2 yhi = max(*themode) + yrng*0.15 xrng = max(*thetime) xlo = 0.0 - xrng*0.05 xhi = xrng*1.05 ; 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 for outcount=0,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=[14,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 ; Plot the result plot, *thetime, *themode, $ TITLE='Mode of Signal vs. Time: '+filename, $ XTITLE='time (seconds)', $ YTITLE='Signal (ADU)', $ 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) 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 + $ ' ADU/s', 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, 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 plots, 0, _A + _B * 0, color=black, thick=lthick plots, max((*thetime)[goodPoints]), $ _A + _B * max((*thetime)[goodPoints]), /CONTINUE, $ color = black, thick = lthick ; 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]) ; calculate time in units of hours from 1/1/2001 ; leap month calculation valid until 2099 jumbo_hour = ((year-2001)*365+(year-2001)/4)*24. + $ (month2days(month,year)+day-1.)*24. + hour + minu/60. + sec/3600. print,filename, dettemp, jumbo_hour, _B printf, djunit, filename, dettemp, jumbo_hour, _B ; read blank line separator datasep = '' readf, dflun, datasep ; free up dynamic memory ptr_free, combo, theTime, theMode endwhile free_lun, dflun free_lun, djunit end