;+ ; ; NAME:wavpol ; ; MODIFICATION HISTORY:Written By Chris Chaston, 30-10-96 ; :Modified by Vassilis, 2001-07-11 ; :Modified by Olivier Le Contel, 2008-07 ; to be able to change nopfft and steplength ;PURPOSE:To perform polarisation analysis of three orthogonal component time ; series data. ; ;EXAMPLE: wavpol,ct,Bx,By,Bz,timeline,freqline,powspec,degpol,waveangle,elliptict,helict ; ;CALLING SEQUENCE: wavpol,ct,Bx,By,Bz,timeline,freqline,powspec,degpol,waveangle,elliptict,helict ; ;INPUTS:ct,Bx,By,Bz, are IDL arrays of the time series data; ct is cline time ; ; Subroutine assumes data are in righthanded fieldaligned ; coordinate system with Z pointing the direction ; of the ambient magnetic field. ; ; threshold:-if this keyword is set then results for ellipticity, ; helicity and wavenormal are set to Nan if below 0.6 deg pol ; ;;Keywords: ; nopfft(optional) = Number of points in FFT ; ; steplength(optional) = The amount of overlap between successive FFT intervals ; ;OUTPUTS: The program outputs five spectral results derived from the ; fourier transform of the covariance matrix (spectral matrix) ; These are follows: ; ; Wave power: On a linear scale, at this stage no units ; ; Degree of Polarisation: ; This is similar to a measure of coherency between the input ; signals, however unlike coherency it is invariant under ; coordinate transformation and can detect pure state waves ; which may exist in one channel only.100% indicates a pure ; state wave. Less than 70% indicates noise. For more ; information see J. C. Samson and J. V. Olson 'Some comments ; on the description of the polarization states ; of waves' Geophys. J. R. Astr. Soc. (1980) v61 115-130 ; ; Wavenormal Angle: ; the angle between the direction of minimum ; variance calculated from the complex off diagonal ; elements of the spectral matrix and the Z direction ; of the input ; ac field data. For magnetic field data in ; field aligned coordinates this is the ; wavenormal angle assuming a plane wave. ; ; Ellipticity:The ratio (minor axis)/(major axis) of the ; ellipse transcribed by the field variations of the ; components transverse to the Z direction. The sign ; indicates the direction of rotation of the field vector in ; the plane. Negative signs refer to left-handed ; rotation about the Z direction. In the field ; aligned coordinate system these signs refer to ; plasma waves of left and right handed ; polarisation. ; ; Helicity:Similar to Ellipticity except defined in terms of the ; direction of minimum variance instead of Z. Stricltly the Helicity ; is defined in terms of the wavenormal direction or k. ; However since from single point observations the ; sense of k cannot be determined, helicity here is ; simply the ratio of the minor to major axis transverse to the ; minimum variance direction without sign. ; ; ;RESTRICTIONS:-If one component is an order of magnitude or more greater than ; the other two then the polarisation results saturate and erroneously ; indicate high degrees of polarisation at all times and ; frequencies. Time series should be eyeballed before running the program. ; For time series containing very rapid changes or spikes ; the usual problems with Fourier analysis arise. ; Care should be taken in evaluating degree of polarisation results. ; For meaningful results there should be significant wave power at the ; frequency where the polarisation approaches ; 100%. Remembercomparing two straight lines yields 100% polarisation. ; ; $LastChangedBy: jwl $ ; $LastChangedDate: 2012-09-17 10:32:23 -0700 (Mon, 17 Sep 2012) $ ; $LastChangedRevision: 10921 $ ; $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/ssl_general/tags/tdas_8_00/science/wavpol/wavpol.pro $ ;- pro wavpol,ct,Bx,By,Bz,timeline,freqline,powspec,degpol,waveangle,elliptict,helict,$ nopfft=nopfft_input,steplength = steplength_input nopoints=n_elements(Bx) If size(nopfft_input, /type) ne 0 then nopfft = nopfft_input else nopfft = 256 If size(steplength_input, /type) ne 0 then steplength = steplength_input else steplength = nopfft/2 ; ;steplength =128 ;overlap between successive FFT intervals ;nopfft=256 ;no. of points in FFT nosteps=(nopoints-nopfft)/steplength ;total number of FFTs leveltplot=0.000001 ;power rejection level 0 to 1 lam=dblarr(2) nosmbins=7 ;No. of bins in frequency domain ;!p.charsize=2 ;to include in smoothing (must be odd) aa=[0.024,0.093,0.232,0.301,0.232,0.093,0.024] ;smoothing profile based on Hanning ; ;ARRAY DEFINITIONS fields=make_array(3, nopoints,/double) power=Make_Array(nosteps,nopfft/2) specx=Make_Array(nosteps,nopfft,/dcomplex) specy=Make_Array(nosteps,nopfft,/dcomplex) specz=Make_Array(nosteps,nopfft,/dcomplex) wnx=DBLARR(nosteps,nopfft/2) wny=DBLARR(nosteps,nopfft/2) wnz=DBLARR(nosteps,nopfft/2) vecspec=Make_Array(nosteps,nopfft/2,3,/dcomplex) matspec=Make_Array(nosteps,nopfft/2,3,3,/dcomplex) ematspec=Make_Array(nosteps,nopfft/2,3,3,/dcomplex) matsqrd=Make_Array(nosteps,nopfft/2,3,3,/dcomplex) matsqrd1=Make_Array(nosteps,nopfft/2,3,3,/dcomplex) trmatspec=Make_Array(nosteps,nopfft/2,/double) powspec=Make_Array(nosteps,nopfft/2,/double) trmatsqrd=Make_Array(nosteps,nopfft/2,/double) degpol=Make_Array(nosteps,nopfft/2,/double) alpha=Make_Array(nosteps,nopfft/2,/double) alphasin2=Make_Array(nosteps,nopfft/2,/double) alphacos2=Make_Array(nosteps,nopfft/2,/double) alphasin3=Make_Array(nosteps,nopfft/2,/double) alphacos3=Make_Array(nosteps,nopfft/2,/double) alphax=Make_Array(nosteps,nopfft/2,/double) alphasin2x=Make_Array(nosteps,nopfft/2,/double) alphacos2x=Make_Array(nosteps,nopfft/2,/double) alphasin3x=Make_Array(nosteps,nopfft/2,/double) alphacos3x=Make_Array(nosteps,nopfft/2,/double) alphay=Make_Array(nosteps,nopfft/2,/double) alphasin2y=Make_Array(nosteps,nopfft/2,/double) alphacos2y=Make_Array(nosteps,nopfft/2,/double) alphasin3y=Make_Array(nosteps,nopfft/2,/double) alphacos3y=Make_Array(nosteps,nopfft/2,/double) alphaz=Make_Array(nosteps,nopfft/2,/double) alphasin2z=Make_Array(nosteps,nopfft/2,/double) alphacos2z=Make_Array(nosteps,nopfft/2,/double) alphasin3z=Make_Array(nosteps,nopfft/2,/double) alphacos3z=Make_Array(nosteps,nopfft/2,/double) gamma=Make_Array(nosteps,nopfft/2,/double) gammarot=Make_Array(nosteps,nopfft/2,/double) upper=Make_Array(nosteps,nopfft/2,/double) lower=Make_Array(nosteps,nopfft/2,/double) lambdau=Make_Array(nosteps,nopfft/2,3,3,/dcomplex) lambdaurot=Make_Array(nosteps,nopfft/2,2,/dcomplex) thetarot=Make_Array(nopfft/2,/double) thetax=DBLARR(nosteps,nopfft) thetay=DBLARR(nosteps,nopfft) thetaz=DBLARR(nosteps,nopfft) aaa2=DBLARR(nosteps,nopfft/2) helicity=Make_Array(nosteps,nopfft/2,3) ellip=Make_Array(nosteps,nopfft/2,3) waveangle=Make_Array(nosteps,nopfft/2) halfspecx=Make_Array(nosteps,nopfft/2,/dcomplex) halfspecy=Make_Array(nosteps,nopfft/2,/dcomplex) halfspecz=Make_Array(nosteps,nopfft/2,/dcomplex) ; ; DEFINE ARRAYS ; xs=Bx & ys=By & zs=Bz sampfreq=1/(ct(1)-ct(0)) endsampfreq=1/(ct(nopoints-1)-ct(nopoints-2)) if sampfreq NE endsampfreq then dprint, 'Warning: file sampling ' + $ 'frequency changes',sampfreq,'Hz to',endsampfreq,'Hz' else dprint, 'ac ' + $ 'file sampling frequency',sampfreq,'Hz' ; print,' ' dprint, 'Total number of steps',nosteps print,' ' counter_start = 0 for j=0,(nosteps-1) do begin if 10*double(j)/(nosteps-1) gt (counter_start+1) then begin dprint, strtrim(100*double(j)/(nosteps-1),2) + ' % Complete ' dprint, ' Processing step no. :'+ strtrim(j+1,2) counter_start++ endif ; ;FFT CALCULATION smooth=0.08+0.46*(1-cos(2*!DPI*findgen(nopfft)/nopfft)) tempx=smooth*xs(0:nopfft-1) tempy=smooth*ys(0:nopfft-1) tempz=smooth*zs(0:nopfft-1) specx(j,*)=(fft(tempx,/double));+Complex(0,j*steplength*3.1415/32)) specy(j,*)=(fft(tempy,/double));+Complex(0,j*steplength*3.1415/32)) specz(j,*)=(fft(tempz,/double));+Complex(0,j*steplength*3.1415/32)) halfspecx(j,*)=specx(j,0:(nopfft/2-1)) halfspecy(j,*)=specy(j,0:(nopfft/2-1)) halfspecz(j,*)=specz(j,0:(nopfft/2-1)) xs=shift(xs,-steplength) ys=shift(ys,-steplength) zs=shift(zs,-steplength) ;CALCULATION OF THE SPECTRAL MATRIX matspec(j,*,0,0)=halfspecx(j,*)*conj(halfspecx(j,*)) matspec(j,*,1,0)=halfspecx(j,*)*conj(halfspecy(j,*)) matspec(j,*,2,0)=halfspecx(j,*)*conj(halfspecz(j,*)) matspec(j,*,0,1)=halfspecy(j,*)*conj(halfspecx(j,*)) matspec(j,*,1,1)=halfspecy(j,*)*conj(halfspecy(j,*)) matspec(j,*,2,1)=halfspecy(j,*)*conj(halfspecz(j,*)) matspec(j,*,0,2)=halfspecz(j,*)*conj(halfspecx(j,*)) matspec(j,*,1,2)=halfspecz(j,*)*conj(halfspecy(j,*)) matspec(j,*,2,2)=halfspecz(j,*)*conj(halfspecz(j,*)) ;CALCULATION OF SMOOTHED SPECTRAL MATRIX for k=(nosmbins-1)/2, (nopfft/2-1)-(nosmbins-1)/2 do begin ematspec(j,k,0,0)=TOTAL(aa(0:(nosmbins-1))*matspec(j,(k-(nosmbins-1)/2):(k+(nosmbins-1)/2),0,0)) ematspec(j,k,1,0)=TOTAL(aa(0:(nosmbins-1))*matspec(j,(k-(nosmbins-1)/2):(k+(nosmbins-1)/2),1,0)) ematspec(j,k,2,0)=TOTAL(aa(0:(nosmbins-1))*matspec(j,(k-(nosmbins-1)/2):(k+(nosmbins-1)/2),2,0)) ematspec(j,k,0,1)=TOTAL(aa(0:(nosmbins-1))*matspec(j,(k-(nosmbins-1)/2):(k+(nosmbins-1)/2),0,1)) ematspec(j,k,1,1)=TOTAL(aa(0:(nosmbins-1))*matspec(j,(k-(nosmbins-1)/2):(k+(nosmbins-1)/2),1,1)) ematspec(j,k,2,1)=TOTAL(aa(0:(nosmbins-1))*matspec(j,(k-(nosmbins-1)/2):(k+(nosmbins-1)/2),2,1)) ematspec(j,k,0,2)=TOTAL(aa(0:(nosmbins-1))*matspec(j,(k-(nosmbins-1)/2):(k+(nosmbins-1)/2),0,2)) ematspec(j,k,1,2)=TOTAL(aa(0:(nosmbins-1))*matspec(j,(k-(nosmbins-1)/2):(k+(nosmbins-1)/2),1,2)) ematspec(j,k,2,2)=TOTAL(aa(0:(nosmbins-1))*matspec(j,(k-(nosmbins-1)/2):(k+(nosmbins-1)/2),2,2)) endfor ;CALCULATION OF THE MINIMUM VARIANCE DIRECTION AND WAVENORMAL ANGLE aaa2(j,*)=SQRT(IMAGINARY(ematspec(j,*,0,1))^2+IMAGINARY(ematspec(j,*,0,2))^2+IMAGINARY(ematspec(j,*,1,2))^2) wnx(j,*)=ABS(IMAGINARY(ematspec(j,*,1,2))/aaa2(j,*)) wny(j,*)=-ABS(IMAGINARY(ematspec(j,*,0,2))/aaa2(j,*)) wnz(j,*)=IMAGINARY(ematspec(j,*,0,1))/aaa2(j,*) waveangle(j,*)=ATAN(Sqrt(wnx(j,*)^2+wny(j,*)^2),abs(wnz(j,*))) ;CALCULATION OF THE DEGREE OF POLARISATION ;calc of square of smoothed spec matrix matsqrd(j,*,0,0)=ematspec(j,*,0,0)*ematspec(j,*,0,0)+ematspec(j,*,0,1)*ematspec(j,*,1,0)+ematspec(j,*,0,2)*ematspec(j,*,2,0) matsqrd(j,*,0,1)=ematspec(j,*,0,0)*ematspec(j,*,0,1)+ematspec(j,*,0,1)*ematspec(j,*,1,1)+ematspec(j,*,0,2)*ematspec(j,*,2,1) matsqrd(j,*,0,2)=ematspec(j,*,0,0)*ematspec(j,*,0,2)+ematspec(j,*,0,1)*ematspec(j,*,1,2)+ematspec(j,*,0,2)*ematspec(j,*,2,2) matsqrd(j,*,1,0)=ematspec(j,*,1,0)*ematspec(j,*,0,0)+ematspec(j,*,1,1)*ematspec(j,*,1,0)+ematspec(j,*,1,2)*ematspec(j,*,2,0) matsqrd(j,*,1,1)=ematspec(j,*,1,0)*ematspec(j,*,0,1)+ematspec(j,*,1,1)*ematspec(j,*,1,1)+ematspec(j,*,1,2)*ematspec(j,*,2,1) matsqrd(j,*,1,2)=ematspec(j,*,1,0)*ematspec(j,*,0,2)+ematspec(j,*,1,1)*ematspec(j,*,1,2)+ematspec(j,*,1,2)*ematspec(j,*,2,2) matsqrd(j,*,2,0)=ematspec(j,*,2,0)*ematspec(j,*,0,0)+ematspec(j,*,2,1)*ematspec(j,*,1,0)+ematspec(j,*,2,2)*ematspec(j,*,2,0) matsqrd(j,*,2,1)=ematspec(j,*,2,0)*ematspec(j,*,0,1)+ematspec(j,*,2,1)*ematspec(j,*,1,1)+ematspec(j,*,2,2)*ematspec(j,*,2,1) matsqrd(j,*,2,2)=ematspec(j,*,2,0)*ematspec(j,*,0,2)+ematspec(j,*,2,1)*ematspec(j,*,1,2)+ematspec(j,*,2,2)*ematspec(j,*,2,2) Trmatsqrd(j,*)=matsqrd(j,*,0,0)+matsqrd(j,*,1,1)+matsqrd(j,*,2,2) Trmatspec(j,*)=ematspec(j,*,0,0)+ematspec(j,*,1,1)+ematspec(j,*,2,2) degpol(j,(nosmbins-1)/2:(nopfft/2-1)-(nosmbins-1)/2)=(3*Trmatsqrd(j,(nosmbins-1)/2:(nopfft/2-1)-(nosmbins-1)/2)-Trmatspec(j,(nosmbins-1)/2: (nopfft/2-1)-(nosmbins-1)/2)^2)/(2*Trmatspec(j,(nosmbins-1)/2: (nopfft/2-1)-(nosmbins-1)/2)^2) ;CALCULATION OF HELICITY, ELLIPTICITY AND THE WAVE STATE VECTOR alphax(j,*)=Sqrt(ematspec(j,*,0,0)) alphacos2x(j,*)=Double(ematspec(j,*,0,1))/Sqrt(ematspec(j,*,0,0)) alphasin2x(j,*)=-Imaginary(ematspec(j,*,0,1))/Sqrt(ematspec(j,*,0,0)) alphacos3x(j,*)=Double(ematspec(j,*,0,2))/Sqrt(ematspec(j,*,0,0)) alphasin3x(j,*)=-Imaginary(ematspec(j,*,0,2))/Sqrt(ematspec(j,*,0,0)) lambdau(j,*,0,0)=alphax(j,*) lambdau(j,*,0,1)=Complex(alphacos2x(j,*),alphasin2x(j,*)) lambdau(j,*,0,2)=Complex(alphacos3x(j,*),alphasin3x(j,*)) alphay(j,*)=Sqrt(ematspec(j,*,1,1)) alphacos2y(j,*)=Double(ematspec(j,*,1,0))/Sqrt(ematspec(j,*,1,1)) alphasin2y(j,*)=-Imaginary(ematspec(j,*,1,0))/Sqrt(ematspec(j,*,1,1)) alphacos3y(j,*)=Double(ematspec(j,*,1,2))/Sqrt(ematspec(j,*,1,1)) alphasin3y(j,*)=-Imaginary(ematspec(j,*,1,2))/Sqrt(ematspec(j,*,1,1)) lambdau(j,*,1,0)=alphay(j,*) lambdau(j,*,1,1)=Complex(alphacos2y(j,*),alphasin2y(j,*)) lambdau(j,*,1,2)=Complex(alphacos3y(j,*),alphasin3y(j,*)) alphaz(j,*)=Sqrt(ematspec(j,*,2,2)) alphacos2z(j,*)=Double(ematspec(j,*,2,0))/Sqrt(ematspec(j,*,2,2)) alphasin2z(j,*)=-Imaginary(ematspec(j,*,2,0))/Sqrt(ematspec(j,*,2,2)) alphacos3z(j,*)=Double(ematspec(j,*,2,1))/Sqrt(ematspec(j,*,2,2)) alphasin3z(j,*)=-Imaginary(ematspec(j,*,2,1))/Sqrt(ematspec(j,*,2,2)) lambdau(j,*,2,0)=alphaz(j,*) lambdau(j,*,2,1)=Complex(alphacos2z(j,*),alphasin2z(j,*)) lambdau(j,*,2,2)=Complex(alphacos3z(j,*),alphasin3z(j,*)) ;HELICITY CALCULATION for k=0, nopfft/2-1 do begin for xyz=0,2 do begin upper(j,k)=Total(2*double(lambdau(j,k,xyz,0:2))*(Imaginary(lambdau(j,k,xyz,0:2)))) lower(j,k)=Total((Double(lambdau(j,k,xyz,0:2)))^2-(Imaginary(lambdau(j,k,xyz,0:2)))^2) if (upper(j,k) GT 0.00) then gamma(j,k)=ATAN(upper(j,k),lower(j,k)) else gamma(j,k)=!DPI+(!DPI+ATAN(upper(j,k),lower(j,k))) lambdau(j,k,xyz,*)=exp(Complex(0,-0.5*gamma(j,k)))*lambdau(j,k,xyz,*) helicity(j,k,xyz)=1/(SQRT(Double(lambdau(j,k,xyz,0))^2+Double(lambdau(j,k,xyz,1))^2+Double(lambdau(j,k,xyz,2))^2)/SQRT(Imaginary(lambdau(j,k,xyz,0))^2+Imaginary(lambdau(j,k,xyz,1))^2+Imaginary(lambdau(j,k,xyz,2))^2)) ;ELLIPTICITY CALCULATION uppere=Imaginary(lambdau(j,k,xyz,0))*Double(lambdau(j,k,xyz,0))+Imaginary(lambdau(j,k,xyz,1))*Double(lambdau(j,k,xyz,1)) lowere=-Imaginary(lambdau(j,k,xyz,0))^2+Double(lambdau(j,k,xyz,0))^2-Imaginary(lambdau(j,k,xyz,1))^2+Double(lambdau(j,k,xyz,1))^2 if uppere GT 0 then gammarot(j,k)=ATAN(uppere,lowere) else gammarot(j,k)=!DPI+!DPI+ATAN(uppere,lowere) lam=lambdau(j,k,xyz,0:1) lambdaurot(j,k,*)=exp(complex(0,-0.5*gammarot(j,k)))*lam(*) ellip(j,k,xyz)=Sqrt(Imaginary(lambdaurot(j,k,0))^2+Imaginary(lambdaurot(j,k,1))^2)/Sqrt(Double(lambdaurot(j,k,0))^2+Double(lambdaurot(j,k,1))^2) ellip(j,k,xyz)=-ellip(j,k,xyz)*(Imaginary(ematspec(j,k,0,1))*sin(waveangle(j,k)))/abs(Imaginary(ematspec(j,k,0,1))*sin(waveangle(j,k))) endfor endfor endfor ; end of main body ;AVERAGING HELICITY AND ELLIPTICITY RESULTS elliptict=(ellip(*,*,0)+ellip(*,*,1)+ellip(*,*,2))/3 helict=(helicity(*,*,0)+helicity(*,*,1)+helicity(*,*,2))/3 ; CREATING OUTPUT STRUCTURES ; timeline=ct(0)+ABS(nopfft/2)/sampfreq+findgen(nosteps)*steplength/sampfreq binwidth=sampfreq/nopfft freqline=binwidth*findgen(nopfft/2) ;scaling power results to units with meaning ;W=nopfft*Total(smooth^2); original Chaston ; redefining W ; 8Sep2012 jhl W=Total(smooth^2) / double(nopfft) ; switch to divide by nopfft powspec(*,1:nopfft/2-2)=1/W*2*trmatspec(*,1:nopfft/2-2)/binwidth powspec(*,0)=1/W*trmatspec(*,0)/binwidth powspec(*,nopfft/2-1)=1/W*trmatspec(*,nopfft/2-1)/binwidth return end