;events
;
;2011/02/15T01:56:00  -->2011046
;2011/03/09T23:23:00  -->2011068
;2011/07/30T02:09:00  -->2011211



fi='/scr/s07/kawate/CDAW13/EVE_data/EVL_L2_2011046_00_003_01.fit.gz'
d=eve_read_whole_fits(fi)

w=string(d.linesmeta[*].wave_center)
n=d.linesmeta[*].name
t=string(d.linesmeta[*].logt)
b=d.linesmeta[*].blends
u=string(d.linesmeta[*].wave_max)
l=string(d.linesmeta[*].wave_min)



line=[w[0],n[0],t[0],b[0],u[0],l[0]]
for i=1,n_elements(w)-1 do $
  line=[[line],[w[i],n[i],t[i],b[i],u[i],l[i]]]


;EVLラインリスト
      9.39260 Fe XVIII       6.81000                     9.43000       9.33000
      13.1240 Fe VIII        5.57000 Fe VIII             13.1700       13.0400
      13.2850 Fe XX          6.97000 Fe XXIII            13.3200       13.2300
      17.1070 Fe IX          5.81000                     17.2400       17.0200
      17.7243 Fe X           5.99000 Fe VII              17.8300       17.6300
      18.0407 Fe XI          6.07000 Fe X, Fe VII        18.1500       17.9600
      19.5120 Fe XII         6.13000                     19.6100       19.4300
      20.2044 Fe XIII        6.19000 Fe XI, Fe XII       20.3200       20.1400
      21.1331 Fe XIV         6.27000 Fe XII              21.2000       21.0700
      25.6317 He II          4.75000                     25.6800       25.5500
      28.4150 Fe XV          6.30000                     28.5000       28.3000
      30.3783 He II          4.70000                     30.5000       30.2500
      33.5410 Fe XVI         6.43000                     33.6100       33.4700
      36.0758 Fe XVI         6.43000                     36.2000       36.0200
      36.8076 Mg IX          5.99000 Fe VIII-XIV         36.8900       36.7100
      46.5221 Ne VII         5.71000                     46.7400       46.3200
      49.9406 Si XII         6.29000                     50.0400       49.8400
      52.5795 O III          4.92000 O II                52.7200       52.4200
      55.4370 O IV           5.19000 O IV                55.6400       55.2000
      58.4334 He I           4.16000                     58.6800       58.2200
      59.9598 O III          4.92000                     60.1400       59.8400
      62.4943 Mg X           6.05000                     62.6800       62.2800
      62.9730 O V            5.37000                     63.1800       62.7400
      71.8535 O II           4.48000 O II                72.0000       71.7200
      77.0409 Ne VIII        5.81000                     77.1800       76.9000
      79.0199 O IV           5.19000 O III, O IV         79.1400       78.9000
      97.2537 H I            3.84000                     97.4400       97.0800
      97.7030 C III          4.84000                     97.8600       97.5600
      102.572 H I            3.84000                     102.700       102.420
      103.190 O VI           5.47000                     103.320       103.100


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;flux積分

fi='/scr/s07/kawate/CDAW13/EVE_data/EVS_L2_2011046_00_003_01.fit.gz'
d=eve_read_whole_fits(fi)

ti=0
lnum=n_elements(line[0,*])

flux=fltarr(lnum)
eflux=fltarr(lnum)
hwhm=fltarr(lnum)


wl=d.spectrummeta.wavelength
dwl=wl[1:n_elements(wl)-1]-wl[0:n_elements(wl)-2]


for i=0,lnum-1 do begin
  if line[1,i] eq 'O IV    ' then continue
  win=where((wl ge line[5,i]) and (wl le line[4,i]))
  
  flux[i]=total(d.spectrum[ti].irradiance[win]*dwl[win])
  eflux[i]=total(d.spectrum[ti].irradiance[win]*d.spectrum[ti].precision[win]*dwl[win])
  hm=max(d.spectrum[ti].irradiance[win])*0.5
  hwhm[i]=(max(wl[win[where(d.spectrum[ti].irradiance[win] ge hm)]]) $
           - min(wl[win[where(d.spectrum[ti].irradiance[win] ge hm)]]))*0.5


endfor

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;単位変換

;W/m2 @1AU -> erg/s/cm2/sr
flux=flux*1e3*(1.5e11)^(-2)*(6.96e10^2*2*!pi) *1e6
eflux=eflux*1e3*(1.5e11)^(-2)*(6.96e10^2*2*!pi) *1e6


;nm -> A
hwhm=hwhm;*10
lind=float(line[0,*])*10

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;ファイル作成

obs=[string(lind[0]),string(flux[0]),string(eflux[0]),string(hwhm[0]),line[1,0]]
for i=1,lnum-1 do begin
  if hwhm[i] eq 0 then continue
  obs=[[obs],[string(lind[i]),string(flux[i]),string(eflux[i]),string(hwhm[i]),line[1,i]]]
