;+
;  PROCEDURE :  rbsp_poynting_flux
;
;  PURPOSE  : calculates Poynting flux (ergs/s/cm^2). Separates into field-aligned and perp components
;			  Returns results as tplot variables. Includes
;			  values with and w/o spinaxis component
;
;
;  REQUIRES:  tplot library
;
;  KEYWORDS:
;       Bw -> tplot name of the [n,3] magnetic field
;             waveform (nT) in MGSE
;			  Ew -> tplot name of the [n,3] electric field
;			        waveform (mV/m) in MGSE
;			  Tshort, Tlong -> short and long period of waveform to use.
;
;       Bo -> (optional keyword) array of DC magnetic field directions in MGSE.
;					Use, for example, if Bw is from AC-coupled data.
;         If not included then Bw is downsampled and used as background field
;
;       method1 -> uses smoothing instead of a bandpass_filter which
;                  results in a sharper freq rolloff than smoothing method.
;                  Both methods give very similar results for test chorus waves
;       method2 (default -- seems to work much better for higher freqs) ->
;                  The waveform data are first downsampled to 1/Tshort to avoid the wasted
;			             computation of having to run program with data at unnecessarily high
;			             sample rate. The waveform is then run through
;			             bandpass_filter to the frequency range flow=1/Tlong...fhigh=1/Tshort
;
;   NOTES:     DO NOT INPUT DATA WITH SIGNIFICANT GAPS
;
;********************************************************************
;             Tested on Van Allen Probes chorus from EFW's B1
;             data on 2014-08-27 at ~07:42 on probe A. A is at +20
;             mlat and the Pflux indicates propagation away from eq
;             with magnitude values consistent with those in Li et
;             al., 2013 and Santolik et al., 2010
;********************************************************************
;			 Poynting flux coord system
;   		 	P1mgse = Bmgse x xhat_mgse  (xhat_mgse is spin axis component)
;				  P2mgse = Bmgse x P1mgse
;  		   	P3mgse = Bmgse
;
;
;
;			 The output tplot variables are:
;
;			 	These three output variables contain a mix of spin axis and spin plane components:
;		 		  pflux_perp1  -> Poynting flux in perp1 direction
;			 		pflux_perp2  -> Poynting flux in perp2 direction
; 			 	pflux_para   -> Poynting flux along Bo
;
;			 	These partial Poynting flux calculations contain only spin plane Ew.
;			 		pflux_nospinaxis_perp
;			 		pflux_nospinaxis_para
;
;       All of the above variables projected
;       to the ionosphere
;
;
;   CREATED:  11/28/2012
;   CREATED BY:  Aaron W. Breneman
;    LAST MODIFIED:  MM/DD/YYYY   v1.0.0
;    MODIFIED BY:
;
;-


pro rbsp_poynting_flux,Bw,Ew,Tshort,Tlong,Bo=Bo,method1=method1


  get_data,Bw,data=Bw_test
  get_data,Ew,data=Ew_test

  if ~is_struct(Bw_test) or ~is_struct(Ew_test) then begin
    print,'**************************************************'
    print,'ONE OR BOTH OF THE TPLOT VARIABLES CONTAIN NO DATA'
    print,'....RETURNING....'
    print,'**************************************************'
    return
  endif

  if is_struct(Bw_test) then begin

  ;Get DC magnetic field and use to define P1,P2,P3 directions
  if ~keyword_set(Bo) then begin
    rbsp_downsample,Bw,suffix='_DC',1/40.
    Bdc = Bw + '_DC'
  endif else begin
    tinterpol_mxn,Bo,Bw,newname='Mag_mgse_DC',/spline
    Bdc = 'Mag_mgse_DC'
  endelse


  ;Interpolate to get MagDC and Ew data to be on the same times as the Bw data
  get_data,Bw,data=goo
  times = goo.x
  tinterpol_mxn,Ew,times,/spline
  Ew = Ew + '_interp'
  tinterpol_mxn,Bdc,times,/spline
  Bdc = Bdc + '_interp'


  ;Define new coordinate system
  nelem = n_elements(times)

  ;P1 unit vector
  p1mgse = double([[replicate(0,nelem)],$
  [replicate(0,nelem)],$
  [replicate(0,nelem)]])
  get_data,Bdc,data=Bmgse_dc
  for i=0L,nelem-1 do p1mgse[i,*] = crossp(Bmgse_dc.y[i,*],[1.,0.,0.])
  ;normalize p1gse
  p1mag = fltarr(nelem)
  for i=0L,nelem-1 do p1mag[i] = sqrt(p1mgse[i,0]^2 + p1mgse[i,1]^2 + p1mgse[i,2]^2)
  for i=0L,nelem-1 do p1mgse[i,*] = p1mgse[i,*]/p1mag[i]


  ;P2 unit vector
  p2mgse = p1mgse
  for i=0L,nelem-1 do p2mgse[i,*] = crossp(Bmgse_dc.y[i,*],p1mgse[i,*])
  ;normalize p2mgse
  p2mag = fltarr(nelem)
  for i=0L,nelem-1 do p2mag[i] = sqrt(p2mgse[i,0]^2 + p2mgse[i,1]^2 + p2mgse[i,2]^2)
  for i=0L,nelem-1 do p2mgse[i,*] = p2mgse[i,*]/p2mag[i]



