C.IRATCU........D.......................................................
C
C   FUNCTION            - RATIONAL CHEBYCHEV APPROXIMATION TO
C                           THE FUNCTION F
C   USAGE               - CALL IRATCU (F,PHI,G,AA,BB,L,M,P,Q,WK,IER)
C   PARAMETERS   F      - USER SUPPLIED EXTERNAL FUNCTION
C                           SUBPROGRAM F(X). F IS TO BE APPROXIMATED
C                           OVER THE INTERVAL (AA,BB)
C                PHI    - VARIABLE TRANSFORMATION WHICH MUST BE
C                           CONTINUOUS AND MONOTONIC. IT IS SPECIFIED
C                           BY A USER SUPPLIED EXTERNAL FUNCTION
C                           SUBPROGRAM
C                G      - WEIGHT FUNCTION WHICH MUST BE CONTINUOUS
C                           AND NON VANISHING IN (AA,BB).
C                           IT IS SPECIFIED BY A USER SUPPLIED EXTERNAL
C                           FUNCTION SUBPROGRAM
C                AA     - LOWER END POINT OF INTERVAL (AA,BB)
C                BB     - UPPER END POINT OF INTERVAL (AA,BB)
C                L      - DEGREE OF NUMERATOR OF RATIONAL APPROXIMATION
C                M      - DEGREE OF DENOMINATOR OF RATIONAL
C                           APPROXIMATION. M MUST NOT BE GREATER THAN L
C                P      - OUTPUT ARRAY P(1),P(2),...,P(L+1) OF
C                           NUMERATOR COEFFICIENTS
C                Q      - OUTPUT ARRAY Q(1),Q(2),...,Q(M+1) OF
C                           DENOMINATOR COEFFICIENTS
C                WK     - WORK AREA OF DIMENSION .GE. (NP1+6)*NP1
C                           WHERE NP1=L+M+2
C                IER    - ERROR PARAMETER
C                         TERMINAL ERROR = 128+N
C                           N = 1 INDICATES A POLE IN THE TRIAL
C                             RATIONAL APPROXIMANT DUE TO A SEQ. OF
C                             CRITICAL POINTS WHICH IS NOT MONOTONIC
C                           N = 2 INDICATES A POLE IN THE TRIAL
C                             RATIONAL APPROXIMANT DUE TO A LARGE
C                             DEL(X(I)) WHERE X(I) IS A TRIAL CRITICAL
C                             POINT
C                         WARNING ERROR = 32+N
C                           N = 1 INDICATES ALGORITHM DID NOT CONVERGE
C                             WITHIN 10 ITERATIONS
C   PRECISION           - DOUBLE
C   REQ'D IMSL ROUTINES - UERTST
C   AUTHOR/IMPLEMENTOR  - O.G. JOHNSON/C.L. SMITH
C   LANGUAGE            - FORTRAN
C.......................................................................
C   LATEST REVISION     - MARCH 29, 1971
C
      SUBROUTINE IRATCU  (F,PHI,G,AA,BB,L,M,P,Q,WK,IER)
C
      DIMENSION          P(1),Q(1),WK(1)
      REAL*8             P,Q,A,WK,AA,BB,XLB,DN,DNP1,XXI
      REAL*8             BB1,XXM1,ZZ,TEST,TEMPL,U,Z1,Y2,H1,Y3,PZ,T1,T2
      REAL*8             T3,T4,DEL,Z2,Z,Y,TEMP1,TEMPL1,SUM,PHI,F,G
      INTEGER            V
      DATA               PI,EPS,EQUIB/3.1415926535897932,.0001,1.E-12/
      IER = 0
      XLB = 0.
      ITNO = 1
      LP1 = L+1
      LPM = L+M
      N = LPM+1
      DN = N
      DNP1 = DN+1.0
      RTEST=0.
      NP1 = N+1
      MP1 = M+1
C                                  SET UP INDEXING FOR WORK AREA
      IB1 = NP1
      IB2 = IB1+NP1
      IB3 = IB2 + NP1
      IB4 = IB3 + NP1
      IB5 = IB4 + NP1
      IB6 = IB5 + NP1
      DO 5  I=1,LP1
         P(I) = 0.
    5 CONTINUE
      DO 10  I=1,MP1
         Q(I) = 0.0
   10 CONTINUE
      WK(1) = 1.0
      DO 15  I=1,NP1
         WK(I+1) = -WK(I)
   15 CONTINUE
      Q(1) = 1.0
