Program Vlasov 
  implicit None
  real(8),dimension(:,:,:),allocatable::Force_x,Force_y,Force_z
  real(8) M(500,500,500)
  
  real::eps,delta=0.001
  integer::i,j,N,k,N_particules,lx,ly,lz,A,t_simu,u,v,w,file_iterator,file_iterator_particules
  integer, external :: sphere
  real(8)::Density, Mean_Field,Fx,Fy,Fz
  real(8),dimension(3):: Force
  real(8) x,y,z,Normalisation, imp_cdm_1,pos_cdm_1,imp_cdm_1,pos_cdm_2
  real(8),dimension(:,:),allocatable:: R_i, L_i
  real(8),dimension(:,:),allocatable:: P_i
  real(8), parameter::  rho_0 = 0.17d0, m_part = 939.d0, hbarc = 197.3d0, L= 0.5
  real(8) pi, r0, pf,r,init_imp
  pi = acos(-1.d0)
  r0 = (3.d0/(4.d0*pi*rho_0))**(1.d0/3)
  A = 100
  pf = hbarc*(3.d0/2*pi**2*rho_0)**(1.d0/3)
  r=r0*A**(1.d0/3)
   write(*,*) 'Extension Fermi', r
  N_particules =100000
  N = int(100/L)
  allocate(R_i(3,N_particules))
  allocate(L_i(3,N_particules))
  allocate(P_i(3,N_particules))
  allocate(Force_x(N,N,N))
  allocate(Force_y(N,N,N))
  allocate(Force_z(N,N,N))  
i = 1
  open(5,file='Rand')
  open(11,file='Trajectories')
  do while ( i <= N_particules)
    call random_number(x)
    call random_number(y)
    call random_number(z)
    R_i(1,i) = (2*x-1)*r
    R_i(2,i) = (2*y-1)*r 
    R_i(3,i) = (2*z-1)*r 
    if (R_i(1,i)**2+R_i(2,i)**2+R_i(3,i)**2 <= r**2) then
      write(5,*) i,R_i(:,i)+N*L/2
      write(5,*) 
      i = i+1
    end if 
   enddo
   

 i =1  
 do while ( i <= N_particules)
    call random_number(x)
    call random_number(y)
    call random_number(z)
    P_i(1,i) = (2*x-1)*pf
    P_i(2,i) = (2*y-1)*pf
    P_i(3,i) = (2*z-1)*pf
    if (P_i(1,i)**2+P_i(2,i)**2+P_i(3,i)**2 <= pf**2) then
      !P_i(:,i) = P_i(:,i) + 5 * R_i(:,i)
      !P_i(1,i) = P_i(1,i) + 25 *R_i(2,i)
      !P_i(2,i) = P_i(2,i) + 25 * R_i(1,i)
      i = i+1
    end if 
   enddo
write(*,*) 'Initial Impulsion (MeV)'
read(*,*) init_imp

  forall (i = 1:N_particules/2)
   R_i(1,i) = R_i(1,i) - 25.d0
   P_i(1,i)=P_i(1,i) + init_imp
   end forall


  
   open(1,file='Space.mat')
   open(2,file='Density.mat')
   open(3,file='Meanfield.mat')
   open(4,file='Force.mat')
write(*,*) 'Taille Boite', N
   Normalisation = N_particules * L**3

M = 0
do i=1,N_particules
  call near_index(N_particules,R_i(:,i),N,lx,ly,lz)
if (lx> 0 .and. ly>0 .and. lz>0 .and. lx<=n .and. ly<=n .and. lz<=n) then
  M(lx,ly,lz) = M(lx,ly,lz) + 1
endif
  enddo
   
   do i = 1,N
      do j=1,N
	write(1,*) i ,j , N/2, M(i,j,N/2)   
 	call calcul_density (A,i,j,N/2, M, Density, Mean_Field, Normalisation,N)
     	write(2,*) i,j,N/2, Density
     	write(3,*) i,j,N/2, Mean_Field
     	call calcul_force(A,i,j,N/2, M, Density, Mean_Field, Normalisation,Force,N)
        Force_x(i,j,N/2) = Force(1)
	Force_y(i,j,N/2) = Force(2)
	Force_z(i,j,N/2) = Force(3)
     	write(4,'(3(i4,1x),3(1x,f9.2))') i,j,N/2, Force(1), Force(2), Force(3)
     enddo
       write(1,*)
       write(2,*)
       write(3,*)
       write(4,*)
   enddo
   
   