;*********************************************
;Test to make sure unit vectors are orthogonal
;*********************************************

;for i=0,3000 do print,acos(total(p1mgse[i,*]*p2mgse[i,*])/(p1mag[i]*p2mag[i]))/!dtor   ;perp!
;for i=0,3000 do print,acos(total(p1mgse[i,*]*Bmgse.y[i,*])/(p1mag[i]*Bmag[i]))/!dtor   ;perp!
;for i=0,3000 do print,acos(total(p2mgse[i,*]*Bmgse.y[i,*])/(p2mag[i]*Bmag[i]))/!dtor   ;perp!

;**************************************************
;**************************************************


  ;----------------------------------------------------------------------------
  ;Now we've defined our Poynting flux unit vectors as P1mgse,P2mgse,
  ;P3mgse=Bmgse_uvec. Project the Ew, Bw and Bo data into these three directions
  ;----------------------------------------------------------------------------

  ;Background magnetic field in Pflux coord
  Bmag_dc = sqrt(Bmgse_dc.y[*,0]^2 + Bmgse_dc.y[*,1]^2 + Bmgse_dc.y[*,2]^2)
  Bmgse_dc_uvec = Bmgse_dc.y
  Bmgse_dc_uvec[*,0] = Bmgse_dc.y[*,0]/Bmag_dc
  Bmgse_dc_uvec[*,1] = Bmgse_dc.y[*,1]/Bmag_dc
  Bmgse_dc_uvec[*,2] = Bmgse_dc.y[*,2]/Bmag_dc

  P3mgse = Bmgse_dc_uvec


  get_data,Ew,data=Emgse
  get_data,Bw,data=Bmgse

  Emgse = Emgse.y
  Bmgse = Bmgse.y


  Ep1 = fltarr(nelem) & Ep2 = Ep1 & Ep3 = Ep1
  for i=0L,nelem-1 do Ep1[i] = total(reform(Emgse[i,*])*reform(P1mgse[i,*]))
  for i=0L,nelem-1 do Ep2[i] = total(reform(Emgse[i,*])*reform(P2mgse[i,*]))
  for i=0L,nelem-1 do Ep3[i] = total(reform(Emgse[i,*])*reform(P3mgse[i,*]))
  Ep = [[Ep1],[Ep2],[Ep3]]


  Bp1 = fltarr(nelem) & Bp2 = Bp1 & Bp3 = Bp1
  for i=0L,nelem-1 do Bp1[i] = total(reform(Bmgse[i,*])*reform(P1mgse[i,*]))
  for i=0L,nelem-1 do Bp2[i] = total(reform(Bmgse[i,*])*reform(P2mgse[i,*]))
  for i=0L,nelem-1 do Bp3[i] = total(reform(Bmgse[i,*])*reform(P3mgse[i,*]))
  Bp = [[Bp1],[Bp2],[Bp3]]


  ;Determine data rate
  goo = rbsp_sample_rate(times,out_med_avg=medavg)
  rate = medavg[0]


