program initial_distribution_osc

use precision_def
use random_mod
use nr
implicit none
integer, parameter:: n=1000, nsamples=100
real(rp) x(n),xave(nsamples), harvest/1./, y(n)
real(rp) x_average, x_rms
integer i,j,k,l
integer m
integer ind(n)/n*0/, ixx(n)
integer ix, idum
integer sample_size
integer nturns/100/

! initial

do i=1,n

 call ran_gauss(harvest)
 x(i) = harvest

end do

do ix=1,nturns

do i=1,n
   ind(i) = ran(idum)*n*10
end do

call indexx(ind, ixx)

do i=1,n
  y(i) = x(ixx(i))
end do


  x(:) = y(:)

sample_size = n/nsamples
do j=1,nsamples-1
 k=(j-1)*sample_size+1
 l=j*sample_size
 xave(j) = sum(x(k:l))/sample_size
end do

 x_average = sum(x)/n
 x_rms = sqrt(sum((x-x_average)**2)/n)
write(6,'(a1,1x,a,es12.4,a,es12.4)')'#', 'average = ', x_average,' rms = ',x_rms
write(11,*)
write(11,'(a1,1x,a10,a10,2a12)')'#', 'turn','sample','x','x_average'
write(11,*)

do i=1,n
write (11,'(2i10,2es12.4)')i,(i-1)/sample_size+1,x(i),xave((i-1)/sample_size+1)
end do

forall(i=1:n)x(i) = x(i) - xave((i-1)/sample_size+1)

write(13,*)
write(13,'(a1,1x,a10,a10,2a12)')'#', 'turn','sample','x','x_average'
write(13,*)

do j=1,nsamples-1
 k=(j-1)*sample_size+1
 l=j*sample_size
 xave(j) = sum(x(k:l))/sample_size
end do

do i=1,n
write (13,'(2i10,2es12.4)')i,(i-1)/sample_size+1,x(i),xave((i-1)/sample_size+1)
end do

 x_average = sum(x)/n
 x_rms = sqrt(sum((x-x_average)**2)/n)
if(ix == 1)write(14,'(a1,1x,a12,a12)')'#', 'average = ',' rms = '
write(14,'(i10,2es12.4)')ix, x_average,x_rms

end do

end