; This script computes and plots the response of our fast DEM method to a set of log-normal ; distributions, with specified width and total em ('width' and 'emtot' below), at temperatures ; ranging from 10^5.5 to 10^7.0 Kelvin. .compile ../estimate_quantile.pro .compile ../dist_median.pro .compile ../firdem.pro .compile ../firdem_get_pixel.pro .compile aia_firdem_wrapper.pro .compile aia_firdem_test_2.pro logtmin = 5.5 logtmax = 7.0 tmin = 10^logtmin tmax = 10^logtmax nt = 100 nclvl = 8 nx = 25 ny = 40 emtot = 5.0e28 width = 0.2 tempsin = tmin*(tmax/tmin)^(findgen(nt)/(nt-1)) demarr = aia_firdem_test_2(tempsin(0), tempsout, dem_in, tr_struct=tr_struct, logtmin=logtmin-0.5, logtmax=logtmax+0.5, sigma_logn=width*0.5, nx=nx, ny=ny ) itout_min = value_locate(tempsout,tmin) itout_max = value_locate(tempsout,tmax) nt_out = n_elements(tempsout) demarr = dblarr(nt,nt_out) chi2arr = dblarr(nt,nx,ny) times = dblarr(nt) for i=0,nt-1 do begin &$ demarr(i,*) = aia_firdem_test_2(tempsin(i),tempsout,dem_in, tr_struct=tr_struct, chi2s=chi2s, nx=nx, ny=ny, times=time, emtot=emtot, logtmin=logtmin-0.5, logtmax=logtmax+0.5, sigma_logn=width*0.5) &$ times(i) = time &$ chi2arr(i,*,*) = chi2s &$ print,i, nt &$ endfor logtavgs = dblarr(nt) logtout = alog10(tempsout) for i=0,nt-1 do logtavgs(i) = dist_median(logtout,demarr(i,*)>0) nt_out = n_elements(tempsout) logtin = alog10(tempsin) logtout = alog10(tempsout) ia = lindgen(nt_out) ishifts = lonarr(nt_out) for i=0,nt-1 do ishifts(i) = ia(where(abs(tempsout-tempsin(i)) eq min(abs(tempsout-tempsin(i))))) negdemarr = dblarr(nt,nt_out) posdemarr = dblarr(nt,nt_out) clrdemarr = dblarr(nt,nt_out,3) ipos = where(demarr gt 0,count) if count gt 0 then posdemarr(ipos) = demarr(ipos) clrdemarr(*,*,0) = posdemarr negs = where(demarr lt 0,count) if count gt 0 then negdemarr(negs) = - demarr(negs) clrdemarr(*,*,2) = negdemarr xscale = alog10(tempsin(1))-alog10(tempsin(0)) yscale = alog10(tempsout(1))-alog10(tempsout(0)) x0 = alog10(tempsin(0)) y0 = alog10(tempsout(0)) yrange = [alog10(min(tempsout)),alog10(max(tempsout))] yt = 'Extracted Log!D10!N(T) (Kelvin)' xt = 'Injected Log!D10!N(T) (Kelvin)' set_plot,'ps' emtot_str=string(format='(%"%7.1e")',emtot) width_str=string(format='(%"%g")',width) device,file="plots/aia_gresp_width"+width_str+"_"+emtot_str+"EM.eps",/encapsulated !p.font=0 device,/inches,xsize=5.6,ysize=1.6,/times !p.multi=[0,3,1] cmax = emtot/sqrt(4*acos(0)*(0.5*width)^2) ;cmax = max(smooth(max(posdemarr,dimension=2),round(0.25*nt))) plot_image,posdemarr, min=0, max=cmax, scale=[xscale,yscale], origin=[x0,y0], xtitle = xt, ytitle = yt, title='AIA DEM Gaussian Response', yrange=yrange,/nosquare clevels = min(demarr)+(cmax-min(demarr))*findgen(nclvl)/(nclvl-1) contour, demarr,alog10(tempsin),alog10(tempsout),levels=clevels,/overplot,c_colors=64 oplot,alog10(tempsin),alog10(tempsin),linestyle=2,color=255,thick=2 oplot,logtin,logtavgs,linestyle=0,color=0,thick=6 oplot,logtin,logtavgs,linestyle=0,color=255,thick=2 legend,'Extracted EMWMT',linestyle=0,color=0,textcolors=255,charsize=0.45,box=0,thick=6 legend,'Extracted EMWMT',linestyle=0,color=255,textcolors=255,charsize=0.45,box=0,thick=2 title = 'EM=' + emtot_str +' cm!E-5!N' + ', Width='+width_str + ' Dex' qlvls = [0.95,0.9,0.75,0.5] qlvllabels = ['95th percentile','90th percentile','75th percentile','50th percentile'] nqlvls = n_elements(qlvls) qlvlcurves = dblarr(nt, nqlvls) for i=0,nt-1 do for j=0,nqlvls-1 do qlvlcurves(i,j) = estimate_quantile(chi2arr(i,*,*),qlvls(j)) plot,alog10(tempsin),qlvlcurves(*,0),linestyle=0,title = title,xtitle=xt,ytitle="Chi squared",yrange=[0,20] for i=0,nqlvls-1 do oplot,alog10(tempsin),qlvlcurves(*,i),linestyle=i legend,qlvllabels,linestyle=indgen(nqlvls),charsize=0.55,box=0 plot,alog10(tempsin),times/nx/ny,title="Average Time per DEM",xtitle=xt,ytitle="Time (Seconds)" device,/close set_plot,'x'