;--------------------------------------------------
;Method 1 - use IDL's smooth function
;--------------------------------------------------

  if keyword_set(method1) then begin


    ;Find lower and upper periods for smoothing
    detren = floor(Tlong * rate)
    smoo = floor(Tshort * rate)

    dE1=Ep1-smooth(Ep1,detren,/nan)
    dE2=Ep2-smooth(Ep2,detren,/nan)
    dE3=Ep3-smooth(Ep3,detren,/nan)
    dB1=Bp1-smooth(Bp1,detren,/nan)
    dB2=Bp2-smooth(Bp2,detren,/nan)
    dB3=Bp3-smooth(Bp3,detren,/nan)

    fE1=smooth(dE1,smoo,/nan)
    fE2=smooth(dE2,smoo,/nan)
    fE3=smooth(dE3,smoo,/nan)
    fB1=smooth(dB1,smoo,/nan)
    fB2=smooth(dB2,smoo,/nan)
    fB3=smooth(dB3,smoo,/nan)

    B1bg=smooth((Bmgse_dc_uvec[*,0]),smoo,/nan)
    B2bg=smooth((Bmgse_dc_uvec[*,1]),smoo,/nan)
    B3bg=smooth((Bmgse_dc_uvec[*,2]),smoo,/nan)

    fE1 = fE1/1000. & fE2 = fE2/1000. & fE3 = fE3/1000.
    fB1 = fB1/1d9 & fB2 = fB2/1d9 & fB3 = fB3/1d9


  endif else begin

  ;--------------------------------------------------
  ;Method 2 - use bandpass filter
  ;--------------------------------------------------

    ;----------------------------------------------------
    ;Downsample both the Bw and Ew based on Tshort.
    ;No need to have the cadence at a higher rate than 1/Tshort
    ;----------------------------------------------------

    sr = 2/Tshort
    nyquist = sr/2.

    rbsp_downsample,[Bw,Ew],suffix='_DS_tmp',sr

    Bw = Bw +  '_DS_tmp'
    Ew = Ew +  '_DS_tmp'


    ;--------------------------------------------------
    ;At this point Ep, Bp and Bdc have a sample rate of
    ;2/Tshort and are sampled at the same times.
    ;We want to bandpass so that the lowest possible
    ;frequency is 1/Tlong Samples/sec
    ;--------------------------------------------------


    ;Define frequencies as a fraction of Nyquist
    flow = (1/Tlong)/nyquist
    fhigh = (1/Tshort)/nyquist


    ;Zero-pad these arrays to speed up FFT
    fac = 1
    nelem = n_elements(Ep[*,0])
    while 2L^fac lt n_elements(Ep[*,0]) do fac++

    addarr = fltarr(2L^fac - nelem) ;array of zeros
    Ep2 = [Ep,[[addarr],[addarr],[addarr]]]
    Bp2 = [Bp,[[addarr],[addarr],[addarr]]]

    Epf = BANDPASS_FILTER(Ep2,flow,fhigh) ;,/gaussian)
    Bpf = BANDPASS_FILTER(Bp2,flow,fhigh) ;,/gaussian)

    smoo = floor(Tshort * nyquist)
    B1bg=smooth((Bmgse_dc_uvec[*,0]),smoo,/nan)
    B2bg=smooth((Bmgse_dc_uvec[*,1]),smoo,/nan)
    B3bg=smooth((Bmgse_dc_uvec[*,2]),smoo,/nan)


    ;Remove the padded zeros
    fE = Epf[0:nelem-1,*]/1000.   ;V/m
    fB = Bpf[0:nelem-1,*]/1d9     ;Tesla

    fE1 = fE[*,0] & fE2 = fE[*,1] & fE3 = fE[*,2]
    fB1 = fB[*,0] & fB2 = fB[*,1] & fB3 = fB[*,2]

  endelse



