|
[Fortran] 纯文本查看 复制代码 program main
dimension o(2),x(2)
real*4 o,x,c,d
data ntrap /20/
integer i
common o,x,c,d
do i=1,ntrap,1
call f1()
write(*,*) 'x(1)=',x(1),'x(2)=',x(2)
write(*,*) 'o(1)=',o(1),'o(2)=',o(2)
c=x(1)
d=x(2)
o(1)=o(1)-2.0*c/1.0*vdotn(o,x/1.0,2)
o(2)=o(2)-2.0*d/1.0*vdotn(o,x/1.0,2)
end do
stop
end
block data
implicit none
dimension o(2),x(2)
real*4 o,x,c,d
common o,x,c,d
data o(1) /-0.7071/
data o(2) /0.7071/
data c /0.0/
data d /-1.0/
end block data
subroutine f1()
dimension x(2),y(2)
common o,x,c,d
data x/2*0.0/
b=2.0
n=2
m=10
eps=1.0e-5
call dnmtc(x,n,b,m,eps,f,y) !第39行
write(*,*)
do 10 i=1,n
10 write(*,100) i,x(i)
write(*,*)
100 format(5x,'x(',i2,1x,')=',e13.6)
end subroutine
function f(x,n)
dimension x(n),o(n)
real*4 o,x,c,d
common o,x,c,d !第49行
f1=x(1)*x(1)+x(2)*x(2)-1.0
f2=o(2)*x(1)-o(2)*c-o(1)*x(2)+o(1)*d
f=sqrt(f1*f1+f2*f2)
end
subroutine dnmtc(x,n,b,m,eps,f,y)
dimension x(n),y(n)
real*4 o,x,c,d
common o,x,c,d !第58行
double precision r
real nrnd1
a=b
k=1
r=1.0d0
z=f(x,n)
10 if (a.ge.eps) then
l=l+1
do 20 i=1,n
20 y(i)=-a+2.0*a*nrnd1(r)+x(i)
z1=f(y,n)
k=k+1
if (z1.ge.z) then
if (k.gt.m) then
k=1
a=a/2.0
end if
goto 10
else
k=1
do 30 i=1,n
30 x(i)=y(i)
z=z1
if (z.ge.eps) goto 10
end if
end if
return
end subroutine
real function nrnd1(r)
double precision s,u,v,r
s=65536.0
u=2053.0
v=13849.0
m=r/s
r=r-m*s
r=u*r+v
m=r/s
r=r-m*s
nrnd1=r/s
return
end 我写了一个代码,主要是由一个初始点和初始方向确定的直线和一个圆相交,得到的下一个点继续作为初始点,又有一个新的方向作为新方向,继续和圆相交,一直循环。循环的次数为ntrap.我运行不知道为什么没有执行循环,输出in(1),in(2)也没有执行。
x(1)为x,x(2)为y,c为x初坐标,d为y初坐标,in(1)为直线斜率的x坐标,in(2)为直线斜率的y坐标。直线方程为:(x(1)-c)/in(1)=(x(2)-c)/in(2)
我想实现在自己设定的function f 中的o,c,d是从主程序里的o,c,d得来的,每次循环产生新的o,c,d,传递给function f,用来计算,不知道怎么实现。用全局变量不知道写的对不对。程序运行报错如图片所示,还请大神指点迷津!
|
|