; NAME:
;			variance_analysis.pro (v1.0)
;
; PURPOSE:
;			To read variance vs. time values from an input data file
;			and produce plots from the data
;
; CALLING SEQUENCE:
;			variance_analysis, testdir[, infile[, outdir]]
;
; INPUTS:
;			testdir		the name of a directory in which the input data file resides
;
; KEYWORD PARAMETERS:
;			infile		the name of the input data file.  defaults to 'variancevals.txt' if not set
;			outdir		the name of the directory where output plots are written; defaults to the same directory as testdir
;
; EXAMPLE:
;			variance_analysis, '.', outdir = '.\ReadNoise_Results'
;
; REFERENCE:
;
;
; MODIFICATION HISTORY:
;       Written by:  Ernie Morse, IDTL, April 27, 2003
;

Pro variance_analysis, testdir, infile=infile, outdir=outdir

	;begin by setting defaults for keyword values not provided by the calling context
	if (not keyword_set(outdir)) then outdir = testdir
	if (not keyword_set(infile)) then infile = 'variancevals.txt'

	;set directory path delimiter based on operating system in use
	if (!version.os_family eq 'unix') then begin
		pathdelim = '/'
	endif else begin
		pathdelim = '\'
	endelse

	;check to see if input and outdirectory names end in path delimiter--if not, add it to the end
	if (strmid(testdir, strlen(testdir)-1, 1) ne pathdelim) then testdir = testdir + pathdelim
	if (strmid(outdir, strlen(outdir)-1, 1) ne pathdelim) then outdir = outdir + pathdelim

	; open the input data file and extract the header values
	openr, inlun, testdir+infile, /get_lun

	; declare some variables to hold results from file
	detector = (ref_ver = (ref_type = (gdate = (region = ''))))

	readf, inlun, detector, format='(A)'
	readf, inlun, ref_ver, format='(A)'
	readf, inlun, ref_type, format='(A)'
	readf, inlun, gdate, format='(A)'
	readf, inlun, region, format='(A)'


	; declare variables to hold the values from the file
	dettemp = 0
	numvals = 0

	while (not eof(inlun)) do begin

		; read the detector temperature and number of reads
		readf, inlun, dettemp, numvals

		; create 2D array to hold read times in column 0 and variances in column 1
		time_var = fltarr(2, numvals)

		; read the time-variance values
		readf, inlun, time_var

		; there is a blank line separating blocks of data.  It doesn't have to be read in explicitly, but sometimes
		; there is an extra blank line at the end of the file that does need to be read.  Reading it here prevents an EOF error.
		datasep = ''
		readf, inlun, datasep

		; now that the data is present, make PostScript and JPEG plots

		; make file names by extracting name of infile up to the extension and adding new
		; extensions of the appropriate type
		jpgfile = outdir + (strsplit(infile, '.', /extract))[0] + '_' + strtrim(string(dettemp),2) + 'K.jpg'
		psfile = outdir + (strsplit(infile, '.', /extract))[0] + '_' + strtrim(string(dettemp),2) + 'K.ps'

		 ; set axis ranges for plots--the actual data range takes up the central 65% of the plot height
		 yrng = (max(time_var[1,1:*]) - min(time_var[1,1:*]))/0.65
		 ylo = min(time_var[1,1:*]) - yrng*0.2
		 yhi = max(time_var[1,1:*]) + yrng*0.15
		 xrng = max(time_var[0,*])
		 xlo = 0.0 - xrng*0.05
		 xhi = xrng*1.05

		 ; set plot colors
		 black = 0
         white = 255

		for outcount = 0,1 do begin
			;first set up for landscape Postscript plots
			if outcount eq 0 then begin
				set_plot, 'ps'
				!p.font = -1
				device, filename=psfile, /landscape
				fsub = psfile                            ;filename to be written at bottom of plot

				!Y.MINOR=10
				!y.margin=[6,2]
				!x.margin=[14,4]

				; set line and text characteristics for publication quality
				!P.CHARSIZE=1.2
				!P.CHARTHICK=4.
				!X.THICK=8.
				!Y.THICK=8.
				!P.THICK=8.
				!P.MULTI=0
				!P.TICKLEN=0.03
				symsize=1.5
				filenamecharsize=1.1

			endif else begin
				set_plot, 'z'
				device, set_resolution=[8000,6000]
				device, set_font='Courier'
				fsub = jpgfile

				!Y.MINOR=10
				!y.margin=[6,4]
				!x.margin=[14,4]

				; set line and text characteristics for publication quality
				!P.CHARSIZE=10
				!P.CHARTHICK=15.
				!X.THICK=20.
				!Y.THICK=20.
				!P.THICK=20.
				symsize=10.
				filenamecharsize=8.
			endelse

			; Define a special plotting symbol to plot big, fat points
			_radius = !P.charsize ;; 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                   ;; Position angle
			_delta_theta = 2 * !PI / 10 ;; Angle step size
			for i=0, 9 do begin
				_X[i] = _radius * cos(_theta) ;; Project radius onto x axis
				_Y[i] = _radius * sin(_theta) ;; Project radius onto y axis
				_theta = _theta + _delta_theta ;; Increment angle
			endfor
			usersym, _X, _Y, /FILL    ;; Define the symbol

           ; Plot the result
           plot, time_var[0,1:*], time_var[1,1:*], $
                TITLE='Variance vs. Time at ' + strtrim(string(dettemp),2) + 'K', $
                XTITLE='Time (seconds)', $
                YTITLE='Variance', $
                psym=8, $
                XRANGE=[xlo,xhi], $
                XSTYLE=1, $
                YRANGE=[ylo,yhi], $
                YSTYLE=1, $
                BACKGROUND=white, $
                COLOR=black


    		; write plot annotations
			xyouts, xlo+xrng*0.025, ylo + (yrng*0.12), 'Analysis Region: ' + region, color=black
			xyouts, xlo+xrng*0.025, ylo + (yrng*0.09), 'Refsub Version: ' + ref_ver, color=black
			xyouts, xlo+xrng*0.025, ylo + (yrng*0.06), 'Refsub Type: ' + ref_type, color=black
			xyouts, xhi-xrng*0.025, ylo + (yrng*0.03), 'Date: ' + gdate, color=black, alignment=1
			xyouts, xlo+xrng*0.025, yhi - (yrng*0.045), 'Detector:  ' + detector, color=black


		   ; write the file name outside the margin of the plot
		   xyouts, 1.0, 40, fsub, charsize=filenamecharsize, color=black, /device

           ;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

	endwhile

	free_lun, inlun


End