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

km = 1d5
distobs = 4506172416.d * km

;---------------------------------------
; Set up parameters for ey92 model calculation/test case
;--------------------------------------
table = mrdfits('testphi.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

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

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

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

deltaymid = oclc_epq03_deltaymid(ymid)

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

stop

;---------------------------------------
; 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
;---------------------------------------

snubot = oclc_epq03_snu(rbot, rtop, deltathetamid, ib)
bnubot = oclc_epq03_bnu(rbot, rtop, thetamid, deltathetamid, ib)
snubot1 = oclc_epq03_snu(rbot, rtop, deltathetamid, 1)
bnubot1 = oclc_epq03_bnu(rbot, rtop, thetamid, deltathetamid, 1)
nubot = snubot + bnubot
nubot_model = oclc_ey92_nu(ey92_r0,ey92_nu0,ey92_lam0,$
                           ey92_a,ey92_b, rbot)

nubot_atmint = dblarr(imax)
for i = 0L, imax - 1 do $
  nubot_atmint[i] = atmint(rbot, -thetabot/rbot, rbot[i], rtop[0], 0., /ex)
nubot_atmint = nubot_atmint/!pi

snubot_atmint = dblarr(imax)
for i = ib, imax-1 do $
  snubot_atmint[i] = atmint(rbot, -thetabot/rbot, rbot[i], rtop[ib], 0., /ex)
snubot_atmint = snubot_atmint/!pi

lam = oclc_ey92_lam(ey92_r0,ey92_lam0, ey92_a, ey92_b, rbot_model)
y1 = (sqrt( rbot_model[0]^2 - rbot_model^2 ) / rbot_model) * sqrt(lam/2.)
y = dblarr(imax)
y[ib:imax-1] = ( sqrt( rbot_model[ib]^2 - rbot_model[ib:imax-1]^2 ) / $
                 rbot_model[ib:imax-1]) * $
  sqrt(lam[ib:imax-1]/2.)

plot, rbot/km, nubot_model, /yl, xr=rbot[[imax,ib]]/km
;oplot, rbot/km, snubot, /li
oplot, rbot/km, snubot_atmint, li=2
oplot, rbot/km, nubot_model * erf(y)

plot, rbot/km, (snubot_atmint/erf(y) - nubot_model)/nubot_model, $
  xr=rbot[[imax,ib]]/km

end