;+
;PROCEDURE:	mvn_sta_prod_cal,all=all,units=units
;PURPOSE:	
;	Calibrates STATIC data products generated by mvn_sta_apid_handler.pro
;INPUT:		
;
;KEYWORDS:
;	all	0/1		if not set, then raw tplot structures are deleted
;	units	string		select the units for generated tplot structures
;	apids	strarr		if set, selectes subset of apids to generate tplot structures
;
;
;CREATED BY:	J. McFadden	2012/10/04
;VERSION:	1
;LAST MODIFICATION:  2013/06/11
;MOD HISTORY:
;
;NOTES:	  
;	
;-
pro mvn_sta_prod_cal,all=all,units=units,apids=apids

if keyword_set(units) then begin
	un = units
	if units eq 'cnts' then un=0
	if units eq 'eflx' then un=1
endif else un=0

cols=get_colors()

; decompression array

decomp19=dblarr(256)
decomp19[0:127]=$
[0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0,$
16.0,17.0,18.0,19.0,20.0,21.0,22.0,23.0,24.0,25.0,26.0,27.0,28.0,29.0,30.0,31.0,$
32.5,34.5,36.5,38.5,40.5,42.5,44.5,46.5,48.5,50.5,52.5,54.5,56.5,58.5,60.5,62.5,$
65.5,69.5,73.5,77.5,81.5,85.5,89.5,93.5,97.5,101.5,105.5,109.5,113.5,117.5,121.5,125.5,$
131.5,139.5,147.5,155.5,163.5,171.5,179.5,187.5,195.5,203.5,211.5,219.5,227.5,235.5,243.5,251.5,$
263.5,279.5,295.5,311.5,327.5,343.5,359.5,375.5,391.5,407.5,423.5,439.5,455.5,471.5,487.5,503.5,$
527.5,559.5,591.5,623.5,655.5,687.5,719.5,751.5,783.5,815.5,847.5,879.5,911.5,943.5,975.5,1007.5,$
1055.5,1119.5,1183.5,1247.5,1311.5,1375.5,1439.5,1503.5,1567.5,1631.5,1695.5,1759.5,1823.5,1887.5,1951.5,2015.5]

decomp19[128:255]=$
[2111.5,2239.5,2367.5,2495.5,2623.5,2751.5,2879.5,3007.5,3135.5,3263.5,3391.5,3519.5,3647.5,3775.5,3903.5,4031.5,$
4223.5,4479.5,4735.5,4991.5,5247.5,5503.5,5759.5,6015.5,6271.5,6527.5,6783.5,7039.5,7295.5,7551.5,7807.5,8063.5,$
8447.5,8959.5,9471.5,9983.5,10495.5,11007.5,11519.5,12031.5,12543.5,13055.5,13567.5,14079.5,14591.5,15103.5,15615.5,16127.5,$
16895.5,17919.5,18943.5,19967.5,20991.5,22015.5,23039.5,24063.5,25087.5,26111.5,27135.5,28159.5,29183.5,30207.5,31231.5,32255.5,$
33791.5,35839.5,37887.5,39935.5,41983.5,44031.5,46079.5,48127.5,50175.5,52223.5,54271.5,56319.5,58367.5,60415.5,62463.5,64511.5,$
67583.5,71679.5,75775.5,79871.5,83967.5,88063.5,92159.5,96255.5,100351.5,104447.5,108543.5,112639.5,116735.5,120831.5,124927.5,129023.5,$
135167.5,143359.5,151551.5,159743.5,167935.5,176127.5,184319.5,192511.5,200703.5,208895.5,217087.5,225279.5,233471.5,241663.5,249855.5,258047.5,$
270335.5,286719.5,303103.5,319487.5,335871.5,352255.5,368639.5,385023.5,401407.5,417791.5,434175.5,450559.5,466943.5,483327.5,499711.5,516095.5]

; Calibration values and arrays -- nominal values used as place holders
									; ESA Entrance Aperture Hex Grids: 91.9% transmission
									; ESA Exit Grid pair at spyder plate: 89.2% transmission
									; TOF Start/Stop suppression grids:   83.5% & 82.6% transmission
									; TOF foil frame grids 333 lpi: 65% transmission
	geom = 	0.016							; nominal ESA geometric factor 0.001057* R*R*E (from CODIF simulation), R=3.895 cm (dia = 3.071 inches)
	geom = geom*(0.919)^3*0.835*0.62				; grid attenuation (3 entrance grids - 91.9%, 1 foil suppresion grid - 83.5%, 1 foil frame grid - 62%)
	geom = geom*(0.7556)*(0.9)					; post attenuation (ESA entrance and exit posts) 

; we need to include an energy dependent transmission -- impact high energy ions, low energy ions see leakage fields
	geom = geom*(0.892)^2						; ESA exit grids are 89.2% transmission -- impact high energy ions, low energy ions see leakage fields
	geom = geom*(0.9)						; TOF ribs attenuate high energy ions (~0.9)

	scale = 1.046							; energy correction, produces the proper energy scaling for SLUT tables
	dgf45 = [.55,.75,.9,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,.9,.75,.55]	; deflector attentuation for 45 deg deflection, determined Jasper Halekas simulations
	e_att = 0.1  ; 0.09133							; electrostatic grid nominal attenuation factor
	e_ano = replicate(1.,16)					; electrostatic grid attenuator variation with anode, assumes no variation 
	t0    = -2.15							; TOF offset in ns, note that this may drift over time, will need to track it 

; TDC peaks in MH hot  +56C:  TDC2 222,466,709	TDC1 228,472,715	; 20121105_011833_fm_tvac_lpt_hot3,  243.5, 24MHz --> 5.844bins/ns , offsets=21.5,15.5
; TDC peaks in MH room +29C:  TDC2 222,465,709	TDC1 228,471,715	; 20121102_021952_fm_tvac_prod4,     243.5, 24MHz --> 5.844bins/ns , offsets=21.5,15.5
; TDC peaks in MH cold -26C:  TDC2 218,457,698	TDC1 224,463,703	; 20121105_064622_fm_tvac_lpt_cold3, 240.0, 24MHz --> 5.760bins/ns , offsets=22.0,16.0
 	b_ns  = 5.80 							; TOF bins/ns for TOF timing (not sure if this is 5.844 or 5.855 for flight unit) 
	mec = [1.,1.,1.,.3,.04,.01,.01,.01,.01,.01,.02,.1,1.,1.,1.,1.]	; mech gf attenuator variation with anode
	agf = [1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.]		; anode dependent gf that depends on ESA grid transmissions, foil area, etc.

;************************************************************************************************ 
;************************************************************************************************ 

; We need a correction for qual/total  -- a form of dead time due to start droop

;************************************************************************************************ 
;************************************************************************************************ 

; Declare time dependent efficiency arrays -- assume 10 changes in these arrays over time

	ef_time_boundary=dblarr(10)
	ef_time_boundary(0)   = time_double('2001-01-01')
	ef_time_boundary(1)   = time_double('2013-11-01')
	ef_time_boundary(2:9) = time_double('2020-11-01')

	ef1 = fltarr(16,10)						; start tof efficiency versus anode & time, Rate(TimeA and TimeB)/Rate(TimeRST)
	ef2 = fltarr(16,10)						; stop  tof efficiency versus anode & time, Rate(TimeC and TimeD)/Rate(TimeRST)
	ef3 = fltarr(64,10)						; start mass dependent efficiency of foils, assume this is independent of foil  *may not be needed
	ef4 = fltarr(64,10)						; stop  mass dependent efficiency of foils, assume this is independent of foil, *may not be needed
	ef5 = fltarr(64,10)						; start energy dependent efficiency of foils, assume this is independent of foil
	ef6 = fltarr(64,10)						; stop  energy dependent efficiency of foils, assume this is independent of foil

