Part II

@rinit ;initialize the IDL session for RADYN
readradyn ; reads RADYN output with the default name 'radyn_out.cdf
readradyn, 'radyn_out.F11_d4_Ec25_10s.cdf' ; reads RADYN output file with specific name
readradyn,'radyn_out.cdf', dt_interval=1.0 ; reads that RADYN output file but only reading every 1 s.
readradyn
help ; list all the RADYN variables 
help, z1t, timet, ndep
print,interpol(taut[*,0],z1t[*,0],0.)
help,z1t,cmass1t,tg1t,d1t,ne1t,vz1t,n1t,pg1t  ; some of the most useful hydrodynamic 
help,bheat1t,bmom1t,f20t ; a few flare specific variables. 
plot,timet,f20t ; plots the injected beam flux as a function of time.
plot,timet,tg1t[0,*] ; Plot the temperature at the top of the corona as a function of time

*Quesiton 9*

mxmovie, z1t*1e-5, tg1t, /ylog, id = timet ; animated plot of temperature as a function of height
***If this previous command does not work, and if you are running XQuartz 2.7.11 run these commands in the Terminal:
sudo mv /opt/X11/lib/libXt.6.dylib{,.bak}
sudo cp /opt/X11/lib{/flat_namespace,}/libXt.6.dylib
mxmovie, z1t*1e-5, tg1t, yr=[3000,100000], id = timet ; zoom-in on chromospheric temperatures
mxmovie, z1t*1e-5,alog10(tg1t),z1t*1e-5,bheat1t/2000 + 3.5,/ps, id=timet ; plot temperature and beam heating rate as a function height
print,timet[78]
interpol(tg1t[*,78],z1t[*,78]/1e5, 900.)
interpol(tg1t[*,0],z1t[*,0]/1e5, 900.)
mxmovie, z1t*1e-5, n1t[*,8,2,*]/totnt[*,2,*], id = timet ; plot the ionization fraction He III / (He I + He II + He III)

*Question 15*

rad2plot
rad2plot, z1t*1e-5, alog10(tg1t), bheat1t, 6.5 ; plot temperature and beam heating rate at 6.5 s on same plot
rad2plot, z1t*1e-5, alog10(tg1t), bheat1t, 6.5, /oplot_t0,xr=[-100,2000]
rad2plot, z1t/1e5, trl1t, bheat1t, 6.5, xr=[-100,1500],yr=[0,1.5e4],y2r=[0,1.5e4]
rad2plot, z1t/1e5, trl1t + coolt1t, bheat1t, 6.5, xr=[-100,1500],yr=[0,1.5e4],y2r=[0,1.5e4]
for i =0,40 do print,irad[i],jrad[i],ielrad[i],alamb[i],coolt[114,i,78]
print,tg1t[114,78], ne1t[114,78]
mxmovie, z1t*1e-5, vz1t*1e-5,id=timet    ;; look at evolution of gas
rad2plot,z1t/1e5,vz1t/1e5,alog10(pg1t),1.,xr=[-100,1500],yr=[-5,150],y2r=[0,6],/oplot_t0
rad2plot,z1t/1e5,vz1t/1e5,alog10(tg1t),1.,xr=[-100,1500],yr=[-20,150],y2r=[3,7]

rad2plot,z1t/1e5,vz1t/1e5,alog10(tg1t),1.,xr=[-100,1500],yr=[-20,150],y2r=[3,7]
rad2plot,z1t/1e5,vz1t/1e5,alog10(tg1t),6.5,xr=[-100,1500],yr=[-20,150],y2r=[3,7]
rad2plot,z1t/1e5,vz1t/1e5,alog10(pg1t),9.9,xr=[-100,1500],yr=[-5,150],y2r=[1,6],/oplot_t0

*Question 21*
rad2plot,z1t/1e5,vz1t/1e5,alog10(pg1t),0.2,xr=[-100,1500],yr=[-5,120],y2r=[.1,6],/oplot_t0
ad2plot,z1t/1e5,vz1t/1e5,alog10(tg1t),0.1,xr=[800,1500],yr=[-5,120],y2r=[3.5,6],/oplot_t0

*Question 23a*
rad2plot,z1t/1e5,vz1t/1e5,alog10(pg1t),9.9,xr=[-100,1500],yr=[-5,100],y2r=[1,6],/oplot_t0
rad2plot,z1t/1e5,vz1t/1e5,alog10(d1t),6.5,xr=[-100,1500],yr=[-5,100],y2r=[-15,-7],/oplot_t0
rad2plot, z1t*1e-5, alog10(tg1t), bheat1t, 6.5,xr=[-100,1500] ; plot temperature and beam heating rate at 6.5 s on same plot

