      SUBROUTINE CUBRUL(F, VEC, P, IDATA, RDATA)
C
C       BASIC CUBATURE RULE PAIR OVER A TRIANGLE
C
C       PARAMETERS
C         F  - EXTERNAL FUNCTION - SEE COMMENTS TO CUBTRI
C         VEC- MATRIX OF BASE VECTORS AND ORIGIN (INPUT)
C         P  - TRIANGLE DESCRIPTION VECTOR OF DIMENSION 6
C               P(1) - TRANSFORMED X COORDINATE OF ORIGIN VERTEX(INPUT)
C               P(2) - TRANSFORMED Y COORDINATE OF ORIGIN VERTEX(INPUT)
C               P(3) - DISTANCE OF OTHER VERTICES IN THE DIRECTIONS
C                     OF THE BASE VECTORS (INPUT)
C               P(4) - LESS ACCURATE ESTIMATED INTEGRAL (OUTPUT)
C               P(5) - MORE ACCURATE ESTIMATED INTEGRAL (OUTPUT)
C               P(6) - ABS(P(5)-P(4))   (OUTPUT)
C
C       CUBRUL EVALUATES A LINEAR COMBINATION OF BASIC INTEGRATION
C       RULES HAVING D3 SYMMETRY.  THE AREAL COORDINATES PERTAINING TO
C       THE J-TH RULE ARE STORED IN W(I,J),I=1,2,3.  THE CORRESPONDING
C       WEIGHTS ARE W(4,J) AND W(5,J), WITH W(5,J) BELONGING TO THE
C       MORE ACCURATE FORMULA.  IF W(1,J).EQ.W(2,J), THE INTEGRATION
C       POINT IS THE CENTROID, ELSE IF W(2,J).EQ.W(3,J), THE EVALUATION
C       POINTS ARE ON THE MEDIANS.  IN BOTH CASES ADVANTAGE IS TAKEN OF
C       SYMMETRY TO AVOID REPEATING FUNCTION EVALUATIONS.
C
C       THE FOLLOWING DOUBLE PRECISION VARIABLES ARE USED TO AVOID
C       UNNECESSARY ROUNDING ERRORS IN FLOATING POINT ADDITION.
C       THEY MAY BE DECLARED SINGLE PRECISION IF DOUBLE PRECISION IS
C       NOT AVAILABLE AND FULL ACCURACY IS NOT NEEDED.
      DOUBLE PRECISION A1, A2, S, SN, DZERO, DONE, DTHREE, DSIX
      REAL AREA, ORIGIN(2), P(6), RDATA(1), TVEC(2,3), VEC(2,3), W(5,6)
      INTEGER IDATA(1)
      EXTERNAL F
C
C       W CONTAINS POINTS AND WEIGHTS OF THE INTEGRATION FORMULAE
C       NQUAD - NUMBER OF BASIC RULES USED
C
C       THIS PARTICULAR RULE IS THE 19 POINT EXTENSION (DEGREE 8) OF
C       THE FAMILIAR 7 POINT RULE (DEGREE 5).
C
C     SIGMA=SQRT(7)
C     PHI=SQRT(15)
C     W(1,1),W(2,1),W(3,1) = 1/3
C     W(4,1) = 9/40
C     W(5,1) = 7137/62720 - 45*SIGMA/1568
C     W(1,2) = 3/7 + 2*PHI/21
C     W(2,2),W(3,2) = 2/7 - PHI/21
C     W(4,2) = 31/80 - PHI/400
C     W(5,2) = - 9301697/4695040 - 13517313*PHI/23475200
C            + 764885*SIGMA/939008 + 198763*PHI*SIGMA/939008
C     W(*,3) = W(*,2) WITH PHI REPLACED BY -PHI
C     W(1,5) = 4/9 + PHI/9 + SIGMA/9 - SIGMA*PHI/45
C     W(2,5),W(3,5) = 5/18 - PHI/18 - SIGMA/18 + SIGMA*PHI/90
C     W(4,5) = 0
C     W(5,5) = 102791225/59157504 + 23876225*PHI/59157504
C            - 34500875*SIGMA/59157504 - 9914825*PHI*SIGMA/59157504
C     W(*,4) = W(*,5) WITH PHI REPLACED BY -PHI
C     W(1,6) = 4/9 + SIGMA/9
C     W(2,6) = W(2,4)
C     W(3,6) = W(2,5)
C     W(4,6) = 0
C     W(5,6) = 11075/8064 - 125*SIGMA/288
C
      DATA NQUAD /6/, W(1,1), W(2,1), W(3,1) /3*.33333333333333333333333
     * 33E0/, W(4,1), W(5,1) /.225E0,.3786109120031468330830822E-1/,
     * W(1,2), W(2,2), W(3,2) /.7974269853530873223980253E0,2*
     * .1012865073234563388009874E0/, W(4,2), W(5,2)
     * /.3778175416344814577870518E0,.1128612762395489164329420E0/,
     * W(1,3), W(2,3), W(3,3) /.5971587178976982045911758E-1,2*
     * .4701420641051150897704412E0/, W(4,3), W(5,3)
     * /.3971824583655185422129482E0,.2350720567323520126663380E0/
      DATA W(1,4), W(2,4), W(3,4) /.5357953464498992646629509E0,2*
     * .2321023267750503676685246E0/, W(4,4), W(5,4)
     * /0.E0,.3488144389708976891842461E0/, W(1,5), W(2,5), W(3,5)
     * /.9410382782311208665596304E0,2*.2948086088443956672018481E-1/,
     * W(4,5), W(5,5) /0.E0,.4033280212549620569433320E-1/, W(1,6),
     * W(2,6), W(3,6) /.7384168123405100656112906E0,
     * .2321023267750503676685246E0,.2948086088443956672018481E-1/,
     * W(4,6), W(5,6) /0.E0,.2250583347313904927138324E0/