C                                  COMPUTE INITIAL GUESSES
C                                  OF CRITICAL POINTS
      XXI = 0.0
      KK = N/2
      BB1 = (BB-AA)*.5
      XXM1 = (BB+AA)*.5
      DO 20  I=1,KK
         XXI = XXI + 1.
         ZZ = DCOS(PI*(1.-XXI/N))*BB1
         IJK4 = IB4+I+1
         WK(IJK4) = ZZ+XXM1
         J = N-I+1
         IJK4 = IB4+J
         WK(IJK4) = XXM1-ZZ
   20 CONTINUE
      WK(IB4+1) = AA
      WK(IB5) = BB
   25 DO 30  I=1,NP1
         IJK4 = IB4+I
         IJK5 = IB5+I
         WK(IJK5) = PHI(WK(IJK4))
   30 CONTINUE
C                                  SET UP AND SOLVE LINEAR EQUATIONS
      KRET = 1
      GO TO 45
   35 XLB = (WK(IB2)+XLB*DN)/DNP1
      IF (M.EQ.0) GO TO 135
      DO 120 I=1,3
         KRET = 2
         GO TO 45
   40    TEST = DABS(XLB-WK(IB2))
         XLB = (WK(IB2)+XLB*DN)/DNP1
         IF (TEST .LE. .01*DABS(XLB)) GO TO 125
         GO TO 120
   45    KK = L+2
C                                  SET UP LINEAR EQUATIONS
         DO 70  V=1,NP1
            IJK = IB6+V
            WK(IJK) = 1.0
            IJK5 = IB5+V
            IF (L.LE.0) GO TO 55
            DO 50  J=2,LP1
               IC1 = IB6+(J-1)*NP1+V
               IC2 = IC1-NP1
               WK(IC1)=WK(IC2)*WK(IJK5)
   50       CONTINUE
   55       IJK1 = IB1 + V
            IJK4 = IB4+V
            WK(IJK1) = F(WK(IJK4))
            IF (M.LE.0) GO TO 65
            IJK0 = V+1
            TEMPL = XLB*G(WK(IJK4))*WK(IJK0)-WK(IJK1)
            IC1 = IB6+(KK-1)*NP1+V
            WK(IC1) = WK(IJK5)*TEMPL
            DO 60 J=KK,LPM
               IC1 = IB6+J*NP1+V
               IC2 = IC1-NP1
               WK(IC1) = WK(IC2)*WK(IJK5)
   60       CONTINUE
            IC1 = IB6+NP1*(NP1-1)+V
            WK(IC1) = WK(IJK0)*G(WK(IJK4))
   65       CONTINUE
   70    CONTINUE
