module modulo1
implicit none
private
public :: interazione,initpot,kr
!public :: fdft 
integer,parameter::kr=selected_real_kind(12)
contains

subroutine initpot(vext, nx, ny, nz, g)
integer, intent(in) :: nx, ny, nz
integer :: i,j,k
real(kind=kr), intent(in) ::g 
real(kind=kr),intent(out),dimension(nx,ny,nz) :: vext
vext=0.
do i=1,nz
  vext(:,:,i)=g*i
end do
end subroutine initpot

subroutine interazione(pos,nbody,f,upot,vext,rho,nx,ny,nz,dx,temp,fakemass,nrho,ds)
real(kind=kr), intent(in), dimension(:,:) :: pos
real(kind=kr), intent(inout), dimension(nx,ny,nz) :: rho
real(kind=kr),dimension(nx,ny,nz) :: force 
real(kind=kr), intent(in) :: temp, dx
integer, intent(in) :: nbody, nx,ny,nz,nrho
real(kind=kr), intent(in) :: fakemass,ds
real(kind=kr), intent(in), dimension(nx,ny,nz) :: vext
real(kind=kr), intent(out) :: upot
real(kind=kr) :: free
real(kind=kr), intent(out), dimension(:,:) :: f
real(kind=kr), dimension(size(pos,1)) :: posij
real(kind=kr) :: rij
integer :: i,j

force = -temp*log(rho)-vext
do i=1,nrho
  rho = rho + force/fakemass * ds**2
  force = -temp*log(rho)-vext
  free = temp*sum(rho*(log(rho)-1))+sum(rho*vext)
  free = free*dx
write(13,*) i, free
end do
upot = 0
f = 0
do i=1,nbody
do j=1,nbody
  if( i==j ) cycle
  posij = pos(:,i)-pos(:,j)
  rij=sqrt( dot_product(posij,posij) )
  upot = upot + rij**(-1) 
  f(:,i) = f(:,i) + rij**(-3)*posij!-fdft(i,pos,rho,ngrid,dx)
end do
end do
upot = upot/2 + free

end subroutine interazione

end module modulo1

program corpi3d
use modulo1, rk=>kr
implicit none
integer :: nstep,it,nbody,nsave,ngrid,ios,i,j, jz
real(kind=rk) :: dt,mepot,mekin,massa=1.,alfa=1.,vmax
real(kind=rk),dimension(:,:),allocatable :: pos
real(kind=rk),dimension(:),allocatable :: velcm
real(kind=rk),dimension(:,:),allocatable :: ekin,vel,f
logical :: dyn,anneal,cont
integer, parameter ::  nx=10,ny=10,nz=10,nrho=10000
real(kind=rk) :: temp=1, dx=1., fakemass=1, ds=0.001, g=1.
real(kind=rk), dimension(nx,ny,nz) :: rho=1./(nx*ny*nz), vext

write(unit=*,fmt="(a)",advance="no")"dynamics?"
read*, dyn
if (dyn) then
  write(unit=*,fmt="(a)",advance="no")"annealing?"
  read*, anneal
  if (anneal) then
    write(unit=*,fmt="(a)",advance="no")"scaling parameter?"
    read*, alfa
  end if
end if
write(unit=*,fmt="(a)",advance="no")"n of bodies?"
read*, nbody
write(unit=*,fmt="(a)",advance="no")"how many it in output?"
read*, nsave

allocate(vel(3,nbody))
allocate(ekin(3,nbody))
allocate(pos(3,nbody))
allocate(f(3,nbody))
allocate(velcm(3))

write(unit=*,fmt="(a)",advance="no")"time step: "
read*,dt
write(unit=*,fmt="(a)",advance="no")"n.step: "
read*,nstep
!print*,"mass of the bodies: "
!read*,massa
write(unit=*,fmt="(a)",advance="no") "Continuation run?"
read*, cont

vel=0.
if (cont) then
  read(unit=11) pos,vel
else
  do
    read(unit=10,fmt=*,iostat=ios) pos 
    if (ios<=0) exit
  end do
  if (dyn) then
!    do
!      read(unit=9,fmt=*,iostat=ios) vel 
!      if (ios<=0) exit
!    end do
    call random_number(vel)
    write(unit=*,fmt="(a)",advance="no")"vmax? "
    read*, vmax 
    vel=2*vmax*(vel-0.5)
    do i=1,3
      velcm(i) = sum(vel(i,:))/nbody
      vel(i,:) = vel(i,:) - velcm(i)
    end do
  end if
end if

call initpot(vext, nx, ny, nz, g)
call interazione(pos,nbody,f,mepot,vext,rho,nx,ny,nz,dx,temp,fakemass,nrho,ds)
 do jz=1,nz
 write(unit=12,fmt=*) jz, rho(1,1,jz), rho(3,3,jz), rho (5,5,jz)
end do
do it = 1,nstep
  if (dyn) then
    do i = 1,nbody
      pos(:,i) = pos(:,i) + vel(:,i) * dt + 0.5* f(:,i)/massa * dt**2
      vel(:,i) = vel(:,i) + 0.5 * dt * f(:,i)/massa
    end do
    call interazione(pos,nbody,f,mepot,vext,rho,nx,ny,nz,dx,temp,fakemass,nrho,ds)
    do i=1,nbody
      vel(:,i) = vel(:,i) + 0.5 * dt * f(:,i)/massa
      ekin(:,i) = 0.5 * massa * (vel(:,i))**2
    end do
    mekin = sum(ekin)
    if (mod(it,nstep/nsave).eq.0) then
      write(unit=1,fmt=*)it,it*dt,pos,vel
      write(unit=2,fmt=*)it,mekin,mepot,mekin+mepot
      write(unit=12,fmt=*)it, rho
    end if
    if (anneal) vel=alfa*vel
  else
    do i = 1,nbody
      pos(:,i) = pos(:,i) + f(:,i)/massa * dt**2
    end do
    call interazione(pos,nbody,f,mepot,vext,rho,nx,ny,nz,dx,temp,fakemass,nrho,ds)
    if (mod(it,nstep/nsave).eq.0) write(unit=1,fmt=*)it,pos
    if (mod(it,nstep/nsave).eq.0) write(unit=2,fmt=*)it,mepot
  end if
end do

do i=1,nbody
do j=i+1,nbody
write(unit=3,fmt=*) i,j,sqrt(dot_product(pos(:,i)-pos(:,j),pos(:,i)-pos(:,j)))
end do
end do

write(unit=4) pos,vel

end program corpi3d