C     RE277A TRAFICORPORATION
C     TRANSFER FUNCTION EXPRESSED AS A RATIO OF TWO COMPLEX POLYNOMIALS
C     WHERE THE NUMERATOR IS OF ONE DEGREE LESS THAN THE DENOMINATOR
C
      DIMENSION OMEGA(250),PH(250),DBM(250),R(250),UR(250),W(250),AMDA(5
     10),V(50),S(25),T(25),U(50),A(25,25),B(25),EM(250),P(15),Q(15),QABS
     21(250),QABS2(250),QABSQ(250),QABS(250),PABS1(250),PABS2(250),PABSQ
     3(250),PABS(250),CEM(250),CDBM(250),ERORM(250),CPH(250),THETA1(250)
     4,THETA2(250),ERORP(250),ALPA(250),BETA(250),PHE(250),CPHR(250)
C
C     READ INPUT
C
 1005 FORMAT(6E12.5)
 1006 FORMAT(3I6,4E12.5)
 1000 READ 1006,IDNO,NOMEGA,N,CURNER,WALT
      DO 1010 J = 1,NOMEGA
 1010 READ 1005,OMEGA(J),DBM(J),PHE(J)
 1500 READ 1005,WX,SCALE,CONV1,CONV2
      ERMAG1=0.
      ERFAZ1=0.
      MAXPT1=0
      MAXPT3=0
      ITER=1
      DO 1510 J = 1,NOMEGA
 1510 W(J) = WX
C
C     CONVERSION OF EXPERIMENTAL G(J) INTO THE REAL AND IMAGINARY PARTS
C
      DO 500 J = 1,NOMEGA
      PH(J) = PHE(J)*0.0174533
      EM(J) = EXPF(DBM(J)/8.6858896)
      R(J) = EM(J)*COSF(PH(J))
  500 UR(J) = EM(J)*SINF(PH(J))
      N2 = N*2
      QZERO = 1.
      NOLAMD = N*2-1
C
C     COMPUTATION OF LAMBDA,U,T AND S
C
    7 DO 15 J = 1,N2,2
    5 TEMP = 0.
      DO 10 K = 1,NOMEGA
      J1 = J-1
   10 TEMP = OMEGA(K)**J1*W(K)+TEMP
      AMDA(J) = TEMP
   15 AMDA(J+1) = 0.
   20 DO 35 J = 1,N2,2
   25 TEMP = 0.
      DO 30 K = 1,NOMEGA
      J1 = J+1
   30 TEMP = OMEGA(K)**(J1)*(R(K)**2+UR(K)**2)*W(K)+TEMP
      U(J) = TEMP
   35 U(J+1) = 0.
   40 DO 55 J = 1,N
   45 TEMP = 0.
      DO 50 K = 1,NOMEGA
      J2 = J*2-1
   50 TEMP = OMEGA(K)**(J2)*UR(K)*W(K)+TEMP
   55 T(J) = TEMP
      L=N-1
   60 DO 75 J= 1,L
   65 TEMP = 0.
C=======================================================================
      DO 70 K = 1,NOMEGA
      J2 = J*2
   70 TEMP = OMEGA(K)**(J2)*R(K)*W(K)+TEMP
   75 S(J) = TEMP
   80 TEMP = 0.
      DO 81 K = 1,NOMEGA
   81 TEMP = R(K)*W(K)+TEMP
      SZERO = TEMP
C
C     SUBSTITUTION OF PROPER MAGNITUDES IN THE MATRIX
C
      DO 2000 J = 1,N
      DO 2000 I = 1,N
      K = I+J-1
 2000 A(I,J) = AMDA(K)
      MST1 = N+1
      DO 2020 J = MST1,N2
      DO 2020 I = MST1,N2
      K = I + J-1-N2
 2020 A(I,J) = U(K)
      DO 6000 I = 1,N
      IV = 2*I
      V(IV-1) = T(I)
 6000 V(IV) = S(I)
      DO 6005 I = MST1,N2
      DO 6005 J = 1,N
      K = I + J-N-1
 6005 A(I,J) = V(K)
      DO 6010 I = 1,N
      DO 6010 J = MST1,N2
      K = I+J-N-1
 6010 A(I,J) = V(K)
