;Function LLA_to_XYZ
;  Convert a given location relative to a spheroid in lat,lon,and alt,
;  into an equivalent cartesian vector
;
;Input
;  LLA - A grid of vectors. First element is geodetic latitude in radians,
;        second is longitude in radians, third is altitude in same units
;        as ellipsoid a (km by default)
;Keyword input
;  a   - Equatorial radius of ellipsoid. Defaults to WGS-84 parameter, in km
;  flat- Ellipsoid flatness ratio 1-(b/a). Defaults to WGS-84 parameter
;  sid - Greenwich Mean Siderial time. If set appropriately, returned vector
;        will be in ECI rather than ECEF. Defaults to zero, giving ECEF
;Returns
;  A grid of xyz vectors, in the same length unit as ellipsoid a (km by default)
function lla_to_xyz,lla, a=a, flat=flat, sid=sid

  ;Default ellipsoid if none is specified
  if n_elements(a)    eq 0 then a=get_wgs84_const(/a);
  if n_elements(flat) eq 0 then flat=get_wgs84_const(/f);
  if n_elements(sid)  eq 0 then sid=0.0d;

  resolve_grid,lla,x=lat,y=lon,z=alt

  ;We really only need a and e^2. Do it this way, rather than
  ;use get_wgs84_const(/e2), so that a custom ellipsoid works also
  b = a*(1.0d - flat)           ;semi minor axis
  e1sqr = (a*a - b*b) / (a*a)   ;first  numerical eccentricity
  N = a / sqrt(1.0d - e1sqr * sin(lat)^2d);
  x = (N                  + alt) * cos(lat) * cos(lon);
  y = (N                  + alt) * cos(lat) * sin(lon);
  z = (N * (1.0d - e1sqr) + alt) * sin(lat);

  return,compose_grid(x,y,z);

end