subroutine int_theta2(lat, ix_start, ix_end, datum_value_a, datum_value_b)
  use bmad
  implicit none
  type (lat_struct) lat
  integer i
  integer ix_start, ix_end
  real(rp) u0, datum_value_a, datum_value_b, f1a, f2a, f1b, f2b
  u0 = 0
  datum_value_a = 0
  datum_value_b = 0
  do i=ix_start+1,ix_end
!    f1a = ( cos(lat%ele(i-1)%a%phi+u0) - lat%ele(i-1)%a%alpha*cos(lat%ele(i-1)%a%phi+u0) )**2/lat%ele(i-1)%a%beta
!    f2a = ( cos(lat%ele(i  )%a%phi+u0) - lat%ele(i  )%a%alpha*cos(lat%ele(i  )%a%phi+u0) )**2/lat%ele(i  )%a%beta
!    f1b = ( cos(lat%ele(i-1)%b%phi+u0) - lat%ele(i-1)%b%alpha*cos(lat%ele(i-1)%b%phi+u0) )**2/lat%ele(i-1)%b%beta
!    f2b = ( cos(lat%ele(i  )%b%phi+u0) - lat%ele(i  )%b%alpha*cos(lat%ele(i  )%b%phi+u0) )**2/lat%ele(i  )%b%beta
     f1a = (lat%ele(i-1)%a%alpha**2+1.)/lat%ele(i-1)%a%beta
     f2a = (lat%ele(i  )%a%alpha**2+1.)/lat%ele(i  )%a%beta
     f1b = (lat%ele(i-1)%b%alpha**2+1.)/lat%ele(i-1)%b%beta
     f2b = (lat%ele(i  )%b%alpha**2+1.)/lat%ele(i  )%b%beta

!    datum_value_a = datum_value_a + lat%ele(i)%value(l$)*(f1a+f2a)/2.0d0
!    datum_value_b = datum_value_b + lat%ele(i)%value(l$)*(f1b+f2b)/2.0d0
     datum_value_a = datum_value_a + lat%ele(i)%value(l$)*(f1a+f2a)/2.0d0
     datum_value_b = datum_value_b + lat%ele(i)%value(l$)*(f1b+f2b)/2.0d0
 enddo
 return
 end