pro euv_dimming_disp,dir,iev,vers,io
;+
; Project     : SDO/AIA,  
;
; Name        : EUV_DIMMING_DISP 
;
; Category    : data analysis of EUV dimming and CME mass, speed, and energy
;
; Explanation : Reads AIA images and performs EUV dimming analysis
;
; Syntax      : IDL>euv_dimming_disp,dir,iev,io
;
; Inputs      : dir,iev,io
;
; Outputs     ; postscrpt plot
;
; History     : 24-Feb-2016, Version 1 written by Markus J. Aschwanden
;
; Contact     : aschwanden@lmsal.com
;-

;________________READ SAVEFILE____________________________
string0,3,iev,iev_str
savefile='dimming_'+iev_str+vers+'.sav'
restore,dir+savefile	
print,'read parameters from = ',dir+savefile
;save,filename=dir+savefile,input,para,exptime_,flux_grid,x_grid,y_grid,$
;       flux_prof,emtot_prof,telog_prof,tsig_prof,chi_prof,wave_,time,dem_grid,$
;       l_cm,rho,angle,area,vol,em_tot,dem_t,nel,te_mk,telog,mass,em_model,$
;       d0,x_lasco,v_lasco,t_lasco,a_model,v_model,x_model,ekin_model,$
;       q_t,q_obs,em_obs,a0_model,t0_model,ta_model,chi_model,ind_fit,h0
;	e_grav,v_esc
nt	=n_elements(time)
tmargin1=input.tmargin1                 ;s
tmargin2=input.tmargin2
rise    =input.rise
decay   =input.decay
duration=rise+decay
cadence =input.cadence
it_start=long(tmargin1/cadence) < (nt-1)
it_peak =long((tmargin1+rise)/cadence) < (nt-1)
it_end  =long((tmargin1+rise+decay)/cadence) < (nt-1)
tstart  =time(it_start)
tpeak   =time(it_peak)
tend    =time(it_end)
r0	=6.96e10 			;cm solar radius
dt	=cadence
em1     =max(emtot_prof,imax)
n_model =n_elements(chi_model) 
indv  =where(v_model gt 0)
statistic,v_model(indv),v_bm,dv_bm            		;mean velocity of CME bulk mass
statistic,ekin_model(indv),e_bm,de_bm
statistic,chi_model(indv),chi_avg,chi_sig
statistic,t_lasco,t_bm,dt_bm

;________________READ LASCO CATALOG__________________________
lasco_file='~/work_global/ecme_lasco.dat'
readcol,lasco_file,flare,dateobs_lasco,start_goes,start_lasco,dt_lasco,h_lasco,$
 v1_lasco,v2_lasco,v3_lasco,m_lasco,e_lasco,skipline=3,format='(I,A,A,A,F,F,I,I,I,F,F)'
ind0	=where(flare eq iev,n)
i	=ind0(0)
x_le	=h_lasco(i)*r0	;distance of leading edge detection at LASCO/C2 
v_le3   =[v1_lasco(i),v2_lasco(i),v3_lasco(i)]*1.e5 ;velocities of leading edge
ind3	=where(v_le3 ne 0)
statistic,v_le3(ind3),v_le,dv_le
hh_le   =strmid(start_lasco(i),0,2)
mm_le   =strmid(start_lasco(i),3,2)
ss_le   =strmid(start_lasco(i),6,2)
t_le    =hh_le*3600.+mm_le*60.+ss_le
t0_le   =t_le-(x_le-d0)/v_le
ind	=where(time ge t0_le)
x_le_	=fltarr(nt)+d0
x_le_(ind)=((time(ind)-t0_le)*v_le)+d0
v_le_	=fltarr(nt)
v_le_(ind)=v_le 
m_le    =m_lasco(i)
e_le    =e_lasco(i)

;_________________DISPLAY________________________________________
form    =1      ;0=landscape, 1=portrait
char	=1.
char2	=1.2
if (io ne 0) then char=0.8    ;character size
plotname=dir+'dimming_'
ref     =''     ;label at bottom of Fig.
unit    =0      ;window number
ct	=5
nt	=n_elements(time)
nwave	=n_elements(wave_)
r0	=6.96e10 ;cm
fignr	=iev_str+vers