C
C     ASSIGNING CORRECT SIGNS TO THE MATRIX ELEMENTS
C
      LL=N+2
      L=N+1
      J=2
 8000 DO 8010 I=L,N2,2
 8010 A(I,J)=-A(I,J)
      J=J+1
      IF(N-J)8100,8020,8020
 8020 DO 8030 I=1,N2
 8030 A(I,J)=-A(I,J)
      J=J+1
      IF(N-J)8100,8035,8035
 8035 DO 8040 I=2,N,2
 8040 A(I,J)=-A(I,J)
      DO 8050 I=LL,N2,2
 8050 A(I,J)=-A(I,J)
      J=J+2
      IF(N-J)8100,8000,8000
 8100 J=N+1
 8105 DO 8110 I=2,N,2
 8110 A(I,J)=-A(I,J)
      J=J+2
      IF(N2-J)8500,8120,8120
 8120 DO 8130 I=1,N,2
 8130 A(I,J)=-A(I,J)
      DO 8140 I=L,N2,2
 8140 A(I,J)=-A(I,J)
      J=J+1
C=======================================================================
      IF(N2-J)8500,8150,8150
 8150 DO 8160 I=1,N2
 8160 A(I,J)=-A(I,J)
      J=J+1
      IF(N2-J)8500,8105,8105
 8500 B(1) = SZERO
      DO 8501 I=2,N
 8501 B(I) = V(I-1)
      B(N+1)=0.
      NTWO = N+2
      DO 8505 I=NTWO,N2
      K=I-N-1
 8505 B(I) = U(K)
C
C     PRINT MATRIX ELEMENTS ROW BY ROW IF SENSE SWITCH 6 IS DEPRESSED
C
      IF(SENSE SWITCH 5)8510,8520
 8510 PRINT 3755,IDNO,ITER
      IF(SENSE SWITCH 6)8513,8520
 8513 PRINT 8511
 8511 FORMAT(42H    THE MATRIX ELEMENTS PRINTED ROW BY ROW/)
      DO 8512 I=1,N2
 8512 PRINT3792,(A(I,J),J=1,N2)
C
C     SCALE MATRIX ELEMENTS
C
 8520 DO 8700 I=1,N2
      DO 8700 J=1,N2
 8700 A(I,J)=A(I,J)*SCALE
      DO 8701 I=1,N2
 8701 B(I)=B(I)*SCALE
C
C     SOLUTION OF MATRIX EQUATION
C
      CALL MATINV(A,N2,B,1,DETRM)
      IF ACCUMULATOR OVERFLOW 8540,8530
 8530 IF QUOTIENT OVERFLOW 8550,8531
 8540 PRINT 8545,IDNO,ITER
 8545 FORMAT(9H1PROBLEM I6,42H HAD AN ACCUMULATOR OVERFLOW IN ITERATION
     1I3)
      GO TO 1000
 8550 PRINT 8555,IDNO,ITER
 8555 FORMAT(9H1PROBLEM I6,38H HAD A QUOTIENT OVERFLOW IN ITERATION I3)
      GO TO 1000
 8531 PZERO = B(1)
      N1=N-1
      DO 2500 I=1,N1
 2500 P(I) = B(I+1)
      N1=N+1
      DO 2505 I=N1,N2
      IMN=I-N
 2505 Q(IMN) = B(I)
C
C     CALCULATION OF THE MAGNITUDE,PHASE AND ERROR
C
      IF(XMODF(N,2))2510,2510,2509
 2509 NSTOP = (N-1)/2
      GO TO 2516
 2510 NSTOP = N/2
 2516 DO 2521 J = 1,NOMEGA
 2515 TEMP = 1.0
      DO 2520 I = 1,NSTOP
