pro euv_dimming_data,dir,iev,vers,test ;+ ; Project : SDO/AIA, ; ; Name : EUV_DIMMING ; ; 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,dir,iev,vers,test ; ; Inputs : dir,iev,vers,test ; ; Outputs ; IDL save file: *.sav ; ; History : 24-Feb-2016, Version 1 written by Markus J. Aschwanden ; ; Contact : aschwanden@lmsal.com ;- print,'__________________EUV DIMMING_____________________________ string0,3,iev,iev_str catfile ='dimming_'+iev_str+vers+'.txt' 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 ;________________________TIME SEQUENCE____________________________ nt =n_elements(event_) &print,'Number of time steps nt= ',nt for i=0,nt-1 do begin print,'==============> TIME STEP ',i+1,'/',nt,' <================' it =i+1 string0,3,it,it_str savefile_t='dimming_'+iev_str+'_'+it_str+vers+'.sav' restore,dir+savefile_t ;--> INPUT,PARA dateobs =strmid(dateobs_(i),0,15) instr =input.instr fov0 =input.fov0 ;width of position angle range nbin =input.nbin helpos =input.helpos ;flare position heliographic coordinates wave_disp=4 ;193 A channel wave_ =para.wave_ x1_sun =para.x1_sun x2_sun =para.x2_sun y1_sun =para.y1_sun y2_sun =para.y2_sun fluxmax =para.fluxmax wave8 =input.wave8 indw =where(wave8 ne 0,nwave) wave_ =strtrim(string(wave8[indw],'(I5)'),2) nwave =n_elements(wave_) ;print,'__________TIME WINDOW______________________________________ cadence =12 ;[seconds] AIA cadence for 7 wavelength fitlers year =strmid(dateobs,0,4) month =strmid(dateobs,4,2) day =strmid(dateobs,6,2) hour =strmid(dateobs,9,2) minu =strmid(dateobs,11,2) seco =strmid(dateobs,13,2) t1_file =year+'-'+month+'-'+day+' '+hour+':'+minu+':'+seco t1_sec =anytim(t1_file,/sec) t2_sec =t1_sec+cadence ;12 s AIA cadence for 7 wavelength fitlers t1 =anytim(t1_sec,/hxrbs) t2 =anytim(t2_sec,/hxrbs) print,'Start time = ',t1 print,'End time = ',t2 ;READ AIA DATA_______________________________________________ if (instr eq 'AIA') then begin aia=sdo_time2files(t1,t2,/aia,waves=wave_,img_type='LIGHT',level=1.5) if (aia(0) eq '') then begin nwave =0 wave_ =['none'] print,'No corresponding AIA image found' goto,end_proc endif nwave =n_elements(aia) ;reduce to existing files wave_ =strarr(nwave) read_sdo,aia,index,data_wave ;-->data(4096,4096,7) help,index,data_wave for iwave=0,nwave-1 do begin wave=long(index[iwave].wavelnth) if (wave le 99) then ndigit=2 if (wave ge 100) then ndigit=3 if (wave ge 1000) then ndigit=4 string0,ndigit,wave,wave_str wave_(iwave)=wave_str endfor endif ;DISPLAY IMAGE______________________________________________ if (test eq 1) then begin for iwave=0,nwave-1 do begin print,'read file = ',aia(iwave),' wavelength=',wave_(iwave) window,0,xsize=1024,ysize=1024 loadct,0 statistic,data_wave(*,*,iwave),zavg,zsig tv,bytscl(rebin(data_wave(*,*,iwave),1024,1024),min=0,max=3*zsig) endfor endif ;________________WAVELENGTH LOOP_____________________________ nwave =n_elements(wave_) exptime_=fltarr(nwave) first_wave=0 for iwave=0,nwave-1 do begin wave =wave_[iwave] ;print,'__________EXTRACT FULL EUV IMAGE_______________________________ first_wave=first_wave+1 data =data_wave(*,*,iwave) dateobs=index(iwave).date_obs naxis1 =index(iwave).naxis1 naxis2 =index(iwave).naxis2 cdelt1 =index(iwave).cdelt1 cdelt2 =index(iwave).cdelt2 crpix1 =index(iwave).crpix1 &if (crpix1 gt naxis1) then crpix1=naxis1/2 crpix2 =index(iwave).crpix2 &if (crpix2 gt naxis2) then crpix2=naxis2/2 crval1 =index(iwave).crval1 crval2 =index(iwave).crval2 exptime=index(iwave).exptime exptime_(iwave)=exptime eph =get_sun(dateobs) rsun =eph(1) ;rsun=index.rsun_obs in arcseconds rpix =rsun/cdelt1 ;solar radius in number of pixels dpix_euv=1./rpix ;pixel size in units of solar radius i0_sun =long(crpix1-crval1/cdelt1) j0_sun =long(crpix2-crval2/cdelt2) if (first_wave eq 1) then begin nx =long(1+(x2_sun-x1_sun)/dpix_euv) ny =long(1+(y2_sun-y1_sun)/dpix_euv) endif image1 =fltarr(nx,ny) i1 =long(i0_sun+x1_sun/dpix_euv) > 0 i2 =(i1+nx) < (naxis1-1) j1 =long(j0_sun+y1_sun/dpix_euv) > 0 j2 =(j1+ny) < (naxis2-1) image1 =data(i1:i2,j1:j2) ;_________________DISPLAY EUV SCANS_______________________________ if (test ge 1) and (wave eq wave_disp) then begin nx_disp=1024 ny_disp=1024 window,0,xsize=nx_disp,ysize=ny_disp loadct,1 x1_ =0.0 x2_ =1.0 y1_ =0.0 y2_ =1.0 !p.position=[x1_,y1_,x2_,y2_] r0 =1.3 !x.range=[-1.,1.]*r0 !y.range=[-1.,1.]*r0 !x.style=1 !y.style=1 plot,[x1_sun,x2_sun,x2_sun,x1_sun,x1_sun],[y1_sun,y1_sun,y2_sun,y2_sun,y1_sun] oplot,[0,0],[-1.3,1.3],linestyle=2 oplot,[-1.3,1.3],[0,0],linestyle=2 nphi =361 phi =2.*!pi*findgen(nphi)/float(nphi-1) oplot,cos(phi),sin(phi) rsun =1 pos =0. crota2 =0. blong =0. blat =0. dlat =10. coordsphere,rsun,pos,crota2,blat,blong,dlat !noeras=1 x1_ =(r0+x1_sun)/(2.*r0) x2_ =(r0+x2_sun)/(2.*r0) y1_ =(r0+y1_sun)/(2.*r0) y2_ =(r0+y2_sun)/(2.*r0) !p.position=[x1_,y1_,x2_,y2_] !x.range=[x1_sun,x2_sun] !y.range=[y1_sun,y2_sun] loadct,3 plot,[x1_sun,x2_sun,x2_sun,x1_sun,x1_sun],[y1_sun,y1_sun,y2_sun,y2_sun,y1_sun] nxw =long(!d.x_vsize*(x2_-x1_)+0.5) nyw =long(!d.y_vsize*(y2_-y1_)+0.5) z =congrid(image1,nxw,nyw) ind0 =where(image1 ne 0) statistic,image1(ind0),zavg,zsig c1 =zavg-3*zsig c2 =zavg+3*zsig tv,bytscl(z,min=c1,max=c2),x1_,y1_,xsize=(x2_-x1_),ysize=(y2_-y1_),/normal endif ;print,'__________FLUX GRID_______________________________________ if (first_wave eq 1) then begin nxx =long(nx/nbin) nyy =long(ny/nbin) dpix_bin=dpix_euv*nbin x_grid=x1_sun+dpix_bin*findgen(nxx) y_grid=y1_sun+dpix_bin*findgen(nyy) flux_grid=fltarr(nxx,nyy,nwave) endif image2 =float(image1)/exptime ;normalized to exposure time of 1 s if (nbin eq 1) then flux_grid(*,*,iwave)=image2(0:nx-1,0:ny-1) if (nbin ge 2) then $ flux_grid(*,*,iwave)=rebin(image2(0:nxx*nbin-1,0:nyy*nbin-1),nxx,nyy) END_WAVE: endfor ;________________SAVE PARAMETERS______________________________ para.dpix_euv=dpix_euv ;pixel size in units of solar radius rsun_cm =6.96e10 ;cm dx_grid_cm=nbin*dpix_euv*rsun_cm ;grid pixel size [cm] save,filename=dir+savefile_t,input,para,exptime_,flux_grid,x_grid,y_grid,dx_grid_cm print ,'parameters saved in file = ',savefile_t endfor END_PROC: end