!Listing 4.1. Начало 
!1) головной модуль NR
        external SG13, SPLAIN, rkf4
  ! Head module      
        DIMENSION BEXPER(42),HEXPER(42), ASPLA(42), DIFKUR(42), p(4), Rsi(42), COFs(42),Reri(42), COFer(42), w(10)
      integer i, m, nexper, n
      real  LA1,LB1,LC1,La2,Lb2,Lc2,Lan,Lbn,Lcn
    !  double precision UAB, UBC, ULMAB, ULMBC,Ua2x,Ub2x,Uc2x 
      COMMON /EU1/ ULMAB, ULMBC,UAB,UBC,Ua2x,Ub2x,Uc2x, Ea1 &
        /EU1a/ UfAm,UfBm,UfCm, UfA,UfB,UfC, UfAmm, UfBmm, UfCmm, UfAm1,UfBm1,UfCm1 &
                /EU2/ WA1,WB1,WC1,Wa2,Wb2,Wc2,Wa2x,Wb2x,Wc2x  &
            /EU3/ RA1,RB1,RC1,Ra2,Rb2,Rc2,Ran,Rbn,Rcn , Rd   &
            /EU4/ DLs,DLer, Ss,Ser,ks   &
            /EU5/ UF,ii   &
            /EU6/ NEXPER,BEXPER,HEXPER   &
            /EU7/ Rsi,COFs, Reri, COFer    &
            /EU8/ Bst,fi,bst3  &
            /EU10/  LA1,LB1,LC1,La2,Lb2,Lc2,Lan,Lbn,Lcn
              					open(unit=5,file="TR.dat")
      						open(unit=6,file="dat.txt")
                               open(unit=10,file="B-1-2-3-xx.dat")
                                 open(unit=12,file="Fi_1-xx.dat")
                                    open(unit=16,file="U_1-xx.dat")
                ks = 4 ! 4 loop magnetic fluxes 
                NEXPER=42
       read(5,*) (BEXPER(I),HEXPER(I),I=1,NEXPER)
 !   write(6,*) "B=f(H) - experiment"
  !     write(6,*) (BEXPER(I),HEXPER(I),I=1,NEXPER)
  ! External procedure SG13 - smoothing the experimental magnetization curve
 CALL SG13(BEXPER,HEXPER,NEXPER)
  !    write(6,*) "B=f(h), sg13"
 !      write(6,*) (BEXPER(I),HEXPER(I),I=1,NEXPER)

                                                BEXPER(1)=0.
                                                HEXPER(1)=0.
  ! External procedure SPLAIN - for determining the coefficients of a cubic spline                                             
      CALL SPLAIN(BEXPER,HEXPER,NEXPER,ASPLA)
  !                  write(6,*) "Aspa - splain"
 ! write(6,*) (ASPLA(M),M=1,NEXPER)
                 
       !!                   4     FORMAT('aspla x4',12G12.4)
                      NEXPM1=NEXPER-1  
 DO  M=1,NEXPER
              DO  I=1,NEXPM1
      IF(BEXPER(M)>BEXPER(I+1)) GOTO 101
                                             NL=I
          GOTO 120
 101 CONTINUE
 Enddo
 120  CONTINUE
                                             NR=NL+1
      DB =BEXPER(NR)-BEXPER(NL)
      DB2=0.5/DB
      BBT=BEXPER(NR)-BEXPER(M)
      BTB=BEXPER(M) -BEXPER(NL)
      DH =HEXPER(NR)-HEXPER(NL)
      DASP=ASPLA(NR)-ASPLA(NL)
      DIFKUR(M)=(-ASPLA(NL)*(BBT*BBT))*DB2+ (ASPLA(NR)*(BTB*BTB))*DB2+(DH)/DB-((DASP)*0.16666)*DB
 Enddo
    
 !write(6,*) "B, dH/dB, H"
 ! write(6,*) (bexper(m),DIFKUR(M),hexper(m),M=1,NEXPER)
   
       READ(5,*) n, UfA,UfB,Ufc,  WA1,WB1,WC1,Wa2,Wb2,Wc2,  RA1,RB1,RC1,Ra2,Rb2,Rc2,Ran,Rbn,Rcn, &
                 Rd,   LA1,LB1,LC1,La2,Lb2,Lc2,Lan,Lbn,Lcn, DLs,DLer, Ss,Ser, F 

      write(6,*) "n, UfA,UfB,Ufc,  WA1,WB1,WC1,Wa2,Wb2,Wc2,  RA1,RB1,RC1,Ra2,Rb2,Rc2,Ran,Rbn,Rcn" 
      write(6,*)     " Rd,   LA1,LB1,LC1,La2,Lb2,Lc2,Lan,Lbn,Lcn, DLs,DLer, Ss,Ser, F "
           write(6,*) n, UfA,UfB,Ufc,  WA1,WB1,WC1,Wa2,Wb2,Wc2,  RA1,RB1,RC1,Ra2,Rb2,Rc2,Ran,Rbn,Rcn, &
      				  Rd,   LA1,LB1,LC1,La2,Lb2,Lc2,Lan,Lbn,Lcn, DLs,DLer, Ss,Ser, F 
      
        read(5,*) (p(i),i=1,4)
  				  write(6,*) "t_bedin, t_end, dt_initial step, accuracy"
                  write(6,*)(P(I),I=1,4)
  
