program Fortran_Figure6_4_a implicit none save integer, parameter :: long = selected_real_kind(15,50) integer, parameter :: ranmax=1000000 integer, parameter :: maxocc=100 integer, parameter :: maxK=50 integer, parameter :: numberofrealizations=10000000 real(kind=long), parameter :: D=10.0_long real(kind=long), parameter :: k1=1.0_long real(kind=long), parameter :: k2=8.0_long real(kind=long), parameter :: L=1.0_long real(kind=long), dimension(ranmax) :: RanNumUnif integer, dimension(maxocc,4) :: positions ! compartment indices real(kind=long), dimension(maxK) :: averageA real(kind=long) :: time,h,Aone,k1h,diff,a0,savetime,ss integer :: ii,iru,maxused,nr,total,K,jj,pom1,pom2,pom3,pp,ll logical :: pombool,pombool2 integer size, seed(2) data seed /10, 10/ size = 2 call random_seed(SIZE=size) call random_seed(PUT=seed(1:size)) call random_number(RanNumUnif) do ii=1,maxK averageA(ii)=0.0_long end do iru=0 K=1 time=0.0_long savetime=0.01_long Aone=real(int(sqrt(k2/(2.0_long*k1)))) a0=k1*Aone*(Aone-1.0_long)+k2 iru=iru+1 if (iru>ranmax) then call random_number(RanNumUnif) iru=1 end if time=time+(1.0_long/a0)*log(1/RanNumUnif(iru)) do nr=1,numberofrealizations do while (timeranmax) then call random_number(RanNumUnif) iru=1 end if if (k2>RanNumUnif(iru)*a0) then Aone=Aone+1.0_long else Aone=Aone-2.0_long end if a0=k1*Aone*(Aone-1.0_long)+k2 iru=iru+1 if (iru>ranmax) then call random_number(RanNumUnif) iru=1 end if time=time+(1.0_long/a0)*log(1/RanNumUnif(iru)) end do averageA(1)=averageA(1)+Aone end do averageA(1)=averageA(1)/real(numberofrealizations) write(*,*) K,averageA(K) do K=2,maxK time=0.0_long h=L/real(K) diff=D/(h*h) k1h=k1/(h*h*h) do ii=1,2 iru=iru+1 if (iru>ranmax) then call random_number(RanNumUnif) iru=1 end if pom1=int(RanNumUnif(iru)*K)+1 iru=iru+1 if (iru>ranmax) then call random_number(RanNumUnif) iru=1 end if pom2=int(RanNumUnif(iru)*K)+1 iru=iru+1 if (iru>ranmax) then call random_number(RanNumUnif) iru=1 end if pom3=int(RanNumUnif(iru)*K)+1 if (ii.EQ.1) then maxused=1 positions(maxused,1)=pom1 positions(maxused,2)=pom2 positions(maxused,3)=pom3 positions(maxused,4)=1 else pombool=.TRUE. jj=0 do while (pombool.AND.(jj1)) then a0=a0+diff*positions(ii,4) endif if ((positions(ii,2)1)) then a0=a0+diff*positions(ii,4) endif if ((positions(ii,3)1)) then a0=a0+diff*positions(ii,4) endif enddo iru=iru+1 if (iru>ranmax) then call random_number(RanNumUnif) iru=1 end if time=time+(1.0_long/a0)*log(1/RanNumUnif(iru)) do nr=1,numberofrealizations do while (timeranmax) then call random_number(RanNumUnif) iru=1 end if ss=k2 if (ss>RanNumUnif(iru)*a0) then iru=iru+1 if (iru>ranmax) then call random_number(RanNumUnif) iru=1 end if pom1=int(RanNumUnif(iru)*K)+1 iru=iru+1 if (iru>ranmax) then call random_number(RanNumUnif) iru=1 end if pom2=int(RanNumUnif(iru)*K)+1 iru=iru+1 if (iru>ranmax) then call random_number(RanNumUnif) iru=1 end if pom3=int(RanNumUnif(iru)*K)+1 if (maxused.EQ.0) then maxused=1 positions(maxused,1)=pom1 positions(maxused,2)=pom2 positions(maxused,3)=pom3 positions(maxused,4)=1 else pombool=.TRUE. jj=0 do while (pombool.AND.(jj1)) then a0=a0+diff endif if ((pom21)) then a0=a0+diff endif if ((pom31)) then a0=a0+diff endif else pombool=.TRUE. ii=0 do while (pombool.AND.(ii1) then ss=ss+k1h*real(positions(ii,4)*(positions(ii,4)-1)) if (ss>RanNumUnif(iru)*a0) then a0=a0-k1h*real(positions(ii,4)*(positions(ii,4)-1)) positions(ii,4)=positions(ii,4)-2 a0=a0+k1h*real(positions(ii,4)*(positions(ii,4)-1)) a0=a0-6.0_long*diff if ((positions(ii,1)1)) then a0=a0-2.0_long*diff endif if ((positions(ii,2)1)) then a0=a0-2.0_long*diff endif if ((positions(ii,3)1)) then a0=a0-2.0_long*diff endif if (positions(ii,4).EQ.0) then if (ii1) then ss=ss+diff*real(positions(ii,4)) if (ss>RanNumUnif(iru)*a0) then pom1=positions(ii,1)-1 pom2=positions(ii,2) pom3=positions(ii,3) pombool2=.TRUE. jj=0 do while (pombool2.AND.(jjRanNumUnif(iru)*a0) then pom1=positions(ii,1)+1 pom2=positions(ii,2) pom3=positions(ii,3) pombool2=.TRUE. jj=0 do while (pombool2.AND.(jj1) then ss=ss+diff*real(positions(ii,4)) if (ss>RanNumUnif(iru)*a0) then pom1=positions(ii,1) pom2=positions(ii,2)-1 pom3=positions(ii,3) pombool2=.TRUE. jj=0 do while (pombool2.AND.(jjRanNumUnif(iru)*a0) then pom1=positions(ii,1) pom2=positions(ii,2)+1 pom3=positions(ii,3) pombool2=.TRUE. jj=0 do while (pombool2.AND.(jj1) then ss=ss+diff*real(positions(ii,4)) if (ss>RanNumUnif(iru)*a0) then pom1=positions(ii,1) pom2=positions(ii,2) pom3=positions(ii,3)-1 pombool2=.TRUE. jj=0 do while (pombool2.AND.(jjRanNumUnif(iru)*a0) then pom1=positions(ii,1) pom2=positions(ii,2) pom3=positions(ii,3)+1 pombool2=.TRUE. jj=0 do while (pombool2.AND.(jjranmax) then call random_number(RanNumUnif) iru=1 end if time=time+(1.0_long/a0)*log(1/RanNumUnif(iru)) enddo total=0 do ii=1,maxused total=total+positions(ii,4) enddo averageA(K)=averageA(K)+real(total) end do averageA(K)=averageA(K)/real(numberofrealizations) write(*,*) K,averageA(K) open(unit=1, file='data_Figure6_4_a.dat', status='replace') 100 format(i6,X,e12.5) do ii=1,K write(1,100) ii,averageA(ii) end do close(unit=1) end do end program Fortran_Figure6_4_a