;+
; Name: aia_firdem_test_2
;
; Performs AIA test inversions for a DEM consisting of a log-normal distribution with central
; temperature set by temp_in. Used by aia_firdem_test_script_4.pro to test the 
; DEM inversion over a range of temperatures.
;
; Inputs:
;	-temp_in: The central temperature of the lognormal
;	-niter: Maximum number of iterations for the DEM inversion (optional)
; Outputs:
;	-t2: Temperature grid of the returned output DEM
;	-dem_in: The input DEM generated by the procedure
;
; Optional Arguments:
;
;	-tr_struct: Temperature response, as returned by aia_get_response.
;	-chi2s: Array containing the final chi-squared for each DEM inversion performed.
;	-emtot: Total EM of the input DEM. Default is 1.0e29.
;	-nx,ny: nx*ny determines the total number of runs (each with a different noise
;			realization) to run on the input DEM. Defaults are 40 each.
;	-times: Total time taken to compute the DEMs.
;	-sigma_logn: Half width of input lognormal DEM. Default is 0.15.
;	-logtmin,logtmax: Minimum and maximum of temperature grid.
;
; Joseph Plowman (plowman@physics.montana.edu) 09-17-12
;-
function aia_firdem_test_2, temp_in, t2, dem_in, tr_struct=tr_struct, chi2s=chi2s, emtot=emtot, nx=nx, ny=ny, times = times, sigma_logn=sigma_logn, logtmin=logtmin, logtmax=logtmax

	common aia_firdem_test_2_static_variables, seed

	if(keyword_set(tr_struct) eq 0) then tr_struct= aia_get_response(/temp,/dn,/evenorm)	
	ichan = where(tr_struct.channels ne "A304")
	channels = tr_struct.channels(ichan)
	tr0 = tr_struct.all(*,ichan)
	logt0=tr_struct.logte
	t0 = 10^logt0
	
	if(n_elements(sigma_logn) eq 0) then sigma_logn = 0.15
	if(n_elements(emtot) eq 0) then emtot = 1.0e29
	
	if(n_elements(nx) eq 0) then nx=40
	if(n_elements(ny) eq 0) then ny=40

	nt = 1000
	nchan = n_elements(channels)
	exptimes = 1.0+dblarr(nchan)
	tr = dblarr(nt,nchan)
	if(n_elements(logtmin) eq 0) then logtmin = 5.5
	if(n_elements(logtmax) eq 0) then logtmax = 7.5
	tmin = 10.0^logtmin
	tmax = 10.0^logtmin

	logt = logtmin+(logtmax-logtmin)*findgen(nt)/(nt-1)
	t=10^logt
	temps=t
	logt_in = alog10(temp_in)

	images0 = dblarr(nx,ny,nchan)
	images = dblarr(nx,ny,nchan)

;	itarr = where(abs(temp_in-temps) eq min(abs(temp_in-temps)))
;	it = itarr(0)

	dem_in = dblarr(nt)
	dem_in = emtot*exp(-(logt-logt_in)^2/(2*(sigma_logn)^2))/sqrt(4*acos(0)*(sigma_logn)^2)

	images = aia_firdem_synthesize_data(logt, dem_in, channels, exptimes, tr_struct=tr_struct, $
			nx=nx, ny=ny, errs=errs, data0=images0)		
		
;	for i=0,nchan-1 do begin
;		tr02 = spl_init(t0,tr0(*,i),/double)
;		tr(*,i) = spl_interp(t0,tr0(*,i),tr02,temps)
;		integrand = dem_in*tr(*,i)
;		images0(*,*,i) += int_tabulated(logt,integrand)
;	endfor
	
;	G=dnperphot(channels)
;	errs = dblarr(nchan)
;	errors_array = images0
;	for i=0,nchan-1 do begin
;		G = dnperphot(channels(i))
;		G=G(0)
;		rdn = read_noise(channels(i))
;		rdn = rdn(0)
;		for j=0,nx-1 do begin
;			for k=0,ny-1 do begin
;				photons = randomu(seed,poisson=images0(j,k,i)/G)
;				images(j,k,i) = photons*G+randomn(seed)*rdn
;			endfor
;;			for i=0,nchan-1 do images(j,k,i) = randomu(seed,poisson=images0(j,k,i))
;;		endfor
;		photons = randomu(seed,200,200,poisson=images0(0,0,i)/G)		
;		errs(i) = stddev(photons*G+randomn(seed,200,200)*rdn)
;		errors_array(*,*,i)=errs(i)
;;		message,"Breakpoint!"
;	endfor
	
	estimated_errors = aia_firdem_estimate_error(channels,images)
	
	error_estimate_average_deviation = dblarr(nchan)
	error_estimate_average_spread = dblarr(nchan)
	
	for i=0,nchan-1 do begin
		error_estimate_average_deviation(i) = mean(estimated_errors(*,*,i)-errs(i))/errs(i)
		error_estimate_average_spread(i) = stddev(estimated_errors(*,*,i)-errs(i))/errs(i)
	endfor
	
;	print,error_estimate_average_deviation
;	print,error_estimate_average_spread
;	print,errs
	
	tstart=systime(1)
	aia_firdem_wrapper, channels, exptimes, images, dem_out, tr_struct=tr_struct, chi2s = chi2s, logtlo=logtmin, logthi=logtmax
	times=systime(1)-tstart

	firdem_get_pixel,dem_out,0,0,t2,dem2
	
	return,dem2	
	
end