!	read(5,*) Wa2x,Wb2x,Wc2x
!	write(6,*) "Wa2x,Wb2x,Wc2x"
!	write(6,*)Wa2x,Wb2x,Wc2x   
                  
  !    54                 FORMAT(' PR(I)=',4G13.5)
  !     10  FORMAT (i3,2x,2G13.4/  3G13.4  /  3G13.4/  3G13.4/ 2G13.4/ G13.4)
                 
               DO   I=1,NEXPER
    Rsi(I)=DIFKUR(I)*DLs/Ss
    Reri(I)=DIFKUR(I)*Dler/Ser
Enddo
                    !   write(6,*) "Rst_i, i"
                     !                 write(6,*)(RSTI(I),I=1,NEXPER)
                                     
               DO   KK=1,NEXPM1
      BTEKT=(BEXPER(KK+1)-BEXPER(KK))
      BTEKT1=1./BTEKT
      COFs(KK)=(Rsi(KK+1)-Rsi(KK))*BTEKT1
      COFer(KK)=(Reri(KK+1)-Reri(KK))*BTEKT1
Enddo
   !   ULMAB=SQRT(2.)*UAB
    !  ULMBC=SQRT(2.)*UBC
        UfAm1=SQRT(2.)*UfA
        UfBm1=SQRT(2.)*UfB
        UfCm1=SQRT(2.)*UfC

      UF=2.0*3.14159*F
DO   I=1,n
                   w(I)=0.0
Enddo
!External procedure RKf4(Runge-Kutta method)- determination of instantaneous values
! of currents and magnetic fluxes
      CALL RKf4(P,w,N,f1)
            Write(*,*) 
      Write(*,*) "The program has finished work"
      STOP
  END
!2) модуль SG13,
         subroutine sg13(b,h,kt)
! Smoothing the experimental magnetization curve of electrical steel  
         DIMENSION H(kt),B(kt),HR(kt)
         DO  J=0,KT-3
         K=J+2
         B3=0.
         H3=0.
         DO  I=1+J,3+J
         B3=B3+B(I)
       H3=H3+H(I)
 End do
         B3S=B3/3.
         H3S=H3/3.
         CH=0.
         ZN=0.
         DO  I=1+J,3+J
         CH=CH+(B(I)-B3S)*(H(I)-H3S)
        ZN=ZN+(B(I)-B3S)**2
 End do
         A=CH/ZN
         HR(K)=A*(B(K)-B3S)+H3S
         IF(J==0)HR(1)=A*(B(1)-B3S)+H3S
         IF(J==(KT-3))HR(KT)=A*(B(KT)-B3S)+H3S
End do
			do  ii=1,kt
 	 h(ii)=hr(ii)
 End do
         return
         END