;________________________AIA WAVELENGTH FLUXES_________________
fig_open,io,form,char,fignr,plotname,unit
loadct,ct
x1_	=0.1
x2_	=0.5
y2_	=0.96
y1_	=0.80
!p.position=[x1_,y1_,x2_,y2_]
!p.title=' '
ymax0=max(flux_prof)
!x.range=minmax(time/3600.)+[-0.1,+0.1]			;FOV coordinates
!y.range=[0,ymax0]
!x.style=1
!y.style=0
!x.title=''
!y.title='AIA Emission measure [cm!U-5!N]'
plot,time/3600.,fltarr(nt)
xyouts,tpeak/3600.-0.01,0,' Flare peak',size=char,orient=90
xyouts,tstart/3600.-0.01,0,' Flare start',size=char,orient=90
xyouts,tend/3600.-0.01,0,' Flare end',size=char,orient=90
color	=50+150*findgen(6)/6.
for iwave=0,nwave-1 do begin
 wave=wave_(iwave)
 ymax=max(flux_prof(*,iwave))
 col	=color(iwave)
 factor	=1.0
 if (ymax lt ymax0/10.) then factor=10.0
 oplot,(time/3600.),flux_prof(*,iwave)*factor,color=col
 oplot,(tstart/3600.)*[1,1],!y.crange
 oplot,(tpeak/3600.)*[1,1],!y.crange
 oplot,(tend/3600.)*[1,1],!y.crange
 str=string(wave_(iwave),'(I4)')+' A x'+string(factor,'(I2)')
 xyouts_norm,0.7,0.9-0.1*iwave,str,char,0,col
 !noeras=1
endfor
qnorm	=!y.crange(1)*0.9/max(emtot_prof)
oplot,time/3600.,emtot_prof*qnorm,linestyle=2

;________________________ACCELERATION_________________
y2_	=0.77
y1_	=0.61
!p.position=[x1_,y1_,x2_,y2_]
!y.range=[0,max(a_model)*1.2] 
plot,time/3600.,a_model(*,0),ytitle='a [cm/s!U2!N]'
for m=0,n_model-1 do oplot,time/3600.,a_model(*,m)	

;________________________VELOCITY_________________
y2_	=0.58
y1_	=0.42
!p.position=[x1_,y1_,x2_,y2_]
!p.title=' '
!x.title=' '
!y.range=[0,(max(v_model)/1.e5) > (v_le/1.e5)]*1.2
;..........................AIA BULK MASS............................
plot,time/3600.,v_model/1.e5,ytitle='v [km/s]'
for m=1,n_model-1 do oplot,time/3600.,v_model(*,m)/1.e5
str1='v!DBM!N = '+string(v_bm/1.e5,'(I4)')
str2=string(dv_bm/1.e5,'(I4)')+' km/s'
xyouts_norm,0.02,0.8,str1+'+',char
xyouts_norm,0.02,0.8,str1+'_'+str2,char
if (t_bm lt time(nt-1)) then begin
 xyouts,t_bm/3600.,v_bm/1.e5,' BM',size=char
 oplot,(t_bm/3600.)*[1,1],(v_bm*[-1,1])/1.e5,psym=4,thick=3
 oplot,(t_bm/3600.)*[1,1],(v_bm+dv_bm*[-1,1])/1.e5
 oplot,(t_bm/3600.)*[1,1],!y.crange,linestyle=1
endif
;...........................LASCO/C2 LEADING EDGE..................
if (v_le gt 0) then begin
 oplot,time/3600.,v_le_/1.e5,color=128
 str1='v!DLE!N = '+string(v_le/1.e5,'(I4)')
 str2=string(dv_le/1.e5,'(I4)')+' km/s'
 xyouts_norm,0.02,0.9,str1+'+',char,0,128
 xyouts_norm,0.02,0.9,str1+'_'+str2,char,0,128
 xyouts,t_le/3600.,v_le/1.e5,' LE',size=char,color=128
 oplot,(t_le/3600.)*[1,1],(v_le*[-1,1])/1.e5,psym=4,thick=3,color=128
 oplot,(t0_le/3600.)*[1,1],[0.,0.],psym=4,thick=3,color=128
 oplot,(t_le/3600.)*[1,1],(v_le+dv_le*[-1,1])/1.e5,color=128
 oplot,(t_le/3600.)*[1,1],!y.crange,linestyle=1,color=128
 oplot,(t0_le/3600.)*[1,1],!y.crange,color=128,linestyle=2
 oplot,(tstart/3600.)*[1,1],!y.crange
 xyouts,tstart/3600.-0.01,0,' Flare start',size=char,orient=90
 xyouts,t0_le/3600.-0.01,0,' LASCO onset',size=char,orient=90,color=128
