;+
;FUNCTION: cal_bsn2(sx,sy,sz,bx,by,bz,vmag,bow=bow)
;INPUT:
; sx,sy,sz: x, y, and z coodinates of spacecraft in GSE with
; units of Re
; bx,by,bz: magnetic field direction
;OPTIONS:
; bow: structure containing {L,ecc,x0}
; L: standoff parameter with units of Re
; ecc: eccentricity of shock
; x0: focus location in Re
;
; hyperbolic bow shock, see JGR 1981, p.11401, Slavin Fig.7
; r = L/(1+ecc*cos(theta))
; 1 = (x-x0-c)^2/a^2 - (y^2+z^2)/b^2
; default L = b^2/a = 23.5
; ecc = c/a = 1.15
; xoffset= 3
; c = (a^2 + b^2)^.5 = L*e/(e^2-1) = 83.8
; a = L/(e^2-1) = 72.87
; b = L/(e^2-1)^.5 = 41.38
;PURPOSE:
; Function returns a structure that describes the interaction of a B-field
; with a hyperboloid bow shock model.
; Returned structure elements include:
; th_bn: the angle (in degrees) between shock normal and the
; field line that passes through the spacecraft
; l1: the distance along the field line to the shock
; l2: the distance from a point that is missdist from the
; spacecraft in x to the tangent point
; d,m: the distance along x from the spacecraft to a point
; where the B field line would be tangent to the
; bow shock. A postive d means that the field
; line has already intersected the shock. A positive
; m indicates that the field line has not yet
; intersected the shock.
; intpos: position vector in GSE of point where field line
; originating at spacecraft intersects bow shock
; All distances are in Re.
;CREATED BY:
; P. Schroeder
;LAST MODIFICATION:
; L. Winslow 9/15/00
;LAST MODIFICATION: %W% %E%
;-
function cal_bsn2, sx, sy, sz, bx, by, bz, vmag, bow=bow
if keyword_set(bow) then begin
if data_type(bow) ne 8 then $
bow={standoff:23.5,eccentricity:1.15,x_offset:3.0}
endif else begin
bow={standoff:23.5,eccentricity:1.15,x_offset:3.0}
endelse
L=bow.standoff
ecc=bow.eccentricity
xoffset=bow.x_offset
a = L/(ecc^2-1.)
b = L/(ecc^2-1.)^.5
c = L*ecc/(ecc^2-1.)
x0 = sx-xoffset-c
y0 = sy
z0 = sz
aa = b^2*bx^2-a^2*by^2-a^2*bz^2
bb = 2.*b^2*bx*x0-2.*a^2*by*y0-2.*a^2*bz*z0
cc = b^2*x0^2-a^2*y0^2-a^2*z0^2-a^2*b^2
crittest=bb^2-4.*aa*cc
haveint = float(crittest ge 0.)
crittest = haveint*crittest
t1 = (-bb+sqrt(crittest))/2./aa
t2 = (-bb-sqrt(crittest))/2./aa
p1x = x0+t1*bx
p1y = y0+t1*by
p1z = z0+t1*bz
p2x = x0+t2*bx
p2y = y0+t2*by
p2z = z0+t2*bz
p1l0 = float(p1x lt 0.)
p2l0 = float(p2x lt 0.)
t1gt2 = float(abs(t1) gt abs(t2))
havenint = float((p1l0+p2l0) gt 0)
haveint = havenint*haveint
mypx = p1x*p1l0 + p2x*p2l0 - p1x*t1gt2*p1l0*p2l0 - p2x*(not t1gt2)*p1l0*p2l0 + $
1.*(not p1l0)*(not p2l0)
mypy = p1y*p1l0 + p2y*p2l0 - p1y*t1gt2*p1l0*p2l0 - p2y*(not t1gt2)*p1l0*p2l0 + $
1.*(not p1l0)*(not p2l0)
mypz = p1z*p1l0 + p2z*p2l0 - p1z*t1gt2*p1l0*p2l0 - p2z*(not t1gt2)*p1l0*p2l0 + $
1.*(not p1l0)*(not p2l0)
l1pllb = float((bx*(mypx-x0)+by*(mypy-y0)+bz*(mypz-z0)) ge 0.)
magb = sqrt(bx^2+by^2+bz^2)
magdelf = sqrt(4.*mypx^2/a^4+4.*mypy^2/b^4+4.*mypz^2/b^4)
costheta = -(2.*mypx/a^2*bx-2.*mypy/b^2*by-2.*mypz/b^2*bz)/magb/magdelf
th_bn = replicate(!values.f_nan,n_elements(sx))
l1 = replicate(!values.f_nan,n_elements(sx))
intpos = replicate(!values.f_nan,3,n_elements(sx))
normal = intpos
indxint = where(haveint,cnt)
if cnt ne 0 then begin
th_bn(indxint) = acos(costheta(indxint))*180./!pi
l1(indxint) =sqrt((mypx(indxint)-x0(indxint))^2+ $
(mypy(indxint)-y0(indxint))^2+(mypz(indxint)-z0(indxint))^2)
intpos(0,indxint) = mypx(indxint)+xoffset+c
intpos(1,indxint) = mypy(indxint)
intpos(2,indxint) = mypz(indxint)
normal(0,indxint) = 2.*mypx(indxint)/a^2
normal(1,indxint) = 2.*mypy(indxint)/b^2
normal(2,indxint) = 2.*mypz(indxint)/b^2
endif
ajunk = (a*by/b/bx)^2-1.
bjunk = 2.*(a/b/bx)^2*by*bz
cjunk = (a*bz/b/bx)^2-1.
a_t = ajunk*by^2+bjunk*by*bz+cjunk*bz^2
b_t = 2.*ajunk*by*y0+bjunk*(bz*y0+by*z0)+2.*cjunk*bz*z0
c_t = ajunk*y0^2+bjunk*y0*z0+cjunk*z0^2-b^2
tpointcrit = b_t^2-4.*a_t*c_t
havetpoint = float(tpointcrit ge 0.)
tpointcrit = tpointcrit*havetpoint
tt1 = (-b_t+sqrt(tpointcrit))/2./a_t
tt2 = (-b_t-sqrt(tpointcrit))/2./a_t
ytang1 = y0 + tt1*by
ytang2 = y0 + tt2*by
ztang1 = z0 + tt1*bz
ztang2 = z0 + tt2*bz
xtang1 = (a/b)^2/bx*(ytang1*by+ztang1*bz)
xtang2 = (a/b)^2/bx*(ytang2*by+ztang2*bz)
xt1lt0 = float(xtang1 lt 0.)
xtang = xtang1*xt1lt0 + xtang2*(not xt1lt0)
ytang = ytang1*xt1lt0 + ytang2*(not xt1lt0)
ztang = ztang1*xt1lt0 + ztang2*(not xt1lt0)
x0tang = xtang - tt1*bx*xt1lt0 - tt2*bx*(not xt1lt0)
d = replicate(!values.f_nan,n_elements(sx))
l2 = replicate(!values.f_nan,n_elements(sx))
indxtpoint = where(havetpoint,cnt)
if cnt ne 0 then begin
d(indxtpoint) = x0tang(indxtpoint)-x0(indxtpoint)
l2(indxtpoint) = sqrt((x0tang(indxtpoint)-xtang(indxtpoint))^2+ $
(y0(indxtpoint)-ytang(indxtpoint))^2+ $
(z0(indxtpoint)-ztang(indxtpoint))^2)
endif
outstruct = {th_bn: 0., l: 0., lc: !values.f_nan, lm: !values.f_nan,$
d: 0., m: 0., vc: 0.,intpos: replicate(0.,3),normal: replicate(0., 3)}
outarr = replicate(outstruct,n_elements(sx))
if n_elements(sx) eq 1 then begin
outarr.th_bn = th_bn(0)
if l1pllb(0) then begin
outarr.l = l1(0)
l2 = l2(0)
endif else begin
outarr.l = -l1(0)
l2 = -l2(0)
endelse
if d(0) ge 0. then begin
outarr.d = d(0)
outarr.m = !values.f_nan
outarr.lc = l2
endif else begin
outarr.d = !values.f_nan
outarr.m = -d(0)
outarr.lm = l2
endelse
outarr.vc = vmag(0)*outarr.lc/outarr.d
outarr.intpos = intpos
outarr.normal = normal
endif else begin
outarr.th_bn = th_bn
outarr.l = l1*l1pllb - l1*(not l1pllb)
l2 = l2*l1pllb - l2*(not l1pllb)
dge0 = where(d ge 0.,cnt1)
if cnt1 eq 1 then dge0 = dge0(0)
dlt0 = where(d lt 0.,cnt2)
if cnt2 eq 1 then dlt0 = dlt0(0)
outarr.d = d
outarr.m = -d
if cnt2 ne 0 then begin
outarr(dlt0).d = !values.f_nan
outarr(dlt0).lm = l2(dlt0)
endif
if cnt1 ne 0 then begin
outarr(dge0).m = !values.f_nan
outarr(dge0).lc = l2(dge0)
endif
outarr.vc = vmag*outarr.lc/outarr.d
outarr.intpos = intpos
outarr.normal = normal
endelse
return,outarr
end