;+
;
; NAME:wavpol
;
; MODIFICATION HISTORY:Written By Chris Chaston, 30-10-96
; :Modified by Vassilis, 2001-07-11
;
;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
;
;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: pcruce $
; $LastChangedDate: 2007-11-11 16:40:33 -0800 (Sun, 11 Nov 2007) $
; $LastChangedRevision: 2025 $
; $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/ssl_general/trunk/cotrans/special/fac/fac_matrix_make.pro $
;-
pro wavpol,ct,Bx,By,Bz,timeline,freqline,powspec,degpol,waveangle,elliptict,helict
nopoints=n_elements(Bx)
;
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 print,'Warning: file sampling ' + $
'frequency changes',sampfreq,'Hz to',endsampfreq,'Hz' else print,'ac ' + $
'file sampling frequency',sampfreq,'Hz'
;
print,' '
print, 'Total number of steps',nosteps
print,' '
for j=0,(nosteps-1) do begin
print,'processing step no.',j+1
;
;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)
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