;+
;#############################################
; DEPRECATED --> Please use moka_mms_pad_fpi
;#############################################
;
;Procedure:
; moka_mms_pad
;
;Purpose:
; Returns a pitch-angle-distribution from MMS FPI data (angle vs energy plot)
; as well as energy spectrum in the omni, para, perp and anti-para directions.
; One-count-level is also returned.
;
;Calling Sequence:
; structure = moka_mms_pad(bname, tname [,index] [,trange=trange] [,units=units],[,/norm],
; [,nbin=nbin], [,vname=vname] [,/structure])
;
; INPUT:
; bname: magnetic field, tplot-variable name, use burst data
; tname: FPI data, tplot-variable such as "mms?_des_dist_brst"
; index: (NOW DEPRECATED! after a struggle with apj2016_egyspec.pro)
; trange: Two element time range to constrain the requested data (See also mms_get_fpi_dist)
; nbin: number of bins in the pitch-angle direction
; vname: bulk flow velocity for frame transformation, tplot-variable name,
; vname & tname should have the same data_rate
; norm: Set this keyword for normalizing the data at each energy bin
; units: units for both the pitch-angle-distribution (pad) and energy spectrum.
; Options are 'eflux' [eV/(cm!U2!N s sr eV)] or
; 'df' [s!U3!N / km!U6!N']
; The default is 'eflux'. The return structure contains a tag "UNITS".
; pr___0: pitch angle range for the "para" spectrum, default = [0,45]
; pr__90: pitch angle range for the "perp" spectrum, default = [45,135]
; pr_180: pitch angle range for the "anti-para" spectrum, default = [135,180]
;
;Output:
; a structure containing the result
;
;Example:
; MMS> trange = '2015-11-04/'+['04:57:49','04:57:50']
; MMS> tname = 'mms3_des_dist_brst'
; MMS> bname = 'mms3_fgm_b_dmpa_brst_l2_bvec'
; MMS> vname = 'mms3_des_bulk_dbcs_brst'
; MMS> pad = moka_mms_pad(bname, tname, trange, vname=vname)
; MMS> plotxyz,pad.PA, pad.EGY, pad.DATA,/noisotropic,/ylog,zlog=1,$
; xrange=[0,180],zrange=[1e+5,1e+9],xtitle='pitch angle',ytitle='energy'
;
;History:
; Created by Mitsuo Oka on 2016-05-15
; Fixed energy bin mistake 2017-01-28
; Fixed para and anti-para mistake (thanks to R. Mistry) 2017-03-14
; Fixed eflux calculation 2017-05-12
; Added SUBTRACT_ERROR keyword 2017-10-17
; Removed unnecessary vname check 2017-10-19
;
;$LastChangedBy: egrimes $
;$LastChangedDate: 2018-04-03 15:14:57 -0700 (Tue, 03 Apr 2018) $
;$LastChangedRevision: 24992 $
;$URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/spdsoft/tags/spedas_5_0/projects/mms/particles/deprecated/moka_mms_pad.pro $
;-
FUNCTION moka_mms_pad, bname, tname, trange, units=units, nbin=nbin, vname=vname, $
norm=norm, ename=ename, pr___0 = pr___0, pr__90 = pr__90, pr_180=pr_180,$
daPara=da,daPerp=da2,oclreal=oclreal, single_time=single_time, $
subtract_bulk = subtract_bulk, subtract_error=subtract_error, $
vlm=vlm; An additional frametransformation
compile_opt idl2
;------------
; initialize
;------------
if undefined(tname) then message,'Please specify FPI data.'
if undefined(bname) then message,'Please specify mag data.'
if undefined(trange) then message, 'Please specify trange.'
if undefined(nbin) then nbin = 18 ; Number of pitch-angle bins
if undefined(da) then da = 45.
if undefined(da2) then da2= 45.
if undefined(pr___0) then pr___0 = [0., da]; para pitch-angle range
if undefined(pr__90) then pr__90 = [90.-da2, 90.+da2]; perp pitch-angle range
if undefined(pr_180) then pr_180 = [180.-da, 180.]; anti-para pitch-angle range
if ~undefined(vname) then begin
if undefined(subtract_bulk) then begin
msg = "moka_mms_pad now requires the keyword 'subtract_bulk' = 1 (in addition to specifying "
msg += "vname) when subtracting bulk velocity"
result = dialog_message(msg,/center)
return, -1
endif
endif
if size(tname,/type) eq 7 then begin; if 'tname' is string...
distime = mms_get_dist(tname,/time); start time
if (time_double(trange[1]) lt distime[0]) or (max(distime) lt time_double(trange[0])) then begin
print, 'trange=',time_string(trange,prec=4)
print, 'distime=',time_string([distime[0], max(distime)],prec=4)
print,'trange is out of '+tname+' time range.'
return, -1
endif
dist = mms_get_dist(tname,index,trange=trange,/structure,subtract_error=subtract_error,error=ename)
if (n_tags(dist) eq 0) then begin
print,'FPI data could not be extracted from the specified time period.'
return, -1
endif
spc = strmid(tname,6,1)
endif else stop
USEERR = 0
if ~undefined(ename) then begin
if size(ename,/type) eq 7 then begin; if 'ename' is string...
distErr = mms_get_dist(ename,index,trange=trange,/structure)
if n_tags(distErr) gt 0 then begin
USEERR = 1
endif else begin
print, 'WARNING: distErr not loaded properly.'
endelse
endif else stop
endif
units_lc = undefined(units) ? 'eflux' : strlowcase(units)
units_gl = strmatch(units_lc,'eflux') ? 'eV/(cm!U2!N s sr eV)' : 's!U3!N / km!U6!N'
nmax = n_elements(dist)
tr = [dist[0].TIME, dist[nmax-1].END_TIME]
dr = !dpi/180.d0
rd = 1.d0/dr
if strmatch(units_lc,'eflux') then begin
species_lc = strlowcase(dist[0].species)
;get mass of species
case species_lc of
'i': A=1;H+
'hplus': A=1;H+
'heplus': A=4;He+
'heplusplus': A=4;He++
'oplus': A=16;O+
'oplusplus': A=16;O++
'e': A=1d/1836;e-
else: message, 'Unknown species: '+species_lc
endcase
;scaling factor between df_km and flux units
flux_to_df = A^2 * 0.5447d * 1d6
;convert between km^6 and cm^6 for df_km
cm_to_km = 1d30
; See 'mms_convert_flux_units'
in = [2,-1,0]; units_in = 'df_km'
out = [0,0,0]; units_out = 'eflux'
exp = in + out
endif
;-----------------------------
; Pitch-angle Bins
;-----------------------------
kmax = nbin
pamin = 0.
pamax = 180.
dpa = (pamax-pamin)/double(kmax); bin size
wpa = pamin + findgen(kmax)*dpa + 0.5*dpa; bin center values
pa_bin = [pamin + findgen(kmax)*dpa, pamax]
;-----------------------------
; Energy Bins
;-----------------------------
wegy = dist[0].ENERGY[*,0,0]
if keyword_set(subtract_bulk) then begin
wegy = [2*wegy[0]- wegy[1], wegy]
endif
jmax = n_elements(wegy)
egy_bin = 0.5*(wegy + shift(wegy,-1))
egy_bin[jmax-1] = 2.*wegy[jmax-1]-egy_bin[jmax-2]
egy_bin0 = 2.*wegy[ 0]-egy_bin[ 0] > 0
egy_bin = [egy_bin0, egy_bin]
;----------------
; Prep
;----------------
pad = fltarr(jmax, kmax)
mmax = 4
f_dat = fltarr(jmax,mmax); Four elements for para, perp, anti-para and omni directions
f_psd = fltarr(jmax,mmax)
f_err = fltarr(jmax,mmax)
f_cnt = fltarr(jmax,mmax)
count_pad = lonarr(jmax, kmax); number of events in each bin
count_dat = lonarr(jmax, mmax)
;------------------------------------
; Magnetic field & Bulk Velocity
;------------------------------------
bnrm = double(moka_tplot_average(bname, tr,norm=1)); .... Normalized!
Vbulk = [0.d0,0.d0,0.d0]
Vbulk_para = 0.d0
Vbulk_perp = 0.d0
Vbulk_vxb = 0.d0
Vbulk_exb = 0.d0
if ~undefined(vname) then begin
get_data,vname,data=DV
endif
;----------------
; Main Loop
;----------------
iecl = 0L
iecm = 0L
for n=0,nmax-1 do begin
;----------------
;Sanitize Data.
;----------------
if USEERR then begin
moka_mms_clean_data,dist[n],output=data,units=units_lc,disterr=distErr[n]
endif else begin
moka_mms_clean_data,dist[n],output=data,units=units_lc
endelse
;------------------------
; Bulk Velocity
;------------------------
if ~undefined(vname) then begin
result = min(DV.x-data.TIME,m,/nan,/abs)
Vbulk = double(reform(DV.y[m,*]))
endif
vbpara = bnrm[0]*Vbulk[0]+bnrm[1]*Vbulk[1]+bnrm[2]*Vbulk[2]
vbperp = Vbulk - vbpara
vbperp_abs = sqrt(vbperp[0]^2+vbperp[1]^2+vbperp[2]^2)
vxb = [-Vbulk[1]*bnrm[2]+Vbulk[2]*bnrm[1],-Vbulk[2]*bnrm[0]+Vbulk[0]*bnrm[2],-Vbulk[0]*bnrm[1]+Vbulk[1]*bnrm[0]]
vxbabs = sqrt(vxb[0]^2+vxb[1]^2+vxb[2]^2)
vxbnrm = vxb/vxbabs
exb = [vxb[1]*bnrm[2]-vxb[2]*bnrm[1],vxb[2]*bnrm[0]-vxb[0]*bnrm[2],vxb[0]*bnrm[1]-vxb[1]*bnrm[0]]
exbabs = sqrt(exb[0]^2+exb[1]^2+exb[2]^2)
exbnrm = exb/exbabs
Vbulk_para += vbpara
Vbulk_perp += vbperp_abs
Vbulk_vxb += vxbnrm[0]*Vbulk[0]+vxbnrm[1]*Vbulk[1]+vxbnrm[2]*Vbulk[2]
Vbulk_exb += exbnrm[0]*Vbulk[0]+exbnrm[1]*Vbulk[1]+exbnrm[2]*Vbulk[2]
;------------------------------------
; Particle Velocities & Pitch Angles
;------------------------------------
; Spherical to Cartesian
erest = data.mass * !const.c^2 / 1e6; convert mass from eV/(km/s)^2 to eV
vabs = !const.c * sqrt( 1 - 1/((data.ENERGY/erest + 1)^2) ) / 1000.;velocity in km/s
sphere_to_cart, vabs, data.theta, data.phi, vx, vy, vz
; Frame transformation
; if undefined(vlmpot) then begin
if keyword_set(subtract_bulk) then begin; Plasma rest-frame
vx -= Vbulk[0]
vy -= Vbulk[1]
vz -= Vbulk[2]
endif
; endif else begin; Another frame (under development)
; bvec = double(moka_tplot_average(bname, tr,norm=0))
; moka_mms_vlm, vx, vy, vz, Vbulk, bvec, vlm
; endelse
; Pitch angles
dp = (bnrm[0]*vx + bnrm[1]*vy + bnrm[2]*vz)/sqrt(vx^2+vy^2+vz^2)
idx = where(dp gt 1.d0, ct) & if ct gt 0 then dp[idx] = 1.d0
idx = where(dp lt -1.d0, ct) & if ct gt 0 then dp[idx] = -1.d0
pa = rd*acos(dp)
; Cartesian to Spherical
cart_to_sphere, vx, vy, vz, vnew, theta, phi, /ph_0_360
data.energy = erest*(1.d0/sqrt(1.d0-(vnew*1000.d0/!const.c)^2)-1.d0); eV
data.phi = phi
data.theta = theta
data.pa = pa
;------------------------
; DISTRIBUTE
;------------------------
imax = n_elements(data.data_dat)
for i=0,imax-1 do begin; for each particle
; Find energy bin
result = min(egy_bin-data.energy[i],j,/abs)
if egy_bin[j] gt data.energy[i] then j -= 1
if j eq jmax then j -= 1
; Find pitch-angle bin
result = min(pa_bin-data.pa[i],k,/abs)
if pa_bin[k] gt data.pa[i] then k -= 1
if k eq kmax then k -= 1
if j ge 0 then begin
;-----------------------
; Find new eflux
;-----------------------
; If shifted to plasma rest-frame, 'eflux' should be re-evaluated
; from 'psd' because 'eflux' depends on the particle energy. We don't need to
; worry about this if we want the output in 'psd'.
newenergy = wegy[j]
if strmatch(units_lc,'eflux') then begin; If plasma rest-frame AND efluxs
; See 'mms_convert_flux_units'
newdat = double(data.data_psd[i]) * double(newenergy)^exp[0] * (flux_to_df^exp[1] * cm_to_km^exp[2])
newpsd = newdat
newerr = double(data.data_err[i]) * double(newenergy)^exp[0] * (flux_to_df^exp[1] * cm_to_km^exp[2])
endif else begin
newdat = data.data_dat[i]
newpsd = data.data_psd[i]
newerr = data.data_err[i]
endelse
; pitch-angle distribution
pad[j,k] += newdat
count_pad[j,k] += 1L
; energy spectrum (para, perp, anti-para)
m = -1
if (pr__90[0] le data.pa[i]) and (data.pa[i] le pr__90[1]) then begin
m = 1
endif else begin
if (pr___0[0] le data.pa[i]) and (data.pa[i] le pr___0[1]) then m=0
if (pr_180[0] le data.pa[i]) and (data.pa[i] le pr_180[1]) then m=2
endelse
if (m ge 0) and (m le 2) then begin
f_dat[j,m] += newdat
f_psd[j,m] += newpsd
f_err[j,m] += newerr
f_cnt[j,m] += data.data_cnt[i]
count_dat[j,m] += 1L
endif
; energy spectrum (omni-direction)
m = 3
f_dat[j,m] += newdat
f_psd[j,m] += newpsd
f_err[j,m] += newerr
f_cnt[j,m] += data.data_cnt[i]
count_dat[j,m] += 1L
endif else begin; if j ge 0
iecl += 1L
endelse
endfor; for each particle
iecm += imax
endfor; for n=0,nmax-1
if iecl gt 0 then begin
print, '---'
print, strtrim(string(iecl),2)+' out of '+strtrim(string(iecm),2)+' particles were not'
print,'used in the result because of the bulk-speed subtraction'
print,'and their energies have moved outside the defined energy bins.'
print,'---'
endif
pad /= float(count_pad)
f_dat /= float(count_dat)
f_psd /= float(count_dat)
f_err /= float(count_dat)
f_cnt /= float(count_dat)
vbulk_para /= float(nmax)
vbulk_perp /= float(nmax)
vbulk_vxb /= float(nmax)
vbulk_exb /= float(nmax)
idx = where(~finite(pad),ct)
if ct gt 0 then pad[idx] = 0.
;---------------
; ANGLE PADDING
;---------------
padnew = fltarr(jmax, kmax+2)
padnew[0:jmax-1,1:kmax] = pad
padnew[0:jmax-1, 0] = padnew[0:jmax-1, 1]
padnew[0:jmax-1,kmax+1] = padnew[0:jmax-1,kmax]
pad = padnew
wpa_new = [wpa[0]-dpa,wpa,wpa[kmax-1]+dpa]
wpa = wpa_new
;---------------
; NORMALIZE
;---------------
padnorm = pad
for j=0,jmax-1 do begin; for each energy
peak = max(pad[j,0:kmax+1],/nan); find the peak
if peak eq 0 then begin
padnorm[j,0:kmax+1] = 0.
endif else begin
padnorm[j,0:kmax+1] /= peak
endelse
endfor
if keyword_set(norm) then pad = padnorm
;----------------------------
; Effective one-count-level
;----------------------------
; 'f_psd' is the PSD averaged over time and angular ranges.
; 'f_cnt' is the counts averaged over time and angular ranges.
; 'count_dat is the total number of time and angular bins.
f_ocl = f_psd/(f_cnt*float(count_dat))
; Real one-count-level of the instrument in a sampling time and in an angular bin
if keyword_set(oclreal) then begin
f_ocl = f_psd/f_cnt
endif
;---------------
; OUTPUT
;---------------
return, {egy:wegy, pa:wpa, data:transpose(pad), datanorm:transpose(padnorm), $
numSlices: nmax, nbin:kmax, units:units_gl,subtract_bulk:subtract_bulk,$
egyrange:[min(wegy),max(wegy)], parange:[min(wpa),max(wpa)], $
spec___0:f_psd[*,0], spec__90:f_psd[*,1], spec_180:f_psd[*,2], spec_omn:f_psd[*,3], $
cnts___0:f_cnt[*,0], cnts__90:f_cnt[*,1], cnts_180:f_cnt[*,2], cnts_omn:f_cnt[*,3], $
oclv___0:f_ocl[*,0], oclv__90:f_ocl[*,1], oclv_180:f_ocl[*,2], oclv_omn:f_ocl[*,3], $
eror___0:f_err[*,0], eror__90:f_err[*,1], eror_180:f_err[*,2], eror_omn:f_err[*,3], $
trange:tr, vbulk_para:vbulk_para, vbulk_perp_abs:vbulk_perp, vbulk_vxb:vbulk_vxb, $
vbulk_exb:vbulk_exb, bnrm:bnrm, Vbulk:Vbulk}
END