!3) модуль SPLAIN,
      SUBROUTINE SPLAIN(B,H,N,Z)
      ! Определение коэффициентов кубического сплайна
  !    external gaus
      integer i, j, n
    DIMENSION B(N),H(N),A(N,N),Z(N),DB(100),DH(100),C(2), x(n)
!write(*,*) "splain-tect simq n = ", n
!write(6,*) "splain-tect simq n= ", n
    DO 11 I=1,N
      DO 10 J=1,N
  10  A(I,J)=0.0
  11 continue
      NM1=N-1
      DO 20 I=2,N
      DB(I)=B(I)-B(I-1)
  20  DH(I)=H(I)-H(I-1)
      DO 50 I=2,NM1
      A(I,I-1)=DB(I)/6.0
      A(I,I+1)=DB(I+1)/6.0
      A(I,I)=(DB(I)+DB(I+1))/3.0
      Z(I)=DH(I+1)/DB(I+1)-DH(I)/DB(I)
  50  CONTINUE
      BS=(B(1)+B(2))/2.0
      DO  65 I=1,2
  65  C(I)=((BS-B(I))/DB(2))**2
      A(1,1)=3.0*(1.0-C(2))
      A(1,1)=3.0*(1.0-C(2))
      A(1,2)=3.0*C(1)
!     Z(1)=DH(2)/DB(2)*(6.0/DB(2)-1.0)
      Z(1)=0.
      BS=(B(N-1)+B(N))/2.0
      DO 75 I=1,2
  75  C(I)=((BS-B(N-2+I))/DB(N))**2
      A(N,N-1)=3.0*C(2)
      A(N,N)=3.0*(1.0-C(1))
!     Z(N)=DH(N)/DB(N)*(1.0-6.0/DB(N))
      Z(N)=0.
     !                          write(*,1) (z(M),M=1,n)
     ! Метод прогонки: расчёт коэффициентов кубического сплайна
      CALL Prog(A,Z,N,X)
       !                         write(*,2) (z(M),M=1,n)
            !                     write(*,3) (x(M),M=1,n)
                               
                     !   1     FORMAT(' z1',12G12.4)      
                      !       2     FORMAT(' z2',12G12.4)   
                           !       3     FORMAT(' x3',12G12.4)  
      RETURN

   CONTAINS 
! Решение системы ЛАУ с трёхдиагональной матрицей коэффициентов.
         subroutine Prog(A,B,N,X)
   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  
          ip1=i+1  
        x(i)=(b(i)-x(ip1)*a(i,ip1))/a(i,i)
        end do
 
 end SUBROUTINE
  END
!4) модуль GAUS2,
        subroutine GAUS2(A,B,N,X)
! External procedure GAUS2(Gaussian method) - determination of derivatives of currents and magnetic fluxes
              DIMENSION a(N,N),b(N),x(N)
   !                                     write(6,*) "Gaus2-test1 rkf4 N=" , N      
   !                                     write(*,*) "Gaus2-test1 rkf4 N=" , N   
              NM1=N-1
              DO 303 K=1,NM1
!          aa
              KP1=K+1
              L=K
              DO 400 I=KP1,N
              IF(ABS(A(I,K))-ABS(A(L,K))) 400,400,401
 401          L=I
 400          CONTINUE
              IF(L-K) 500,500,405
 405          DO 410 J=K,N
              TEMP=A(K,J)
              A(K,J)=A(L,J)
 410          A(L,J)=TEMP
  !                                      write(6,*) "Gaus-test2"       

              TEMP=B(K)
              B(K)=B(L)
              B(L)=TEMP
!           xx
!           yy
 500                    DO 300  I=KP1,N
                        F=A(I,K)/A(K,K)
                                      !  write(6,*) "Gaus-test3", F, i, k, a(i,k), a(k,k)
                        A(I,K)=0.
                                      !  write(6,*) "Gaus-test4"       
                        DO 301 J=KP1,N
                                       ! write(6,*) "Gaus-test5", F ,i, j, A(I,J), k, a(k,j)
                    A(I,J)=A(I,J)-F*A(K,J)
   301           continue
                                      !  write(6,*) "Gaus-test6", F ,i, j, A(I,J)     
                    B(I)=B(I)-F*B(K)
                                      !  write(6,*) "Gaus-test7", F, i, b(i)      
   300               continue
                                       ! write(6,*) "Gaus-test8", F, i, b(i)  
     303               continue
                                       ! write(6,*) "Gaus-test9", F, i, b(i)                                         

