function deprob, x, y, a, b, lat_se, lon_se $
               , radians=radians, planetocentric=planetocentric $
               , npangle=npangle, empty=empty
 ;+
 ; NAME:
 ;     deprob
 ;     vvvv           vv 
 ;     DEPRoject onto OBlate spheroid
 ;     ^^^^           ^^             
 ;
 ; PURPOSE:
 ;     Convert x/y coordinates (typically positions on an image) to
 ;     latitude/longitude coordinates on a rotationally symmetrical, 
 ;     spheroidal planet.
 ;     See green notebook pages 160-162, 164 for more details.
 ;     NOTE: everything gets done as double, regardless of input.
 ;     NOTE: x/y's outside of the planet will be set to
 ;           lat/lon/emission angle = -666.
 ;
 ; INPUTS:
 ;     x,y - coordinates of point(s) for which to calculate lat/lon.
 ;           Dimensions of x, y must be identical, units of
 ;           x, y must be same as units of a, b.
 ;     a,b - ellipsoid axes. b is along the rotational axis, a is
 ;           in the equatorial plane. same units as x and y.
 ;     lat_se - latitude of the sub-observer point on the  surface.
 ;           (0 at equator, 90 at north pole, -90 at south pole)
 ;     lon_se - longitude of the sub-observer point on the  surface.
 ;           (0 to 360 west longitude).
 ;
 ; KEYWORDS:
 ;     radians - set this keyword if the angle input/outputs are in radians
 ;           (default is angles input/output in degrees).
 ;     planetocentric - indicate that lat_se is planetocentric (default 
 ;           is planetographic). note that lon_se is not affected because
 ;           the planet is assumed to have rotational symmetry.
 ;     NPangle - Angle of the North Pole measured CCW w.r.t. y-axis
 ;           of the input image coordinate system. (Same convention as 
 ;           JPL horizons, if north is up in the image.) default is 
 ;           npangle=0, or input image already rotated so that the 
 ;           planet's rotational axis is vertical.
 ;     empty - does no math but returns a 1-element structure of-666d
 ;           the right output type.
 ;
 ; OUTPUTS:
 ;     function returns a structure with all vaules set to -666
 ;     if there was an error. otherwise, output is a structure
 ;     ob, where ob.lat is the latitude (planetographic), ob.lon
 ;     is the longitude, ob.lat_pc is the planetocentric
 ;     latitude, and ob.emi is the emission angle. each of these
 ;     has type double and same dimensions as x and y.
 ;
 ; EXAMPLE: obi = deprob(x_km, y_km, a_km, b_km, lat_se, lon_se, npang=17.)
 ;
 ; MODIFICATION HISTORY:
 ;
 ;     revised 2005-11-29 by mikewong: (v2.0) added /empty option, changed output
 ;                                     struct to include emission angle and both
 ;                                     p-centric and p-graphic latitudes, changed
 ;                                     meaning of /planetocentric, fixed some square
 ;                                     roots of negative roundoff errors.
 ;
 ;     written 2005-11-11 by mikewong: (v1.0)
 ;-



  ; VALIDATE AND PARSE INPUT
  ; 
  null_result = {lat:-666d, lon:-666d, emi:-666d, lat_pc:-666d}
  if keyword_set(empty) then return, null_result

  sx = size(x)
  sy = size(y)
  if min(sx eq sy) eq 0 then begin
      print, 'DEPROB: X and Y size mismatch'
      return, null_result
  endif

  if n_elements(NPangle) eq 0 then npangle =0d
  if not keyword_set(radians) then begin ; make sure we're using radians
      se_lat = lat_se  * !dpi / 180
      se_lon = lon_se  * !dpi / 180
      NP_ang = npangle * !dpi / 180
  endif else begin
      se_lat = lat_se + 0d
      se_lon = lon_se + 0d
      NP_ang = npangle + 0d
  endelse

  if keyword_set(planetocentric) then begin ; make sure lat is planetocentric
      se_lat_pg = atan((a/b)^2 * tan(se_lat))
  endif else begin
      se_lat_pg = se_lat
      se_lat    = atan((b/a)^2 * tan(se_lat_pg))
  endelse





  ; CREATE A CLIPPING CIRCLE (R = a) AND ONLY CALCULATE INSIDE IT.
  ; 
  clip_c = where(((x/a)^2 + (y/a)^2) le 1, nclip_c)
  if nclip_c eq 0 then begin
      print, 'DEPROB: no planet !! [nclip_c < 1]'
      return, null_result
  endif




  ; ROTATION AROUND z-AXIS (xhat to yhat)
  ; 
  if Np_Ang ne 0 then begin
      yrot = y[clip_c]*cos(np_ang) - x[clip_c]*sin(np_ang)
      xrot = x[clip_c]*cos(np_ang) + y[clip_c]*sin(np_ang)
  endif else begin
      yrot = y[clip_c] * 1d
      xrot = x[clip_c] * 1d
  endelse






  ; ROTATION AROUND xrot-AXIS (zhat to yrothat)
  ;
  ; yp can be found by solving a quadratic equation. have to also
  ; set the signs of yp and zp because they each have 2 possible solutions.
  if se_lat ne 0 then begin

      ; setup quadratic coeffs to solve for yp
      ; 
      Aq = (cos(se_lat))^2d0 + (a/b)^2 * (sin(se_lat))^2
      Bq = -2d0 * yrot * cos(se_lat)
      Cq = yrot^2d0 + (sin(se_lat))^2 * (xrot^2 - a^2)

      fourAC = 4 * Aq * Cq
      B2     = Bq^2


      ; now that we know yp, we can make a new clipping ellipse and 
      ; drop points outside it
      ; 
      clip_ce = where(fourAC le B2, nclip_ce)
      if nclip_ce eq 0 then begin
          print, 'DEPROB: no planet !! [nclip_ce < 1]'
          return, null_result
      endif
      clip_e = clip_c[clip_ce]


      ; choose proper root for yp
      ; 
      if se_lat ge 0 then $
        yp = (-1 * Bq[clip_ce] + $
              (B2[clip_ce] - fourAC[clip_ce])^0.5) / (2 * Aq) else $
        yp = (-1 * Bq[clip_ce] - $
              (B2[clip_ce] - fourAC[clip_ce])^0.5) / (2 * Aq)
  

      ; use eqn. of ellipsoid to find |zp|, use dyp/dyrot to choose sign of root for zp
      ; 
      dyp_dyrot = (yp * cos(se_lat) - yrot[clip_ce]) / $
        ( ((cos(se_lat))^2 + (a *sin(se_lat)/b)^2) * yp $
          - yrot[clip_ce]*cos(se_lat) )
      zsign = dyp_dyrot / abs(dyp_dyrot)

      root = (1d0 - (yp/b)^2 - (xrot[clip_ce]/a)^2) > 0d
      zp = a * root^0.5  * zsign 

  endif else begin ; for se_lat = 0 special case:

      clip_ce = where(((xrot/a)^2 + (yrot/b)^2) le 1, nclip_ce)
      if nclip_ce eq 0 then begin
          print, 'DEPROB: no planet !! [nclip_ce < 1]'
          return, null_result
      endif
      clip_e = clip_c[clip_ce]
      
      yp = yrot[clip_ce]

      root = (1d0 - (yp/b)^2 - (xrot[clip_ce]/a)^2) > 0d
      zp = a * root^0.5  * zsign 

  endelse







  ; ROTATION AROUND y'-AXIS (z'hat to xrothat)
  ; 
  ypp = yp
  zpp = zp * cos(se_lon)           + xrot[clip_ce] * sin(se_lon)
  xpp = xrot[clip_ce] * cos(se_lon) - zp* sin(se_lon)








  ; CALCULATE LATITUDE AND LONGITUDE AND EMISSION ANGLE
  ; 
  ; longitude
  ; 
  lon = x * 0d0 -666
  lon[clip_e] = (12*!dpi - atan(xpp,zpp)) mod (2*!dpi)


  ; latitude
  ; 
  lat    = x * 0d0 -666
  lat_pc = x * 0d0 -666
  ns = ypp / abs(ypp)

  lat_pc[clip_e] = ns * atan(b / (a* ((B/ypp)^2 -1)^0.5))
  lat[clip_e] = ns * atan(a / (b* ((B/ypp)^2 -1)^0.5))


  ; emission angle (law of cosines)
  ; 
  emi = x * 0d0 -666

  A_angle= abs(lon[clip_e] - se_lon)
  wrap = where(A_angle ge !dpi, n_wrap)
  if n_wrap ge 1 then A_angle[wrap] = (2*!dpi - A_angle[wrap]) mod (2*!dpi)

  b_side = !dpi/2 - se_lat ; use planetocentric to get line of sight
  c_side = !dpi/2 - lat[clip_e] ; use planetographic to get the normal

  emi[clip_e] = acos(cos(b_side)*cos(c_side) + sin(b_side)*sin(c_side)*cos(A_angle))






  ; SET OUTPUT UNITS
  ; 
  if not keyword_set(radians) then begin 
      lat[clip_e]    = lat[clip_e]  * 180 / !dpi
      lon[clip_e]    = lon[clip_e]  * 180 / !dpi
      emi[clip_e]    = emi[clip_e]  * 180 / !dpi
      lat_pc[clip_e] = lat_pc[clip_e]*180 / !dpi
  endif ;else keep in radians





  return, {lat:lat, lon:lon, emi:emi, lat_pc:lat_pc}
end