      SUBROUTINE CUBTRI(F, T, EPS, MCALLS, ANS, ERR, NCALLS, W, NW,
     * IDATA, RDATA, IER)
C
C       ADAPTIVE CUBATURE OVER A TRIANGLE
C
C       PARAMETERS
C          F     - USER SUPPLIED EXTERNAL FUNCTION OF THE FORM
C                  F(X,Y,IDATA,RDATA)
C                  WHERE X AND Y ARE THE CARTESIAN COORDINATES OF A
C                  POINT IN THE PLANE, AND IDATA AND RDATA ARE INTEGER
C                  AND REAL VECTORS IN WHICH DATA MAY BE PASSED.
C          T     - ARRAY OF DIMENSION (2,3) WHERE T(1,J) AND T(2,J)
C                  ARE THE X AND Y COORDINATES OF THE J-TH VERTEX OF
C                  THE GIVEN TRIANGLE (INPUT)
C          EPS   - REQUIRED TOLERANCE (INPUT).  IF THE COMPUTED
C                  INTEGRAL IS BETWEEN-1 AND 1, AN ABSOLUTE ERROR
C                  TEST IS USED, ELSE A RELATIVE ERROR TEST IS USED.
C          MCALLS- MAXIMUM PERMITTED NUMBER OF CALLS TO F (INPUT)
C          ANS   - ESTIMATE FOR THE VALUE OF THE INTEGRAL OF F OVER
C                  THE GIVEN TRIANGLE (OUTPUT)
C          ERR   - ESTIMATED ABSOLUTE ERROR IN ANS (OUTPUT)
C          NCALLS- ACTUAL NUMBER OF CALLS TO F (OUTPUT).  THIS
C                  PARAMETER MUST BE INITIALIZED TO 0 ON THE FIRST
C                  CALL TO CUBTRI FOR A GIVEN INTEGRAL (INPUT)
C          W     - WORK SPACE.  MAY NOT BE DESTROYED BETWEEN CALLS TO
C                  CUBTRI IF RESTARTING IS INTENDED
C          NW    - LENGTH OF WORK SPACE (INPUT).
C                  IF NW .GE. 3*(19+3*MCALLS)/38, TERMINATION DUE TO
C                  FULL WORK SPACE WILL NOT OCCUR.
C          IER   - TERMINATION INDICATOR (OUTPUT)
C                  IER=0   NORMAL TERMINATION, TOLERANCE SATISFIED
C                  IER=1   MAXIMUM NUMBER OF CALLS REACHED
C                  IER=2   WORK SPACE FULL
C                  IER=3   FURTHER SUBDIVISION OF TRIANGLES IMPOSSIBLE
C                  IER=4   NO FURTHER IMPROVEMENT IN ACCURACY IS
C                        POSSIBLE DUE TO ROUNDING ERRORS IN FUNCTION
C                        VALUES
C                  IER=5   NO FURTHER IMPROVEMENT IN ACCURACY IS
C                        POSSIBLE BECAUSE SUBDIVISION DOES NOT
C                        CHANGE THE ESTIMATED INTEGRAL. MACHINE
C                        ACCURACY HAS PROBABLY BEEN REACHED BUT
C                        THE ERROR ESTIMATE IS NOT SHARP ENOUGH.
C
C       CUBTRI IS DESIGNED TO BE CALLED REPEATEDLY WITHOUT WASTING
C       EARLIER WORK.  THE PARAMETER NCALLS IS USED TO INDICATE TO
C       CUBTRI AT WHAT POINT TO RESTART, AND MUST BE RE-INITIALIZED
C       TO 0 WHEN A NEW INTEGRAL IS TO BE COMPUTED.  AT LEAST ONE OF
C       THE PARAMETERS EPS, MCALLS AND NW MUST BE CHANGED BETWEEN
C       CALLS TO CUBTRI, ACCORDING TO THE RETURNED VALUE OF IER. NONE
C       OF THE OTHER PARAMETERS MAY BE CHANGED IF RESTARTING IS DONE.
C       IF IER=3 IS ENCOUNTERED, THERE PROBABLY IS A SINGULARITY
C       SOMEWHERE IN THE REGION.  THE ERROR MESSAGE INDICATES THAT
C       FURTHER SUBDIVISION IS IMPOSSIBLE BECAUSE THE VERTICES OF THE
C       SMALLER TRIANGLES PRODUCED WILL BEGIN TO COALESCE TO THE
C       PRECISION OF THE COMPUTER.  THIS SITUATION CAN USUALLY BE
C       RELIEVED BY SPECIFYING THE REGION IN SUCH A WAY THAT THE
C       SINGULARITY IS LOCATED AT THE THIRD VERTEX OF THE TRIANGLE.
C       IF IER=4 IS ENCOUNTERED, THE VALUE OF THE INTEGRAL CANNOT BE
C       IMPROVED ANY FURTHER. THE ONLY EXCEPTION TO THIS OCCURS WHEN A
C       FUNCTION WITH HIGHLY IRREGULAR BEHAVIOUR IS INTEGRATED (E.G.
C       FUNCTIONS WITH JUMP DISCONTINUITIES OR VERY HIGHLY OSCILLATORY
C       FUNCTIONS). IN SUCH A CASE THE USER CAN DISABLE THE ROUNDING
C       ERROR TEST BY REMOVING THE IF STATEMENT SHORTLY AFTER STATEMENT
C       NUMBER 70.
C
      EXTERNAL F
      INTEGER IDATA(1), IER, MCALLS, NCALLS, NW
      REAL ALFA, ANS, ANSKP, AREA, EPS, ERR, ERRMAX, H, Q1, Q2, R1, R2,
     * RDATA(1), D(2,4), S(4), T(2,3), VEC(2,3), W(6,NW), X(2)
