phi Function

  real function phi(rp)
    implicit none
    real, intent (in) :: rp
    integer, parameter :: nimax = 200
    real :: pbar, pb(nimax), dpb, qi, qnext
    integer :: i

    pbar = psifun(rp)
    ! Construct a uniformly space grid from 0 to pbar
    dpb = pbar / float(nimax - 1)
    do i = 1, nimax
       pb(i) = dpb * float(i - 1)
    end do

    ! Integrate area under q(pb) curve from pb = 0 to pbar
    phi = 0. ; qi = geom%qfun(pb(1))
    do i = 1, nimax - 1
       qnext = geom%qfun(pb(i+1))
       phi = phi + 0.5 * dpb * (qnext + qi)
       qi = qnext
    end do
    phi = phi * (rpmax - rpmin)
  end function phi