endif

;________________________DISTANCE_________________
y2_	=0.39
y1_	=0.23
!p.position=[x1_,y1_,x2_,y2_]
!y.range=[0,(max(x_model)/r0) > (x_le/r0)]*1.2
;..........................AIA BULK MASS............................
plot,time/3600.,x_model(*,0)/r0,ytitle='x [R!Dsun!N]'
for m=1,n_model-1 do oplot,time/3600.,x_model(*,m)/r0
if (t_bm lt time(nt-1)) then begin
 str1	='x!DBM!N = '+string(x_le/r0,'(F5.2)')+' R!Dsun!N'
 xyouts_norm,0.02,0.9,str1,char,0,128
 xyouts,t_bm/3600.,x_le/r0,' BM',size=char
 oplot,(t_bm/3600.)*[1,1],(x_le*[-1,1])/r0,psym=4,thick=3
 oplot,(t_bm/3600.)*[1,1],!y.crange,linestyle=1
endif
oplot,!x.crange,(d0/r0)*[1,1],linestyle=1
;...........................LASCO/C2 LEADING EDGE..................
if (x_le gt 0) then begin
 oplot,time/3600.,x_le_/r0,color=128
 oplot,!x.crange,(x_le/r0)*[1,1],linestyle=1,color=128
 xyouts,t_le/3600.,x_le/r0,' LE',size=char,color=128
 oplot,(t_le/3600.)*[1,1],!y.crange,linestyle=1,color=128
 oplot,(t_le/3600.)*[1,1],(x_le*[-1,1])/r0,psym=4,thick=3,color=128
 oplot,(t0_le/3600.)*[1,1],!y.crange,color=128,linestyle=2
 oplot,(tstart/3600.)*[1,1],!y.crange
 xyouts,tstart/3600.-0.01,0,' Flare start',size=char,orient=90
 xyouts,t0_le/3600.-0.01,0,' LASCO onset',size=char,orient=90,color=128
endif

;________________________EMISSION MEASURE______________
y2_	=0.20
y1_	=0.04
!p.position=[x1_,y1_,x2_,y2_]
!x.title='Time  t[hrs]'
!y.range=[0,max(em_obs)*1.2]
plot,time(min(ind_fit))*[1,1]/3600.,!y.crange,linestyle=2
oplot,time/3600.,em_model(*,0),thick=3
oplot,time/3600.,em_obs(*,0),psym=10,thick=1
for m=1,n_model-1 do begin
 oplot,time/3600.,em_model(*,m),color=50+25*m
 oplot,time/3600.,em_obs,psym=10,color=50+25*m
endfor
itm1	=max(ind_fit) < (nt-1)
oplot,time(itm1)*[1,1]/3600.,!y.crange,linestyle=1
itm2	=min(ind_fit) < (nt-1)
oplot,time(itm2)*[1,1]/3600.,!y.crange,linestyle=1
oplot,(tpeak/3600.)*[1,1],!y.crange
xyouts,tpeak/3600.-0.01,0,' Flare peak',size=char,orient=90
nfit	=n_elements(ind_fit)
fit_str =' Fitting range n!Dfit!N='+string(nfit,'(I3)')
xyouts,time(ind_fit(0))/3600.,!y.crange(1)*0.9,fit_str,size=char

;________________________DEM_____________________________
x1_	=0.60
x2_	=0.95
y2_	=0.96
y1_	=0.61
!p.position=[x1_,y1_,x2_,y2_]
demlog=alog10(dem_t > 1.)
plot,telog,demlog(*,0),xrange=minmax(telog),yrange=minmax(demlog),$
     xtitle='log(T)',ytitle='log(EM [cm!U-5!N K!U-1!N])'
for it=1,nt do oplot,telog,demlog(*,it-1)
telog_med=median(telog_prof)
oplot,telog_med*[1,1],!y.crange

