function near_sph_line_f,t,r=r,v=v,_extra=extra
  return,(xyz_to_lla(r+v*t,_extra=extra))[2];Just return the altitude
end

;Function near_sph_line
;  Calculate the nearest point on the ray R+V*t to a spheroid
;  centered on the origin
;
;  If ray hits surface, then R+V*t is on the surface, dist is zero,
;  lat and lon are coords on ellipse surface
;
;  Since it is a ray, only positive t are considered. If calculated
;  t would be negative, this means that the calculated closest point
;  is behind R. Closest point is then moved to R.
;
;Input
;  r - Origin of the input ray
;  v - direction vector of the input ray, must not be [0,0,0]
;  /neg_dist_ok - If set, finds the point on the ray with the minimum
;                 altitude, even if this altitude is negative
;Input keywords
;  (a,flat, etc) keywords used by xyz_to_lla to define the spheroid
;Output keywords
;  xyz - Point in space on line which is selected
;  lat - Geodetic latitude of selected point, in radians
;  lon - Longitude of selected point, in radians
;  dist - Altitude of selected point, in same length units as a, km by default
;Returns
;  t of selected point
function near_sph_line,r,v,lat=lat,lon=lon,dist=dist,xyz=xyz,neg_dist_ok=neg_dist_ok,_extra=extra
  ;Check for actual spheroid impact
  if keyword_set(neg_dist_ok) then begin
    ;Force this to be a NaN so that the next if works as if there were no itersection
    xyz=[1d,1,1]*!values.d_nan
  end else begin
    ;Check if there really is an intersection
    xyz=raysurfintersect(R,V,0,_extra=extra)
  end
  if finite(xyz[0],/nan) then begin
    ;If xyz is nan, there was no intersection (or neg_dist_ok was set)
    ;if the keyword is set, get the minimum dist even if it is inside the surface
    dist=golden(-2*norm(r),0,2*norm(r),'near_sph_line_f',1d-6,r=r,v=v,xopt=t,_extra=extra)
    xyz=r+v*(t>0)
    lla=xyz_to_lla(xyz,_extra=extra)
    lat=lla[0]
    lon=lla[1]
    if t lt 0 then begin
      dist=lla[2]
    end
    return,(t>0)
  end else begin
    ;The ray does intersect the surface
    junk=max(abs(v),w)   ;This will only fail if all 3 components are 0, in which case
    t=(xyz[w]-r[w])/v[w] ;it will fail long before here
    lla=xyz_to_lla(xyz,_extra=extra)
    lat=lla[0]
    lon=lla[1]
    dist=0
    return,t
  end
end