!           zz
                              X(N)=B(N)/A(N,N)
                              I=NM1
 710                          IP1=I+1
                              S=0
                              DO 700 J=IP1,N
 700                          S=S+A(I,J)*X(J)
                              X(I)=(B(I)-S)/A(I,I)
                              I=I-1
                             IF(I) 800,800,710
    800          RETURN
                  END
!5) модуль RKF4.
        SUBROUTINE RKf4(P,W,N,f)
    !    external fct, outp
    !Implementation of the Runge-Kutta method with Felberg's modification:
 !   determination of the required quantities - instantaneous values 
 !   of currents and magnetic fluxes as a function of time
            DIMENSION y(n), w(n), a(n),b(n),c(n),d(n),e(n),u(n), f(n), P(5)
            integer e3
  !           write(*,*) "RKtest1"  , n     
            X =P(1) ! begin X = 0 c
            XK=P(2) ! end X = 0.2 c
            H =P(3) ! initial value of step h = 0.0002
            E1=P(4) ! specified accuracy -  0.00001
            do 13 i=1,n
 13         Y(i)=W(i)
		C1=1./9.
		C2=1./64.
		C3=1./2.4
		C4=1./15.
		C5=1./144.
		C6=1./6.
        C7=1./45.
        C8=1./75.
		C9=1./150.
        C10=1./12.
        C11=1./432.
		C13=1./27
 !                    write(*,*) "RKtest2"       
            CALL OUTP(X)
 !                        write(*,*) "RKtest3"       
  110           e3 = 0           
!          write(*,*) "RKtest4"         
            CALL FCT(X,Y,F,N)
            N1=0
            DO 1 J=1,N
            A(J)=F(J)*H
 1          Y(J)=W(J)+2.*A(J)*c1
            X=X+2.*H*c1
 !                        write(*,*) "RKtest5"       
            CALL  FCT(X,Y,F,N)
            DO 2 J=1,N
            B(J)=H*F(J)
 2          Y(J)=W(J)+A(J)*c10+b(J)*0.25
            x=x+h*c1
            CALL FCT(X,Y,F,N)
            DO 3 J=1,N
            C(J)=F(J)*H
 3          Y(J)=W(J)+(34.5*a(j)-121.5*b(j)+135.*c(j))*c2
            X=X+H*c3
            CALL FCT(X,Y,F,N)
            DO 4 J=1,N
            D(J)=F(J)*H
            Q=W(J)-17.*A(J)*C10+6.75*B(J)
 4          Y(J)=Q-5.4*C(J)+16.*D(J)*C4
            X=X+H*0.25
            CALL FCT(X,Y,F,N)
            DO 5 J=1,N
            E(J)=F(J)*H
            Q=W(J)+65.*A(J)*C11-0.3125*B(J)
 5           Y(J)=Q+0.8125*C(J)+4.*D(J)*C13+5.*E(J)*C5
            X=X-H*C6
            CALL FCT(X,Y,F,N)
            DO 260 J=1,N
            U(J)=H*F(J)
            Y(J)=W(J)+A(J)*C1+0.45*C(J)+16.*D(J)*C7+E(J)*C10
                     Q=A(J)*C9-0.03*C(J)+16.*D(J)*C8+E(J)*0.05
            E2=ABS(Q-0.24*U(J))
            IF(E2<=E1) GOTO 250
            E3=1
            GOTO 260
 250        IF(E2<(E1*0.05))  N1=N1+1
  260        CONTINUE
            X=X+H*C6
            CALL FCT(X,Y,F,N)
            IF(E3==0)  GOTO  290
            X=X-H
            DO 6 J=1,N
 6          Y(J)=W(J)
            H=H*0.5
            GOTO  110
  290        IF(N1==n)  H=H+H
            if(h>p(3))  h=p(3)
            CALL OUTP(X)
            do 15 j=1,n
 15         W(j)=Y(j)
            IF(X>=XK)  GOTO 7
  goto  110        
  7         RETURN
 CONTAINS 
 
   SUBROUTINE OUTP(X)
   !Procedure OUTP - Output of results: instantaneous values
   ! of currents, voltages, magnetic inductions  
      DIMENSION Bst(3),fi(10),bst3(3)
      COMMON /EU1/ ULMAB, ULMBC,UAB,UBC,Ua2x,Ub2x,Uc2x,Ea1 &
         /EU1a/ UfAm,UfBm,UfCm, UfA,UfB,UfC, UfAmm, UfBmm, UfCmm,UfAm1,UfBm1,UfCm1 &    
      /EU2/ WA1,WB1,WC1,Wa2,Wb2,Wc2,Wa2x,Wb2x,Wc2x   &
             /EU8/ Bst,fi,bst3
 !           WRITE(*,*) "RKtest2a"
