;+
; NAME: 
;  invertlc
; PURPOSE: 
;  Invert an occultation lightcurve
; DESCRIPTION:
;  Invert an occultation lightcurve assuming no ray crossing, spherical 
;  symmetry, large planet, and no far-limb contributions, using the Abel transform.
;
;  Calls atmint for the line-of-sight integral
;  Calls iso to calculate the bending angle at the top of the last bin
;
; CATEGORY:
;  Occultation analysis
; CALLING SEQUENCE:
;  pro invertlc, rho, flux, dist, h, rhoh, theta, r, refr
;
; INPUTS:
;  rho - shadow radius, in increasing order
;  flux   - normalized stellar flux
;  dist - Observer-atmosphere distance
;  h - best guess at the scale height (to compute angles above the last datum)
;  rhoh - best guess at shadow radius of half-light (to compute angles above the last datum)
;
; OPTIONAL INPUT PARAMETERS:
; KEYWORD INPUT PARAMETERS:
; OUTPUTS:
;  theta - bending angle
;  r - planet radius
;  refr - refractivity
;  
; KEYWORD OUTPUT PARAMETERS:
; COMMON BLOCKS:
; SIDE EFFECTS:
; RESTRICTIONS:
;  flux must be positive.
; PROCEDURE:
; MODIFICATION HISTORY:
;  01/11/26 - Written by Leslie A. Young, Southwest Research Institute
;-

pro invertlc, rho, flux, dist, h, rhoh, theta, r, refr

  ; Set up arrays.  The arrays r, theta, and refractivity are indexed
  ; on the bin boundaries of the rho and flux arrays.
  ;
  n = n_elements(rho)
  drho = rho[1]-rho[0]
  theta = fltarr(n+1)
  r = dblarr(n+1)
  refr = fltarr(n+1)
  
  ; get the top theta and r, based on the passed parameters rhoh and H
  fluxN = iso(rho[n-1]+drho/2-rhoh,0,h)
  theta[n] = (H/dist)*(1 - 1/fluxN)
  r[n] = rho[n-1]+drho/2 - dist*theta[n]-rhoh
  
  ; Calculate the rest of the r's and theta's
  for i=n-1,0,-1 do begin &$
    theta[i] = theta[i+1] - (1-flux[i])*drho/dist &$
    r[i] = rho[i]-drho/2-rhoh - dist*theta[i] &$
  end
  r=r+rhoh
  
  ; pad r and theta above the last (highest) datum prior to calling atmint
  rh = rhoh + h
  rpad = r[n] + (findgen(20)+1)*0.5 * H
  thetapad = -(h/dist) * exp(-(rpad-rh)/h)  
  rr = [r,rpad]
  tt = [theta, thetapad]
  
  ; integrate the bending angle to get the refractivity
  for i=n,0,-1 do begin 
    refr[i] = 2*atmint(reverse(rr), -reverse(tt), rr[i], max(rr), 0.) 
  end
  refr = refr/(2.*!pi*rh)
end


;+
; NAME: 
;  test_invertlc
; PURPOSE: 
;  Test the proceedure invertlc vs. the analytic (isothermal) solution
; DESCRIPTION:
;  Calculate an exponential (e.g., isothermal) refractivity, and the
;  lightcurve that would make.  Invert the lightcurve w/ invertlc. 
;  Plot both refractivities (original as a line, inverted as diamonds)
;
; CATEGORY:
;  Occultation analysis
; CALLING SEQUENCE:
;  pro test_invertlc
;
; INPUTS:
; OPTIONAL INPUT PARAMETERS:
; KEYWORD INPUT PARAMETERS:
; OUTPUTS:
; KEYWORD OUTPUT PARAMETERS:
; COMMON BLOCKS:
; SIDE EFFECTS:
;   Resizes window 0 and plots to it.
; RESTRICTIONS:
; PROCEDURE:
; MODIFICATION HISTORY:
;  01/11/26 - Written by Leslie A. Young, Southwest Research Institute
;-
pro test_invertlc

  H=1.
  dist=1
  rhoh=1e4*h-h
  rho = rhoh + findgen(200.)*0.1-15.
  flux = iso(rho,rhoh,h)
  invertlc, rho, flux, dist, h, rhoh, theta, r, refr

  rh = rhoh+h
window,0,xs=700,ys=350
;plot, r-rh, -theta, /ylo,xtit='r',ytit='-theta',/ys
;oplot, r-rh, (h/dist)*exp( -(r-rh)/h), ps=4

plot,  r-rh, refr, /ylo,xtit='r',ytit='refr',/ys
oplot,  r-rh, (H/dist)/sqrt(2.*!pi*rh/H)*exp( -(r-rh)/h), ps=4

end