thm_crib_poynting_flux.pro (Documentation for /home/davin/idl/socware/)



; This crib sheet shows how to correct the THEMIS-EFI high-frequency
; data to compensate for the transfer function, and then shows how to calculate
; Poynting flux from the EFI and SCM data.
; 2012-05-23, jmm, changed input to have user prompted for test case.

; Set the time and satellite for the event
start_input:
print, 'Choose test case:'
print, 'Input 1 for Whistler sampled at 8 kHz, THEMIS D, 2008-06-05 12:48:47 to 2008-06-05 12:48:54'
print, 'Input 2 for Oblique whistler sampled at 16 kHz, THEMIS D, 2010-01-01 00:10:02 to 2010-01-01 00:10:07.9'
print, 'Input 3 for Whistler sampled at 16 kHz, THEMIS E, 2009-07-26 12:26:40 to 2009-07-26 12:27:00'
read, tc
case tc of
  1: begin ; Whistler sampled at 8 kHz
       timespan,'2008-06-05'
       probe='d'
       interval=[str2time('2008-06-05 12:48:47'),str2time('2008-06-05 12:48:54')]
       trange=interval+[2.5,-2.5]
       power_cutoff=0.1 ; lowest power shown on 1D timeseries plots of phase, etc.
       ave_phase_yrange=[0,180]
     end
  2: begin ; Oblique whistler sampled at 16 kHz.
       timespan,'2010-01-01'
       probe='d'
       interval=[str2time('2010-01-01 00:10:02'),str2time('2010-01-01 00:10:07.9')]
       trange=interval+[3.0,-0.5]
       power_cutoff=0.1 ; lowest power shown on 1D timeseries plots of phase, etc.
       ave_phase_yrange=[180,360]
     end
  3: begin ; Whistler sampled at 16 kHz.
       timespan,'2009-07-26'
       probe='e'
       interval=[str2time('2009-07-26 12:26:40'),str2time('2009-07-26 12:27:00')]
       trange=interval+[4.5,-13.5]
       power_cutoff=0.03 ; lowest power shown on 1D timeseries plots of phase, etc.
       ave_phase_yrange=[0,180]
   end
   else: begin
       print, 'Please input 1, 2, or 3'
       goto, start_input
   end
endcase

;===============================
; Load the data
;===============================
; Load the data
thm_load_state,probe=probe[0], /get_support_data
thm_load_efi,probe=probe[0],datatype='efw',coord='dsl',trange=interval+[-10,10]
thm_load_scm,probe=probe[0],datatype='scw',coord='dsl',trange=interval+[-10,10],nk=256,/edge_truncate,fmin=10.0,fcut=10.0,despin=0
thm_load_fgm,probe=probe[0],datatype='fgs',coord='dsl',level='l2'

; rotate to FAC
tsmooth2, 'th'+probe[0]+'_fgs_dsl', 5, newname = 'smooth_B'
thm_fac_matrix_make, 'smooth_B', other_dim='Zdsl', newname = 'fac_mat'
tvector_rotate, 'fac_mat', 'th'+probe[0]+'_efw', newname = 'th'+probe[0]+'_efw'
tvector_rotate, 'fac_mat', 'th'+probe[0]+'_scw', newname = 'th'+probe[0]+'_scw'


;===============================
; Polarization analysis
;===============================
twavpol,'th'+probe[0]+'_scw',nopfft=256,steplength=64
get_data,'th'+probe[0]+'_scw_waveangle',lim=lim,dlim=dlim,data=angle
angle.y/=!dtor
store_data,'th'+probe[0]+'_scw_waveangle',lim=lim,dlim=dlim,data=angle
options,'th'+probe[0]+'_scw_waveangle',zrange=[0,40],yrange=[0,4096],ystyle=1
options,'th'+probe[0]+'_scw_powspec',zlog=1,yrange=[0,4096],ystyle=1
options,'th'+probe[0]+'_scw_degpol',yrange=[0,4096],ystyle=1


;===============================
; Correct for EFI transfer function in time domain
; Corrections include:
;   plasma-probe response
;   anti-aliasing filter response
;   ADC interleaving timing
;   DFB digital filter response
;===============================

; Set up kernels
get_data,'th'+probe[0]+'_efw',lim=lim,dlim=dlim,data=efw
fsample=double(round(1/(efw.x[1]-efw.x[0])))
print,'Sampling frequency for EFI:',fsample
kernel_length=1024
df=fsample/double(kernel_length)
f=dindgen(kernel_length)*df
f[kernel_length/2+1:*] -= double(kernel_length)*df
thm_comp_efi_response, 'SPB', f, SPB_resp,rsheath=5d6,/complex_response
thm_comp_efi_response, 'AXB', f, AXB_resp,rsheath=5d6,/complex_response
if fsample eq 16384 then begin
  print,'Assuming AC coupling.'
  E12_resp=1/(SPB_resp*thm_eac_filter_resp(f)*thm_adc_resp('E12AC',f)*thm_dfb_dig_filter_resp(f, fsample))
  E34_resp=1/(SPB_resp*thm_eac_filter_resp(f)*thm_adc_resp('E34AC',f)*thm_dfb_dig_filter_resp(f, fsample))
  E56_resp=1/(AXB_resp*thm_eac_filter_resp(f)*thm_adc_resp('E56AC',f)*thm_dfb_dig_filter_resp(f, fsample))