!                      IF(ABS(X-A)-0.000001) 8,8,6
   !   WRITE(7,99) X,(bst(I),I=1,3)
   !   WRITE(8,*) X,(bst3(I),I=1,3)
  !    WRITE(9,*) X,bst9
          
      WRITE(10,*) X,bst(1),bst(2),bst(3)
      WRITE(12,*) X,fi(5),fi(7) ,fi(9)
   !   WRITE(14,*) X,fi(6),fi(8),fi(10)
   !   WRITE(15,*) X,Ua2x,Ub2x,Uc2x
      WRITE(16,*) X,UfAmm, UfBmm, UfCmm
   !    WRITE(17,*) X,Ea1
          RETURN
 !   99   FORMAT(10G16.6)
  !  5 FORMAT(T2,F9.7,1X,9(F6.3,1X))
  !  7 FORMAT(T2,F9.7,1X,3(F9.4,1X))
!    9 FORMAT('+',T35,3(F8.2,1X))
 !  20     FORMAT(' ',F9.7,1X,-4P4F9.1)
!   30     FORMAT('+',T55,4(F7.3,1X))
          END SUBROUTINE
          
      SUBROUTINE FCT(X,Y,ff,n)
 ! Procedure FCT - formation of a matrix of coefficients and a column vector of free members
external gaus2
         REAL       LA1,LB1,LC1,La2,Lb2,Lc2,Lan,Lbn,Lcn
      DIMENSION Y(10),FI(10),PC(4),BST(3),BerV(3),BerN(3), B(10),bst3(3),&
               BEXPER(42),HEXPER(42),Rsi(42), COFs(42),Reri(42), COFer(42), &
               A(10,10),UMGN(2),BTEK(10),ff(10),x3(10),Rs(3),RerV(3),RerN(3)
      COMMON /EU1/ ULMAB, ULMBC,UAB,UBC,Ua2x,Ub2x,Uc2x,Ea1 &
         /EU1a/ UfAm,UfBm,UfCm, UfA,UfB,UfC, UfAmm, UfBmm, UfCmm,UfAm1,UfBm1,UfCm1 &
             /EU2/ WA1,WB1,WC1,Wa2,Wb2,Wc2,Wa2x,Wb2x,Wc2x   &
            /EU3/ RA1,RB1,RC1,Ra2,Rb2,Rc2,Ran,Rbn,Rcn,Rd   &
            /EU4/ DLs,DLer, Ss,Ser,ks  &
            /EU5/ UF,ii   &
            /EU6/ NEXPER,BEXPER,HEXPER   &
            /EU7/ Rsi,COFs, Reri, COFer    &
            /EU8/ Bst,fi,bst3   &
            /EU10/  LA1,LB1,LC1,La2,Lb2,Lc2,Lan,Lbn,Lcn   
        APRLIN(X2,X1,Y1,X2Y2)=Y1+(X2Y2)*(X2-X1)
 !                                write(*,*) "FCTtest1" , n      

        ii=1
      T=X
    !                                    write(*,666) T,(Y(I),I=1,n)
               DO 2001 I=ks+1,n
 2001 FI(I)=Y(I)
               DO 3001 I=1,ks
 3001                      PC(I)=Y(I)
