program Fortran_Figure3_2b

implicit none

save

 integer, parameter :: long = selected_real_kind(15,50)
 integer, parameter :: ranpom=1000000
 integer, parameter :: ranpom2=500000

 real(kind=long), parameter :: k1=0.001_long
 real(kind=long), parameter :: k2=0.75_long
 real(kind=long), parameter :: k3=165.0_long
 real(kind=long), parameter :: k4=10000.0_long
 real(kind=long), parameter :: k5=200.0_long
 real(kind=long), dimension(2,ranpom2) :: RanNum 
 real(kind=long), dimension(ranpom) :: RanNumNorm 
 real(kind=long), dimension(600) :: histogram 
 real(kind=long) :: X,time,nsave,pom1,pom2,dt
 integer :: ii,zz,irn,tt,oldtt
 
 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

 time=0.0_long
 nsave=0.0_long
 dt=0.001_long
 X=100.0_long
   
 do ii=1,600
    histogram(ii)=0.0_long
 end do  
 
 oldtt=0
   
 do while (time<10000000.0_long)
          
    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 
	 
    X=X+dt*(-k1*X*X*X+k2*X*X-k3*X+k4)+k5*sqrt(dt)*RanNumNorm(irn)
    time=time+dt
     
    tt=nint(10.0*time)
    if (tt>oldtt) then
       ii=nint(X)
       histogram(ii)=histogram(ii)+1.0_long           
       nsave=nsave+1.0_long
       oldtt=tt
    end if
 end do   

 do ii=1,600
    histogram(ii)=histogram(ii)/nsave
 end do  
 
 open(unit=1, file='data_Figure3_2b.dat', status='replace')

100 format(e12.6)
 
 do ii=1,600
    write(1,100) histogram(ii)
 end do

 close(unit=1)
 
end program Fortran_Figure3_2b