pro euvi_loopdem,imagefiles,loopfile,iloop,backgr,w_pix,outfile ;+ ; Project : STEREO ; ; Name : EUVI_LOOPDEM ; ; Category : Data analysis ; ; Explanation : This routine extracts the flux F(NW,NS,NSC,NF) ; in curved subarrays ; from stereopair images saved in files IMAGEFILES, ; with NW=width of strip, NS=length of stripe, NF=3 ; filters (171, 195, 298), and NSC=2 spacecraft (A,B), ; with the coordinates X(NW,NS,NSC) and Y(NW,NS,NSC) ; defined by loop coordinates [x(s),y(s)] in loopfile ; under loop number index (ILOOP-1). A background ; B(NW,NS,NSC,NF) is evaluated according to method ; BACKGR=1,2,3,4. A DEM-fit is then performed to the ; background-subtracted fluxes, yielding a temperature ; profile T(s,2), temperature (gaussian DEM) width dT(s,2), ; and density profile n(s,2), with average width w. ; All the output parameters are stored in outfile ; ; Syntax : IDL>imagefiles,loopfile,iloop,backgr,nw,outfile ; ; Inputs : imagefiles ;containing [image_pair,para] for 171, 195, 284 ; loopfile ;file with loop coordinates ; iloop ;selected loop nr ; backgr ;background subtraction method ; nw ;width of extracted curved array ; ; Output outfile ;IDL savefile containing output parameters ; (XS,YS,F,B,S,T,DT,N,NS,W) ; ; History : Dec 2007, Version 1 written by Markus J. Aschwanden ; ; Contact : aschwanden@lmsal.com ;- common powell_para,f1,f2,f3,dfit,dx0,dxg nw =long(2*(2+w_pix)+1) &print,'nw=',nw ib1 =2 ib2 =nw-3 dfit =nw/2-ib1 dx0 =0.5 dxg =1. xw =findgen(nw) x0 =float(nw/2) nww =(ib2-ib1)*10 xww =ib1+0.1*findgen(nww) ;WAVELENGTH IMAGES_____________________________________________________ nfilter =n_elements(imagefiles) if (nfilter ne 3) then stop,'EUVI_LOOPDEM: input has not 3 imagefiles' for iw=0,2 do begin restore,imagefiles(iw) ;-->image_pair(2,nx,ny),para,index_a,index_b dim =size(image_pair) nx =dim(2) ny =dim(3) ;LOOP COORDINATES______________________________________________________ if (iw eq 0) then begin readcol,loopfile,iloop_,ix_,iy_,iz_,ih_,ih_raw,dh_,nsm_,poly_,chi2_,q1_,q2_ loop_index=iloop-1 ind =where(iloop_ eq loop_index,np) xp =ix_(ind) yp =iy_(ind) zp =iz_(ind) euvi_spline,para,xp,yp,zp,x,y,z,d,r,r_n euvi_stereo_rot,x,y,para,r_n,x_b,y_b,z_b ns =n_elements(x) endif ;LOOP LENGTH COORDINATE AND ALTITUDE_________________________________ if (iw eq 0) then begin s =fltarr(ns) h =fltarr(ns) for is=1,ns-1 do s(is)=s(is-1)+$ sqrt((x(is)-x(is-1))^2+(y(is)-y(is-1))^2+(z(is)-z(is-1))^2) cdelt1=para.cdelt1 s =s*cdelt1*0.725 ;Mm h =(r_n-1.)*696. ;Mm l =max(s) print,'loop length (3D) = ',l ,' Mm' s2 =fltarr(ns) for is=1,ns-1 do s2(is)=s2(is-1)+sqrt((x(is)-x(is-1))^2+(y(is)-y(is-1))^2) s2 =s2*cdelt1*0.725 ;Mm l2 =max(s2) print,'loop length (2D) = ',l2,' Mm' endif ;CURVED ARRAY_________________________________________________________ if (iw eq 0) then begin xs =fltarr(ns,nw,2) ys =fltarr(ns,nw,2) f =fltarr(ns,nw,2,3) b =fltarr(ns,nw,2,3) g =f endif imagea =reform(image_pair(1,*,*)) imageb =reform(image_pair(0,*,*)) array_curve,imagea,x ,y ,nw,xwa,ywa,fwa array_curve,imageb,x_b,y_b,nw,xwb,ywb,fwb xs(*,*,0)=xwa xs(*,*,1)=xwb ys(*,*,0)=ywa ys(*,*,1)=ywb f(*,*,0,iw)=fwa f(*,*,1,iw)=fwb ;BACKGROUND___________________________________________________________ backgr_ =['No background','Bright-Dark','Peak-Valley','Linear Background',$ 'Quadratic Background','Cubic Background','Cospatial Gaussian-Fit',$ 'Smoothed (3x3) Background','Smoothed (5x5) Background',$ 'Smoothed (7x7) Background'] method =backgr_(backgr+2) if (backgr ge 5) and (backgr le 7) then begin if (backgr eq 5) then nsm=3 if (backgr eq 6) then nsm=5 if (backgr eq 7) then nsm=7 imagea_s=smooth(imagea,nsm) imageb_s=smooth(imageb,nsm) array_curve,imagea_s,x ,y ,nw,xwa,ywa,bwa array_curve,imageb_s,x_b,y_b,nw,xwb,ywb,bwb b(*,*,0,iw)=bwa b(*,*,1,iw)=bwb endif for isc=0,1 do begin if (backgr eq -2) then begin ;no background if (isc eq 0) then b(*,*,isc,iw)=0. if (isc eq 1) then b(*,*,isc,iw)=0. endif if (backgr eq -1) then begin ;bright-dark method i1 =long(min(x)-nw)>0 i2 =long(max(x)+nw)0 j2 =long(max(y)+nw)',dev x0 =(coeff(0)>(nw/2-dx0))<(nw/2+dx0) xg =(coeff(1)>0.)0 y2 =coeff(3)>0 y3 =coeff(4)>0 b1 =coeff(5)>0 b2 =coeff(6)>0 b3 =coeff(7)>0 db1 =coeff(8) db2 =coeff(9) db3 =coeff(10) coeff =[x0,xg,y1,y2,y3,b1,b2,b3,db1,db2,db3] b(is,*,isc,0)=b1+db1*(xw-x0) b(is,*,isc,1)=b2+db2*(xw-x0) b(is,*,isc,2)=b3+db3*(xw-x0) g(is,*,isc,0)=b1+db1*(xww-x0)+y1*exp(-(xww-x0)^2/(2.*xg^2)) g(is,*,isc,1)=b2+db2*(xww-x0)+y2*exp(-(xww-x0)^2/(2.*xg^2)) g(is,*,isc,2)=b3+db3*(xww-x0)+y3*exp(-(xww-x0)^2/(2.*xg^2)) endfor endfor endif b =b < f ;no negative flux ;EUVI FILTER RATIOS AND DEM-FIT_________________________________________ dir ='/net/socrates/Users/aschwand/idl_stereo/' respfile=concat_dir('$SSW_MJASTEREO','idl/euvi_response.dat') readcol,respfile,t_mk,r1_a,r2_a,r3_a,r4_a,r1_b,r2_b,r3_b,r4_b nt =n_elements(t_mk) dt_mk =shift(t_mk,-1)-t_mk dt_mk(nt-1)=dt_mk(nt-2) filterfile=dir+'euvi_filterratio.dat' readcol,filterfile,t0,sigt,q1a,q2a,q1b,q2b t =fltarr(ns,2) dt =fltarr(ns,2) n =fltarr(ns,2) flux =fltarr(ns,2,3) frac =fltarr(ns,2,3) wid =fltarr(ns,2,3) tot =fltarr(ns,2,3) tfilter =[1.0,1.5,2.0] pix8 =cdelt1*0.725 w8 =fltarr(2) for isc=0,1 do begin if (isc eq 0) then begin ;STEREO A q1 =q1a q2 =q2a r1 =r1_a r2 =r2_a r3 =r3_a endif if (isc eq 1) then begin ;STEREO B q1 =q1b q2 =q2b r1 =r1_b r2 =r2_b r3 =r3_b endif for is=0,ns-1 do begin floop1 =f(is,ib1:ib2,isc,0)-b(is,ib1:ib2,isc,0) floop2 =f(is,ib1:ib2,isc,1)-b(is,ib1:ib2,isc,1) floop3 =f(is,ib1:ib2,isc,2)-b(is,ib1:ib2,isc,2) tot(is,isc,0)=total(floop1) ;total loop flux 171 A (s=is) tot(is,isc,1)=total(floop2) ;total loop flux 195 A (s=is) tot(is,isc,2)=total(floop3) ;total loop flux 284 A (s=is) if (max(floop1) gt 0) then wid(is,isc,0)=total(floop1)/max(floop1) ;equiv width if (max(floop2) gt 0) then wid(is,isc,1)=total(floop2)/max(floop2) ;equiv width if (max(floop3) gt 0) then wid(is,isc,2)=total(floop3)/max(floop3) ;equiv width endfor tot1 =total(tot(*,isc,0)) tot2 =total(tot(*,isc,1)) tot3 =total(tot(*,isc,2)) ind1 =where(wid(*,isc,0) gt 0.) ind2 =where(wid(*,isc,1) gt 0.) ind3 =where(wid(*,isc,2) gt 0.) wid_pix=avg(wid(ind1,isc,0)) if (tot2 gt tot1) then wid_pix=avg(wid(ind2,isc,1)) if (tot3 gt tot1) and (tot3 gt tot2) then wid_pix=avg(wid(ind3,isc,2)) w8(isc)=wid_pix*pix8 for is=0,ns-1 do begin im =nw/2 b1 =b(is,nw/2,isc,0) b2 =b(is,nw/2,isc,1) b3 =b(is,nw/2,isc,2) f1 =f(is,nw/2,isc,0)-b1 f2 =f(is,nw/2,isc,1)-b2 f3 =f(is,nw/2,isc,2)-b3 frac(is,isc)=(f1+f2+f3)/(b1+b2+b3+f1+f2+f3) if (f1 gt 0) and (f2 gt 0) then begin q1_obs=f2/f1 q2_obs=f3/f2 dev =sqrt((q1_obs-q1)^2+(q2_obs-q2)^2) devmin =min(dev,icase) t(is,isc) =t0(icase) dt(is,isc)=sigt(icase) em_mod=exp(-(t_mk-t(is,isc))^2/(2.*dt(is,isc)^2)) f2_mod=total(em_mod*r2*dt_mk) ;r2=>195 A em44 =f2/f2_mod n(is,isc)=1.e10*sqrt(em44/(wid_pix*pix8^3)) ;cm-3 em =em44*exp(-(t_mk-t(is,isc))^2/(2.*dt(is,isc)^2)) flux(is,isc,0)=(total(em*r1*dt_mk)+b1) flux(is,isc,1)=(total(em*r2*dt_mk)+b2) flux(is,isc,2)=(total(em*r3*dt_mk)+b3) endif endfor endfor ;SAVE OUTPUT____________________________________________________________ save,filename=outfile,xs,ys,f,g,b,s,h,t,dt,n,ns,w8,wid,l,l2,flux,xww,ib1,ib2,nw,frac,method print,'output saved in file = ',outfile end