; +
;PROCEDURE:      cal_bs_param,pos,magf,BOW=bow,VEL=vel,VSW=vsw
;PURPOSE:
;       Procedure returns parameters that describe intersection w/ Bow Shock.
;INPUT: 
;       pos:    vector(3),  [Xs,Ys,Zs] s/c position in GSE with units of Re
;       magf:   vector(3),  [Bx,By,Bz] magnetic field vector.
;OPTIONS 
;       BOW:    structure       Structure containing {L,ecc,X0}
;               L:      scalar,         standoff parameter with units of Re
;               ecc:    scalar,         eccentricity of shock
;               X0:     scalar,         focus location in Re 
;       VEL:    scalar,         parallel particle velocity along magnetic field
;       VSW:    vector,         solar wind velocity, same units as VEL
;
;       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
;               x0= 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
;OUTPUTS:    (All through keywords)
;       BSN is the angle between the shock normal and the
;              field line that passes through the s/c.
;       LSN is the distance along the field line to the shock.
;       NSN is the distance from shock nose to field line crossing.
;       SHPOS is the position at the shock-B intersection.
;       SHNORM is the shock normal vector.
;       CONNECT is 0 if unconnected,  1 if connected.
;       STRUCT:  a structure to which all these parameters are added.
;
;NOTES: 
;       see "get_bsn", "add_bsn"
;
;CREATED BY:
;       J.McFadden      95-9-20
;LAST MODIFICATION:  95/10/17
;       95-10-5 McFadden    Added BOW keyword.
;       95-10-0 Larson   Major Mods   changed to procedure. Added shock normal
; -
pro cal_bs_param,pos,magf,bow=bow,vpar=vpar,vsw=vsw,  $
   shpos   = shpos, $
   shnorm  = shnorm, $
   BSN     = bsn,  $
   LSN     = lsn,   $
   NSN     = nsn,   $
   CONNECT = connect,  $
   STRUCT  = struct
   

if data_type(bow) ne 8 then $
    bow={standoff:23.3,eccentricity:1.16,x_offset:3.0}

missing = !values.f_nan
shpos = replicate(missing,3)
shnorm = replicate(missing,3)
bsn = missing
lsh = missing
nsn = missing
connect = missing
 
L=bow.standoff
ecc=bow.eccentricity
x0=bow.x_offset

if not keyword_set(vel) then vel=3.e5
if not keyword_set(vsw) then vsw=[0.,0.,0.]
bfld = magf/sqrt(total(magf^2))    ; normalize
x0=1.d*x0
y0=0.d
z0=0.d
a = L/(ecc^2-1)
b = L/(ecc^2-1)^.5
c = L*ecc/(ecc^2-1)
nsh = [c+x0+a,0,0]

x1=pos(0)-x0-c
y1=pos(1)-y0
z1=pos(2)-z0
a1 = b^2*bfld(0)^2 - a^2*(bfld(1)^2+bfld(2)^2)
b1 = 2*b^2*bfld(0)*x1 - 2*a^2*(bfld(1)*y1+bfld(2)*z1)
c1 = b^2*x1^2 - a^2*(y1^2+z1^2+b^2)

b24ac= b1^2-4*a1*c1
if b24ac ge 0 then begin

   l1 = (-b1+(b24ac)^.5)/(2.*a1)
   l2 = (-b1-(b24ac)^.5)/(2.*a1)

   sh1 = pos + l1 * bfld
   sh2 = pos + l2 * bfld
;   sh1=[pos(0)+l1*bfld(0),pos(1)+l1*bfld(1),pos(2)+l1*bfld(2)]
;   sh2=[pos(0)+l2*bfld(0),pos(1)+l2*bfld(1),pos(2)+l2*bfld(2)]
   connect = 1

   if sh1(0) gt nsh(0) and sh2(0) gt nsh(0) then begin 
;       print,'B field line does not cross shock'
       connect = 0
   endif else if sh1(0) gt nsh(0) then begin
       shpos=sh2
       lsh=l2
   endif else if sh2(0) lt nsh(0) and abs(l1) gt abs(l2) then begin
       shpos=sh2
       lsh=l2
   endif else begin
       shpos=sh1
       lsh=l1
   endelse

 if abs(1-(shpos(0)-x0-c)^2/a^2+(shpos(1)^2+shpos(2)^2)/b^2) gt .0001 then begin
        print,' Shock calculation error!!!'
 endif

; use gradient to dermine shock normal
   shnorm = [2/a^2 , -2/b^2, 2/b^2] * shpos
   shnorm = shnorm/sqrt(total(shnorm^2))        ; normalize it


   bsn= !radeg*acos(total(shnorm*bfld))
endif

add_str_element,struct,'connect',connect
add_str_element,struct,'shpos',shpos
add_str_element,struct,'shnorm',shnorm
add_str_element,struct,'bsn',bsn
add_str_element,struct,'lsn',lsh

return

end