call calcul_density(A,N/2,N/2,N/2,M,Density,Mean_Field,Normalisation,N)
  write(*,*) 'Densité au centre', Density
  
    
file_iterator = 0
 do t_simu = 1, 6000
	write(*,*) t_simu
	imp_cdm = 0.
	pos_cdm = 0.
     	call equation_of_motion(R_i,P_i,N_particules,t_simu,M,Density,Mean_Field,Normalisation,Force_x,Force_y,Force_z,N)
 M=0.

do i=1,N_particules/2
call near_index(N_particules,R_i(:,i),N,lx,ly,lz)
if (lx> 0 .and. ly>0 .and. lz>0 .and. lx<=n .and. ly<=n .and. lz<=n) then
	M(lx,ly,lz) = M(lx,ly,lz) + 1
	pos_cdm_1 = (pos_cdm_1 + R_i(1,i)**2 + R_i(2,i)**2+R_i(3,i)**2)
	imp_cdm_1 = (imp_cdm_1 + P_i(1,i)**2 + P_i(2,i)**2 + P_i(3,i)**2)
endif
enddo

do i=N_particules/2+1,N_particules
call near_index(N_particules,R_i(:,i),N,lx,ly,lz)
if (lx> 0 .and. ly>0 .and. lz>0 .and. lx<=n .and. ly<=n .and. lz<=n) then
	M(lx,ly,lz) = M(lx,ly,lz) + 1
	pos_cdm_2 = (pos_cdm_2 + R_i(1,i)**2 + R_i(2,i)**2+R_i(3,i)**2)
	imp_cdm_2 = (imp_cdm_2 + P_i(1,i)**2 + P_i(2,i)**2 + P_i(3,i)**2)
endif
enddo

    do u= 1,N
    	do v= 1,N
    		call calcul_density(A,u,v,N/2,M,Density,Mean_Field,Normalisation,N)
		if (mod(t_simu,10) /= 0) then
   		write(300+file_iterator,'(3(i4),3(1x,f9.2))') u,v, N/2, M(u,v,N/2), Density, Mean_field
		endif
   	enddo
        write(300+file_iterator,*) 
 enddo
 close(300+file_iterator)
file_iterator = file_iterator + 1
write(25,*) t_simu, pos_cdm_1, imp_cdm_1,imp_cdm_2,pos_cdm_2
enddo


    
deallocate(R_i)
deallocate(P_i)
end Program Vlasov

    
subroutine equation_of_motion(R,P,N_particules,t_simu,T,Density,Mean_Field,Normalisation,Force_x,Force_y,Force_z,N)
    implicit none 
    integer::N,i,j,k,l,N_particules,t_simu,lx,ly,lz
    real(8)::Density, Mean_Field,delta_t,Fx,Fy,Fz
    real(8) :: Normalisation
    real(8),dimension(3):: Force
    real(8),dimension(N,N,N)::Force_x,Force_y,Force_z
    real(8),dimension(3,N_Particules) :: R, P
    real(8),dimension(N,N,N)::T
    real(8),parameter :: masse = 939.d0
    delta_t = 0.1

  do i=1,N_particules  

	call near_index(N_particules,R(:,i),N,lx,ly,lz)
if (lx> 0 .and. ly>0 .and. lz>0 .and. lx<=n .and. ly<=n .and. lz<=n) then
	call calcul_force(100,lx,ly,lz,T,Density,Mean_Field,Normalisation,Force,N)
 
 
    R(1,i)=R(1,i)+P(1,i)*delta_t/masse+0.5d0*Force(1)*(delta_t**2)/masse
    R(2,i)=R(2,i)+P(2,i)*delta_t/masse+0.5d0*Force(2)*(delta_t**2)/masse
    R(3,i)=R(3,i)+P(3,i)*delta_t/masse+0.5d0*Force(3)*(delta_t**2)/masse
  
    P(1,i)=P(1,i)+0.5d0*Force(1)*delta_t
    P(2,i)=P(2,i)+0.5d0*Force(2)*delta_t
    P(3,i)=P(3,i)+0.5d0*Force(3)*delta_t

