N = 32;
  theta = 2*pi*(0:N-1)/N;
  x = 0:0.1:10;
  f = cos(sin(theta')*x);
  J0num = mean(f);
  J0ex = besselj(0,x);
  err = J0num- J0ex;
  subplot(2,2,1)
  plot(x,J0num,'x',x,J0ex,'-')
  xlabel('x')
  ylabel('J_0(x)')
  legend('Per Trap N=32','Exact')
  subplot(2,2,2)
  plot(x,err)
  xlabel('x')
  ylabel('N=32 error')