;--------------------------------------------------
;Calculate Poynting flux
;--------------------------------------------------

  muo = 4d0*!DPI*1d-7        ; -Permeability of free space (N/A^2)

  ;J/m^2/s
  S1=(fE2*fB3-fE3*fB2)/muo
  S2=(fE3*fB1-fE1*fB3)/muo
  S3=(fE1*fB2-fE2*fB1)/muo

  ;No calculate Pflux for the versions without the "spinaxis" component
  ;Sp below stands for spinplane. Note that Sp1-hat is a different direction than p1-hat
  Sp1 = (fE1*fB3)/muo  ;remaining perp (to Bo) component while ignoring less reliable direction (in -1*p2-hat direction)
  Sp2 = (fE1*fB2)/muo  ;projected onto parallel (Bo) direction while ignoring less reliable direction (in p3-hat direction)


  ;erg/s/cm2
  S1 = S1*(1d7/1d4) & S2 = S2*(1d7/1d4) & S3 = S3*(1d7/1d4)
  Sp1 = Sp1*(1d7/1d4) & Sp2 = Sp2*(1d7/1d4)

  Svec_pfluxcoord = [[S1],[S2],[S3]]
  S1tmp = S1
  S1tmp[*] = 0.
  Svec_nospinaxis_pfluxcoord = [[S1tmp],[-1*Sp1],[Sp2]]


;--------------------------------------------------
;Find angle b/t Poynting flux and Bo
;--------------------------------------------------

  Bbkgnd = [[B1bg],[B2bg],[B3bg]]

  S_mag = sqrt(S1^2 + S2^2 + S3^2)

  ;Bo defined to be along the third component
  angle_SB = acos((S3*Bbkgnd[*,2])/S_mag)/!dtor
  store_data,'angle_pflux_Bo',data={x:times,y:angle_SB}


  ;Smooth the wave normal angle calculation
  stime = 100.*(1/rate)
  if stime lt (times[n_elements(times)-1]-times[0]) then rbsp_detrend,'angle_pflux_Bo',stime

  tplot,[Bw,Ew,'pflux_para','pflux_perp1','pflux_perp2','pflux_Ew','pflux_Bw',$
  'angle_pflux_Bo','angle_pflux_Bo_smoothed','Mag_mgse_DC_interp']




  ;Define the P unit vectors in terms of input coord system (xhat,yhat,zhat)
  Svecx = fltarr(n_elements(S1)) & Svecy = Svecx & Svecz = Svecx
  for i=0L,n_elements(S1)-1 do Svecx[i] = S1[i]*p1mgse[i,0] + S2[i]*p2mgse[i,0] + S3[i]*p3mgse[i,0]
  for i=0L,n_elements(S1)-1 do Svecy[i] = S1[i]*p1mgse[i,1] + S2[i]*p2mgse[i,1] + S3[i]*p3mgse[i,1]
  for i=0L,n_elements(S1)-1 do Svecz[i] = S1[i]*p1mgse[i,2] + S2[i]*p2mgse[i,2] + S3[i]*p3mgse[i,2]
  Svec_mgse = [[Svecx],[Svecy],[Svecz]]

  ;Same as above but with the "nospinaxis" version. From above, Sp1 is in the
  ;-1*p2-hat direction and Sp2 is in the p3-hat direction
  Svecx2 = fltarr(n_elements(Sp1)) & Svecy2 = Svecx2 & Svecz2 = Svecx2
  for i=0L,n_elements(S1)-1 do Svecx2[i] = Sp1[i]*(-1)*p2mgse[i,0] + Sp2[i]*p3mgse[i,0]
  for i=0L,n_elements(S1)-1 do Svecy2[i] = Sp1[i]*(-1)*p2mgse[i,1] + Sp2[i]*p3mgse[i,1]
  for i=0L,n_elements(S1)-1 do Svecz2[i] = Sp1[i]*(-1)*p2mgse[i,2] + Sp2[i]*p3mgse[i,2]
  Svec_mgse2 = [[Svecx2],[Svecy2],[Svecz2]]






;------------------------------------
;Estimate mapped Poynting flux
;From flux tube conservation B1*A1 = B2*A2  (A=cross sectional area of flux tube)
;B2/B1 = A1/A2
;P = EB/A ~ 1/A
;P2/P1 = A1/A2 = B2/B1
;Assume an ionospheric magnetic field of 45000 nT at 100km. This value shouldn't change too much
;and is good enough for a rough mapping estimate.
;------------------------------------


  S1_ion = 45000d * S1/Bmag_dc
  S2_ion = 45000d * S2/Bmag_dc
  S3_ion = 45000d * S3/Bmag_dc
  Sp1_ion = 45000d * Sp1/Bmag_dc
  Sp2_ion = 45000d * Sp2/Bmag_dc