C=======================================================================
      K = 2*I
 2520 TEMP = (-1.0)**I*OMEGA(J)**K*Q(K)+TEMP
 2521 QABS1(J) = TEMP
 2525 IF (XMODF(N,2))2530,2530,2529
 2529 NSTOP = (N+1)/2
      GO TO 2535
 2530 NSTOP = N/2
 2535 DO 2550 J = 1,NOMEGA
 2540 TEMP = 0.
      DO 2545 I = 1,NSTOP
      K = 2*I-1
 2545 TEMP = -(-1.)**I*OMEGA(J)**K*Q(K)+TEMP
 2550 QABS2(J) = TEMP
      DO 2560 J = 1,NOMEGA
      QABSQ(J) = QABS1(J)**2+QABS2(J)**2
      W(J) = 1.0/QABSQ(J)
 2560 QABS(J) = SQRTF(QABSQ(J))
 2600 IF (XMODF(N,2) )2605,2605,2610
 2605 NSTOP = (N-2)/2
      GO TO 2620
 2610 NSTOP = (N-1)/2
 2620 DO 2626 J = 1,NOMEGA
 2621 TEMP = PZERO
      DO 2625 I = 1,NSTOP
      K = 2*I
 2625 TEMP = (-1.)**I*OMEGA(J)**K*P(K)+TEMP
 2626 PABS1(J) = TEMP
 2630 IF (XMODF(N,2))2640,2640,2641
 2640 NSTOP = N/2
      GO TO 2645
 2641 NSTOP = (N-1)/2
 2645 DO 2655 J = 1,NOMEGA
 2646 TEMP = 0.
      DO 2650 I = 1,NSTOP
      K = 2*I-1
 2650 TEMP = -(-1.)**I*OMEGA(J)**K*P(K)+TEMP
 2655 PABS2(J) = TEMP
 2660 DO 2670 J = 1,NOMEGA
      PABSQ(J) = PABS1(J)**2+PABS2(J)**2
      PABS(J) = SQRTF(PABSQ(J))
      CEM(J) = PABS(J)/QABS(J)
      CDBM(J)=8.6858896*LOGF(CEM(J))
      ERORM(J)=ABSF(DBM(J)-CDBM(J))
      ALTER1 = ABSF(CDBM(J))
      IF (ERORM(J)-CURNER*ALTER1)2661,2661,2663
 2663 IF(ITER-1)2661,2661,2665
 2665 DBM(J) = CDBM(J)
      PRINT 2664,J
 2664 FORMAT(///66H    EXPERIMENTAL MAGNITUDE CHANGED TO CALCULATED MAGN
     1ITUDE AT J = I4/)
 2661 ALPA(J) = ABSF(PABS2(J)/PABS1(J))
      IF  QUOTIENT   OVERFLOW 7300,7301
 7301 THETA1(J) = ATANF(ALPA(J))
 7304 BETA(J) = ABSF(QABS2(J)/QABS1(J))
      IF QUOTIENT    OVERFLOW 7302,7303
 7503 THETA2(J) = ATANF(BETA(J))
      GO TO 7305
 7300 THETA1(J)=1.570796425
      GO TO 7304
 7302 THETA2(J)=1.570796425
 7305 IF(QABS2(J))2700,2701,2701
 2700 THETA2(J) = -THETA2(J)
C=======================================================================
 2701 IF (QABS1(J))2710,2711,2711
 2710 THETA2(J) = 3.1415927-THETA2(J)
 2711 IF (PABS2(J))2720,2721,2721
 2720 THETA1(J) = -THETA1(J)
 2721 IF (PABS1(J))2730,2731,2731
 2730 THETA1(J) = 3.1415927-THETA1(J)
 2731 CPHR(J) = THETA1(J)-THETA2(J)
      CPH(J) = CPHR(J)*57.2957795
      ERORP(J)=ABSF(CPH(J)-PHE(J))
      ALTER2 = ABSF(CPH(J))
      IF (ERORP(J)-WALT*ALTER2)2670,2670,2671
 2671 IF(ITER-1)2670,2670,2673
 2673 PHE(J) = CPH(J)
      PRINT 2672,J
 2672 FORMAT(///58H    EXPERIMENTAL PHASE CHANGED TO CALCULATED PHASE AT
     1 J = I4/)
 2670 CONTINUE
      ERMAG2= ERORM(1)
      MAXPT2 = 1
      DO 3700 I = 2,NOMEGA
      IF(ERMAG2-ERORM(I))3710,3700,3700
 3710 ERMAG2= ERORM(I)
      MAXPT2 = I
 3700 CONTINUE
      ERFAZ2 = ERORP(1)
      MAXPT4 = 1
      DO 3720 I = 2,NOMEGA
      IF (ERFAZ2-ERORP(I))3715,3720,3720
 3715 ERFAZ2 = ERORP(I)
      MAXPT4 = I
 3720 CONTINUE
