fortran代码运行出错请教

夏沫有约

[Fortran] 纯文本查看 复制代码

 program main
           dimension in(2),x(2)
           data ntrap /20/
           integer i
       real*8 in,c,d
           in(1)=-sqrt(2.0)/2.0
           in(2)=sqrt(2.0)/2.0
           c=0.0
           d=-1.0
        do i=1,ntrap,1
           call external f()
                   c=x(1)
                   d=x(2)
           in(1)=in(1)-2.0*c/1.0*vdotn(in(1),x(1),1) ！新方向的x坐标表达式
           in(2)=in(2)-2.0*d/1.0*vdotn(in(1),x(1),1) ！新方向的y坐标表达式
                   write(*,*) 'x(1)=',x(1),'x(2)=',x(2)
                   write(*,*) 'in(1)=',in(1),'in(2)=',in(2)
            end do
          stop
          end 
           
          subroutine external f()
          dimension x(2),y(2)
          data x/2*0.0/
          b=2.0
          n=2
          m=10
          eps=1.0e-5
          call dnmtc(x,n,b,m,eps,y)
          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),in(n)
          f1=x(1)*x(1)+x(2)*x(2)-1.0
          f2=in(2)*x(1)-in(2)*c-in(1)*x(2)+in(1)*d
          f=sqrt(f1*f1+f2*f2)
          end

          subroutine dnmtc(x,n,b,m,eps,y)
          dimension x(n),y(n)
          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)
新人，还不太懂，代码可能写的有问题，没有报错，但结果不对，还请大神们指点！感谢感谢！

夏沫有约

我是在Windows下的Linux系统操作的，可以运行出来，得到的只有x(1),x(2)的值，感觉循环那块没有起作用，不知道要怎么改。

li913

代码语法错误，根本就不能执行。

[Fortran] 纯文本查看 复制代码

program main
           dimension in(2),x(2)
           data ntrap /20/
           integer i
       real*8 in,c,d
           in(1)=-sqrt(2.0)/2.0
           in(2)=sqrt(2.0)/2.0
           c=0.0
           d=-1.0
        do i=1,ntrap,1
           !call external f()
           call  f1()
                   c=x(1)
                   d=x(2)
           in(1)=in(1)-2.0*c/1.0*vdotn(in(1),x(1),1) !新方向的x坐标表达式
           in(2)=in(2)-2.0*d/1.0*vdotn(in(1),x(1),1) !新方向的y坐标表达式
                   write(*,*) 'x(1)=',x(1),'x(2)=',x(2)
                   write(*,*) 'in(1)=',in(1),'in(2)=',in(2)
            end do
          stop
          end 
           
          subroutine  f1()
          dimension x(2),y(2)
          data x/2*0.0/
          b=2.0
          n=2
          m=10
          eps=1.0e-5
          call dnmtc(x,n,b,m,eps,y)
          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),in(n)
          f1=x(1)*x(1)+x(2)*x(2)-1.0
          f2=in(2)*x(1)-in(2)*c-in(1)*x(2)+in(1)*d
          f=sqrt(f1*f1+f2*f2)
          end

          subroutine dnmtc(x,n,b,m,eps,y)
          dimension x(n),y(n)
          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