舍选法举例 subroutine acc rej PAW>call acc rej.f real rvec(1) call hbook1(10,'x(r)',100,0.,10.,0) PAW zone 13;h/pl10;h/pl 20;h/pl 30 call hbook1(20,'x(r)',100,0.,10.,0.) call hbook230,'fx)vs.x(r)',100,0,10.,100,0.,1.1,0.) if(z.gt.fmax)then fmax=-999. fmax=z*1.2 doi=1,100 write(6,*)'z greater than fmax' call ranmar(rvec,1) end if r=0+rvec(1)*(10.-0.) call hfill(10,r,0.,1.0) f=0.5*exp(-r/2.) 找f(x)的最大值 call ranmar(rvec,1) if(fmax.It.f)fmax=f u=rvec(1)*fmax end do if(u.lt.z)then fmax=1.2*fmax call hfill(20,r,0.,1.0) ntot=0 call hfill(30,r,u,1.0) doi=1,10000 ntot=ntot+1 call ranmar(rvec,1) end if r=0+vec(1)*(10.-0.) end do z=0.5*exp(-r/2.) write(6,*)'ntot=',ntot return end舍选法举例 subroutine acc_rej real rvec(1) call hbook1(10,'x(r)',100,0.,10.,0.) call hbook1(20,'x(r)',100,0.,10.,0.) call hbook2(30,'f(x) vs. x(r)',100,0.,10.,100,0.,1.1,0.) fmax=-999. do i=1,100 call ran mar(rvec,1) r=0+rvec(1)*(10.-0.) f=0.5*exp(-r/2.) if(fmax.lt.f)fmax=f end do fmax=1.2*fmax ntot=0 do i=1,10000 call ran mar(rvec,1) r=0+rvec(1)*(10.-0.) z=0.5*exp(-r/2.) if(z.gt.fmax)then fmax=z*1.2 writ e(6,*)'z greater than fmax' end if call hfill(10,r,0.,1.0) call ranmar(rvec,1) u=rvec(1)*fmax if(u.lt.z)then call hfill(20,r,0.,1.0) call hfill(30,r,u,1.0) ntot=ntot+1 end if end do write(6,*)'ntot=',ntot return end 找 f ( x )的最大值 PAW zone13; h/pl10; h/pl 20; h/pl 30 PAW call acc_rej.f > >