;+
; PROCEDURE:
; mvn_swia_dirEt
; PURPOSE:
; Makes directional E-t spectrograms in the specified frame from SWIA Coarse data.
; a tplot variable will be generated for only the specified phi
; and theta bounds in the chosen coordinate frame.
; CALLING SEQUENCE:
; mvn_swia_diret_any
; INPUTS:
; None (SWIA data and SPICE kernels need to have been loaded)
; KEYWORDS:
; all optional
; FRAME: specifies the frame (Def: 'MSO')
; UNITS: specifies the units ('eflux', 'counts', etc.) (Def: 'eflux')
; ARCHIVE: uses archive data instead of survey
; ATTVEC: generates tplot variables showing SWIA XYZ vectors in the specified frame
; TRANGE: time range to compute directional spectra (Def: all)
; THETA_bounds: polar angle range in degrees [MAX: 0 TO 180] over which data is wanted
; (Def: 45 to 135)
; PHI_bounds: azimuth angle range in degrees [MAX: 0 to 360] over which data is wanted
; (Def: 315 to 45)
; CREATED BY:
; Rob Lillis on 2016-06-22, modified from Yuki Harada's mvn_swia_dirEt
;
; $LastChangedBy: rlillis3 $
; $LastChangedDate: 2016-06-23 18:15:58 -0700 (Thu, 23 Jun 2016) $
; $LastChangedRevision: 21359 $
; $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/spdsoft/tags/spedas_3_1/projects/maven/swia/mvn_swia_diret_any.pro $
;-
pro mvn_swia_diret_any, frame=frame, units=units, archive=archive, attvec=attvec, trange=trange, verbose=verbose, theta_bounds = theta_bounds,phi_bounds = phi_bounds
common mvn_swia_data
if not keyword_set(frame) then frame='MSO'
if not keyword_set(units) then units='eflux'
if not keyword_set (theta_bounds) then theta_bounds = [45.0, 135.0]
if not keyword_set (phi_bounds) then phi_bounds = [315,45.0]
; very important, change the bounds to the same coordinate system used
; by Yuki, i.e. latitude instead of co-latitude and longitude from
; -180 to 180
theta_bounds_new = reverse (90.0 - theta_bounds)
phi_bounds_new = convert_azimuth (phi_bounds)
if keyword_set(archive) then time = swica.time_unix else time = swics.time_unix
if keyword_set(thld_theta) then thld_theta = abs(thld_theta) else thld_theta = 45
if keyword_set(trange) then begin
idx = where(time ge trange[0] and time le trange[1], idx_cnt)
if idx_cnt gt 0 then time = time[idx] else begin
dprint,dlevel=1,verbose=verbose,'No data in the specified time range.'
return
endelse
endif
center_time = dblarr(n_elements(time))
energy = fltarr(n_elements(time),48)
eflux = fltarr(n_elements(time),48)
pX_new = fltarr(n_elements(time),3)
pY_new = fltarr(n_elements(time),3)
pZ_new = fltarr(n_elements(time),3)
for i=0ll,n_elements(time)-1 do begin ;- time loop
if i mod 1000 eq 0 then dprint,dlevel=1,verbose=verbose,i,' /',n_elements(time)
d = mvn_swia_get_3dc(time[i],archive=archive)
d = conv_units(d,units)
center_time[i] = (d.time+d.end_time)/2.d
energy[i,*] = d.energy[*,0]
sphere_to_cart,1.,d.theta,d.phi,vx,vy,vz
q = spice_body_att('MAVEN_SWIA',frame,center_time[i], $
/quaternion,check_object='MAVEN_SPACECRAFT', $
verbose=-1)
t2 = q[0]*q[1] ;- cf. quaternion_rotation.pro
t3 = q[0]*q[2]
t4 = q[0]*q[3]
t5 = -q[1]*q[1]
t6 = q[1]*q[2]
t7 = q[1]*q[3]
t8 = -q[2]*q[2]
t9 = q[2]*q[3]
t10 = -q[3]*q[3]
vxnew = 2*( (t8 + t10)*vx + (t6 - t4)*vy + (t3 + t7)*vz ) + vx
vynew = 2*( (t4 + t6)*vx + (t5 + t10)*vy + (t9 - t2)*vz ) + vy
vznew = 2*( (t7 - t3)*vx + (t2 + t9)*vy + (t5 + t8)*vz ) + vz
thetanew = 90. - acos(vznew)*!radeg
phinew = atan(vynew,vxnew)*!radeg
if keyword_set(attvec) then begin
pX_new[i,0] = 2*( (t8 + t10)*1. + (t6 - t4)*0. + (t3 + t7)*0. ) + 1.
pX_new[i,1] = 2*( (t4 + t6)*1. + (t5 + t10)*0. + (t9 - t2)*0. ) + 0.
pX_new[i,2] = 2*( (t7 - t3)*1. + (t2 + t9)*0. + (t5 + t8)*0. ) + 0.
pY_new[i,0] = 2*( (t8 + t10)*0. + (t6 - t4)*1. + (t3 + t7)*0. ) + 0.
pY_new[i,1] = 2*( (t4 + t6)*0. + (t5 + t10)*1. + (t9 - t2)*0. ) + 1.
pY_new[i,2] = 2*( (t7 - t3)*0. + (t2 + t9)*1. + (t5 + t8)*0. ) + 0.
pZ_new[i,0] = 2*( (t8 + t10)*0. + (t6 - t4)*0. + (t3 + t7)*1. ) + 0.
pZ_new[i,1] = 2*( (t4 + t6)*0. + (t5 + t10)*0. + (t9 - t2)*1. ) + 0.
pZ_new[i,2] = 2*( (t7 - t3)*0. + (t2 + t9)*0. + (t5 + t8)*1. ) + 1.
endif
; now find the indices where the angles are within the ranges we want
if phi_bounds_new [0] gt phi_bounds_new [1] then idx = $
where( (thetanew) ge theta_bounds_new[0] and $
(thetanew) le theta_bounds_new[1] and $
((phinew) ge phi_bounds_new [0] or $
(phinew) le phi_bounds_new [1]) , idx_cnt ) else idx = $
where( (thetanew) ge theta_bounds_new[0] and $
(thetanew) le theta_bounds_new[1] and $
(phinew) ge phi_bounds_new [0] and $
(phinew) le phi_bounds_new [1] , idx_cnt )
if idx_cnt gt 0 then begin
w = d.data * 0.
w[idx] = 1.
if strlowcase(units) ne 'counts' then $
eflux[i,*] = total(d.data*d.domega*w,2)/total(d.domega*w,2) $
else $
eflux[i,j] = total(d.data*w,2)
endif else eflux[i,*] = !values.f_nan
endfor ;- time loop end
if keyword_set(archive) then type = 'swica' else type = 'swics'
store_data,'mvn_'+type+'_en_'+units+'_'+frame+'_theta'+roundst(theta_bounds[0]) +'_' + $
roundst(theta_bounds[1])+'_phi'+roundst(phi_bounds[0]) +'_' + $
roundst(phi_bounds[1]),$
data={x:center_time,y:eflux,v:energy}, $
dlim={spec:1,zlog:1,ylog:1,yrange:minmax(energy),ystyle:1, $
ytitle:type+'!c'+frame+'theta'+roundst(theta_bounds[0]) +'_' + $
roundst(theta_bounds[1])+'!cphi'+roundst(phi_bounds[0]) +'_' + $
roundst(phi_bounds[1])+'!cEnergy [eV]', $
ztitle:units,datagap:180},verbose=verbose
if keyword_set(attvec) then begin
store_data,'mvn_'+type+'_'+frame+'_Xvec', $
data={x:center_time,y:pX_new}, $
dlim={yrange:[-2,2],ystyle:1,labflag:1, $
ytitle:'Xswia!c'+frame, $
labels:['x','y','z'],colors:'bgr'},verbose=verbose
store_data,'mvn_'+type+'_'+frame+'_Yvec', $
data={x:center_time,y:pY_new}, $
dlim={yrange:[-2,2],ystyle:1,labflag:1, $
ytitle:'Yswia!c'+frame, $
labels:['x','y','z'],colors:'bgr'},verbose=verbose
store_data,'mvn_'+type+'_'+frame+'_Zvec', $
data={x:center_time,y:pZ_new}, $
dlim={yrange:[-2,2],ystyle:1,labflag:1, $
ytitle:'Zswia!c'+frame, $
labels:['x','y','z'],colors:'bgr'},verbose=verbose
endif
end