 module global 
   double precision, parameter :: pi = 4.d0 * atan(1.d0) 
   integer, parameter :: N_samp_max = 1000 
   integer N_realiz         ! # of stochastic trajectories
   integer N_samples        ! # of time samples
   double precision a,s     ! some parameters 
   double precision time    ! the time
   double precision dt      ! the time step
   double precision dth     ! the time step**0.5
   double precision, dimension (N_samp_max) :: xm, xm2, ts
   character (len=10) seed_char
 end module global
!========================================================  
 program Milshtein 
!
! This code implement the Milshtein method to 
! integrate a one-dimensional stochastic flow 
! of the form 
!
!   dx
!   -- = f(x,t) + g(x) xi(t)
!   dt
!
! where xi(t) is a Gaussian stochastic process 
! with delta correlation and unit variance
!
! This code implements the linear equation with 
! multiplicative noise:
!
!   f(x,t) = a*x
!   g(x,t) = s*x
!
!======================================================== 
!
 use global

 integer k                ! floating sample counter
 integer N_steps          ! the total number of steps
 
 double precision x       ! the state variable
 double precision f,g,dg  ! the model functions
 double precision x0      ! the initial conditions
 double precision T_tot   ! the total integration time

 character (len=20) arg

!======================================================== 
 if ( iargc() == 1 ) then
     call getarg(1,arg)
     read (arg,*) N_realiz
 else 
     print*, ' Usage: milshtein.gfx [ # realizations ]'
     STOP
 end if
 write (seed_char,'(i10.10)') N_realiz
!========================================================
! Model parameters

 a = -0.1d0
 s = 1.0d0
 
!======================================================== 
 
 dt = 0.1d0
 dth = dsqrt(dt)
 N_samples = 100        
 T_tot = 20.d0
 N_steps = int(T_tot/dt)
 i_samp = int ( dfloat(N_steps) / dfloat(N_samples) )
 do i=1, N_samples
    xm(i) = 0.d0 
    xm2(i) = 0.d0
    ts(i) = i_samp*i*dt
 end do
 x0 = 0.01d0
 
!########################################################
! Loop over many stochastic realizations 

 DO is=1,N_realiz 
 
!======================================================== 
! Integration

 x = x0
 k = 0
 do i=1, N_steps
    call mils_tstep(x)
    if ( mod(i,i_samp) == 0 ) then 
         k = k + 1 
         print*, k, x, is
         call do_averages(x,k)
    end if     
 end do
 
!########################################################
 END DO   ! end loop over many stochastic realizations 

 call print_final
 
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 contains
 
 subroutine do_averages(x,i)
 use global 
 integer i
 double precision x
 
 xm(i)  = xm(i)  + x 
 xm2(i) = xm2(i) + x*x  
 
 end subroutine 
 
!--------------------------------------------------------
 
 subroutine print_final
 use global

 open (1,file='xm_Nreal_'//seed_char//'.dat', &
         status='unknown')
 do i=1, N_samples
    write (1,'(3f20.8)') ts(i), xm(i)/N_realiz,     &
    dsqrt((xm2(i)/N_realiz - (xm(i)/N_realiz)**2)/  &
           N_realiz)
 end do
 close (1)
 
 end subroutine 

!---------------------------------------------------------

 subroutine mils_tstep(x)
 use global 
 double precision x,r1,r2,u,f,g,dg
 external model 
 
 call random_number(r1)
 call random_number(r2)
 u = sqrt(-2.0d0 * dlog(r1)) * cos (2.d0 * pi * r2)
 
 call model(f,g,dg,x)
 x = x + g*dth*u + ( f + g*dg*u*u/2.0d0 )*dt
 
 end subroutine

!---------------------------------------------------------

 subroutine model(f,g,dg,x)
 use global 
 double precision x,f,g,dg
 
 f  = a*x 
 g  = s*x 
 dg = s
 
 end subroutine

 end program  
 
 
 
 
 
 
 