!  Ua2x = PC(1)* Wa2x * 0.001
!  Ub2x = (PC(2)-PC(1))* Wb2x * 0.001
 ! Uc2x = -PC(2)* Wc2x * 0.001
 ! Ea1 = PC(1)* Wa1 * 0.0001

 	!			do 1 i=1,11              
                 Bst(1)= PC(1)/Ss
                 Bst(2)=(PC(2)-PC(1))/Ss
                 Bst(3)= -PC(2)/Ss
             BerV(1)=(PC(1)+PC(3))/Ser
             BerV(2)=(PC(2)+PC(3))/Ser
             BerV(3)= PC(3)/Ser
             BerN(1)=(PC(1)+PC(4))/Ser
             BerN(2)=(PC(2)+PC(4))/Ser
             BerN(3)= PC(4)/Ser
 				do 1 i=1,3              
                   Btek(i)=Bst(i)
                   Btek(i+3)=BerV(i)
 1                 Btek(i+6)=BerN(i)

 			nn=0              
33        continue
						 DO 3 K=1,3
                         knn=k+nn
                          if(BTEK(Knn)>=0.) GOTO 15
                 BTEK(Knn)=ABS(BTEK(Knn))
   15 CONTINUE
                             NNEXP=21
      IF(BTEK(Knn)>BEXPER(NNEXP)) GOTO 300
                           NL=1
                           NR=11
                           ND=10
      GOTO 299
  300 CONTINUE
                           NL=22
                           NR=32
                           ND=10
      IF(BTEK(Knn)<BEXPER(22)) GOTO 38
  299 CONTINUE
                   DO 301 I=NL,NR,ND
      IF(BTEK(Knn)>BEXPER(I+ND)) GOTO 301
                                               NL=I
      GOTO 302
  301 CONTINUE
                           NL=41
!                          NR=42
          GOTO 295
  302 CONTINUE
                           NR=NL+ND
               IF(ND==5) GOTO 296
               IF(ND==1) GOTO 295
                           ND=5
                           NR=NR-5
      GOTO 299
  296 CONTINUE
                           ND=1
                           NR=NR-1
      GOTO 299
  295 CONTINUE
                           goto 19
38                         nL=21
19      continue
					if(nn.eq.3) goto 101
          			if(nn.eq.6) goto 102
          Rs(K)=APRLIN(BTEK(K),BEXPER(NL),Rsi(NL),COFs(NL))
       !           write(*,779) NL,BEXPER(NL),K,BTEK(K),Rsi(NL),RST(K)
!                            write(1,777) NL,BEXPER(NL),k,BTEK(K)
 					if(nn.eq.0) goto 3
101      continue
         RerV(K) =APRLIN(BTEK(Knn),BEXPER(NL),Reri(NL),COFer(NL))
          if(nn.eq.3) goto 3
102      RerN(K) =APRLIN(BTEK(Knn),BEXPER(NL),Reri(NL),COFer(NL))
3        continue
          if(nn.eq.3) goto 202
          if(nn.eq.6) goto 203
          nn=3
          goto 33
202       nn=6
          goto 33
203       continue
!    3 CONTINUE
                DO 18 I=1,n
                DO 18 J=1,n
   18             A(I,J)=0.