;________________________PARAMETERS________________________
!p.position=[0.6,0.05,0.95,0.55]
!x.range=[0,1]
!y.range=[-1.0,2.0]
!x.style=5
!y.style=5
plot,[0,0],[0,0]
xyouts,0.0,2.0,'Flare # '+iev_str,size=char2
xyouts,0.0,1.9,'Version = '+vers,size=char2
tstart_ut =strmid(input.tstart,0,20)
tstart_sec=anytim(input.tstart,/sec) mod 86400.
tpeak_sec=rise+anytim(input.tstart,/sec) mod 86400.
xyouts,0.0,1.8,'t!Dstart!N='+tstart_ut+' UT',size=char2
rise	=(input.rise)/60.
xyouts,0.0,1.7,'t!Dpeak!N-t!Dstart!N='+string(rise,'(F5.1)')+' min',size=char2
dur	=(input.rise+input.decay)/60.
xyouts,0.0,1.6,'t!Dend!N-t!Dpeakt!N='+string(dur,'(F5.1)')+' min',size=char2
t_dimm	=time(ind_fit(0))
dt_dimm =(t_dimm-tpeak_sec)/60.
xyouts,0.0,1.5,'t!Ddimm!N-t!Dpeak!N='+string(dt_dimm,'(F6.1)')+' min',size=char2
xyouts,0.0,1.4,'dt='+string(dt/60.,'(F4.1)')+' min',size=char2
xyouts,0.0,1.3,'GOES = '+input.goes,size=char2
xyouts,0.0,1.2,'Heliogr. '+input.helpos+' ('+$
       string(rho*180./!pi,'(F4.1)')+' deg)',size=char2
xyouts,0.0,1.1,'FOV ='+string(input.fov0,'(F4.2)')+' R!Dsun!N',size=char2
xyouts,0.0,1.0,'NBIN='+string(input.nbin,'(I4)'),size=char2
xyouts,0.0,0.9,'L='+string(l_cm/1.e8,'(I4)')+' Mm, ('+$
       string(angle,'(I3)')+'!U0!N)',size=char2
xyouts,0.0,0.8,'EM='+string(em_tot,'(E7.1)')+' cm!U-3!N',size=char2
xyouts,0.0,0.7,'n!De!N='+string(nel,'(E7.1)')+' cm!U-3!N',size=char2
xyouts,0.0,0.6,'T='+string(te_mk,'(F5.2)')+' MK',size=char2
xyouts,0.0,0.5,'m='+string(mass,'(E7.1)')+' g [AIA]',size=char2
xyouts,0.0,0.4,'v!DBM!N='+string(v_bm/1.e5,'(I5)')+' +'+$
          string(dv_bm/1.e5,'(I5)')+' km/s [AIA]',size=char2
xyouts,0.0,0.3,'v!Desc!N='+string(v_esc/1.e5,'(I5)')+' km/s',size=char2
statistic,ekin_model/1.e30,ekin_avg,ekin_sig
xyouts,0.0,0.2,'E!Dkin!N='+string(e_bm/1.e30,'(F7.2)')+' +'+$
          string(de_bm/1.e30,'(F6.2)')+' 10!U30!N erg [AIA]',size=char2
xyouts,0.0,0.1,'E!Dgrav!N='+string(e_grav/1.e30,'(F7.2)')+' 10!U30!N erg',size=char2
xyouts,0.0, 0.0,'q!Ddimm!N='+string(1.-em_obs(nt-1)/max(em_obs),'(f5.2)'),size=char2
statistic,chi_model,chi_avg,chi_sig
xyouts,0.0,-0.1,'chi  ='+string(chi_avg,'(f4.1)'),size=char2

;______________________LASCO LEADING EDGE______________________________
if (v_le gt 0) or (m_le gt 0) or (e_le gt 0) then begin
 xyouts,0.25,-0.4,'LASCO :',size=char2,color=128
 dt_lasco=(t0_le-tstart_sec)/60.
 xyouts,0.0,-0.5,'t!DLASCO!N-t_!Dstart!N=  '+string(dt_lasco,'(F5.1)')+' min',size=char2,color=128
 xyouts,0.0,-0.6,'m='+string(m_le,'(E7.1)')+' g',size=char2,color=128
 xyouts,0.0,-0.7,'x!LE!N='+string(x_le/r0,'(F5.2)'),size=char2,color=128
 xyouts,0.0,-0.8,'v!DLE!N='+string(v_le/1.e5,'(I5)')+'+'+$
	string(dv_le/1.e5,'(I5)')+' km/s',size=char2,color=128
 xyouts,0.0,-0.9,'E!Dkin!N='+string(e_le/1.e30,'(F7.2)')+$
	' 10!U30!N erg',size=char2,color=128
 xyouts,0.0,-1.0,'v-ratio='+string(v_bm/v_le,'(f5.2)'),size=char2,color=128
endif
fig_close,io,fignr,ref
end