;---------------------------------------
; common parameters
;---------------------------------------

km = 1d5
distobs = 4506172416.d * km

;---------------------------------------
; Set up parameters for ey92 model calculation/test case
;--------------------------------------
table = mrdfits('testphi_v2.fits',1,head, /silent)
ey92_r0 = double( sxpar(head,'R0')) * km
ey92_nu0 = double( sxpar(head,'NU0'))
ey92_lam0 = double( sxpar(head,'LAM0'))
ey92_a = double( sxpar(head,'A'))
ey92_b = double( sxpar(head,'B'))
ey92_r1 = double( sxpar(head,'R1')) * km
ey92_kappa1 = double( sxpar(head,'KAPPA1'))
ey92_htau1 = double( sxpar(head,'HTAU1')) * km
ey92_rsurf = double( sxpar(head,'RSURF')) * km
ey92_order = 4
;delvarx, table, head  ; cleanup

ey92_mu0 = double( sxpar(head,'MU0'))
ey92_Mp = double( sxpar(head,'MP'))
ey92_nustp = double( sxpar(head,'NUSTP'))


;---------------------------------------
; set up the y,phi pairs
;---------------------------------------

;imax = 4001L - 100
imax = 3457
ymid = reverse( (dindgen(imax)+600)*km )

kappa1 = 0.d & farside = 0 & haze = 1
phimid = oclc_ey92_phi_of_rho(ey92_r0, ey92_nu0, ey92_lam0, ey92_a, ey92_b,$
                              ey92_r1, ey92_kappa1, ey92_htau1, distobs, ey92_rsurf,$
                              farside,haze,ey92_order,ymid)
stop
deltaymid = oclc_epq03_deltaymid(ymid)

ib = ( where(phimid lt 0.5) ) [0]

;---------------------------------------
; calculate the model r at bin boundaries
;---------------------------------------

ybot = ymid + deltaymid/2.
ytop = ymid - deltaymid/2.
rbot_model = oclc_ey92_r_of_rho(ey92_r0,ey92_nu0,ey92_lam0,$
                                ey92_a,ey92_b, distobs, ey92_order, ybot)
rtop_model = oclc_ey92_r_of_rho(ey92_r0,ey92_nu0,ey92_lam0,$
                                ey92_a,ey92_b, distobs, ey92_order, ytop)
rmid_model = oclc_ey92_r_of_rho(ey92_r0,ey92_nu0,ey92_lam0,$
                                ey92_a,ey92_b, distobs, ey92_order, ymid)

;---------------------------------------
; make input/output r
;---------------------------------------

rbot = rbot_model
rbot[ib:imax - 1] = 0

oclc_epq03_r,ymid, deltaymid, phimid, ib, imax, rbot

;---------------------------------------
; calculate the theta at bin boundaries
;---------------------------------------

thetabot_model = oclc_ey92_theta(ey92_r0,ey92_nu0,ey92_lam0,$
                                 ey92_a,ey92_b,ey92_order,rbot_model)

thetatop_model = oclc_ey92_theta(ey92_r0,ey92_nu0,ey92_lam0,$
                                 ey92_a,ey92_b,ey92_order,rtop_model)

thetabot = oclc_epq03_theta(ybot,rbot, distobs)

plot, ymid/km, (thetabot - thetabot_model)/thetabot_model,chars=1.4

;---------------------------------------
; calculate the change in theta at bin center
;---------------------------------------

deltathetamid = oclc_epq03_deltathetamid(deltaymid, thetabot)

deltathetamid_model = thetabot_model - thetatop_model

deltarmid = deltaymid - distobs * deltathetamid

deltarmid_model = (rbot_model - rtop_model)

rmid = rbot - deltarmid/2.

rtop = rbot - deltarmid

dthetamid = deltathetamid/deltarmid

thetamid = thetabot - deltathetamid/2.

dthetamid_model = oclc_ey92_dtheta(ey92_r0,ey92_nu0,ey92_lam0,$
                                 ey92_a,ey92_b,ey92_order,rmid_model)

dthetamid_model2 = deltathetamid_model/deltarmid_model

; dthetamid_model differs from dthetamid by 0.01 at 4000 km.

;---------------------------------------
; calculate the refractivity
;---------------------------------------

nubot_model = oclc_ey92_nu(ey92_r0,ey92_nu0,ey92_lam0,$
                           ey92_a,ey92_b, rbot)

;stop

snubot = oclc_epq03_snu(rbot, rtop, deltathetamid, ib)
bnubot = oclc_epq03_bnu(rbot, rtop, thetamid, deltathetamid, ib,0)
nubot = snubot + bnubot
nubot[0:ib-1] = nubot_model[0:ib-1]

;---------------------------------------
; calculate the number density
;---------------------------------------

nbot_model = oclc_ey92_n(ey92_r0,ey92_nu0,ey92_lam0,$
                         ey92_a,ey92_b,ey92_nustp, rbot)

nbot =  oclc_epq03_n(rbot, rtop, thetamid, deltathetamid, $
                     ib, 0, ey92_nustp)
    
;---------------------------------------
; calculate the pressure
;---------------------------------------

pbot_model = oclc_ey92_p(ey92_r0,ey92_nu0,ey92_lam0,$
                         ey92_a,ey92_b,ey92_mp, ey92_mu0, ey92_nustp, rbot)

pbot = oclc_epq03_p(rbot, rtop, thetamid, deltathetamid, $
                         ib, 0, ey92_mp, ey92_mu0, ey92_nustp)

;---------------------------------------
; calculate the temperature
;---------------------------------------

tbot_model = oclc_ey92_t(ey92_r0,ey92_lam0,ey92_b,ey92_mp, ey92_mu0, rbot)

tbot = fltarr(imax)
tbot[ib:imax-1] = pbot[ib:imax-1] / (!phys.k * nbot[ib:imax-1])




end