f_trap Function

  pure real function f_trap(b_mag, theta, jacob)
    implicit none
    real, dimension(-nth:nth), intent (in) :: b_mag, theta, jacob
    real :: ftu, ftl, havg, h2avg, h(-nth:nth)
    ! use expressions from Y. R. Lin-Liu and Bob Miller, coding from Q. P. Liu
    h = min(b_mag / maxval(b_mag),1.0)
    havg = fluxavg(h, theta, jacob)
    h2avg = fluxavg(h**2, theta, jacob)
    ftu = 1.0 - (h2avg / havg**2) * (1.0 - sqrt(1.0 - havg) * (1.0 + 0.5 * havg))
    ftl = 1.0 - h2avg * fluxavg((1.0 - sqrt(1.0 - h) * (1.0 + 0.5 * h)) / h**2, theta, jacob)
    f_trap = 0.75 * ftu + 0.25 * ftl
  end function f_trap