; This script demonstrates how to create synthetic AIA data from a random test DEM, reconstruct ; the DEM using the firdem code, and plot the results. Repeatedly running the script can ; provide invaluable insight into how well input DEMs can generally be recovered. .compile ../estimate_quantile.pro .compile ../firdem_get_pixel.pro .compile ../firdem.pro .compile aia_firdem_wrapper.pro ; Generate AIA Temperature Response functions: if(n_elements(tr) eq 0) then tr=aia_get_response(/temp,/dn,/evenorm) npeaks=5 logt = tr.logte nt = n_elements(tr.logte) pi = 2*acos(0) ; Generate some Gaussian peaks: tcen_min = 6.0 ; Minimum log_10(T) for Gaussian peaks tcen_max = 7.0 ; Maximum log_10(T) for Gaussian peaks tsig_min = 0.025 ; Minimum standard deviation (in log_10(T)) for Gaussian peaks tsig_max = 0.150 ; Maximum standard deviation (in log_10(T)) for Gaussian peaks etot_min = 0.5e28 ; Minimum total emission for Gaussian peaks etot_max = 1.5e28 ; Maximum total emission for Gaussian peaks cens = tcen_min+(tcen_max-tcen_min)*randomu(seed,npeaks) sigs = tsig_min+(tsig_max-tsig_min)*randomu(seed,npeaks) amps = etot_min+(etot_max-etot_min)*randomu(seed,npeaks) dem = dblarr(nt) for i=0,npeaks-1 do dem += amps(i)*exp(-(logt-cens(i))^2/(2*sigs(i)^2))/sqrt(2*pi*sigs(i)^2) ; Produce AIA data numbers. The subroutine aia_firdem_synthesize_data in aia_firdem_wrapper.pro ; will do the same thing: channels_good = where(tr.channels ne 'A304') ; This channel is unsuitable for DEM analysis channels = tr.channels(channels_good) nchan = n_elements(channels) tresps = tr.all(*,channels_good) dns0 = dblarr(1,1,nchan) for i=0,nchan-1 do dns0(*,*,i) = int_tabulated(logt,dem*tresps(*,i)) G = aia_firdem_dnperphot(channels) rdn = aia_firdem_readnoise(channels) phots0 = dns0/G > 0; Convert to photons nx=50 ; The total number of random noise realizations to produce is nx*ny. ny=50 dns=dblarr(nx,ny,nchan) ; Generate data numbers with noise: for i=0,nchan-1 do dns(*,*,i) = G(i)*randomu(seed,nx,ny,poisson=phots0(i))+rdn(i)*randomn(seed) exptimes = 1.0+dblarr(nchan) ; Compute the DEM: aia_firdem_wrapper, channels, exptimes, dns, dem_struc, tr_struct=tr, telapsed=t ; Plot The DEMs: plot,logt,dem,title=string(t/nx/ny)+" Seconds per DEM" for i=0,4 do begin &$ for j=0,4 do begin &$ firdem_get_pixel,dem_struc,i,j,t,dem_out &$ oplot,alog10(t),dem_out,linestyle=1 &$ endfor &$ endfor legend,["Input","Output"], linestyle=[0,1]