;+
;PROCEDURE: swe_engy_snap
;PURPOSE:
; Plots omnidirectional energy spectrum snapshots in a separate window for times
; selected with the cursor in a tplot window. Hold down the left mouse button
; and slide for a movie effect.
;
;USAGE:
; swe_engy_snap
;
;INPUTS:
;
;KEYWORDS:
; UNITS: Plot the data in these units. See mvn_swe_convert_units.
; Default = 'eflux'.
;
; TIMES: Make a plot for these times. (Placeholder only.)
;
; TPLOT: Get energy spectra from tplot variable instead of SWEA
; common block.
;
; FIXY: Use a fixed y-axis range. Default = 1 (yes).
;
; KEEPWINS: If set, then don't close the snapshot window(s) on exit.
;
; ARCHIVE: If set, show shapshots of archive data (A5).
;
; BURST: Synonym for ARCHIVE.
;
; SPEC: Named variable to hold the energy spectrum at the last time
; selected.
;
; SUM: If set, use cursor to specify time ranges for averaging.
;
; TSMO: Smoothing interval, in seconds. Default is no smoothing.
;
; ERROR_BARS: If set, plot energy spectra with error bars. Does not work
; when the TPLOT option (see above) is set, because statistical
; uncertainties are not stored in the tplot variable.
;
; POT: Overplot an estimate of the spacecraft potential. Must run
; mvn_swe_sc_pot first.
;
; SCP: Temporarily override any other estimates of the spacecraft
; potential and force it to be this value.
;
; SHIFTPOT: Correct for the spacecraft potential. Spectra are shifted from
; the spacecraft frame to the plasma frame.
;
; DEMAX: Maximum width of spacecraft potential signature.
;
; PEPEAKS: Overplot the nominal energies of the photoelectron energy peaks
; at 23 and 27 eV.
;
; PEREF: Overplot photoelectron reference spectra
;
; BCK: Plot background level (Potassium-40 decay and penetrating
; particles only).
;
; MAGDIR: Print magnetic field geometry (azim, elev, clock) on the plot.
;
; PDIAG: Plot potential estimator in a separate window.
;
; PXLIM: X limits (Volts) for diagnostic plot.
;
; MB: Perform a Maxwell-Boltzmann fit to determine density and
; temperature. Uses a moment calculation to determine the
; halo density, which is defined as the high energy residual
; after subtracting the best-fit Maxwell-Boltzmann.
;
; KAP: Instead of the halo moment calculation, fit the halo with
; a kappa function to estimate halo density.
;
; MOM: Instead of fitting the core with a Maxwell-Boltzmann, use
; a moment calculation for all energies above the spacecraft
; potential.
;
; ERANGE: Energy range for computing the moment. Only effective when
; keyword MOM is set.
;
; FLEV: Calculate the signal level at this energy, using interpolation
; as needed.
;
; SCAT: Plot the scattered photoelectron population, which is defined
; as the low-energy residual after subtracting the best-fit
; Maxwell-Boltzmann.
;
; SEC: Calculate secondary electron spectrum using McFadden's
; semi-empirical approach.
;
; SSCALE: Scale factor for secondary electron spectrum. Default = 5.
; This scales the secondary electron production efficiency.
;
; DDD: Create an energy spectrum from the nearest 3D spectrum and
; plot for comparison.
;
; ABINS: Anode bin mask (16 elements: 0=off, 1=on). Default = all on.
;
; DBINS: Deflector bin mask (6 elements: 0=off, 1=on). Default = all on.
;
; OBINS: 3D solid angle bin mask (96 elements: 0=off, 1=on).
; Default = reform(ABINS # DBINS).
;
; MASK_SC: Mask solid angle bins that view the spacecraft. Default = yes.
; This masking is in addition to OBINS.
;
; NOERASE: Overplot all spectra after the first.
;
; VOFFSET: Vertical offset when overplotting spectra.
;
; RAINBOW: With NOERASE, overplot spectra using up to 6 different colors.
;
; POPEN: Set this to the name of a postscript file for output.
;
; WSCALE: Window size scale factor.
;
; CSCALE: Character size scale factor.
;
; XRANGE: Override the default horizontal axis range with this.
;
; YRANGE: Override the default vertical axis range with this.
;
; TRANGE: Plot snapshot for this time range. Can be in any
; format accepted by time_double. (This disables the
; interactive time range selection.)
;
; $LastChangedBy: dmitchell $
; $LastChangedDate: 2017-05-08 17:30:59 -0700 (Mon, 08 May 2017) $
; $LastChangedRevision: 23286 $
; $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/spdsoft/tags/spedas_2_00/projects/maven/swea/swe_engy_snap.pro $
;
;CREATED BY: David L. Mitchell 07-24-12
;-
pro swe_engy_snap, units=units, keepwins=keepwins, archive=archive, spec=spec, ddd=ddd, $
abins=abins, dbins=dbins, obins=obins, sum=sum, pot=pot, pdiag=pdiag, $
pxlim=pxlim, mb=mb, kap=kap, mom=mom, scat=scat, erange=erange, $
noerase=noerase, scp=scp, fixy=fixy, pepeaks=pepeaks, $
burst=burst, rainbow=rainbow, mask_sc=mask_sc, sec=sec, $
bkg=bkg, tplot=tplot, magdir=magdir, bck=bck, shiftpot=shiftpot, $
xrange=xrange,yrange=frange,sscale=sscale, popen=popen, times=times, $
flev=flev, pylim=pylim, k_e=k_e, peref=peref, error_bars=error_bars, $
trange=tspan, tsmo=tsmo, wscale=wscale, cscale=cscale, voffset=voffset
@mvn_swe_com
@swe_snap_common
common swe_pot_com, Espan, thresh, dEmax, minflux
mass = 5.6856297d-06 ; electron rest mass [eV/(km/s)^2]
c1 = (mass/(2D*!dpi))^1.5
c2 = (2d5/(mass*mass))
c3 = 4D*!dpi*1d-5*sqrt(mass/2D) ; assume isotropic electron distribution
tiny = 1.e-31
if (size(snap_index,/type) eq 0) then swe_snap_layout, 0
if not keyword_set(archive) then aflg = 0 else aflg = 1
if keyword_set(burst) then aflg = 1
if not keyword_set(units) then units = 'eflux'
if keyword_set(sum) then npts = 2 else npts = 1
if keyword_set(tsmo) then begin
npts = 1
dosmo = 1
delta_t = double(tsmo)/2D
endif else dosmo = 0
if (size(error_bars,/type) eq 0) then ebar = 1 else ebar = keyword_set(error_bars)
if keyword_set(ddd) then dflg = 1 else dflg = 0
if keyword_set(noerase) then oflg = 0 else oflg = 1
if (size(scp,/type) eq 0) then scp = !values.f_nan else scp = float(scp[0])
if (size(fixy,/type) eq 0) then fixy = 1
if keyword_set(fixy) then fflg = 1 else fflg = 0
if keyword_set(rainbow) then rflg = 1 else rflg = 0
if keyword_set(sec) then dosec = 1 else dosec = 0
if not keyword_set(sscale) then sscale = 5D else sscale = double(sscale[0])
if keyword_set(bkg) then dobkg = 1 else dobkg = 0
if keyword_set(shiftpot) then spflg = 1 else spflg = 0
if (n_elements(xrange) ne 2) then xrange = [1.,1.e4]
if not keyword_set(wscale) then wscale = 1.
if not keyword_set(cscale) then cscale = 1.
if not keyword_set(voffset) then voffset = 1.
if (size(popen,/type) eq 7) then begin
psflg = 1
psname = popen[0]
csize1 = 1.2*cscale
csize2 = 1.0*cscale
endif else begin
psflg = 0
csize1 = 1.2*cscale
csize2 = 1.4*cscale
endelse
tflg = 0
if keyword_set(tplot) then begin
get_data,'swe_a4',data=dat,limits=lim,index=i
if (i gt 0) then begin
str_element,lim,'ztitle',units_name,success=ok
if (ok) then units_name = strlowcase(units_name) else units_name = 'unknown'
tspec = {time:dat.x, data:dat.y, energy:dat.v, units_name:units_name}
tflg = 1
endif else print,'No SPEC data found in tplot.'
endif
if (tflg and ebar) then begin
print,"Can't plot error bars with TPLOT option."
ebar = 0
endif
get_data,'alt',data=alt
if (size(alt,/type) eq 8) then begin
doalt = 1
get_data,'sza',data=sza
get_data,'lon',data=lon
get_data,'lat',data=lat
endif else doalt = 0
domag = 0
if keyword_set(magdir) then begin
get_data,'mvn_B_1sec_iau_mars',data=mag
if (size(mag,/type) eq 8) then domag = 1 else domag = 0
endif
if (n_elements(abins) ne 16) then abins = replicate(1B, 16)
if (n_elements(dbins) ne 6) then dbins = replicate(1B, 6)
if (n_elements(obins) ne 96) then begin
obins = replicate(1B, 96, 2)
obins[*,0] = reform(abins # dbins, 96)
obins[*,1] = obins[*,0]
endif else obins = reform(byte(obins)) # [1B,1B]
if (size(mask_sc,/type) eq 0) then mask_sc = 1
if keyword_set(mask_sc) then obins = swe_sc_mask * obins
if (size(pot,/type) eq 0) then dopot = 1 else dopot = keyword_set(pot)
if keyword_set(pepeaks) then dopep = 1 else dopep = 0
if (size(peref,/type) eq 8) then doper = 1 else doper = 0
if keyword_set(scat) then scat = 1 else scat = 0
if keyword_set(mb) then begin
mb = 1
dopot = 1
endif else mb = 0
if keyword_set(kap) then kap = 1 else kap = 0
if keyword_set(mom) then begin
mom = 1
dopot = 1
mb = 0
kap = 0
endif else mom = 0
if keyword_set(pdiag) then begin
get_data,'df',data=df
get_data,'d2f',data=d2f,index=i
if (i gt 0) then pflg = 1 else pflg = 0
endif else pflg = 0
if (not tflg) then begin
if (size(mvn_swe_engy,/type) ne 8) then mvn_swe_makespec
if (aflg) then begin
if (size(mvn_swe_engy_arc,/type) ne 8) then begin
print,"No SPEC archive data."
return
endif
endif else begin
if (size(mvn_swe_engy,/type) ne 8) then begin
print,"No SPEC survey data."
return
endif
endelse
endif
if (fflg) then begin
case strupcase(units) of
'COUNTS' : drange = [1e0, 1e5]
'RATE' : drange = [1e1, 1e6]
'CRATE' : drange = [1e1, 1e6]
'FLUX' : drange = [1e1, 3e8]
'EFLUX' : drange = [1e4, 3e9]
'E2FLUX' : drange = [1e6, 1e11]
'DF' : drange = [1e-18, 1e-8]
else : drange = [0,0]
endcase
endif
if (size(swe_hsk,/type) eq 8) then begin
if (n_elements(swe_hsk) gt 2L) then hflg = 1 else hflg = 0
endif else hflg = 0
if keyword_set(keepwins) then kflg = 0 else kflg = 1
if keyword_set(archive) then aflg = 1 else aflg = 0
case n_elements(tspan) of
0 : tsflg = 0
1 : begin
tspan = time_double(tspan)
tsflg = 1
kflg = 0
end
else : begin
tspan = minmax(time_double(tspan))
tsflg = 1
kflg = 0
end
endcase
; Put up snapshot window(s)
Twin = !d.window
window, /free, xsize=wscale*Eopt.xsize, ysize=wscale*Eopt.ysize, xpos=Eopt.xpos, ypos=Eopt.ypos
Ewin = !d.window
if (hflg) then begin
window, /free, xsize=Hopt.xsize, ysize=Hopt.ysize, xpos=Hopt.xpos, ypos=Hopt.ypos
Hwin = !d.window
endif
if (pflg) then begin
window, /free, xsize=Sopt.xsize, ysize=Sopt.ysize, xpos=Sopt.xpos, ypos=Sopt.ypos
Pwin = !d.window
endif
; Get the spectrum closest to the selected time
ok = 1
print,'Use button 1 to select time; button 3 to quit.'
wset,Twin
if (~tsflg) then begin
trange = 0
ctime,trange,npoints=npts,/silent
if (npts gt 1) then cursor,cx,cy,/norm,/up ; Make sure mouse button released
endif else trange = tspan
if (size(trange,/type) ne 5) then begin
wdelete,Ewin
if (hflg) then wdelete,Hwin
if (pflg) then wdelete,Pwin
wset,Twin
return
endif
if (dosmo) then begin
tmin = min(trange, max=tmax)
trange = [(tmin - delta_t), (tmax + delta_t)]
endif
if (tflg) then begin
if (finite(scp)) then pot = scp else pot = 0.
if ((npts eq 1) and (~dosmo)) then begin
dt = min(abs(trange[0] - tspec.time), i)
spec = {time:tspec.time[i], data:reform(tspec.data[i,*]), $
energy:tspec.energy, units_name:tspec.units_name, $
sc_pot:pot}
endif else begin
tmin = min(trange, max=tmax)
i = where((tspec.time ge tmin) and (tspec.time le tmax), count)
if (count gt 0L) then begin
spec = {time:mean(tspec.time[i]), data:average(tspec.data[i,*],1,/nan), $
energy:tspec.energy, units_name:tspec.units_name, sc_pot:pot}
endif
endelse
endif else begin
spec = mvn_swe_getspec(trange, /sum, archive=aflg, units=units, yrange=yrange)
if (finite(scp)) then pot = scp $
else if (finite(spec.sc_pot)) then pot = spec.sc_pot else pot = 0.
spec.sc_pot = pot
endelse
if (fflg) then yrange = drange
if keyword_set(frange) then yrange = frange
; Correct for spacecraft potential. For instrumental units (COUNTS, RATE, or
; CRATE) only shift in energy. For flux units (FLUX, EFLUX), shift in energy
; and correct the signal level to ensure conservation of phase space density.
if (spflg) then begin
if (stregex(units,'flux',/boo,/fold)) then begin
mvn_swe_convert_units, spec, 'df'
spec.energy -= pot
mvn_swe_convert_units, spec, units
endif else spec.energy -= pot
endif
case strupcase(units) of
'COUNTS' : ytitle = 'Raw Counts'
'RATE' : ytitle = 'Raw Count Rate'
'CRATE' : ytitle = 'Count Rate'
'EFLUX' : ytitle = 'Energy Flux (eV/cm2-s-ster-eV)'
'E2FLUX' : ytitle = 'Energy Flux (eV/cm2-s-ster)'
'FLUX' : ytitle = 'Flux (1/cm2-s-ster-eV)'
'DF' : ytitle = 'Dist. Function (1/cm3-(km/s)3)'
else : ytitle = 'Unknown Units'
endcase
if (hflg) then dt = min(abs(swe_hsk.time - spec.time), jref) ; closest HSK
nplot = 0
xs = 0.68
dys = 0.03
yoff = 1.
while (ok) do begin
x = spec.energy
y = spec.data * yoff
if (ebar) then dy = sqrt(spec.var) * yoff
phi = spec.sc_pot
ys = 0.90
if (~oflg) then yoff *= voffset
if (psflg) then popen, psname + string(nplot,format='("_",i2.2)') $
else wset, Ewin
; Put up an Energy Spectrum with (optionally) model fit, background, scattering, etc.
psym = 10
if ((nplot eq 0) or oflg) then plot_oo,x,y,yrange=yrange,/ysty,xrange=xrange, $
xtitle='Energy (eV)', ytitle=ytitle,charsize=csize2,psym=psym,title=time_string(spec.time), $
xmargin=[10,3] $
else oplot,x,y,psym=psym
if (ebar) then errplot,x,(y-dy)>tiny,y+dy,width=0
if (rflg) then begin
col = (nplot mod 6) + 1
oplot,x,y,psym=psym,color=col
errplot,x,(y-dy)>tiny,y+dy,width=0,color=col
endif
if (dflg) then begin
if (npts gt 1) then ddd = mvn_swe_get3d(trange,/all,/sum) $
else ddd = mvn_swe_get3d(spec.time)
mvn_swe_convert_units, ddd, units
dt = min(abs(swe_hsk.time - ddd.time), kref)
if (ddd.time gt t_mtx[2]) then boom = 1 else boom = 0
indx = where(obins[*,boom] eq 1B, onorm)
omask = replicate(1.,64) # obins[*,boom]
onorm = float(onorm)
spec3d = total(ddd.data*omask,2,/nan)/onorm
oplot,ddd.energy[*,0],spec3d,psym=psym,color=4
endif
; Secondary electrons produced by primary electron impact inside the instrument.
; Method from McFadden, adapted from Dave Evans.
if (dosec) then begin
odat = spec
mvn_swe_convert_units, odat, 'crate'
energy = odat.energy
nenergy = odat.nenergy
sec_spec = dblarr(nenergy)
if (finite(scp)) then pot = scp $
else if (finite(phi)) then pot = phi else pot = 0.
kndx = where(energy gt pot)
kmax = max(kndx)
alpha = 1.35
Tmax = 2.283
Emax = 325.
k = 2.2
scale = sscale
Vbias = 0. ; primaries not passing through exit grid
Erat = (energy + Vbias)/Emax ; effect of V0 cancels when using swe_swp
arg = Tmax*(Erat^alpha) < 80. ; avoid underflow
delta = (Erat^(1. - alpha))*(1. - exp(-arg))/(1. - exp(-Tmax))
eff = scale*(1. - exp(-k*delta))/(1. - exp(-k))
for k=1,kmax-1 do sec_spec[k:kmax] += eff[k]*odat.data[k]/energy[k:kmax]^2.0
odat.data = sec_spec
sec_dat = odat
mvn_swe_convert_units, sec_dat, units
dif_dat = spec
dif_dat.data = (spec.data - sec_dat.data) > 1.
oplot,sec_dat.energy[kndx],sec_dat.data[kndx],color=5,line=2
oplot,dif_dat.energy,dif_dat.data,color=5,psym=10
str_element, spec, 'sec_spec', sec_dat.data, /add ; secondary spectrum
str_element, spec, 'dif_spec', dif_dat.data, /add ; primary spectrum
endif
; Background counts resulting from the wings of the energy response function.
; Empirical: based on fits to solar wind core
; Experimental.
if (dobkg) then begin
odat = spec
mvn_swe_convert_units, odat, 'crate'
energy = odat.energy
nenergy = odat.nenergy
if (finite(scp)) then pot = scp $
else if (finite(phi)) then pot = phi else pot = 0.
kndx = where(energy le pot, kcnt)
if (kcnt gt 0L) then begin
bkg_spec = dblarr(nenergy)
kmax = nenergy - 1
scale = 1.5d-1
for k=1,(nenergy-2) do begin
denergy = energy - energy[k]
bkg_spec[0:k-1] += scale*odat.data[k]/denergy[0:k-1]^2.0
bkg_spec[k+1:kmax] += scale*odat.data[k]/denergy[k+1:kmax]^2.0
endfor
bkg_spec[kndx] = !values.f_nan
odat.data = bkg_spec
bkg_dat = odat
mvn_swe_convert_units, bkg_dat, units
dif_dat = spec
dif_dat.data = (spec.data - bkg_dat.data) > 1.
oplot,bkg_dat.energy,bkg_dat.data,color=5,line=2
oplot,dif_dat.energy,dif_dat.data,color=5,psym=10
endif
endif
if keyword_set(bck) then begin
bck = spec
eff0 = 1./mvn_swe_crosscal('2014-11-15',/silent)
eff = 1./mvn_swe_crosscal(bck.time,/silent)
brate = (0.6*(average(eff,/nan)/eff0))[0]
bck.data[*] = brate ; background crate per anode at periapsis
bck.units_name = 'crate'
mvn_swe_convert_units, bck, units
oplot, bck.energy, bck.data, line=2, color=4 ; periapsis
brate = ((2./3.)*(0.97/0.63) + (1./3.)) ; GCR's + Potassium-40
oplot, bck.energy, bck.data*brate, line=2, color=4 ; apoapsis
bck.data[*] = (1./swe_min_dtc - 1.)/swe_dead ; saturation CRATE
bck.units_name = 'crate'
mvn_swe_convert_units, bck, units
oplot, bck.energy, bck.data, line=2, color=4
endif
if (dopot) then begin
apot = abs(pot)
if (spflg) then oplot,[apot,apot],yrange,line=2,color=4 $
else oplot,[apot,apot],yrange,line=2,color=6
endif
if (dopep) then begin
oplot,[23.,23.],yrange,line=2,color=1
oplot,[27.,27.],yrange,line=2,color=1
oplot,[60.,60.],yrange,line=2,color=1
endif
if (doper) then begin
oplot, peref.x, peref.y
oplot, peref.x, peref.y2
endif
if (doalt) then begin
dt = min(abs(alt.x - spec.time), aref)
xyouts,xs,ys,string(round(alt.y[aref]), format='("ALT = ",i5)'),charsize=csize1,/norm
ys -= dys
if (~mb and ~mom) then begin
xyouts,xs,ys,string(round(sza.y[aref]), format='("SZA = ",i5)'),charsize=csize1,/norm
ys -= dys
endif
endif
if (domag) then begin
dt = min(abs(mag.x - spec.time), mref)
str_element, mag, 'azim', success=ok
if (ok) then begin
xyouts,xs,ys,string(round(mag.azim[mref]), format='("Baz = ",i5)'),charsize=csize1,/norm
ys -= dys
endif
str_element, mag, 'elev', success=ok
if (ok) then begin
xyouts,xs,ys,string(round(mag.elev[mref]), format='("Bel = ",i5)'),charsize=csize1,/norm
ys -= dys
endif
str_element, mag, 'clock', success=ok
if (ok) then begin
xyouts,xs,ys,string(round(mag.clock[mref]), format='("Bclk = ",i5)'),charsize=csize1,/norm
ys -= dys
endif
endif
if (mb) then begin
E1 = spec.energy
F1 = spec.data - spec.bkg
counts = spec
mvn_swe_convert_units, counts, 'counts'
cnts = counts.data
sig2 = counts.var ; variance w/ digitization noise
sdev = F1 * (sqrt(sig2)/(cnts > 1.))
p = swe_maxbol()
if (finite(scp)) then p.pot = scp $
else if (finite(phi)) then p.pot = phi else p.pot = 0.
psep = 2.0
indx = where(E1 gt psep*p.pot)
Fpeak = max(F1[indx],k,/nan)
Epeak = E1[indx[k]]
p.t = Epeak/2.
p.n = Fpeak/(4.*c1*c2*sqrt(p.t)*exp((p.pot/p.t) - 2.))
Elo = Epeak*0.7 < ((Epeak/2.) > (psep*phi))
imb = where((E1 gt Elo) and (E1 lt Epeak*2.))
if (n_elements(erange) gt 1) then begin
Emin = min(erange, max=Emax)
imb = where((E1 ge Emin) and (E1 le Emax))
endif
fit,E1[imb],F1[imb],dy=sdev[imb],func='swe_maxbol',par=p,names='N T',/silent
N_core = p.n
if (kap) then begin
Fh = F1 - swe_maxbol(E1,par=p)
ikap = where(E1 gt Epeak*3.)
Fhmax = max(Fh[ikap],k,/nan)
Ehmax = E1[ikap[k]]
Th = Ehmax/2.
p.k_n = Fhmax/(4.*c1*c2*sqrt(Th)*exp((p.pot/Th) - 2.))
p.k_vh = sqrt(Ehmax/mass)
ikap = where((E1 gt Epeak*0.8) and (E1 lt Epeak*25.))
fit,E1[ikap],F1[ikap],func='swe_maxbol',par=p,names='N T K_N',/silent
fit,E1[ikap],F1[ikap],func='swe_maxbol',par=p,names='N T K_N K_VH',/silent
; fit,E1[indx],F1[indx],func='swe_maxbol',par=p,names='N T K_N K_VH K_K',/silent
N_tot = p.n + p.k_n
pk = p
pk.n = 0.
oplot,E1[ikap],swe_maxbol(E1[ikap],par=pk),color=3
endif else begin
dE = E1
dE[0] = abs(E1[1] - E1[0])
for i=1,62 do dE[i] = abs(E1[i+1] - E1[i-1])/2.
dE[63] = abs(E1[63] - E1[62])
j = where(E1 gt Epeak*2., n_e)
E_halo = E1[j]
F_halo = F1[j] - swe_maxbol(E_halo, par=p)
oplot,E_halo,F_halo,color=1,psym=10
prat = (p.pot/E_halo) < 1.
N_halo = c3*total(dE[j]*sqrt(1. - prat)*(E_halo^(-1.5))*F_halo)
N_tot = N_core + N_halo
endelse
if (spflg) then jndx = indgen(64) else jndx = where(E1 gt p.pot)
col = 4
oplot,E1[jndx],swe_maxbol(E1[jndx],par=p),thick=2,color=col,line=1
oplot,E1[imb],swe_maxbol(E1[imb],par=p),color=col,thick=2
xyouts,xs,ys,string(N_tot,format='("N = ",f5.2)'),color=col,charsize=csize1,/norm
ys -= dys
xyouts,xs,ys,string(p.T,format='("T = ",f5.2)'),color=col,charsize=csize1,/norm
ys -= dys
xyouts,xs,ys,string(p.pot,format='("V = ",f5.2)'),color=col,charsize=csize1,/norm
ys -= dys
if (kap) then begin
xyouts,xs,ys,string(p.k_n,format='("Nh = ",f5.2)'),color=3,charsize=csize1,/norm
ys -= dys
xyouts,xs,ys,string(p.k_vh,format='("Vh = ",f6.0)'),color=3,charsize=csize1,/norm
ys -= dys
xyouts,xs,ys,string(p.k_k,format='("k = ",f5.2)'),color=3,charsize=csize1,/norm
ys -= dys
endif else begin
xyouts,xs,ys,string(N_halo,format='("Nh = ",f5.2)'),color=1,charsize=csize1,/norm
ys -= dys
endelse
if (scat) then begin
kndx = where((E1 gt phi) and (E1 lt Epeak), count)
if (count gt 0L) then begin
x_scat = E1[kndx]
y_scat = F1[kndx] - swe_maxbol(E1[kndx], par=p)
kndx = where(E1 le phi, count)
if (count gt 0L) then begin
x_scat = [x_scat, E1[kndx]]
y_scat = [y_scat, F1[kndx]]
endif
oplot,x_scat,y_scat,color=3,psym=10
endif
endif
endif
if (mom) then begin
eflux = spec
mvn_swe_convert_units, eflux, 'eflux'
E1 = eflux.energy
F1 = eflux.data - eflux.bkg
S1 = sqrt(eflux.var)
dE = E1
dE[0] = abs(E1[1] - E1[0])
for i=1,62 do dE[i] = abs(E1[i+1] - E1[i-1])/2.
dE[63] = abs(E1[63] - E1[62])
if (n_elements(erange) gt 1) then begin
Emin = min(erange, max=Emax)
j = where((E1 ge Emin) and (E1 le Emax), n_e)
endif else begin
if (finite(scp)) then pot = scp $
else if (finite(phi)) then pot = phi else pot = 0.
j = where(E1 gt pot, n_e)
j = j[0:(n_e-2)] ; one channel cushion from s/c potential
n_e--
endelse
oplot,spec.energy[j],spec.data[j],color=1,psym=10
prat = (pot/E1[j]) < 1.
Enorm = c3*dE[j]*sqrt(1. - prat)*(E1[j]^(-1.5))
N_j = Enorm*F1[j]
S_j = Enorm*S1[j]
N_tot = total(N_j)
N_sig = sqrt(total(S_j^2.))
print,"Density = ",N_tot," +/- ",N_sig,format='(a,f6.3,a,f6.3)'
Enorm = (2./3.)*c3*dE[j]*((1. - prat)^1.5)*(E1[j]^(-0.5))
P_j = Enorm*F1[j]
S_j = Enorm*S1[j]
pres = total(P_j)
psig = sqrt(total(S_j^2.))
temp = pres/N_tot ; temperature corresponding to kinetic energy density
tsig = temp*sqrt((N_sig/N_tot)^2. + (psig/pres)^2.)
print,"Temperature = ",temp," +/- ",tsig,format='(a,f6.3,a,f6.3)'
xyouts,xs,ys,string(N_tot,format='("N = ",f6.3)'),color=1,charsize=csize1,/norm
ys -= dys
xyouts,xs,ys,string(temp,format='("T = ",f6.2)'),color=1,charsize=csize1,/norm
ys -= dys
xyouts,xs,ys,string(pot,format='("V = ",f6.2)'),color=col,charsize=csize1,/norm
ys -= dys
endif
if (~mb and ~mom) then begin
xyouts,xs,ys,string(pot,format='("V = ",f6.2)'),color=col,charsize=csize1,/norm
ys -= dys
endif
if keyword_set(flev) then begin
logE0 = alog10((float(flev) > min(spec.energy)) < max(spec.energy))
logE = alog10(spec.energy)
logF = alog10(spec.data)
E0 = 10.^logE0
F0 = 10.^interpol(logF, logE, logE0)
oplot,[E0],[F0],psym=4,color=6,symsize=1.5,thick=2
xyouts,xs,ys,string(F0,format='("F = ",e8.2)'),color=6,charsize=csize1,/norm
ys -= dys
Emsg = strtrim(string(E0,format='(f7.1)'),2)
print,"Flux (",Emsg," eV) = ",F0,format='(a,a,a,e8.2)'
endif
if (dflg) then begin
xyouts,xs,ys,'3D',charsize=csize1,/norm,color=4
ys -= dys
endif
if keyword_set(k_e) then begin
xyouts,xs,ys,string(swe_Ke[0],format='("Ke = ",f4.2)'),charsize=csize1,/norm
ys -= dys
endif
if (psflg) then pclose
if (pflg) then begin
wset, Pwin
xs = 0.71
ys = 0.90
dys = 0.03
if not keyword_set(pxlim) then xlim = Espan else xlim = minmax(pxlim)
if not keyword_set(pylim) then begin
indx = where((df.v ge xlim[0]) and (df.v le xlim[1]))
ymin = min(df.y[indx],/nan) < min(d2f.y[indx],/nan)
ymax = max(df.y[indx],/nan) > max(d2f.y[indx],/nan)
ylim = [floor(100.*ymin), ceil(100.*ymax)]/100.
endif else ylim = minmax(pylim)
dt = min(abs(d2f.x - trange[0]), kref)
px = reform(d2f.v[kref,*])
py = reform(df.y[kref,*])
py2 = reform(d2f.y[kref,*])
n_e = n_elements(py)
zcross = py2 * shift(py2,1)
zcross[0] = 1.
indx = where((zcross lt 0.) and (py gt thresh[0]), ncross)
title = string(spec.sc_pot,format='("Potential = ",f5.1," V")')
plot,px,py,xtitle='Potential (V)',ytitle='dF and d2F',$
xrange=xlim,/xsty,yrange=ylim,/ysty,title=title,charsize=csize2
oplot,[spec.sc_pot,spec.sc_pot],ylim,line=2,color=6
oplot,px,py2,color=4
oplot,xlim,[0,0],line=2
oplot,xlim,[thresh,thresh],line=2,color=5
if (ncross gt 0L) then begin
k = max(indx) ; lowest energy feature above threshold
pymax = py[k]
; pymax = max(py,k0) ; largest slope feature above threshold
; k = k0
pymin = pymax/2.
while ((py[k] gt pymin) and (k lt n_e-1)) do k++
kmax = k
k = max(indx)
while ((py[k] gt pymin) and (k gt 0)) do k--
kmin = k
dE = px[kmin] - px[kmax]
oplot,[px[kmin],px[kmax]],[pymin,pymin],color=6
if ((kmax eq (n_e-1)) or (kmin eq 0)) then dE = 2.*dEmax
if (dE lt dEmax) then k = max(indx) else k = -1 ; only accept narrow features
; if (dE lt dEmax) then k = k0 else k = -1 ; only accept narrow features
for j=0,(ncross-1) do oplot,[px[indx[j]],px[indx[j]]],ylim,color=2
if (k gt 0) then begin
xyouts,xs,ys,string(px[k],format='("V = ",f6.2)'),color=col,charsize=csize1,/norm
ys -= dys
oplot,[px[k],px[k]],ylim,color=6,line=2
endif else begin
xyouts,xs,ys,"V = NaN",color=col,charsize=csize1,/norm
ys -= dys
endelse
xyouts,xs,ys,string(dE,format='("dE = ",f6.2)'),charsize=csize1,/norm
ys -= dys
xyouts,xs,ys,string(dEmax,format='("dEmax = ",f6.2)'),charsize=csize1,/norm
ys -= dys
xyouts,xs,ys,string(thresh,format='("thresh = ",f6.2)'),charsize=csize1,/norm
ys -= dys
endif
endif
; Print out housekeeping in another window
if (hflg) then begin
wset, Hwin
csize = 1.4
x1 = 0.05
x2 = 0.75
x3 = x2 - 0.12
y1 = 0.95 - 0.035*findgen(28)
fmt1 = '(f7.2," V")'
fmt2 = '(f7.2," C")'
fmt3 = '(i2)'
j = jref
k = swe_hsk[j].ssctl
if (k lt 4) then tabnum = mvn_swe_tabnum(swe_hsk[j].chksum[k]) else tabnum = -1
erase
xyouts,x1,y1[0],/normal,"SWEA Housekeeping",charsize=csize
xyouts,x1,y1[1],/normal,time_string(swe_hsk[j].time),charsize=csize
xyouts,x1,y1[3],/normal,"P28V",charsize=csize
xyouts,x2,y1[3],/normal,string(swe_hsk[j].P28V,format=fmt1),charsize=csize,align=1.0
xyouts,x1,y1[4],/normal,"MCP28V",charsize=csize
xyouts,x2,y1[4],/normal,string(swe_hsk[j].MCP28V,format=fmt1),charsize=csize,align=1.0
xyouts,x1,y1[5],/normal,"NR28V",charsize=csize
xyouts,x2,y1[5],/normal,string(swe_hsk[j].NR28V,format=fmt1),charsize=csize,align=1.0
xyouts,x1,y1[6],/normal,"MCPHV",charsize=csize
xyouts,x2,y1[6],/normal,string(sigfig(swe_hsk[j].MCPHV,3),format=fmt1),charsize=csize,align=1.0
xyouts,x1,y1[7],/normal,"NRV",charsize=csize
xyouts,x2,y1[7],/normal,string(swe_hsk[j].NRV,format=fmt1),charsize=csize,align=1.0
xyouts,x1,y1[9],/normal,"P12V",charsize=csize
xyouts,x2,y1[9],/normal,string(swe_hsk[j].P12V,format=fmt1),charsize=csize,align=1.0
xyouts,x1,y1[10],/normal,"N12V",charsize=csize
xyouts,x2,y1[10],/normal,string(swe_hsk[j].N12V,format=fmt1),charsize=csize,align=1.0
xyouts,x1,y1[11],/normal,"P5AV",charsize=csize
xyouts,x2,y1[11],/normal,string(swe_hsk[j].P5AV,format=fmt1),charsize=csize,align=1.0
xyouts,x1,y1[12],/normal,"N5AV",charsize=csize
xyouts,x2,y1[12],/normal,string(swe_hsk[j].N5AV,format=fmt1),charsize=csize,align=1.0
xyouts,x1,y1[13],/normal,"P5DV",charsize=csize
xyouts,x2,y1[13],/normal,string(swe_hsk[j].P5DV,format=fmt1),charsize=csize,align=1.0
xyouts,x1,y1[14],/normal,"P3P3DV",charsize=csize
xyouts,x2,y1[14],/normal,string(swe_hsk[j].P3P3DV,format=fmt1),charsize=csize,align=1.0
xyouts,x1,y1[15],/normal,"P2P5DV",charsize=csize
xyouts,x2,y1[15],/normal,string(swe_hsk[j].P2P5DV,format=fmt1),charsize=csize,align=1.0
xyouts,x1,y1[17],/normal,"ANALV",charsize=csize
xyouts,x2,y1[17],/normal,string(swe_hsk[j].ANALV,format=fmt1),charsize=csize,align=1.0
xyouts,x1,y1[18],/normal,"DEF1V",charsize=csize
xyouts,x2,y1[18],/normal,string(swe_hsk[j].DEF1V,format=fmt1),charsize=csize,align=1.0
xyouts,x1,y1[19],/normal,"DEF2V",charsize=csize
xyouts,x2,y1[19],/normal,string(swe_hsk[j].DEF2V,format=fmt1),charsize=csize,align=1.0
xyouts,x1,y1[20],/normal,"V0V",charsize=csize
xyouts,x2,y1[20],/normal,string(swe_hsk[j].V0V,format=fmt1),charsize=csize,align=1.0
xyouts,x1,y1[22],/normal,"ANALT",charsize=csize
xyouts,x2,y1[22],/normal,string(swe_hsk[j].ANALT,format=fmt2),charsize=csize,align=1.0
xyouts,x1,y1[23],/normal,"LVPST",charsize=csize
xyouts,x2,y1[23],/normal,string(swe_hsk[j].LVPST,format=fmt2),charsize=csize,align=1.0
xyouts,x1,y1[24],/normal,"DIGT",charsize=csize
xyouts,x2,y1[24],/normal,string(swe_hsk[j].DIGT,format=fmt2),charsize=csize,align=1.0
xyouts,x1,y1[26],/normal,"SWEEP TABLE",charsize=csize
xyouts,x2,y1[26],/normal,string(tabnum,format=fmt3),charsize=csize,align=1.0
endif
; Get the next button press
nplot++
if (~tsflg) then begin
wset,Twin
trange = 0
ctime,trange,npoints=npts,/silent
if (npts gt 1) then cursor,cx,cy,/norm,/up ; make sure mouse button is released
if (size(trange,/type) eq 5) then begin
if (dosmo) then begin
tmin = min(trange, max=tmax)
trange = [(tmin - delta_t), (tmax + delta_t)]
endif
if (tflg) then begin
if (finite(scp)) then pot = scp else pot = 0.
if ((npts eq 1) and (~dosmo)) then begin
dt = min(abs(trange[0] - tspec.time), i)
spec = {time:tspec.time[i], data:reform(tspec.data[i,*]), $
energy:tspec.energy, units_name:tspec.units_name, $
sc_pot:pot}
endif else begin
tmin = min(trange, max=tmax)
i = where((tspec.time ge tmin) and (tspec.time le tmax), count)
if (count gt 0L) then begin
spec = {time:mean(tspec.time[i]), data:average(tspec.data[i,*],1,/nan), $
energy:tspec.energy, units_name:tspec.units_name, sc_pot:pot}
endif
endelse
endif else begin
spec = mvn_swe_getspec(trange, /sum, archive=aflg, units=units, yrange=yrange)
if (finite(scp)) then pot = scp $
else if (finite(spec.sc_pot)) then pot = spec.sc_pot else pot = 0.
spec.sc_pot = pot
endelse
if (spflg) then begin
if (stregex(units,'flux',/boo,/fold)) then begin
mvn_swe_convert_units, spec, 'df'
spec.energy -= pot
mvn_swe_convert_units, spec, units
endif else spec.energy -= pot
endif
if (fflg) then yrange = drange
if keyword_set(frange) then yrange = frange
if (hflg) then dt = min(abs(swe_hsk.time - trange[0]), jref)
endif else ok = 0
endif else ok = 0
endwhile
if (kflg) then begin
wdelete, Ewin
if (hflg) then wdelete, Hwin
if (pflg) then wdelete, Pwin
endif
wset, Twin
return
end