求路过的大神看一下为什么这样算不出来abc1,是定义的不够合理吗
[Fortran] 纯文本查看 复制代码 program main
implicit real(kind=8)(a-h,o-z)
pi=3.1416926d0
vm_ini=2.0d0
abc1=dsqrt(-dlog(rf(0)))
abc2=2.0d0*pi*rf(0)
p=abc1*dsin(abc2)*vm_ini
end
function rf(idum)
implicit real(kind=8)(a-h,o-z)
save ma,inext,inextp
parameter (mbig=1000000000,mseed=161803398,mz=0,fac=1.e-9)
dimension ma(55)
data iff/0/
if(idum<0.or.iff==0) then
iff=1
mj=mseed-iabs(idum)
mj=mod(mj,mbig)
ma(55)=mj
mk=1
do 50 i=1,54
ii=mod(21*i,55)
ma(ii)=mk
mk=mj-mk
if(mk<mz) mk=mk+mbig
mj=ma(ii)
50 continue
do 100 k=1,4
do 60 i=1,55
ma(i)=ma(i)-ma(1+mod(i+30,55))
if(ma(i)<mz) ma(i)=ma(i)+mbig
60 continue
100 continue
inext=0
inextp=31
endif
200 inext=inext+1
if(inext==56) inext=1
inextp=inextp+1
if(inextp==56) inextp=1
mj=ma(inext)-ma(inextp)
if(mj<mz) mj=mj+mbig
ma(inext)=mj
rf=mj*fac
if(rf>1.e-8.and.rf<0.999999999) return
goto 200
pause
end
|