; NAME:
;   darkjumbo2.pro (v2.4)
;
; PURPOSE:
;   This program analyzes dark current images and produces a plot of
;   dark current per image as a function of time.
;
; CALLING SEQUENCE:
;   darkjumbo2, testdir [,outdir [, infile]]
;
; INPUTS:
;   testdir     directory containing input file that contains image file names,
;          detector temperature and dark current values for the data set.
;
; KEYWORD PARAMETERS:
;   outdir       directory where plots will be written.  The default is to use
;          the directory specified by testdir
;   infile       name of input file.  Default value is 'darkjumbo.txt'
;
;
; EXAMPLE
;       darkjumbo2,'/mydirectory/results'
;
;
; REFERENCE:
;
; MODIFICATION HISTORY:
;   Written by:
;     Sungsoo Kim & Don Figer, IDTL, August 16, 2002
;
;   Modified by:
;     Ernie Morse, IDTL, Nov. 26, 2002
;         Changed name of main function to darkjumbo2 to match file name
;
;     Don Figer, IDTL, April 23, 2003
;      Cleaned up code, modified plot characteristics
;
;     Ernie Morse, IDTL, April 23, 2003 (v2.3)
;      Changed determination of y-axis range and changed to printing
;      detector temperatures as integers
;
;     Ernie Morse, IDTL, November 5, 2003 (v2.4)
;      Updated header and cleaned up code to bring it into compliance with
;      IDTL standards.
;
;     Don Figer, IDTL, June 13, 2004 (v2.5)
;      Added option to add an extension (ext) on to the output file names
;      so that one can produce plots for initial dark data without overwriting
;      other plots.
;
;       Ernie Morse, IDTL, July 16, 2004 (v2.6)
;           Changed plot command parameters to force exponential notation, modified
;           plot margin and added !y.charsize to make y-axis text a little smaller
;
;		Don Figer, RIDL, July 30, 2007
;			Added capability to accept a conversion gain (e) and pass it to rnplot

FUNCTION getdatestring, hours
    ; leap month calculation valid until 2099
    days = ceil(hours / 24)
    leapdays = days / 1461
    year = (days - leapdays) / 365
    dayofyear = (days - (leapdays * 1461)) mod 365
    leap = 0
    if((year mod 4) eq 0) then leap = 1
    CASE 1 OF
       (dayofyear ge 0) AND (dayofyear le 31): month = '1'
       (dayofyear ge 32) AND (dayofyear le (59+leap)): month = '2'
       (dayofyear ge (60+leap)) AND (dayofyear le (90+leap)): month = '3'
       (dayofyear ge (91+leap)) AND (dayofyear le (120+leap)): month = '4'
       (dayofyear ge (121+leap)) AND (dayofyear le (151+leap)): month = '5'
       (dayofyear ge (152+leap)) AND (dayofyear le (181+leap)): month = '6'
       (dayofyear ge (182+leap)) AND (dayofyear le (212+leap)): month = '7'
       (dayofyear ge (213+leap)) AND (dayofyear le (243+leap)): month = '8'
       (dayofyear ge (244+leap)) AND (dayofyear le (273+leap)): month = '9'
       (dayofyear ge (278+leap)) AND (dayofyear le (304+leap)): month = '10'
       (dayofyear ge (305+leap)) AND (dayofyear le (334+leap)): month = '11'
       (dayofyear ge (335+leap)) AND (dayofyear le (365+leap)): month = '12'
    ELSE:  begin
       print,'Problem with month2days: ',month
       stop
    end
    ENDCASE

    CASE 1 OF
       (month eq '1'): day = strtrim(string(dayofyear),2)
       (month eq '2'): day = strtrim(string(dayofyear mod 31),2)
       (month eq '3'): day = strtrim(string(dayofyear mod (59+leap)),2)
       (month eq '4'): day = strtrim(string(dayofyear mod (90+leap)),2)
       (month eq '5'): day = strtrim(string(dayofyear mod (120+leap)),2)
       (month eq '6'): day = strtrim(string(dayofyear mod (151+leap)),2)
       (month eq '7'): day = strtrim(string(dayofyear mod (181+leap)),2)
       (month eq '8'): day = strtrim(string(dayofyear mod (212+leap)),2)
       (month eq '9'): day = strtrim(string(dayofyear mod (243+leap)),2)
       (month eq '10'): day = strtrim(string(dayofyear mod (273+leap)),2)
       (month eq '11'): day = strtrim(string(dayofyear mod (304+leap)),2)
       (month eq '12'): day = strtrim(string(dayofyear mod (334+leap)),2)
    ELSE:  begin
       print,'Problem with month2days: ',month
       stop
    end
    ENDCASE

    return, month+'/'+day+'/'+strtrim(string(year+2001),2)