*Question 26*

rad_thinline_intensity,'Fe_XXI_1354', intensity1354, contrib1354
plot,timet,intensity1354
rad2plot,z1t*1e-5,contrib1354,vz1t*1e-5,6.5
radxdetailed
radxdetailed, user_mu=4, continuum, line1, line2
help,line1,/str
mxmovie, line1.lam, line1.int, line1.int0, xr= [6560,6566]

device,decom = 0
loadct,39
nlam = n_elements(line1.lam[*,0])
time_index = 78 ; t=6.5s
print,timet[78]
contour,line1.contribf[*,*,time_index],z1t[*,time_index]*1e-5,line1.lam[*,time_index],/fill,nlev = 256, xr = [0,2000],yr = [6566,6560],ytitle='Wavelength (A)',xtitle='Height (km)',xmargin=[10,10], ystyle=1+8
loadct,0,/silent
oplot,line1.int[*,time_index]*2000./(max(line1.int[*,time_index])),line1.lam[*,time_index],color=150,linestyle=2,thick=3
oplot,[0,2000],6562.95 + [0,0],linestyle=1
loadct,39,/silent
oplot,z1t[*,time_index]/1e5,vz1t[*,time_index]/cc * 6562.95 * (-1.0) + 6562.95,color=255,thick=3,linestyle=0
axis,yaxis=1,yr=[3,5],/save,ytitle='log_10 Temperature'
oplot,z1t[*,time_index]*1e-5,alog10(tg1t[*,time_index]),thick=3,color=255,linestyle=3

mxmovie,continuum.lam,continuum.int,continuum.int0,/xlog,/ylog,id=timet
plot,continuum.lam[*,78],continuum.int[*,78]-continuum.int0[*,78],xr=[3200,5200],psym=-1,yr=[0,1e6],/ystyle,/xstyle,xtitle=’wavelength
(angstroms)’, ytitle=’intensity mu=0.95 (erg/s/cm2/sr/Ang)’
print,interpol(continuum.int[*,78]-continuum.int0[*,78],continuum.lam[*,78],3600.)/interpol(continuum.int[*,78]-continuum.int0[*,78],continuum.lam[*,78],4170.)

***Part III***
readradyn,'radyn_out.3.5f11.cdf',dt_int=0.1,/ebal
mxmovie, z1t*1e-5,alog10(tg1t),z1t*1e-5,bheat1t/6000 + 3.5,/ps, id=timet 

loadct,39 & device,decomp=0
emovie,id=timet,xr=[-0.1,1.5],yr=[-1e15,1e15]
mxmovie,z1t/1e5,vz1t/1e5,id=timet,yr=[-100,100]

rad2plot,z1t/1e5,vz1t/1e5,alog10(d1t),1.,xr=[-100,1500],yr=[-100,100],y2r=[-14,-6],/oplot_t0
rad2plot,z1t/1e5,vz1t/1e5,alog10(d1t),6.5,xr=[-100,1500],yr=[-100,100],y2r=[-14,-6],/oplot_t0
rad2plot,z1t/1e5,alog10(tg1t),alog10(d1t),1,xr=[-100,1500],yr=[3,7],y2r=[-14,-6],/oplot_t0
rad2plot,z1t/1e5,alog10(tg1t),alog10(d1t),6.5,xr=[-100,1500],yr=[3,7],y2r=[-14,-6],/oplot_t0
rad2plot,z1t/1e5,alog10(tg1t),alog10(pg1t),6.5,xr=[-100,1500],yr=[3,7],y2r=[1,4],/oplot_t0

plot,z1t[*,10]/1e5,tg1t[*,10],yr=[2500,1e7],/ylog,xr=[-100,1500],/xstyle,/ystyle
oplot,z1t[*,65]/1e5,tg1t[*,65],linestyle=2
print,interpol(alog10(cmass1t[*,10]), tg1t[*,10],9500) ; t=1s
print,interpol(alog10(cmass1t[*,65]), tg1t[*,65],9500)  ; t=6.5s
print,interpol(z1t[*,65]/1e5, tg1t[*,65],1e5) 
print,interpol(z1t[*,65]/1e5, tg1t[*,65],9500.)        
print,interpol(z1t[*,10]/1e5, tg1t[*,10],1e5) 
print,interpol(z1t[*,10]/1e5, tg1t[*,10],9500.)
loadct,39 & device,decomp=0
emovie







Part II:  The high beam flux simulation:

readradyn,'radyn_out.f11.cdf'

readradyn

