;+ ; Name: aia_firdem_curve_test ; ; This procedure tests AIA DEM inversion of arbitrary curves provided by the calling ; procedure. The test is repeated 2500 times for benchmarking and to examine variability ; due to shot and read noise. 25 Of these results are plotted, and one is returned to the ; calling procedure. ; ; Inputs: ; ; temps_in0: Temperature grid for the input DEM. ; dem_in0: The input DEM curve to be tested. ; ; Outputs: ; ; temps_out: Temperature grid for the output DEM. ; dem_out: The extracted output DEM for one of the noise realizations generated. ; ; Optional Arguments: ; ; niter: The maximum number of iterations in the DEM inversion (otherwise the default ; in aia_firdem_wrapper is used). ; tr_struct: The temperature response structure returned by aia_get_response. ; its: The number of iterations taken for each noise realization. ; chi2s: The chi squared for each noise realization. ; data: AIA channel intensities for each noise realization. ; seed: Optional input seed. ; name: String to add to plot title. ; ; Joseph Plowman (plowman@physics.montana.edu) 09-17-12 ;- pro aia_firdem_curve_test, temps_in0, dem_in0, temps_out, dem_out, niter=niter, tr_struct=tr_struct,its=its,chi2s=chi2s, data=data, seed=seed,name=name common aia_firdem_curve_test, seed2 if(n_elements(seed) eq 0) then seed=34985739857398 if(n_elements(tr_struct) eq 0) then tr_struct = aia_get_response(/temp,/dn) ichan = where(tr_struct.channels ne "A304") channels = tr_struct.channels(ichan) nchan = n_elements(channels) tr0 = tr_struct.all(*,ichan) logt0=tr_struct.logte t0 = 10^logt0 nx=50 ny=50 exptimes = 1.0+dblarr(nchan) nt = 481 tr = dblarr(nt,nchan) logtmin = 5.5 if(logtmin lt min(temps_in0)) then logtmin = min(temps_in0) logtmax = 7.5 if(logtmax gt max(temps_in0)) then logtmax = max(temps_in0) tmin = 10.0^logtmin tmax = 10.0^logtmin logt = logtmin+(logtmax-logtmin)*findgen(nt)/(nt-1) t=10^logt dem_in = dblarr(nt) dem_in02 = spl_init(temps_in0,dem_in0,/double) dem_in = spl_interp(temps_in0,dem_in0,dem_in02,logt) images = aia_firdem_synthesize_data(logt, dem_in, channels, exptimes, tr_struct=tr_struct, $ nx=nx, ny=ny, errs=errs, data0=images0) sigmas = dblarr(nchan) for i=0,nchan-1 do begin sigmas(i)=stddev(images(*,*,i)) endfor aia_firdem_wrapper, channels, exptimes, images, dem_out, niter=niter, tr_struct=tr_struct, $ its=its,chi2s=chi2s, telapsed=telapsed firdem_get_pixel,dem_out,0,0,t2,dem2,/zeronegs qlvls = [0.25,0.5,0.75] nqlvls = n_elements(qlvls) ymin = min([min(dem2),min(dem_in)]) ymax = 1.1*max([max(dem2),max(dem_in)]) yrange = [ymin,ymax] chi2_95 = estimate_quantile(chi2s,0.95) title = string(format='(%"%7.1e%s%7.1e")',telapsed/nx/ny," Seconds per DEM; 95% Chi Squared=",chi2_95) if(n_elements(name) gt 0) then title = name+title plot,alog10(t),dem_in,linestyle=0, title=title, xmargin=[14,3], yrange=yrange,xrange=[tmin,tmax],xtitle="Log10(T)",ytitle="DEM [cm^-5/Log10(T)]" nxp = min([nx,5]) nyp = min([ny,5]) ii=1 print,"Plotted chi squareds: " plotted_chi2s=dblarr(nxp*nyp) for i=0,nxp-1 do begin for j=0,nyp-1 do begin ii++ firdem_get_pixel,dem_out,i,j,t2,dem2,/zeronegs oplot,alog10(t2),dem2,linestyle=1,color=128 plotted_chi2s(i*nyp+j) = chi2s(i,j) endfor endfor print,plotted_chi2s legend,["Input","Output"], linestyle=[0,1] oplot,alog10(t),dem_in,linestyle=0 negs = where(dem2 lt 0,count) if (count gt 0) then dem2(negs) = 0 temps_out = t2 dem_out = dem2 temps_in = t end