C
      DATA DZERO /0.D0/, DONE /1.D0/, DTHREE /3.D0/, DSIX /6.D0/,
     * POINT5 /.5E0/
C
C       SCALE BASE VECTORS AND OBTAIN AREA
      DO 20 I=1,2
        ORIGIN(I) = VEC(I,3) + P(1)*VEC(I,1) + P(2)*VEC(I,2)
        DO 10 J=1,2
          TVEC(I,J) = P(3)*VEC(I,J)
   10   CONTINUE
   20 CONTINUE
      AREA = POINT5*ABS(TVEC(1,1)*TVEC(2,2)-TVEC(1,2)*TVEC(2,1))
      A1 = DZERO
      A2 = DZERO
C
C       COMPUTE ESTIMATES FOR INTEGRAL AND ERROR
      DO 40 K=1,NQUAD
        X = ORIGIN(1) + W(1,K)*TVEC(1,1) + W(2,K)*TVEC(1,2)
        Y = ORIGIN(2) + W(1,K)*TVEC(2,1) + W(2,K)*TVEC(2,2)
        S = DBLE(F(X,Y,IDATA,RDATA))
        SN = DONE
        IF (W(1,K).EQ.W(2,K)) GO TO 30
        X = ORIGIN(1) + W(2,K)*TVEC(1,1) + W(1,K)*TVEC(1,2)
        Y = ORIGIN(2) + W(2,K)*TVEC(2,1) + W(1,K)*TVEC(2,2)
        S = S + DBLE(F(X,Y,IDATA,RDATA))
        X = ORIGIN(1) + W(2,K)*TVEC(1,1) + W(3,K)*TVEC(1,2)
        Y = ORIGIN(2) + W(2,K)*TVEC(2,1) + W(3,K)*TVEC(2,2)
        S = S + DBLE(F(X,Y,IDATA,RDATA))
        SN = DTHREE
        IF (W(2,K).EQ.W(3,K)) GO TO 30
        X = ORIGIN(1) + W(1,K)*TVEC(1,1) + W(3,K)*TVEC(1,2)
        Y = ORIGIN(2) + W(1,K)*TVEC(2,1) + W(3,K)*TVEC(2,2)
        S = S + DBLE(F(X,Y,IDATA,RDATA))
        X = ORIGIN(1) + W(3,K)*TVEC(1,1) + W(1,K)*TVEC(1,2)
        Y = ORIGIN(2) + W(3,K)*TVEC(2,1) + W(1,K)*TVEC(2,2)
        S = S + DBLE(F(X,Y,IDATA,RDATA))
        X = ORIGIN(1) + W(3,K)*TVEC(1,1) + W(2,K)*TVEC(1,2)
        Y = ORIGIN(2) + W(3,K)*TVEC(2,1) + W(2,K)*TVEC(2,2)
        S = S + DBLE(F(X,Y,IDATA,RDATA))
        SN = DSIX
   30   S = S/SN
        A1 = A1 + W(4,K)*S
        A2 = A2 + W(5,K)*S
   40 CONTINUE
      P(4) = SNGL(A1)*AREA
      P(5) = SNGL(A2)*AREA
      P(6) = ABS(P(5)-P(4))
      RETURN
      END
