subroutine integrate_charge(theta1, theta2, q, theta, pt, button)

     use precision_def
     implicit none

     real(rp) q(:), theta(:)
     real(rp) theta1, theta2, button
     real(rp) sum, dTheta, dQp, dT, ds

     integer np
     integer i

     logical pt(:)

     np = size(q)

! print *
     sum=0.
     do i=2,np-1
      if(.not. pt(i-1) .and. pt(i))then
       if(theta(i) < theta1 .or. theta1 < theta(i-1))then
        print *,' problem with boundaries '
        print '(i,3e12.4)',i, theta(i-1), theta1, theta(i)
        stop
       endif

       dTheta = theta(i)-theta(i-1)
       dQp = (Q(i) - Q(i-1))/dTheta
       dT = theta1-theta(i-1)
!       ds = (0.5 * dQp * dT**2 + Q(i-1) * dT)/dTheta + 0.5*Q(i)
       ds = 0.5 * (dQp *(theta(i)**2-2*theta(i)*theta(i-1)-theta1**2+2*theta1*theta(i-1)) &
                  + 2*Q(i-1)*(theta(i)-theta1))/dTheta  + 0.5*Q(i)
       sum = sum + ds
!       print '(i,7e12.4)',i,theta1, theta2, theta(i-1),theta(i), q(i-1), q(i), ds
      endif
      if(pt(i) .and. pt(i-1) .and. pt(i+1))then
        sum = sum + Q(i)
!        print '(i,6e12.4)',i,theta1, theta2, theta(i-1),theta(i), theta(i+1),Q(i)
      endif                     
      if(pt(i) .and. .not. pt(i+1))then
       if(theta(i+1) < theta1 .or. theta2 < theta(i))then
        print *,' problem with boundaries '
        print '(i,3e12.4)',i, theta(i), theta2, theta(i+1)
        stop
       endif
       dTheta = theta(i+1)-theta(i)
       dQp = (Q(i+1) - Q(i))/dTheta
       dT = theta2-theta(i)
       ds =  (0.5 * dQp * dT**2 + Q(i) * dT)/dTheta + 0.5*Q(i)
       ds = 0.5 * (dQp *(theta(i)**2+theta2**2-2*theta2*theta(i)) &
                  + 2*Q(i)*(theta2-theta(i)))/dTheta  + 0.5*Q(i)
       sum = sum + ds
!       print '(i,7e12.4)',i,theta1, theta2, theta(i),theta(i+1),q(i), q(i+1), ds
      endif
     end do
!print *
      button = sum
         

     end