;--------------------------------------------------
;Store all as tplot variables
;--------------------------------------------------

  store_data,'pflux_para',data={x:times,y:S3}
  store_data,'pflux_perp1',data={x:times,y:S1}
  store_data,'pflux_perp2',data={x:times,y:S2}
  store_data,'pflux_nospinaxis_para',data={x:times,y:Sp2}
  store_data,'pflux_nospinaxis_perp',data={x:times,y:-1*Sp1}
  store_data,'pflux_Ew',data={x:times,y:1000.*[[fE1],[fE2],[fE3]]} ;change back to mV/m
  store_data,'pflux_Bw',data={x:times,y:1d9*[[fB1],[fB2],[fB3]]}   ;change back to nT
  store_data,'pflux_para_iono',data={x:times,y:S3_ion}
  store_data,'pflux_perp1_iono',data={x:times,y:S1_ion}
  store_data,'pflux_perp2_iono',data={x:times,y:S2_ion}
  store_data,'pflux_nospinaxis_para_iono',data={x:times,y:Sp2_ion}
  store_data,'pflux_nospinaxis_perp_iono',data={x:times,y:-1*Sp1_ion}
  store_data,'angle_pflux_Bo',data={x:times,y:angle_SB}
  ;store_data,'pflux_para_mgse',data={x:times,y:pflux_para_mgse}
  ;store_data,'pflux_perp1_mgse',data={x:times,y:pflux_perp1_mgse}
  ;store_data,'pflux_perp2_mgse',data={x:times,y:pflux_perp2_mgse}
  store_data,'pflux_pfluxcoord',data={x:times,y:Svec_pfluxcoord}
  store_data,'pflux_nospinaxis_pfluxcoord',data={x:times,y:Svec_nospinaxis_pfluxcoord}
  store_data,'pflux_mgse',data={x:times,y:Svec_mgse}
  store_data,'pflux_nospinaxis_mgse',data={x:times,y:Svec_mgse2}




  options,'pflux_pfluxcoord','ytitle','Poynting flux in P1,P2,P3 coord system!C[erg/cm^2/s]'
  options,'pflux_nospinaxis_pfluxcoord','ytitle','Poynting flux in P1,P2,P3 coord system!Ccalculated without spinaxis component!C[erg/cm^2/s]'
  options,'pflux_mgse','ytitle','Poynting flux in mgse coord system!C[erg/cm^2/s]'
  options,'pflux_nospinaxis_mgse','ytitle','Poynting flux in mgse coord system!Ccalculated without spinaxis component!C[erg/cm^2/s]'
  options,'pflux_Ew','ytitle','Ew!Cpflux coord!C[mV/m]'
  options,'pflux_Bw','ytitle','Bw!Cpflux coord!C[nT]'
  options,'pflux_Ew','ytitle','Ew!Cpflux!C[mV/m]'
  options,'pflux_Bw','ytitle','Bw!Cpflux!C[nT]'
  options,'pflux_nospinaxis_perp','ytitle','pflux!Ccomponent!Cperp to Bo!C[erg/cm^2/s]'
  options,'pflux_nospinaxis_para','ytitle','pflux!Ccomponent!Cpara to Bo!C[erg/cm^2/s]'
  options,'pflux_perp1','ytitle','pflux!Cperp1 to Bo!C[erg/cm^2/s]'
  options,'pflux_perp2','ytitle','pflux!Cperp2 to Bo!C[erg/cm^2/s]'
  options,'pflux_para','ytitle','pflux!Cparallel to Bo!C[erg/cm^2/s]'
  options,'pflux_nospinaxis_perp_iono','ytitle','pflux!Cmapped!Cto ionosphere!Cperp to Bo!C[erg/cm^2/s]'
  options,'pflux_nospinaxis_para_iono','ytitle','pflux!Cmapped!Cto ionosphere!Cpara to Bo!C[erg/cm^2/s]'
  options,'pflux_perp1_iono','ytitle','pflux!Cmapped!Cperp1 to Bo!C[erg/cm^2/s]'
  options,'pflux_perp2_iono','ytitle','pflux!Cmapped!Cperp2 to Bo!C[erg/cm^2/s]'
  options,'pflux_para_iono','ytitle','pflux!Cmapped!Cparallel to Bo!C[erg/cm^2/s]'
  options,'angle_pflux_Bo','ytitle','Angle (deg) b/t!CBo and Pflux'
  options,'pflux_nospinaxis_perp','labels','No spin axis!C comp'
  options,'pflux_nospinaxis_para','labels','No spin axis!C comp!C+ along Bo'
  options,'pflux_nospinaxis_perp_iono','labels','No spin axis!C comp'
  options,'pflux_nospinaxis_para_iono','labels','No spin axis!C comp!C+ along Bo'
  options,'pflux_Ew','labels','Red=parallel!C  to Bo'
  options,'pflux_Bw','labels','Red=parallel!C  to Bo'
  options,'pflux_nospinaxis_para','colors',2
  options,'pflux_nospinaxis_perp','colors',1
  options,'pflux_nospinaxis_para_iono','colors',2
  options,'pflux_nospinaxis_perp_iono','colors',1

