pro ipc_reduce, indir=indir, param=param

if keyword_set (indir) then begin
	cd, indir
endif else begin
	print, 'directory not specified'
	stop
endelse

if keyword_set (param) then parfile=param else parfile='IPC_VSUB.param'

openr, U, parfile, /GET_LUN ;U= LUN unit

; Call parseline to read the temps for the experiment
parseline_new, U, temps, num_temps
print, temps

; Call parseline to read vsub voltage array
parseline_new, U, vsub_volts, num_volts
print, vsub_volts

num_temps=n_elements (temps)
num_vsub_volts=n_elements(vsub_volts)

ipc_temp_vsub=fltarr(num_temps,num_vsub_volts)

file_mkdir,'IPC_results'
results_dir=indir+'\IPC_results'

for i=0, (num_temps-1) do begin

	for j=0, (num_vsub_volts-1) do begin
		cd, indir

		darklist = 'dark_' + strtrim(string(temps[i]), 2) + '_' + strtrim(string(vsub_volts[j]), 2) + '.lst'
		ipclist =  'ipc_' + strtrim(string(temps[i]), 2) + '_' + strtrim(string(vsub_volts[j]), 2) + '.lst'
		print, darklist
		print, ipclist

		readcol, darklist, darkfile, format='(A)'
		readcol, ipclist, ipcfile, format='(A)'

		ipc=float(readfits(ipcfile[0],header,nslice=1))
		dark=float(readfits(darkfile[0],header,nslice=1))

		diff=ipc-dark

		writefits, 'diff_' + strtrim(string(temps[i]), 2) + '_' + strtrim(string(vsub_volts[j]), 2) + '.fits', diff,header

		cd,results_dir
		y0=215
		x0=1319.
		num=450
		ipcs=fltarr(5,5,15,450)
		sumipcs=fltarr(15,450)

		; loop over the 15 repeated patterns across the array
		for k=0,10 do begin
			; loop over 450 central pixels going down the array
			for l=0,num-1 do begin
				x=x0+256.*k
				y=y0+8.*l
				ipcs(0:2,*,k,l)=(diff(x-2:x,y-2:y+2))/(diff(x,y))
				x=x0+256.*k-79
				ipcs(3:4,*,k,l)=(diff(x+1:x+2,y-2:y+2))/(diff(x,y))
				sumipcs(k,l)=round((total(ipcs(*,*,k,l))-1.)*10000.)/100.
			endfor
		endfor
	;stop
		; collapse all 5x5 maps into a single median map
		medipc=median(median(ipcs,dimension=4),dimension=3)
		medipc_percent=round(medipc*10000.)/100.
		medsumipcs=median(sumipcs)
		print,medipc_percent

		openw,1,'ipc_'+ strtrim(string(temps[i]), 2) + '_' + strtrim(string(vsub_volts[j]), 2) + '.txt'
		printf,1,medipc_percent
		close,1

		set_plot,'win'
		surface,medipc_percent,charsize=3,/lego,ztitle='percent'
		set_plot,'ps'
		device, filename = 'ipc_'+ strtrim(string(temps[i]), 2) + '_' + strtrim(string(vsub_volts[j]), 2) + '.ps',/landscape
		!p.font = -1    ; use vector fonts for PS
		!P.MULTI = 0    ; Set up for one plot per page
		cd,current=curdir
		fsub = curdir+'\ipc_'+ strtrim(string(temps[i]), 2) + '_' + strtrim(string(vsub_volts[j]), 2) + '.ps'   ; set text to print at bottom of plot

		; set some size values
		asize = 1.3    ; charsize of annotations
		athick = 4.0    ; charthick of annotations
		lthick = 8.0    ; line thickness
		!x.thick=lthick
		!y.thick=lthick
		!z.thick=lthick

		surface,medipc_percent,charsize=asize*2,/lego,ztitle='percent',thick=lthick,charthick=athick,xtitle='x pixel',ytitle='y pixel'
		; Add labels & fitted line
		; write the file name outside the margin of the plot
		black=0
		xyouts, 1.0, 40, fsub, charsize = asize*0.75, charthick = athick, color = black, /device

		device,/close

		; work on sumipcs
		outliers=where(sumipcs gt median(sumipcs)+stddev(sumipcs,/nan)*2. or sumipcs lt median(sumipcs)-stddev(sumipcs,/nan)*2.,noutliers )
		if (noutliers ne 0) then sumipcs(outliers)=!VALUES.F_NAN
		outliers=where(sumipcs gt median(sumipcs)+stddev(sumipcs,/nan)*2. or sumipcs lt median(sumipcs)-stddev(sumipcs,/nan)*2.,noutliers )
		if (noutliers ne 0) then sumipcs(outliers)=!VALUES.F_NAN
		device, filename = 'sumipc_'+ strtrim(string(temps[i]), 2) + '_' + strtrim(string(vsub_volts[j]), 2) + '.ps',/landscape
		fsub = curdir+'\sumipc_'+ strtrim(string(temps[i]), 2) + '_' + strtrim(string(vsub_volts[j]), 2) + '.ps'   ; set text to print at bottom of plot
		surface,sumipcs,[indgen(15)*256.+x0],[indgen(450)*8.+y0],charsize=asize*2,/lego,ztitle='percent',thick=lthick,charthick=athick,xtitle='x pixel',ytitle='y pixel'
		xyouts, 1.0, 40, fsub, charsize = asize*0.75, charthick = athick, color = black, /device
		device,/close

		; add IPC to array to be plotted
		ipc_temp_vsub(i,j)=medsumipcs
	endfor
endfor
openw,1,'sumipc_total.txt'
printf,1,ipc_temp_vsub
close,1
stop
end