C                                  SOLVE LINEAR EQUATIONS USING GAUSIAN
C                                  ELIMINATION WITH PARTIAL PIVOTING
         DO 105 NN=1,NP1
            IC3 = IB6+(NN-1)*NP1+NN
            TEMP1 = DABS(WK(IC3))
            II = NN
            IJK1 = IB1+NN
            IF (NN.EQ.NP1) GO TO 85
            NNP1 = NN+1
            DO 75  KK=NNP1,NP1
               IC1 = IB6+(NN-1)*NP1+KK
               TEMPL1 = DABS(WK(IC1))
               IF (TEMP1.GE.TEMPL1) GO TO 75
               II = KK
               TEMP1 = TEMPL1
   75       CONTINUE
            IF (II.EQ.NN) GO TO 85
            DO 80  MM=NN,NP1
               IC1 = IB6+(MM-1)*NP1+NN
               IC2 = IC1-NN+II
               TEMP1 = WK(IC1)
               WK(IC1) = WK(IC2)
               WK(IC2) = TEMP1
   80       CONTINUE
            TEMP1 = WK(IJK1)
            IJK2 = IB1+II
            WK(IJK1) = WK(IJK2)
            WK(IJK2) = TEMP1
   85       WK(IJK1) = WK(IJK1)/WK(IC3)
            DO 90  JJ=NN,NP1
               JK = NP1-JJ+NN
               IC1 = IB6+(JK-1)*NP1+NN
               WK(IC1) = WK(IC1)/WK(IC3)
   90       CONTINUE
            IF (NN.EQ.NP1) GO TO 105
            DO 100 II=NNP1,NP1
               IC1 = IB6+(NN-1)*NP1+II
               DO 95  JJ=NNP1,NP1
                  IC2 = IB6+(JJ-1)*NP1+II
                  IC4 = IC2-II+NN
                  WK(IC2) = WK(IC2)-WK(IC1)*WK(IC4)
   95          CONTINUE
               IJKN = IB1+II
               WK(IJKN) = WK(IJKN)-WK(IC1)*WK(IJK1)
  100       CONTINUE
  105    CONTINUE
         DO 115 II=1,NP1
            JJ = NP1-II+1
            IJK1 = IB1+JJ
            IC1 = IB6+(JJ-1)*NP1+JJ
            WK(IJK1) = WK(IJK1)/WK(IC1)
            JJM1 = JJ-1
            IF (JJ.EQ.1) GO TO 115
            DO 110 KK=1,JJM1
               KJ = JJM1-KK+1
               IJKN = IB1+KJ
               IC2 = IC1-JJ+KJ
               WK(IJKN) = WK(IJKN)-WK(IC2)*WK(IJK1)
  110       CONTINUE
  115    CONTINUE
         GO TO (35,40), KRET
  120 CONTINUE
  125 DO 130 I=1,M
         IJK1 = IB1+LP1+I
         Q(I+1) = WK(IJK1)
  130 CONTINUE
  135 CONTINUE
      DO 140  I=1,LP1
         IJK1 = IB1+I
         P(I) = WK(IJK1)
  140 CONTINUE
  145 U = DSIGN(1.D0,XLB)
      Z1 = 0.0
C                                  SEARCH FOR NEW CRITICAL POINTS
      IF (N.GT.1) GO TO 150
      WK(IB3+1) = .015D0*(WK(IB4+2)-WK(IB4+1))
      WK(IB3+2) = -WK(IB3+1)
      GO TO 160
  150 DO 155  I=2,N
         IJK3 = IB3+I
         IJK4 = IB4+I+1
         IJKN = IB4+I-1
         WK(IJK3) = .015D0*(WK(IJK4)-WK(IJKN))
  155 CONTINUE
      WK(IB3+1) = WK(IB3+2)*.5
      IJK3 = IB3+N
      WK(IB4) = -WK(IJK3)*.5
  160 DO 240  I=1,NP1
         IJK4 = IB4+I
         Y2 = WK(IJK4)
         IJK3 = IB3+I
         H1 = WK(IJK3)
         Y3 = Y2+H1
         PZ = Y2
         KRET = 1
C                                  COMPUTE APPROXIMATION ERROR AT PZ
  165    T1 = F(PZ)
         T2 = PHI(PZ)
C                                  ESTIMATE MAXIMUM VALUE OF F
         T5=DABS(T1)
         IF(T5.GT.RTEST) RTEST=T5
         IJK1 = IB1+LP1
         T3 = WK(IJK1)
         DO 170  II=1,L
            KD = LP1-II
            IJK1 = IB1+KD
            T3 = T3*T2+WK(IJK1)
  170    CONTINUE
         T4 = 0.
         DO 175  II=1,M
            KD = NP1-II
            IJK1 = IB1+KD
            T4 = (T4+WK(IJK1))*T2
  175    CONTINUE
         T4 = T4+1.0
         T2 = T3/T4
         DEL = (T2-T1)/G(PZ)
C                                  DEL IS THE APPROXIMATION ERROR AT PZ
C                                  RETURN TO THE CALLING SEGEMENT
         GO TO (180,185,205,225,250,275),KRET
C                                  DIRECT SEARCH FOR ADJUSTED CRITICAL
C                                  POINTS
  180    Z2 = DEL*U
         PZ = Y3
         KRET = 2
         GO TO 165
  185    Z3 = DEL*U
         IF (Z3.GT.Z2) GO TO 190
         H1 = -H1
         Z = Z3
         Z3 = Z2
         Z2 = Z
         Y = Y3
         Y3 = Y2
         Y2 = Y
  190    Y = Y3+H1
         IF (Y.GE.AA) GO TO 195
         Y = AA
         GO TO 220
  195    IF (Y.LE.BB) GO TO 200
         Y = BB
         GO TO 220
  200    PZ = Y
         KRET = 3
         GO TO 165
  205    Z = DEL*U
         IF (Z.LE.Z3) GO TO 210
         Y2 = Y3
         Y3 = Y
         Z2 = Z3
         Z3 = Z
         GO TO 190
  210    Y = -Z3-Z3+Z+Z2
         IF(Y.EQ.0) GO TO 215
         Y = (Y2+Y3)*.5+H1*(Z2-Z3)/Y
         GO TO 220
  215    Y = Y3
  220    WK(IJK4) = Y
         PZ = Y
         KRET = 4
         GO TO 165
  225    IJK2 = IB2+I
         WK(IJK2) = DEL
         U = -U
