;+
; NAME:
;  rd_aat
; PURPOSE: (one line)
;  read the AAT lightcurve
; DESCRIPTION:
;  read the AAT lightcurve
; CATEGORY:
;  
; CALLING SEQUENCE:
;  rd_aat, y, phi, param
; INPUTS:
; OPTIONAL INPUT PARAMETERS:
;   inputdir - directory that contains kao.utts
;   im : return only immersion
;   em : return only emersion
;   sort : sort in order of DECREASING y
;   yrange : return only values in the given yrange
;   f0 : adjustment to eastward position (km). Default = -221.007
;   g0 : adjustment to northwardward position (km). Default = -1110.653
;   upper : upper baseline.  Default = 1.0112830 (from fit_-1100.dat)
;   lower : lower baseline.  Default = 0.0248204 (from fit_-1100.dat)
; 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:
;  
; MODIFICATION HISTORY:
;  Written 2006 Sep 9, Leslie Young
;  2006 Oct 5: change default f0, g0 to -220.980, -1111.47
;  2006 Oct 10: change upper to 1.01127, lower to 0.0211205,
;    
;  2006 Oct 10: change parameters to reflect change in 
;    reference radius from 1250 to 1275 km.
;-

pro rd_aat, y, phi, param, inputdir=inputdir, im=im, em=em, $
            sort=sort, yrange=yrange, distobs=distobs, newp=newp, $
            f0=f0, g0=g0, uppper=upper, lower=lower, $
            hazefrac=hazefrac

    km = 1d5

    distobs = 4506172416.00d * km

    if not keyword_set(inputdir) then inputdir = '../20060916_forefy/INPUT'

    if not keyword_set(f0) then f0 =  -220.980 ; was -221.007
    if not keyword_set(g0) then g0 =  -1111.470 ; -1110.653
    
    fn = inputdir + '/fit_-1100.dat'
    readcol, fn, name, val, err, /skip, for='(A,F,F)', /sil
    if not keyword_set(newp) then begin
        param = double(val[ [1,2,4,5,6,3] ])
    endif else begin
        if n_elements(lower) eq 0 then begin
            param = [2.2174736, 21.962135, 1263.2126, $
                     5.5097857e-10, 15.794196, -0.63491595]
        endif else begin
            param =  [2.11114, 21.4060, 1241.13,  $
                      4.15557e-09, 33.5809, -0.841211]
        endelse
    endelse

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

    fn = inputdir + '/lc_AAT_6.fits'
    table = mrdfits(fn,1,head)
    y=rho_of_fg(table, f0, g0) * km
    if not keyword_set(upper) then upper = 1.01127   ; was val[8]
    if n_elements(lower) eq 0 then lower = 0.0211205 ; was val[7] 
    ; don't use keyword_set to allow lower=0. to stay 0.
    phi = (table.lc - lower) / (upper - lower)
    
    n = n_elements(y)

    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
    
    ; fit haze with specified f0, g0 if we're removing haze
    if n_elements(hazefrac) eq 0 then hazefrac = 0.
    if float(hazefrac) ne 0. then begin
        rm_haze, y, phi, hazefrac
    end
    
                                ; 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(sort) then begin
        indx = reverse(sort(y))
        y = y[indx]
        phi = phi[indx]
    end
    
end