!  norma ritorna un numero random estratto da una distribuzione normale con 
!  media zero e  varianza 2. Il metodo e'  approssimato ma sufficiente per 
!  le necessita'  del moto browniano
subroutine norma(x)
implicit none
integer :: i,ncalls
real,dimension(12) :: harvest12
real,intent(out)   :: x
call random_number(harvest12)
x= sqrt(2.0)*sum(harvest12-0.5)
end subroutine norma


PROGRAM Brown
 IMPLICIT NONE
 INTEGER, PARAMETER  :: kr=selected_int_kind(6), npart=100
 REAL(kind=kr),DIMENSION(2,npart)   :: pos,pos0,vel,f
 REAL,DIMENSION(2)                  :: harvest  ! variabile per i n. random
 REAL                               :: dt,gamma,t,w,msq
 REAL,DIMENSION(npart)              :: mass
 INTEGER                            :: it,nit,i,j
 
 WRITE(unit=*,fmt="(a)",advance="no")"Inserisci l'incremento di tempo : "
 READ*,dt
 WRITE(unit=*,fmt="(a)",advance="no")"Inserisci il numero di iterazioni che vuoi fare : "
 READ*,nit 
 WRITE(unit=*,fmt="(a)",advance="no")"Inserisci la massa  : "
 READ*,mass(1)
 mass(2:npart)=mass(1)
 WRITE(unit=*,fmt="(a)",advance="no")"Inserisci gamma e kT  : "
 READ*,gamma,t
! WRITE(unit=*,fmt="(a)",advance="no")"Inserisci la posizione iniziale : "
! READ*,pos0(1,1),pos0(2,1)

 vel=0    !   qui si parte dall'  origine con vel. nulla
 pos0=0
 
 it=0       
 pos = pos0
 
! CALL interazione(pos,f)  ! se mai si volesse mettere un'  interazione
f = 0                      ! qui non c'e':  solo attrito e forza random
WRITE(unit=1,fmt=*)"# iteraz., tempo, pos_x, pos_y, vel_x, vel_y della prima particella"
WRITE(unit=1,fmt=*)it,it*dt,pos(1,1),pos(2,1),vel(1,1),vel(2,1)

 DO it=1,nit
    DO j=1,npart
       DO i=1,2
          call norma(w)
          vel(i,j) = vel(i,j)*( 1 - gamma*dt/mass(j))  + dt  * f(i,j)/mass(j) &
                     + w*sqrt(gamma*t*dt)/mass(j)
          pos(i,j) = pos(i,j) + vel(i,j) * dt 
       END DO
       
    END DO 
    
    !CALL interazione(pos,f)
    msq = sum( (pos - pos0)**2 )/npart
        
    WRITE(unit=1,fmt=*)it,it*dt,pos(1,1),pos(2,1),vel(1,1),vel(2,1) 
    WRITE(unit=2,fmt=*)it,it*dt,msq
    
 END DO

END PROGRAM Brown