;+
; NAME: 
;  iso
; PURPOSE: 
;  Calculate an isothermal lightcurve
; DESCRIPTION:
;  The isothermal occultation equations are
;
;    theta(r) = - (H/D) exp(-(r-r0)/h)
;    r = rho + D theta(r)
;    rh = rhoh + H
;    phi = 1/(1 + exp(-(r-rh)/H)
;  where D = distance, H = scale height, theta=bending angle,
;  rho=shadow radius, rhoh = half-light shadow radius
;  r = planet radius, rh = helf-light planet radius
;  phi = flux.
;
;  Define x = (r-rh)/H and y = (rho-rhoh)/H, then this becomes the
;  transendental equation
;
;   x = y + exp(-x) - 1
;
;  Recast as f = y + exp(-x) - 1 - x = 0,
;  solve using Newton's method [x = x - f / (df/fx) ]
;
; CATEGORY:
;  Occultation analysis
; CALLING SEQUENCE:
;  phi = iso(rho, rhoh, h)
;
; INPUTS:
;  rho - shadow radius (array)
;  rhoh - shadow radius of half-light
;  h - scale height
;
; OPTIONAL INPUT PARAMETERS:
; KEYWORD INPUT PARAMETERS:
; OUTPUTS:
;  flux   - normalized stellar flux
;  
; KEYWORD OUTPUT PARAMETERS:
; COMMON BLOCKS:
; SIDE EFFECTS:
; RESTRICTIONS:
; PROCEDURE:
; MODIFICATION HISTORY:
;  01/11/26 - Written by Leslie A. Young, Southwest Research Institute
;-
function iso, rho, rho0, h

  y=(rho-rho0)/h
  
  ; set up initial guesses for x
  x=double(y-1)
  g = where(y le -2.1528622055,count)
  if count gt 0 then x(g) = -alog(1-y(g))
  g = where(y gt -2.1528622055 and y le 2, count)
  if count gt 0 then x(g) = y(g)/2
  
  ; iterate to convergence
  for i = 0, 10 do begin 
    phi=1/(1.+exp(-x)) 
    x = x + (y-1.+exp(-x)-x)*phi
  end
  
  return,phi

end