; docformat = 'rst'
;+
;Returns information on the intersection between a ray and a spheroid
;  - The ellipsoid is centered at the origin, and has its polar axis along
;    the Z axis.
;  - Spheroid may be either oblate, prolate, or perfectly spherical.
;  - Rays are defined by vector parametric equation R=R0+V*t.
;  
;Valid intersections always have a positive t parameter. A ray may fail to 
;have a valid intersection due to either:
;  - The full line the ray is part of never intersects the spheroid
;  - The full line does intersect the spheroid either once or twice, but all
;    such intersections are on the wrong half of the line (negative t)
;  
; :Author:
;   Chris Jeppesen : 07-15-2005
;   
; :Examples:
;   result = raysurfintersect(ray_origin, ray_vec, altitude)
;-
;+
; :Params:
;   R0 : in, required, type=vector or array of vectors
;     3-element vector describing the origin of the ray, in distance units
;     matching those of the reference spheroid (km by default). May be a 
;     grid of vectors.
;   V : in, required, type=vector or array of vectors
;     3-element vector describing the direction of the ray. May be a grid 
;     of vectors.
;     NOTE! R0 and V grids must be compatible in the IDL sense. One of 
;     these must be true:
;     1) Both R0 and V are single vectors
;     2) One of R0 or V are single vectors, the other is a grid
;     3) R0 and V0 are grids of the same shape and size
;   Alt : in, required, type=double or array of double
;     Altitude of surface of interest above reference spheroid, in 
;     distance units matching those of the reference spheroid (km by 
;     default). May be a grid of scalars compatible in the IDL sense with
;     R0 and V.
;     
; :Keywords:
;   a : in, optional, type=double
;     Reference spheroid equatorial radius. Defaults to WGS-84 spheroid
;     value in km. Must be scalar.
;   Re : in, optional, type=double
;     Same as a, above. If both are set, this one takes precedence.
;   flat : in, optional, type=double
;     Spheroid flatness = 1-(Rpolar/Requatorial). Defaults to WGS-84
;     ellipsoid value. Positive for oblate spheroids, zero for perfect
;     spheres, negative for prolate spheroids. Must be scalar.
;   Rp : in, optional, type=double
;     Reference spheroid polar radius. If not set, it is calculated from
;     Re and flat above. Must be scalar.   
;   t : out, optional, type=boolean
;     The t parameter for a valid intersection (see above) or NaN if there
;     isn't one. If either R0 or V are grids, this is an array of the same
;     shape and size. 
;     
; :Returns:
;   A 3-element vector of closest valid intersection to R0, or 
;   [NaN,NaN,NaN] if there is no valid intersection. If either R0 or V are 
;   grids, return value is a grid. If some elements of the grid have valid 
;   intersections and some do not, result will have NaN only at grid cells 
;   where there is no valid intersection.
;-
function raysurfintersect,R0,V,Alt,a=a,Re=Re,flat=flat,Rp=Rp,t=t
  ;Default reference spheroid is WGS-84
  if n_elements(a)    eq 0 then a=get_wgs84_const(/re);
  if n_elements(Re)   eq 0 then Re=a;
  if n_elements(flat) eq 0 then flat=get_wgs84_const(/f);
  if n_elements(Rp)   eq 0 then Rp=Re*(1-flat);

  ;Calculate size of spheroid of interest
  cloud_re=Re+Alt;
  cloud_rp=Rp+Alt;

  ;Transform the input vectors such that
  ;this becomes a ray-sphere intersection problem.
  resolve_grid,R0,x=x0,y=y0,z=z0
  x0n=x0/cloud_re;
  y0n=y0/cloud_re;
  z0n=z0/cloud_rp;

  resolve_grid,V,x=xv,y=yv,z=zv
  xvn=xv/cloud_re
  yvn=yv/cloud_re
  zvn=zv/cloud_rp

  ;Quadratic equation coefficients for quadratic equation
  ;describing intersection point. A*t^2+b*t+C=0 at intersection points
  A=xvn^2+yvn^2+zvn^2;
  B=2*(x0n*xvn+y0n*yvn+z0n*zvn);
  C=x0n^2+y0n^2+z0n^2-1;

  ;Quadratic formula discriminant
  D=B^2-4*A*C;
  ;Set the size of the t arrays
  t=xvn*x0n*0;
  t1=t;
  t2=t;

  ;Separate the valids from the invalids
  HasIntersect=D ge 0;
  wHasIntersect=where(HasIntersect,nHasIntersect,complement=wNoIntersect,ncomplement=nNoIntersect);
  if(nHasIntersect ne 0) then begin
    ;For these, the ray intersects the surface, calculate t at the point(s) of intersection
    ;using the rest of the quadratic formula
    t1[wHasIntersect]=(-B[wHasIntersect]+sqrt(D[wHasIntersect]))/(2*A[wHasIntersect]);
    t2[wHasIntersect]=(-B[wHasIntersect]-sqrt(D[wHasIntersect]))/(2*A[wHasIntersect]);

    ;Find the places where t1 is the smaller of the positive values
    T1IsIt=HasIntersect and (((t1 gt 0) and (t1 lt t2)) or ((t1 gt 0) and (t2 lt 0)))
    wT1IsIt=where(T1IsIt,nT1IsIt,complement=wNotT1IsIt,ncomplement=nnotT1IsIt);
    if(nT1IsIt ne 0) then begin
      t[wT1IsIt]=t1[wT1IsIt];
    end

    ;For these, either t1 is negative or t2 is less than t1. Check if t2 is positive
    if (nnotT1IsIt ne 0) then begin
      ;Find the places where t2 is the smaller, but positive
      notT2IsIt=T1IsIt or (t2 le 0)
      wnotT2IsIt=where(notT2IsIt,nnotT2IsIt,complement=wT2IsIt,ncomplement=nT2IsIt);
      if(nT2IsIt ne 0) then begin
        t[wT2IsIt]=t2[wT2IsIt];
      end
      ;For these, t2 is negative also. Both solutions are on the wrong side of the line
      if (nnotT2IsIt ne 0) then begin
        t[wnotT2IsIt]=!values.F_NaN;
      end
    end
  end
  ;For these, the ray never intersects the surface
  if(nNoIntersect ne 0) then begin
    t[wNoIntersect]=!values.F_NaN;
  end

  ;Now we have the correct t parameter(s), calculate the intersection point(s).
  rxn=x0n+t*xvn;
  ryn=y0n+t*yvn;
  rzn=z0n+t*zvn;

  ;Reverse the transformation from above, and compose answer into a grid of vectors
  return,compose_grid(rxn*cloud_re,ryn*cloud_re,rzn*cloud_rp)

end