;FUNCTION GetPixelRayCam
;  Purpose
;    Get the pixel rays for all pixels in a camera, in camera coordinates.
;  Input
;    cam       - Camera number. 0=PX, 1=PY, 2=MX, 3=MY
;    caminfo   - camera info structure as created by LoadCamInfo()
;    /grid     - If not set (default) treats each ray as the center of a pixel. The
;                outside rays will be set in from the edge of the FOV by 1/2 pixel.
;                If set, treats each ray as a node on a grid covering the entire FOV.
;                The outer rays on all four sides will be exactly on the edge of 
;                the FOV.
;  Return
;    a 2D grid of vectors (3D array, [col,row,vec_component]). Each is
;    a normalized vector pointing in the direction a particular pixel looks,
;    in Camera coordinates.
;    Camera frame has axis pointing at +z, increasing col pointing at +y
;    Increasing row pointing at -x
;
;  Distortion model: Each vector is pointed at
;    [kx*tan(thetax)+poly(d),ky*tan(thetay)+poly(d),1]
;   where
;    thetax is angle along horizontal from camera axis to pixel column
;    thetay is angle along vertical   from camera axis to pixel row
;    kx and ky are the appropriate scale factors for each axis
;    d=sqrt(thetax^2+thetay^2) in pixels
;    poly(theta) is some polynomial correction factor
function GetPixelRayCAM,Cam,caminfo,grid=grid
  dist=[0,0,0,0,0] ;[constant,Linear,quadratic,cubic,4th order,etc...]
                   ;This is a null correction - IE camera is perfectly
                   ; represented by tan(theta)
  cols=(*(caminfo.cols))[cam];
  rows=(*(caminfo.rows))[cam];
  pixI_=dindgen(cols)
  pixI_=rebin(pixI_,cols,rows)
  pixJ_=transpose(dindgen(rows))
  pixJ_=rebin(pixJ_,cols,rows)
  PixI=PixI_-(cols/2d)+0.5d
  PixJ=PixJ_-(rows/2d)+0.5d
  if ~keyword_set(grid) then begin
    CamFovPerPixX=(*(caminfo.thfovx))[cam]/((cols)/2d);
    CamFovPerPixY=(*(caminfo.thfovy))[cam]/((rows)/2d);
  end else begin
    CamFovPerPixX=(*(caminfo.thfovx))[cam]/((cols-1)/2d);
    CamFovPerPixY=(*(caminfo.thfovy))[cam]/((rows-1)/2d);
  end
  d=sqrt(pixI^2+pixJ^2)

  result=dblarr(cols,rows,3)
  result[*,*,0]=CamFovPerPixY*PixJ+poly(d,dist)
  result[*,*,1]=CamFovPerPixX*PixI+poly(d,dist)
  result[*,*,2]=1;
  return,normalize_grid(result);

end