博文
高斯和随机数生成及验证 random.f
||
program gauss_pro
doubleprecision b,c,d,min,max,x
doubleprecision RANF,gauss,RAND3,rand0,rand
doubleprecision a(12000),g(12000),k(12000),f(12000)
integer g1(12000),k1(12000),a1(12000),f1(12000)
c min=15.0
c max=-10.0
x=0.002
open(UNIT=88,file='gauss.dat',status='unknown')
do i=1,200
write(*,*) rand0()
a(i)=GAUSS(1,x,0)
g(i)=GAUSS(2,x,-3)
k(i)=GAUSS(3,x,0)
f(i)=Gauss(4,x,0)
c if(min.gt.a(i)) min=a(i)
c if(max.lt.a(i)) max=a(i)
enddo
b=-4.0-0.02
do i=1,80
c=b+0.1*i
d=c+0.1
f1(i)=0
g1(i)=0
k1(i)=0
a1(i)=0
do j=1,200
if(a(j).ge.c.and.a(j).lt.d) a1(i)=a1(i)+1
if(g(j).ge.c.and.g(j).lt.d) g1(i)=g1(i)+1
if(k(j).ge.c.and.k(j).lt.d) k1(i)=k1(i)+1
if(f(j).ge.c.and.f(j).lt.d) f1(i)=f1(i)+1
enddo
write(88,'(2x,f10.5,2x,4(i5,2x))') c,a1(i),g1(i),k1(i),f1(i)
enddo
c write(*,*) max,min
close(88)
do i=1,100
write(22,*) i,gauss00(2.0d-3,0.5d0,1.0d-3,1.0d-3)
enddo
end
FUNCTION GAUSS (n,x,size)
C *******************************************************************
C ** RANDOM VARIATE FROM THE STANDARD NORMAL DISTRIBUTION. **
C ** **
C ** THE DISTRIBUTION IS GAUSSIAN WITH ZERO MEAN AND UNIT VARIANCE.**
C ** **
C ** REFERENCE: **
C ** **
C ** KNUTH D, THE ART OF COMPUTER PROGRAMMING, (2ND EDITION **
C ** ADDISON-WESLEY), 1978 **
C ** **
C ** ROUTINE REFERENCED: **
C ** **
C ** REAL FUNCTION RANF ( DUMMY ) **
C ** RETURNS A UNIFORM RANDOM VARIATE ON THE RANGE ZERO TO ONE **
C *******************************************************************
doubleprecision A1, A3, A5, A7, A9
doubleprecision A10, A30, A50, A70, A90
PARAMETER ( A10 = 3.949846138, A30 = 0.252408784 )
PARAMETER ( A50 = 0.076542912, A70 = 0.008355968 )
PARAMETER ( A90 = 0.029899776 )
doubleprecision SUM, R, R2
doubleprecision rand,ranf,rand3,rand0,x
INTEGER I,n,size
a1=a10*10**real(size)
a3=a30*10**real(size)
a5=a50*10**real(size)
a7=a70*10**real(size)
a9=a90*10**real(size)
SUM = 0.0
select case(n)
case(1)
DO I = 1, 12
SUM = SUM +rand()
enddo
case(2)
DO I = 1, 12
SUM = SUM +ranf()
enddo
case(3)
DO I = 1, 12
SUM = SUM +rand0()
enddo
case(4)
DO I = 1, 12
SUM = SUM +rand3()
enddo
end select
R = ( SUM - 6.0 ) / 4.0
R2 = R * R
GAUSS = (((( A9 * R2 + A7 ) * R2 + A5 ) * R2 + A3 ) * R2 +A1 )
+ * R
gauss=gauss+x
RETURN
END
C-----------------------------------------------------------------------C
C RANDOM NUMBER
C-----------------------------------------------------------------------C
doubleprecision FUNCTION RANF ()
INTEGER L, C, M
PARAMETER ( L = 1029, C = 221591, M = 1048576 )
INTEGER SEED
SAVE SEED
DATA SEED / 0 /
SEED = MOD ( SEED * L + C, M )
RANF = REAL ( SEED ) / M
RETURN
END
C---------------------------------------------------------------C
c Randdom number c
C---------------------------------------------------------------C
doubleprecision FUNCTION RAND3 ()
INTEGER*4 IR1,JJ1,KK1,MM1
character*20 time_now,hh,mm,ss
integer*4 i_h,i_m,i_s,K
SAVE IR1,K
JJ1=5243
KK1=55397
MM1=262139
if(K.ne.10001)K=0
if(K.eq.0)then
call time(time_now)
hh=time_now(1:2)
mm=time_now(4:5)
ss=time_now(7:8)
read(hh,*)I_h
read(mm,*)I_m
read(ss,*)I_s
IR1=i_h*3600+i_m*60+i_s
K=10001
endif
IR1=Iabs(MOD(IR1*JJ1+KK1,MM1))
RAND3=(REAL(IR1)+0.5)/REAL(MM1)
END
C---------------------------------------------------------------C
c Randdom number c
C---------------------------------------------------------------C
doubleprecision FUNCTION RAND0()
doubleprecision x
call random_seed ()
c 系统根据日期和时间随机地提供种子
call random_number (x)
c 每次的随机数就都不一样了
rand0=x
end function
C-----------------------------------------------------------------------C
C-----------------------------------------------------------------------C
doubleprecision FUNCTION GAUSS00(a,b,c,delta)
implicit doubleprecision (a-h,o-z)
parameter (pi=3.14159 )
gauss00=b*delta/(pi*(c*ranf()-a)**2+delta*delta)
end FUNCTION
https://blog.sciencenet.cn/blog-588243-481917.html
上一篇：将题录从endnote导到noteexpress中 endnote2ne.f90
下一篇：dmol 结构沿虚频方向移动 dmol_ima_freq_car.f90
