;+ ;Procedure: ; moka_mms_pad_fpi ; ;Purpose: ; To process MMS FPI data and return ; (1) pitch-angle-distribution (angle vs energy plot) ; (2) energy spectrum in the omni, para, perp and anti-para directions. ; (3) One-count-level is also returned. ; ;Calling Sequence: ; ; (Similar to spd_slice2d. See also moka_mms_pad_fpi_crib.pro) ; ; structure = moka_mms_pad_fpi(dist [,disterr] $ ; [,time=time [,window=window | samples=samples]] ; [trange=trange] ... ) ; ; INPUT: ; DIST : A pointer to 3D data structure. ; DISTERR: A pointer to 3D data error structure ; TRANGE : Two-element time range over which data will be averaged. (string or double) ; TIME : Time at which the pad will be computed. (string or double) ; SAMPLES: Number of nearest samples to TIME to average. (int/double) ; If neither SAMPLES nor WINDOW are specified then default=1. ; WINDOW: Length in seconds from TIME over which data will be averaged. (int/double) ; CENTER_TIME: Flag denoting that TIME should be midpoint for window instead of beginning. ; ; MAG_DATA: Name of tplot variable containing magnetic field data or 3-vector. ; This will be used for pitch-angle calculation and must be in the ; same coordinates as the particle data. ; VEL_DATA: Name of tplot variable containing the bulk velocity data or 3-vector. ; This will be used for pitch-angle calculation and must be in the ; same coordinates as the particle data. ; ; nbin: number of bins in the pitch-angle direction ; ; 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_km' [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 results ; ;History: ; 2016-05-15 Created by Mitsuo Oka ; 2017-01-28 Fixed energy bin mistake ; 2017-03-14 Fixed para and anti-para mistake (thanks to R. Mistry) ; 2017-05-12 Fixed eflux calculation ; 2017-10-17 Added SUBTRACT_ERROR keyword ; 2017-10-17 Changed the interface so that it works like spd_slice2d ; ;$LastChangedBy: moka $ ;$LastChangedDate: 2017-09-30 11:03:14 -0700 (Sat, 30 Sep 2017) $ ;$LastChangedRevision: 24073 $ ;$URL: svn+ssh://ambrosia.ssl.berkeley.edu/repos/spdsoft/trunk/projects/mms/particles/moka/moka_mms_pad.pro $ ;- FUNCTION moka_mms_pad_fpi, input, input_err, $ ;------ Time options --------- time=time_in, $ window=window, $ center_time=center_time, $ trange=trange_in, $ samples=samples, $ ;----- Support Data ----- mag_data=mag_data, $ vel_data=vel_data, $ ;------ Output options ---------- nbin=nbin, $; Number of pitch-angle bins da=da, da2=da2, $ pr___0=pr___0, pr__90=pr__90, pr_180=pr_180, $ units=units,$ norm=norm, $ subtract_bulk=subtract_bulk, $ _extra = _extra compile_opt idl2 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(subtract_bulk) then subtract_bulk = 0 ;----------------- ; Constant ;----------------- invalid = 0b fail = '' dr = !dpi/180.d0 rd = 1.d0/dr ;----------------- ; CHECK ;----------------- if ~ptr_valid(input[0]) and ~is_struct(input[0]) then begin fail = 'Invalid data. Input must be pointer or structure array.' dprint, dlevel=1, fail return, invalid endif ;-------------------------------- ; Aggregate inputs ; -multiple inputs are allowed to make things easier on the user ; -pointer arrays are used as a way of supporting dissimilar ; structures, which cannot be concatenated ;------------------------------------------------------------ ;copy any structs to new pointer var, other inputs are copied and checked next switch n_params() of 2: if is_struct(input_err) then pDerr = ptr_new(input_err) else pDerr = input_err 1: if is_struct(input) then pD = ptr_new(input) else pD = input endswitch USEERR = ptr_valid(pDerr) nmax = n_elements(*pD) ;--------------------- ; Time range ;--------------------- ; CASE A: time and window ; CASE B: time and samples ; CASE C: trange ; if not Case C, make sure we have time and window/samples if undefined(trange_in) then begin if undefined(time_in) then begin fail = 'Please specifiy a time or time range over which to compute the pad. For example: '+ $ ssl_newline()+' "TIME=t, WINDOW=w" or "TRANGE=tr" or "TIME=t, SAMPLES=n"' dprint, dlevel=1, fail return, invalid endif else begin if undefined(window) && undefined(samples) then begin samples = 1 ;use single closest distribution by default endif endelse endif ; get the time range if one was specified if ~undefined(time_in) then time = time_double(time_in[0]) ; [CASE C] ; time range already provided if ~undefined(trange_in) then trange = minmax(time_double(trange_in)) ; [CASE A] ; get the time range if a time & window were specified instead if undefined(trange) && keyword_set(window) then begin if keyword_set(center_time) then begin trange = [time - window/2d, time + window/2d ] endif else begin trange = [time, time + window ] endelse endif ; [CASE B] ; if no time range or window was specified then get a time range ; from the N closest samples to the specied time ; (defaults to 1 if SAMPLES is not defined) if undefined(trange) then begin trange = spd_slice2d_nearest(pD, time, samples) endif ;-------------------------------------- ; get 'n_samples' ;-------------------------------------- ; While the user may specify 'samples, 'n_samples' here is the actual number of ; samples to be used. ; ; check that there is data in range before proceeding for i=0, n_elements(pD)-1 do begin times_ind = spd_slice2d_intrange(pD[i], trange, n=ndat); the indices of all samples in the specified trange n_samples = array_concat(ndat,n_samples) endfor n_samples = total(n_samples) if n_samples lt 1 then begin fail = 'No particle data in the time window: '+strjoin(time_string(trange),', ')+ $ '. Time samples may be at low cadence; try adjusting the time window.' dprint, dlevel=1, fail return, invalid endif dprint, dlevel=3, strtrim(n_samples,2) + ' samples in time window' if keyword_set(fail) then return, invalid nns = min(times_ind,/nan) nne = max(times_ind,/nan) ;------------------------------------------------------------ ; Check Support data ; ; Just to check if mag_data and vel_data exist or not. ; The values 'bfield' and 'vbulk' will not be used in the main loop because ; they are taken from the specified trange Instead, we use the values at each sample. ; Thus, the pitch angles are calculated at each sample. ;------------------------------------------------------------ spd_slice2d_get_support, mag_data, trange, output=bfield spd_slice2d_get_support, vel_data, trange, output=vbulk if undefined(bfield) then begin fail = 'Magnetic field data needed to calculate pitch-angles' dprint, dlevel=1, fail return, invalid endif if undefined(vbulk) then begin if keyword_set(subtract_bulk) then begin fail = 'Velocity data needed to subtract bulk velocity.' dprint, dlevel=1, fail return, invalid endif endif else begin get_data,vel_data,data=DV endelse ;------------------------------------------------------------ ; Unit ;--------------------------------------------- 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' species_lc = strlowcase((*pD)[0].species) if strmatch(units_lc,'eflux') then begin ;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 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] ;----------------------------- ; Azimuthal Bins ;----------------------------- amax = 16 azmin = 0. azmax = 180. daz = (azmax-azmin)/double(amax); bin size waz = azmin + findgen(amax)*daz + 0.5*daz; bin center values az_bin = [azmin + findgen(amax)*daz, azmax] ;----------------------------- ; Polar Bins ;----------------------------- pmax = 16 polmin = 0. polmax = 360. dpol = (polmax-polmin)/double(pmax); bin size wpol = polmin + findgen(pmax)*dpol + 0.5*dpol; bin center values pol_bin = [polmin + findgen(pmax)*dpol, polmax] ;----------------------------- ; PA (az x pol) ;----------------------------- pa_azpol = dblarr(amax, pmax) pa_azpol[*, *] = !values.d_nan ;----------------------------- ; Energy Bins ;----------------------------- wegy = (*pD)[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_avg = [0., 0., 0.] babs_avg = 0. Vbulk_avg = [0.d0,0.d0,0.d0] Vbulk_para = 0.d0 Vbulk_perp = 0.d0 Vbulk_vxb = 0.d0 Vbulk_exb = 0.d0 ;---------------- ; Main Loop ;---------------- iecl = 0L iecm = 0L for n=nns,nne do begin ;---------------- ;Sanitize Data. ;---------------- if USEERR then begin moka_mms_clean_data,(*pD)[n],output=data,units=units_lc,disterr=(*pDerr)[n] endif else begin moka_mms_clean_data,(*pD)[n],output=data,units=units_lc endelse ;------------------------ ; Magnetic field direction ;------------------------ tr = [(*pD)[n].TIME, (*pD)[n].END_TIME] bfield = spd_tplot_average(mag_data, tr) babs = sqrt(bfield[0]^2+bfield[1]^2+bfield[2]^2) bnrm = bfield/babs bnrm_avg += bnrm babs_avg += babs ;------------------------ ; Bulk Velocity ;------------------------ if keyword_set(subtract_bulk) then begin result = min(DV.x-data.TIME,m,/nan,/abs) Vbulk = double(reform(DV.y[m,*])) endif else begin Vbulk = [0., 0., 0.] endelse 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] Vbulk_avg += Vbulk ;------------------------------------ ; 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 azang = 90.-data.theta ;------------------------ ; 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 ; Find azimuthal bin result = min(az_bin-azang[i],a,/abs) if az_bin[a] gt azang[i] then a -= 1 if a eq amax then a -= 1 ; Find polar bin result = min(pol_bin-data.phi[i],p,/abs) if pol_bin[p] gt data.phi[i] then p -= 1 if p eq pmax then p -= 1 if j ge 0 then begin pa_azpol[a, p] = data.pa[i] ;----------------------- ; 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=nns, nne if iecl gt 0 then begin dprint,dlevel=1, '---' dprint,dlevel=1,strtrim(string(iecl),2)+' out of '+strtrim(string(iecm),2)+' particles were not' dprint,dlevel=1,'used in the result because of the bulk-speed subtraction' dprint,dlevel=1,'and their energies have moved outside the defined energy bins.' dprint,dlevel=1,'---' 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(n_samples) vbulk_perp /= float(n_samples) vbulk_vxb /= float(n_samples) vbulk_exb /= float(n_samples) Vbulk_avg /= float(n_samples) bnrm_avg /= float(n_samples) babs_avg /= float(n_samples) 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, {trange:trange,$ egy:wegy, pa:wpa, data:pad, datanorm:padnorm, $ numSlices: n_samples, 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], $ vbulk_para:vbulk_para, vbulk_perp_abs:vbulk_perp, vbulk_vxb:vbulk_vxb, $ vbulk_exb:vbulk_exb, bnrm:bnrm_avg, Vbulk:Vbulk_avg, bfield:bnrm_avg*babs_avg,$ species:species_lc, pa_azpol:pa_azpol, wpol: wpol, waz: waz} END