C       ACTUAL DIMENSION OF W IS (6,NW/6)
C
      DOUBLE PRECISION TANS, TERR, DZERO
      COMMON /CUBSTA/ TANS, TERR
C       THIS COMMON IS REQUIRED TO PRESERVE TANS AND TERR BETWEEN CALLS
C       AND TO SAVE VARIABLES IN FUNCTION RNDERR
      DATA NFE /19/, S(1), S(2), S(3), S(4) /3*1E0,-1E0/, D(1,1),
     * D(2,1) /0.0,0.0/, D(1,2), D(2,2) /0.0,1.0/, D(1,3), D(2,3)
     * /1.0,0.0/, D(1,4), D(2,4) /1.0,1.0/
C       NFE IS THE NUMBER OF FUNCTION EVALUATIONS PER CALL TO CUBRUL.
      DATA ZERO /0.E0/, ONE /1.E0/, DZERO /0.D0/, POINT5 /.5E0/
C
C      CALCULATE DIRECTION VECTORS, AREA AND MAXIMUM NUMBER
C      OF SUBDIVISIONS THAT MAY BE PERFORMED
      DO 20 I=1,2
        VEC(I,3) = T(I,3)
        DO 10 J=1,2
          VEC(I,J) = T(I,J) - T(I,3)
   10   CONTINUE
   20 CONTINUE
      MAXC = (MCALLS/NFE+3)/4
      IER = 1
      MAXK = MIN0(MAXC,(NW/6+2)/3)
      IF (MAXC.GT.MAXK) IER = 2
      AREA = ABS(VEC(1,1)*VEC(2,2)-VEC(1,2)*VEC(2,1))*POINT5
      K = (NCALLS/NFE+3)/4
      MW = 3*(K-1) + 1
      IF (NCALLS.GT.0) GO TO 30
C
C       TEST FOR TRIVIAL CASES
      TANS = DZERO
      TERR = DZERO
      IF (AREA.EQ.ZERO) GO TO 90
      IF (MCALLS.LT.NFE) GO TO 100
      IF (NW.LT.6) GO TO 110
C
C       INITIALIZE DATA LIST
      K = 1
      MW = 1
      W(1,1) = ZERO
      W(2,1) = ZERO
      W(3,1) = ONE
      CALL CUBRUL(F, VEC, W(1,1), IDATA, RDATA)
      TANS = W(5,1)
      TERR = W(6,1)
      NCALLS = NFE
C
C       TEST TERMINATION CRITERIA
   30 ANS = TANS
      ERR = TERR
      IF (ERR.LT.AMAX1(ONE,ABS(ANS))*EPS) GO TO 90
      IF (K.EQ.MAXK) GO TO 120
C
C       FIND TRIANGLE WITH LARGEST ERROR
      ERRMAX = ZERO
      DO 40 I=1,MW
        IF (W(6,I).LE.ERRMAX) GO TO 40
        ERRMAX = W(6,I)
        J = I
   40 CONTINUE
C
C       SUBDIVIDE TRIANGLE INTO FOUR SUBTRIANGLES AND UPDATE DATA LIST
      DO 50 I=1,2
        X(I) = W(I,J)
   50 CONTINUE
      H = W(3,J)*POINT5
      IF (RNDERR(X(1),H,X(1),H).NE.ZERO) GO TO 130
      IF (RNDERR(X(2),H,X(2),H).NE.ZERO) GO TO 130
      ANSKP = SNGL(TANS)
      TANS = TANS - DBLE(W(5,J))
      TERR = TERR - DBLE(W(6,J))
      R1 = W(4,J)
      R2 = W(5,J)
      JKP = J
      Q1 = ZERO
      Q2 = ZERO
      DO 70 I=1,4
        DO 60 L=1,2
          W(L,J) = X(L) + H*D(L,I)
   60   CONTINUE
        W(3,J) = H*S(I)
        CALL CUBRUL(F, VEC, W(1,J), IDATA, RDATA)
        Q2 = Q2 + W(5,J)
        Q1 = Q1 + W(4,J)
        J = MW + I
   70 CONTINUE
      ALFA = 1E15
      IF (Q2.NE.R2) ALFA = ABS((Q1-R1)/(Q2-R2)-ONE)
      J = JKP
      DO 80 I=1,4
        W(6,J) = W(6,J)/ALFA
        TANS = TANS + W(5,J)
        TERR = TERR + W(6,J)
        J = MW + I
   80 CONTINUE
      MW = MW + 3
      NCALLS = NCALLS + 4*NFE
      K = K + 1
C
C       IF ANSWER IS UNCHANGED, IT CANNOT BE IMPROVED
      IF (ANSKP.EQ.SNGL(TANS)) GO TO 150
C
C       REMOVE THIS IF STATEMENT TO DISABLE ROUNDING ERROR TEST
      IF (K.GT.3 .AND. ABS(Q2-R2).GT.ABS(Q1-R1)) GO TO 140
      GO TO 30
C
C       EXITS FROM SUBROUTINE
   90 IER = 0
      GO TO 120
  100 IER = 1
      GO TO 120
  110 IER = 2
  120 ANS = TANS
      ERR = TERR
      RETURN
  130 IER = 3
      GO TO 120
  140 IER = 4
      GO TO 120
  150 IER = 5
      GO TO 120
      END