endfor


openw,1,'dem_in20110215.obs'
printf,1,obs
close,1


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;DEM計算

chianti_dem,output='dem_in20110215',file_input='dem_in20110215.obs',pressure=3.e15, mesh_points=[4,4.5,5,5.5,6,6.5,7]

;ionization equiilbrium file選択


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;保留;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;white light flareにおけるEUV continuum
;H I recombination  --> ~95nm


fi='/scr/s07/kawate/CDAW13/EVE_data/EVS_L2_2011046_01_003_01.fit.gz'
d=eve_read_whole_fits(fi)
plot,d.spectrummeta.wavelength,d.spectrum[150].irradiance>0, $
  /ylog,yrange=[1e-7,1e-3],xtitle='wavelength [nm]', $
  ytitle='irradiance [W/m2/nm]',charsize=1.5

oplot,d.spectrummeta.wavelength,d.spectrum[333].irradiance>0,color=cred

xyouts,5,5,'2011/02/15 01:25:00',/dev,charsize=1.5
xyouts,5,25,'2011/02/15 01:55:30',/dev,charsize=1.5,color=cred

write_png,'20110215eve_sp.png',tvrd(true=1)

;;;;;

plot,d.spectrummeta.wavelength,d.spectrum[333].irradiance/ $
  d.spectrum[150].irradiance,/ylog,xtitle='wavelength [nm]', $
  ytitle='irradiance ratio of flare peak to quiet',charsize=1.5

plots,[0,120],[1,1],color=cred