C                                  CHECK FOR PROPER ORDERING OF
C                                  ADJUSTED CRITICAL POINTS
         IF (.NOT.(I.GT.1.AND.WK(IJK4).LE.WK(IJK4-1))) GO TO 230
         IER = 129
         GO TO 9000
  230    IJK2 = IB2+I
         Z = DABS(WK(IJK2))
C                                  CHECK FOR POLES IN TRIAL
C                                  RATIONAL APPROXIMANT
         IF (Z.LT.10.) GO TO 235
         IER = 130
         GO TO 9000
  235    Y = DABS(XLB)
         Z = DABS(Z-Y)/Y
         IF (Z1.LT.Z) Z1=Z
  240 CONTINUE
C                                  ADJUST END CRITICAL POINTS SO THAT
C                                  THEY DO NOT LIE OUTSIDE THE
C                                  INTERVAL (AA,BB)
      IJK4 = IB4+1
      IF (WK(IJK4).LE.AA) GO TO 265
      H1 = (WK(IJK4)-AA)*.0625
      U = -DSIGN(1.D0,XLB)
      Z3 = 0.0
      Y = AA
      I = 1
  245 PZ = Y
      KRET = 5
      GO TO 165
  250 Z = DEL*U
      IF (Z.LE.Z3) GO TO 255
      Z3 = Z
      Z2 = Y
  255 Y = Y+H1
      I = I+1
      IF (I.LE.16) GO TO 245
      Z = DABS(XLB)
      IF (Z3.LE.Z) GO TO 265
      DO 260  I=2,NP1
         J = NP1-I+2
         IJK2 = IB2+J
         WK(IJK2) = WK(IJK2-1)
         IJK4 = IB4+J
         WK(IJK4) = WK(IJK4-1)
  260 CONTINUE
      IJK4 = IB4+1
      WK(IJK4) = Z2
      IJK2 = IB2+1
      WK(IJK2) = Z3*U
      GO TO 290
  265 IF (BB.LE.WK(IB5)) GO TO 295
      H1 = .0625*(BB-WK(IB5))
      Z3 = 0.0
      Y = BB
      U = -DSIGN(1.D0,WK(IB3))
      I = 1
  270 PZ = Y
      KRET = 6
      GO TO 165
  275 Z = DEL*U
      IF (Z3.GE.Z) GO TO 280
      Z3 = Z
      Z2 = Y
  280 Y = Y-H1
      I = I+1
      IF (I.LE.16) GO TO 270
      Z = DABS(XLB)
      IF (Z.GE.Z3) GO TO 295
      DO 285  I=1,N
         IJK2 = IB2+I
         WK(IJK2) = WK(IJK2+1)
         IJK4 = IB4+I
         WK(IJK4) = WK(IJK4+1)
  285 CONTINUE
      WK(IB5) = Z2
      WK(IB3) = Z3*U
  290 XLB = -XLB
      Z = DABS(Z3-Z)/Z
      IF (Z1.LT.Z) Z1=Z
  295 CONTINUE
C                                  CHECK FOR NO CONVERGENCE
      IF (ITNO.GE.10) GO TO 315
C                                  CHECK FOR CONVERGENCE
      IF (Z1.LE.EPS) GO TO 9005
C                                  TEST FOR SIGNIFICIENT DIGITS IN Z1
      IF(DABS(XLB).LT. RTEST*EQUIB) GO TO 9005
      SUM=0.
      DO 310  I=1,NP1
         IJK2 = IB2+I
         SUM = SUM+WK(I)*WK(IJK2)
  310 CONTINUE
      XLB = SUM/DNP1
      ITNO = ITNO+1
      GO TO 25
  315 IER = 33
 9000 CONTINUE
      CALL UERTST(IER,'IRATCU')
 9005 RETURN
      END