!                  do 58 j=1,ks
!58                  a(j,j)=Rst(j)
!----------------
     !                if(t>0.1) Lan= 0.00032
        !             if(t>0.1) Lbn=3200.0
           !          if(t>0.1) Lcn=3200.0  
                    
      !     Wa2 = 0.0
          ! Wb2 = 0.0
         !  Wc2 = 0.0
 !                     if(t>=0.1) Wa2=60.0
  !                    if(t>=0.1) Wb2=60.0
   !                   if(t>=0.1) Wc2=60.0
 A(1,1)= Rs(1)+Rs(2)+RerV(1)+RerN(1)+4*Rd
  A(1,2)=-(Rs(2)+2*Rd)
  A(1,3)= RerV(1)
  A(1,4)= RerN(1)
  A(1,5)= WA1
  A(1,6)=-Wa2
  A(1,7)=-WB1
  A(1,8)= Wb2
  	 A(2,1)=-(Rs(2)+2*Rd)
     A(2,2)= Rs(2)+Rs(3)+RerV(2)+RerN(2)+4*Rd
     A(2,3)= RerV(2)
     A(2,4)= RerN(2)
 	 A(2,7)= WB1
 	 A(2,8)=-Wb2
   	 A(2,9)=-WC1
  	 A(2,10)=Wc2
		A(3,1)=  RerV(1)
                                A(3,2)=  RerV(2)  		
		A(3,3)=  RerV(1)+ RerV(2)+ RerV(3)
			A(4,1)= RerN(1)
      	                                A(4,2)= RerN(2)  		
	  		A(4,4)= RerN(1)+RerN(2)+RerN(3)
        		        A(5,1)=-(WA1+WB1)
   		        A(5,2)= WB1
 		        A(5,5)= LA1
 		        A(5,7)=-LB1
  					A(6,1)= Wa2+Wb2
  					A(6,2)=-Wb2
 		                                                A(6,6)=  La2+Lan
 		                                                A(6,8)=-(Lb2+Lbn)
  						A(7,1)=  WB1
  						A(7,2)=-(WB1+WC1)
 		       		                                 A(7,7)= LB1
 		      		                                 A(7,9)=-LC1
  							A(8,1)= -Wb2
  							A(8,2)=  Wb2+Wc2
 		            		                                                A(8,8)=   Lb2+Lbn
 		       			                                A(8,10)=-(Lc2+Lcn)
                                 DO 77 L=5,9,2
                                   A(9,L)= 1.
 77                               A(10,L+1)= 1.
!-----------------------------------
                   !     write(*,111)ii
        !      if(ii.le.9)    WRITE(*,88)((A(I,J),J=1,n),I=1,n)
            DO 50 I=1,n
  50           b(I)=0.0
!  xx
                ! if(t>=0.04) ulmab=sqrt(2.)*Uab*1.03
                ! if(t>=0.04) ulmbc=sqrt(2.)*Ubc*1.03
                ! yy
                
                 !   if(t<0.06) UfAm=UfAm1
                  !  if(t<0.06) UfBm=UfBm1
                   ! if(t<0.06) UfCm=UfCm1
                  !  if(t>=0.06) UfAm=UfAm1*1.1
                  !  if(t>=0.06) UfBm=UfBm1*1.1
                  !  if(t>=0.06) UfCm=UfCm1*1.1
                 !   if(t>0.14) UfAm=UfAm1
                 !   if(t>0.14) UfBm=UfBm1
                 !  if(t>0.14) UfCm=UfCm1
                 
                  !   if(t<0.06) Ran=50000.0
                  !   if(t<0.06) Rbn=50000.0
                  !   if(t<0.06) Rcn=50000.0
                  !   if(t>=0.06) Ran=50000.0*1.1
                   !  if(t>=0.06) Rbn=50000.0*1.1
                  !   if(t>=0.06) Rcn=50000.0*1.1
                  !   if(t>0.14) Ran=50000.0
                 !    if(t>0.14) Rbn=50000.0
                  !   if(t>0.14) Rcn=50000.0
! zz
             !    if(t>=0.14) ulmab=sqrt(2.)*UAB
              !   if(t>=0.14) ulmbc=sqrt(2.)*UBC
 ! kz                        ulmab=sqrt(2.)*UAB
            !     if(t>=0.08) ulmbc=sqrt(2.)*UBC*0
                                             TUFa= UF*T
                                  TUFb=TUFa-2.094395
                                  TUFc=TUFa-4.18879
                !     if(t>=0.1) UfAm1 = sqrt(2.0)*6.3005027E+7
             !                          	 if(t>=0.1) UfBm1 = sqrt(2.0)*5.7735027E+7
				! 	if(t>=0.1) UfCm1 = sqrt(2.0)*5.1005027E+7

                  UfAmm= UfAm1*sin(tufa)
                  UfBmm= UfBm1*sin(tufb)
                  UfCmm= UfCm1*sin(tufc)
               UMGN(1)=(UfAmm-UfBmm)
               UMGN(2)=(UfBmm-UfCmm)

