program photon_spline_test

use photon_init_mod
use photon_init_spline_mod
use nr

implicit none

type (photon_init_splines_struct) splines

real(rp) gamma, g_bend, prob, prob2, photon_energy
real(rp) e_factor, E_rel, gamma_phi, phi, phi2, prob_target
integer i, j, n
logical ok
character(100) spline_dir

!

call cesr_getarg(1, spline_dir)
call photon_read_spline (trim(spline_dir), splines)

gamma = 5e9 / m_electron
g_bend = 1.0_rp / 30
e_factor = 3 * h_bar_planck * c_light * gamma**3 * g_bend / 2 

do i = 0, 5
  photon_energy = 10**i
  prob = photon_energy_integ_prob (photon_energy, g_bend, gamma)
  call spline_evaluate (splines.energy_prob, photon_energy, ok, prob2)
  print *, photon_energy, prob, prob2
enddo

!-----------------------------------------------

do i = 1, size(splines%v_angle), 2
  photon_energy = splines%energy_prob(i)%x
  E_rel = photon_energy / e_factor
  n = size(splines%v_angle(i)%prob)
  print *
  do j = 0, 5
    prob_target = j / 5.0_rp
    call photon_vert_angle_init (E_rel, gamma_phi, prob_target)
    phi  = gamma_phi / gamma
    phi2 = zbrent (v_angle_func, splines%v_angle(i)%prob(1)%x, splines%v_angle(i)%prob(n)%x, 1d-6) 
    print '(4es12.4)', photon_energy, prob_target, phi, phi2
  enddo
enddo

!-------------------------------------------------
contains

function v_angle_func(v_angle) result (rel_prob)

real(rp), intent(in) :: v_angle
real(rp) integ_prob, rel_prob

call spline_evaluate (splines%v_angle(i)%prob, v_angle, ok, integ_prob)
rel_prob = integ_prob - prob_target

end function

end program