endif
enddo

 T=0.
do i=1,N_particules
call near_index(N_particules,R(:,i),N,lx,ly,lz)
	
	if (lx> 1 .and. ly>1 .and. lz>1 .and. lx<=n .and. ly<=n .and. lz<=n) then
	T(lx,ly,lz) = T(lx,ly,lz) + 1
	
	endif
enddo

do i= 1,N_particules
	call near_index(N_particules,R(:,i),N,lx,ly,lz)
if (lx> 0 .and. ly>0 .and. lz>0 .and. lx<=n .and. ly<=n .and. lz<=n) then
 call calcul_force(100,lx,ly,lz,T,Density,Mean_Field,Normalisation,Force,N)

    P(1,i)=P(1,i)+0.5d0*Force(1)*delta_t
    P(2,i)=P(2,i)+0.5d0*Force(2)*delta_t
    P(3,i)=P(3,i)+0.5d0*Force(3)*delta_t
endif
enddo

end 



subroutine near_index(N_particules,R_i,N,lx,ly,lz)
implicit none
integer::N,i,N_particules,lx,ly,lz
real(8),dimension(3) :: R_i, L_i
real(8), parameter::L=0.5
        L_i(:) = R_i(:) + N*L/2
	lx = nint(L_i(1)/L) + 1
	ly = nint(L_i(2)/L) + 1
	lz = nint(L_i(3)/L) + 1 
end
 

  subroutine calcul_density(A,x,y,z,M,Density,Mean_Field,Normalisation,n)
    implicit none
    real(8) M(500,500,500)
    integer::x,y,z,n,k1,k2,k3,A
    real(8),dimension(3)::poids
    real(8)::rho, Density, Mean_Field
    real(8):: Normalisation
    real(8), parameter::rho_0=0.17d0
    poids(1) = 0.25
    poids(2) = 0.5
    poids(3) = 0.25
    rho = 0.
    Mean_Field = 0.
  if (x> 1 .and. y>1 .and. z>1 .and. x<=n-1 .and. y<=n-1 .and. z<=n-1) then
    do k1 = -1,1
      do k2 =-1,1
	do k3 = -1,1
	   rho = rho + poids(k1+2)*poids(k2+2)*poids(k3+2)*M(x+k1,y+k2,z+k3)    
	enddo
      enddo
    enddo
endif
    Density = A*rho/(Normalisation)
    Mean_Field = -356.* Density/rho_0 + 303.*(Density/rho_0)**(7.d0/6.d0)
    end
    

    
  subroutine calcul_force (A,x,y,z,T,Density,Mean_Field,Normalisation,Force,n)
      implicit none
      real(8) T(n,n,n)
      integer::x,y,z,n,A
      real(8)::Density, Mean_Field
      real(8) :: Normalisation
      real(8),dimension(2)::Ux,Uy,Uz
      real(8), parameter::L=0.5
      real(8),dimension(3):: Force
      if (x>1 .and. y>1 .and. z>1 .and. x<=n-1 .and. y<=n-1 .and. z<=n-1) then
      call calcul_density(A,x+1,y,z,T,Density,Mean_Field,Normalisation,n)
      Ux(1) = Mean_Field
      call calcul_density(A,x-1,y,z,T,Density,Mean_Field,Normalisation,n)
      Ux(2) = Mean_Field
      Force(1) =-1.d0*(Ux(1) - Ux(2))/(2*L)
      
      call calcul_density(A,x,y+1,z,T,Density,Mean_Field,Normalisation,n)
      Uy(1) = Mean_Field
      call calcul_density(A,x,y-1,z,T,Density,Mean_Field,Normalisation,n)
      Uy(2) = Mean_Field
      Force(2) = -1.d0*(Uy(1) - Uy(2))/(2*L)
      
      call calcul_density(A,x,y,z+1,T,Density,Mean_Field,Normalisation,n)
      Uz(1) = Mean_Field
      call calcul_density(A,x,y,z-1,T,Density,Mean_Field,Normalisation,n)
      Uz(2) = Mean_Field
      Force(3) = -1.d0*(Uz(1) - Uz(2))/(2*L)
      
      call calcul_density(A,x,y,z,T,Density,Mean_Field,Normalisation,n)
endif
      end
      
  
  
