-
Notifications
You must be signed in to change notification settings - Fork 0
/
wavpol.pro
480 lines (420 loc) · 24.7 KB
/
wavpol.pro
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
;+
;
; 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
; :Modified by O. Le Contel, 2016-03
; to be able to change frequency averaging parameter by adding bin_freq keyword
; :Modified by O. Le Contel, LPP, 2016-07, in order to manage data gaps in the waveform
; using test written by K. Bromund (in thm_cal_scm)
; :Modified by egrimes, merging OLE's changes with SPEDAS's wavpol:
; gamma -> gammay (avoids confusion with IDL's gamma function)
; redefined W (fixed bug reported by Justin Lee) -> W=Total(smooth^2) / double(nopfft)
; added pspec3 input, pspec[x,y,z], returns pspec3 (changes from Justin Lee) - added to original 10/10/2013
; converted () to [], fixed tabbing
; updated documentation with changes from Justin Lee - added to original 9/23/2014
; :Modified by egrimes, now checking for 0s in the output time series, setting those
; data values to NaNs
;
;
;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
;
; bin_freq (optional): No. of bins in frequency domain
;
;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 (units of nT^2/Hz if input Bx, By, Bz are in nT)
;
; 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. See:
; Means, J. D. (1972), Use of the three-dimensional covariance
; matrix in analyzing the polarization properties of plane waves,
; J. Geophys. Res., 77(28), 5551-5559,
; doi:10.1029/JA077i028p05551.
;
; Ellipticity:
; The ratio (minor axis)/(major axis) of the ellipse transcribed
; by the field variations of the components transverse to the
; Z direction (Samson and Olson, 1980). The sign indicates
; the direction of rotation of the field vector in the plane (cf.
; Means, (1972)).
; 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 polarization.
;
; 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: egrimes $
; $LastChangedDate: 2022-08-18 14:05:08 -0700 (Thu, 18 Aug 2022) $
; $LastChangedRevision: 31024 $
; $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/spdsoft/trunk/general/science/wavpol/wavpol.pro $
;-
pro wavpol,ct,Bx,By,Bz,timeline,freqline,powspec,degpol,waveangle,elliptict,helict,pspec3,$
samp_per=samp_per_input, nopfft=nopfft_input,steplength = steplength_input, bin_freq = bin_freq_input, err_flag = err_flag
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
If size(bin_freq_input, /type) ne 0 then bin_freq = bin_freq_input else bin_freq = 3
If size(err_flag, /type) ne 0 then err_flag = err_flag else err_flag = 0
nopoints=n_elements(Bx)
iano = intarr(nopoints)
dt = ct[1:*]-ct
beginsampfreq=1d/(ct[1]-ct[0])
endsampfreq=1d/(ct[nopoints-1]-ct[nopoints-2])
if beginsampfreq NE endsampfreq then dprint,'Warning: file sampling ' + $
'frequency changes',beginsampfreq,'Hz to',endsampfreq,'Hz' else dprint,'ac ' + $
'file sampling frequency',beginsampfreq,'Hz'
samp_freq = beginsampfreq
samp_per = 1d/samp_freq
; ========= Time anomaly detection from test written by K. Bromund (in thm_cal_scm)
; ========= time reversal detection
reverse = where(dt lt 0, n_reverse)
if n_reverse gt 0 then iano[reverse] = 16 ; time reverse
; ========= discontinuity in time detection
accuracy = 0.01 ; the accuracy of the sampling frequency should be about 1 percent.
discont_trigger = accuracy*samp_per
;discontinuity = where(abs(dt-1.0/samp_freq) gt 2.0d-5, n_discont)
discontinuity = where(abs(dt-1d/samp_freq) gt discont_trigger, n_discont)
if n_discont gt 0 then iano[discontinuity] = 17 ; discontinuity in time
; ========= end of file
iano[nopoints-1] = 22
;; discontinuity in sample rate (Minor bug: some changes in sample rate
;; can be mislabled as discontinuities in time)
;iano[ind_r[n_elements(ind_r)-1]] = 22
errs = where(iano ge 15, n_batches)
dprint,''
dprint, 'Number of continuous batches: ', n_batches
if n_batches gt 8d4 then begin
dprint,''
dprint, 'Large number of batches. Returning to avoid memory runaway.'
err_flag = 1;
return
endif
nbp_fft_batches = lonarr(n_batches)*0.
; Total numbers of FFT calculations including 1 leap frog for each batch
ind_batch0 = 0
nosteps = 0
for batch = 0,n_batches-1 do begin
nosteps = nosteps+floor((errs[batch] - ind_batch0)/steplength,/L64)
ind_batch0 = errs[batch]
endfor
nosteps= nosteps+n_batches
;nosteps=(nopoints-nopfft)/steplength ;total number of FFTs
leveltplot=0.000001 ;power rejection level 0 to 1
lam=dblarr(2)
nosmbins=bin_freq ;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
timeline = dblarr(nosteps)*0.
;
;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=Make_Array(nosteps,nopfft/2);DBLARR(nosteps,nopfft/2)
wny=Make_Array(nosteps,nopfft/2);DBLARR(nosteps,nopfft/2)
wnz=Make_Array(nosteps,nopfft/2);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)
xrmatspec=Make_Array(nosteps,nopfft/2,/double) ; added 10Sep2013 jhl
yrmatspec=Make_Array(nosteps,nopfft/2,/double) ; added 10Sep2013 jhl
zrmatspec=Make_Array(nosteps,nopfft/2,/double) ; added 10Sep2013 jhl
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)
gammay=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)
pspecx=Make_Array(nosteps,nopfft/2,/double)
pspecy=Make_Array(nosteps,nopfft/2,/double)
pspecz=Make_Array(nosteps,nopfft/2,/double)
pspec3=Make_Array(nosteps,nopfft/2,3,/double)
;
; DEFINE ARRAYS
;
;xs=Bx & ys=By & zs=Bz
;
dprint,' '
dprint, 'Total number of steps',nosteps
dprint,' '
ind0 = 0L
KK = 0L
for batch = 0l, n_batches-1 do begin
ind1 = errs[batch]
nbp_batch = ind1-ind0+1L
ind1_ref = ind1
KK_batch_start = KK
xs = Bx[ind0:ind1]
ys = By[ind0:ind1]
zs = Bz[ind0:ind1]
good_data = where(finite(xs), ngood,/L64)
if(ngood Gt nopfft) then begin
;=== Number of fft spectra
nbp_fft_batches[batch] = floor(ngood/steplength,/L64)
dprint,'Total number of possible FFT in the batch n°', batch,' is:',nbp_fft_batches[batch]
ind0_fft=0L
for j=0L,nbp_fft_batches[batch]-1L do begin
ind1_fft = nopfft*(j+1L)-1L
ind1_ref_fft = ind1_fft
;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[KK,*]=(fft(tempx,/double));+Complex(0,j*steplength*3.1415/32))
specy[KK,*]=(fft(tempy,/double));+Complex(0,j*steplength*3.1415/32))
specz[KK,*]=(fft(tempz,/double));+Complex(0,j*steplength*3.1415/32))
halfspecx[KK,*]=specx[KK,0:(nopfft/2-1)]
halfspecy[KK,*]=specy[KK,0:(nopfft/2-1)]
halfspecz[KK,*]=specz[KK,0:(nopfft/2-1)]
xs=shift(xs,-steplength)
ys=shift(ys,-steplength)
zs=shift(zs,-steplength)
;CALCULATION OF THE SPECTRAL MATRIX
matspec[KK,*,0,0]=halfspecx[KK,*]*conj(halfspecx[KK,*])
matspec[KK,*,1,0]=halfspecx[KK,*]*conj(halfspecy[KK,*])
matspec[KK,*,2,0]=halfspecx[KK,*]*conj(halfspecz[KK,*])
matspec[KK,*,0,1]=halfspecy[KK,*]*conj(halfspecx[KK,*])
matspec[KK,*,1,1]=halfspecy[KK,*]*conj(halfspecy[KK,*])
matspec[KK,*,2,1]=halfspecy[KK,*]*conj(halfspecz[KK,*])
matspec[KK,*,0,2]=halfspecz[KK,*]*conj(halfspecx[KK,*])
matspec[KK,*,1,2]=halfspecz[KK,*]*conj(halfspecy[KK,*])
matspec[KK,*,2,2]=halfspecz[KK,*]*conj(halfspecz[KK,*])
;CALCULATION OF SMOOTHED SPECTRAL MATRIX
for k=(nosmbins-1)/2, (nopfft/2-1)-(nosmbins-1)/2 do begin
ematspec[KK,k,0,0]=TOTAL(aa[0:(nosmbins-1)]*matspec[KK,(k-(nosmbins-1)/2):(k+(nosmbins-1)/2),0,0])
ematspec[KK,k,1,0]=TOTAL(aa[0:(nosmbins-1)]*matspec[KK,(k-(nosmbins-1)/2):(k+(nosmbins-1)/2),1,0])
ematspec[KK,k,2,0]=TOTAL(aa[0:(nosmbins-1)]*matspec[KK,(k-(nosmbins-1)/2):(k+(nosmbins-1)/2),2,0])
ematspec[KK,k,0,1]=TOTAL(aa[0:(nosmbins-1)]*matspec[KK,(k-(nosmbins-1)/2):(k+(nosmbins-1)/2),0,1])
ematspec[KK,k,1,1]=TOTAL(aa[0:(nosmbins-1)]*matspec[KK,(k-(nosmbins-1)/2):(k+(nosmbins-1)/2),1,1])
ematspec[KK,k,2,1]=TOTAL(aa[0:(nosmbins-1)]*matspec[KK,(k-(nosmbins-1)/2):(k+(nosmbins-1)/2),2,1])
ematspec[KK,k,0,2]=TOTAL(aa[0:(nosmbins-1)]*matspec[KK,(k-(nosmbins-1)/2):(k+(nosmbins-1)/2),0,2])
ematspec[KK,k,1,2]=TOTAL(aa[0:(nosmbins-1)]*matspec[KK,(k-(nosmbins-1)/2):(k+(nosmbins-1)/2),1,2])
ematspec[KK,k,2,2]=TOTAL(aa[0:(nosmbins-1)]*matspec[KK,(k-(nosmbins-1)/2):(k+(nosmbins-1)/2),2,2])
endfor
;CALCULATION OF THE MINIMUM VARIANCE DIRECTION AND WAVENORMAL ANGLE
aaa2[KK,*]=SQRT(IMAGINARY(ematspec[KK,*,0,1])^2+IMAGINARY(ematspec[KK,*,0,2])^2+IMAGINARY(ematspec[KK,*,1,2])^2)
wnx[KK,*]=ABS(IMAGINARY(ematspec[KK,*,1,2])/aaa2[KK,*])
wny[KK,*]=-ABS(IMAGINARY(ematspec[KK,*,0,2])/aaa2[KK,*])
wnz[KK,*]=IMAGINARY(ematspec[KK,*,0,1])/aaa2[KK,*]
waveangle[KK,*]=ATAN(Sqrt(wnx[KK,*]^2+wny[KK,*]^2),abs(wnz[KK,*]))
;CALCULATION OF THE DEGREE OF POLARISATION
;calc of square of smoothed spec matrix
matsqrd[KK,*,0,0]=ematspec[KK,*,0,0]*ematspec[KK,*,0,0]+ematspec[KK,*,0,1]*ematspec[KK,*,1,0]+ematspec[KK,*,0,2]*ematspec[KK,*,2,0]
matsqrd[KK,*,0,1]=ematspec[KK,*,0,0]*ematspec[KK,*,0,1]+ematspec[KK,*,0,1]*ematspec[KK,*,1,1]+ematspec[KK,*,0,2]*ematspec[KK,*,2,1]
matsqrd[KK,*,0,2]=ematspec[KK,*,0,0]*ematspec[KK,*,0,2]+ematspec[KK,*,0,1]*ematspec[KK,*,1,2]+ematspec[KK,*,0,2]*ematspec[KK,*,2,2]
matsqrd[KK,*,1,0]=ematspec[KK,*,1,0]*ematspec[KK,*,0,0]+ematspec[KK,*,1,1]*ematspec[KK,*,1,0]+ematspec[KK,*,1,2]*ematspec[KK,*,2,0]
matsqrd[KK,*,1,1]=ematspec[KK,*,1,0]*ematspec[KK,*,0,1]+ematspec[KK,*,1,1]*ematspec[KK,*,1,1]+ematspec[KK,*,1,2]*ematspec[KK,*,2,1]
matsqrd[KK,*,1,2]=ematspec[KK,*,1,0]*ematspec[KK,*,0,2]+ematspec[KK,*,1,1]*ematspec[KK,*,1,2]+ematspec[KK,*,1,2]*ematspec[KK,*,2,2]
matsqrd[KK,*,2,0]=ematspec[KK,*,2,0]*ematspec[KK,*,0,0]+ematspec[KK,*,2,1]*ematspec[KK,*,1,0]+ematspec[KK,*,2,2]*ematspec[KK,*,2,0]
matsqrd[KK,*,2,1]=ematspec[KK,*,2,0]*ematspec[KK,*,0,1]+ematspec[KK,*,2,1]*ematspec[KK,*,1,1]+ematspec[KK,*,2,2]*ematspec[KK,*,2,1]
matsqrd[KK,*,2,2]=ematspec[KK,*,2,0]*ematspec[KK,*,0,2]+ematspec[KK,*,2,1]*ematspec[KK,*,1,2]+ematspec[KK,*,2,2]*ematspec[KK,*,2,2]
Trmatsqrd[KK,*]=matsqrd[KK,*,0,0]+matsqrd[KK,*,1,1]+matsqrd[KK,*,2,2]
Trmatspec[KK,*]=ematspec[KK,*,0,0]+ematspec[KK,*,1,1]+ematspec[KK,*,2,2]
degpol[KK,(nosmbins-1)/2:(nopfft/2-1)-(nosmbins-1)/2]=(3*Trmatsqrd[KK,(nosmbins-1)/2:(nopfft/2-1)-(nosmbins-1)/2]-Trmatspec[KK,(nosmbins-1)/2: (nopfft/2-1)-(nosmbins-1)/2]^2)/(2*Trmatspec[KK,(nosmbins-1)/2: (nopfft/2-1)-(nosmbins-1)/2]^2)
Xrmatspec[KK,*]=ematspec[KK,*,0,0] ; added 10Sep2013 jhl
Yrmatspec[KK,*]=ematspec[KK,*,1,1] ; added 10Sep2013 jhl
Zrmatspec[KK,*]=ematspec[KK,*,2,2] ; added 10Sep2013 jhl
;CALCULATION OF HELICITY, ELLIPTICITY AND THE WAVE STATE VECTOR
alphax[KK,*]=Sqrt(ematspec[KK,*,0,0])
alphacos2x[KK,*]=Double(ematspec[KK,*,0,1])/Sqrt(ematspec[KK,*,0,0])
alphasin2x[KK,*]=-Imaginary(ematspec[KK,*,0,1])/Sqrt(ematspec[KK,*,0,0])
alphacos3x[KK,*]=Double(ematspec[KK,*,0,2])/Sqrt(ematspec[KK,*,0,0])
alphasin3x[KK,*]=-Imaginary(ematspec[KK,*,0,2])/Sqrt(ematspec[KK,*,0,0])
lambdau[KK,*,0,0]=alphax[KK,*]
lambdau[KK,*,0,1]=Complex(alphacos2x[KK,*],alphasin2x[KK,*])
lambdau[KK,*,0,2]=Complex(alphacos3x[KK,*],alphasin3x[KK,*])
alphay[KK,*]=Sqrt(ematspec[KK,*,1,1])
alphacos2y[KK,*]=Double(ematspec[KK,*,1,0])/Sqrt(ematspec[KK,*,1,1])
alphasin2y[KK,*]=-Imaginary(ematspec[KK,*,1,0])/Sqrt(ematspec[KK,*,1,1])
alphacos3y[KK,*]=Double(ematspec[KK,*,1,2])/Sqrt(ematspec[KK,*,1,1])
alphasin3y[KK,*]=-Imaginary(ematspec[KK,*,1,2])/Sqrt(ematspec[KK,*,1,1])
lambdau[KK,*,1,0]=alphay[KK,*]
lambdau[KK,*,1,1]=Complex(alphacos2y[KK,*],alphasin2y[KK,*])
lambdau[KK,*,1,2]=Complex(alphacos3y[KK,*],alphasin3y[KK,*])
alphaz[KK,*]=Sqrt(ematspec[KK,*,2,2])
alphacos2z[KK,*]=Double(ematspec[KK,*,2,0])/Sqrt(ematspec[KK,*,2,2])
alphasin2z[KK,*]=-Imaginary(ematspec[KK,*,2,0])/Sqrt(ematspec[KK,*,2,2])
alphacos3z[KK,*]=Double(ematspec[KK,*,2,1])/Sqrt(ematspec[KK,*,2,2])
alphasin3z[KK,*]=-Imaginary(ematspec[KK,*,2,1])/Sqrt(ematspec[KK,*,2,2])
lambdau[KK,*,2,0]=alphaz[KK,*]
lambdau[KK,*,2,1]=Complex(alphacos2z[KK,*],alphasin2z[KK,*])
lambdau[KK,*,2,2]=Complex(alphacos3z[KK,*],alphasin3z[KK,*])
;HELICITY CALCULATION
for k=0, nopfft/2-1 do begin
for xyz=0,2 do begin
upper[KK,k]=Total(2*double(lambdau[KK,k,xyz,0:2])*(Imaginary(lambdau[KK,k,xyz,0:2])),/NAN) ; Add /NAN OLe 2016
lower[KK,k]=Total((Double(lambdau[KK,k,xyz,0:2]))^2-(Imaginary(lambdau[KK,k,xyz,0:2]))^2,/NAN) ; Add /NAN OLE 2016
if (upper[KK,k] GT 0.00) then gammay[KK,k]=ATAN(upper[KK,k],lower[KK,k]) else gammay[KK,k]=!DPI+(!DPI+ATAN(upper[KK,k],lower[KK,k]))
lambdau[KK,k,xyz,*]=exp(Complex(0,-0.5*gammay[KK,k]))*lambdau[KK,k,xyz,*]
helicity[KK,k,xyz]=1/(SQRT(Double(lambdau[KK,k,xyz,0])^2+Double(lambdau[KK,k,xyz,1])^2+Double(lambdau[KK,k,xyz,2])^2)/SQRT(Imaginary(lambdau[KK,k,xyz,0])^2+Imaginary(lambdau[KK,k,xyz,1])^2+Imaginary(lambdau[KK,k,xyz,2])^2))
;ELLIPTICITY CALCULATION
uppere=Imaginary(lambdau[KK,k,xyz,0])*Double(lambdau[KK,k,xyz,0])+Imaginary(lambdau[KK,k,xyz,1])*Double(lambdau[KK,k,xyz,1])
lowere=-Imaginary(lambdau[KK,k,xyz,0])^2+Double(lambdau[KK,k,xyz,0])^2-Imaginary(lambdau[KK,k,xyz,1])^2+Double(lambdau[KK,k,xyz,1])^2
if uppere GT 0 then gammarot[KK,k]=ATAN(uppere,lowere) else gammarot[KK,k]=!DPI+!DPI+ATAN(uppere,lowere)
lam=lambdau[KK,k,xyz,0:1]
lambdaurot[KK,k,*]=exp(complex(0,-0.5*gammarot[KK,k]))*lam[*]
ellip[KK,k,xyz]=Sqrt(Imaginary(lambdaurot[KK,k,0])^2+Imaginary(lambdaurot[KK,k,1])^2)/Sqrt(Double(lambdaurot[KK,k,0])^2+Double(lambdaurot[KK,k,1])^2)
ellip[KK,k,xyz]=-ellip[KK,k,xyz]*(Imaginary(ematspec[KK,k,0,1])*sin(waveangle[KK,k]))/abs(Imaginary(ematspec[KK,k,0,1])*sin(waveangle[KK,k]))
endfor
endfor
binwidth=samp_freq/nopfft
;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[KK,1:nopfft/2-2]=1/W*2*trmatspec[KK,1:nopfft/2-2]/binwidth
powspec[KK,0]=1/W*trmatspec[KK,0]/binwidth
powspec[KK,nopfft/2-1]=1/W*trmatspec[KK,nopfft/2-1]/binwidth
; added 10Sep2013 jhl
pspecx[KK,1:nopfft/2-2]=1/W*2*xrmatspec[KK,1:nopfft/2-2]/binwidth
pspecx[KK,0]=1/W*xrmatspec[KK,0]/binwidth
pspecx[KK,nopfft/2-1]=1/W*xrmatspec[KK,nopfft/2-1]/binwidth
pspecy[KK,1:nopfft/2-2]=1/W*2*yrmatspec[KK,1:nopfft/2-2]/binwidth
pspecy[KK,0]=1/W*yrmatspec[KK,0]/binwidth
pspecy[KK,nopfft/2-1]=1/W*yrmatspec[KK,nopfft/2-1]/binwidth
pspecz[KK,1:nopfft/2-2]=1/W*2*zrmatspec[KK,1:nopfft/2-2]/binwidth
pspecz[KK,0]=1/W*zrmatspec[KK,0]/binwidth
pspecz[KK,nopfft/2-1]=1/W*zrmatspec[KK,nopfft/2-1]/binwidth
pspec3[KK,*,0]=pspecx[KK, *]
pspec3[KK,*,1]=pspecy[KK, *]
pspec3[KK,*,2]=pspecz[KK, *]
ind0_fft = ind0_fft+ steplength
KK_batch_stop = KK
KK = KK+1L
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[KK_batch_start:KK_batch_stop]=ct[ind0]+ABS(nopfft/2)/samp_freq+findgen(nbp_fft_batches[batch])*steplength/samp_freq
; avoid a crash when we hit the end of timeline
if KK eq n_elements(timeline) then continue
; time tag of the leap frog between batch
timeline[KK_batch_stop+1]=ct[ind0]+ABS(nopfft/2)/samp_freq+(nbp_fft_batches[batch]+1)*steplength/samp_freq
KK = KK+1L
Endif Else Begin
binwidth=samp_freq/nopfft
dprint, 'Fourier Transform is not possible'
;print, 'Nbp = ', nbp_fft
dprint, 'Ngood = ', ngood
dprint, 'Required number of points for FFT = ', nopfft
if KK ge n_elements(timeline) then continue
timeline[KK]=ct[ind0]+ABS(nopfft/2)/samp_freq+steplength/samp_freq
powspec[KK,1:nopfft/2-2]=!values.d_nan
powspec[KK,0]=!values.d_nan
powspec[KK,nopfft/2-1]=!values.d_nan
KK = KK+1L
Endelse
ind0 = ind1_ref+1L
endfor ; loop end on batches
freqline=binwidth*findgen(nopfft/2)
; make sure there aren't any missing data points at the end of the output
wherezero = where(timeline eq 0, zerocount)
if zerocount ne 0 then begin
timeline[wherezero] = !values.d_nan
powspec[wherezero, *] = !values.d_nan
if size(elliptict,/type) eq 5 then elliptict[wherezero, *] = !values.d_nan else elliptict = make_array(n_elements(powspec[*,0]),n_elements(powspec[0,*]),value=!values.d_nan)
if size(helict,/type) eq 5 then helict[wherezero, *] = !values.d_nan else helict = make_array(n_elements(powspec[*,0]),n_elements(powspec[0,*]),value=!values.d_nan)
pspecx[wherezero, *] = !values.d_nan
pspecy[wherezero, *] = !values.d_nan
pspecz[wherezero, *] = !values.d_nan
pspec3[wherezero, *, *] = !values.d_nan
endif
return
end