;  ylim,['pflux_perp1','pflux_perp2','pflux_para'],-0.2,0.2
;  ylim,['pflux_nospinaxis_perp','pflux_nospinaxis_para'],-0.2,0.2
;  ylim,['pflux_nospinaxis_perp_iono','pflux_nospinaxis_para_iono'],-10,10
;  ylim,['pflux_perp1_iono','pflux_perp2_iono','pflux_para_iono'],0,0



;;Define the Poynting flux perp and parallel values in terms of MGSE coord
;get_data,'pflux_para',data=ppara
;get_data,'pflux_perp1',data=pperp1
;get_data,'pflux_perp2',data=pperp2
;pflux_para_mgse = [[ppara.y * Bmgse_dc.y[*,0]],[ppara.y * Bmgse_dc.y[*,1]],[ppara.y * Bmgse_dc.y[*,2]]]
;pflux_perp1_mgse = [[pperp1.y * p1mgse[*,0]],[pperp1.y * p1mgse[*,1]],[pperp1.y * p1mgse[*,2]]]
;pflux_perp2_mgse = [[pperp2.y * p2mgse[*,0]],[pperp2.y * p2mgse[*,1]],[pperp2.y * p2mgse[*,2]]]

;pf1mag = sqrt(pflux_para_mgse[*,0]^2 + pflux_para_mgse[*,1]^2 + pflux_para_mgse[*,2]^2)
;pf2mag = sqrt(pflux_perp1_mgse[*,0]^2 + pflux_perp1_mgse[*,1]^2 + pflux_perp1_mgse[*,2]^2)
;pf3mag = sqrt(pflux_perp2_mgse[*,0]^2 + pflux_perp2_mgse[*,1]^2 + pflux_perp2_mgse[*,2]^2)

;for i=0,300 do print,acos(total(pflux_para_mgse[i,*]*pflux_perp1_mgse[i,*])/(pf1mag[i]*pf2mag[i]))/!dtor   ;perp!
;for i=0,300 do print,acos(total(pflux_para_mgse[i,*]*pflux_perp2_mgse[i,*])/(pf1mag[i]*pf3mag[i]))/!dtor   ;perp!
;for i=0,300 do print,acos(total(pflux_perp1_mgse[i,*]*pflux_perp2_mgse[i,*])/(pf2mag[i]*pf3mag[i]))/!dtor   ;perp!


;options,'pflux_para_mgse','ytitle','pflux para!Cin MGSE'
;options,'pflux_perp1_mgse','ytitle','pflux perp1!Cin MGSE'
;options,'pflux_perp2_mgse','ytitle','pflux perp2!Cin MGSE'


;ylim,['pflux_para','pflux_nospinaxis_para','pflux_perp1','pflux_nospinaxis_perp1','pflux_perp2','pflux_nospinaxis_perp2'],-2d-5,1.5d-4
;ylim,['pflux_perp1','pflux_nospinaxis_perp1','pflux_perp2','pflux_nospinaxis_perp2'],-1d-4,1d-4
;tplot,['pflux_para','pflux_nospinaxis_para','pflux_Ew','pflux_Bw']
;tplot,['pflux_perp1','pflux_nospinaxis_perp1','pflux_perp2','pflux_nospinaxis_perp2']
;tplot,['pflux_para','pflux_perp1','pflux_perp2']
;tplot,['pflux_para','pflux_perp1','pflux_perp2']+'_iono'


  endif
end
