function far_theta,theta,e_pot,alpha
rad = !dpi/180.d
st2 = sin(theta*rad)
sg = (st2 gt 0) * 2 -1
st2 = st2^2
e = sqrt(1+4.d*e_pot*(e_pot-1)*st2)
th0 = acos(-1/e) - acos( (2*e_pot*st2-1)/e )
return,th0/rad*sg
end
function el_response,volts,energy,denergy,nsteps=nsteps
common el_response_com,resp0,volts0,energy0,denergy0,nsteps0
if keyword_set(nsteps0) then begin
if nsteps eq nsteps0 and total(volts ne volts0) eq 0 then begin
energy=energy0
denergy=denergy0
return,resp0
endif
endif
message,/info,'Computing instrument response.'
volts0=volts
w = where(volts lt 0,nw)
if nw ne 0 then begin
volts[w] = max(volts)
page,'Voltage correction'
endif
k_an = 6.42
fwhm = .22
if keyword_set(nsteps) then begin
erange = k_an*minmax(volts)*[1-fwhm,1+fwhm]^2
energy = dgen(nsteps,/log,range=erange)
i = lindgen(nsteps)
denergy = abs(energy[i+1]-energy[i-1])/2
endif
nn = n_elements(energy)
nv = n_elements(volts)
dim = dimen(volts)
es = replicate(1.,nv) # energy
kvs = reform(k_an*volts,nv) # replicate(1.,nn)
sigma = fwhm/(2*sqrt(2*alog(2))) * kvs
resp = exp( -((es-kvs)/sigma)^2 /2) / sqrt(2*!pi) / sigma
resp = resp * (replicate(1.,nv) # denergy)
if ndimen(volts) eq 2 then begin
resp = reform(resp,dim[0],dim[1],nn)
resp = total(resp,1)/dim[0]
endif
resp0=resp
energy0=energy
denergy0=denergy
nsteps0=nsteps
return,resp
end
function el_df,energy,theta,phi,parameters=p
mass=5.6856591e-06
nrg_inf = energy-p.sc_pot
vel = sqrt(2*(nrg_inf > 0.01)/mass)
sphere_to_cart,vel,theta,phi,vx,vy,vz
magf = p.magf
bdir = magf/sqrt(total(magf^2))
bx = bdir[0]
by = bdir[1]
bz = bdir[2]
vsw = p.vsw
nn2 = dimen1(energy)
if nn2 gt 2 then begin
if ndimen(energy) ge 2 then begin
e0 = shift(energy,1,0) & e0[0,*] = e0[1,*]
e2 = shift(energy,-1,0) & e2[nn2-1,*] = e2[nn2-2,*]
endif else begin
e0 = shift(energy,1) & e0[0] = e0[1]
e2 = shift(energy,-1) & e2[nn2-1] = e2[nn2-2]
endelse
de = abs(e2-e0)
wght = 0. > (energy+e0-2.*p.sc_pot)/de < 1.
endif else wght = energy gt p.sc_pot
f = 0.
a = 2./mass^2/1e5
photo = p.ph
if keyword_set(photo.eflux[0]) then begin
eflx2 = spl_init(alog(photo.nrg),alog(photo.eflux),/double)
fphoto = exp(spl_interp(alog(photo.nrg),alog(photo.eflux),eflx2,alog(energy))) /energy^2 /a
f = f+ (1.-wght) * fphoto
endif
if p.core.n ne 0 then begin
vcore = vsw + bdir * p.core.v
vsx = vx - vcore[0]
vsy = vy - vcore[1]
vsz = vz - vcore[2]
vtot2 = vsx*vsx + vsy*vsy + vsz*vsz
cos2a = (vsx*bx + vsy*by + vsz*bz)^2/vtot2
e = exp( -0.5*mass*vtot2 / p.core.t * (1.d + cos2a*p.core.tdif) )
k = (mass/2/!dpi)^1.5 * 1e10
fcore = (k * p.core.n * sqrt((1+p.core.tdif)/p.core.t^3) ) * e
f= f+ wght * fcore
endif
if p.halo.n ne 0 then begin
vhalo = vsw + bdir *p.halo.v
vsx = vx - vhalo(0)
vsy = vy - vhalo(1)
vsz = vz - vhalo(2)
vtot2 = vsx*vsx + vsy*vsy + vsz*vsz
cos2a = (vsx*bx + vsy*by + vsz*bz)^2/vtot2
vh2 = (p.halo.k-1.5)*p.halo.vth^2
kc = (!dpi*vh2)^(-1.5) * gamma(p.halo.k+1)/gamma(p.halo.k-.5) *1e10
fhalo = p.halo.n*kc*(1+(vtot2/vh2))^(-p.halo.k-1)
f = f+ wght*fhalo
endif
return,f
end
function st_swea_distfit, x, $
set=set, $
type = type, $
parameters=p
if not keyword_set(p) then begin
pnrg = [5.18234, 5.92896, 7.25627, 9.41315, 12.8144, 18.2895, 27.2489, 41.7663, 65.2431, 103.320, 164.957, 264.837, 426.769, 689.161, 1112.99]
eflux = [3.39420e+08, 2.77170e+08, 2.05402e+08, 1.49702e+08, 9.07756e+07, 3.88106e+07, 2.41517e+07, 1.01327e+07, 2.40907e+06, 400220., 368695., 271865., 145748., 59322.0, 26367.6]
photo = {eflux:eflux,nrg:pnrg }
core = {n:10.d, t:10.d, tdif:0.d, v:0.d}
halo = {n:.1d, vth:5000.d, k:4.3d, v:0.d}
p13 = {ph:photo,core:core, halo:halo, $
e_shift:0.d, v_shift:0.31d, sc_pot:5.d, vsw:[500.d,0.d,0.d], $
magf:[0.,0.,1.], $
ntot:10.d, $
dflag:0, $
deadtime:0.6d-6,expand:8 }
p=p13
endif
if p.dflag then p.core.n = p.ntot-p.halo.n else p.ntot=p.core.n+p.halo.n
if size(/type,x) ne 8 then begin
if 0 then begin
p.sc_pot = 1. > p.sc_pot < 50.
p.ntot = .05 > p.ntot < 200
p.core.n = .05 > p.core.n < 200
p.core.t = 0.5 > p.core.t < 50
p.core.tdif= -.6 > p.core.tdif < .3
p.core.v = -500 > p.core.v < 500
p.halo.n = 0.00 > p.halo.n < 1
p.halo.vth = 2000 > abs(p.halo.vth) < 6000
p.halo.v = -5000 > p.halo.v < 5000
p.halo.k = 3 > p.halo.k < 5
p.vsw = [-900.,-200.,-200.] > p.vsw < [200.,200.,200.]
endif
return,0
endif
energy = x.energy
theta = x.theta
phi = x.phi
if keyword_set(p.expand) then expand=16
str_element,x,'volts',volts
if not keyword_set(volts) then expand=0
if keyword_set(expand) then begin
nn = x.nenergy
resp = el_response(x.volts+p.v_shift,nrg2,nsteps=expand*nn)
nb = x.nbins
i = replicate(1.,expand)
nn2 = nn*expand
energy = nrg2[*] # replicate(1.,nb)
phi = reform(i # phi[*],nn2,nb)
theta = reform(i # theta[*],nn2,nb)
esteps = resp # nrg2
endif
mass = x.mass
a = 2./mass^2/1e5
f= el_df(energy,theta,phi,param=p)
eflux = f* energy^2 * a
if keyword_set(expand) then begin
eflux = resp # eflux
energy = resp # energy
theta = resp # theta
phi = resp # phi
endif
units = x.units_name
case strlowcase(units) of
'df' : data = eflux/energy^2/a
'flux' : data = eflux/energy
'eflux' : data = eflux
else : begin
crate = x.geomfactor *x.gf * eflux
anode = byte((90 - x.theta)/22.5)
deadtime = (p.deadtime/[1.,1.,2.,4.,4.,2.,1.,1.])(anode)
rate = crate/(1+ deadtime *crate)
bkgrate = 0
str_element,p,'bkgrate',bkgrate
rate = rate + bkgrate
case strlowcase(units) of
'crate' : data = crate
'rate' : data = rate
'counts' : data = rate * x.dt
endcase
end
endcase
if keyword_set(set) then begin
x.data = data
x.energy = energy
str_element,/add,x,'sc_pot',p.sc_pot
str_element,/add,x,'deadtime', deadtime
endif
str_element,x,'bins',value = bins
if n_elements(bins) gt 0 then begin
ind = where(bins)
data = data(ind)
endif else data = reform(data,n_elements(data),/overwrite)
if keyword_set(set) and keyword_set(bins) then begin
w = where(bins eq 0,c)
if (c ne 0) and (set eq 2) then x.data(w) = !values.f_nan
endif
return,data
end