write_png,'20110215eve_sp_r.png',tvrd(true=1
;;;;;

!p.multi=[0,1,2]
window,xs=600,ys=800
for i=150,359 do begin
  plot,d.spectrummeta.wavelength,d.spectrum[150].irradiance>0, $
    /ylog,yrange=[1e-7,1e-3],xtitle='wavelength [nm]', $
    ytitle='irradiance [W/m2/nm]',charsize=1.5, $
    title=anytim(tai2utc(d.spectrum[i].tai),/ccsds)

  oplot,d.spectrummeta.wavelength,d.spectrum[i].irradiance>0,color=cred

  plot,d.spectrummeta.wavelength,d.spectrum[i].irradiance/ $
    d.spectrum[150].irradiance,/ylog,xtitle='wavelength [nm]', $
    ytitle='irradiance ratio to quiet',charsize=1.5,yrange=[1e-2,1e2]

  plots,[0,120],[1,1],color=cred

  write_png,'20110215eve_sp'+string(i,format='(i3.3)')+'.png',tvrd(true=1)

endfor



obs_summ = hsi_obs_summary() 
obs_summ->set,obs_time_interval=['2011/02/15 1:45', '2011/02/15 2:00']
counts = obs_summ -> getdata() 
obs_summ -> plotman ,/eclipse, /flare,/log












;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;GOESからtemp., emを求める。

;;event 20110215
o=ogoes()
o->set,tstart='15-feb-11 01:00:00', tend='15-feb-11 02:30:00'
da=o->getdata()
ti=o->getdata(/times)
ut=o->get(/utbase)

goes_tem,yclean=da,tempr=temp,emis=em

window,0
utplot,ti,da[*,0],ut,/ylog,yrange=[1e-8,1e-3],title='GOES low and high WL',ytitle='GOES flux W/m2',charsize=1.5
outplot,ti,da[*,1],ut,line=1

window,1
utplot,ti,temp,ut,title='temperature',ytitle='1e6 Kelvin',charsize=1.5

window,2
utplot,ti,em,ut,title='emission measure',/ylog,ytitle='1e49 cm-3',charsize=1.5

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;background

IDL> print,anytim(anytim(ut,/second)+ti[800],/ccsds)
2011-02-15T01:27:18.389
IDL> print,anytim(anytim(ut,/second)+ti[920],/ccsds)
2011-02-15T01:31:24.145


IDL> print,mean(da[800:920,0])
  7.53229e-07
IDL> print,mean(da[800:920,1])
  1.69097e-08


lbk=mean(da[800:920,0])
hbk=mean(da[800:920,1])

da2=da
da2[*,0]=da[*,0]-lbk
da2[*,1]=da[*,1]-hbk


goes_tem,yclean=da2,tempr=temp,emis=em



window,1
utplot,ti,temp,ut,title='temperature',ytitle='1e6 Kelvin',charsize=1.5
plots,[ti[910],ti[910]],[0,25],color=cred
xyouts,ti[910],2,anytim(anytim(ut,/seconds)+ti[910],/ccsds),color=cred

write_png,'20110215temp.png',tvrd(true=1)


window,2
utplot,ti,em,ut,title='emission measure',/ylog,ytitle='1e49 cm-3',charsize=1.5
plots,[ti[910],ti[910]],[1e-3,1e2],color=cred
xyouts,ti[910],50,anytim(anytim(ut,/seconds)+ti[910],/ccsds),color=cred



write_png,'20110215em.png',tvrd(true=1)


save,ti,ut,temp,em,da,file='20110215goes_tem.sav'

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;event 20110804

o=ogoes()
o->set,tstart='04-aug-11 03:00:00', tend='04-aug-11 04:30:00'
da=o->getdata()
ti=o->getdata(/times)
ut=o->get(/utbase)

goes_tem,yclean=da,tempr=temp,emis=em

window,0
utplot,ti,da[*,0],ut,/ylog,yrange=[1e-9,1e-3],title='GOES low and high WL',ytitle='GOES flux W/m2',charsize=1.5
outplot,ti,da[*,1],ut,line=1

write_png,'20110804goesHL.png',tvrd(true=1)


window,1
utplot,ti,temp,ut,title='temperature',ytitle='1e6 Kelvin',charsize=1.5

window,2
utplot,ti,em,ut,title='emission measure',/ylog,ytitle='1e49 cm-3',charsize=1.5

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;background


IDL> print,anytim(anytim(ut,/second)+ti[1050],/ccsds)
2011-08-04T03:35:50.472
IDL> print,anytim(anytim(ut,/second)+ti[1150],/ccsds)
2011-08-04T03:39:15.269


IDL> print,mean(da[1050:1150,0])
  5.23137e-07
IDL> print,mean(da[1050:1150,1])
  1.69261e-09


lbk=mean(da[1050:1150,0])
hbk=mean(da[1050:1150,1])

da2=da
da2[*,0]=da[*,0]-lbk
da2[*,1]=da[*,1]-hbk


goes_tem,yclean=da2,tempr=temp,emis=em



window,1
utplot,ti,temp,ut,title='temperature',ytitle='1e6 Kelvin',charsize=1.5
plots,[ti[1170],ti[1170]],[0,20],color=cred
xyouts,ti[1170],2,anytim(anytim(ut,/seconds)+ti[1170],/ccsds),color=cred

write_png,'20110804temp.png',tvrd(true=1)


window,2
utplot,ti,em,ut,title='emission measure',/ylog,ytitle='1e49 cm-3',charsize=1.5
plots,[ti[1170],ti[1170]],[1e-3,1e1],color=cred
xyouts,ti[1170],5e-3,anytim(anytim(ut,/seconds)+ti[1170],/ccsds),color=cred



write_png,'20110804em.png',tvrd(true=1)


save,ti,ut,temp,em,da,file='20110804goes_tem.sav'



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;event 20120705

o=ogoes()
o->set,tstart='05-jul-12 03:00:00', tend='05-jul-12 04:00:00'
da=o->getdata()
ti=o->getdata(/times)
ut=o->get(/utbase)

goes_tem,yclean=da,tempr=temp,emis=em

window,0
utplot,ti,da[*,0],ut,/ylog,yrange=[1e-7,1e-3],title='GOES low and high WL',ytitle='GOES flux W/m2',charsize=1.5
outplot,ti,da[*,1],ut,line=1

write_png,'20120705goesHL.png',tvrd(true=1)


window,1
utplot,ti,temp,ut,title='temperature',ytitle='1e6 Kelvin',charsize=1.5

window,2
utplot,ti,em,ut,title='emission measure',/ylog,ytitle='1e49 cm-3',charsize=1.5

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;background


IDL> print,anytim(anytim(ut,/second)+ti[700],/ccsds)
2012-07-05T03:23:53.689
IDL> print,anytim(anytim(ut,/second)+ti[750],/ccsds)
2012-07-05T03:25:36.089


IDL> print,mean(da[700:750,0])
  5.42263e-06
IDL> print,mean(da[700:750,1])
  2.87472e-07


lbk=mean(da[730:750,0])
hbk=mean(da[730:750,1])

da2=da
da2[*,0]=da[*,0]-lbk
da2[*,1]=da[*,1]-hbk


goes_tem,yclean=da2,tempr=temp,emis=em



window,1
utplot,ti,temp,ut,title='temperature',ytitle='1e6 Kelvin',charsize=1.5
plots,[ti[780],ti[780]],[0,40],color=cred
xyouts,ti[780],2,anytim(anytim(ut,/seconds)+ti[780],/ccsds),color=cred

write_png,'20120705temp.png',tvrd(true=1)


window,2
utplot,ti,em,ut,title='emission measure',/ylog,ytitle='1e49 cm-3',charsize=1.5
plots,[ti[780],ti[780]],[1e-4,1e1],color=cred
xyouts,ti[780],1e-3,anytim(anytim(ut,/seconds)+ti[780],/ccsds),color=cred



write_png,'20120705em.png',tvrd(true=1)


save,ti,ut,temp,em,da,file='20120705goes_tem.sav'







;chianti_dem --> [line, flux, eflux, hwhm, name]
;ラインの選択基準 --> EVLでラインブレンドのないライン
;hwhm --> 0.5A固定？