pro tanhgt,pv,Rs,th,tlat,tp,flat=flat,Req=Req,notan=notan,pp=pp
; given a line of sight unit vector (pv) and 
; the observers position vector (Rs) calculate the
; tangent height (th) and the geodetic latitude tlat and 
; the cartesian position of the taangent point (tp) in ECI. Defaults
; for the equitorial radius and flatness of the earth
; are assumed, but ay be input by the user with the
; Req and flat keywords.
; If th is negative, the pierce point (the point at which the
; line of sight intersects the earth) longitude and latitude will
; be returned in the keyword pp.
; if no tangent point is found, the keyword notan will be set on return

notan=0.
f=(1d)/298.25722356
re=double(6378.1370)
if keyword_set(flat) then f=flat
if keyword_set(req) then re=req
one=1d
e=(one/(one-f)^2)-one
al=[1d,1d,e+1d]

Rsu=Rs/sqrt(total(Rs^2))
cv=-crossp(pv,crossp(Rsu,pv))
cv=cv/sqrt(total(cv^2))

pop=total(al*pv^2)
coc=total(al*cv^2)
RoR=total(al*Rs^2)
poc=total(al*pv*cv)
Roc=total(al*Rs*cv)
Rop=total(al*Rs*pv)

a=pop*poc^2-pop^2*coc
b=pop*(Roc*poc-Rop*coc)
c=(Re^2-RoR)*poc^2+((2d)*Roc*poc-coc*Rop)*Rop

rad0=b^2-a*c
if rad0 lt 0 then begin
	notan=1
	return
endif
L=[(-b+sqrt(rad0))/a,(-b-sqrt(rad0))/a]

ah=[coc,coc]
bh=(Roc+L*poc)
ch=pop*L^2+(2d)*L*Rop+RoR-Re^2
h=[[(-bh+sqrt(bh^2-ah*ch))/ah],[(-bh-sqrt(bh^2-ah*ch))/ah]]
th=min(abs(h),ic)
th=h(ic)
LL=L(0)
if ic mod 2 then LL=L(1)

tp=Rs+LL*pv+th*cv
tpu=tp/sqrt((total(tp*tp)))
glat=90.-acos(tpu(2))*!radeg

if th lt 0 and keyword_set(pp) then begin
   ap=pop
   bp=Rop
   cp=(RoR-Re^2)
   Lp=min([(-bp+sqrt(bp^2-ap*cp))/ap,(-bp-sqrt(bp^2-ap*cp))/ap])
   pp=Rs+Lp*pv
endif
return
end