Subroutine analyze (beam, n_bunch, n_particles)
	!Subroutine takes in beam and analyzes it (rms, mean, deviation...)

	use bmad
        use sim_utils	

	implicit none

	type (beam_struct) :: beam	
	
	integer :: n_bunch, n_particles
	
	integer i, j, k, l
	REAL :: x_avg=0, px_avg=0, y_avg=0, py_avg=0, z_avg=0, pz_avg=0, x_dev=0, px_dev=0, y_dev=0, py_dev=0, z_dev=0, pz_dev=0, x_var=0, px_var=0, y_var=0, py_var=0, z_var=0, pz_var=0
        logical first/.false./

	!REAL, DIMENSION(1:n_particles*n_bunch, 1:6) :: vectors
	
	!do i=1, n_bunch
	!   do j=1, n_particles
	!      vectors((i-1)*n_particles + j) = beam%bunch(i)%particle(j)%vec(1:6)
	!   end do
	!end do

	do i=1, n_bunch
	   do k=1, n_particles
	      x_avg = x_avg + beam%bunch(i)%particle(k)%vec(1)/(n_particles*n_bunch)	      
	      px_avg = px_avg + beam%bunch(i)%particle(k)%vec(2)/(n_particles*n_bunch)	      
	      y_avg = y_avg + beam%bunch(i)%particle(k)%vec(3)/(n_particles*n_bunch)	      
	      py_avg = py_avg + beam%bunch(i)%particle(k)%vec(4)/(n_particles*n_bunch)	      
	      z_avg = z_avg + beam%bunch(i)%particle(k)%vec(5)/(n_particles*n_bunch)	      
	      pz_avg = pz_avg + beam%bunch(i)%particle(k)%vec(6)/(n_particles*n_bunch)
	   end do
	end do

	do j=1, n_bunch
	   do l=1, n_particles
	      x_var = x_var + (beam%bunch(j)%particle(l)%vec(1)-x_avg)**2
	      px_var = px_var + (beam%bunch(j)%particle(l)%vec(2)-px_avg)**2
	      y_var = y_var + (beam%bunch(j)%particle(l)%vec(3)-y_avg)**2
	      py_var = py_var + (beam%bunch(j)%particle(l)%vec(4)-py_avg)**2
	      z_var = z_var + (beam%bunch(j)%particle(l)%vec(5)-z_avg)**2
	      pz_var = pz_var + (beam%bunch(j)%particle(l)%vec(6)-pz_avg)**2
	   end do
	end do

	!write(71, *) x_avg, ' '
	!write(72, *) px_avg, ' '
	!write(73, *) y_avg, ' '
	!write(74, *) py_avg, ' '
	!write(75, *) z_avg, ' '
	!write(76, *) pz_avg, ' '
	!write(81, *) sqrt(x_var/(n_particles*n_bunch)), ' '
	!write(82, *) sqrt(px_var/(n_particles*n_bunch)), ' '
	!write(83, *) sqrt(y_var/(n_particles*n_bunch)), ' '
	!write(84, *) sqrt(py_var/(n_particles*n_bunch)), ' '
	!write(85, *) sqrt(z_var/(n_particles*n_bunch)), ' '
	!write(86, *) sqrt(pz_var/(n_particles*n_bunch)), ' ' 
        if(first)write(66,'(12a12)')'x_avg',  'px_avg',  'y_avg',  'py_avg',  'z_avg',  'pz_avg', &
           'x_rms',  'px_rms','y_rms',  'py_var',  'z_var',  'pz_var'
	write(66,'(12es12.4)' ) x_avg,  px_avg,  y_avg,  py_avg,  z_avg,  pz_avg,  sqrt(x_var/(n_particles*n_bunch)),  sqrt(px_var/(n_particles*n_bunch)), &
            sqrt(y_var/(n_particles*n_bunch)),  sqrt(py_var/(n_particles*n_bunch)),  sqrt(z_var/(n_particles*n_bunch)),  sqrt(pz_var/(n_particles*n_bunch))
		


end subroutine