program Fortran_Figure8_2_a

implicit none

save

 integer, parameter :: long = selected_real_kind(15,50)

 integer, parameter :: numberofrealizations=10
 real(kind=long), parameter :: timestep=0.01_long
 integer, parameter :: kernellength=1201
 integer, parameter :: km=1200  ! km=kernellength-1
 integer, parameter :: numberoftimepoints=20001
 real(kind=long), dimension(numberoftimepoints) :: Xvartraj,Vvartraj,Rtraj
 real(kind=long), dimension(numberoftimepoints,numberofrealizations) :: Xvartrajsave
 real(kind=long), dimension(kernellength) :: kappa
 real(kind=long), dimension(km) :: czero
 real(kind=long), dimension(km,km) :: Amatrix,cvector
 integer, parameter :: ranpom=1000000
 integer, parameter :: ranpom2=500000
   
 real(kind=long), parameter :: alpha1=1.0_long
 real(kind=long), parameter :: alpha2=1.0_long
 real(kind=long), parameter :: alpha3=1.0_long
 real(kind=long), parameter :: alpha4=1.0_long

 real(kind=long) :: pom1,pom2,sigma,alpha5
 integer :: ii,jj,ll,zz,nr,irn,tt
 real(kind=long), dimension(2,ranpom2) :: RanNum 
 real(kind=long), dimension(ranpom) :: RanNumNorm 
 
 integer size, seed(2)
 data seed /100, 100/ 
 size = 2

 call random_seed(SIZE=size)
 call random_seed(PUT=seed(1:size))
 call random_number(RanNum)
  
 do zz=1,ranpom2
    pom1=sqrt(2.0_long*log(1.0_long/RanNum(1,zz)))
    pom2=2.0_long*3.14159265358979_long*RanNum(2,zz)
    RanNumNorm(2*zz-1)=pom1*cos(pom2)
    RanNumNorm(2*zz)=pom1*sin(pom2)
 end do   
 
 irn=0
                
 do ii=1,numberoftimepoints
    do jj=1,numberofrealizations
       Xvartrajsave(ii,jj)=0.0_long
    end do
 end do
                
 sigma=(alpha3/sqrt(2.0_long*alpha2))
 
 if (alpha4>(4.0_long*alpha2)) then
    alpha5=sqrt(alpha4*abs(alpha4-4.0_long*alpha2))
    do ii=1,kernellength
       pom1=real(ii-1)*timestep/2.0_long
       kappa(ii)=alpha1*exp(-alpha4*pom1)*(alpha4*sinh(alpha5*pom1)/alpha5+cosh(alpha5*pom1))
    end do
 else
    alpha5=sqrt(alpha4*abs(alpha4-4.0_long*alpha2))
    do ii=1,kernellength
       pom1=real(ii-1)*timestep/2.0_long
       kappa(ii)=alpha1*exp(-alpha4*pom1)*(alpha4*sin(alpha5*pom1)/alpha5+cos(alpha5*pom1))
    end do
 end if

 ! computation of c vector 
 
 tt=1
    cvector(1,tt)=kappa(2)
    Amatrix(1,1)=kappa(1)
    cvector(1,tt)=cvector(1,tt)/Amatrix(1,1)
    czero(tt)=cvector(1,tt)*kappa(tt+1)
    czero(tt)=sigma*sqrt(1-czero(tt)/kappa(1))
 
 do tt=2,km
    do ii=1,tt
       cvector(ii,tt)=kappa(ii+1)
    end do
    do ii=1,tt
       do jj=1,tt
          Amatrix(ii,jj)=kappa(abs(ii-jj)+1)
       end do      
    end do
    do jj=1,tt-1
       do ii=jj+1,tt
          pom1=Amatrix(ii,jj)/Amatrix(jj,jj)
          do ll=jj,tt
	     Amatrix(ii,ll)=Amatrix(ii,ll)-Amatrix(jj,ll)*pom1
	  end do
	  cvector(ii,tt)=cvector(ii,tt)-cvector(jj,tt)*pom1
       end do
    end do
    cvector(tt,tt)=cvector(tt,tt)/Amatrix(tt,tt)
    czero(tt)=cvector(tt,tt)*kappa(tt+1)
    do ii=1,tt-1
       do jj=(tt-ii+1),tt
          cvector(tt-ii,tt)=cvector(tt-ii,tt)-Amatrix(tt-ii,jj)*cvector(jj,tt)
       end do
       cvector(tt-ii,tt)=cvector(tt-ii,tt)/Amatrix(tt-ii,tt-ii) 
       czero(tt)=czero(tt)+cvector(tt-ii,tt)*kappa(tt-ii+1)
    end do
    czero(tt)=sigma*sqrt(1-czero(tt)/kappa(1))
 end do
 
 do nr=1,numberofrealizations
    do ii=1,numberoftimepoints
       Xvartraj(ii)=0.0_long
       Vvartraj(ii)=0.0_long
    end do   

    irn=irn+1
    if (irn>ranpom) then
       call random_number(RanNum)
       do zz=1,ranpom2
          pom1=sqrt(2.0_long*log(1.0_long/RanNum(1,zz)))
          pom2=2.0_long*3.14159265358979_long*RanNum(2,zz)
	  RanNumNorm(2*zz-1)=pom1*cos(pom2)
	  RanNumNorm(2*zz)=pom1*sin(pom2)
       end do   
       irn=1
    end if
    Rtraj(1)=sigma*RanNumNorm(irn)
    Xvartraj(2)=Xvartraj(1)+Vvartraj(1)*timestep
    Vvartraj(2)=Vvartraj(1)-timestep*timestep*kappa(1)*Vvartraj(1)+Rtraj(1)*timestep

    do tt=2,km
       irn=irn+1
       if (irn>ranpom) then
          call random_number(RanNum)
          do zz=1,ranpom2
             pom1=sqrt(2.0_long*log(1.0_long/RanNum(1,zz)))
             pom2=2.0_long*3.14159265358979_long*RanNum(2,zz)
	     RanNumNorm(2*zz-1)=pom1*cos(pom2)
	     RanNumNorm(2*zz)=pom1*sin(pom2)
          end do   
          irn=1
       end if
       Rtraj(tt)=czero(tt-1)*RanNumNorm(irn)
       do ii=1,tt-1
          Rtraj(tt)=Rtraj(tt)+cvector(ii,tt-1)*Rtraj(tt-ii)
       end do	  
       Xvartraj(tt+1)=Xvartraj(tt)+Vvartraj(tt)*timestep
       pom1=0.0_long
       do ii=1,tt
          pom1=pom1+kappa(tt)*Vvartraj(tt-ii+1)
       end do       
       Vvartraj(tt+1)=Vvartraj(tt)-timestep*timestep*pom1+Rtraj(tt)*timestep
    end do   
     
    do tt=(km+1),numberoftimepoints-1
       irn=irn+1
       if (irn>ranpom) then
          call random_number(RanNum)
          do zz=1,ranpom2
             pom1=sqrt(2.0_long*log(1.0_long/RanNum(1,zz)))
             pom2=2.0_long*3.14159265358979_long*RanNum(2,zz)
	     RanNumNorm(2*zz-1)=pom1*cos(pom2)
	     RanNumNorm(2*zz)=pom1*sin(pom2)
          end do   
          irn=1
       end if
       Rtraj(tt)=czero(km)*RanNumNorm(irn)
       do ii=1,km
          Rtraj(tt)=Rtraj(tt)+cvector(ii,km)*Rtraj(tt-ii)
       end do	  
       Xvartraj(tt+1)=Xvartraj(tt)+Vvartraj(tt)*timestep
       pom1=0.0_long
       do ii=1,km+1
          pom1=pom1+kappa(ii)*Vvartraj(tt-ii+1)
       end do       
       Vvartraj(tt+1)=Vvartraj(tt)-timestep*timestep*pom1+Rtraj(tt)*timestep
    end do   

    do ii=1,numberoftimepoints
       Xvartrajsave(ii,nr)=Xvartraj(ii)
    end do
       
 end do

 100 format(e14.6)
	  
 open(unit=1, file='data_Figure8_2_a.dat', status='replace')

 do ii=1,numberofrealizations
    do tt=1,numberoftimepoints        
       write(1,100) Xvartrajsave(tt,ii)
    end do	    
 end do

 close(unit=1)	         	  

   
end program Fortran_Figure8_2_a