END
;---------------------------------------------------------------------------------
PRO darkjumbo2,testdir, outdir=outdir, infile=infile,ext=ext,egain=egain

    if not(keyword_set(ext)) then ext=''

    ; set plot margins
    !y.margin=[6,2]
    !x.margin=[14,4]

    ; set directory path delimiter
    delim=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 delim) then begin
        testdir = testdir + delim
    endif

    ; set outdir equal to testdir if outdir is not set
    if (not keyword_set(outdir)) then outdir = testdir

    ; 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 delim) then begin
       outdir = outdir + delim
    endif

    ; set plotting devices
    devicesave=!d.name
    set_plot,'ps'

    ; define the default file that contains the list of filenames of images
    if (not keyword_set(infile)) then begin
       infile = testdir+'darkjumbo.txt'
    endif else begin
       infile = testdir+infile
    endelse

    rdmode = (region = (scaname = (refcor_version = (refcor_type = ''))))

    ; read image parameters from input file
    openr, lun1, infile, /get_lun
    readf, lun1, scaname, format='(A)'
    readf, lun1, rdmode, format='(A)'
    readf, lun1, refcor_version, format='(A)'
    readf, lun1, refcor_type, format='(A)'
    readf, lun1, region, format='(A)'
    free_lun, lun1

    ; read list of files from input file
    readcol, infile, files, dettemp, hour, slope, format='A,F,F,F', skipline=5
    units=' (ADU/s)'
    if keyword_set(egain) then begin
    	slope=slope*egain
 		units=' (e-/s)'
    endif
