ERF函数Fortran代码
下面贴出的是Fortran版本的误差函数，代码源于徐士良的那本书，但做了一定修改并长期使用过。建议大家复制后再测试一下，编程一定要谨慎！！！为避免与新版本编译器自带ERF函数名字冲突，代码中子程序并没有直接用ERF，而是f_merrf，另外还有两个依赖的子函数也贴出来了。
!误差函数
subroutine f_merrf(x,merrf)
real merrf,x,mgam2,tmp
tmp=x*x
call f_mgam2(0.5,tmp,mgam2)
if (x.ge.0.0) then
merrf=mgam2
else
merrf=-mgam2
end if
return
end subroutine
!伽马函数
subroutine f_mgam1(x,mgam1)
real mgam1,x
real y,t,s,u,a(11)
data a /0.0000677106,-0.0003442342,0.0015397681,&
-0.0024467480,0.0109736958,-0.0002109075,&
0.0742379071,0.0815782188,0.4118402518,&
0.4227843370,1.0/
integer i
if (x.le.0.) then
mgam1=-1.0
return
end if
y=x
if (y.le.1.0) then
t=1.0/(y*(y+1.0))
y=y+2.0
else if (y.le.2.0) then
t=1.0/y
y=y+1.0
else if (y.le.3.0) then
t=1.0
else
t=1.0
do while(y.gt.3.0)
y=y-1.0
t=t*y
end do
end if
s=a(1)
u=y-2.0
do i=1,10
s=s*u+a(i+1)
end do
s=s*t
mgam1=s
return
end subroutine
!不完全伽马函数
subroutine f_mgam2(a,x,mgam2)
real mgam2,a,x
real mgam1,p,q,d,s,s1,p0,q0,p1,q1,qq
integer n
if ((a.le.0.0).or.(x.lt.0.0)) then
mgam2=-1.0
end if
if (x+1.0.eq.1.0) then
mgam2=0.0
return
end if
if (x.gt.1.0d+35) then
mgam2=1.0
return
end if
q=log(x)
q=a*q
qq=exp(q)
if (x.lt.1.0+a) then
p=a
d=1.0/a
s=d
do n=1,100
p=1.0+p
d=d*x/p
s=s+d
if (abs(d).lt.abs(s)*1.0d-07) then
call f_mgam1(a,mgam1)
s=s*exp(-x)*qq/mgam1
mgam2=s
return
end if
end do
else
s=1.0/x
p0=0.0
p1=1.0
q0=1.0
q1=x
do n=1,100
p0=p1+(n-a)*p0
q0=q1+(n-a)*q0
p=x*p0+n*p1
q=x*q0+n*q1
if (abs(q)+1.0.ne.1.0) then
s1=p/q
p1=p
q1=q
if (abs((s1-s)/s1).lt.1.0d-07) then
call f_mgam1(a,mgam1)
s=s1*exp(-x)*qq/mgam1
mgam2=1.0-s
return
end if
s=s1
end if
p1=p
q1=q
end do
end if
call f_mgam1(a,mgam1)
s=1.0-s*exp(-x)*qq/mgam1
mgam2=s
return
end subroutine
https://www.google.com/favicon.icohttp://www.baidu.com/favicon.ico
页:
[1]