;+
; Name: eis_dem_test_2
;
; Performs EIS test inversions for a DEM consisting of a log-normal distribution with central
; temperature set by temp_in. Used by eis_test_dem_script_2.pro to test the 
; DEM inversion over a range of temperatures.
;
; Inputs:
;
;	temp_in: The central temperature of the lognormal
;
; Outputs:
;
;	t2: Temperature grid of the returned output DEM
;	dem_in: The input DEM generated by the procedure
;
; Returns the dem inversion for the first noise realization tried.
;
; 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.
;	exptimes: Exposure times to assume. Default is 30 seconds for each line.
;	niter: Maximum number of iterations for the DEM inversion. Default set in fastdem_iterate.
;	other_noise: Estimated fractional error due to unmodeled noise sources. Default: 0.1
;
; Joseph Plowman (plowman@physics.montana.edu) 09-20-12
;-
function eis_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, exptimes=exptimes, niter=niter, logtmin=logtmin, logtmax=logtmax, pix_scales=pix_scales, other_noise=other_noise

	common eis_dem_iterative_test_2_static_variables, seed

	if(n_elements(sigma_logn) eq 0) then sigma_logn = 0.15
	if(n_elements(emtot) eq 0) then emtot = 5.0e28
	
	if(n_elements(nx) eq 0) then nx=40
	if(n_elements(ny) eq 0) then ny=40
	
	if(n_elements(exptime) eq 0) then 	exptime = 30
	
;	if(n_elements(niter) eq 0) then niter = 5000

;	channels_all = tr_struct.eis_windows
;	channels = ["fe_09_188.497", "fe_10_184.537", "fe_12_195.119", "fe_15_284.163", "fe_16_262.976"]
	channels = ["fe_09_188.497", "fe_10_184.537", "fe_12_195.119", "fe_15_284.163", "fe_16_262.976","mg_05_276.579","mg_06_270.394","mg_07_280.737","fe_09_197.862","fe_11_188.216","fe_11_180.401","s_10_264.233","fe_12_192.394","fe_13_202.044","fe_13_203.826","fe_14_270.519","s_13_256.686","ca_14_193.874","ca_15_200.972","ca_16_208.604","ca_17_192.858","si_07_275.361","si_10_258.371","fe_14_264.79"]

	nchan = n_elements(channels)
	if(n_elements(sigma_logn) eq 0) then sigma_logn = 0.15
	if(n_elements(pix_scales) eq 0) then pix_scales = [2.0,1.0]
	if(n_elements(exptimes) eq 0) then exptimes = 30+dblarr(nchan)
	if(n_elements(emtot) eq 0) then emtot = 5.0e28

	if(n_elements(logtmin) eq 0) then logtmin = 5.0
	if(n_elements(logtmax) eq 0) then logtmax = 7.5
	tmin = 10.0^logtmin
	tmax = 10.0^logtmax

	nt = 10*(logtmax-logtmin)/sigma_logn
	logt = logtmin+(logtmax-logtmin)*findgen(nt)/(nt-1)
	logt_in = alog10(temp_in)
	temps = 10^logt
	
	dem_in = emtot*exp(-(logt-logt_in)^2/(2*sigma_logn^2))/sqrt(4*acos(0)*sigma_logn^2)
	
	images = eis_firdem_synthesize_data(logt, dem_in, channels, exptimes, tr_struct=tr_struct, $
		nx=nx, ny=ny, errs=errs, data0=data0, pix_scales=pix_scales, other_noise=other_noise)	
	
	tstart=systime(1)

	tstart=systime(1)
	eis_firdem_wrapper, channels, images, errs, dem_out, tr_struct=tr_struct, chi2s = chi2s,niter=niter,logtlo=logtmin,logthi=logtmax;, mincond2=0.01
	times=systime(1)-tstart
	firdem_get_pixel,dem_out,0,0,t2,dem2
;	plot,temps,dem_in,linestyle=2,/xlog
;	oplot,t2,dem2,linestyle=0
;	plot,t2,dem2,linestyle=0,/xlog, yrange = [min([dem2,dem_in*temps*alog(10)]),max([dem2,dem_in*temps*alog(10)])],xrange=[tmin,tmax]
;	plot,t2,dem2,linestyle=0,/xlog, yrange = [min([dem2,dem_in]),max([dem2,dem_in])],xrange=[tmin,tmax]
;	oplot,temps,dem_in*temps*alog(10),linestyle=2
;	oplot,temps,dem_in,linestyle=2
;	oplot,t2,dem2, linestyle=0
;	legend,["Input","Output"], linestyle=[2,0]
;	oplot,temps,dem_in*temps*alog(10)		
	
;	images2 = dblarr(nchan)
;	tr2 = dblarr(n_elements(t2),nchan)
;	for i=0,nchan-1 do begin
;		tr02 = spl_init(temps,tr(*,i),/double)
;		tr2(*,i) = spl_interp(temps,tr(*,i),tr02,t2)
;		integrand = dem2*tr2(*,i)
;		print,dem2
;		print,tr2(*,i)
;		images2(i) += integrate(alog10(t2),integrand)
;		images2(i) += integrate(t2,integrand)
;	endfor
;	print,"Input Intensities = "
;	print,transpose(images0(0,0,*))
;	print,"Output Intensities= "
;	print,images2
;	print,"Fractional Error = ", (transpose(images0(0,0,*))-images2)/transpose(images0(0,0,*))
;	print,"chi squared = ", total(((transpose(images0(0,0,*))-images2)/transpose(errs(0,0,*)))^2)
;	temps_out = t2
;	dem_out = dem2
;	temps_in = temps
	
	return,dem2
	
end

pro plot_basis,dem_out
	plot, dem_out.t, dem_out.basis(*,0), /xlog, yrange=[0,max(dem_out.basis)], xtitle='Temperature(Kelvin)',title='Normalized EIS Emissivities/Basis Functions'
	for i=0,n_elements(dem_out.basis(0,*))-1 do begin
		oplot,dem_out.t,dem_out.basis(*,i)
	endfor
end