;   startdate = getdatestring(hour[0])
    caldat, (hour[0] / 24.0d) + julday(1, 1, 2001, 0, 0, 0), month, day, year
    dateparts = [month, day, year]
    dateparts = strtrim(string(dateparts), 2)
    startdate = dateparts[0] + '/' + dateparts[1] + '/' + dateparts[2]

    ; round temperatures to nearest integer and make an array of the
    ; unique temperature values
    dettemp = round(dettemp)
    sortedtemp = dettemp[sort(dettemp)]
    utemp = sortedtemp[uniq(sortedtemp)]

    ; define plot limits
    hour = hour - min(hour)
    xmax = max(hour)
    ymax = max(slope)
    xmin = min(hour)
    nonzero=where(slope gt 0,count)
    ymin=0
    if (count ne 0) then ymin = min(slope(nonzero))
    xrng = (xmax - xmin) / 0.8
    yrng = (ymax - ymin) / 0.7
    ylo = (0.0 - yrng*0.15) < (ymin - yrng*0.15)
    yhi = ymax + yrng*0.15
    xlo = xmin - xrng*0.1
    xhi = xmax + xrng*0.1

    ; define plot symbols
    symbols = [2, 4, 5, 6, 7, 1, 8, 8]
    nplots = ceil(n_elements(utemp)/float(n_elements(symbols)))

    ; draw curves
    for np=0,nplots-1 do begin

       ; define forground and background colors
       white = 255
       black = 0

       ; define plot file names
       psfile = outdir+'darkjumbo_'+scaname + '_' + strtrim(string(np),2) + ext+'.ps'
       jpgfile = outdir+'darkjumbo_'+scaname + '_' + strtrim(string(np),2) + ext+'.jpg'

       for plotcount=0,1 do begin
         if (plotcount eq 0) then begin

          ;set up for landscape Postscript plotting
          set_plot, 'ps'
          !p.font=-1
          device, filename= psfile, /landscape
          fsub = psfile     ;filename to be written at bottom of plot
          !x.minor = 10
          !y.minor = 10
          !y.margin=[6,2]
          !x.margin=[12,4]

          ; set line and text characteristics for publication quality
          !P.CHARSIZE=1.2
          !P.CHARTHICK=4.
          !y.charsize = 1.0
          !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 up to plot to memory buffer.  This plot will later be
          ; written from memory to a JPEG file.
          set_plot, 'z'
          device, set_resolution=[8000,6000]
          device, set_font='Courier'
          fsub = jpgfile   ; filename to be written at bottom of plot

          ;set plot margin & tick marks
          !Y.MINOR=10
          !x.minor=10
          !y.margin=[6,4]
          !x.margin=[15,4]

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

         endelse

         ; The series will be represented by big filled circles
         _radius = 1.0 ; 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

		; make plot axes
		if count ne 0 then begin
			ymin=10^float(round(alog10(min(slope(nonzero))*0.9)-0.5))
			ymax=10^float(round(alog10(max(slope(nonzero))*1.05)+0.5))
	         ; make plot axes
	         plot,[0],[0],/nodata,xrange=[xlo,xhi],xstyle=1,ystyle=1, $
	          xtitle='Time (hours)',ytitle='Dark Current'+units, $
	          title='Dark Current for '+scaname, $
	          color=black, background=white,charsize=!P.CHARSIZE*1.4, $
	;          yrange=[0,ymax*1.05], $
	;          yrange=[1e-4,(1 > (ymax*2))], $
	          /ylog,yrange=[ymin,ymax], $
	          ytickformat = '(f7.3)'

	         ; add line through zero
	         oplot,[!X.CRANGE(0),!X.CRANGE(1)],[0,0],linestyle=0,color=black

	         starttemp = n_elements(symbols)*np
	         endtemp = starttemp + n_elements(symbols) - 1
	         if (endtemp ge n_elements(utemp)) then begin
	          endtemp = n_elements(utemp)-1
	         endif

	         tempslice = utemp[starttemp:endtemp]

	         xloc = xlo + xrng*0.05
	         yloc = yhi

	         ; add data points for each temperature to plot
	         for i=0, n_elements(tempslice)-1 do begin
	          itemp = where(dettemp eq tempslice[i])
	          xdata = hour[itemp]
	          ydata = slope[itemp]
	          ;sort data for plotting
	          ptorder = sort(xdata)
	          xdata = xdata[ptorder]
	          ydata = ydata[ptorder]
	          case i mod n_elements(symbols) OF
	              6: begin
	                  usersym, _X, _Y, /FILL    ; Define the symbol
	              end
	              7: begin
	                  usersym, [-1, -0.5, 0.5, 1, 0.5, -0.5, -1], $
	                      [0, 0.5, 0.5, 0, -0.5, -0.5, 0], /fill
	              end
	          else:
	          endcase

	          ; draw curves
	          oplot,xdata,ydata,psym=symbols[i mod n_elements(symbols)], $
	              color=black,symsize=symsize

	          ; make legend
	          offset=0.15
	          oplot,[!X.CRANGE(0)+(!x.CRANGE(1)-!x.CRANGE(0))*0.86], $
	              [1.5*10^(!Y.CRANGE(0)+offset*i)], $
	              psym=symbols[i mod n_elements(symbols)], color=black,$
	              symsize=symsize
	          xyouts,[!X.CRANGE(0)+(!x.CRANGE(1)-!x.CRANGE(0))*0.88], $
	              [1.4*10^(!Y.CRANGE(0)+offset*i)], $
	              strtrim(string(tempslice[i]),2) + ' K', color=black
	         endfor

			; add legend
			offset=(!y.crange(1)-!y.crange(0))*0.05
			xyouts, [!X.CRANGE(0)+(!x.CRANGE(1)-!x.CRANGE(0))*0.05], [10^(!Y.CRANGE(0)+16*offset)], 'Start date: ' + startdate, color=black
			xyouts, [!X.CRANGE(0)+(!x.CRANGE(1)-!x.CRANGE(0))*0.05], [10^(!Y.CRANGE(0)+17*offset)], 'Refsub Type: ' + refcor_type, color=black
			xyouts, [!X.CRANGE(0)+(!x.CRANGE(1)-!x.CRANGE(0))*0.05], [10^(!Y.CRANGE(0)+18*offset)], 'Refsub Version: ' + refcor_version, color=black
			xyouts, [!X.CRANGE(0)+(!x.CRANGE(1)-!x.CRANGE(0))*0.05], [10^(!Y.CRANGE(0)+19*offset)], 'Analysis Region: ' + region,  color=black

	         ; add file name to bottom of plot
	         xyouts, 1.0, 40.0, fsub, color=black, charsize = filenamecharsize, /device

	         ; write JPG image file
	         if (plotcount eq 1) then begin
	          xyouts, 1.0, 0.0, '_', /device
	          jpgimg = tvrd()
	          write_jpeg, jpgfile, $
	              congrid(jpgimg, 1280, 960, /interp, /center), quality=100
	         endif

		endif

         device, /close

       endfor

    endfor

	; now make plot versus temperature

	; choose the last slope per temperature for plotting so that the value should be settled
	maxdettemp=fltarr(n_elements(utemp))
	maxslope=fltarr(n_elements(utemp))
	for utempindex=0,n_elements(utemp)-1 do begin
		maxdettemp(utempindex)=dettemp(max(where(utemp(utempindex) eq dettemp)))
		maxslope(utempindex)=slope(max(where(utemp(utempindex) eq dettemp)))
	endfor

	; define forground and background colors
	white = 255
	black = 0

	; define plot file names
	psfile = outdir+'darkjumbo_versustemp_'+scaname + '_' + ext+'.ps'
	jpgfile = outdir+'darkjumbo_versustemp_'+scaname + '_' + ext+'.jpg'

	for plotcount=0,1 do begin
		if (plotcount eq 0) then begin

			;set up for landscape Postscript plotting
			set_plot, 'ps'
			!p.font=-1
			device, filename= psfile, /landscape
			fsub = psfile     ;filename to be written at bottom of plot
			!x.minor = 10
			!y.minor = 10
			!y.margin=[6,2]
			!x.margin=[12,4]

			; set line and text characteristics for publication quality
			!P.CHARSIZE=1.2
			!P.CHARTHICK=4.
			!y.charsize = 1.0
			!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 up to plot to memory buffer.  This plot will later be
			; written from memory to a JPEG file.
			set_plot, 'z'
			device, set_resolution=[8000,6000]
			device, set_font='Courier'
			fsub = jpgfile   ; filename to be written at bottom of plot

			;set plot margin & tick marks
			!Y.MINOR=10
			!x.minor=10
			!y.margin=[6,4]
			!x.margin=[15,4]

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

		endelse

		; The series will be represented by big filled circles
		_radius = 1.0 ; 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

		; make plot axes
		gtzero=where(maxslope gt 0,countgtzero)
		if countgtzero gt 0 then begin
			ymin=10^float(round(alog10(min(maxslope(gtzero))*0.9)-0.5))
			ymax=10^float(round(alog10(max(maxslope(gtzero))*1.05)+0.5))
			plot,[0],[0],/nodata,xrange=[min(maxdettemp(where(maxslope ne 0)))*0.9,max(maxdettemp(where(maxslope ne 0)))*1.05],$
				xstyle=1,ystyle=1, xtitle='Temperature (K)',ytitle='Dark Current'+units, $
				title='Dark Current for '+scaname, color=black, background=white,charsize=!P.CHARSIZE*1.4, $
				yrange=[ymin,ymax], /ylog,ytickformat = '(f8.3)'

			; add line through zero
			oplot,[!X.CRANGE(0),!X.CRANGE(1)],[0,0],linestyle=0,color=black

			starttemp = n_elements(symbols)*np
			endtemp = starttemp + n_elements(symbols) - 1
			if (endtemp ge n_elements(utemp)) then begin
				endtemp = n_elements(utemp)-1
			endif

			xloc = xlo + xrng*0.05
			yloc = yhi

			; Make a vector of 16 points, A[i] = 2pi/16:
			A = FINDGEN(17) * (!PI*2/16.)
			; Define the symbol to be a unit circle with 16 points,
			; and set the filled flag:
			USERSYM, COS(A), SIN(A), /FILL

			; add data points for each temperature to plot
	;		itemp = where(dettemp eq tempslice[i])
			xdata = dettemp[itemp]
			ydata = slope[itemp]
			;sort data for plotting
			ptorder = sort(xdata)
			xdata = xdata[ptorder]
			ydata = ydata[ptorder]

			; draw curves
			oplot,maxdettemp(where(maxslope ne 0)),maxslope(where(maxslope ne 0)),psym=8, color=black,symsize=symsize

			; add legend
			offset=(!y.crange(1)-!y.crange(0))*0.05
			xyouts, [!X.CRANGE(0)+(!x.CRANGE(1)-!x.CRANGE(0))*0.05], $
				[10^(!Y.CRANGE(0)+16*offset)], 'Start date: ' + startdate, $
				color=black
			xyouts, [!X.CRANGE(0)+(!x.CRANGE(1)-!x.CRANGE(0))*0.05], $
				[10^(!Y.CRANGE(0)+17*offset)], 'Refsub Type: ' + $
				refcor_type, color=black
			xyouts, [!X.CRANGE(0)+(!x.CRANGE(1)-!x.CRANGE(0))*0.05], $
				[10^(!Y.CRANGE(0)+18*offset)], 'Refsub Version: ' + $
				refcor_version, color=black
			xyouts, [!X.CRANGE(0)+(!x.CRANGE(1)-!x.CRANGE(0))*0.05], $
				[10^(!Y.CRANGE(0)+19*offset)], 'Analysis Region: ' + $
				region,  color=black

			; add file name to bottom of plot
			xyouts, 1.0, 40.0, fsub, color=black, charsize = filenamecharsize, /device

			; write JPG image file
			if (plotcount eq 1) then begin
				xyouts, 1.0, 0.0, '_', /device
				jpgimg = tvrd()
				write_jpeg, jpgfile, $
				congrid(jpgimg, 1280, 960, /interp, /center), quality=100
			endif

		endif

		device, /close

	endfor

    ; restore plotting device to original
    set_plot,devicesave

END