1             subroutine metod(nnod,x,y,value,ktype)

2             include 'type'

3             integer ktype(NDIM)

4             real*8 x(NDIM),y(NDIM),value(NDIM)

5             real*8 pi,angle,angle1,angle2

6             complex*16 z(NDIM),v_a(NDIM)

7             complex*16 iq,tmp,tmp2,v1,v2

8             iq=(0,1)

9             pi=datan(1.d0)*4.d0

10            do i=1,nnod

11              z(i)=dcmplx(x(i),y(i))

12            end do

13            do i=1,i_shag

14              c(i)=0

15              do j=1,nnod

16                b(i,j)=0

17              end do

18            end do

19            do i=1,i_shag

20              do v_i=1,nnod

21                v_a(v_i)=0

22              end do

23              kt1=i+i_step-1

24              if (kt1.eq.0) kt1=nnod

25              kt2=i+i_step+1

26              if (kt2.gt.nnod) kt2=1

27              tmp2=cdlog((z(kt2)-z(i+i_step))/(z(kt1)-z(i+i_step)))

28              v1=z(kt1)-z(i+i_step)

29              v2=z(kt2)-z(i+i_step)

30              call cauch5(dreal(v1),dimag(v1),angle1)

31              call cauch5(dreal(v2),dimag(v2),angle2)

32              angle=angle1-angle2

33              if (angle.lt.0) angle=angle+2.d0*pi

34              v_a(i+i_step)=dcmplx(dreal(tmp2),-angle)

35              v_a(kt2)=v_a(kt2)+1

36              v_a(kt1)=v_a(kt1)-1

37              do j=1,nnod-2

38                m=i+i_step+j

39                if(m.gt.nnod) m=m-nnod

40                n=m+1

41                if(n.gt.nnod) n=n-nnod

42                tmp=(cdlog((z(n)-z(i+i_step))/

43            *       (z(m)-z(i+i_step))))/(z(n)-z(m))

44                v_a(m)=v_a(m)-1-(z(i+i_step)-z(n))*tmp

45                v_a(n)=v_a(n)+1+(z(i+i_step)-z(m))*tmp

46              end do

47              if (ktype(i+i_step).eq.0) then

48                b(i,i+i_step)=-dimag(v_a(i+i_step))

49                c(i)=c(i)-dreal(v_a(i+i_step))*value(i+i_step)

50                do j=1,nnod-1

51                  m=i+j+i_step

52                  if(m.gt.nnod) m=m-nnod

53                  if(ktype(m).eq.0) then

54                    b(i,m)=-dimag(v_a(m))

55                    c(i)=c(i)-dreal(v_a(m))*value(m)

56                  else

57                    b(i,m)=dreal(v_a(m))

58                    c(i)=c(i)+dimag(v_a(m))*value(m)

59                  end if

60                end do

61              else

62                b(i,i+i_step)=dimag(v_a(i+i_step))

63                c(i)=c(i)-dreal(v_a(i+i_step))*value(i+i_step)

64                do j=1,nnod-1

65                  m=i+j+i_step

66                  if(m.gt.nnod) m=m-nnod

67                  if(ktype(m).eq.0) then

68                    b(i,m)=dreal(v_a(m))

69                    c(i)=c(i)-dimag(v_a(m))*value(m)

70                  else

71                    b(i,m)=dimag(v_a(m))

72                    c(i)=c(i)-dreal(v_a(m))*value(m)

73                  end if

74                 end do

75              end if

76            end do

77            return

78            end

79        

80            subroutine cauch5(x,y,angle)

81            implicit double precision (A-H,O-Z)

82            pi=datan(1.d0)*4.d0

83            if (x.eq.0.d0 .and. y.gt.0.d0) angle=pi/2.d0

84            if (x.eq.0.d0 .and. y.lt.0.d0) angle=3*pi/2.d0

85            if (x.gt.0.d0 .and. y.ge.0.d0) angle=datan(y/x)

86            if (x.lt.0.d0 .and. y.ge.0.d0) angle=pi-datan(-y/x)

87            if (x.lt.0.d0 .and. y.lt.0.d0) angle=pi+datan(y/x)

88            if (x.gt.0.d0 .and. y.lt.0.d0) angle=2*pi-datan(-y/x)

89            return

90            end

Рис. 67. Листинг подпрограммы для реализации алгоритма формирования СЛАУ (общий для стационарной и нестационарной задачи)