;+
; :Params:
;   Cam_name : in, required, type=string
;     Camera name, must be a scalar. Examples: "px" - +X (Solar), "mx" -
;     -X (Antisolar), "py" - +Y (Nadir), "my" - -Y (Nadir)
;   time : in, required, type=double
;     GPS second of the image, must be scalar.
;   alt : in, required, type=double
;     Altitude of surface of interest, in km, must be scalar. Can be
;     anything, but typical values are 0=surface of the earth, 83=cloud
;     deck. May even be above camera viewpoint.
;   caminfo : in, required, type=caminfo structure
;     Structure describing cameras. Returned by loadcaminfo().
;   quat : in, required, type=array of double
;     Quaternion of the spacecraft of the form [3D vector, scalar].
;
; :Keywords:
;   view : in, optional, type=boolean
;     Set to calculate viewing geometry as well as geolocation.
;   common_volume : in, optional, type=boolean
;     A reference to a CSG object. If present, and /view is set, those
;     pixels whose rays' intersection with the cloud deck is inside this
;     CSG will be marked. These are the common volume "red pixels".
;   ray_peak : in, optional, type=boolean
;     Do latitdue and longitude at Rayleigh peak altitude.
;   _extra : in, optional, type=boolean
;     Passed into xyz2llagrid.
;
; :Returns:
;  An anonymous structure, many of whose fields are 2D arrays. For these,
;  the two dimensions correspond to pixel row and column number.
;    lat - 2D float array of geodetic latitude of pixel, rad. Ranges from
;    -pi/2 to pi/2 (90degS to 90degN).
;    lon - 2D float array of longitude of pixel, rad. Ranges from -pi to
;    pi (180degW to 180degE).
;
;  If viewing geometry is calculated, the following fields are also
;  present:
;    sat_lat - scalar float, spacecraft geodetic latitude, rad.
;    sat_lon - scalar float, spacecraft longitude, rad.
;    sat_alt - scalar float, spacecraft altitude above ellipsoid, km.
;    track_azimuth - scalar float, spacecraft track azimuth, rad, measured
;    from true north towards the east.
;    ssa - 2D float array, solar scattering angle, supplement of angle
;    from Sun to pixel geolocation to spacecraft, rad.
;    sva - 2D float array, spacecraft view angle, angle from center of
;    Earth to spacecraft to pixel geolocation, in radians.
;    sza - 2D float array, solar zenith angle, angle from Sun to center
;    of Earth to pixel geolocation, in radians.
;    lit - 2D byte array, lighting flag - 1 indicates that this pixel is
;    in sunlight at the cloud deck, 0 indicates shade. Shading is
;    calculated using the earth surface ellipsoid expanded by the ozone
;    altitude.
;    cv - 2D byte array - 2 if the pixel intersection with the cloud deck
;    is in the common volume "red pixel". Zero otherwise. May eventually
;    take on other values, like 1 for a pixel which views the common
;    volume not at the cloud deck "blue pixel". If common volume is not
;    calculated, this field is a grid of zeros of the appropriate size.
;-
function cam_to_lla, Cam_name, time, alt, caminfo, quat,view=view, $
  common_volume=common_volume, ray_peak=ray_peak, _extra=extra

  if n_elements(time) gt 1 then message,"Time must be a scalar!"

  thiscaminfo=getcaminfo(caminfo,Cam_name)
  ;Get a quaternion which transforms from the camera optical axis frame to the
  ;spacecraft body frame
  qcam2sc=cam_to_sc(thiscaminfo);

  ;Read telemetry to get a quaternion which transforms from the spacecraft body
  ;frame to one parallel to ECEF, centered at the spacecraft
  qsc2global=sc2global(time,quat,/ecef, quat_time)
  ; *****  NOTE THE EARLY RETURN UNDER ERROR CONDITIONS *****
  if n_elements(qsc2global) eq 1 && qsc2global eq -1 then return, -1
  spacecraft_z=transformgrid(qsc2global,[0,0,1])

  ;Combine these two to get a quaternion which transforms from the camera optical
  ;axis frame to one parallel to ECEF
  qcam2global=combine_rotation(qcam2sc,qsc2global);
  cam_ref_x=transformgrid(qcam2global,thiscaminfo.cam_ref_x);
  ; Read telemetry to get the spacecraft state vector (position part) for the correct
  ;   moment in time
  RV=scpv(time, drv=drv, /ecef);
  R=RV[0:2]
  dr=drv[0:2]
  V=RV[3:5]
  nadir=-R/norm(R)
  alongtrack=V/norm(V)
  crosstrack=crossp_grid(alongtrack,nadir)
  crosstrack/=norm(crosstrack)
  alongtrack_err=abs(comp(dr,alongtrack))
  crosstrack_err=abs(comp(dr,crosstrack))
  radial_err=abs(comp(dr,nadir))
  spacecraft_nadir_deflection=vangle(nadir,spacecraft_z)*!radeg

  ;Transform the camera pixel vector grid as a group from camera frame to frame parallel
  ;to ECEF. These become the V vector for the rays we will use.
  VGrid=transformgrid(qcam2global,thiscaminfo.pixelvec);

  ;Find the points of intersection between the rays and the cloud deck. The rays
  ;are in the form r(t)=R+VGrid*t, where R is constant for all rays (it's the spacecraft
  ;location) and VGrid is different for each pixel in the image. This returns a grid
  ;of cartesian vectors the same shape as the VGrid, ie one xyz vector per pixel
  xyz=raysurfintersect(R,VGrid,alt);

  ;Convert all the cartesian coordinates to lat, lon, alt.
  lla=xyz_to_lla(xyz);

  ; If there is no need for view geometry
  if ~keyword_set(view) then return, {lat: lla[*,*,0], lon: lla[*,*,1] }

  ;Calculate spacecraft lat, lon, alt
  satlla=xyz_to_lla(R);
  sat_lat=satlla[0];
  sat_lon=satlla[1];
  sat_alt=satlla[2];

  viewGeometry,r,xyz,time,lla=lla,sva=sva,sza=sza,sun=sun,cva=cva,_extra=extra
  ;Unit vector from center of Earth towards Sun
  sun=sunvec(time,/ecef);
  ;calculate solar scattering angle
  ssa=vangle(sun ,VGrid);

  ;calculate polarization angle
  if n_elements(cam_ref_x) gt 0 then begin
    n_sca=crossp_grid(VGrid,sun)
    n_sca=normalize_grid(n_sca)
    n_ref=crossp_grid(cam_ref_x,VGrid)
    n_ref=normalize_grid(n_ref)
    pol=vangle(n_sca,n_ref)
  endif else pol = []

  ;Calculate rayleigh peak geometry
  ray_peak_alt=sza_alt(sza*!radeg)   ;Altitude for rayleigh_peak from sza
  xyz_ray_peak=raysurfintersect(R,vgrid,ray_peak_alt)
  lla_ray_peak=xyz_to_lla(xyz_ray_peak);
  viewGeometry,r,xyz_ray_peak,time,sza=sza_ray_peak,cva=cva_ray_peak,_extra=extra


  ; If there is a valid common_volume check to see if it intersects the image.
  ; Otherwise return a correctly dimensioned array of zeros.
  if keyword_set(common_volume) then begin
    ; Calculate common volume red pixels. Those intersection points which are
    ; inside the common volume CSG are red pixels.
    cv_red=common_volume->is_inside(xyz)
    ; The multiplication by 2 is a legacy feature for which there is at present
    ;   (2016/09/22) no documented reason.
    cv=byte(cv_red*2)
  endif else begin
    cv=bytarr(shape_grid(lla))
  endelse

  if keyword_set(ray_peak) then lla=lla_ray_peak
  
  return,{lat: lla[*,*,0], $
    lon: lla[*,*,1], $
    vgrid: vgrid, $
    sat_lat: sat_lat,    $
    sat_lon: sat_lon,    $
    sat_alt: sat_alt,    $
    track_azimuth:track_azimuth(pv=RV), $
    alongtrack_error:alongtrack_err, $
    crosstrack_error:crosstrack_err, $
    radial_error:radial_err, $
    ssa: ssa,        $
    sza: sza,        $
    sza_ray_peak:sza_ray_peak,     $
    sva: sva,        $
    cva: cva,        $
    cva_ray_peak:cva_ray_peak,      $
    pol: pol,        $
    quat: qcam2global,   $
    quat_time : quat_time, $
    cv:  cv              $
  };
end