dem = rdem(t)
plot, alog10(t), dem[*,78],/ylog, xr = [4.5,5.1]
dem_tarr =dem*0.
for i =0,n_elements(timet)-1 do dem_tarr[*,i]=t
mxmovie,  alog10(dem_tarr),dem,yr=[1e22,1e24],xr=[4,8]

radxdetailed
radxdetailed, user_mu=4, continuum, line1, line2
help,line1,/str
mxmovie, continuum.lam, continuum.int,/xlog,/ylog
mxmovie, continuum.lam, continuum.int - continuum.int0, xr= [0,8000]
mxmovie, line1.lam, line1.int, line1.int0, xr= [6560,6566]
bj_ratio =fltarr(n_elements(timet)) 
for i =0,n_elements(timet)-1 do begin bj_ratio[i]=interpol(continuum.int[*,i]-continuum.int0[*,i],continuum.lam[*,i],3500.) /interpol(continuum.int[*,i]-continuum.int0[*,i],continuum.lam[*,i],4170.)

halpha_lc = fltarr(n_elements(timet)) & heii304_lc = fltarr(n_elements(timet))
for i =0,n_elements(timet)-1 do halpha_lc[i] = tsum(line1.lam[*,i],line1.int[*,i])
for i =0,n_elements(timet)-1 do heii304_lc[i] = tsum(line2.lam[*,i],line2.int[*,i])
plot,timet,heii304_lc/max(heii304_lc)
oplot,timet,(halpha_lc-halpha_lc[0])/max(halpha_lc-halpha_lc[0]),linestyle=2

;;plot,line1.lam[*,30],line1.int[*,30],xr = [6560,6566] & oplot,line1.lam[*,30],line1.int0[*,30],linesty = 2
device,decom = 0
loadct,39
nlam = n_elements(line1.lam[*,0])
time_index = 78 ; t=6.5s
print,timet[78]
;contour,line1.contribf[*,*,time_index]*(dzt[*,time_index]#replicate(1,nlam)),z1t[*,time_index]*1e-5,line1.lam[*,time_index],/fill,nlev = ;256, xr = [0,2000],yr = [6560,6566],ytitle='Wavelength (A)',xtitle='Height (km)',xmargin=[10,10], ystyle=1+8
;axis,yaxis=1,yr=[3,5],/save,ytitle='log_10 Temperature'
;oplot,z1t[*,time_index]*1e-5,alog10(tg1t[*,time_index]),thick=3,color=255

contour,line1.contribf[*,*,time_index],z1t[*,time_index]*1e-5,line1.lam[*,time_index],/fill,nlev = 256, xr = [0,2000],yr = [6560,6566],ytitle='Wavelength (A)',xtitle='Height (km)',xmargin=[10,10], ystyle=1+8
loadct,0,/silent
oplot,line1.int[*,time_index]*2000./(max(line1.int[*,time_index])),line1.lam[*,time_index],color=150,linestyle=2,thick=3
oplot,[0,2000],6562.95 + [0,0],linestyle=1
loadct,39,/silent
oplot,z1t[*,time_index]/1e5,vz1t[*,time_index]/cc * 6562.95 * (-1.0) + 6562.95,color=255,thick=3,linestyle=0
axis,yaxis=1,yr=[3,5],/save,ytitle='log_10 Temperature'
oplot,z1t[*,time_index]*1e-5,alog10(tg1t[*,time_index]),thick=3,color=255,linestyle=3




Now, load in the 3.5F11 model (higher heating rate):

readradyn,'../runs/radyn_out.3.5F11_d4_Ec25.finished.cdf',dt_int=0.1,/ebal


loadct,39 & device,decompose=0
rad2plot, z1t/1e5, trl1t+coolt1t, bheat1t, 0.4, xr=[-100,1500],yr=[0,2e4],y2r=[0,2e4]


3.5F11:  Identify where beam heating is greater than radiative lossess
at early time
use emovie to identify Which detailed element important?
There's one that clearly dominates the others.

What downward velocity then develops?  What is its maximum downward
velecity (when, what temperature, and what density?)
What is the column mass at 10,000 K?

Follow the evolution of the downward velocity).  When does it reach
~10,000 K?  What is the density now?
What is the column mass at 10,000 K?

Look at Ca II K.  What is the

Over the summer:  Calculate the FeXXI spectra as a function of time
for each model and compare to observations.


Look at detailed radation in CC, what are the dominant radiative
losses in detail?
LTE pops in CC?
n1t vs. nstart


The beam hits at what height?
Does beam heating go to photosphere?

Then look at photosphere with emovie.  What is contributing to the
heating in teh photosphere.  This is called radiative backwarming.