;+ ;PROCEDURE: mvn_sta_scpot_load ;PURPOSE: ; Loads scpot into static apid common blocks for science data products - only works if scpot is negative ; ;KEYWORDS: ; check: 0,1 if set, prints diagnostic data ; ;CREATED BY: J. McFadden 15/01/?? ;VERSION: 1 ;LAST MODIFICATION: 16/02/22 ;MOD HISTORY: ; ;NOTES: ; Program assumes that "mvn_sta_l0_load" or "mvn_sta_l2_load,/test" has been run ; Program will use some combination of sta, lpw, swe data to estimate scpot ; Assumption 1 - O2+ ram velocity only valid below 300 km ; Assumption 2 - low energy cutoff only valid when (Ni>100/cc and alt>180km) or in shadow ; Assumption 3 - when the above do not apply, s/c potential set to 0V ; ; TBDs - scpot is not valid 20151127/1:00UT - why is it invalid? ; ; "cutoff" keyword forces the use of c0 cutoff when pot_cutoff_valid=1 which requires enough counts - needed if winds impact potential ; ; Arrays ; pot - potential determined from Vo2-Vram in c6 time base, assume no wind, valid for altitude < alt_vo2 (default=300km) ; pot1 - pot array interpolated to c0 time2 base, used in validity test to determine whether to use heavy ions in pot2 calculation ; pot2 - potential determined from ion cutoff, 3 sweeeps used to increase sensitivity, valid for alt>(200.<(min_alt+35.)>180.), p_only=1 means only protons used ; pot3 - pot2 interpolated to c6 time base ; pot22 - potential determined from total ion cutoff with no averaging - occassionally more accurate than pot2 ; pot33 - pot22 interpolated to c6 time base ; p_only - set to 1 if only protons were used in pot2 determination ; po_only - set to 1 if pot22 is used ; ; pot_alg 0 if invalid ; 1 if Vo2-Vram used (ie. pot) ; 2 if proton only cutoff used (ie. pot2*p_only) ; 3 if all ion cutoff used (ie. pot2*(1-p_only) ; 4 ; ; the complexity of the algorith is due to: ; background coincident events must be ignored ; scattered/trapped ions below s/c potential (especially at potentials <-10V) must be ignored - these are primarily heavies ; lack of protons at <180km to estimate cutoff (assume Vo2=Vram to estimate potential) ; lack of adequate CO2 counts at periapsis for O2-CO2 energy comparison ; marginal statistical counts of protons near periapsis ; transition from negative to positive potential at densities of about 8/cc ; s/c potential is always negative in shadow or in high densities ; ; apid c0 is used for cutoff because ti has ; ;- pro mvn_sta_scpot_load,$ tplot=tplot,$ max_alt=max_alt,$ max_nrg=max_nrg,$ ram_min=ram_min,$ max_den=max_den,$ max_ec=max_ec,$ skip=skip,$ kk2=kk2,$ swe=swe,$ pot_err=pot_err,$ alt_vo2=alt_vo2,$ cutoff=cutoff,$ no_ramdir=no_ramdir common mvn_c6,mvn_c6_ind,mvn_c6_dat common mvn_c0,mvn_c0_ind,mvn_c0_dat common mvn_ca,mvn_ca_ind,mvn_ca_dat if size(mvn_c6_dat,/type) ne 8 or size(mvn_c0_dat,/type) ne 8 or size(mvn_ca_dat,/type) ne 8 then begin print,'Error - c6, c0 and ca data must be loaded' return endif ; mvn_c6_dat.sc_pot[*] = 0. ; mvn_c0_dat.sc_pot[*] = 0. if not keyword_set(max_alt) then max_alt=10000. ; below this altitude s/c potential is assumed negative and algorithm valid when in sunlight ; potential is also assumed negative and algorithm valid in eclipse at higher altitudes ; as a default, we assume the algorithm works everywhere if not keyword_set(max_nrg) then max_nrg=30. ; -1*max_nrg is the maximum negative s/c potential determined by ram flow if not keyword_set(ram_min) then ram_min=0.12 ; -1*ram_min is the minimum negative s/c potential determined by ram flow if not keyword_set(max_den) then max_den=8. ; if in sunlight and >250km and densitymax_ec then s/c potential is set to 0 if not keyword_set(max_ec) then max_ec=50. ; if in sunlight and >250km and densitymax_ec then s/c potential is set to 0 if not keyword_set(pot_err) then pot_err=3. ; ion cutoff pot < pot_err indicates valid potential, when in sunlight and >250km and ec indicates sheath if not keyword_set(alt_vo2) then alt_vo2=300. ; max alt where vel of O2+ can be used to estimate scpot - varies with local time time = (mvn_c6_dat.time + mvn_c6_dat.end_time)/2. npts = n_elements(time) qf_valid = (mvn_c6_dat.quality_flag and 3) eq 0 c6_mode = mvn_c6_dat.mode cols=get_colors() if not keyword_set(kk2) then kk2 = mvn_sta_get_kk2(time[0]) ; correcting for ion suppression introduces errors ; kk2=0. ;print,'kk2=',kk2 ;********************************************************** ; get swe data if swe keyword is set if keyword_set(swe) then begin get_data,'mvn_swe_spec_dens',data=tmp if size(tmp,/type) ne 8 then begin mvn_swe_load_l2,/spec mvn_swe_sumplot,/loadonly mvn_swe_sc_pot,/over,badval=0 mvn_swe_n1d,/mom,background=1e3,minden=1e-4 get_data,'mvn_swe_spec_dens',data=tmp if size(tmp,/type) eq 8 then begin swe_den = interp(tmp.y,tmp.x,time) get_data,'mvn_swe_spec_temp',data=tmp2 swe_te = interp(tmp2.y,tmp2.x,time) swe_den = swe_den*(10./(swe_te>10.)) endif else begin print,'Error - swe data not available, sc potential not loaded' return endelse endif else begin swe_den = interp(tmp.y,tmp.x,time) get_data,'mvn_swe_spec_temp',data=tmp2 swe_te = interp(tmp2.y,tmp2.x,time) swe_den = swe_den*(10./(swe_te>10.)) endelse if keyword_set(tplot) then store_data,'swe_den_scaled_te',data={x:time,y:swe_den} if keyword_set(tplot) then ylim,'swe_den_scaled_te',1,100,1 endif ;********************************************************** ; estimate density and characteristic energy without s/c potential corrections get_data,'mvn_sta_test_density3',data=tmp if size(tmp,/type) ne 8 or (not keyword_set(skip)) then begin get_4dt,'ec_4d','mvn_sta_get_c6',mass=[0,1.7],name='mvn_sta_test_ec_p',energy=[0,30000.] ylim,'mvn_sta_test_ec_p',1,1.e4,1 options,'mvn_sta_test_ec_p',colors=cols.green,ytitle='sta!Cc6!Ctest!Cec_p',labels='H+He+',labpos=1000 get_4dt,'ec_4d','mvn_sta_get_c6',mass=[10,100],name='mvn_sta_test_ec_o',energy=[0,30000.] ylim,'mvn_sta_test_ec_o',1,1.e4,1 options,'mvn_sta_test_ec_o',colors=cols.red,ytitle='sta!Cc6!Ctest!Cec_o',labels='O+O++',labpos=10 get_4dt,'n_4d','mvn_sta_get_c6',mass=[.3,5],name='mvn_sta_test_density1',m_int=1.,energy=[0,30000.] options,'mvn_sta_test_density1',colors=cols.green,labels='H+He+',labpos=100 get_4dt,'n_4d','mvn_sta_get_c6',mass=[10.,100.],name='mvn_sta_test_density2',energy=[0,30000.] options,'mvn_sta_test_density2',colors=cols.red,labels='O+O++',labpos=1 get_4dt,'c_4d','mvn_sta_get_c6',mass=[.3,5],name='mvn_sta_test_cnts_p_low',m_int=1.,energy=[0,12.] options,'mvn_sta_test_cnts_p_low',ytitle='sta!Cc6!CH+!C<12eV!Ccnts',psym=-1,colors=cols.green,labels='H+He+',labpos=100 ylim,'mvn_sta_test_cnts_p_low',1,1000,1 get_4dt,'c_4d','mvn_sta_get_c6',mass=[12,100],name='mvn_sta_test_cnts_o_med',m_int=32,energy=[0,30.] ; this energy range may be too high options,'mvn_sta_test_cnts_o_med',colors=cols.cyan,ytitle='sta!Cc6!CO+,O++!C<30eV!Ccnts',psym=-1 ylim,'mvn_sta_test_cnts_o_med',1,1000,1 get_4dt,'c_4d','mvn_sta_get_c6',mass=[12,100],name='mvn_sta_test_cnts_o_low',m_int=32,energy=[0,12.] ; this energy range may be too high options,'mvn_sta_test_cnts_o_low',ytitle='sta!Cc6!CO+,O++!C<12eV!Ccnts',psym=-1,colors=cols.red,labels='O+O++',labpos=10 ylim,'mvn_sta_test_cnts_o_low',1,1000,1 get_4dt,'c_4d','mvn_sta_get_c6',mass=[12,100],name='mvn_sta_test_cnts_o',m_int=32,energy=[30,40000.] ; this energy range may be too high options,'mvn_sta_test_cnts_o',colors=cols.cyan,ytitle='sta!Cc6!CO+,O++!C>30eV!Ccnts',psym=-1 ylim,'mvn_sta_test_cnts_o',1,10000,1 get_4dt,'v_4d','mvn_sta_get_ca',name='mvn_sta_ca_vel3',m_int=1.,energy=[0,40000.] get_data,'mvn_sta_ca_vel3',data=tmp store_data,'mvn_sta_ca_vel',data={x:tmp.x,y:(total(tmp.y*tmp.y,2))^.5} options,'mvn_sta_ca_vel','sta!Cca!CH+!Cvel!Ckm/s' ylim,'mvn_sta_ca_vel',10,500,1 get_data,'mvn_sta_test_density1',data=tmp1 get_data,'mvn_sta_test_density2',data=tmp2 store_data,'mvn_sta_test_density3',data={x:tmp1.x,y:tmp1.y+tmp2.y} ylim,'mvn_sta_test_density3',1,100.,1 endif get_data,'mvn_sta_test_density1',data=tmp1 get_data,'mvn_sta_test_density2',data=tmp2 get_data,'mvn_sta_test_density3',data=tmp ind = where(tmp.y gt 0.,count) if count gt 1 then begin den = interp(tmp.y[ind],tmp.x[ind],time) p_den = interp(tmp1.y[ind],tmp1.x[ind],time) o_den = interp(tmp2.y[ind],tmp2.x[ind],time) endif else begin print,'Error - STATIC c6 data not loaded' return endelse get_data,'mvn_sta_test_ec_p',data=tmp p_ec = interp(tmp.y[ind],tmp.x[ind],time) get_data,'mvn_sta_test_ec_o',data=tmp o_ec = interp(tmp.y[ind],tmp.x[ind],time) get_data,'mvn_sta_test_cnts_p_low',data=tmp p_low = interp(tmp.y[ind]+.1,tmp.x[ind],time) get_data,'mvn_sta_ca_vel',data=tmp p_vel = interp(tmp.y,tmp.x,time) get_data,'mvn_sta_test_cnts_o',data=tmp o_cnt=interp(tmp.y,tmp.x,time) get_data,'mvn_sta_test_cnts_o_med',data=tmp o_med=interp(tmp.y,tmp.x,time) get_data,'mvn_sta_test_cnts_o_low',data=tmp o_low=interp(tmp.y,tmp.x,time) store_data,'mvn_sta_test_density',data=['mvn_sta_test_density1','mvn_sta_test_density2'] options,'mvn_sta_test_density',ytitle='sta!C!Ctest!C!Cdensity' ylim,'mvn_sta_test_density',.1,1000,1 store_data,'mvn_sta_test_ec',data=['mvn_sta_test_ec_p','mvn_sta_test_ec_o'] ylim,'mvn_sta_test_ec',1,2000,1 options,'mvn_sta_test_ec',ytitle='sta!Cc6!CH+,O+!C!CEc' store_data,'mvn_sta_test_cnts_low',data=['mvn_sta_test_cnts_o_low','mvn_sta_test_cnts_p_low'] ylim,'mvn_sta_test_cnts_low',1,2000,1 options,'mvn_sta_test_cnts_low',ytitle='sta!Cc6!C<10eV!C!CCnts' ;********************************************************** ; determine when sc is in shadow and wake mars_radius = 3386. pos = mvn_c6_dat.pos_sc_mso if min((pos[*,0]^2+pos[*,1]^2+pos[*,2]^2)^.5) lt mars_radius then begin print,'Error - missing ephemeris data in c6' maven_orbit_tplot,/current,result=foo,/LOADONLY,eph=eph if size(foo,/type) eq 8 then begin pos[*,0] = interp(foo.x*mars_radius,foo.t,time) pos[*,1] = interp(foo.y*mars_radius,foo.t,time) pos[*,2] = interp(foo.z*mars_radius,foo.t,time) print,'Alert: mvn_c6_dat.pos_sc_mso has values < mars_radius, using maven_orbit_tplot to determine s/c position' endif else begin print,'Error: sc position not determined - aborting mvn_sta_scpot_load!!!!' return endelse endif shadow = ((pos[*,1]^2+pos[*,2]^2)^.5 - 200.) lt mars_radius and pos[*,0] lt 0. if keyword_set(tplot) then store_data,'mvn_shadow',data={x:time,y:shadow} if keyword_set(tplot) then ylim,'mvn_shadow',-1,2,0 wake = ((pos[*,1]^2+pos[*,2]^2)^.5) lt (mars_radius+500.-pos[*,0]*sin(10./!radeg)) and (pos[*,0] lt 0.) if keyword_set(tplot) then store_data,'mvn_wake',data={x:time,y:wake} if keyword_set(tplot) then ylim,'mvn_wake',-1,2,0 if keyword_set(tplot) then store_data,'mvn_yz',data={x:time,y:(pos[*,1]^2+pos[*,2]^2)^.5 - mars_radius} if keyword_set(tplot) then ylim,'mvn_yz',-4000,4000,0 alt = (pos[*,0]^2+pos[*,1]^2+pos[*,2]^2)^.5 - mars_radius ; min_alt = min(alt) > 150. min_alt = min(alt) ;********************************************************** ; determine sc velocity relative to planet - this is ram velocity get_data,'V_sc_MAVEN_APP',data=tmp81 trange=timerange() count=1 if size(tmp81,/type) eq 8 then if n_elements(tmp81.x) gt 10 then ind = where((tmp81.x lt trange[0]) or (tmp81.x gt trange[1]),count) ; 10 is sort of arbitrary if size(tmp81,/type) eq 8 and count eq 0 then begin vtot = (total(tmp81.y*tmp81.y,2))^.5 vel = interp(vtot,tmp81.x,mvn_c6_dat.time+2.) endif else begin mkernels = mvn_spice_kernels(/load) if not keyword_set(no_ramdir) then mvn_ramdir,trange,frame=['MAVEN_APP'],dt=1,/polar get_data,'V_sc_MAVEN_APP',data=tmp81 if size(tmp81,/type) eq 8 then begin vtot = (total(tmp81.y*tmp81.y,2))^.5 vel = interp(vtot,tmp81.x,mvn_c6_dat.time+2.) endif else begin ; approximate velocity from position motion - used prior to 20160805 vel = total((pos[1:npts-1,*] - pos[0:npts-2,*])^2,2)^.5/(time[1:npts-1]-time[0:npts-2]) vel = [vel[0],[vel]] print,'ERROR - mvn_ramdir did not run!!! - using maven_orbit_tplot to estimate velocity ignoring planet rotation' print,'ERROR - this approximation will introduce sc_pot error at periapsis!!!' endelse endelse if keyword_set(tplot) then store_data,'mvn_sta_scpot_vel',data={x:time,y:vel} nrg_offset=+.0 ; this is included incase sensor HV has an offset ;********************************************************** ; estimate sc potential using c6 O2+ velocity data ; this estimate is only valid below about 300 km pot = fltarr(npts) vd = fltarr(npts) vth = fltarr(npts) ano = fltarr(npts) ms = mvn_c6_dat.mass*32. pd_arr = fltarr(npts,6) mass_o2 = [28,40.] m_o2 = 32. min_o2 = 0 nbad=0 for i=0l,npts-1 do begin if mvn_c6_dat.energy[mvn_c6_dat.swp_ind[i],31,0] lt 2. then begin ; the lower energy limit of the sweep must be less then 2 eV ; the following was used prior to 20170502 - replaced by more accurate vb_pot_4d.pro if 0 then begin cnts=reform(mvn_c6_dat.data[i,*,*]) data=reform(mvn_c6_dat.eflux[i,*,*]) mass=reform(mvn_c6_dat.mass_arr[mvn_c6_dat.swp_ind[i],*,*]) ind = where(mass lt 25 or mass gt 40,count) ; ind = where(mass lt 30 or mass gt 34,count) ; this was used during some testing to see if it worked better if count ne 0 then data[ind]=0. if count ne 0 then cnts[ind]=0. nrg = (reform(mvn_c6_dat.energy[mvn_c6_dat.swp_ind[i],*,*]) + nrg_offset) > 0.01 ; nrg_offset is currently 0, but might be needed ; anode array may be useful in the future for eliminating data when flow not in ram direction min_ca = min(abs(mvn_ca_dat.time +2. - time[i]),ind_ca) anode = reform(replicate(1.,2)#reform(total(reform(mvn_ca_dat.data[ind_ca,*,*],16,4,16),2),16*16),32,16) ; the following kluge to max_nrg was needed for cases with nightside vdis at 10 eV - example 20160107 0317UT max_nrg2 = max_nrg peak_eflx_ca = max(anode[*,7],ind_peak_ca) peak_eflx_c6 = max(total(data,2),ind_peak_c6) if ind_peak_ca gt (ind_peak_c6+1) then begin mintmp = min(abs(nrg[ind_peak_ca,0]-nrg[*,0]),minind) maxnrg = nrg(minind-2,0) if alt[i] le 180. and total(cnts[((minind-2)>1):31,*]) ge 10 then max_nrg2=nrg(((minind-2)>1),0) endif ; bin10ev = min(abs(nrg[*,0]-10.),ind10ev) ; if alt[i] le 180. and total(cnts[ind10ev:31,*]) ge 10 then max_nrg2=10. ; bin6ev = min(abs(nrg[*,0]-6.),ind6ev) ; if alt[i] le 180. and total(cnts[ind6ev:31,*]) ge 10 then max_nrg2=6. ind = where(nrg gt max_nrg2,count) if count ne 0 then data[ind]=0. if count ne 0 then cnts[ind]=0. ind = where(nrg lt 1.,count) if count ne 0 then data[ind]=0. if count ne 0 then cnts[ind]=0. v = (2.*nrg/ms)^.5 maxcnt = max(total(data,2),mind) data[0:((mind-4)>0),*]=0. ; data[((mind+4)<31):31,*]=0. en_peak=nrg[mind,0] data = data*(exp((kk2/nrg)^2) < 10.) ; this conservative correction may be inadequate between 20150201 and 20151001 vd[i] = total(data)/(total(data/v)>1.e-20) vth[i] = (total(data*(v-vd[i])^2/v)/(total(data/v)>1.e-20))^.5 pot[i] = 0.5*ms*(vd[i]^2-vel[i]^2) > ram_min ; ram_min=0.12 the minimum allowed negative potential endif else begin ; pot is not valid above alt_vo2 so set the default if alt[i] gt alt_vo2 then begin cnts = 6. ; this avoids an error with code 38 lines down pot[i] = max_nrg vd[i] = -1. endif else begin dat6 = mvn_sta_get_c6(mvn_c6_dat.time[i]) if dat6.valid then begin if dat6.mode eq 7 then engy_o2 = [0,11] else engy_o2 = [0,max_nrg+5.] ; this is somewhat arbitrary cnts = cb_4d(dat6,mass=mass_o2,m_int=m_o2,energy=engy_o2,mincnt=min_o2) ; not quite consistent with vb_pot_4d.pro ; 3 iterations give an accurate potential ; these iterations avoid small errors between calculated potential from v-drift dat7 = dat6 & dat7.sc_pot=0. vd7 = vb_pot_4d(dat7,mass=mass_o2,m_int=m_o2,energy=engy_o2,mincnt=min_o2,/no_co2) pd7 = 0.5*ms*(vd7^2-vel[i]^2) vd[i] = vd7 ; this works better for high scpot potentials if pd7 gt 5. then begin pd7 = .8*pd7 dat8 = dat6 & dat8.sc_pot = -pd7 vd8 = vb_pot_4d(dat8,mass=mass_o2,m_int=m_o2,energy=engy_o2,mincnt=min_o2,/no_co2) pd8 = 0.5*ms*(vd8^2-vel[i]^2) pd9 = 0. pd10 = 0 pd11 = 0 endif else begin dat8 = dat6 & dat8.sc_pot = -pd7 vd8 = vb_pot_4d(dat8,mass=mass_o2,m_int=m_o2,energy=engy_o2,mincnt=min_o2,/no_co2) pd8 = 0.5*ms*(vd8^2-vel[i]^2) dat9 = dat6 & dat9.sc_pot = -pd7-pd8 vd9 = vb_pot_4d(dat9,mass=mass_o2,m_int=m_o2,energy=engy_o2,mincnt=min_o2,/no_co2) pd9 = 0.5*ms*(vd9^2-vel[i]^2) dat10 = dat6 & dat10.sc_pot = -pd7-pd8-pd9 vd10 = vb_pot_4d(dat10,mass=mass_o2,m_int=m_o2,energy=engy_o2,mincnt=min_o2,/no_co2) pd10 = 0.5*ms*(vd10^2-vel[i]^2) dat11 = dat6 & dat11.sc_pot = -pd7-pd8-pd9-pd10 vd11 = vb_pot_4d(dat11,mass=mass_o2,m_int=m_o2,energy=engy_o2,mincnt=min_o2,/no_co2) pd11 = 0.5*ms*(vd11^2-vel[i]^2) endelse ; these are gathered for testing pd_arr[i,0]=pd7 pd_arr[i,1]=pd8 pd_arr[i,2]=pd9 pd_arr[i,3]=pd10 pd_arr[i,4]=pd11 pd_arr[i,5]=pd7+pd8+pd9+pd10+pd11 ; use of !values.f_nan in the following lines can causes problems ; if finite(pd7+pd8+pd9) then pot[i] = ((pd7+pd8+pd9) > ram_min) else pot[i] = !values.f_nan ; ram_min=0.12 the minimum allowed negative potential if finite(pd7+pd8+pd9) then pot[i] = ((pd7+pd8+pd9) > ram_min) else pot[i] = !values.f_nan ; ram_min=0.12 the minimum allowed negative potential endif else begin vd[i] = !values.f_nan pot[i] = !values.f_nan endelse endelse endelse ; not sure what these lines do -- need to document if total(cnts) lt 5. and alt[i] le 180. then pot[i]=pot[(i-1)>0] if total(cnts) lt 5. and alt[i] gt 180. then pot[i]=pot[(i-1)>0] ; if total(cnts) lt 5. and alt[i] gt 180. and den[i] lt 1. then pot[i]=pot[(i-1)>0] ; if total(cnts) lt 5. and alt[i] gt 180. and den[i] ge 1. then pot[i]=max_nrg ; if alt[i] gt 180. then pot[i]=pot[i] > .8 if (mvn_c6_dat.quality_flag[i] and 192) gt 0 then pot[i]=pot[(i-1)>0] if (mvn_c6_dat.quality_flag[i] and 3) gt 0 then pot[i]=0 endif endfor if keyword_set(tplot) then store_data,'mvn_sta_c6_pot_vd',data={x:time,y:vd} if keyword_set(tplot) then options,'mvn_sta_c6_pot_vd',colors=cols.green if keyword_set(tplot) then store_data,'mvn_sta_c6_pot_vth',data={x:time,y:vth} if keyword_set(tplot) then store_data,'mvn_sta_c6_pot_ec',data={x:time,y:0.5*ms*vd^2} if keyword_set(tplot) then store_data,'mvn_sta_c6_pd_arr',data={x:time,y:pd_arr} if keyword_set(tplot) then options,'mvn_sta_c6_pd_arr',colors=[1,2,3,4,6,0] if keyword_set(tplot) then store_data,'mvn_sta_c6_vo2_pot',data={x:time,y:pot} if keyword_set(tplot) then options,'mvn_sta_c6_vo2_pot',ylog=1,yrange=[.1,10] ;print,minmax(vth/vd) ;********************************************************** ; estimate sc potential using c0 data time2 = (mvn_c0_dat.time + mvn_c0_dat.end_time)/2. npts2 = n_elements(time2) pot2 = fltarr(npts2) pot22 = fltarr(npts2) pot2_cutoff_cnts = fltarr(npts2) pot2_p_cutoff_cnts = fltarr(npts2) pot2_p_cutoff_cnts0 = fltarr(npts2) pot2_p_cutoff_cnts1 = fltarr(npts2) pot2_p_cutoff_cnts2 = fltarr(npts2) pot2_p_cutoff_cnts3 = fltarr(npts2) p_cutoff = intarr(npts2) modep_reject = fltarr(npts2) modem_reject = fltarr(npts2) pot_bad = fltarr(npts2) cnts_lt_pot = fltarr(npts2) swp_ind = mvn_c0_dat.swp_ind ; pot0 = interp(pot,time,time2) ; pot0 is redundant, same as pot1 get_data,'mvn_sta_test_density3',data=tmp3 den2 = interp(tmp3.y,tmp3.x,time2) ind_cutoff = fltarr(npts2) p_only = fltarr(npts2) oxgt30eV=mvn_c0_dat ; get_data,'alt',data=tmp ; alt2 = interp(tmp.y,tmp.x,time2) alt2 = interp(alt,time,time2) pot1 = interp(pot,time,time2) ; pot1 determined from c6 O2+ velocity data sha2 = interp(shadow,time,time2) o_cnt2 = interp(o_cnt,time,time2) ; 30-40,000 eV p_low2 = interp(p_low,time,time2) ; 0-12 eV o_low2 = interp(o_low,time,time2) ; 0-12 eV o_med2 = interp(o_med,time,time2) ; 0-30 eV ; Determine pot2, the following is stored as 'mvn_sta_c6_O2+_sc_pot_h+' ; Designed to not use oxygen during s/c charging events with scpot<-5 eV since sputtered-trapped O2+ ions at 3000 for i=1l,npts2-2 do begin if mvn_c0_dat.energy[mvn_c0_dat.swp_ind[i],63,0] lt 2. then begin ; the lower energy limit of the sweep must be less then 2 eV modep = (mvn_c0_dat.mode[i] eq mvn_c0_dat.mode[i+1]) and (((mvn_c0_dat.quality_flag[i+1] and 64)/64) eq 0) modem = (mvn_c0_dat.mode[i] eq mvn_c0_dat.mode[i-1]) and (((mvn_c0_dat.quality_flag[i-1] and 64)/64) eq 0) ; if (pot1[i] lt 3.) or ((alt2[i] gt 400.) and sha2[i] and (o_cnt2[i] lt 25)) then begin ; used prior to 20161130, checked in on ? ; if (pot1[i] lt 3.) or ((alt2[i] gt 400.) and sha2[i] and (o_cnt2[i] lt 25) and (p_low2[i] lt 20)) then begin ; o_cnt2 >30eV, changed to only use protons when there are enough proton counts, fixed problem on 20160915/2307UT ; if (pot1[i] lt 3.) or ((alt2[i] gt 400.) and sha2[i] and (o_cnt2[i] lt 25) and (p_low2[i] lt 20)) $ ; or ((alt2[i] gt 400.) and not sha2[i] and (o_low2[i] gt 5)) then begin ; o_cnt2 >30eV, changed to only use protons when there are enough proton counts, fixed problem on 20160915/2307UT ; this had problems on 20170223/1508-1528UT ; if ((pot1[i] lt 3.) and (o_cnt2[i] lt 100)) or ((alt2[i] gt 400.) and sha2[i] and (o_cnt2[i] lt 25) and (p_low2[i] lt 20)) then begin ; o_cnt2 >30eV, changed to only use protons when there are enough proton counts, fixed problem on 20160915/2307UT ; changed 20170228 -- needs testing, if ((pot1[i] lt 3.) and (o_cnt2[i] lt 100)) $ or ((alt2[i] gt 400.) and sha2[i] and (o_cnt2[i] lt 25) and (p_low2[i] lt 20)) $ ; or ((alt2[i] gt (min_alt+30.)) and (den2[i] lt 3000.)) then begin ; changed to use all ions when the den<3000., modified 20170330 to use (min_alt+30.) - didn't seem to matter or ((alt2[i] gt 150.) and (den2[i] lt 3000.)) then begin ; changed to use all ions when the den<3000. dat1c=total(reform(mvn_c0_dat.data[i,*,*]),2) dat1p=total(reform(mvn_c0_dat.data[i+1,*,*]),2) dat1m=total(reform(mvn_c0_dat.data[i-1,*,*]),2) endif else begin ; this eliminates sputtereed O2+ when s/c charges to >5V, and just uses protons when there are enough counts dat1c=reform(mvn_c0_dat.data[i,*,0]) dat1p=reform(mvn_c0_dat.data[i+1,*,0]) dat1m=reform(mvn_c0_dat.data[i-1,*,0]) p_only[i]=1 endelse dat3c=reform(mvn_c0_dat.data[i,*,0]) dat3p=reform(mvn_c0_dat.data[i+1,*,0]) dat3m=reform(mvn_c0_dat.data[i-1,*,0]) dat7c=reform(mvn_c0_dat.data[i,*,1]) dat7p=reform(mvn_c0_dat.data[i+1,*,1]) dat7m=reform(mvn_c0_dat.data[i-1,*,1]) ; if p_low2[i] gt 30 then begin ; old code ; modep=0 ; modem=0 ; endif dat2c = shift(dat1c,1) dat2p = shift(dat1p,1) dat2m = shift(dat1m,1) dat4c = shift(dat3c,1) dat4p = shift(dat3p,1) dat4m = shift(dat3m,1) dat5c = shift(dat1c,2) dat6c = shift(dat3c,2) dat8c = shift(dat7c,1) dat8p = shift(dat7p,1) dat8m = shift(dat7m,1) dat9c = shift(dat7c,2) data1 = dat1c+modep*dat1p+modem*dat1m data2 = dat2c+modep*dat2p+modem*dat2m data3 = dat3c+modep*dat3p+modem*dat3m data7 = dat7c+modep*dat7p+modem*dat7m energy0 = reform(mvn_c0_dat.energy[swp_ind[i],*,0]) energy = (energy0 + nrg_offset) > .01 denergy = reform(mvn_c0_dat.denergy[swp_ind[i],*,0])*energy/energy0 ; determine the ion cutoff energy ind0 = 0 ; ind0 = where (data1 ge 1. and data2 ge 1. and (data1+data2) ge 4.,count) ; changed 20170228 ; ind0 = where ( ( (data1 ge 1.) and (data2 ge 1.) and ((data1+data2) ge 4.) and (dat5c ge 1.) ) $ ; or ( ((dat3c+dat4c) ge 3.) and (dat5c ge 1.) ) $ ; or ( ((dat3c+dat4c) ge 2.) and (dat5c ge 2.) ) ,count) ; changed 20170412 ind0 = where ( ( ((dat1c ge 1.) or (data1 ge 2.)) and ((dat2c ge 1.) or (data2 ge 2.)) and ((dat1c+dat2c) ge 2.) and ((data1+data2) ge 4.) and ((data1+data2+dat5c) ge 5.) ) $ or ( ((dat1c ge 1.) or (data1 ge 2.)) and ((dat2c ge 1.) or (data2 ge 2.)) and ((dat1c+dat2c) ge 2.) and ((data1+data2) ge 3.) and ((dat1c+dat2c+dat5c) ge 4.) ) $ or ( (dat3c ge 1.) and ((dat3c+dat4c) ge 3.) and ((dat3c+dat4c+dat6c) ge 4.) ) $ or ( (dat3c ge 1.) and ((dat3c+dat4c) ge 2.) and ((dat3c+dat4c+dat6c) ge 5.) ) ,count) ; if count gt 0 then mind0=max(ind0) ; would crashed, changed 20170812 if count gt 0 then mind0=max(ind0) else mind0=1 mind0 = mind0 > 1 ; 20200528 note: the above algorithm is attempting to pick out the proton cutoff at the statistical noise level ; coincident background from O2+ can produce enough 'c0' noise counts in mass_bin=0 create a false proton cutoff ; under these circumstances the vo2-pot (pot) should be selected, but may not be if the keyword "cutoff" is set. ; don't average over time if dat1c >10 counts, added 20170331 ; if count gt 0 and (dat1c[mind0]+dat2c[mind0]+dat2c[(mind0-1)>0]) gt 10. then begin ; old method if count gt 0 and (dat1c[mind0]+dat2c[mind0]) ge 10. then begin ; if dat1c[mind0] eq 0 then mind0=(mind0-1)>1 ; old method data1=dat1c data2=dat2c data3=dat3c data7=dat7c modep = 0 modem = 0 endif ; interpolate ; if energy[mind0] lt max_nrg and ((alt2[i] gt 180.) or ((data1[mind0]+data2[mind0]) gt 20.)) then begin ; if energy[mind0] lt max_nrg and ((alt2[i] gt (min_alt+25.)) or ((data1[mind0]+data2[mind0]) ge 20.)) then begin ; changed 20200817 to the following line if energy[mind0] lt max_nrg and ((alt2[i] gt (200.<(min_alt+35.)<180.)) or ((data1[mind0]+data2[mind0]) ge 20.)) then begin ; changed 20200817 ; e0 = energy[mind0] - denergy[mind0]/2. ; d1 = data1[mind0] ; e1 = energy[mind0] - e0 ; d2 = data2[mind0-1] ; e2 = energy[mind0-1] - e0 ; scale = 0. > (1-(d1/d2)*(e2/e1)^1) < 1. ; pot2[i]= (energy[mind0] - denergy[mind0]/2. + denergy[mind0]*scale) < max_nrg e0 = energy[mind0] d0 = data1[mind0] e1 = energy[mind0-1] d1 = data1[mind0-1] pot2[i] = e0 + .5*(e1-e0)*(d1-d0)/(d1+d0+3.) ; the "3" here just skews it toward mind0 for low count rates if d1 le 3. then pot2[i]=energy[mind0] ; don't interpolate if d1 is small p_cutoff[i] = (data3[mind0]+data3[mind0-1]) gt (data7[mind0]+data7[mind0-1]) ; this determines whether protons dominate the ion cutoff endif else pot2[i] = max_nrg ; endif else pot2[i] = 2.5 ; if alt2[i] ge 300. and den2[i] lt 1. then pot2[i] = 2.5 if (mvn_c0_dat.quality_flag[i] and 192) gt 0 then pot2[i]=pot2[(i-1)>0] if (mvn_c0_dat.quality_flag[i] and 3) gt 0 then pot2[i]=max_nrg ; if (alt2[i] lt (min_alt+25.)) then pot2[i]=max_nrg ; changed 20200817 to following line, since periapsis is no longer determined by density corridor if (alt2[i] lt (200.<(min_alt+35.)>180.)) then pot2[i]=max_nrg ; changed 20200817, this basically duplicates the "pot_cutoff_valid =" line far below endif else pot2[i] = max_nrg ;************************************** ; this added 20170330 this is used to determine validity of low energy cutoff for along track wind calculations if pot2[i] lt (max_nrg*.9 > 10.) then begin nrg = reform(mvn_c0_dat.energy[swp_ind[i],*,0]) ; minval = min(abs(nrg-pot2[i]),ind77) ind77=mind0 ; mind0 is the bin number of the closest energy step to the potential determined above ; minval = min(abs(nrg-(2.7+pot0[i])*.70),ind87) ; pot0 is redundant, same as pot1 minval = min(abs(nrg-(2.7+pot1[i])*.70),ind87) if alt2[i] lt alt_vo2 then ind77=ind77>ind87 ; this makes sure we don't count o2+ stragglers pot2_cutoff_cnts[i] = total(dat1c[0>(ind77-2):(ind77+1)<63]+modep*dat1p[0>(ind77-2):(ind77+1)<63]+modem*dat1m[0>(ind77-2):(ind77+1)<63]) pot2_p_cutoff_cnts[i] = total(dat3c[0>(ind77-1):(ind77+1)<63]+modep*dat3p[0>(ind77-1):(ind77+1)<63]+modem*dat3m[0>(ind77-1):(ind77+1)<63]) ; why is this more restrictive? why "-1" and not "-2" pot2_p_cutoff_cnts0[i] = (total(dat3c[0>(mind0-2):mind0])+modep*total(dat3p[0>(mind0-1):mind0])+modem*total(dat3m[0>(mind0-1):mind0]))/(1+modep+modem) pot2_p_cutoff_cnts1[i] = total(dat3c[63<(mind0+1):63]+modep*dat3p[63<(mind0+1):63]+modem*dat3m[63<(mind0+1):63]) pot2_p_cutoff_cnts2[i] = total(dat3c[0>(mind0-1):mind0]) pot2_p_cutoff_cnts3[i] = total(dat3c[63<(mind0+1):63]) modep_reject[i] = modep modem_reject[i] = modem ind_cutoff[i]=ind77 endif ;************************************** endfor if keyword_set(tplot) then store_data,'mvn_sta_c0_p_cutoff_pot2',data={x:time2,y:pot2} ; for testing only if keyword_set(tplot) then store_data,'modep_reject',data={x:time2,y:modep_reject} if keyword_set(tplot) then options,'modep_reject',yrange=[-1,2] if keyword_set(tplot) then store_data,'modem_reject',data={x:time2,y:modem_reject} if keyword_set(tplot) then options,'modem_reject',yrange=[-1,2] if keyword_set(tplot) then store_data,'pot2_p_cutoff_cnts0',data={x:time2,y:pot2_p_cutoff_cnts0} if keyword_set(tplot) then options,'pot2_p_cutoff_cnts0',ylog=1,yrange=[1,100],psym=-1 if keyword_set(tplot) then store_data,'pot2_p_cutoff_cnts1',data={x:time2,y:pot2_p_cutoff_cnts1} if keyword_set(tplot) then options,'pot2_p_cutoff_cnts1',ylog=1,yrange=[1,100],psym=-1 if keyword_set(tplot) then store_data,'pot2_p_cutoff_cnts2',data={x:time2,y:pot2_p_cutoff_cnts2} if keyword_set(tplot) then options,'pot2_p_cutoff_cnts2',ylog=1,yrange=[1,100],psym=-1 if keyword_set(tplot) then store_data,'pot2_p_cutoff_cnts3',data={x:time2,y:pot2_p_cutoff_cnts3} if keyword_set(tplot) then options,'pot2_p_cutoff_cnts3',ylog=1,yrange=[1,100],psym=-1 if keyword_set(tplot) then store_data,'ind_cutoff',data={x:time2,y:ind_cutoff} if keyword_set(tplot) then store_data,'p_only',data={x:time2,y:p_only} ylim,'p_only',-1,2,0 if keyword_set(tplot) then store_data,'p_cutoff',data={x:time2,y:p_cutoff} ylim,'p_cutoff',-1,2,0 ; tplot,['pot2_p_cutoff_cnts0','pot2_p_cutoff_cnts1','pot2_p_cutoff_cnts2','pot2_p_cutoff_cnts3','mvn_sta_c6_O2+_lpw_sc_pot_all','mvn_sta_c6_pot_cutoff_valid','mvn_sta_c0_P1A_L_E_pot','mvn_sta_c0_P1A_H_E_pot'] ;********************************************************** ; If counts <10eV exceed 10, density>8., and alt>400 then use the total count cutoff without averaging ; Determine pot22, the following is stored as 'mvn_sta_c6_O2+_sc_pot_o+h+' for i=1l,npts2-2 do begin if mvn_c0_dat.energy[mvn_c0_dat.swp_ind[i],63,0] lt 2. then begin ; the lower energy limit of the sweep must be less then 2 eV if ((alt2[i] gt 400.) and den2[i] gt 8. and ((o_low2[i]+p_low2[i]) gt 10)) then begin ; o_cnt2 >30eV, changed to only use protons when there are enough proton counts, fixed problem on 20160915/2307UT energy0 = reform(mvn_c0_dat.energy[swp_ind[i],*,0]) dat22=total(reform(mvn_c0_dat.data[i,*,*]),2) minval = min(abs(energy0-max_nrg),max_nrg_ind) dat22[0:max_nrg_ind]=0. dat23=shift(dat22,1) mincnts = 3. > (o_low2[i]+p_low2[i])/20. > (o_cnt2[i]+o_med2[i])/40. ; higher threshold to eliminate sputtered trapped O+ background ind22 = where ((dat22 ge 1.) and (dat23 ge 1.) and ((dat22+dat23) ge mincnts),count) if count ge 1 then begin ; interpolate between energy bins if 1 then begin mind22 = max(ind22) d1 = dat22[mind22] d2 = dat22[mind22-1] e1 = energy0[mind22] e2 = energy0[mind22-1] pot22[i] = e1 + .5*(e2-e1)*(1.-(d1>3.)/(d2>d1>3.)) - .5*(e2-e1)*(1.-3./(d1>3.)) ; pot22[i] = e1 endif if 0 then begin mind22 = max(ind22) energy = (energy0 + nrg_offset) > .01 denergy = reform(mvn_c0_dat.denergy[swp_ind[i],*,0])*energy/energy0 e0 = energy[mind22] - denergy[mind22]/2. d1 = dat22[mind22] e1 = energy[mind22] - e0 d2 = dat22[mind22-1] e2 = energy[mind22-1] - e0 scale = 0. > (1-(d1/d2)*(e2/e1)^1) < 1. pot22[i] = (energy[mind22] - denergy[mind22]/2. + denergy[mind22]*scale) < max_nrg endif endif else pot22[i]=max_nrg endif else pot22[i] = max_nrg if (mvn_c0_dat.quality_flag[i] and 192) gt 0 then pot22[i]=pot22[(i-1)>0] if (mvn_c0_dat.quality_flag[i] and 3) gt 0 then pot22[i]=max_nrg endif else pot22[i] = max_nrg endfor ;********************************************************** ; eliminate data when there are significant counts below scpot (almost never) or when the mode does not sweep to energies lt 2. for i=1l,npts2-2 do begin if mvn_c0_dat.energy[mvn_c0_dat.swp_ind[i],63,0] lt 2. then begin ; the lower energy limit of the sweep must be less then 2 eV energy0 = reform(mvn_c0_dat.energy[swp_ind[i],*,0]) pot_tmp = pot2[i] < pot22[i] ind = where(energy0 lt pot_tmp/1.1,count) if count gt 1 then cnts_lt_pot[i] = total(mvn_c0_dat.data[i,ind,*]) if cnts_lt_pot[i] gt 5. and alt2[i] gt 400. and sha2[i] then pot_bad[i] = 1 endif else pot_bad[i] = 1 endfor if keyword_set(tplot) then store_data,'ion_cutoff_pot_bad',data={x:time2,y:pot_bad} & ylim,'ion_cutoff_pot_bad',-1,2,0 if keyword_set(tplot) then store_data,'cnts_lt_pot',data={x:time2,y:cnts_lt_pot} ;********************************************************** pot3 = interp(pot2>max_nrg*pot_bad,time2,time) pot3a = pot3 pot33 = interp(pot22>max_nrg*pot_bad,time2,time) ;pot3 = (pot3 0. pot33_only = pot33 lt pot3 if keyword_set(tplot) then store_data,'pot33_only',data={x:time,y:pot33_only} if keyword_set(tplot) then ylim,'pot33_only',-1,2,0 pot3 = (pot3 0. pot3_cutoff_cnts = interp(pot2_cutoff_cnts,time2,time) pot3_p_cutoff_cnts = interp(pot2_p_cutoff_cnts,time2,time) ;********************************************************** ; estimate sc potential from swea data - this is turned off since it fails too often if 0 then begin get_data,'mvn_swe_sc_pot',data=tmp9 if size(tmp9,/type) ne 8 then begin mvn_swe_load_l2, /spec,/sumplot mvn_swe_sc_pot,/over,badval=0 get_data,'mvn_swe_sc_pot',data=tmp9 endif get_data,'swe_a4',data=tmp6 npts6 = n_elements(tmp6.x) pk6=fltarr(npts6) & ind6=intarr(npts6) pkh=fltarr(npts6) & pkl=fltarr(npts6) pkc=fltarr(npts6) & pkr=fltarr(npts6) en6=fltarr(npts6) & fx6=fltarr(npts6) mpk = min(abs(tmp6.v-10.),ind0) for i=0l,npts6-1 do begin pk6[i]=max(tmp6.y[i,52:59],ind) pkh[i]=tmp6.y[i,52+ind-4] pkl[i]=tmp6.y[i,52+ind+4] pkr[i]=pkh[i]/pkl[i] pkc[i]=2*pk6[i]/(pkh[i]+pkl[i]) en6[i]=tmp6.v[52+ind] fx6[i]=tmp6.y[i,56] endfor ; end cold = (pkr gt .8) and (en6 lt 7.) and (fx6 gt 6.e7) get_data,'mvn_sta_c6_E',data=tmp99 & time=tmp99.x cold2 = fix(interp((1.*cold),tmp6.x,time)+.6) if keyword_set(tplot) then store_data,'cold',data={x:tmp6.x,y:cold} & ylim,'cold',-1,2,0 if keyword_set(tplot) then store_data,'cold2',data={x:time,y:cold2} & ylim,'cold2',-1,2,0 if keyword_set(tplot) then store_data,'swe_id1',data={x:tmp6.x,y:pkr} if keyword_set(tplot) then store_data,'swe_id2',data={x:tmp6.x,y:pkc} if keyword_set(tplot) then store_data,'swe_id3',data={x:tmp6.x,y:pkh} if keyword_set(tplot) then store_data,'swe_id4',data={x:tmp6.x,y:pkl} if keyword_set(tplot) then store_data,'swe_id5',data={x:tmp6.x,y:en6} if keyword_set(tplot) then store_data,'swe_id6',data={x:tmp6.x,y:fx6} ;tplot,['swe_id1','swe_id5','swe_id6','swe_a4_pot','cold','cold2'] if size(tmp9,/type) eq 8 then begin npts9 = n_elements(tmp9.x) pot9 = tmp9.y ; ind = where(pot9[1:npts9-2] ne 0. or (pot9[0:npts9-3] eq 0 and pot9[2:npts9-1] eq 0)) ; pot9 = pot9[ind+1] ; tim9 = tmp9.x[ind+1] ; pot4 = interp(pot9,tim9,time) ; get rid of single/double bad points ind = where(pot9[2:npts9-3] gt 10. and (pot9[0:npts9-5] eq 0 and pot9[1:npts9-4] eq 0 and pot9[3:npts9-2] eq 0 and pot9[4:npts9-1] eq 0),count) if count ne 0 then pot9[ind+2] = 0 ind = where(pot9[2:npts9-3] gt 10. and (pot9[2:npts9-3] gt 2.*pot9[0:npts9-5] and pot9[2:npts9-3] gt 2.*pot9[1:npts9-4] $ and pot9[2:npts9-3] gt 2.*pot9[3:npts9-2] and pot9[2:npts9-3] gt 2.*pot9[4:npts9-1]),count) if count ne 0 then pot9[ind+2] = 0 ; get rid of single/double bad points ind = where(pot9[1:npts9-2] eq 0. and (pot9[0:npts9-3] ne 0 and pot9[2:npts9-1] ne 0),count) if count ne 0 then pot9[ind+1] = (pot9[ind]+pot9[ind+2])/2. ind = where(pot9[1:npts9-2] eq 0. and pot9[0:npts9-3] ne 0,count) if count ne 0 then pot9[ind+1] = pot9[ind] ind = where(pot9[1:npts9-2] eq 0. and pot9[2:npts9-1] ne 0,count) if count ne 0 then pot9[ind+1] = pot9[ind+2] ; do it again - get rid of up to 4 bad points ind = where(pot9[1:npts9-2] eq 0. and (pot9[0:npts9-3] ne 0 and pot9[2:npts9-1] ne 0),count) if count ne 0 then pot9[ind+1] = (pot9[ind]+pot9[ind+2])/2. ind = where(pot9[1:npts9-2] eq 0. and pot9[0:npts9-3] ne 0,count) if count ne 0 then pot9[ind+1] = pot9[ind] ind = where(pot9[1:npts9-2] eq 0. and pot9[2:npts9-1] ne 0,count) if count ne 0 then pot9[ind+1] = pot9[ind+2] pot4 = interp(pot9,tmp9.x,time) pot4a = interp(tmp9.y,tmp9.x,time) ind = where(tmp9.y gt 0.,count) if count eq 0 then print,'Alert: No sc potential data available from swea' endif else begin pot4=replicate(0.,npts) print,'Alert: No sc potential data available from swea - data missing' endelse if keyword_set(tplot) then store_data,'pot4',data={x:time,y:pot4} if keyword_set(tplot) then options,'pot4',colors=1,thick=2 get_data,'mvn_swe_spec_temp,data=tmp8 if size(tmp8,/type) ne 8 then begin mvn_swe_n1d,/mom,background=1e3,minden=1e-4 endif get_data,'mvn_swe_spec_temp',data=tmp8 if size(tmp8,/type) eq 8 then begin tsmooth2,'mvn_swe_spec_temp',3,newname='mvn_swe_spec_temp_3s' ylim,'mvn_swe_spec_temp_3s',1,50.,1 get_data,'mvn_swe_spec_temp_3s',data=tmp8 temp_e = interp(tmp8.y,tmp8.x,time) get_data,'mvn_swe_spec_dens',data=tmp9 tsmooth2,'mvn_swe_spec_dens',3,newname='mvn_swe_spec_dens_3s' ylim,'mvn_swe_spec_dens_3s',0,50.,0 get_data,'mvn_swe_spec_dens_3s',data=tmp9 dens_e = interp(tmp9.y,tmp9.x,time) ; dens_e = 1.2*dens_e*(temp_e gt 13.) ; dens_e = 1.2*dens_e*(temp_e gt 11.) dens_e = 1.2*dens_e ; print,'dens_e minmax = ',minmax(dens_e) endif else begin print,'Alert: No electron temperature data temp_e = replicate(10.,npts) dens_e = replicate(10.,npts) endelse endif ;********************************************************** ;********************************************************** ; combine potential estimates ; if in sunlight and >250km and densitymax_ec then s/c potential is set to 0 ;sc_neg = den gt max_den or alt lt max_alt or shadow ;sc_neg = alt lt max_alt or shadow ;sc_pos = (alt gt max_alt or ((alt gt 250.) and (den lt max_den) and (p_ec gt max_ec)) or ((alt gt 250.) and (p_ec gt 400.) and (p_low lt 2.)) or ((alt gt 400.) and (c6_mode lt 2))) and (not shadow) if keyword_set(swe) then den2 = (den>swe_den) else den2=den ; sc_pos = (((alt gt 250.) and (den2 lt max_den) and (p_ec gt max_ec))) and (not shadow) ; used to define sc_neg, p_ec is proton ec, low density and high proton ec and not shadow ; sc_pos = ( (alt gt 250.) and (den2 lt max_den) and (p_ec gt max_ec) and (pot3 gt 2.5) ) and (not shadow) ; used to define sc_neg, p_ec is proton ec, low density and high proton ec and not shadow ; sc_pos = ( (alt gt 250.) and (den2 lt max_den) and (pot3 gt 2.5 ) ) and (not shadow) ; used to define sc_neg, p_ec is proton ec, low density and high proton ec and not shadow sc_pos = ( (alt gt 250.) and (den2 lt max_den) and (pot3 gt pot_err ) ) and (not shadow) ; changed 20161205 used to define sc_neg, p_ec is proton ec, low density and high proton ec and not shadow scpot_invalid = (alt gt max_alt) or $ ; max altitude for scpot calculation, generally no limit ; ((alt gt 250.) and (p_ec gt 400.) and (p_low lt 4.)) or $ ; eliminate sw or sheath data, p_ec is proton ec, p_low is proton counts less than 10eV - doesn't work well ; ((alt gt 250.) and (p_ec gt 400.) and (p_low lt 4.) and (o_low lt 10)) or $ ; changed 20161205, eliminate sw or sheath data, p_ec is proton ec, p_low is proton counts less than 10eV - doesn't work well ((alt gt 250.) and (p_ec gt 400.) and (p_den gt 2.) and (p_den gt o_den)) or $ ; changed 20170302, eliminate sw or sheath data, p_ec is proton ec, p_den is proton density, modified for feb27,0150UT ((alt gt 400.) and (p_vel gt 120.) and (not shadow) and (p_den gt o_den)) or $ ; changed 20170302, eliminate sheath, 120 is too high for subsolar point, modified for feb27,0150UT ; ((alt gt 250.) and (p_ec gt 400.) and (p_den gt 2.)) or $ ; changed 20161205, eliminate sw or sheath data, p_ec is proton ec, p_den is proton density - doesn't work that well ; ((alt gt 400.) and (p_vel gt 120.) and (not shadow)) or $ ; eliminate sheath, 120 is too high for subsolar point - may need to modify if additional conditions don't work ; ((alt gt 400.) and (p_ec gt 100.) and (o_ec gt 100.) and (not shadow)) or $ ; added 20161130, eliminate sheath, ; ((alt gt 400.) and (p_ec gt 100.) and (o_ec gt 100.) and (pot3 gt 2.5) and (not shadow)) or $ ; added 20161130, eliminate sheath, ; ((alt gt 400.) and (p_ec gt max_ec) and (o_ec gt max_ec) and (pot3 gt pot_err) and (not shadow)) or $ ; added 20161130, eliminate sheath, ; ((alt gt 400.) and (p_ec gt 100.) and (o_med lt 20.) and (not shadow)) or $ ; added 20161130, eliminate sheath, ; ((alt gt 400.) and (p_ec gt max_ec) and (o_med lt 20.) and (not shadow)) or $ ; added 20161130, eliminate sheath, ((alt gt 400.) and (p_low lt 4.) and (o_low lt (10>((o_cnt+o_med)/40.))) and (not shadow)) or $ ; added 20161130, eliminate counts less than background, ((alt gt alt_vo2) and (pot3 ge max_nrg*.99)) or $ ; assume invalid operation ion cutoff exceeds max_nrg ((alt gt 400.) and (c6_mode lt 2) and (not shadow)) or $ ; assume invalid operation when ram and conic mode are used at high altitude ((den lt .5) and (not shadow)) ; assume it fails if density is too low sc_neg = 1-sc_pos if keyword_set(tplot) then store_data,'scpot_invalid',data={x:time,y:scpot_invalid} if keyword_set(tplot) then ylim,'scpot_invalid',-1,2,0 if keyword_set(tplot) then store_data,'sc_neg',data={x:time,y:sc_neg} if keyword_set(tplot) then ylim,'sc_neg',-1,2,0 if keyword_set(tplot) then store_data,'sc_pos',data={x:time,y:sc_pos} if keyword_set(tplot) then ylim,'sc_pos',-1,2,0 if keyword_set(tplot) then store_data,'alt_gt_250',data={x:time,y:(alt gt 250.)} if keyword_set(tplot) then ylim,'alt_gt_250',-1,2,0 if keyword_set(tplot) then store_data,'den_lt_max',data={x:time,y:(den2 lt max_den)} if keyword_set(tplot) then ylim,'den_lt_max',-1,2,0 if keyword_set(tplot) then store_data,'ec_gt_max',data={x:time,y:(p_ec gt max_ec)} if keyword_set(tplot) then ylim,'ec_gt_max',-1,2,0 ;pot_all = -(pot3 < (pot+1000.*(alt gt alt_vo2)))*(1.*sc_neg) + (pot4 > 1.)*(-sc_neg+1.) pot_all = -(pot3 < (pot+1000.*(alt gt alt_vo2)))*(1.*sc_neg)*(1.*qf_valid)*(1.-scpot_invalid) ; sc_pos comes into play, pot3 is ion cutoff, the alt cutoff at alt_vo2=300km should be variable ion_cutoff_used = round((pot3 lt (pot+1000.*(alt gt alt_vo2)))*(1.*sc_neg)*(1.*qf_valid)*(1.-scpot_invalid)) ; true-false for valid ion cutoff if keyword_set(tplot) then store_data,'ion_cutoff_used',data={x:time,y:ion_cutoff_used} ;tplot,/add,'ion_cutoff_used' ; pot_cutoff_valid is used to force selection of pot3 if pot_cutoff_valid=1 ; this is needed if along track winds impact scpot determined by vo2 ;pot_cutoff_valid = (pot3 lt 1.2*(pot+1000.*(alt gt alt_vo2)))*(1.*sc_neg)*(1.*qf_valid)*(1.-scpot_invalid) and (alt gt ((min_alt+25.)>180.)) ; pre-20170331 algorithm pot_cutoff_valid = (pot3_cutoff_cnts gt 5.1)*(1.*sc_neg)*(1.*qf_valid)*(1.-scpot_invalid) and (pot3_p_cutoff_cnts gt 2.1 or alt gt alt_vo2) and (alt gt (200.<(min_alt+35.)>180.)) ; are 35 and 180 the best parameters? ind = where(not pot_cutoff_valid[1:npts-2] and pot_cutoff_valid[0:npts-3] and pot_cutoff_valid[2:npts-1],count) if count gt 0 then pot_cutoff_valid[ind+1] = 1 ind = where(pot_cutoff_valid[1:npts-2] and not pot_cutoff_valid[0:npts-3] and not pot_cutoff_valid[2:npts-1],count) if count gt 0 then pot_cutoff_valid[ind+1] = 0 ; if keyword_set(cutoff) then pot_all = -(pot3*pot_cutoff_valid + pot*(1-pot_cutoff_valid))*(1.*sc_neg)*(1.*qf_valid)*(1.-scpot_invalid) ; used prior to 20180418 if keyword_set(cutoff) then pot_all = -(pot3*pot_cutoff_valid + ((pot+1000.*(alt gt alt_vo2))