program serial
parameter(nmu=20,nsig=20,nb=10,ntz=16,nd=11)
parameter (nsim=2000)
integer k1,k2,k3,k4,k5,i,j,m,n
real xmu,sigma,d,b,tz,beta
data xmu1/.02/,xmu2/.4/
data sig1/.2/,sig2/4./
data b1/0./
data t1/200/,t2/1000/
data d1/1./,d2/2./
data iseed/-9/
open(unit=10,file='s10')
open(unit=20,file='S1T0s.txt')
r=ran2(iseed)
read(20,*) n
read(20,*) m
do k1=1,nmu
xmu=xmu1+(xmu2-xmu1)*(k1-1)/float(nmu-1)
do k2=1,nsig
sigma=sig1+(sig2-sig1)*(k2-1)/float(nsig-1)
do k3=1,nd
d=d1+(d2-d1)*(k3-1)/float(nd-1)
do k4=1,nb
b2=10*d
beta=b1+(b2-b1)*(k4-1)/float(nb-1)
do k5=1,ntz
tz=t1+(t2-t1)*(k5-1)/float(ntz-1)
call varser(nsim,n,m,d,beta,tz,sigma,xmu)
enddo
enddo
enddo
enddo
enddo
end
subroutine varser(nsim,n,m,d,beta,tz,sigma,xmu)
integer nsim,n,m,it,icorr,ntime
real d,beta,tz,sigma,xmu,totcor,totstp,stepsq,pcorr
real perr,avst,errvar,stvar
totcor=0
totstp=0
totsq=0
stepsq=0
do it=1,nsim
call trser(beta,tz,sigma,xmu,d,n,m,icorr,ntime)
totcor=totcor+icorr
totstp=totstp+ntime
stepsq=stepsq+ntime**2
enddo
pcorr=totcor/nsim
perr=1.-pcorr
avst=totstp/nsim
errvar=pcorr*perr/nsim
stvar=(stepsq-avst**2/nsim)/float(nsim-1)
write(10,*) perr,avst,stvar
return
end
subroutine trser(beta,tz,sigma,xmu,d,nel,ntarg,icorr,ntime)
parameter(targb=20)
parameter(nelmax=4)
integer typetr
integer nel,ntarg,icorr,ntime,ndist,iran
real beta,tz,sigma,xmu,d
dimension xmean(nelmax),xran(nelmax)
iran=1
typetr=0
ntime=0
ndist=nel-ntarg
distb0=-d*targb
do i=1,ntarg
xmean(i)=xmu
enddo
do i=1,ndist
xmean(i+ntarg)=-xmu
enddo
do i=1,nel
xran(i)=ran2(iran)
enddo
call sort2(nel,xran,xmean)
do iw=1,nel
distb1=bcalc(nel,tz,beta,ntime,distb0)
call walker(xmean(iw),sigma,distb1,targb,nstep,typetr)
ntime=ntime+nstep
if(typetr .eq. 1 .or. ntime .gt. tz) go to 99
enddo
99 continue
if(ntarg .gt. 0) then
if(typetr .eq. 1) then
icorr=1
else
icorr=0
endif
else
if(typetr .eq. 0) then
icorr=1
else
icorr=0
endif
endif
return
end
function bcalc(nel,tz,beta,ntime,distb0)
if(nel .eq. 4) b1=distb0+2*beta
if(nel .eq. 2) b1=distb0+beta
if(nel .eq. 1) b1=distb0
timfrc=float(ntime)/tz
rem=1-timfrc
bcalc=b1*rem
return
end
subroutine walker(xm,sig,db,tb,n,idec)
iran=1
w=0
n=0
10 continue
step=gasdev(iran)*sig+xm
n=n+1
w=w+step
if(w .ge. tb) then
idec=1
return
endif
if(w .le. db) then
idec=0
return
endif
go to 10
end
FUNCTION gasdev(idum)
INTEGER idum
REAL gasdev
CU USES ran1
INTEGER iset
REAL fac,gset,rsq,v1,v2,ran1
SAVE iset,gset
DATA iset/0/
if (iset.eq.0) then
1 v1=2.*ran2(idum)-1.
v2=2.*ran2(idum)-1.
rsq=v1**2+v2**2
if(rsq.ge.1..or.rsq.eq.0.)goto 1
fac=sqrt(-2.*alog(rsq)/rsq)
gset=v1*fac
gasdev=v2*fac
iset=1
else
gasdev=gset
iset=0
endif
return
END
FUNCTION ran2(idum)
INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
REAL ran2,AM,EPS,RNMX
PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1)
PARAMETER (IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211)
PARAMETER (IR2=3791,NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7)
PARAMETER (RNMX=1.-EPS)
INTEGER idum2,j,k,iv(NTAB),iy
SAVE iv,iy,idum2
DATA idum2/123456789/, iv/NTAB*0/, iy/0/
if (idum.le.0) then
idum=max(-idum,1)
idum2=idum
do 11 j=NTAB+8,1,-1
k=idum/IQ1
idum=IA1*(idum-k*IQ1)-k*IR1
if (idum.lt.0) idum=idum+IM1
if (j.le.NTAB) iv(j)=idum
11 continue
iy=iv(1)
endif
k=idum/IQ1
idum=IA1*(idum-k*IQ1)-k*IR1
if (idum.lt.0) idum=idum+IM1
k=idum2/IQ2
idum2=IA2*(idum2-k*IQ2)-k*IR2
if (idum2.lt.0) idum2=idum2+IM2
j=1+iy/NDIV
iy=iv(j)-idum2
iv(j)=idum
if(iy.lt.1)iy=iy+IMM1
ran2=min(AM*iy,RNMX)
return
END
SUBROUTINE sort2(n,arr,brr)
INTEGER n,M,NSTACK
REAL arr(n),brr(n),b
PARAMETER (M=7,NSTACK=50)
INTEGER i,ir,j,jstack,k,l,istack(NSTACK)
REAL a,temp
jstack=0
l=1
ir=n
1 if(ir-l.lt.M)then
do 12 j=l+1,ir
a=arr(j)
b=brr(j)
do 11 i=j-1,1,-1
if(arr(i).le.a)goto 2
arr(i+1)=arr(i)
brr(i+1)=brr(i)
11 continue
i=0
2 arr(i+1)=a
brr(i+1)=b
12 continue
if(jstack.eq.0)return
ir=istack(jstack)
l=istack(jstack-1)
jstack=jstack-2
else
k=(l+ir)/2
temp=arr(k)
arr(k)=arr(l+1)
arr(l+1)=temp
temp=brr(k)
brr(k)=brr(l+1)
brr(l+1)=temp
if(arr(l+1).gt.arr(ir))then
temp=arr(l+1)
arr(l+1)=arr(ir)
arr(ir)=temp
temp=brr(l+1)
brr(l+1)=brr(ir)
brr(ir)=temp
endif
if(arr(l).gt.arr(ir))then
temp=arr(l)
arr(l)=arr(ir)
arr(ir)=temp
temp=brr(l)
brr(l)=brr(ir)
brr(ir)=temp
endif
if(arr(l+1).gt.arr(l))then
temp=arr(l+1)
arr(l+1)=arr(l)
arr(l)=temp
temp=brr(l+1)
brr(l+1)=brr(l)
brr(l)=temp
endif
i=l+1
j=ir
a=arr(l)
b=brr(l)
3 continue
i=i+1
if(arr(i).lt.a)goto 3
4 continue
j=j-1
if(arr(j).gt.a)goto 4
if(j.lt.i)goto 5
temp=arr(i)
arr(i)=arr(j)
arr(j)=temp
temp=brr(i)
brr(i)=brr(j)
brr(j)=temp
goto 3
5 arr(l)=arr(j)
arr(j)=a
brr(l)=brr(j)
brr(j)=b
jstack=jstack+2
if(jstack.gt.NSTACK)pause 'NSTACK too small in sort2'
if(ir-i+1.ge.j-l)then
istack(jstack)=ir
istack(jstack-1)=i
ir=j-1
else
istack(jstack)=j-1
istack(jstack-1)=l
l=i
endif
endif
goto 1
END
