;+
; NAME:
;  rd_kao
; PURPOSE: (one line)
;  read the KAO lightcurve
; DESCRIPTION:
;  read the KAO lightcurve
; CATEGORY:
;  
; CALLING SEQUENCE:
;  rd_kao, y, phi, param
; INPUTS:
; OPTIONAL INPUT PARAMETERS:
;   im : return only immersion
;   em : return only emersion
;   sort : sort in order of DECREASING y
;   yrange : return only values in the given yrange
;   inputdir - directory that contains kao.utts
; KEYWORD INPUT PARAMETERS:
;  none
; KEYWORD OUTPUT PARAMETERS:
; OUTPUTS:
;  y : shadow radius (aka rho)
;  phi : normalized stellar flux
;  param : [nuref, lamref, r1, kappa1, htau1, b]
; COMMON BLOCKS:
;  None
; SIDE EFFECTS:
; RESTRICTIONS:
;  None
; PROCEDURE:
;  Eliot suggests - use delta rather than dist
; MODIFICATION HISTORY:
;  Written 2006 Sep 9, Leslie Young
;-

pro rd_kao, y, phi, param, inputdir=inputdir, im=im, em=em, sort=sort, yrange=yrange, dist=dist, newp=newp

    km = 1d5

    dist = 4.343d9 * km

    if not keyword_set(newp) then begin
        param = [0.97, 22.4, 1215., 2.22d-8, 29.8, -0.61d]
    endif else begin
        param = [0.93280438, 23.230968, 1214.5133, 2.2716878e-08, 25.699941,-0.82106733]
    endelse

    param[0] = param[0] * exp(-25/55.)  ; move from 1250 to 1275 km ref.
    param[1] = param[1] * (1250/1275.)

    if not keyword_set(inputdir) then inputdir = 'INPUT/'
    
    fn = inputdir + 'kao.utts'
    rdfloat,fn,flux, skip=62
    n = n_elements(flux)
    time = findgen(n) * 0.2 + 0.1 ; seconds after 1988 June 9 10:35:50
    timemid = 96.907            ; s after 1988 June 9 10:35:50
    v = 18.475911d               ; km/s
    rho_min = 865.69            ; km
    y_km = sqrt( rho_min^2 + (time-timemid)^2 * v^2)
    y = y_km * 1d5
    upper = 668.86              ; Pluto + P8 = Pluto * (1 + P8/Pluto)
    lower = upper / (1 + 1/0.24)
    phi = (flux - lower) / (upper - lower)
    
    n = n_elements(y)

                                ; select on immersion or emersion
    imset = keyword_set(im)
    emset = keyword_set(em)
    if imset and emset then begin
        imset = 0b
        emset = 0b
    end
    if imset or emset then begin
        dy = shift(y,-1) - shift(y,1)
        dy[0] = dy[1]
        dy[n-1] = dy[n-2]
        if imset then begin
            indx = where(dy lt 0, n)
            y = y[indx]
            phi = phi[indx]
        endif else begin
            indx = where(dy gt 0, n)
            y = y[indx]
            phi = phi[indx]
        endelse
    endif 
    
    if keyword_set(yrange) then begin
        indx = where(y gt yrange[0] and y lt yrange[1], n)
        y = y[indx]
        phi = phi[indx]
    end
    
    if keyword_set(sort) then begin
        indx = reverse(sort(y))
        y = y[indx]
        phi = phi[indx]
    end
    
end