! A - матрица коэффициентов
! B - вектор-столбец свободных членов
! X - вектор-столбец искомых 
! n -  порядок системы
          DIMENSION A(7,7),B(7),X(7), Be(7),He(7), Ha(7) 
          integer i,j,n
! SLAU3.dat - файл исходных данных
! results - файл результатов
       open(unit=5,file="SLAU3.dat")
       open(unit=12,file="result.txt")
          read(5,*) n,((a(i,j),j=1,n),i=1,n)
           read(5,*)    (b(i),i=1,n) 
            write(12,*) "Исходная система ЛАУ (4.52)"
           write(12,10) ((a(i,j),j=1,n),i=1,n)
           write(12,10) (b(j),j=1,n)
         call Prog(a,b,n,x)
          write(12,*) "Преобразованная система ЛАУ (4.52)"
         write(12,10) ((a(i,j),j=1,n),i=1,n)
         write(12,10) (b(j),j=1,n)
         write(12,*) "Искомые корни системы ЛАУ (4.52):"
          write(12,*) "коэффициенты кубического сплайна"
         write(12,10) (x(j),j=1,n)
 !     Ввод значений индукции Be и напряжённости He в узлах
!       на экспериментальной КН марки 1511 (см. табл. 4.4)
                       read(5,*) (Be(i),He(i), i= 1,n)
           write(12,*)  "Значения напряженности Ha в узлах КН, рассчитанные по (4.32)"
	     do i=2,n
        bi =Be(i-1) 
        h1 =x(i-1)*(Be(i)-bi)*(Be(i)-bi)*(Be(i)-bi)/1.8 
        h2 =x(i)  *(bi-Be(i-1))*(bi-Be(i-1))*(bi-Be(i-1))/1.8 
        h3 =(He(i-1)-(x(i-1)*0.09/6.0))*((Be(i)-bi)/0.3) 
        h4 =(He(i)-(x(i)*0.09/6.0))*((bi-Be(i-1))/0.3) 
        Ha(i) =h1+h2+h3+h4;
         end do
       write(12,10) (Ha(i), i = 2,n)
         stop
10       format(7(1x,f10.2))
         end
         
          subroutine Prog(A,B,N,X)
! в подпрограмме Prog - реализован иетод прогонки
   DIMENSION a(N,N),b(N),x(N)
   real f              
   integer i,im1,ip1,nm1   
! прямой ход - вычисление прогоночных коэффициентов
              do i=2,n
              im1=i-1  
               f=a(i,im1)/a(im1,im1) 
              a(i,i)=a(i,i)-f*a(im1,i)               
              b(i)=b(i)-f*b(im1)
              end do
! обратных ход - определение неизвестных
         x(n)=b(n)/a(n,n)
        nm1=n-1
        do i=nm1,1,-1  
          ip1=i+1  
        x(i)=(b(i)-x(ip1)*a(i,ip1))/a(i,i)
        end do
  return 
 end
