;+
;
; 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: pcruce $
; $LastChangedDate: 2009-06-17 15:35:34 -0700 (Wed, 17 Jun 2009) $
; $LastChangedRevision: 6244 $
; $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/ssl_general/trunk/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 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,' '
counter_start = 0
for j=0,(nosteps-1) do begin

if 10*double(j)/(nosteps-1) gt (counter_start+1) then begin
  print,strtrim(100*double(j)/(nosteps-1),2) + ' % Complete '
  print,' 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)
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