endif else begin
  print,'Assuming DC coupling.'
  E12_resp=1/(SPB_resp*bessel_filter_resp(f,4096,4)*thm_adc_resp('E12DC',f)*thm_dfb_dig_filter_resp(f, fsample))
  E34_resp=1/(SPB_resp*bessel_filter_resp(f,4096,4)*thm_adc_resp('E34DC',f)*thm_dfb_dig_filter_resp(f, fsample))
  E56_resp=1/(AXB_resp*bessel_filter_resp(f,4096,4)*thm_adc_resp('E56DC',f)*thm_dfb_dig_filter_resp(f, fsample))
endelse

; Filter out frequencies above SCM Nyquist (before downsampling)
if fsample eq 16384 then begin
  E12_resp[kernel_length/4:3*kernel_length/4-1]=0
  E34_resp[kernel_length/4:3*kernel_length/4-1]=0
  E56_resp[kernel_length/4:3*kernel_length/4-1]=0
  E12_resp[0]=0
  E34_resp[0]=0
  E56_resp[0]=0
endif

; Transfer kernel into time domain: take inverse FFT and center
E12_resp=shift((fft(E12_resp,1)),kernel_length/2)/kernel_length
E34_resp=shift((fft(E34_resp,1)),kernel_length/2)/kernel_length
E56_resp=shift((fft(E56_resp,1)),kernel_length/2)/kernel_length

; Extend data to account for edge wrap
ndata=n_elements(efw.x)
efw_y=[[replicate(efw.y[0,0],kernel_length/2),efw.y[*,0],replicate(efw.y[ndata-1,0],kernel_length/2)],$
       [replicate(efw.y[0,1],kernel_length/2),efw.y[*,1],replicate(efw.y[ndata-1,1],kernel_length/2)],$
       [replicate(efw.y[0,2],kernel_length/2),efw.y[*,2],replicate(efw.y[ndata-1,2],kernel_length/2)]]

; Deconvolve transfer function
b_length = 8 * kernel_length
efw_y[*,0] = shift(blk_con(E12_resp, efw_y[*,0], b_length=b_length),-kernel_length/2)
efw_y[*,1] = shift(blk_con(E34_resp, efw_y[*,1], b_length=b_length),-kernel_length/2)
efw_y[*,2] = shift(blk_con(E56_resp, efw_y[*,2], b_length=b_length),-kernel_length/2)

; Align data to SCM (downsampling to 8 kHz if needed) and store it
get_data,'th'+probe[0]+'_scw',data=scw
efw={x:scw.x,y:[[interpol(efw_y[kernel_length/2:ndata+kernel_length/2-1,0],efw.x,scw.x)],$
                [interpol(efw_y[kernel_length/2:ndata+kernel_length/2-1,1],efw.x,scw.x)],$
                [interpol(efw_y[kernel_length/2:ndata+kernel_length/2-1,2],efw.x,scw.x)]]}
store_data,'th'+probe[0]+'_efw_corrected',lim=lim,dlim=dlim,data=efw


;===============================
; Poynting flux analysis
;===============================
nfft=128

; Get data
get_data,'th'+probe[0]+'_efw_corrected',lim=lim,dlim=dlim,data=efw
get_data,'th'+probe[0]+'_scw',lim=lim1,dlim=dlim1,data=scw

; Calculate Poynting flux (bandpass filtered, time domain)
ndata = n_elements(efw.x)       ;jmm, 22-may-2012 efw was interpolated earlier
S=dblarr(ndata,3)
E=efw.y
B=scw.y
filter=digital_filter(200./4096,1,50,nfft,/double)
for i=0,2 do E[*,i]=convol(E[*,i],filter,/center)
for i=0,2 do B[*,i]=convol(B[*,i],filter,/center)
S[*,0]= E[*,1]*B[*,2]-E[*,2]*B[*,1]
S[*,1]=-E[*,0]*B[*,2]+E[*,2]*B[*,0]
S[*,2]= E[*,0]*B[*,1]-E[*,1]*B[*,0]
S_conversion=1d-3*1d-9*1d6/(4d-7*!dpi)  ; mV->V, nT->T, W->uW, divide by mu_0
S*=S_conversion
for i=0,2 do S[*,i]=smooth(S[*,i],nfft)
store_data,'S_timeseries',data={x:efw.x,y:S},lim={ytitle:'Poynting flux 200-4096 Hz!C!C[!4l!3W/m2]'}
store_data,'S_tot1',    data={x:efw.x,y:sqrt(total(S*S,2))}

