swe_engy_snap.pro (Documentation for /home/davin/idl/socware/)

;+
;PROCEDURE:   swe_engy_snap
;PURPOSE:
;  Plots 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.  This procedure depends on running swe_plot_dpu first, which unpacks the
;  A4 packets, creating 16 energy spectra per packet.
;
;  If housekeeping data exist (almost always the case), then they are displayed 
;  as text in a small separate window.
;
;USAGE:
;  swe_engy_snap
;
;INPUTS:
;
;KEYWORDS:
;       UNITS:         Plot the data in these units.  See mvn_swe_convert_units.
;                      Default = 'eflux'.
;
;       KEEPWINS:      If set, then don't close the snapshot window(s) on exit.
;
;       ARCHIVE:       If set, show shapshots of archive data (A5).
;
;       SPEC:          Named variable to hold the energy spectrum at the last time
;                      selected.
;
;       SUM:           If set, use cursor to specify time ranges for averaging.
;
;       POT:           Overplot an estimate of the spacecraft potential.  Must run
;                      mvn_swe_sc_pot first.
;
;       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.
;
;       SCAT:          Plot the scattered photoelectron population, which is defined
;                      as the low-energy residual after subtracting the best-fit
;                      Maxwell-Boltzmann.
;
;       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.
;
;       NOERASE:       Overplot all spectra after the first.
;
; $LastChangedBy: dmitchell $
; $LastChangedDate: 2014-11-02 14:44:54 -0800 (Sun, 02 Nov 2014) $
; $LastChangedRevision: 16114 $
; $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/spdsoft/trunk/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, sum=sum, pot=pot, pdiag=pdiag, $
                   pxlim=pxlim, mb=mb, kap=kap, mom=mom, scat=scat, erange=erange, $
                   noerase=noerase, thresh=thresh, scp=scp

  @mvn_swe_com
  common snap_layout, Dopt, Sopt, Popt, Nopt, Copt, Eopt, Hopt

  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

  if not keyword_set(archive) then aflg = 0 else aflg = 1
  if not keyword_set(units) then units = 'eflux'
  if keyword_set(sum) then npts = 2 else npts = 1
  if keyword_set(ddd) then dflg = 1 else dflg = 0
  if not keyword_set(abins) then abins = replicate(1, 16)
  if not keyword_set(dbins) then dbins = replicate(1, 6)
  if keyword_set(noerase) then oflg = 0 else oflg = 1
  if not keyword_set(scp) then scp = 0. else scp = float(scp[0])

  case n_elements(thresh) of
    0    : thresh = [0.05, 0.025, 1.8]
    1    : thresh = [float(thresh), 0.025, 1.8]
    2    : thresh = [float(thresh), 1.8]
    else : thresh = float(thresh[0:2])
  endcase

  get_data,'mvn_swe_shape_par',data=par
  if (size(par,/type) ne 8) then begin
    mvn_swe_shape_par
    get_data,'mvn_swe_shape_par',data=par
  endif

  get_data,'alt',data=alt
  if (size(alt,/type) eq 8) then doalt = 1 else doalt = 0
  
  obins = reform(abins # dbins, 96)
  indx = where(obins eq 1, onorm)
  omask = replicate(1.,64) # obins
  onorm = float(onorm)
  
  if keyword_set(pot) then dopot = 1 else dopot = 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 (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

  if (size(swe_hsk,/type) ne 8) then hflg = 0 else hflg = 1
  if keyword_set(keepwins) then kflg = 0 else kflg = 1
  if keyword_set(archive) then aflg = 1 else aflg = 0

; Put up snapshot window(s)

  Twin = !d.window

  if (size(Dopt,/type) ne 8) then swe_snap_layout, 0

  window, /free, xsize=Eopt.xsize, ysize=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 the selected time
  
  ok = 1

  print,'Use button 1 to select time; button 3 to quit.'

  wset,Twin
  ctime2,trange,npoints=npts,/silent,button=button

  if (size(trange,/type) ne 5) then begin
    wdelete,Ewin
    if (hflg) then wdelete,Hwin
    if (pflg) then wdelete,Pwin
    wset,Twin
    return
  endif
  
  spec = mvn_swe_getspec(trange, /sum, archive=aflg, units=units, yrange=yrange)
  
  case strupcase(spec.units_name) of
    'COUNTS' : ytitle = 'Raw Counts'
    'RATE'   : ytitle = 'Raw Count Rate'
    'CRATE'  : ytitle = 'Count Rate'
    'EFLUX'  : ytitle = 'Energy Flux (eV/cm2-s-ster-eV)'
    '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
  
  first = 1
  xs = 0.71
  dys = 0.03

  while (ok) do begin

    x = spec.energy
    y = spec.data
    phi = spec.sc_pot
    ys = 0.90

    wset, Ewin

; Put up an Energy Spectrum

    psym = 10

    if (first or oflg) then plot_oo,x,y,yrange=yrange,/ysty,xtitle='Energy (eV)',ytitle=ytitle, $
            charsize=1.4,psym=psym,title=time_string(spec.time) $
                       else oplot,x,y,psym=psym

    if (dflg) then begin
      if (npts gt 1) then ddd = mvn_swe_get3d([tmin,tmax],/all,/sum) $
                     else ddd = mvn_swe_get3d(spec.time)
      mvn_swe_convert_units, ddd, spec.units_name
      dt = min(abs(swe_hsk.time - ddd.time), kref)
      spec3d = total(ddd.data*omask,2)/onorm
      oplot,ddd.energy[*,0],spec3d,psym=psym,color=4
    endif

    if (dopot) then begin
      if (finite(phi)) then pot = phi else pot = scp
      oplot,[pot,pot],yrange,line=2,color=6
    endif
    
    if (doalt) then begin
      dt = min(abs(alt.x - spec.time), aref)
      xyouts,xs,ys,string(alt.y[aref], format='("ALT = ",f6.1)'),charsize=1.2,/norm
      ys -= dys
    endif
    
    if (mb) then begin
      E1 = spec.energy
      F1 = spec.data
      
      counts = conv_units(spec,'counts')
      cnts = counts.data
      sig2 = counts.var  ; variance w/ digitization noise
      sdev = F1 * (sqrt(sig2)/(cnts > 1.))

      p = swe_maxbol()
      if (finite(phi)) then p.pot = phi else p.pot = scp
      indx = where(E1 gt 2.*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.8 < ((Epeak/2.) > (2.*phi))
      imb = where((E1 gt Elo) and (E1 lt Epeak*3.))

      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

      jndx = where(E1 gt p.pot)
      col = 4
      oplot,E1[jndx],swe_maxbol(E1[jndx],par=p),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=1.2,/norm
      ys -= dys
      xyouts,xs,ys,string(p.T,format='("T = ",f5.2)'),color=col,charsize=1.2,/norm
      ys -= dys
      xyouts,xs,ys,string(p.pot,format='("V = ",f5.2)'),color=6,charsize=1.2,/norm
      ys -= dys
      if (kap) then begin
        xyouts,xs,ys,string(p.k_n,format='("Nh = ",f5.2)'),color=3,charsize=1.2,/norm
        ys -= dys
        xyouts,xs,ys,string(p.k_vh,format='("Vh = ",f6.0)'),color=3,charsize=1.2,/norm
        ys -= dys
        xyouts,xs,ys,string(p.k_k,format='("k = ",f5.2)'),color=3,charsize=1.2,/norm        
        ys -= dys
      endif else begin
        xyouts,xs,ys,string(N_halo,format='("Nh = ",f5.2)'),color=1,charsize=1.2,/norm
        ys -= dys
      endelse

      if (scat) then begin
        kndx = where((E1 gt phi) and (E1 lt Epeak))
        oplot,E1[kndx],(F1[kndx] - swe_maxbol(E1[kndx], par=p)),color=5,psym=10
      endif
    endif
      
    if (mom) then begin
      E1 = spec.energy
      F1 = spec.data

      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(phi) then pot = phi else pot = scp
        j = where(E1 gt pot, n_e)
        j1 = max(j)
        Fmax = max(F1[0:j1],jmax,/nan)
        Fmin = min(F1[jmax:j1],jmin,/nan)
        j = where(E1 ge E1[jmin+jmax], n_e)
      endelse

      oplot,E1[j],F1[j],color=1,psym=10
      prat = (pot/E1[j]) < 1.

      N_tot = c3*total(dE[j]*sqrt(1. - prat)*(E1[j]^(-1.5))*F1[j])

      Fmax = max(F1[j],k,/nan)
      Emax = E1[j[k]]
      temp = Emax/2.

      xyouts,xs,ys,string(N_tot,format='("N = ",f6.2)'),color=1,charsize=1.2,/norm
      ys -= dys
      xyouts,xs,ys,string(temp,format='("T = ",f6.2)'),color=1,charsize=1.2,/norm
      ys -= dys
      xyouts,xs,ys,string(pot,format='("V = ",f6.2)'),color=6,charsize=1.2,/norm
      ys -= dys
    endif
    
    if (pflg) then begin
      wset, Pwin
      
      if not keyword_set(pxlim) then xlim = [0.,30.] else xlim = minmax(pxlim)

      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.

      dt = min(abs(d2f.x - trange[0]), kref)
      px = reform(d2f.v[kref,*])
      py = reform(df.y[kref,*])
      py2 = reform(d2f.y[kref,*])
      
      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=1.4
      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[0],thresh[0]],line=2,color=5      
      oplot,xlim,[thresh[1],thresh[1]],line=2,color=5      

      if (ncross gt 0L) then begin
        pymax = max(py[indx],k)
        py2max = max(py2[0L:indx[k]])
        
        dt = min(abs(spec.time - par.x), sndx)
        par_y = par.y[sndx]
        
        if ((py2max gt thresh[1]) and (par_y gt thresh[2])) then k = indx[k] else k = -1

        for j=0,(ncross-1) do oplot,[px[indx[j]],px[indx[j]]],ylim,color=2

        if (k gt 0) then begin
          xyouts,xs,0.90,string(px[k],format='("V = ",f6.2)'),color=6,charsize=1.2,/norm
          oplot,[px[k],px[k]],ylim,color=6,line=2
        endif

      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
    
      tabnum = mvn_swe_tabnum(swe_hsk[j].chksum[swe_hsk[j].ssctl])
          
      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

    first = 0

    wset,Twin
    ctime2,trange,npoints=npts,/silent,button=button

    if (size(trange,/type) eq 5) then begin
      spec = mvn_swe_getspec(trange, /sum, archive=aflg, units=units, yrange=yrange)
      if (hflg) then dt = min(abs(swe_hsk.time - trange[0]), jref)
    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