!                                 TUF= UF*T
 !             UMGN(1)=ULMAB*SIN(TUF)
   !                              TUFP=TUF-2.094395
     !         UMGN(2)=ULMBC*SIN(TUFP)
   Rm = 0.    !   Rm =500000.  
                 !    if(t>=0.1) Ran=7.0E+04
                   !  if(t>=0.1) Rbn=5.0E+04
                    ! if(t>=0.1) Rcn=1.0E+04
                     
        b(5)= UMGN(1)-FI(5)*RA1 +FI(7)*RB1 - Rm*fI(5)+Rm*fI(7)
        b(6)=        -FI(6)*Ra2 +FI(8)*Rb2 - Ran*fI(6)+Rbn*fI(8)
        b(7)= UMGN(2)-FI(7)*RB1 +FI(9)*RC1 - Rm*fI(7)+Rm*fI(9)
        b(8)=        -FI(8)*Rb2 +FI(10)*Rc2 - Rbn*fI(8)+Rcn*fI(10)
     !     if(ii.le.9)          write(1,99) (b(I),I=1,n)
 !                        write(*,*) "FCTtest2"       
! n = 10
!External Procedure gaus2 (Gaussian method) - determination of values of derivatives
! of currents and magnetic fluxes
           CALL gaus2(A,B,n,x3)
         !  if(ii.le.9)             write(*,66) (x3(I),I=1,n)
                 !        write(*,*) "FCTtest3"   
!  if(ii.le.9)    WRITE(*,88)((A(I,J),J=1,n),I=1,n)                            

           DO 407 I=1,n
  407     ff(i)   =x3(i)
!               ii=ii+1
      RETURN
!   666                            FORMAT(' T=',F7.5,' K T ',12G11.4)
 ! 779 FORMAT(' NL=',I2,' BEXP=',F5.2,' K=',I2,' BTEK=',F5.2,2X,2(G11.4,2X))
       
 !  777                        FORMAT(' NL=',I2,' BEXP=',F5.2,' K=',I2,' BTEK=',F6.3)
                          
 !  88                      FORMAT(' A ',10G11.4)
 ! 99                                             FORMAT(10G13.6)
! 66                                               FORMAT(' x3',12G12.4)
! 111             format(' ii=',i3)
      END SUBROUTINE
            END
! Listing 4.1. Окончание 

! Файл TR.txt c исходными данными (cм. головной модуль TR, операторы read)
! 0.     0.     0.6    0.81   0.7    1.1   0.8   1.3   0.9   1.52   1.     1.82
! 1.1    2.13   1.2    2.43   1.3    2.75  1.4   3.2   1.45  3.5    1.5    3.9
! 1.55   4.5    1.57   4.78   1.6    5.2   1.63  5.88  1.65  6.32   1.67   6.98
! 1.7    8.     1.73   9.4    1.75  10.4   1.77 12.24  1.8  15.     1.83  19.22
! 1.84  21.44   1.85  23.66   1.86  25.88  1.87 28.2   1.9  38.25   1.91  42.
! 1.92  46.     1.93  52.     1.94  58.    1.95 70.    1.96 82.     1.97  94.
! 1.98 109.     1.99 134.     2.   160.    2.01 200.   2.02 250.    2.03 300.
!10
! 57735027.00  57735027.00  57735027.00  600.  600.  600. 60.  60.  60.
!    70000.00000    70000.00000    70000.00000
! 70.00000    70.00000    70.00000
! 50000.00000    50000.00000    50000.00000 10.
!       63.0000       63.0000       63.0000
! 0.0630000       0.063000       0.0630000
! 32.0000       32.00000       32.0000
! 77. 55. 270. 270.
! 50.
! 0.  0.8  0.0025  0.01
