program Fortran_Figure7_5c implicit none save integer, parameter :: long = selected_real_kind(15,50) integer, parameter :: ranpom=1000000 integer, parameter :: ranpom2=500000 integer, parameter :: nomesh=100 real(kind=long), parameter :: DiffConst=1.0_long real(kind=long), parameter :: a1=0.5_long real(kind=long), parameter :: a2=0.0_long real(kind=long), parameter :: L=4.0_long real(kind=long), parameter :: radius=0.3_long real(kind=long), dimension(2,ranpom2) :: RanNum real(kind=long), dimension(ranpom) :: RanNumNorm real(kind=long), dimension(ranpom) :: RanNumUnif real(kind=long), dimension(2) :: xpos,ypos,zpos real(kind=long), dimension(nomesh,nomesh) :: pdfxy integer :: irn,iru,ii,jj,zz,bouncross real(kind=long) :: alpha,time,dt,pom1,pom2,speedaver,dx,angle real(kind=long) :: fx,fy,fz,distancecube,r 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(RanNum) call random_number(RanNumUnif) 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 iru=0 bouncross=0 alpha=sqrt(2.0_long*DiffConst) xpos(1)=0.0_long xpos(2)=1.0_long xpos(1)=0.0_long xpos(2)=1.0_long do ii=1,nomesh do jj=1,nomesh pdfxy(ii,jj)=0.0_long end do end do time=0.0_long dt=0.001_long do while (time<100000.0_long) time=time+dt ! interaction term between ions distancecube=(xpos(1)-xpos(2))*(xpos(1)-xpos(2)) distancecube=distancecube+(ypos(1)-ypos(2))*(ypos(1)-ypos(2)) distancecube=distancecube+(zpos(1)-zpos(2))*(zpos(1)-zpos(2)) distancecube=sqrt(distancecube)*distancecube fx=a2*(xpos(1)-xpos(2))/distancecube fy=a2*(ypos(1)-ypos(2))/distancecube fz=a2*(zpos(1)-zpos(2))/distancecube ! evolve first molecule over one time step 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 xpos(1)=xpos(1)+(fx+a1)*dt+alpha*RanNumNorm(irn)*sqrt(dt) 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 ypos(1)=ypos(1)+fy*dt+alpha*RanNumNorm(irn)*sqrt(dt) 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 zpos(1)=zpos(1)+fz*dt+alpha*RanNumNorm(irn)*sqrt(dt) ! evolve first molecule over one time step 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 xpos(2)=xpos(2)+(-fx+a1)*dt+alpha*RanNumNorm(irn)*sqrt(dt) 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 ypos(2)=ypos(2)-fy*dt+alpha*RanNumNorm(irn)*sqrt(dt) 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 zpos(2)=zpos(2)-fz*dt+alpha*RanNumNorm(irn)*sqrt(dt) ! boundary conditions do jj=1,2 if (xpos(jj)>L) then xpos(jj)=xpos(jj)-L bouncross=bouncross+1 end if if (xpos(jj)<0.0_long) then xpos(jj)=xpos(jj)+L bouncross=bouncross-1 end if r=sqrt(ypos(jj)*ypos(jj)+zpos(jj)*zpos(jj)) if (r>radius) then if (r<2.0_long*radius) then ypos(jj)=ypos(jj)*(2.0_long*radius-r)/r zpos(jj)=zpos(jj)*(2.0_long*radius-r)/r else ypos(jj)=0.0_long zpos(jj)=0.0_long end if end if end do ! saving ii=nint(real(nomesh)*xpos(1)/L+0.5_long) jj=nint(real(nomesh)*xpos(2)/L+0.5_long) if (ii>nomesh) then ii=nomesh end if if (ii<1) then ii=1 end if if (jj>nomesh) then jj=nomesh end if if (jj<1) then jj=1 end if pdfxy(ii,jj)=pdfxy(ii,jj)+1.0_long end do dx=L/real(nomesh) do ii=1,nomesh do jj=1,nomesh pdfxy(ii,jj)=pdfxy(ii,jj)*dt/time/dx/dx end do end do speedaver=real(bouncross)*L/time/2.0_long open(unit=1, file='data_Figure7_5_c.dat', status='replace') 100 format(ES15.6E3) do ii=1,nomesh do jj=1,nomesh write(1,100) pdfxy(ii,jj) end do end do close(unit=1) write(*,*) bouncross,speedaver end program Fortran_Figure7_5c