;	for ief1_def = replicate(1.,16

; Independent response arrays
	ano = fltarr(16) & for i=0,15 do ano = (i-7)*22.5		; anode angles




; SLUT parameters for each STATIC mode
	max_md = 8
	slut = fltarr(max_md+1,7)
;	slut[md,*]	  Estart,	 Estop,	defmax,	defctr,	gridon,	gridscale,	maxgrid		; mode name
	slut[0,*] = 	[   50.0,	0.7040,	45.0,	0.,	10.0,	2.0,		25.]		; startup
	slut[1,*] = 	[   50.0,	0.7040,	45.0,	0.,	10.0,	2.0,		25.]		; ram 
	slut[2,*] = 	[  500.0,	1.2340,	45.0,	0.,	 0.0,	2.0,		25.]		; conic 
	slut[3,*] = 	[30000.0,	2.6080,	45.0,	0.,	 0.0,	2.0,		25.]		; pickup 
	slut[4,*] = 	[30000.0,	2.6080,	 0.0,	0.,	 0.0,	2.0,		25.]		; scan
	slut[5,*] = 	[30000.0,	2.6080,	45.0,	0.,	 0.0,	2.0,		25.]		; eclipse
; TBD modes
	slut[6,*] = 	[   28.7,	0.5000,	22.5,	0.,	11.0,	2.0,		25.]		; ram
	slut[7,*] = 	[ 1912.1,	0.5000,	45.0,	0.,	11.3,	2.0,		25.]		; conic
	slut[8,*] = 	[28680.0,	3.0000,	45.0,	0.,	 0.0,	2.0,		25.]		; pickup

; MLUT parameters for each STATIC mode
	mlut = fltarr(max_md+1,64,256) 
	for j=0,8 do mlut[j,*,*] = replicate(1.,64)#reform(replicate(1.,4)#findgen(64),256)			; default MLUT
; 
;	need to include the mlut for pickup mode

; Organize calibration arrays as [mode,energy,def,anode,mass,time] - time is needed for time dependent efficiencies

; Declare arrays

	nrg = fltarr(max_md+1,64)					; energy sweep table, determined from SLUTs
	def = fltarr(max_md+1,16)					; deflector values, determined from SLUTs

	ggf = fltarr(max_md+1,64,16) 					; electrostatic grid gf attenuation vs energy, determined from SLUTs
	dgf = fltarr(max_md+1,16) 					; def gf attenuation vs deflection, determined from SLUTs

	
	agf = fltarr(max_md+1,64,16,16,64)				; total gf as a function of (energy,def,anode,mass)

	tof=fltarr(max_md+1,64)						; tof in ns for product bins (determined from MLUTs, assumes MLUT accounts for changes in tof due to initial energy)
	mss=fltarr(max_md+1,64)						; mass based on TOF assuming no energy straggling
	dms=fltarr(max_md+1,64)						; accounts for variable TOF mass bin range (determined from MLUTs)
	ms2=fltarr(max_md+1,16,64)					; could be used to account for variations in foil thickness to correct mass

; TOF arrays
; 	15keV H+ takes 11.79 ns for 2 cm
;	TOFs (TDC1,TDC2) are converted from 1024 TOF bins to 256 mass bins with a pseudo log compression - see FPGA spec
;		Formula: 0-255 -> 0-127, 256-511 -> 128-191, 512-1024 -> 192-255  
;		highest resolution mass bins average 2 bins, introduces factor 1/2
; 		Typical product mass compression (256->64) introduces factor 1/4 
	dt = 11.79*b_ns/4./2.						; TOF in bins for 15keV proton					

	tof=fltarr(max_md+1,64)						; tof
	mas=fltarr(max_md+1,64)						; mass based on TOF assuming no energy straggling
	dm=fltarr(max_md+1,64)						; accounts for variable TOF mass bin range 
; 		tof_def:  256 bins summed -> 64 bins
		tof_def = findgen(64) & tof_def(32:47)=32.+2.*findgen(16) & tof_def(48:63)=64.+4*findgen(16) 
	tof_def = findgen(64) & tof_def(32:47)=32.+2.*findgen(16) & tof_def(48:63)=64.+4*findgen(16) 




; Generate the energy, deflector, attenuation arrays

	for md=0,max_md do begin
	   ; energy array
		aa = (slut[md,0]/slut[md,1])^(1./63) 
		nrg[md,*] = scale*slut[md,0]/aa^findgen(64) 
	   ; deflector array
		def[md,*] = slut[md,2]*(findgen(16)-7.5)						; this may need a sign change depending upon coordinate transformations
	   ; def attenuation
		if slut[md,2] eq 45. then dgf[md,*]=dgf45 else dgf[md,*]=1.				; this assumes either defmax is either 45. or <30. deg.  
	   ; grid attenuation
		ggf[md,*]=1.										; no electrostatic attenuation unless grid is on
		ind = where(nrg[md,*] lt slut[md,4],cnt) 						; find energy steps when grid is on
		if cnt gt 0 then ggf[md,ind]=e_att							; this assumes uniform grid attenuation

;		, start out assuming no anode dependence

	endfor
 
; Generate tof/mass tables

	for md=0,max_md do begin
		tof(md,*) = tof_def					; this is a default array until I incorporate MLUTs
		mas(md,*) = ((tof_def-t0)/dt)^2				; this is a default array until I incorporate MLUTs
		dm(md,0:32)=1. & dm(31:47)=2. & dm(48:63)=4. 		; this is a default array until I incorporate MLUTs
	endfor	


; Generate the gf tables (0-15 anode number, and energy for electrostatic attenuator)

	gf=fltarr(4,16,64)
;	geom = 



; Generate the attenuator tables (depends 0-3 state of attenuators and energy for electrostatic attenuator)

	mec = fltarr(16) & mec = [1.,1.,1.,.8,.3,.02,.02,.02,.02,.02,.04,.1,1.,1.,1.,1.]
	  
fm=1
if fm then begin
	att =fltarr(max_md+1,4,16,64) & att[*,*,*,*]=1.
	; startup mode SLUT=0, MLUT=0, 
;		att(*,0,*,*) = 0.1
;		att(*,1,*,*) = 0.1
;		att(*,2,*,*) = 1.
;		att(*,3,*,*) = 1.
endif else begin
endelse


; Generate the efficiency tables (should we assume eff is energy independent?)
	eff=fltarr(max_md+1,64)
fm=1
if fm then begin
	; startup mode SLUT=0, MLUT=0, 
		tof(0,*) = tof_def
		mas(0,*) = ((tof_def-t0)/dt)^2
		dm(0,0:32)=1. & dm(31:47)=2. & dm(48:63)=4.  ; accounts for variable TOF mass bin range
endif else begin
endelse


; print, 'made it this far'

;***************************************************************************************************************
;***************************************************************************************************************
;***************************************************************************************************************


if keyword_set(apids) then test = fix((total(apids eq 'c0') + total(apids eq 'C0')) < 1) else test=0
if not keyword_set(apids) or test then begin

	get_data,'mvn_STA_C0_DATA',data=t
	if size(/type,t) eq 8 then begin
		npts=128									; 64Ex1Dx1Ax2M
		ind1 = where(t.x gt 0, nn)
		tt=t.x[0:nn-1]

		get_data,'mvn_STA_C0_MODE',data=md
			md2 = md.y[0:nn-1] mod 128

		get_data,'mvn_STA_C0_AVG',data=cavg							; average state
		avg = 2^(cavg.y and replicate(7,nn))
		sum = (cavg.y and replicate(8,nn))/8
		comp = (cavg.y and replicate(192,nn))/128
		ind = where(avg[0:nn-2] ne avg[1:nn-1],count)					; corrected timing for averaging
		if count ne 0 then begin
			for i=0,count-1 do begin
				a0 = avg[ind[i]]
				a1 = avg[ind[i]+1]
				dt = a1*4.
				nf = abs(round( ( (tt[ind[i]+1]-tt[ind[i]])/4. -(a0+a1)/2. ) /(a1-a0) ))
				if nf gt 7 then nf=0
				for j=ind[i]-nf+1,ind[i] do tt[j]=tt[ind[i]+1]-dt*(ind[i]+1-j)
				for j=ind[i]-nf+1,ind[i] do md2[j]=md2[ind[i]+1]
			endfor
		endif

		store_data,'mvn_sta_C0_mode',data={x:tt,y:md2}					; corrected modes
			ylim,'mvn_sta_C0_mode',-1,6,0

		get_data,'mvn_STA_C0_ATTEN',data=catt						; attenuator state
		att1 = (catt.y[0:nn-1] and replicate(192,nn))/64
		att0 = att1
		store_data,'mvn_STA_C0_att_s',data={x:tt,y:att1}
		att2 = (catt.y[0:nn-1] and replicate(48,nn))/16
		store_data,'mvn_STA_C0_att_e',data={x:tt,y:att2}
		att3 = (catt.y[0:nn-1] and replicate(15,nn))
		store_data,'mvn_STA_C0_att_w',data={x:tt,y:att3}
		ind0 = where(att1[0:nn-2] ne att1[1:nn-1],cnt0)					; find att transitions
		if cnt0 ne 0 then begin
			nmax = 1024/npts
			for i=0,cnt0-1 do begin
				inds = ind0[i]
				if att3[inds] ge nmax then print,'APID C8 packet header error: WWWW too large, wwww=',att3[inds]
				idelay = fix(2/avg[inds])					; observed internal delays of 2 x 4sec, hence the "2"
				iww = att3[inds] - 1
				if iww lt 0 then iww = nmax-1
				att0[inds-nmax+1:inds] = att1[inds+1] 				; replace transition att w/ next value
				att0[inds-nmax+1:inds-nmax+1+iww+idelay] = att1[inds]		; use ww and interal delays to correct att at previous value
			endfor
		endif
		store_data,'mvn_sta_C0_att',data={x:tt,y:att0}
		ylim,'mvn_sta_C0_att',-1,4,0

		energy=nrg[md2,*]
		tmp=decomp19[transpose(reform(t.y[0:nn-1,*],nn,2,64),[0,2,1])]			; [time,en,ma]

		store_data,'mvn_sta_C0_P1A_E',data={x:tt,y:reform(total(tmp,3),nn,64),v:energy}
		store_data,'mvn_sta_C0_P1A_M',data={x:tt,y:reform(total(tmp,2),nn,2),v:findgen(2)}
		store_data,'mvn_sta_C0_P1A_E_M0',data={x:tt,y:reform(tmp[*,*,0],nn,64),v:energy}
		store_data,'mvn_sta_C0_P1A_E_M1',data={x:tt,y:reform(tmp[*,*,1],nn,64),v:energy}
		store_data,'mvn_sta_C0_P1A_tot',data={x:tt,y:total(total(tmp,2),2)}
		store_data,'mvn_sta_C0_P1A_E_M0_tot',data={x:tt,y:total(reform(tmp[*,*,0],nn,64),2),v:energy}
		store_data,'mvn_sta_C0_P1A_E_M1_tot',data={x:tt,y:total(reform(tmp[*,*,1],nn,64),2),v:energy}
		store_data,'mvn_sta_C0_P1A_all',data={x:tt,y:reform(tmp,nn,npts),v:findgen(npts)}


;		store_data,'mvn_sta_C0_E_eflx',data={x:tt,y:reform(total(tmp,3),nn,64),v:energy}

		ylim,'mvn_sta_C0_P1A_tot',0,0,1
		ylim,'mvn_sta_C0_P1A_E',.4,40000.,1
		ylim,'mvn_sta_C0_P1A_M',-1,2,0
		ylim,'mvn_sta_C0_P1A_E_M0',0,64,0
		ylim,'mvn_sta_C0_P1A_E_M1',0,64,0
		zlim,'mvn_sta_C0_P1A_E',1,1.e6,1
		zlim,'mvn_sta_C0_P1A_M',1,1.e6,1
		zlim,'mvn_sta_C0_P1A_E_M0',1,1.e6,1
		zlim,'mvn_sta_C0_P1A_E_M1',1,1.e6,1
		datagap=64.
		options,'mvn_sta_C0_P1A_E',datagap=datagap
		options,'mvn_sta_C0_P1A_M',datagap=datagap
		options,'mvn_sta_C0_P1A_E_M0',datagap=datagap
		options,'mvn_sta_C0_P1A_E_M1',datagap=datagap
		options,'mvn_sta_C0_P1A_tot',datagap=datagap
		options,'mvn_sta_C0_att',datagap=datagap
		options,'mvn_sta_C0_mode',datagap=datagap

		options,'mvn_sta_C0_P1A_E','spec',1
		options,'mvn_sta_C0_P1A_M','spec',1
		options,'mvn_sta_C0_P1A_E_M0','spec',1
		options,'mvn_sta_C0_P1A_E_M1','spec',1
		options,'mvn_sta_C0_P1A_E',ytitle='sta!CP1A-C0!C!CEnergy!C(eV)'
		options,'mvn_sta_C0_P1A_M',ytitle='sta!CP1A-C0!C!CMass!Cbin'
		options,'mvn_sta_C0_P1A_E_M0',ytitle='sta!CP1A-C0!CH+He+!C!CEnergy!Cbin'
		options,'mvn_sta_C0_P1A_E_M1',ytitle='sta!CP1A-C0!CO+O2+!C!CEnergy!Cbin'
	endif
endif

;***************************************************************************************************************
; APID C2

if keyword_set(apids) then test = fix((total(apids eq 'c2') + total(apids eq 'C2')) < 1) else test=0
if not keyword_set(apids) or test then begin

	get_data,'mvn_STA_C2_DATA',data=t
	if size(/type,t) eq 8 then begin
		npts=1024									; 32Ex1Dx1Ax32M
		ind1 = where(t.x gt 0, nn)
		tt=t.x[0:nn-1]

		get_data,'mvn_STA_C2_MODE',data=md
			md2 = md.y[0:nn-1] mod 128
			ylim,'mvn_STA_C2_MODE',-1,6,0

;	Correct tt for mode change boundaries, may want to add similar feature to other packets
;	assume start time in packet is correct and end time reflects incorrect CCxxSNNN bits
;	correct center time by using the start time of next packet

		get_data,'mvn_STA_C2_AVG',data=cavg							; average state
			avg = 2^(cavg.y[0:nn-1] and 7)
			sum = (cavg.y[0:nn-1] and 8)/8
		tt1 = tt - sum*2.*avg								; corrected timing for averaging
		tt2 = tt + sum*2.*avg
		if nn gt 2 then begin
			tt2[0:nn-2] = tt2[0:nn-2] < tt1[1:nn-1]
			tt3 = (tt1 + tt2)/2.
;			print,tt3-tt
			tt = tt3
		endif

		get_data,'mvn_STA_C2_ATTEN',data=catt						; attenuator state
		att0 = (catt.y[0:nn-1] and replicate(192,nn))/64
		if nn gt 2 then begin
			ind0 = where(att0[0:nn-2] ne att0[1:nn-1],cnt0)
			if cnt0 ne 0 then begin
				for i=0,cnt0-1 do begin
					inds = ind0[i]
					idelay = fix(2/avg[inds])
					if idelay gt 0 then att0[inds+1:(inds+idelay)<(nn-1)] = att0[inds]
				endfor
			endif
		endif
		store_data,'mvn_sta_C2_att',data={x:tt,y:att0}
		ylim,'mvn_sta_C2_att',-1,4,0

		get_data,'mvn_STA_C2_SEQ_CNTR',data=tmp7
		store_data,'mvn_sta_C2_SEQ_CNTR2',data={x:tt,y:tmp7.y[0:nn-1]}

		energy=nrg[md2,*]
		energy=total(reform(energy,nn,2,32),2)/2.
		mass=mas[md2,*]
		mass=total(reform(mass,nn,2,32),2)/2.
		dmass=dm[md2,*]
		dmass=total(reform(dmass,nn,2,32),2)
		tmp=decomp19[transpose(reform(t.y[0:nn-1,*],nn,32,32),[0,2,1])]			; [time,en,ma]

		store_data,'mvn_sta_C2_P1B_Ebin',data={x:tt,y:total(tmp,3),v:findgen(32)}
		store_data,'mvn_sta_C2_P1B_Mbin',data={x:tt,y:total(tmp,2),v:findgen(32)}
		store_data,'mvn_sta_C2_P1B_E',data={x:tt,y:total(tmp,3),v:energy}
		store_data,'mvn_sta_C2_P1B_M',data={x:tt,y:total(tmp,2),v:mass}
		store_data,'mvn_sta_C2_P1B_tot',data={x:tt,y:total(total(tmp,2),2)}
		store_data,'mvn_sta_C2_P1B_all',data={x:tt,y:reform(tmp,nn,npts),v:findgen(npts)}

		ylim,'mvn_sta_C2_P1B_tot',0,0,1
		ylim,'mvn_sta_C2_P1B_Ebin',0,32,0
		ylim,'mvn_sta_C2_P1B_Mbin',0,32,0
		ylim,'mvn_sta_C2_P1B_E',.4,40000.,1
		ylim,'mvn_sta_C2_P1B_M',.2,100,1
		zlim,'mvn_sta_C2_P1B_Ebin',1,1.e6,1
		zlim,'mvn_sta_C2_P1B_Mbin',1,1.e6,1
		zlim,'mvn_sta_C2_P1B_E',1,1.e6,1
		zlim,'mvn_sta_C2_P1B_M',1,1.e6,1
		datagap=256.
		options,'mvn_sta_C2_P1B_Ebin',datagap=datagap
		options,'mvn_sta_C2_P1B_Mbin',datagap=datagap
		options,'mvn_sta_C2_P1B_E',datagap=datagap
		options,'mvn_sta_C2_P1B_M',datagap=datagap
		options,'mvn_sta_C2_P1B_tot',datagap=datagap
		options,'mvn_sta_C2_att',datagap=datagap
		options,'mvn_STA_C2_MODE',datagap=datagap


		options,'mvn_sta_C2_P1B_Ebin','spec',1
		options,'mvn_sta_C2_P1B_Mbin','spec',1
		options,'mvn_sta_C2_P1B_E','spec',1
		options,'mvn_sta_C2_P1B_M','spec',1
		options,'mvn_sta_C2_P1B_Ebin',ytitle='sta!CP1B-C2!C!CEnergy!Cbin'
		options,'mvn_sta_C2_P1B_Mbin',ytitle='sta!CP1B-C2!C!CMass!Cbin'
		options,'mvn_sta_C2_P1B_E',ytitle='sta!CP1B-C2!C!CEnergy!C(eV)'
		options,'mvn_sta_C2_P1B_M',ytitle='sta!CP1B-C2!C!CMass!C(amu)'
	endif
endif

;***************************************************************************************************************

if keyword_set(apids) then test = fix((total(apids eq 'c4') + total(apids eq 'C4')) < 1) else test=0
if not keyword_set(apids) or test then begin

	get_data,'mvn_STA_C4_DATA',data=t
	if size(/type,t) eq 8 then begin
		npts=256									; 4Ex1Dx1Ax64M
		ind1 = where(t.x gt 0, nn)
		tt=t.x[0:nn-1]

		get_data,'mvn_STA_C4_MODE',data=md
			md2 = md.y[0:nn-1] mod 128

		get_data,'mvn_STA_C4_AVG',data=cavg							; set avg state
		avg = 2^(cavg.y and replicate(7,nn))
		sum = (cavg.y and replicate(8,nn))/8
		comp = (cavg.y and replicate(192,nn))/128
		ind = where(avg[0:nn-2] ne avg[1:nn-1],count)
		if count ne 0 then begin
			for i=0,count-1 do begin
				a0 = avg[ind[i]]
				a1 = avg[ind[i]+1]
				dt = a1*4.
				nf = abs(round( ( (tt[ind[i]+1]-tt[ind[i]])/4. -(a0+a1)/2. ) /(a1-a0) ))
				if nf gt 15 then nf=0
				for j=ind[i]-nf+1,ind[i] do tt[j]=tt[ind[i]+1]-dt*(ind[i]+1-j)
			endfor
		endif

		store_data,'mvn_sta_C4_mode',data={x:tt,y:md2}
			ylim,'mvn_sta_C4_mode',-1,6,0

		get_data,'mvn_STA_C4_ATTEN',data=catt						; set attenuator state
		att1 = (catt.y[0:nn-1] and replicate(192,nn))/64
		att0 = att1
		store_data,'mvn_STA_C4_att_s',data={x:tt,y:att1}
		att2 = (catt.y[0:nn-1] and replicate(48,nn))/16
		store_data,'mvn_STA_C4_att_e',data={x:tt,y:att2}
		att3 = (catt.y[0:nn-1] and replicate(15,nn))
		store_data,'mvn_STA_C4_att_w',data={x:tt,y:att3}
		ind0 = where(att1[0:nn-2] ne att1[1:nn-1],cnt0)					; find att transitions
		if cnt0 ne 0 then begin
			nmax = 1024/npts
			for i=0,cnt0-1 do begin
				inds = ind0[i]
				if att3[inds] ge nmax then print,'APID C8 packet header error: WWWW too large, wwww=',att3[inds]
				idelay = fix(2/avg[inds])					; observed internal delays of 2 x 4sec, hence the "2"
				iww = att3[inds] - 1
				if iww lt 0 then iww = nmax-1
				att0[inds-nmax+1:inds] = att1[inds+1] 				; replace transition att w/ next value
				att0[inds-nmax+1:inds-nmax+1+iww+idelay] = att1[inds]		; use ww and interal delays to correct att at previous value
			endfor
		endif
		store_data,'mvn_sta_C4_att',data={x:tt,y:att0}
		ylim,'mvn_sta_C4_att',-1,4,0

		energy=nrg[md2,*]
		energy=total(reform(energy,nn,16,4),2)/16.
		tmp=decomp19[transpose(reform(t.y[0:nn-1,*],nn,64,4),[0,2,1])]			; [time,en,ma]

		store_data,'mvn_sta_C4_P1C_E',data={x:tt,y:total(tmp,3),v:energy}
		store_data,'mvn_sta_C4_P1C_M',data={x:tt,y:total(tmp,2),v:indgen(64)}
		store_data,'mvn_sta_C4_P1C_tot',data={x:tt,y:total(total(tmp,2),2)}
		store_data,'mvn_sta_C4_P1C_all',data={x:tt,y:reform(tmp,nn,npts),v:findgen(npts)}

		ylim,'mvn_sta_C4_P1C_tot',0,0,1
		ylim,'mvn_sta_C4_P1C_E',.4,40000.,1
		ylim,'mvn_sta_C4_P1C_M',0,64,0
		zlim,'mvn_sta_C4_P1C_E',1,1.e6,1
		zlim,'mvn_sta_C4_P1C_M',1,1.e6,1
		datagap=64.
		options,'mvn_sta_C4_P1C_E',datagap=datagap
		options,'mvn_sta_C4_P1C_M',datagap=datagap
		options,'mvn_sta_C4_P1C_tot',datagap=datagap
		options,'mvn_sta_C4_att',datagap=datagap
		options,'mvn_sta_C4_mode',datagap=datagap

		options,'mvn_sta_C4_P1C_E','spec',1
		options,'mvn_sta_C4_P1C_M','spec',1
		options,'mvn_sta_C4_P1C_E',ytitle='sta!CP1C-C4!C!CEnergy!C(eV)'
		options,'mvn_sta_C4_P1C_M',ytitle='sta!CP1C-C4!C!CMass!Cbin'
	endif
endif


;***************************************************************************************************************

if keyword_set(apids) then test = fix((total(apids eq 'c6') + total(apids eq 'C6')) < 1) else test=0
if not keyword_set(apids) or test then begin

	get_data,'mvn_STA_C6_DATA',data=t
	if size(/type,t) eq 8 then begin
		npts=2048									; 32Ex1Dx1Ax64M
		np = 2										; np is number of packets per measurement
		ind1 = where(t.x gt 0, nn)
		tt=t.x[0:nn-1]

		get_data,'mvn_STA_C6_DIAG',data=diag
		en = (diag.y[0:nn-1] AND 15)
		dt=0.1
			store_data,'mvn_sta_C6_DIAG_EN',data={x:tt+en*dt,y:en}
			get_data,'mvn_STA_C6_SEQ_CNTR',data=tmp2
			store_data,'mvn_sta_C6_SEQ_CNTR_EN',data={x:tt+en*dt,y:tmp2.y[0:nn-1]}
			store_data,'mvn_sta_C6_DATA_EN',data={x:tt+en*dt,y:total(t.y[0:nn-1,*],2)}

		ind = where(en eq 0,ndis)							; ndis is number of complete distributions
		if ind[ndis-1]+np-1 gt nn-1 then ndis=ndis-1					; eliminate last distribution if not complete


		if ndis ge 1 then begin

			dat = t.y[ind[0]:ind[0]+ndis*np-1,*]
			tdis = tt[ind[0:ndis-1]]

			get_data,'mvn_STA_C6_MODE',data=md
				md2 = md.y[ind[0:ndis-1]] mod 128
			get_data,'mvn_STA_C6_AVG',data=cavg
			avg = 2^(cavg.y[ind[0:ndis-1]] and 7)
			sum = (cavg.y[ind[0:ndis-1]] and 8)/8

			get_data,'mvn_STA_C6_ATTEN',data=catt						; attenuator state
			att0 = (catt.y[ind[0:ndis-1]] and replicate(192,ndis))/64
			if ndis ge 2 then ind0 = where(att0[0:ndis-2] ne att0[1:ndis-1],cnt0) else cnt0=0
			if cnt0 ne 0 then begin
				for i=0,cnt0-1 do begin
					inds = ind0[i]
					idelay = fix(2/avg[inds])
					if idelay gt 0 then att0[inds+1:(inds+idelay)<(nn-1)] = att0[inds]
				endfor
			endif
			store_data,'mvn_sta_C6_att',data={x:tdis,y:att0}
			ylim,'mvn_sta_C6_att',-1,4,0


			energy=nrg[md2,*]
			energy=total(reform(energy,ndis,2,32),2)/2.
			tmp=decomp19[reform(transpose(reform(dat,np,ndis,64,16),[1,3,0,2]),ndis,32,64)]		; [time,en,def,ma]

			store_data,'mvn_sta_C6_P1D_E',data={x:tdis,y:total(tmp,3),v:energy}
			store_data,'mvn_sta_C6_P1D_M',data={x:tdis,y:total(tmp,2),v:indgen(64)}
			store_data,'mvn_sta_C6_P1D_tot',data={x:tdis,y:total(total(tmp,3),2)}
			store_data,'mvn_sta_C6_P1D_all',data={x:tdis,y:reform(tmp,ndis,npts),v:indgen(npts)}

			ylim,'mvn_sta_C6_P1D_tot',0,0,1
			ylim,'mvn_sta_C6_P1D_E',.4,40000.,1
			ylim,'mvn_sta_C6_P1D_M',-1,64,0
			zlim,'mvn_sta_C6_P1D_E',1,1.e6,1
			zlim,'mvn_sta_C6_P1D_M',1,1.e6,1
			options,'mvn_sta_C6_P1D_E',datagap=512.
			options,'mvn_sta_C6_P1D_M',datagap=512.
			options,'mvn_sta_C6_P1D_tot',datagap=512.
			options,'mvn_sta_C6_P1D_E','spec',1
			options,'mvn_sta_C6_P1D_M','spec',1
			options,'mvn_sta_C6_P1D_E',ytitle='sta!CP1D-C6!C!CEnergy!C(eV)'
			options,'mvn_sta_C6_P1D_M',ytitle='sta!CP1D-C6!C!CMass!Cbin'
		endif
	endif
endif


;***************************************************************************************************************

if keyword_set(apids) then test = fix((total(apids eq 'c8') + total(apids eq 'C8')) < 1) else test=0
if not keyword_set(apids) or test then begin

	get_data,'mvn_STA_C8_DATA',data=t
	if size(/type,t) eq 8 then begin
		npts=512									; 32Ex16Dx1Ax1M
		ind1 = where(t.x gt 0, nn)
		tt=t.x[0:nn-1]

		get_data,'mvn_STA_C8_MODE',data=md
			md2 = md.y[0:nn-1] mod 128
		get_data,'mvn_STA_C8_AVG',data=cavg
		avg = 2^(cavg.y and replicate(7,nn))
		sum = (cavg.y and replicate(8,nn))/8
		comp = (cavg.y and replicate(192,nn))/128
		ind = where(avg[0:nn-2] ne avg[1:nn-1],count)
		if count ne 0 then begin
			for i=0,count-1 do begin
				a0 = avg[ind[i]]
				a1 = avg[ind[i]+1]
				dt = a1*4.
				nf = abs(round( ( (tt[ind[i]+1]-tt[ind[i]])/4. -(a0+a1)/2. ) /(a1-a0) ))
				if nf gt 1 then nf=0
				for j=ind[i]-nf+1,ind[i] do tt[j]=tt[ind[i]+1]-dt*(ind[i]+1-j)
			endfor
		endif
		store_data,'mvn_sta_C8_mode',data={x:tt,y:md2}
			ylim,'mvn_sta_C8_mode',-1,6,0

		get_data,'mvn_STA_C8_ATTEN',data=catt						; set attenuator state
		att1 = (catt.y[0:nn-1] and replicate(192,nn))/64
		att0 = att1
		store_data,'mvn_STA_C8_att_s',data={x:tt,y:att1}
		att2 = (catt.y[0:nn-1] and replicate(48,nn))/16
		store_data,'mvn_STA_C8_att_e',data={x:tt,y:att2}
		att3 = (catt.y[0:nn-1] and replicate(15,nn))
		store_data,'mvn_STA_C8_att_w',data={x:tt,y:att3}
		ind0 = where(att1[0:nn-2] ne att1[1:nn-1],cnt0)					; find att transitions
		if cnt0 ne 0 then begin
			nmax = 1024/npts
			for i=0,cnt0-1 do begin
				inds = ind0[i]
				if att3[inds] ge nmax then print,'APID C8 packet header error: WWWW too large, wwww=',att3[inds]
				idelay = fix(2/avg[inds])					; observed internal delays of 2 x 4sec, hence the "2"
				iww = att3[inds] - 1
				if iww lt 0 then iww = nmax-1
				att0[inds-nmax+1:inds] = att1[inds+1] 				; replace transition att w/ next value
				att0[inds-nmax+1:inds-nmax+1+iww+idelay] = att1[inds]		; use ww and interal delays to correct att at previous value
			endfor
		endif
		store_data,'mvn_sta_C8_att',data={x:tt,y:att0}
		ylim,'mvn_sta_C8_att',-1,4,0

		energy=nrg[md2,*]
		energy=total(reform(energy,nn,2,32),2)/2.
		tmp=decomp19[transpose(reform(t.y[0:nn-1,*],nn,16,32),[0,2,1])]			; [time,en,def]

		store_data,'mvn_sta_C8_P2_E',data={x:tt,y:total(tmp,3),v:energy}
		store_data,'mvn_sta_C8_P2_D',data={x:tt,y:total(tmp,2),v:indgen(16)}
		store_data,'mvn_sta_C8_P2_tot',data={x:tt,y:total(total(tmp,2),2)}
		store_data,'mvn_sta_C8_P2_all',data={x:tt,y:reform(tmp,nn,npts),v:findgen(npts)}

		ylim,'mvn_sta_C8_P2_tot',0,0,1
		ylim,'mvn_sta_C8_P2_E',.4,40000.,1
		ylim,'mvn_sta_C8_P2_D',0,16,0
		zlim,'mvn_sta_C8_P2_E',1,1.e6,1
		zlim,'mvn_sta_C8_P2_D',1,1.e6,1
		options,'mvn_sta_C8_P2_E',datagap=64.
		options,'mvn_sta_C8_P2_D',datagap=64.
		options,'mvn_sta_C8_P2_tot',datagap=64.
		options,'mvn_sta_C8_P2_E','spec',1
		options,'mvn_sta_C8_P2_D','spec',1
		options,'mvn_sta_C8_P2_E',ytitle='sta!CP2-C8!C!CEnergy!C(eV)'
		options,'mvn_sta_C8_P2_D',ytitle='sta!CP2-C8!C!CDef!Cbin'
	endif
endif


;***************************************************************************************************************

if keyword_set(apids) then test = fix((total(apids eq 'ca') + total(apids eq 'CA')) < 1) else test=0
if not keyword_set(apids) or test then begin

	get_data,'mvn_STA_CA_DATA',data=t
	if size(/type,t) eq 8 then begin
		npts=1024									; 16Ex4Dx16Ax1M
		ind1 = where(t.x gt 0, nn)
		tt=t.x[0:nn-1]

		get_data,'mvn_STA_CA_MODE',data=md
			md2 = md.y[0:nn-1] mod 128

		get_data,'mvn_STA_CA_AVG',data=cavg							; average state
			avg = 2^(cavg.y[0:nn-1] and 7)
			sum = (cavg.y[0:nn-1] and 8)/8
		tt1 = tt - sum*2.*avg								; corrected timing for averaging
		tt2 = tt + sum*2.*avg

		if nn gt 2 then begin
			tt2[0:nn-2] = tt2[0:nn-2] < tt1[1:nn-1]
			tt3 = (tt1 + tt2)/2.
			tt = tt3
		endif


		get_data,'mvn_STA_CA_ATTEN',data=catt						; attenuator state
		att0 = (catt.y[0:nn-1] and replicate(192,nn))/64
		ind0 = where(att0[0:nn-2] ne att0[1:nn-1],cnt0)
		if cnt0 ne 0 then begin
			for i=0,cnt0-1 do begin
				inds = ind0[i]
				idelay = fix(2/avg[inds])
				if idelay gt 0 then att0[inds+1:(inds+idelay)<(nn-1)] = att0[inds]
			endfor
		endif
		store_data,'mvn_sta_CA_att',data={x:tt,y:att0}
		ylim,'mvn_sta_CA_att',-1,4,0

		energy=nrg[md2,*]
		energy=total(reform(energy,nn,4,16),2)/4.
		tmp=decomp19[transpose(reform(t.y[0:nn-1,*],nn,16,4,16),[0,3,2,1])]		; [time,en,def,an]

		store_data,'mvn_sta_CA_P3_E',data={x:tt,y:total(total(tmp,4),3),v:energy}
		store_data,'mvn_sta_CA_P3_D',data={x:tt,y:total(total(tmp,4),2),v:indgen(4)}
		store_data,'mvn_sta_CA_P3_A',data={x:tt,y:total(total(tmp,3),2),v:indgen(16)}
		store_data,'mvn_sta_CA_P3_tot',data={x:tt,y:total(total(total(tmp,2),2),2)}
		store_data,'mvn_sta_CA_P3_all',data={x:tt,y:reform(tmp,nn,npts),v:findgen(npts)}

		ylim,'mvn_sta_CA_P3_tot',0,0,1
		ylim,'mvn_sta_CA_P3_E',.4,40000.,1
		ylim,'mvn_sta_CA_P3_D',-1,4,0
		ylim,'mvn_sta_CA_P3_A',-1,16,0
		zlim,'mvn_sta_CA_P3_E',1,1.e6,1
		zlim,'mvn_sta_CA_P3_D',1,1.e6,1
		zlim,'mvn_sta_CA_P3_A',1,1.e6,1
		options,'mvn_sta_CA_P3_E',datagap=64.
		options,'mvn_sta_CA_P3_D',datagap=64.
		options,'mvn_sta_CA_P3_A',datagap=64.
		options,'mvn_sta_CA_P3_tot',datagap=64.
		options,'mvn_sta_CA_P3_E','spec',1
		options,'mvn_sta_CA_P3_D','spec',1
		options,'mvn_sta_CA_P3_A','spec',1
		options,'mvn_sta_CA_P3_E',ytitle='sta!CP3-CA!C!CEnergy!C(eV)'
		options,'mvn_sta_CA_P3_D',ytitle='sta!CP3-CA!C!CDef!Cbin'
		options,'mvn_sta_CA_P3_A',ytitle='sta!CP3-CA!C!CAnode!Cbin'
	endif
endif


;***************************************************************************************************************

if keyword_set(apids) then test = fix((total(apids eq 'cc') + total(apids eq 'CC')) < 1) else test=0
if not keyword_set(apids) or test then begin

	get_data,'mvn_STA_CC_DATA',data=t
	if size(/type,t) eq 8 then begin
		npts=8192									; 32Ex8Dx1Ax32M
		np = 8										; np is number of packets per measurement
		ind1 = where(t.x gt 0, nn)
		tt=t.x[0:nn-1]

		get_data,'mvn_STA_CC_DIAG',data=diag
		en = (diag.y[0:nn-1] AND 15)
		dt=0.1
			store_data,'mvn_sta_CC_DIAG_EN',data={x:tt+en*dt,y:en}
			get_data,'mvn_STA_CC_SEQ_CNTR',data=tmp2
			store_data,'mvn_sta_CC_SEQ_CNTR_EN',data={x:tt+en*dt,y:tmp2.y[0:nn-1]}
			store_data,'mvn_sta_CC_DATA_EN',data={x:tt+en*dt,y:total(t.y[0:nn-1,*],2)}

		ind = where(en eq 0,ndis)							; ndis is number of complete distributions
		if ind[ndis-1]+np-1 gt nn-1 then ndis=ndis-1					; eliminate last distribution if not complete

		if ndis ge 1 then begin

			dat = t.y[ind[0]:ind[0]+ndis*np-1,*]
			tdis = tt[ind[0:ndis-1]]

			get_data,'mvn_STA_CC_MODE',data=md
				md2 = md.y[ind[0:ndis-1]] mod 128
			get_data,'mvn_STA_CC_AVG',data=cavg
			avg = 2^(cavg.y[ind[0:ndis-1]] and 7)
			sum = (cavg.y[ind[0:ndis-1]] and 8)/8

			get_data,'mvn_STA_CC_ATTEN',data=catt						; attenuator state
			att0 = (catt.y[ind[0:ndis-1]] and replicate(192,ndis))/64
			if ndis ge 2 then ind0 = where(att0[0:ndis-2] ne att0[1:ndis-1],cnt0) else cnt0=0
			if cnt0 ne 0 then begin
				for i=0,cnt0-1 do begin
					inds = ind0[i]
					idelay = fix(2/avg[inds])
					if idelay gt 0 then att0[inds+1:(inds+idelay)<(nn-1)] = att0[inds]
				endfor
			endif
			store_data,'mvn_sta_CC_att',data={x:tdis,y:att0}
			ylim,'mvn_sta_CC_att',-1,4,0


			energy=nrg[md2,*]
			energy=total(reform(energy,ndis,2,32),2)/2.
			tmp=decomp19[reform(transpose(reform(dat,np,ndis,32,8,4),[1,4,0,3,2]),ndis,32,8,32)]		; [time,en,def,ma]

			store_data,'mvn_sta_CC_P4A_E',data={x:tdis,y:total(total(tmp,4),3),v:energy}
			store_data,'mvn_sta_CC_P4A_D',data={x:tdis,y:total(total(tmp,4),2),v:indgen(8)}
			store_data,'mvn_sta_CC_P4A_M',data={x:tdis,y:total(total(tmp,2),2),v:indgen(32)}
			store_data,'mvn_sta_CC_P4A_tot',data={x:tdis,y:total(total(total(tmp,4),3),2)}
			store_data,'mvn_sta_CC_P4A_all',data={x:tdis,y:reform(tmp,ndis,npts),v:indgen(npts)}

			ylim,'mvn_sta_CC_P4A_tot',0,0,1
			ylim,'mvn_sta_CC_P4A_E',.4,40000.,1
			ylim,'mvn_sta_CC_P4A_D',-1,8,0
			ylim,'mvn_sta_CC_P4A_M',-1,32,0
			zlim,'mvn_sta_CC_P4A_E',1,1.e6,1
			zlim,'mvn_sta_CC_P4A_D',1,1.e6,1
			zlim,'mvn_sta_CC_P4A_M',1,1.e6,1
			options,'mvn_sta_CC_P4A_E',datagap=128.
			options,'mvn_sta_CC_P4A_D',datagap=128.
			options,'mvn_sta_CC_P4A_M',datagap=128.
			options,'mvn_sta_CC_P4A_tot',datagap=128.
			options,'mvn_sta_CC_P4A_E','spec',1
			options,'mvn_sta_CC_P4A_D','spec',1
			options,'mvn_sta_CC_P4A_M','spec',1
			options,'mvn_sta_CC_P4A_E',ytitle='sta!CP4A-CC!C!CEnergy!C(eV)'
			options,'mvn_sta_CC_P4A_D',ytitle='sta!CP4A-CC!C!CDef!Cbin'
			options,'mvn_sta_CC_P4A_M',ytitle='sta!CP4A-CC!C!Mass!Cbin'
		endif
	endif
endif


;***************************************************************************************************************

if keyword_set(apids) then test = fix((total(apids eq 'cd') + total(apids eq 'CD')) < 1) else test=0
if not keyword_set(apids) or test then begin

	get_data,'mvn_STA_CD_DATA',data=t
	if size(/type,t) eq 8 then begin
		npts=8192									; 32Ex8Dx1Ax32M
		np = 8										; np is number of packets per measurement
		ind1 = where(t.x gt 0, nn)
		tt=t.x[0:nn-1]

		get_data,'mvn_STA_CD_DIAG',data=diag
		en = (diag.y[0:nn-1] AND 15)
		dt=0.1
			store_data,'mvn_sta_CD_DIAG_EN',data={x:tt+en*dt,y:en}
			get_data,'mvn_STA_CD_SEQ_CNTR',data=tmp2
			store_data,'mvn_sta_CD_SEQ_CNTR_EN',data={x:tt+en*dt,y:tmp2.y[0:nn-1]}
			store_data,'mvn_sta_CD_DATA_EN',data={x:tt+en*dt,y:total(t.y[0:nn-1,*],2)}

		ind = where(en eq 0,ndis)							; ndis is number of complete distributions
		if ind[ndis-1]+np-1 gt nn-1 then ndis=ndis-1					; eliminate last distribution if not complete

		if ndis ge 1 then begin

			dat = t.y[ind[0]:ind[0]+ndis*np-1,*]
			tdis = tt[ind[0:ndis-1]]

			get_data,'mvn_STA_CD_MODE',data=md
				md2 = md.y[ind[0:ndis-1]] mod 128
			get_data,'mvn_STA_CD_AVG',data=cavg
			avg = 2^(cavg.y[ind[0:ndis-1]] and 7)
			sum = (cavg.y[ind[0:ndis-1]] and 8)/8

			get_data,'mvn_STA_CD_ATTEN',data=catt						; attenuator state
			att0 = (catt.y[ind[0:ndis-1]] and replicate(192,ndis))/64
			if ndis ge 2 then ind0 = where(att0[0:ndis-2] ne att0[1:ndis-1],cnt0) else cnt0=0
			if cnt0 ne 0 then begin
				for i=0,cnt0-1 do begin
					inds = ind0[i]
					idelay = fix(2/avg[inds])
					if idelay gt 0 then att0[inds+1:(inds+idelay)<(nn-1)] = att0[inds]
				endfor
			endif
			store_data,'mvn_sta_CD_att',data={x:tdis,y:att0}
			ylim,'mvn_sta_CD_att',-1,4,0


			energy=nrg[md2,*]
			energy=total(reform(energy,ndis,2,32),2)/2.
			tmp=decomp19[reform(transpose(reform(dat,np,ndis,32,8,4),[1,4,0,3,2]),ndis,32,8,32)]		; [time,en,def,ma]

			store_data,'mvn_sta_CD_P4A_E',data={x:tdis,y:total(total(tmp,4),3),v:energy}
			store_data,'mvn_sta_CD_P4A_D',data={x:tdis,y:total(total(tmp,4),2),v:indgen(8)}
			store_data,'mvn_sta_CD_P4A_M',data={x:tdis,y:total(total(tmp,2),2),v:indgen(32)}
			store_data,'mvn_sta_CD_P4A_tot',data={x:tdis,y:total(total(total(tmp,4),3),2)}
			store_data,'mvn_sta_CD_P4A_all',data={x:tdis,y:reform(tmp,ndis,npts),v:indgen(npts)}

			ylim,'mvn_sta_CD_P4A_tot',0,0,1
			ylim,'mvn_sta_CD_P4A_E',.4,40000.,1
			ylim,'mvn_sta_CD_P4A_D',-1,8,0
			ylim,'mvn_sta_CD_P4A_M',-1,32,0
			zlim,'mvn_sta_CD_P4A_E',1,1.e6,1
			zlim,'mvn_sta_CD_P4A_D',1,1.e6,1
			zlim,'mvn_sta_CD_P4A_M',1,1.e6,1
			options,'mvn_sta_CD_P4A_E',datagap=16.
			options,'mvn_sta_CD_P4A_D',datagap=16.
			options,'mvn_sta_CD_P4A_M',datagap=16.
			options,'mvn_sta_CD_P4A_tot',datagap=16.
			options,'mvn_sta_CD_P4A_E','spec',1
			options,'mvn_sta_CD_P4A_D','spec',1
			options,'mvn_sta_CD_P4A_M','spec',1
			options,'mvn_sta_CD_P4A_E',ytitle='sta!CP4A-CD!C!CEnergy!C(eV)'
			options,'mvn_sta_CD_P4A_D',ytitle='sta!CP4A-CD!C!CDef!Cbin'
			options,'mvn_sta_CD_P4A_M',ytitle='sta!CP4A-CD!C!CMass!Cbin'
		endif
	endif
endif


;***************************************************************************************************************

if keyword_set(apids) then test = fix((total(apids eq 'ce') + total(apids eq 'CE')) < 1) else test=0
if not keyword_set(apids) or test then begin

	get_data,'mvn_STA_CE_DATA',data=t
	if size(/type,t) eq 8 then begin
		npts=16384									; 16Ex4Dx16Ax16M
		np = 16										; np is number of packets per measurement, npts/1024
		ind1 = where(t.x gt 0, nn)
		tt=t.x[0:nn-1]

		get_data,'mvn_STA_CE_DIAG',data=diag
		en = (diag.y[0:nn-1] AND 15)
		dt=0.1
			store_data,'mvn_sta_CE_DIAG_EN',data={x:tt+en*dt,y:en}
			get_data,'mvn_STA_CE_SEQ_CNTR',data=tmp2
			store_data,'mvn_sta_CE_SEQ_CNTR_EN',data={x:tt+en*dt,y:tmp2.y[0:nn-1]}
			store_data,'mvn_sta_CE_DATA_EN',data={x:tt+en*dt,y:total(t.y[0:nn-1,*],2)}

		ind = where(en eq 0,ndis)							; ndis is number of complete distributions
		if ind[ndis-1]+np-1 gt nn-1 then ndis=ndis-1					; eliminate last distribution if not complete

		if ndis ge 1 then begin		
		   if ndis*np gt nn then begin
			print,'Error - APID CE diag bits in header set incorrectly
				print,'data array'
				help,t,/st
				print,'# packets=',nn
				print,'diag bits=',en
				print,'diag bits should cycle from 0 to 16, over ',nn,' packets'
		   endif else begin

			dat = t.y[ind[0]:ind[0]+ndis*np-1,*]
			tdis = tt[ind[0:ndis-1]]

			get_data,'mvn_STA_CE_MODE',data=md
				md2 = md.y[ind[0:ndis-1]] mod 128
			get_data,'mvn_STA_CE_AVG',data=cavg
			avg = 2^(cavg.y[ind[0:ndis-1]] and 7)
			sum = (cavg.y[ind[0:ndis-1]] and 8)/8

			get_data,'mvn_STA_CE_ATTEN',data=catt						; attenuator state
			att0 = (catt.y[ind[0:ndis-1]] and replicate(192,ndis))/64
			if ndis ge 2 then ind0 = where(att0[0:ndis-2] ne att0[1:ndis-1],cnt0) else cnt0=0
			if cnt0 ne 0 then begin
				for i=0,cnt0-1 do begin
					inds = ind0[i]
					idelay = fix(2/avg[inds])
					if idelay gt 0 then att0[inds+1:(inds+idelay)<(nn-1)] = att0[inds]
				endfor
			endif
			store_data,'mvn_sta_CE_att',data={x:tdis,y:att0}
			ylim,'mvn_sta_CE_att',-1,4,0


			energy=nrg[md2,*]
			energy=total(reform(energy,ndis,4,16),2)/4.
			tmp=decomp19[transpose(reform(dat,np,ndis,16,16,4),[1,0,4,2,3])]	; [time,en,def,an,ma]

			store_data,'mvn_sta_CE_P4B_E',data={x:tdis,y:total(total(total(tmp,5),4),3),v:energy}
			store_data,'mvn_sta_CE_P4B_D',data={x:tdis,y:total(total(total(tmp,5),4),2),v:indgen(4)}
			store_data,'mvn_sta_CE_P4B_A',data={x:tdis,y:total(total(total(tmp,5),3),2),v:indgen(16)}
			store_data,'mvn_sta_CE_P4B_M',data={x:tdis,y:total(total(total(tmp,4),3),2),v:indgen(16)}
			store_data,'mvn_sta_CE_P4B_tot',data={x:tdis,y:total(total(total(total(tmp,5),4),3),2)}
			store_data,'mvn_sta_CE_P4B_all',data={x:tdis,y:reform(tmp,ndis,npts),v:indgen(npts)}

			ylim,'mvn_sta_CE_P4B_tot',0,0,1
			ylim,'mvn_sta_CE_P4B_E',.4,40000.,1
			ylim,'mvn_sta_CE_P4B_D',-1,4,0
			ylim,'mvn_sta_CE_P4B_A',-1,16,0
			ylim,'mvn_sta_CE_P4B_M',-1,16,0
			zlim,'mvn_sta_CE_P4B_E',1,1.e6,1
			zlim,'mvn_sta_CE_P4B_D',1,1.e6,1
			zlim,'mvn_sta_CE_P4B_A',1,1.e6,1
			zlim,'mvn_sta_CE_P4B_M',1,1.e6,1
			options,'mvn_sta_CE_P4B_E',datagap=256.
			options,'mvn_sta_CE_P4B_D',datagap=256.
			options,'mvn_sta_CE_P4B_A',datagap=256.
			options,'mvn_sta_CE_P4B_M',datagap=256.
			options,'mvn_sta_CE_P4B_tot',datagap=256.
			options,'mvn_sta_CE_P4B_E','spec',1
			options,'mvn_sta_CE_P4B_D','spec',1
			options,'mvn_sta_CE_P4B_A','spec',1
			options,'mvn_sta_CE_P4B_M','spec',1
			options,'mvn_sta_CE_P4B_E',ytitle='sta!CP4B-CE!C!CEnergy!C(eV)'
			options,'mvn_sta_CE_P4B_D',ytitle='sta!CP4B-CE!C!CDef!Cbin'
			options,'mvn_sta_CE_P4B_A',ytitle='sta!CP4B-CE!C!CAnode!Cbin'
			options,'mvn_sta_CE_P4B_M',ytitle='sta!CP4B-CE!C!CMass!Cbin'

		   endelse
		endif
	endif
endif


;***************************************************************************************************************

if keyword_set(apids) then test = fix((total(apids eq 'cf') + total(apids eq 'CF')) < 1) else test=0
if not keyword_set(apids) or test then begin

	get_data,'mvn_STA_CF_DATA',data=t
	if size(/type,t) eq 8 then begin
		npts=16384									; 16Ex4Dx16Ax16M
		np = 16										; np is number of packets per measurement, npts/1024
		ind1 = where(t.x gt 0, nn)
		tt=t.x[0:nn-1]

		get_data,'mvn_STA_CF_DIAG',data=diag
		en = (diag.y[0:nn-1] AND 15)
		dt=0.1
			store_data,'mvn_sta_CF_DIAG_EN',data={x:tt+en*dt,y:en}
			get_data,'mvn_STA_CF_SEQ_CNTR',data=tmp2
			store_data,'mvn_sta_CF_SEQ_CNTR_EN',data={x:tt+en*dt,y:tmp2.y[0:nn-1]}
			store_data,'mvn_sta_CF_DATA_EN',data={x:tt+en*dt,y:total(t.y[0:nn-1,*],2)}

		ind = where(en eq 0,ndis)							; ndis is number of complete distributions
		if ind[ndis-1]+np-1 gt nn-1 then ndis=ndis-1					; eliminate last distribution if not complete

		if ndis ge 1 then begin

			dat = t.y[ind[0]:ind[0]+ndis*np-1,*]
			tdis = tt[ind[0:ndis-1]]

			get_data,'mvn_STA_CF_MODE',data=md
				md2 = md.y[ind[0:ndis-1]] mod 128
			get_data,'mvn_STA_CF_AVG',data=cavg
			avg = 2^(cavg.y[ind[0:ndis-1]] and 7)
			sum = (cavg.y[ind[0:ndis-1]] and 8)/8

			get_data,'mvn_STA_CF_ATTEN',data=catt						; attenuator state
			att0 = (catt.y[ind[0:ndis-1]] and replicate(192,ndis))/64
			if ndis ge 2 then ind0 = where(att0[0:ndis-2] ne att0[1:ndis-1],cnt0) else cnt0=0
			if cnt0 ne 0 then begin
				for i=0,cnt0-1 do begin
					inds = ind0[i]
					idelay = fix(2/avg[inds])
					if idelay gt 0 then att0[inds+1:(inds+idelay)<(nn-1)] = att0[inds]
				endfor
			endif
			store_data,'mvn_sta_CF_att',data={x:tdis,y:att0}
			ylim,'mvn_sta_CF_att',-1,4,0


			energy=nrg[md.y[ind[0:ndis-1]] mod 128,*]
			energy=total(reform(energy,ndis,4,16),2)/4.
			tmp=decomp19[transpose(reform(dat,np,ndis,16,16,4),[1,0,4,2,3])]	; [time,en,def,an,ma]

			store_data,'mvn_sta_CF_P4B_E',data={x:tdis,y:total(total(total(tmp,5),4),3),v:energy}
			store_data,'mvn_sta_CF_P4B_D',data={x:tdis,y:total(total(total(tmp,5),4),2),v:indgen(4)}
			store_data,'mvn_sta_CF_P4B_A',data={x:tdis,y:total(total(total(tmp,5),3),2),v:indgen(16)}
			store_data,'mvn_sta_CF_P4B_M',data={x:tdis,y:total(total(total(tmp,4),3),2),v:indgen(16)}
			store_data,'mvn_sta_CF_P4B_tot',data={x:tdis,y:total(total(total(total(tmp,5),4),3),2)}
			store_data,'mvn_sta_CF_P4B_all',data={x:tdis,y:reform(tmp,ndis,npts),v:indgen(npts)}

			ylim,'mvn_sta_CF_P4B_tot',0,0,1
			ylim,'mvn_sta_CF_P4B_E',.4,40000.,1
			ylim,'mvn_sta_CF_P4B_D',-1,4,0
			ylim,'mvn_sta_CF_P4B_A',-1,16,0
			ylim,'mvn_sta_CF_P4B_M',-1,16,0
			zlim,'mvn_sta_CF_P4B_E',1,1.e6,1
			zlim,'mvn_sta_CF_P4B_D',1,1.e6,1
			zlim,'mvn_sta_CF_P4B_A',1,1.e6,1
			zlim,'mvn_sta_CF_P4B_M',1,1.e6,1
			options,'mvn_sta_CF_P4B_E',datagap=64.
			options,'mvn_sta_CF_P4B_D',datagap=64.
			options,'mvn_sta_CF_P4B_A',datagap=64.
			options,'mvn_sta_CF_P4B_M',datagap=64.
			options,'mvn_sta_CF_P4B_tot',datagap=64.
			options,'mvn_sta_CF_P4B_E','spec',1
			options,'mvn_sta_CF_P4B_D','spec',1
			options,'mvn_sta_CF_P4B_A','spec',1
			options,'mvn_sta_CF_P4B_M','spec',1
			options,'mvn_sta_CF_P4B_E',ytitle='sta!CP4B-CF!C!CEnergy!C(eV)'
			options,'mvn_sta_CF_P4B_D',ytitle='sta!CP4B-CF!C!CDef!Cbin'
			options,'mvn_sta_CF_P4B_A',ytitle='sta!CP4B-CF!C!CAnode!Cbin'
			options,'mvn_sta_CF_P4B_M',ytitle='sta!CP4B-CF!C!CMass!Cbin'
		endif
	endif
endif


;***************************************************************************************************************

if keyword_set(apids) then test = fix((total(apids eq 'd0') + total(apids eq 'D0')) < 1) else test=0
if not keyword_set(apids) or test then begin

	get_data,'mvn_STA_D0_DATA',data=t
	if size(/type,t) eq 8 then begin
		npts=16384									; 32Ex4Dx16Ax8M
		np = 16										; np is number of packets per measurement, npts/1024
		ind1 = where(t.x gt 0, nn)
		tt=t.x[0:nn-1]

		get_data,'mvn_STA_D0_DIAG',data=diag
		en = (diag.y[0:nn-1] AND 15)
		dt=0.1
			store_data,'mvn_sta_D0_DIAG_EN',data={x:tt+en*dt,y:en}
			get_data,'mvn_STA_D0_SEQ_CNTR',data=tmp2
			store_data,'mvn_sta_D0_SEQ_CNTR_EN',data={x:tt+en*dt,y:tmp2.y[0:nn-1]}
			store_data,'mvn_sta_D0_DATA_EN',data={x:tt+en*dt,y:total(t.y[0:nn-1,*],2)}

		ind = where(en eq 0,ndis)							; ndis is number of complete distributions
		if ind[ndis-1]+np-1 gt nn-1 then ndis=ndis-1					; eliminate last distribution if not complete

		if ndis ge 1 then begin

			dat = t.y[ind[0]:ind[0]+ndis*np-1,*]
			tdis = tt[ind[0:ndis-1]]

			get_data,'mvn_STA_D0_MODE',data=md
				md2 = md.y[ind[0:ndis-1]] mod 128
			get_data,'mvn_STA_D0_AVG',data=cavg
			avg = 2^(cavg.y[ind[0:ndis-1]] and 7)
			sum = (cavg.y[ind[0:ndis-1]] and 8)/8

			get_data,'mvn_STA_D0_ATTEN',data=catt						; attenuator state
			att0 = (catt.y[ind[0:ndis-1]] and replicate(192,ndis))/64
			if ndis ge 2 then ind0 = where(att0[0:ndis-2] ne att0[1:ndis-1],cnt0) else cnt0=0
			if cnt0 ne 0 then begin
				for i=0,cnt0-1 do begin
					inds = ind0[i]
					idelay = fix(2/avg[inds])
					if idelay gt 0 then att0[inds+1:(inds+idelay)<(nn-1)] = att0[inds]
				endfor
			endif
			store_data,'mvn_sta_D0_att',data={x:tdis,y:att0}
			ylim,'mvn_sta_D0_att',-1,4,0


			energy=nrg[md.y[ind[0:ndis-1]] mod 128,*]
			energy=total(reform(energy,ndis,2,32),2)/2.
			tmp=decomp19[reform(transpose(reform(dat,np,ndis,16,8,4,2),[1,5,0,4,2,3]),ndis,32,4,16,8)]	; [time,en,def,an,ma]

			store_data,'mvn_sta_D0_P4C_E',data={x:tdis,y:total(total(total(tmp,5),4),3),v:energy}
			store_data,'mvn_sta_D0_P4C_D',data={x:tdis,y:total(total(total(tmp,5),4),2),v:indgen(4)}
			store_data,'mvn_sta_D0_P4C_A',data={x:tdis,y:total(total(total(tmp,5),3),2),v:indgen(16)}
			store_data,'mvn_sta_D0_P4C_M',data={x:tdis,y:total(total(total(tmp,4),3),2),v:indgen(16)}
			store_data,'mvn_sta_D0_P4C_tot',data={x:tdis,y:total(total(total(total(tmp,5),4),3),2)}
			store_data,'mvn_sta_D0_P4C_all',data={x:tdis,y:reform(tmp,ndis,npts),v:indgen(npts)}

			ylim,'mvn_sta_D0_P4C_tot',0,0,1
			ylim,'mvn_sta_D0_P4C_E',.4,40000.,1
			ylim,'mvn_sta_D0_P4C_D',-1,4,0
			ylim,'mvn_sta_D0_P4C_A',-1,16,0
			ylim,'mvn_sta_D0_P4C_M',-1,8,0
			zlim,'mvn_sta_D0_P4C_E',1,1.e6,1
			zlim,'mvn_sta_D0_P4C_D',1,1.e6,1
			zlim,'mvn_sta_D0_P4C_A',1,1.e6,1
			zlim,'mvn_sta_D0_P4C_M',1,1.e6,1
			options,'mvn_sta_D0_P4C_E',datagap=512.
			options,'mvn_sta_D0_P4C_D',datagap=512.
			options,'mvn_sta_D0_P4C_A',datagap=512.
			options,'mvn_sta_D0_P4C_M',datagap=512.
			options,'mvn_sta_D0_P4C_tot',datagap=512.
			options,'mvn_sta_D0_P4C_E','spec',1
			options,'mvn_sta_D0_P4C_D','spec',1
			options,'mvn_sta_D0_P4C_A','spec',1
			options,'mvn_sta_D0_P4C_M','spec',1
			options,'mvn_sta_D0_P4C_E',ytitle='sta!CP4C-D0!C!CEnergy!C(eV)'
			options,'mvn_sta_D0_P4C_D',ytitle='sta!CP4C-D0!C!CDef!Cbin'
			options,'mvn_sta_D0_P4C_A',ytitle='sta!CP4C-D0!C!CAnode!Cbin'
			options,'mvn_sta_D0_P4C_M',ytitle='sta!CP4C-D0!C!CMass!Cbin'
		endif
	endif
endif


;***************************************************************************************************************

if keyword_set(apids) then test = fix((total(apids eq 'd1') + total(apids eq 'D1')) < 1) else test=0
if not keyword_set(apids) or test then begin

	get_data,'mvn_STA_D1_DATA',data=t
	if size(/type,t) eq 8 then begin
		npts=16384									; 32Ex4Dx16Ax8M
		np = 16										; np is number of packets per measurement, npts/1024
		ind1 = where(t.x gt 0, nn)
		tt=t.x[0:nn-1]

		get_data,'mvn_STA_D1_DIAG',data=diag
		en = (diag.y[0:nn-1] AND 15)
		dt=0.1
			store_data,'mvn_sta_D1_DIAG_EN',data={x:tt+en*dt,y:en}
			get_data,'mvn_STA_D1_SEQ_CNTR',data=tmp2
			store_data,'mvn_sta_D1_SEQ_CNTR_EN',data={x:tt+en*dt,y:tmp2.y[0:nn-1]}
			store_data,'mvn_sta_D1_DATA_EN',data={x:tt+en*dt,y:total(t.y[0:nn-1,*],2)}

		ind = where(en eq 0,ndis)							; ndis is number of complete distributions
		if ind[ndis-1]+np-1 gt nn-1 then ndis=ndis-1					; eliminate last distribution if not complete

		if ndis ge 1 then begin

			dat = t.y[ind[0]:ind[0]+ndis*np-1,*]
			tdis = tt[ind[0:ndis-1]]

			get_data,'mvn_STA_D1_MODE',data=md
				md2 = md.y[ind[0:ndis-1]] mod 128
			get_data,'mvn_STA_D1_AVG',data=cavg
			avg = 2^(cavg.y[ind[0:ndis-1]] and 7)
			sum = (cavg.y[ind[0:ndis-1]] and 8)/8

			get_data,'mvn_STA_D1_ATTEN',data=catt						; attenuator state
			att0 = (catt.y[ind[0:ndis-1]] and replicate(192,ndis))/64
			if ndis ge 2 then ind0 = where(att0[0:ndis-2] ne att0[1:ndis-1],cnt0) else cnt0=0
			if cnt0 ne 0 then begin
				for i=0,cnt0-1 do begin
					inds = ind0[i]
					idelay = fix(2/avg[inds])
					if idelay gt 0 then att0[inds+1:(inds+idelay)<(nn-1)] = att0[inds]
				endfor
			endif
			store_data,'mvn_sta_D1_att',data={x:tdis,y:att0}
			ylim,'mvn_sta_D1_att',-1,4,0


			energy=nrg[md.y[ind[0:ndis-1]] mod 128,*]
			energy=total(reform(energy,ndis,2,32),2)/2.
			tmp=decomp19[reform(transpose(reform(dat,np,ndis,16,8,4,2),[1,5,0,4,2,3]),ndis,32,4,16,8)]	; [time,en,def,an,ma]

			store_data,'mvn_sta_D1_P4C_E',data={x:tdis,y:total(total(total(tmp,5),4),3),v:energy}
			store_data,'mvn_sta_D1_P4C_D',data={x:tdis,y:total(total(total(tmp,5),4),2),v:indgen(4)}
			store_data,'mvn_sta_D1_P4C_A',data={x:tdis,y:total(total(total(tmp,5),3),2),v:indgen(16)}
			store_data,'mvn_sta_D1_P4C_M',data={x:tdis,y:total(total(total(tmp,4),3),2),v:indgen(16)}
			store_data,'mvn_sta_D1_P4C_tot',data={x:tdis,y:total(total(total(total(tmp,5),4),3),2)}
			store_data,'mvn_sta_D1_P4C_all',data={x:tdis,y:reform(tmp,ndis,npts),v:indgen(npts)}

			ylim,'mvn_sta_D1_P4C_tot',0,0,1
			ylim,'mvn_sta_D1_P4C_E',.4,40000.,1
			ylim,'mvn_sta_D1_P4C_D',-1,4,0
			ylim,'mvn_sta_D1_P4C_A',-1,16,0
			ylim,'mvn_sta_D1_P4C_M',-1,8,0
			zlim,'mvn_sta_D1_P4C_E',1,1.e6,1
			zlim,'mvn_sta_D1_P4C_D',1,1.e6,1
			zlim,'mvn_sta_D1_P4C_A',1,1.e6,1
			zlim,'mvn_sta_D1_P4C_M',1,1.e6,1
			options,'mvn_sta_D1_P4C_E',datagap=64.
			options,'mvn_sta_D1_P4C_D',datagap=64.
			options,'mvn_sta_D1_P4C_A',datagap=64.
			options,'mvn_sta_D1_P4C_M',datagap=64.
			options,'mvn_sta_D1_P4C_tot',datagap=64.
			options,'mvn_sta_D1_P4C_E','spec',1
			options,'mvn_sta_D1_P4C_D','spec',1
			options,'mvn_sta_D1_P4C_A','spec',1
			options,'mvn_sta_D1_P4C_M','spec',1
			options,'mvn_sta_D1_P4C_E',ytitle='sta!CP4C-D1!C!CEnergy!C(eV)'
			options,'mvn_sta_D1_P4C_D',ytitle='sta!CP4C-D1!C!CDef!Cbin'
			options,'mvn_sta_D1_P4C_A',ytitle='sta!CP4C-D1!C!CAnode!Cbin'
			options,'mvn_sta_D1_P4C_M',ytitle='sta!CP4C-D1!C!CMass!Cbin'
		endif
	endif
endif


;***************************************************************************************************************

if keyword_set(apids) then test = fix((total(apids eq 'd2') + total(apids eq 'D2')) < 1) else test=0
if not keyword_set(apids) or test then begin

	get_data,'mvn_STA_D2_DATA',data=t
	if size(/type,t) eq 8 then begin
		npts=4096									; 32Ex1Dx16Ax8M
		np = 4										; np is number of packets per measurement, npts/1024
		ind1 = where(t.x gt 0, nn)
		tt=t.x[0:nn-1]

		get_data,'mvn_STA_D2_DIAG',data=diag
		en = (diag.y[0:nn-1] AND 15)
		dt=0.1
			store_data,'mvn_sta_D2_DIAG_EN',data={x:tt+en*dt,y:en}
			get_data,'mvn_STA_D2_SEQ_CNTR',data=tmp2
			store_data,'mvn_sta_D2_SEQ_CNTR_EN',data={x:tt+en*dt,y:tmp2.y[0:nn-1]}
			store_data,'mvn_sta_D2_DATA_EN',data={x:tt+en*dt,y:total(t.y[0:nn-1,*],2)}

		ind = where(en eq 0,ndis)							; ndis is number of complete distributions
		if ind[ndis-1]+np-1 gt nn-1 then ndis=ndis-1					; eliminate last distribution if not complete

		if ndis ge 1 then begin

			dat = t.y[ind[0]:ind[0]+ndis*np-1,*]
			tdis = tt[ind[0:ndis-1]]

			get_data,'mvn_STA_D2_MODE',data=md
				md2 = md.y[ind[0:ndis-1]] mod 128
			get_data,'mvn_STA_D2_AVG',data=cavg
			avg = 2^(cavg.y[ind[0:ndis-1]] and 7)
			sum = (cavg.y[ind[0:ndis-1]] and 8)/8

			get_data,'mvn_STA_D2_ATTEN',data=catt						; attenuator state
			att0 = (catt.y[ind[0:ndis-1]] and replicate(192,ndis))/64
			if ndis ge 2 then ind0 = where(att0[0:ndis-2] ne att0[1:ndis-1],cnt0) else cnt0=0
			if cnt0 ne 0 then begin
				for i=0,cnt0-1 do begin
					inds = ind0[i]
					idelay = fix(2/avg[inds])
					if idelay gt 0 then att0[inds+1:(inds+idelay)<(nn-1)] = att0[inds]
				endfor
			endif
			store_data,'mvn_sta_D2_att',data={x:tdis,y:att0}
			ylim,'mvn_sta_D2_att',-1,4,0


			energy=nrg[md.y[ind[0:ndis-1]] mod 128,*]
			energy=total(reform(energy,ndis,2,32),2)/2.
			tmp=decomp19[reform(transpose(reform(dat,np,ndis,16,8,8),[1,4,0,2,3]),ndis,32,16,8)]	; [time,en,an,ma]

			store_data,'mvn_sta_D2_P4D_E',data={x:tdis,y:total(total(tmp,4),3),v:energy}
			store_data,'mvn_sta_D2_P4D_A',data={x:tdis,y:total(total(tmp,4),2),v:indgen(16)}
			store_data,'mvn_sta_D2_P4D_M',data={x:tdis,y:total(total(tmp,3),2),v:indgen(16)}
			store_data,'mvn_sta_D2_P4D_tot',data={x:tdis,y:total(total(total(tmp,4),3),2)}
			store_data,'mvn_sta_D2_P4D_all',data={x:tdis,y:reform(tmp,ndis,npts),v:indgen(npts)}

			ylim,'mvn_sta_D2_P4D_tot',0,0,1
			ylim,'mvn_sta_D2_P4D_E',.4,40000.,1
			ylim,'mvn_sta_D2_P4D_A',-1,16,0
			ylim,'mvn_sta_D2_P4D_M',-1,8,0
			zlim,'mvn_sta_D2_P4D_E',1,1.e6,1
			zlim,'mvn_sta_D2_P4D_A',1,1.e6,1
			zlim,'mvn_sta_D2_P4D_M',1,1.e6,1
			options,'mvn_sta_D2_P4D_E',datagap=256.
			options,'mvn_sta_D2_P4D_A',datagap=256.
			options,'mvn_sta_D2_P4D_M',datagap=256.
			options,'mvn_sta_D2_P4D_tot',datagap=256.
			options,'mvn_sta_D2_P4D_E','spec',1
			options,'mvn_sta_D2_P4D_A','spec',1
			options,'mvn_sta_D2_P4D_M','spec',1
			options,'mvn_sta_D2_P4D_E',ytitle='sta!CP4D-D2!C!CEnergy!C(eV)'
			options,'mvn_sta_D2_P4D_A',ytitle='sta!CP4D-D2!C!CAnode!Cbin'
			options,'mvn_sta_D2_P4D_M',ytitle='sta!CP4D-D2!C!CMass!Cbin'
		endif
	endif
endif


;***************************************************************************************************************

if keyword_set(apids) then test = fix((total(apids eq 'd3') + total(apids eq 'D3')) < 1) else test=0
if not keyword_set(apids) or test then begin

	get_data,'mvn_STA_D3_DATA',data=t
	if size(/type,t) eq 8 then begin
		npts=4096									; 32Ex1Dx16Ax8M
		np = 4										; np is number of packets per measurement, npts/1024
		ind1 = where(t.x gt 0, nn)
		tt=t.x[0:nn-1]

		get_data,'mvn_STA_D3_DIAG',data=diag
		en = (diag.y[0:nn-1] AND 15)
		dt=0.1
			store_data,'mvn_sta_D3_DIAG_EN',data={x:tt+en*dt,y:en}
			get_data,'mvn_STA_D3_SEQ_CNTR',data=tmp2
			store_data,'mvn_sta_D3_SEQ_CNTR_EN',data={x:tt+en*dt,y:tmp2.y[0:nn-1]}
			store_data,'mvn_sta_D3_DATA_EN',data={x:tt+en*dt,y:total(t.y[0:nn-1,*],2)}

		ind = where(en eq 0,ndis)							; ndis is number of complete distributions
		if ind[ndis-1]+np-1 gt nn-1 then ndis=ndis-1					; eliminate last distribution if not complete

		if ndis ge 1 then begin

			dat = t.y[ind[0]:ind[0]+ndis*np-1,*]
			tdis = tt[ind[0:ndis-1]]

			get_data,'mvn_STA_D3_MODE',data=md
				md2 = md.y[ind[0:ndis-1]] mod 128
			get_data,'mvn_STA_D3_AVG',data=cavg
			avg = 2^(cavg.y[ind[0:ndis-1]] and 7)
			sum = (cavg.y[ind[0:ndis-1]] and 8)/8

			get_data,'mvn_STA_D3_ATTEN',data=catt						; attenuator state
			att0 = (catt.y[ind[0:ndis-1]] and replicate(192,ndis))/64
			if ndis ge 2 then ind0 = where(att0[0:ndis-2] ne att0[1:ndis-1],cnt0) else cnt0=0
			if cnt0 ne 0 then begin
				for i=0,cnt0-1 do begin
					inds = ind0[i]
					idelay = fix(2/avg[inds])
					if idelay gt 0 then att0[inds+1:(inds+idelay)<(nn-1)] = att0[inds]
				endfor
			endif
			store_data,'mvn_sta_D3_att',data={x:tdis,y:att0}
			ylim,'mvn_sta_D3_att',-1,4,0


			energy=nrg[md.y[ind[0:ndis-1]] mod 128,*]
			energy=total(reform(energy,ndis,2,32),2)/2.
			tmp=decomp19[reform(transpose(reform(dat,np,ndis,16,8,8),[1,4,0,2,3]),ndis,32,16,8)]	; [time,en,an,ma]

			store_data,'mvn_sta_D3_P4D_E',data={x:tdis,y:total(total(tmp,4),3),v:energy}
			store_data,'mvn_sta_D3_P4D_A',data={x:tdis,y:total(total(tmp,4),2),v:indgen(16)}
			store_data,'mvn_sta_D3_P4D_M',data={x:tdis,y:total(total(tmp,3),2),v:indgen(16)}
			store_data,'mvn_sta_D3_P4D_tot',data={x:tdis,y:total(total(total(tmp,4),3),2)}
			store_data,'mvn_sta_D3_P4D_all',data={x:tdis,y:reform(tmp,ndis,npts),v:indgen(npts)}

			ylim,'mvn_sta_D3_P4D_tot',0,0,1
			ylim,'mvn_sta_D3_P4D_E',.4,40000.,1
			ylim,'mvn_sta_D3_P4D_A',-1,16,0
			ylim,'mvn_sta_D3_P4D_M',-1,8,0
			zlim,'mvn_sta_D3_P4D_E',1,1.e6,1
			zlim,'mvn_sta_D3_P4D_A',1,1.e6,1
			zlim,'mvn_sta_D3_P4D_M',1,1.e6,1
			options,'mvn_sta_D3_P4D_E',datagap=64.
			options,'mvn_sta_D3_P4D_A',datagap=64.
			options,'mvn_sta_D3_P4D_M',datagap=64.
			options,'mvn_sta_D3_P4D_tot',datagap=64.
			options,'mvn_sta_D3_P4D_E','spec',1
			options,'mvn_sta_D3_P4D_A','spec',1
			options,'mvn_sta_D3_P4D_M','spec',1
			options,'mvn_sta_D3_P4D_E',ytitle='sta!CP4D-D3!C!CEnergy!C(eV)'
			options,'mvn_sta_D3_P4D_A',ytitle='sta!CP4D-D3!C!CAnode!Cbin'
			options,'mvn_sta_D3_P4D_M',ytitle='sta!CP4D-D3!C!CMass!Cbin'
		endif
	endif
endif


	store_data,'mvn_sta_P4_E',data=['mvn_sta_CC_P4A_E','mvn_sta_CE_P4B_E','mvn_sta_D0_P4C_E','mvn_sta_D2_P4D_E']
		ylim,'mvn_sta_P4_E',.4,40000.,1
		options,'mvn_sta_P4_E',ytitle='sta!CP4D !C!C Energy (eV)'
	store_data,'mvn_sta_P4_M',data=['mvn_sta_CC_P4A_M','mvn_sta_CE_P4B_M','mvn_sta_D0_P4C_M','mvn_sta_D2_P4D_M']
		ylim,'mvn_sta_P4_M',-1,32
		options,'mvn_sta_P4_M',ytitle='sta!CP4D !C!C Mass bin'
	store_data,'mvn_sta_P4_D',data=['mvn_sta_CC_P4A_D','mvn_sta_CE_P4B_D','mvn_sta_D0_P4C_D']
		ylim,'mvn_sta_P4_D',-1,8
		options,'mvn_sta_P4_D',ytitle='sta!CP4D !C!C Def bin'
	store_data,'mvn_sta_P4_A',data=['mvn_sta_CE_P4B_A','mvn_sta_D0_P4C_A','mvn_sta_D2_P4D_A']
		ylim,'mvn_sta_P4_A',-1,16
		options,'mvn_sta_P4_A',ytitle='sta!CP4D !C!C Anode'

		
;***************************************************************************************************************

if keyword_set(apids) then test = fix((total(apids eq 'd4') + total(apids eq 'D4')) < 1) else test=0
if not keyword_set(apids) or test then begin

	get_data,'mvn_STA_D4_DATA',data=t
	if size(/type,t) eq 8 then begin
		npts=128									; 1Ex4Dx16Ax2M
		ind1 = where(t.x gt 0, nn)
		tt=t.x[0:nn-1]

		get_data,'mvn_STA_D4_MODE',data=md
			md2 = md.y[0:nn-1] mod 128
		get_data,'mvn_STA_D4_AVG',data=cavg
		avg = cavg.y and replicate(7,nn)
		sum = (cavg.y and replicate(8,nn))/8
		comp = (cavg.y and replicate(192,nn))/128
		ind = where(avg[0:nn-2] ne avg[1:nn-1],count)
		if count ne 0 then begin
			for i=0,count-1 do begin
				a0 = 2^avg[ind[i]]
				a1 = 2^avg[ind[i]+1]
				dt = a1*4.
				nf = abs(round( ( (tt[ind[i]+1]-tt[ind[i]])/4. -(a0+a1)/2. ) /(a1-a0) ))
				if nf gt 7 then nf=0
				for j=ind[i]-nf+1,ind[i] do tt[j]=tt[ind[i]+1]-dt*(ind[i]+1-j)
			endfor
		endif
		store_data,'mvn_sta_D4_mode',data={x:tt,y:md2}
			ylim,'mvn_sta_D4_mode',-1,6,0

		tmp=decomp19[transpose(reform(t.y[0:nn-1,*],nn,16,2,4),[0,3,1,2])]			; [time,def,an,ma]

		store_data,'mvn_sta_D4_P4E_D',data={x:tt,y:reform(total(total(tmp,4),3),nn,4),v:findgen(4)}
		store_data,'mvn_sta_D4_P4E_A',data={x:tt,y:reform(total(total(tmp,4),2),nn,16),v:findgen(16)}
		store_data,'mvn_sta_D4_P4E_M',data={x:tt,y:reform(total(total(tmp,2),2),nn,2),v:findgen(2)}
		store_data,'mvn_sta_D4_P4E_tot',data={x:tt,y:total(total(total(tmp,2),2),2)}
		store_data,'mvn_sta_D4_P4E_all',data={x:tt,y:reform(tmp,nn,npts),v:findgen(npts)}

		ylim,'mvn_sta_D4_P4E_tot',0,0,1
		ylim,'mvn_sta_D4_P4E_D',-1,4,0
		ylim,'mvn_sta_D4_P4E_A',-1,16,0
		ylim,'mvn_sta_D4_P4E_M',-1,2,0
		zlim,'mvn_sta_D4_P4E_D',1,1.e6,1
		zlim,'mvn_sta_D4_P4E_A',1,1.e6,1
		zlim,'mvn_sta_D4_P4E_M',1,1.e6,1
		options,'mvn_sta_D4_P4E_D',datagap=64.
		options,'mvn_sta_D4_P4E_A',datagap=64.
		options,'mvn_sta_D4_P4E_M',datagap=64.
		options,'mvn_sta_D4_P4E_tot',datagap=64.
		options,'mvn_sta_D4_P4E_D','spec',1
		options,'mvn_sta_D4_P4E_A','spec',1
		options,'mvn_sta_D4_P4E_M','spec',1
		options,'mvn_sta_D4_P4E_D',ytitle='sta!CP4E-D4!C!CDef!Cbin'
		options,'mvn_sta_D4_P4E_A',ytitle='sta!CP4E-D4!C!CAnode!Cbin'
		options,'mvn_sta_D4_P4E_M',ytitle='sta!CP4E-D4!C!CMass!Cbin'
	endif
endif


;***************************************************************************************************************
; this section needs a fix

if keyword_set(apids) then test = fix((total(apids eq 'd3') + total(apids eq 'D3')) < 1) else test=0

if 0 and not keyword_set(apids) or test then begin

	get_data,'mvn_STA_D6_DATA',data=t
	if size(/type,t) eq 8 then begin
		npts=49152									; 32Ex8Dx1Ax32M
		np = 48										; np is number of packets per measurement
		ind1 = where(t.x gt 0, nn)
		tt=t.x[0:nn-1]

		get_data,'mvn_STA_D6_MODE',data=md
		get_data,'mvn_STA_D6_AVG',data=avg
		get_data,'mvn_STA_D6_DIAG',data=diag
		en = (diag.y[0:nn-1] AND 63)
		dt=0.1
			store_data,'mvn_sta_D6_DIAG_EN',data={x:tt+en*dt,y:en}
			get_data,'mvn_STA_D6_SEQ_CNTR',data=tmp2
			store_data,'mvn_sta_D6_SEQ_CNTR_EN',data={x:tt+en*dt,y:tmp2.y[0:nn-1]}
;			store_data,'mvn_sta_D6_DATA_EN',data={x:tt+en*dt,y:total(t.y[0:nn-1,*],2)}

		ind = where(en eq 0,ndis)							; ndis is number of complete distributions
		if ind[ndis-1]+np-1 gt nn-1 then ndis=ndis-1					; eliminate last distribution if not complete

		if ndis ge 1 then begin
			events=intarr(340*nn,10)

			for i=0,ndis-1 do begin
			  for j=0,47 do begin
				while ind lt lst do begin

;*************************************************************
;	pntr=0l
;	pktl=0l
;	npkt=0l
;	ntot=0l
;	hdr  = adat(pntr:pntr+hdrl-1)
;	last_apid_cntr= hdr(3)
;	last_time = (1.d*((hdr[6]*256ul+hdr[7])*256ul+hdr[8])*256ul+hdr[9])+(hdr[10]+hdr[11]/256.d)/256.
;while pntr+1 lt data_size do begin
;	hdr  = adat[pntr:pntr+hdrl-1]
;	apid_cntr= hdr(3)
;	spin_num = hdr[15]
;	time = (1.d*((hdr[6]*256ul+hdr[7])*256ul+hdr[8])*256ul+hdr[9])+(hdr[10]+hdr[11]/256.d)/256.
;	time2 = time
;	pktl = hdr[4]*256+hdr[5]+7-hdrl
;	pktp = 0l & datl=0
;	pkt  = adat[pntr+hdrl:pntr+hdrl+pktl-1]
;	sub_spin=0
;	last_mode=mode_decode[pkt[0],pkt[1]]
;	last_apid_cntr= byte(apid_cntr+1)
;	mode=last_mode 
;	if mode ne -1 then datl=dat_len[mode]
;
;
;	while pktp+datl+4 le pktl do begin
;
;		config_1 = pkt[pktp] & config_2 = pkt[pktp+1]
;		mode = mode_decode[config_1,config_2]
;		if mode ne -1 and mode eq last_mode then begin
;
;			datl = dat_len[mode]
;		endif else begin
;			pktp = pktl
;			dprint,dlevel=1,'Error: Invalid ESA config header bytes or mode change in packet - skipping packet'
;			dprint,dlevel=1,mode,config_1,config_2,pntr,pktl,pktp
;			datl=0 & nsp=1
;		endelse
;
;		pktp = pktp + 4 + datl
;
;		sub_spin = sub_spin + nsp
;		last_mode = mode	
;		if ntot-1000l*ceil(ntot/1000l) eq 0 then dprint,dlevel=1,pntr,'  ',pktl,'  ',time_string(time),'  ',adat[pntr+hdrl],'  ',adat[pntr+hdrl+1]
;	endwhile
;	pntr = pntr + pktl + hdrl
;	dprint,dlevel=1,pntr,'  ',ntot
;	last_time=time
;
;endwhile
;
;*******************************************************************************
			  	endwhile
			  endfor	
			endfor

			dat = t.y[ind[0]:ind[0]+ndis*np-1,*]
			tdis = tt[ind[0:ndis-1]]

			energy=nrg[md.y[ind[0:ndis-1]] mod 128,*]
			energy=total(reform(energy,ndis,2,32),2)/2.
			tmp=decomp19[reform(transpose(reform(dat,np,ndis,32,8,4),[1,4,0,3,2]),ndis,32,8,32)]		; [time,en,def,ma]

			store_data,'mvn_sta_D6_P4A_E',data={x:tdis,y:total(total(tmp,4),3),v:energy}
			store_data,'mvn_sta_D6_P4A_D',data={x:tdis,y:total(total(tmp,4),2),v:indgen(8)}
			store_data,'mvn_sta_D6_P4A_M',data={x:tdis,y:total(total(tmp,2),2),v:indgen(32)}
			store_data,'mvn_sta_D6_P4A_tot',data={x:tdis,y:total(total(total(tmp,4),3),2)}
			store_data,'mvn_sta_D6_P4A_all',data={x:tdis,y:reform(tmp,ndis,npts),v:indgen(npts)}

			ylim,'mvn_sta_D6_P4A_tot',0,0,1
			ylim,'mvn_sta_D6_P4A_E',.4,40000.,1
			ylim,'mvn_sta_D6_P4A_D',-1,8,0
			ylim,'mvn_sta_D6_P4A_M',-1,32,0
			zlim,'mvn_sta_D6_P4A_E',1,1.e6,1
			zlim,'mvn_sta_D6_P4A_D',1,1.e6,1
			zlim,'mvn_sta_D6_P4A_M',1,1.e6,1
			options,'mvn_sta_D6_P4A_E',datagap=128.
			options,'mvn_sta_D6_P4A_D',datagap=128.
			options,'mvn_sta_D6_P4A_M',datagap=128.
			options,'mvn_sta_D6_P4A_tot',datagap=128.
			options,'mvn_sta_D6_P4A_E','spec',1
			options,'mvn_sta_D6_P4A_D','spec',1
			options,'mvn_sta_D6_P4A_M','spec',1
			options,'mvn_sta_D6_P4A_E',ytitle='sta!CP4A-D6!C!CEnergy!C(eV)'
			options,'mvn_sta_D6_P4A_D',ytitle='sta!CP4A-D6!C!CDef!Cbin'
			options,'mvn_sta_D6_P4A_M',ytitle='sta!CP4A-D6!C!CMass!Cbin'
		endif
	endif
endif


;***************************************************************************************************************

if keyword_set(apids) then test = fix((total(apids eq 'd7') + total(apids eq 'D7')) < 1) else test=0
if not keyword_set(apids) or test then begin

	get_data,'mvn_STA_D7_DATA',data=t
	if size(/type,t) eq 8 then begin
		npts=512									; 512 points
		dt = 4.d/512.*findgen(512)
		ind1 = where(t.x gt 0, nn)
		tt= reform(replicate(1.d,512)#t.x[0:nn-1] + dt#replicate(1.d,nn),nn*512)
		dd = reform(transpose(t.y[0:nn-1,*]),nn*512)
		fshkpcal = -5./2^15
		store_data,'mvn_sta_D7_data',data={x:tt,y:dd}
		store_data,'mvn_sta_D7_data_cal',data={x:tt,y:dd*fshkpcal}
	endif
endif


;***************************************************************************************************************

if keyword_set(apids) then test = fix((total(apids eq 'd8') + total(apids eq 'D8')) < 1) else test=0
if not keyword_set(apids) or test then begin

	get_data,'mvn_STA_D8_DATA',data=t
	if size(/type,t) eq 8 then begin
		npts=12									; 12R
		ind1 = where(t.x gt 0, nn)
		tt=t.x[0:nn-1]

		get_data,'mvn_STA_D8_MODE',data=md
		get_data,'mvn_STA_D8_AVG',data=cavg
		avg = cavg.y and replicate(7,nn)
		sum = (cavg.y and replicate(8,nn))/8
		comp = (cavg.y and replicate(192,nn))/128
		ind = where(avg[0:nn-2] ne avg[1:nn-1],count)
		if count ne 0 then begin
			for i=0,count-1 do begin
				a0 = 2^avg[ind[i]]
				a1 = 2^avg[ind[i]+1]
				dt = a1*4.
				nf = abs(round( ( (tt[ind[i]+1]-tt[ind[i]])/4. -(a0+a1)/2. ) /(a1-a0) ))
				if nf gt 15 then nf=0
				for j=ind[i]-nf+1,ind[i] do tt[j]=tt[ind[i]+1]-dt*(ind[i]+1-j)
			endfor
		endif
		store_data,'mvn_sta_D8_mode',data={x:tt,y:md.y[0:nn-1]}
			ylim,'mvn_sta_D8_mode',-1,6,0

		tmp=decomp19[t.y[0:nn-1,*]]						; [time,rate]

		store_data,'mvn_sta_D8_R1_Time_ABCD',data={x:tt,y:tmp[*,0:3],v:findgen(4)}
		store_data,'mvn_sta_D8_R1_Time_RST',	data={x:tt,y:tmp[*,4]}
		store_data,'mvn_sta_D8_R1_Time_NoStart',	data={x:tt,y:tmp[*,5]}
		store_data,'mvn_sta_D8_R1_Time_Unqual',	data={x:tt,y:tmp[*,6]}
		store_data,'mvn_sta_D8_R1_Time_Qual',	data={x:tt,y:tmp[*,7]}
		store_data,'mvn_sta_D8_R1_Time_AnRej',	data={x:tt,y:tmp[*,8]}
		store_data,'mvn_sta_D8_R1_Time_MaRej',	data={x:tt,y:tmp[*,9]}
		store_data,'mvn_sta_D8_R1_Time_A&B',	data={x:tt,y:tmp[*,10]}
		store_data,'mvn_sta_D8_R1_Time_C&D',	data={x:tt,y:tmp[*,11]}
		store_data,'mvn_sta_D8_R1_all',data={x:tt,y:tmp,v:findgen(12)}


		options,'mvn_sta_D8_R1_Time_ABCD','colors',[cols.green,cols.blue,cols.red,cols.yellow]

		ylim,'mvn_sta_D8_R1*',0,0,1
		options,'mvn_sta_D8_R1*',datagap=64.
	endif
endif

;***************************************************************************************************************

if keyword_set(apids) then test = fix((total(apids eq 'd9') + total(apids eq 'D9')) < 1) else test=0
if not keyword_set(apids) or test then begin

	get_data,'mvn_STA_D9_DATA',data=t
	if size(/type,t) eq 8 then begin
		npts=768									; 32Ex1Dx1Ax32M
		ind1 = where(t.x gt 0, nn)
		tt=t.x[0:nn-1]

		get_data,'mvn_STA_D9_MODE',data=md

		tmp=decomp19[transpose(reform(t.y[0:nn-1,*],nn,12,64),[0,2,1])]			; [time,en,ma]

		store_data,'mvn_sta_D9_R2_Eeff',data={x:tt,y:reform(tmp[*,*,7]/(tmp[*,*,4]+.1),nn,64),v:findgen(64)}
		store_data,'mvn_sta_D9_R2_Aeff',data={x:tt,y:reform(total(tmp[*,*,7],2)/(total(tmp[*,*,4],2)+.1),nn)}
		store_data,'mvn_sta_D9_R2_RST',data={x:tt,y:reform(total(tmp[*,*,4],2))}
		store_data,'mvn_sta_D9_R2_all',data={x:tt,y:reform(tmp,nn,npts),v:findgen(npts)}

		ylim,'mvn_sta_D9_R2*',0,0,1
		options,'mvn_sta_D9_R2_Eeff',datagap=64.
		options,'mvn_sta_D9_R2_Aeff',datagap=64.
		options,'mvn_sta_D9_R2_RST',datagap=64.
		options,'mvn_sta_D9_R2_all',datagap=64.
		options,'mvn_sta_D9_R2_Eeff','spec',1
		options,'mvn_sta_D9_R2_Eeff',ytitle='sta!CEeff-D9!C!CEnergy!Cbin'
		zlim,'mvn_sta_D9_R2_Eeff',0.,1.,0
	endif
endif


;***************************************************************************************************************

if keyword_set(apids) then test = fix((total(apids eq 'da') + total(apids eq 'DA')) < 1) else test=0
if not keyword_set(apids) or test then begin

	get_data,'mvn_STA_DA_DATA',data=t
	if size(/type,t) eq 8 then begin
		npts=1024									; 32Ex1Dx1Ax32M

		ind1 = where(t.x gt 0, nn)
		tt=t.x[0:nn-1]
		da=decomp19[t.y[0:nn-1,*]]

;	Correct tt for mode change boundaries, code could be simplified
;	assume first measurement in a packet may be an average at old cadence
;	This may not be needed since for flight I plan to not average data for apid DA
 
		nda = indgen(nn)
		for k=0,nn-1 do begin
			itmp = where(da[k,*] ne decomp19[255],count)
			nda[k]=count
		endfor
		store_data,'mvn_sta_DA_ndata',data={x:tt,y:nda}
		store_data,'mvn_sta_DA_ndata2',data={x:tt,y:1025-nda}
			ylim,'mvn_sta_DA_ndata2',.5,1024,1

		get_data,'mvn_STA_DA_AVG',data=avg
			gg = [64,128,256,1024]
			g2 = (avg.y[0:nn-1] and 48)/2^4		; g2 = ind into gg array
			nd = gg[g2]				; nd = # data points per sweep
			nm = nda/nd				; nm = # sweeps in a packet
			cd = 2^(avg.y[0:nn-1] and 7)		; cd = # sweeps between measurements, cadence
		dt3=tt
	; NOTE: corrections like this can screw up if packets are dropped.
		if nn gt 2 then begin
			dt3[0:nn-2] = (tt[1:nn-1] - tt[0:nn-2] - nm[0:nn-2]*cd[0:nn-2]*4.d)
;			print,dt3
			tt[0:nn-2] = tt[0:nn-2] + dt3[0:nn-2]
		endif


; 	Necessary code

		get_data,'mvn_STA_DA_MODE',data=md
			md2 = md.y[0:nn-1] mod 128
		energy=nrg[md2,*]

;	this section makes time series out of DA data
		; the following is correct when cd=1, spreads our measurements in time otherwise 
		get_data,'mvn_STA_DA_AVG',data=avg
			gg = [64,128,256,1024]
			g2 = (avg.y[0:nn-1] and 48)/2^4		; g2 = ind into gg array
			nd = gg[g2]				; nd = # data points per sweep
			nm = 1024/nd				; nm = # sweeps in a packet
			cd = 2^(avg.y[0:nn-1] and 7)		; cd = # sweeps between measurements, cadence
			ss = (avg.y[0:nn-1] and 8)/2^3		; ss = snapshot or average
			n2 = (ss*(cd-1)+1) 			; n2 = # sweeps averaged per measurement, either "cd" or 1
			dt = n2*4.d/nd				; dt = averaging time if  correct when data are averaged over multiple sweeps 
; old			tt1=tt-2.0d*n2*nm			; old code, the -2.0d*n2*nm corrects back to start time from mvn_sta_apid_decom.pro, removed
			tt1=tt

			tp1 = dblarr(1024,4) & tp1[*,0]=findgen(64)#replicate(1.,16) & tp1[*,1]=findgen(128)#replicate(1.,8) & tp1[*,2]=findgen(256)#replicate(1.,4)  & tp1[*,3]=findgen(1024)
			tp2 = dblarr(1024,4) & tp2[*,0]=replicate(1.,64)#findgen(16) & tp2[*,1]=replicate(1.,128)#findgen(8) & tp2[*,2]=replicate(1.,256)#findgen(4)


;			tt3 = reform((findgen(1024)+.5)#dt + tp2[*,g2]*(replicate(1.,1024)#cd) + replicate(1.,1024)#tt1,1024*nn)
			tt3 = reform(tp1[*,g2]*(replicate(1.,1024)#dt) + tp2[*,g2]*(replicate(1.,1024)#cd*4.d) + replicate(1.,1024)#tt1,1024*nn)

; old			tt3 = reform((findgen(1024)+.5)#dt + replicate(1.,1024)#tt1,1024*nn)
			da3 = reform(transpose(da),nn*1024)
			in3 = where(da3 ne decomp19[255],ct3) 

		store_data,'mvn_sta_DA_R3_T',data={x:tt3[in3],y:da3[in3]}
			ylim,'mvn_sta_DA_R3_T',0,0,1
			options,'mvn_sta_DA_R3_T',datagap=64.
			options,'mvn_sta_DA_R3_T',psym=1


		tm = total(nm)
		da2 = fltarr(tm,64) 					
		tt2 = dblarr(tm)					
		en2 = fltarr(tm,64)

;	     this section makes energy sweeps out of DA data
		m=0
		for i=0,nn-1 do begin									; nn packets
			for j=0,nm[i]-1 do begin							; nm sweeps per packet
				jj=j*nd[i]								; jj index into da[0:nn-1,0:1024] array
				if da[i,jj] ne decomp19[255] then begin
					tmp      = total(reform(da[i,jj:jj+nd[i]-1],nd[i]/64,64),1)
					da2[m,*] = tmp
; old					tt2[m]   = tt[i] - 2.0d*n2[i]*(nm[i]-1) + j*4.d*n2[i]		
					tt2[m]   = tt[i] + 2.0d*n2[i]           + j*4.d*cd[i]		
					en2[m,*] = energy[i,*]
					m=m+1
				endif
			endfor
		endfor
		ind2 = where(tt2 gt 0, n) 
		store_data,'mvn_sta_DA_R3_E',data={x:tt2[0:n-1],y:da2[0:n-1,*],v:en2[0:n-1,*]}
			ylim,'mvn_sta_DA_R3_E',.4,40000,1
			zlim,'mvn_sta_DA_R3_E',1,1.e6,1
			options,'mvn_sta_DA_R3_E',datagap=64.
			options,'mvn_sta_DA_R3_E',spec=1

		store_data,'mvn_sta_DA_R3_sum',data={x:tt2[0:n-1],y:total(da2[0:n-1,*],2)}
			ylim,'mvn_sta_DA_R3_sum',0,0,1
			options,'mvn_sta_DA_R3_sum',datagap=64.

	endif
endif


;***************************************************************************************************************

if keyword_set(apids) then test = fix((total(apids eq 'db') + total(apids eq 'DB')) < 1) else test=0
if not keyword_set(apids) or test then begin

	get_data,'mvn_STA_DB_DATA',data=t
	if size(/type,t) eq 8 then begin
		npts=1024									; 32Ex1Dx1Ax32M
		ind = where(t.x gt 0, nn)
		tt=t.x[0:nn-1]
		store_data,'mvn_sta_DB_DATA',data={x:tt,y:t.y[0:nn-1,*],v:findgen(1024)}

		tmp=decomp19[t.y[0:nn-1,*]]							; [time,ma]

			aa = fltarr(nn,64)
			ind1=indgen(8) & ind2=indgen(16) & ind3=indgen(32)
			for i=0,nn-1 do begin 
				for j=0,31 do aa(i,j) = total(tmp[i,j*8+ind1],2)
				for j=32,47 do aa(i,j) = total(tmp[i,(j-16)*16+ind2],2)
				for j=48,63 do aa(i,j) = total(tmp[i,(j-32)*32+ind3],2)
			endfor

		store_data,'mvn_sta_DB_MH_M',data={x:tt,y:tmp,v:findgen(1024)}
		store_data,'mvn_sta_DB_MH_M64',data={x:tt,y:aa,v:indgen(64)}
		store_data,'mvn_sta_DB_MH_tot',data={x:tt,y:total(tmp,2)}

		ylim,'mvn_sta_DB_DATA',0,0,1
		ylim,'mvn_sta_DB_MH_tot',0,0,1
		ylim,'mvn_sta_DB_MH_M',0,1024,0
		ylim,'mvn_sta_DB_MH_M64',0,64,0
		zlim,'mvn_sta_DB_MH_M',1,1.e6,1
		zlim,'mvn_sta_DB_MH_M64',1,1.e6,1
		datagap=512.
		options,'mvn_sta_DB_DATA',datagap=datagap
		options,'mvn_sta_DB_MH_M',datagap=datagap
		options,'mvn_sta_DB_MH_M64',datagap=datagap
		options,'mvn_sta_DB_MH_tot',datagap=datagap
		options,'mvn_sta_DB_MH_M','spec',1
		options,'mvn_sta_DB_MH_M64','spec',1
		options,'mvn_sta_DB_MH_M',ytitle='sta!CMHist-DB!C!CMass!Cbin'
		options,'mvn_sta_DB_MH_M64',ytitle='sta!CMHist-DB!C!CMass!Cbin'

		; zero common then load the data
;			common mvn_sta_mh_raw,mvn_sta_mh_ind,mvn_sta_mh_dat & mvn_sta_mh_dat=0 & mvn_sta_mh_ind=-1l


	endif
endif


		options,'mvn_sta_*','x_no_interp',1
		options,'mvn_sta_*','y_no_interp',1

; form combined plots

	tt1=0 & tt2=0 & tt3=0
	get_data,'mvn_sta_CA_P3_A',data=t1
	get_data,'mvn_sta_D4_P4E_A',data=t2
	get_data,'mvn_sta_D2_P4D_A',data=t3
	if size(/type,t1) eq 8 then tt1=1
	if size(/type,t2) eq 8 then tt2=1
	if size(/type,t3) eq 8 then tt3=1
	if (tt1) then store_data,'mvn_sta_A',data=['mvn_sta_CA_P3_A'] 
	if (tt2) then store_data,'mvn_sta_A',data=['mvn_sta_D4_P4E_A']
	if (tt3) then store_data,'mvn_sta_A',data=['mvn_sta_D2_P4D_A'] 
	if (tt1 and tt2) then store_data,'mvn_sta_A',data=['mvn_sta_CA_P3_A','mvn_sta_D4_P4E_A']
	if (tt1 and tt3) then store_data,'mvn_sta_A',data=['mvn_sta_CA_P3_A','mvn_sta_D2_P4D_A']
	if (tt2 and tt3) then store_data,'mvn_sta_A',data=['mvn_sta_D4_P4E_A','mvn_sta_D2_P4D_A']
	if (tt1 and tt2 and tt3) then store_data,'mvn_sta_A',data=['mvn_sta_CA_P3_A','mvn_sta_D4_P4E_A','mvn_sta_D2_P4D_A']
		ylim,'mvn_sta_A',-1,16,0
		zlim,'mvn_sta_A',1,1.e6,1
		options,'mvn_sta_A','spec',1
		options,'mvn_sta_A',ytitle='sta!C!CAnode!Cbin'

	tt1=0 & tt2=0 
	get_data,'mvn_sta_C8_P2_D',data=t1
	get_data,'mvn_sta_D4_P4E_D',data=t2
	if size(/type,t1) eq 8 then tt1=1
	if size(/type,t2) eq 8 then tt2=1
	if (tt1) then store_data,'mvn_sta_D',data=['mvn_sta_C8_P2_D']
	if (tt2) then store_data,'mvn_sta_D',data=['mvn_sta_D4_P4E_D']
	if (tt1 and tt2) then store_data,'mvn_sta_D',data=['mvn_sta_C8_P2_D','mvn_sta_D4_P4E_D']
		ylim,'mvn_sta_D',-1,16,0
		zlim,'mvn_sta_D',1,1.e6,1
		options,'mvn_sta_D','spec',1
		options,'mvn_sta_D',ytitle='sta!C!CDef!Cbin'

; get rid of raw data tplot structures

if not keyword_set(all) then begin
		store_data,delete='mvn_STA_*'
endif


end