pro aia_loop_autodem,fileset,wave_,fov,ns_step,nw,segmin,vers,teem_table,chimax,test ;+ ; Project : AIA/SDO ; ; Name : AIA_LOOP_AUTODEM ; ; Category : Data analysis ; ; Explanation : Automated EM and Te analysis of coronal loops ; ; Syntax : IDL>aia_loop_autodem,fileset,wave_,fov,ns_step,nw,segmin,vers,teem_table ; ; Inputs : fileset = beginning of filename of FITS datafiles ; wave_[6] = wavelengths ; fov[4] = [x1,y1,x2,y2] field-of-view in solar radii ; ns_step = number of pixels to integrate along loop segment ; nw = loop cross-sectional width in pixels ; segmin = minimum loop length (in solar radii) to display ; vers = label in filename (e.g., vers='a') ; teem_table = DEM loopup table generated with AIA_TEEM_TABLE.PRO ; chimax = maximum chi-square limit ; test = test display (1) or none (0) ; ; Outputs : savefile *_teem_cube.sav ; ; History : 09-Mar-2011, Version 1 written by Markus J. Aschwanden ; 12-May-2011, check level nr, exposure normalization (Ding Yuan) ; 13-May-2011, goodness-of-fit (6 fluxes and cross-sections) ; 13-May-2011, chi-square with Poisson statistics ; 18-May-2011, change units of FOV from pixels to solar radii ; 18-May-2011, change segmin from pixels to solar radii ; 11-Jun-2011, add VERS to specify filename ; ; Contact : aschwanden@lmsal.com ;- ;____________________LOOP SEGMENT COUNTER___________________________ io =0 ct =3 contr =3 nwave =n_elements(wave_) nseg =0 nloop =0 loopfiles=strarr(nwave) for iw=0,nwave-1 do begin searchstring=fileset+vers+'_'+wave_(iw)+'_auto.dat' files =file_search(searchstring,count=nfiles) if (nfiles ge 1) then begin loopfiles[iw]=files[0] readcol,loopfiles[iw],iloop,x,y,z,s_length nl =max(iloop)+1 nloop =nloop+nl for il=0,nl-1 do begin ind =where(il eq iloop) spline_p,x(ind),y(ind),xx,yy nn =n_elements(xx) nseg_loop=long(nn/ns_step) nseg =nseg+nseg_loop endfor print,wave_(iw),' Number of loops=',nl endif endfor cube =fltarr(ns_step,nw,nseg,nwave) loop_par=fltarr(nseg,10) ;____________________READ DATA_______________________________________ texp_ =fltarr(nwave) for iw=0,nwave-1 do begin searchfile=fileset+'*'+wave_(iw)+'.fits' datafiles=file_search(searchfile,count=nfiles) if (nfiles ge 1) then begin datafile =datafiles[0] read_sdo,datafile,index0,data0 aia_prep,index0,data0,index,data texp =index.exptime level =index.lvl_num ;check level 1.0 or 1.5 number cdelt1 =index.cdelt1 crpix1 =index.crpix1 crpix2 =index.crpix2 rsun =index.rsun_obs rpix =rsun/cdelt1 x1 =fov[0] x2 =fov[2] y1 =fov[1] y2 =fov[3] i1 =long(x1*rpix+crpix1+0.5) j1 =long(y1*rpix+crpix2+0.5) if (iw eq 0) then begin i2 =long(x2*rpix+crpix1+0.5) j2 =long(y2*rpix+crpix2+0.5) nx =i2-i1+1 ny =j2-j1+1 if (iw eq 0) then images=fltarr(nx,ny,nwave) endif i2 =i1+nx-1 j2 =j1+ny-1 print,'FOV in pixels [i1,j1,i2,j1]=',i1,j1,i2,j2 if (level ge 1.5) then texp=1.0 ;already normalized texp_[iw]=texp images(*,*,iw)=data(i1:i2,j1:j2)/texp dpix =(fov(2)-fov(0))/float(nx) ;solar radii per pixel i1_fov =fov(0)/dpix ;pixels from Sun center j1_fov =fov(1)/dpix ;pixels from Sun center endif endfor ;for iw=0,nwave-1 ;_____________________READ DEM LOOKUP TABLE__________________________ restore,teem_table ;-->wave_,q94,area,resp_corr,telog,dte,tsig,flux dim =size(flux) nte =dim(1) nsig =dim(2) if (dim(3) ne nwave) then STOP,'DIM(FLUX) NE NWAVE' ;_____________________LOOP SEGMENT___________________________________ iseg =0 il =0 w_ =findgen(nw) w0 =w_(nw/2) nsm =nw/6 for iw=0,nwave-1 do begin ;loop detections in set of 6 filters if (loopfiles[iw] eq '') then goto,skip_read readcol,loopfiles[iw],iloop,x,y,z,s_length nloop =max(iloop)+1 for il=0,nloop-1 do begin ind =where(il eq iloop,ns) length=max(s_length(ind)) spline_p,x(ind),y(ind),xx,yy nn =n_elements(xx) ;number of points per loop nseg =long(nn/ns_step) ;number of loop segments for is=0,nseg-1 do begin is1 =is*ns_step is2 =is1+ns_step-1 i_loop=xx(is1:is2)/dpix+crpix1-i1 j_loop=yy(is1:is2)/dpix+crpix2-j1 for iwave=0,nwave-1 do begin array_curve,images(*,*,iwave),i_loop,j_loop,nw,xw,yw,zw cube(*,*,iseg,iwave)=float(zw) endfor loop_par(iseg,0)=iw ;wavelength set of loop detection loop_par(iseg,1)=il ;loop number if (is eq 0) then iseg1=iseg if (is eq nseg-1) then iseg2=iseg iseg=iseg+1 endfor ;for is=0,nseg-1 ;_____________________DISPLAY__________________________________________ idet =long(loop_par(iseg1,0)) lnr =long(loop_par(iseg1,1)) ;loop number print,'Chan#='+string(idet,'(I1)')+', Loop#'+string(lnr,'(I5)')+$ ', Segm #'+string(iseg,'(I5)')+', Length='+string(length,'(f8.5)') POSTSCRIPT: form =1 ;0=landscape, 1=portrait char =0.8 ;character size char2 =1.2 ;character size string0,5,lnr,fignr plotname='aia_loop_'+wave_(idet)+'_' ref ='' ;label at bottom of Fig. unit =0 ;window number if (length ge segmin) then begin fig_open,io,form,char,fignr,plotname,unit loadct,ct dx =0.45/nwave dy =0.45/(nseg>10) endif flux_obs=fltarr(nw,nwave,nseg) noise =fltarr(nw,nwave,nseg) fwhm =fltarr(nseg) valid =fltarr(nseg) cube3 =fltarr(nw,ns_step,nseg,nwave) cube3min=fltarr(nwave) cube3max=fltarr(nwave) pix_cm=0.6*0.716e8 nfree =3 ib1 =nw/6 ib2 =(nw*5)/6 ;__________________________DISPLAY LOOP PROFILE_____________________ gauss_norm=fltarr(nw,nseg) ib_ =fltarr(2,nseg) for is=iseg1,iseg2 do begin is_ =is-iseg1 prof =fltarr(nw) for i=0,ns_step-1 do prof=prof+reform(cube(i,*,is,idet))/float(ns_step) prof_smo=smooth(prof,3) if (max(prof) le 1) then prof_smo=smooth(prof_smo,3) b1 =avg(prof_smo(0:ib1)) b2 =avg(prof_smo(ib2:nw-1)) back =b1+(b2-b1)*(findgen(nw)-ib1)/(ib2-ib1) diff =(prof-back)>0 c0 =max(diff) c1 =w_(nw/2) c2 =(total(diff)/c0)/2.35482 coeff0=[c0,c1,c2] gauss =gaussfit(w_,diff,coeff,estimates=coeff0,nterms=nfree) peak =coeff(0)>0. w1 =coeff(1) gwid =abs(coeff(2)) wid =gwid*2.35482 ;FWHM gauss_norm(*,is_)=exp(-(w_-w1)^2/(2.*gwid^2)) > 0 gauss =peak*gauss_norm(*,is_) flux_obs(*,idet,is_)=diff fwhm(is_)=wid*pix_cm ib3 =(long(w1-gwid-0.5)>0 )<(nw-1) &ib_(0,is_)=ib3 ib4 =(long(w1+gwid+0.5)>ib3)<(nw-1) &ib_(1,is_)=ib4 fint0=total(diff(ib3:ib4)) > 0.1 if (peak le 0. ) then valid(is_)=-1 if (w1 lt ib1 ) or (w1 gt ib2 ) then valid(is_)=-1 if (wid le 1. ) or (wid ge nw-1) then valid(is_)=-1 for iwave=0,nwave-1 do begin prof=fltarr(nw) for i=0,ns_step-1 do prof=prof+reform(cube(i,*,is,iwave))/float(ns_step) prof_smo=smooth(prof,3) if (max(prof) le 1) then prof_smo=smooth(prof_smo,3) b1 =avg(prof_smo( 0:ib1)) b2 =avg(prof_smo(ib2:nw-1)) back=b1+(b2-b1)*(findgen(nw)-ib1)/(ib2-ib1) diff =(prof-back) > 0. fint =total(diff(ib3:ib4)) > 0.1 peak =max(flux_obs(*,idet,is_))*(fint/fint0) gauss=peak*gauss_norm(*,is_) flux_obs(*,iwave,is_)=diff noise(*,iwave,is_)=sqrt(prof)/sqrt(ns_step*texp_[iwave]) for i=0,ns_step-1 do cube3(*,i,is_,iwave)=reform(cube(i,*,is,iwave))$ -back/ns_step cube3min(iwave)=min(cube3(*,*,*,iwave)) cube3max(iwave)=max(cube3(*,*,*,iwave)) if (length ge segmin) then begin x1_ =0.05+iwave*dx x2_ =x1_+dx y1_ =0.51+(is_)*dy y2_ =y1_+dy y1 =min(prof) y2 =max(prof) dyy =y2/10. !p.position=[x1_,y1_,x2_,y2_] !y.range=[y1-dyy,y2+dyy*1.05] !x.range=[0,nw+1] !x.style=5 !y.style=5 plot,w_,prof,psym=10 oplot,w_,back,psym=10 if (valid(is_) eq 0) then oplot,w_,back+gauss,color=200,thick=2 oplot,[0,nw,nw,0,0],[y1-dyy,y1-dyy,y2+dyy,y2+dyy,y1-dyy] !noeras=1 endif endfor ;for iwave=0,nwave-1 endfor ;for is=iseg1,iseg2 ;__________________________DISPLAY LOOP CUBE_____________________ for is=iseg1,iseg2 do begin is_ =is-iseg1 for iwave=0,nwave-1 do begin if (length ge segmin) then begin x1_ =0.05+iwave*dx x2_ =x1_+dx y1_ =0.05+(is_)*dy y2_ =y1_+dy !p.position=[x1_,y1_,x2_,y2_] !p.title=' ' !x.range=[0,nw+1] !y.range=[0,ns_step] !x.style=5 !y.style=5 image =cube3(*,*,is_,iwave) if (io eq 0) then begin nxw =long(!d.x_vsize*(x2_-x1_)+0.5) nyw =long(!d.y_vsize*(y2_-y1_)+0.5) z =congrid(image,nxw,nyw) endif if (io ne 0) then z=image z2 =cube3max[iwave]>2. z1 =cube3min[iwave]<(-2.) plot,[0,0],[0,0] tv,bytscl(z,min=z1,max=z2),x1_,y1_,xsize=(x2_-x1_),ysize=(y2_-y1_),/normal loopnr =' - Loop #'+wave_(idet)+'/'+string(lnr,'(i4)') if (is eq iseg1) then begin xyouts,x1_,y1_-0.02,wave_(iwave)+' A',size=char2,/normal if (iwave eq idet) then xyouts,x1_,y1_-0.035,'DETECT',size=char,/normal if (iwave eq 0) then xyouts,0.05,0.97,fileset+loopnr,size=char2,/normal endif !noeras=1 endif endfor ;for iwave=0,nwave-1 endfor ;for is=iseg1,iseg2 ;__________________________PARTIAL IMAGE_______________________________ if (length ge segmin) then begin imin =min(xx)/dpix+crpix1-i1 ;pixels of partial image imax =max(xx)/dpix+crpix1-i1 jmin =min(yy)/dpix+crpix2-j1 jmax =max(yy)/dpix+crpix2-j1 range=((imax-imin)>(jmax-jmin))+20 x0 =(imin+imax)/2 y0 =(jmin+jmax)/2 ii1 =long((x0-range/2))>0 ii2 =long((x0+range/2))<(nx-1) jj1 =long((y0-range/2))>0 jj2 =long((y0+range/2))<(ny-1) image=images(ii1:ii2,jj1:jj2,iw) i_loop=xx/dpix+crpix1-i1 j_loop=yy/dpix+crpix2-j1 array_curve,images(*,*,iw),i_loop,j_loop,nw,i_xw,j_yw,zw xw =(i_xw+i1-crpix1)*dpix yw =(j_yw+j1-crpix2)*dpix x1_ =0.55 x2_ =0.95 y1_ =0.67 y2_ =0.97 !p.title=wave_(iw)+' A' !x.range=([ii1,ii2]+i1-crpix1)*dpix !y.range=([jj1,jj2]+j1-crpix2)*dpix !x.style=1 !y.style=1 !p.position=[x1_,y1_,x2_,y2_] if (io eq 0) then begin nxw =long(!d.x_vsize*(x2_-x1_)+0.5) nyw =long(!d.y_vsize*(y2_-y1_)+0.5) z =congrid(image,nxw,nyw) endif if (io ne 0) then z=image statistic,image,z0,dz c1=z0-contr*dz c2=z0+contr*dz plot,[0,0],[0,0] tv,bytscl(z,min=c1,max=c2),x1_,y1_,xsize=(x2_-x1_),ysize=(y2_-y1_),/normal oplot,xw(*,0) ,yw(*,0) ,color=255,thick=3 oplot,xw(*,nw/2),yw(*,nw/2),color=255,thick=1,linestyle=2 oplot,xw(*,nw-1),yw(*,nw-1),color=255,thick=3 oplot,xw(0,*) ,yw(0,*) ,color=255,thick=3 oplot,xw(nn-1,*),yw(nn-1,*),color=255,thick=3 endif ;____________________________DEM FORWARD-FIT___________________________ for is=iseg1,iseg2 do begin is_ =is-iseg1 chimin=9999. sigbest=chimin i3 =ib_(0,is_) i4 =ib_(1,is_) n34 =i4-i3+1 np =(n34*nwave-nfree) for k=0,nte-1 do begin for l=0,nsig-1 do begin peak_dem1 =reform(flux(k,l,*)) flux_dem1 =fltarr(n34,nwave) for iwa=0,nwave-1 do flux_dem1(*,iwa)=peak_dem1(iwa)*gauss_norm(i3:i4,is_) em1 =total(flux_obs(i3:i4,*,is_))/total(flux_dem1(*,*)) flux_mod =flux_dem1*em1 flux_diff =(flux_obs(i3:i4,*,is_)-flux_mod(*,*))>0. flux_sigm =noise(i3:i4,*,is_) chi =sqrt(total(flux_diff^2/flux_sigm^2)/np) if (chi le chimin) then begin chimin =chi em_best =alog10(em1) te_best =telog(k) sig_best =tsig(l) f_mod =flux_mod(*,*) f_obs =flux_obs(i3:i4,*,is_) f_noi =noise(i3:i4,*,is_) endif endfor endfor w_cm =fwhm(is_) log_te1=te_best-sig_best*(2.35482/2.) log_te2=te_best+sig_best*(2.35482/2.) em_tot=10.^em_best*(10.^log_te2-10.^log_te1) ;cm-5 ne_cm3=sqrt(em_tot/w_cm) loop_par(is,2)=em_best ;log(EM) loop_par(is,3)=te_best ;log(T) loop_par(is,4)=sig_best ;log(sig_T) loop_par(is,5)=chimin ;chi2 loop_par(is,6)=w_cm < 1.e9 ;wid[cm] loop_par(is,7)=ne_cm3 ;ne[cm-3] if (test eq 1) then begin window,1,xsize=768,ysize=768 n =n34*nwave x =findgen(n) ym =max(f_obs)*2 !p.position=[0.1,0.1,0.9,0.9] plot,x,f_obs,xstyle=0,ystyle=2,xrange=[0,n],yrange=[-ym,ym],psym=10 oplot,x,x*0.,linestyle=2 oplot,x,f_mod,thick=3 errplot,x,f_obs-f_noi,f_obs+f_noi xyouts_norm,0.1,0.9,'chi='+string(chimin,'(f8.2)'),char*2 read,'continue?',yes if (yes eq 0) then stop window,unit endif endfor ;___________________DISPLAY BEST-FIT PARAMETERS________________________ if (length ge segmin) then begin log_te =loop_par(iseg1:iseg2,3) log_sig=loop_par(iseg1:iseg2,4) log_ne =alog10(loop_par(iseg1:iseg2,7)) wid_mm =loop_par(iseg1:iseg2,6)/1.e8 chi =loop_par(iseg1:iseg2,5) s =findgen(nseg)*ns_step*0.6*0.916 ;loop length s[Mm] x1_ =0.60 x2_ =0.95 y1_ =0.50 y2_ =0.63 !p.title='' !x.range=minmax(s)+1.0*[-1,+1] !y.range=[5.7,8.0] !x.style=1 !y.style=0 !p.position=[x1_,y1_,x2_,y2_] ind=where(valid eq 0,nind) if (nind le 1) then goto,skip_plot plot,s(ind),log_te(ind),thick=2,psym=-4,ytitle='Temperature log(T!De!N)' xgrey=[s(ind),reverse(s(ind))] ygrey=[log_te(ind)+log_sig(ind),reverse(log_te(ind)-log_sig(ind))] polyfill,xgrey,ygrey,color=200 oplot,s(ind),log_te(ind),thick=2,psym=-4 for is=0,nind-1 do oplot,s(ind(is))*[1,1],log_te(ind(is))+log_sig(ind(is))*[-1,1] for j=6,7 do oplot,!x.range,[j,j],linestyle=1 statistic,log_te(ind),n_avg,n_sig xyouts_norm,0.1,0.82,'log(T!De!N)='+string(n_avg,'(f5.2)')+'+',char2 xyouts_norm,0.1,0.82,'log(T!De!N)='+string(n_avg,'(f5.2)')+'_'+$ string(n_sig,'(f5.2)'),char2 statistic,log_sig(ind),s_avg,s_sig xyouts_norm,0.1,0.64,'log(!9s!3!DT!N)='+string(s_avg,'(f5.2)')+'+',char2 xyouts_norm,0.1,0.64,'log(!9s!3!DT!N)='+string(s_avg,'(f5.2)')+'_'+$ string(s_sig,'(f5.2)'),char2 y1_ =0.35 y2_ =0.48 !p.position=[x1_,y1_,x2_,y2_] !y.range=[8>min(log_ne(ind)), (11