pro euv_dimming_calc,input,dir,catfile,run,iflare,test,dem,ifit2 ;+ ; 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,input,catfile,dir,iflare ; ; Inputs : dir,savefile ; ; Outputs ; IDL save file: *.sav ; ; History : 24-Feb-2016, Version 1 written by Markus J. Aschwanden ; ; Contact : aschwanden@lmsal.com ;- common dimming,model,t_hrs,a_t,v_t,x_t,q_t,q_obs,h0,ind_fit print,'__________________EUV_DIMMING_CALC3_________________________ readcol,dir+catfile,event_,dateobs_,noaa_,helpos_,instr_,$ w1_,w2_,w3_,w4_,w5_,w6_,w7_,w8_,fov_,$ format='(I,A,A,A,A,A,A,A,A,A,A,A,A,A)',/silent nt =n_elements(event_) i =0 wave8 =[w1_(i),w2_(i),w3_(i),w4_(i),w5_(i),w6_(i),w7_(i),w8_(i)] indw =where(wave8 ne 0,nwave) wave_ =strtrim(string(wave8(indw),'(I5)'),2) tmargin1=input.tmargin1 ;s tmargin2=input.tmargin2 rise =input.rise decay =input.decay duration=rise+decay cadence =input.cadence it_start=long(tmargin1/cadence) it_peak =long((tmargin1+rise)/cadence) it_end =long((tmargin1+rise+decay)/cadence) helpos =input.helpos lon =long(strmid(helpos,4,2)) lat =long(strmid(helpos,1,2)) l =lon*!pi/180. b =lat*!pi/180. angle =sqrt(l^2+b^2) r0 =6.96e10 ;cm d0 =r0*sin(angle) ;distance of CME lauch from disk center ;________________________AIA RESPONSE FUNCTION____________________ teem_table='teem_table.sav' te_range=[0.5,20]*1.e6 ;[K] valid temperature range of DEM tsig =0.1*(1+findgen(10)) ;values of Gaussian log temp widths tresp =aia_get_response(/temp,/full,/dn) telog_ =tresp.logte telog1 =alog10(te_range(0)) telog2 =alog10(te_range(1)) ind_te =where((telog_ ge telog1) and (telog_ le telog2),nte) telog =telog_(ind_te) ichan_ =fltarr(nwave) resp =fltarr(nte,nwave) for iw=0,nwave-1 do begin filter ='A'+wave_(iw) if (wave_(iw) eq '094') or (wave_(iw) eq '94') then filter='A94' ichan =where(tresp.channels eq filter) resp_ =tresp.tresp(*,ichan) resp(*,iw)=resp_(ind_te) endfor ;_____________________CALCULATES LOOPUP TABLES____________________ if (dem eq 1) then begin dte1 =10.^telog(1:nte-1)-10.^telog(0:nte-2) dte =[dte1(0),dte1] em1 =1. nsig =n_elements(tsig) flux =fltarr(nte,nsig,nwave) for i=0,nte-1 do begin for j=0,nsig-1 do begin em_te =em1*exp(-(telog-telog(i))^2/(2.*tsig(j)^2)) for iw=0,nwave-1 do flux(i,j,iw)=total(resp(*,iw)*em_te*dte) endfor endfor save,filename=teem_table,wave_,resp,telog,dte,tsig,flux print,'Lookup table created : ',teem_table endif restore,teem_table ;-->wave_,resp,telog,dte,tsig,flux ;__________________FLUX GRID DATA_________________________ time =fltarr(nt) flux_prof =fltarr(nt,nwave) emtot_prof=fltarr(nt) telog_prof=fltarr(nt) tsig_prof =fltarr(nt) chi_prof =fltarr(nt) wid_prof =fltarr(nt) for it=0,nt-1 do begin print,'time step = ',it,nt-1 nr_it =it+1 ind =where((event_) eq nr_it,nind) if (nind le 0) then stop,'NLFFF_EVENT : NO MATCHING EVENT FOUND' i =ind(0) ;selected entry string0,3,nr_it,nr_str hh =strmid(dateobs_(i), 9,2) mm =strmid(dateobs_(i),11,2) ss =strmid(dateobs_(i),13,2) time(it)=hh+mm/60.+ss/3600. if (it ge 1) then $ ;midnight jump if (time(it) lt time(it-1)) then time(it)=time(it)+24. savefile =dateobs_(i)+'_'+instr_(i)+'_'+nr_str+'_'+run+'.sav' ;__________________FLUX PROFILES____________________________________ restore,dir+savefile ;--> PARA,FLUX_GRID(NX,NY,NWAVE) for iwave=0,nwave-1 do flux_prof(it,iwave)=avg(flux_grid(*,*,iwave)) if (it eq 0) then begin dim =size(flux_grid) nx =dim(1) ny =dim(2) em_xy =fltarr(nx,ny,nt) telog_xy=fltarr(nx,ny,nt) tsig_xy =fltarr(nx,ny,nt) chi2 =fltarr(nx,ny,nt) endif ;________________DEM ANALYSIS________________________________ dim =size(flux) ;precalculated TE+EM loopup table nte =dim[1] nsig =dim[2] nwave =dim[3] nfree =3 for i=0,nx-1 do begin ;x,y grid for j=0,ny-1 do begin em_best =0. telog_best=0. tsig_best=0. flux_obs=reform(flux_grid(i,j,*)) ;flux per second counts=flux_obs*exptime_ noise =sqrt(counts/exptime_) chi_best=9999. for k=0,nte-1 do begin for l=0,nsig-1 do begin flux_dem1=reform(flux(k,l,*)) em1 =total(flux_obs)/total(flux_dem1) flux_dem=flux_dem1*em1 ind =where(flux_obs ne 0.,nind) chi =sqrt(total((flux_obs(ind)-flux_dem(ind))^2/noise(ind)^2)/(nind-nfree)) if (chi le chi_best) then begin em_best =em1 telog_best =telog(k) tsig_best =tsig(l) chi_best =chi endif endfor endfor em_xy(i,j,it) =em_best telog_xy(i,j,it)=telog_best tsig_xy(i,j,it) =tsig_best chi2(i,j,it) =chi_best endfor ;x location loop endfor ;y location loop endfor ;time step loop ;_________________SPATIAL SYNTHESIS OF DEM____________________________________ dem_t =fltarr(nte,nt) dem_grid=fltarr(nx,ny,nt) ngrid =float(nx*ny) chi_med =fltarr(nt) for it=0,nt-1 do begin for j=0,ny-1 do begin for i=0,nx-1 do begin em_best =em_xy(i,j,it) telog_best =telog_xy(i,j,it) tsig_best =tsig_xy(i,j,it) em_kelvin =em_best*exp(-(telog-telog_best)^2/(2.*tsig_best^2))/ngrid dem_grid(i,j,it)=total(em_kelvin*dte) ;DEM localized x,y dem_t(*,it)=dem_t(*,it)+em_kelvin ;DEM summed x,y endfor endfor emtot_prof(it)=total(dem_t(*,it)*dte) ;EM [cm-5] dem_max=max(dem_t(*,it),imax) ;EM [cm-5 K-1] telog_prof(it)=telog(imax) ;T[MK] chi_med(it)=median(chi2(*,*,it)) endfor ;__________________ELIMINATE WORST DEM FITS_____________________ chi_bad =5.0 ind_bad =where(chi_med ge chi_bad,nbad) if (nbad ge 1) then begin for i=0,nbad-1 do begin it =ind_bad(i) it1 =(it-1) &if (it eq 0 ) then it1=1 it2 =(it+1) &if (it eq nt-1) then it2=nt-2 emtot_prof(it) =(emtot_prof(it1)+emtot_prof(it2))/2. dem_grid(*,*,it)=(dem_grid(*,*,it1)+dem_grid(*,*,it2))/2. dem_t(*,it)=(dem_t(*,it1)+dem_t(*,it2))/2. telog_xy(*,*,it)=(telog_xy(*,*,it1)+telog_xy(*,*,it2))/2. endfor print,'Eliminate time with bad fit it_bad=',it endif te_mk =10.^(median(telog_xy(*,*,0)))/1.e6 h0 =4.7e9*te_mk ;[cm] density scale height ;_________________SELECT DIMMING LOCATIONS______________________ npeak =fltarr(nt) for i=0,nx-1 do begin for j=0,ny-1 do begin em_peak=max(dem_grid(i,j,*),itp) npeak(itp)=npeak(itp)+1 endfor endfor peak =max(npeak(1:it_end),itpeak1) &itpeak=itpeak1+1 window,1,xsize=512,ysize=512 clearplot plot,time,npeak,psym=10 oplot,time(itpeak)*[1,1],!y.crange,thick=2 ;_________________TOTAL EMISSION DECREASES EVOLUTION____________ itol =1 em_decr=fltarr(nt) for i=0,nx-1 do begin for j=0,ny-1 do begin em_peak=max(dem_grid(i,j,*),itp) dem_model=fltarr(nt) if (itp ge itpeak-itol) and (itp le itpeak+itol) then begin for it=nt-2,0,-1 do begin diff=reform(dem_grid(i,j,it)-dem_grid(i,j,it+1)) > 0. dem_model(it)=dem_model(it+1)+diff endfor endif em_decr=em_decr+dem_model endfor endfor q_obs =em_decr/max(em_decr) em_max =max(emtot_prof) em_obs =em_max*q_obs ;_________________CME VELOCITY AND ACCELERATION___________ n_model=4 a_model =fltarr(nt,n_model) v_model =fltarr(nt,n_model) x_model =fltarr(nt,n_model) em_model =fltarr(nt,n_model) a0_model =fltarr(n_model) t0_model =fltarr(n_model) ta_model =fltarr(n_model) qe_model =fltarr(n_model) ekin_model=fltarr(n_model) chi_model =fltarr(n_model) x_lasco =fltarr(n_model) v_lasco =fltarr(n_model) t_lasco =fltarr(n_model) t_hrs =time ;hrs dt_hrs =time(1)-time(0) dt =dt_hrs*3600. t0_hrs =time(imax) ;s it =findgen(nt) v0 =3.0e7 ;cm/s for model=1,4 do begin m =model-1 tau_hrs=0.1 ;hrs a0 =v0/(dt_hrs*3600.) qmin =0.01 coeff =[a0/1.e5,tau_hrs,t0_hrs,qmin] ncoeff =n_elements(coeff) xi =fltarr(ncoeff,ncoeff) for i=0,ncoeff-1 do xi(i,i)=1.0 ftol =1.e-4 nfit =ifit2(1)-ifit2(0)+1 ind_fit=ifit2(0)+findgen(nfit) powell,coeff,xi,ftol,fmin,'euv_dimming_powell' chi =euv_dimming_powell(coeff) ;-->a_model,v_model,x_model,q_t print,coeff,model,chi acc5 =coeff(0) tau_hrs =coeff(1) t0_hrs =coeff(2) ;hrs qmin =coeff(3) a_model(*,m)=a_t v_model(*,m)=v_t x_model(*,m)=x_t+d0 em_model(*,m)=em_max*q_t a0_model(m)=acc5*1.e5 t0_model(m)=t0_hrs*3600. ta_model(m)=(t0_hrs+tau_hrs)*3600. ;s chi_model(m)=chi endfor ;__________________LASCO/C2 TIME AND VELOCITY_____________________ readcol,'lasco.dat',flare,dateobs_lasco,h_lasco,v1_lasco,v2_lasco,v3_lasco,$ m_lasco,e_lasco,c,skipline=1,format='(I,A,F,I,I,I,F,F,A)' ind0 =where(flare eq iflare,n) i =ind0(0) d_lasco =h_lasco(i) x_lasco =r0*h_lasco(i) ;cm for m=0,n_model-1 do begin ind1 =where(x_model(*,m) ge x_lasco,n_lasco) if (max(x_model(*,m)) lt x_lasco) then ind1=[nt-1] v_lasco(m)=v_model(ind1(0),m) ;cm/s t_lasco(m)=time(ind1(0))*3600. ;s endfor ;_________________DIMMING AREA_________________________________ emax =max(emtot_prof,it_max) emin =min(emtot_prof(it_max:nt-1),dit_min) em_tot =emax-emin e1_xy =em_xy(*,*,it_max) e2_xy =em_xy(*,*,nt-1) de_med =median(e1_xy-e2_xy) e3_xy =e1_xy-de_med*2. ind_dimm=where(e2_xy lt e3_xy,npix_dimm) q_dimm =float(npix_dimm)/float(nx*ny) wid_dimm=sqrt(float(npix_dimm)) ;width in units of grid bins l_cm =wid_dimm*dx_grid_cm ;cm area =l_cm^2 ;equivalent area [cm>2] vol =l_cm^3 ;volume [cm3] angle =90.*(l_cm/(r0*!pi/2.)) ;angular width [deg] ;_________________DIMMING MASS AND KINETIC ENERGY________________ nel =sqrt(em_tot/l_cm) mp =1.67e-24 ;proton mass mass =mp*nel*vol ekin_model=(1./2.)*mass*v_lasco^2 ;________________SAVE PARAMETERS______________________________ string0,3,iflare,nr_str savefile2='dimming_'+nr_str+run+'.sav' save,filename=dir+savefile2,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,angle,area,vol,em_tot,dem_t,nel,te_mk,telog,mass,time,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 print,'parameters saved in file = ',savefile2 end