C
C     PRINT RESULTS OF CURRENT ITERATION IF SENSE SWITCH 5 IS DEPRESSED
C
      NLESS1=N-1
      IF (SENSE SWITCH 5)3750,9000
 3755 FORMAT(/////111H11520/RE277   COMPLEX CURVE FITTING ROUTINE
     1                                              PROBLEM NUMBER I4//1
     21H ITERATION I4/)
 3750 PRINT 3760,ERMAG2,MAXPT2,ERFAZ2,MAXPT4
 3760 FORMAT(/18H          ERMAG = F10.5,17H  OCCURED AT J = I4,
     114H     ERPASE = F10.5,17H  OCCURED AT J = I4)
      PRINT 3765
 3765 FORMAT(/50H    NUMERATOR COEFFICIENTS P(0),P(1),P(2),ETC. ARE)
      PRINT 3766,PZERO,(P(I),I=1,NLESS1)
 3766 FORMAT(/1PE20.5,1P6E15.5)
      PRINT 3770
 3770 FORMAT(/52H    DENOMINATOR COEFFICIENTS Q(0),Q(1),Q(2),ETC. ARE)
      PRINT 3766,QZERO,(Q(I),I = 1,N)
      PRINT 3775
 3775 FORMAT(/88H   J       OMEGA           CDBM           CPH
     1ERORM          ERORP            W/)
      PRINT 3780,(J,OMEGA(J),CDBM(J),CPH(J),ERORM(J),ERORP(J),W(J),J = 1
     1,NOMEGA)
 3780 FORMAT(I4,1P6E15.5)
      IF (SENSE SWITCH 6) 3781,9000
 3781 PRINT 3782
 3782 FORMAT(/28H    LAMDAS ARE THE FOLLOWING)
      PRINT 3766,(AMDA(J),J=1,NOLAMD)
      PRINT 3784
 3784 FORMAT(/27H    S(0),S(1),S(2),ETC. ARE)
      PRINT 3766,SZERO,(S(J),J=1,NLESS1)
C=======================================================================
      PRINT 3786
 3786 FORMAT(/27H    T(1),T(2),T(3),ETC. ARE)
      PRINT 3766,(T(J),J=1,N)
      PRINT 3788
 3788 FORMAT(/27H    U(1),U(2),U(3),ETC. ARE)
      PRINT 3766,(U(J),J=1,N2)
 3792 FORMAT(/1P8E15.5/)
C
C     TEST FOR CONVERGENCE
C
 9000 IF(MAXPT1-MAXPT2)9100,9050,9100
 9050 TEST=ABSF(ERMAG1-ERMAG2)
      IF(TEST-CONV1)9051,9051,9100
 9051 IF(MAXPT3-MAXPT4)9100,9052,9100
 9052 TEST=ABSF(ERFAZ1-ERFAZ2)
      IF(TEST-CONV2)9010,9010,9100
 9100 MAXPT1=MAXPT2
      MAXPT3=MAXPT4
      ERMAG1=ERMAG2
      ERFAZ1=ERFAZ2
      ITER=ITER+1
      GO TO 7
C
C     WRITE OUTPUT TAPE 6
C
 9010 WRITE OUTPUT TAPE 6,3755,IDNO,ITER
      WRITE OUTPUT TAPE 6,3760,ERMAG2,MAXPT2,ERFAZ2,MAXPT4
      WRITE OUTPUT TAPE 6,3765
      WRITE OUTPUT TAPE 6,3766,PZERO,(P(I),I=1,NLESS1)
      WRITE OUTPUT TAPE 6,3770
      WRITE OUTPUT TAPE 6,3766,QZERO,(Q(I),I = 1,N)
      WRITE OUTPUT TAPE 6,3775
      WRITE OUTPUT TAPE 6,3780,(J,OMEGA(J),CDBM(J),CPH(J),ERORM(J),ERORP
     1(J),W(J),J=1,NOMEGA)
 9999 GO TO 1000
      END(0,1,0,0,0)