; Perform FFTs for frequency-domain calculation
nfft=128
stride=32
ndata=n_elements(efw.x)
efw_fft=dcomplexarr(long(ndata-nfft)/stride+1,nfft,3)
scw_fft=dcomplexarr(long(ndata-nfft)/stride+1,nfft,3)
win=hanning(nfft,/double)
win/=mean(win^2)  ; preserve energy
i=0L
for j=0L,ndata-nfft-1,stride do begin
  for k=0,2 do efw_fft[i,*,k]=fft(efw.y[j:j+nfft-1,k]*win)
  for k=0,2 do scw_fft[i,*,k]=fft(scw.y[j:j+nfft-1,k]*win)
  i++
endfor
t=scw.x[0]+(dindgen(i-1)*stride+nfft/2)/8192.
freq=(findgen(nfft/2)+0.5)*8192/nfft
bw=8192/nfft
efwlim={spec:1,zlog:1,ylog:0,yrange:[100,4096],ystyle:1,zrange:[1e-8,1e-4]}
scwlim={spec:1,zlog:1,ylog:0,yrange:[100,4096],ystyle:1,zrange:[1e-10,1e-6]}
store_data,'efw_fft_x',data={x:t,y:abs(efw_fft[*,0:nfft/2,0])^2/bw,v:freq},lim=efwlim
store_data,'efw_fft_y',data={x:t,y:abs(efw_fft[*,0:nfft/2,1])^2/bw,v:freq},lim=efwlim
store_data,'efw_fft_z',data={x:t,y:abs(efw_fft[*,0:nfft/2,2])^2/bw,v:freq},lim=efwlim
store_data,'scw_fft_x',data={x:t,y:abs(scw_fft[*,0:nfft/2,0])^2/bw,v:freq},lim=scwlim
store_data,'scw_fft_y',data={x:t,y:abs(scw_fft[*,0:nfft/2,1])^2/bw,v:freq},lim=scwlim
store_data,'scw_fft_z',data={x:t,y:abs(scw_fft[*,0:nfft/2,2])^2/bw,v:freq},lim=scwlim

; Calculate Poynting flux (frequency domain)
Sx= double(efw_fft[*,*,1]*conj(scw_fft[*,*,2])-efw_fft[*,*,2]*conj(scw_fft[*,*,1]))*S_conversion
Sy=-double(efw_fft[*,*,0]*conj(scw_fft[*,*,2])-efw_fft[*,*,2]*conj(scw_fft[*,*,0]))*S_conversion
Sz= double(efw_fft[*,*,0]*conj(scw_fft[*,*,1])-efw_fft[*,*,1]*conj(scw_fft[*,*,0]))*S_conversion
bw=8192/nfft
indx=where(freq ge 200)
Stot=sqrt(total(Sx[*,indx]^2+Sy[*,indx]^2+Sz[*,indx]^2,2))
zrange=max(Stot)*0.1*[-1,1]/bw
store_data,'S_x',     data={x:t,y:Sx/bw,v:freq},lim={spec:1,zlog:0,ylog:0,yrange:[100,4096],ystyle:1,zrange:zrange,ytitle:'Poynting flux S!Bx!N!C!C[!4l!3W/m!A2!N/Hz]'}
store_data,'S_y',     data={x:t,y:Sy/bw,v:freq},lim={spec:1,zlog:0,ylog:0,yrange:[100,4096],ystyle:1,zrange:zrange,ytitle:'Poynting flux S!By!N!C!C[!4l!3W/m!A2!N/Hz]'}
store_data,'S_z',     data={x:t,y:Sz/bw,v:freq},lim={spec:1,zlog:0,ylog:0,yrange:[100,4096],ystyle:1,zrange:zrange,ytitle:'Poynting flux S!Bz!N!C!C[!4l!3W/m!A2!N/Hz]'}
store_data,'S_tot2',data={x:t,y:Stot},lim={color:1}
store_data,'S_tot' ,data=['S_tot1','S_tot2'],lim={ytitle:'|S| 200-4096 Hz!C!C[!4l!3W/m2]'}


;===============================
; Plot
;===============================
window,0,ysize=900
window,1,ysize=900
tplot,['??w_*fft_?'],trange=trange,title='Calibrated data with transfer function corrected',window=0
tplot,['th'+probe[0]+'_scw_'+['powspec','degpol','waveangle'],'S_?','S_timeseries'],trange=trange,title='Poynting flux',window=1

end