1      SUBROUTINE LQM0(F,U,V,RES8,EST)
2C
3C
4C***REVISION HISTORY (YYMMDD)
5C   000601   Changed AMAX1/AMIN1 to generic MAX/MIN
6C
7C      PURPOSE
8C           TO COMPUTE - IF = INTEGRAL OF F OVER THE TRIANGLE
9C           WITH VERTICES (U(1),V(1)),(U(2),V(2)),(U(3),V(3)), AND
10C           ESTIMATE THE ERROR,
11C                      - INTEGRAL OF ABS(F) OVER THIS TRIANGLE
12C
13C      CALLING SEQUENCE
14C           CALL LQM0(F,U,V,RES11,EST)
15C        PARAMETERS
16C           F       - FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
17C                     F(X,Y); THE ACTUAL NAME FOR F NEEDS TO BE
18C                     DECLARED E X T E R N A L IN THE CALLING
19C                     PROGRAM
20C           U(1),U(2),U(3)- ABSCISSAE OF VERTICES
21C           V(1),V(2),V(3)- ORDINATES OF VERTICES
22C           RES6    - APPROXIMATION TO THE INTEGRAL IF, OBTAINED BY THE
23C                     LYNESS AND JESPERSEN RULE OF DEGREE 6, USING
24C                     12 POINTS
25C           RES8   - APPROXIMATION TO THE INTEGRAL IF, OBTAINED BY THE
26C                     LYNESS AND JESPERSEN RULE OF DEGREE 8,
27C                     USING 16 POINTS
28C           EST     - ESTIMATE OF THE ABSOLUTE ERROR
29C           DRESC   - APPROXIMATION TO THE INTEGRAL OF ABS(F- IF/DJ),
30C                     OBTAINED BY THE RULE OF DEGREE 6, AND USED FOR
31C                     THE COMPUTATION OF EST
32C
33C      REMARKS
34C           DATE OF LAST UPDATE : 10 APRIL 1984 O.W. RECHARD NBS
35C
36C           SUBROUTINES OR FUNCTIONS CALLED :
37C                   - F (USER-SUPPLIED INTEGRAND FUNCTION)
38C                   - R1MACH FOR MACHINE DEPENDENT INFORMATION
39C
40C .....................................................................
41C
42      REAL DJ,DF0,DRESC,EMACH,EST,F,FV,F0,
43     *  RES6,RES8,U1,U2,U3,UFLOW,V1,V2,V3,W,W60,W80,X,Y
44     *  ,ZETA1,ZETA2,Z1,Z2,Z3,RESAB6
45      REAL R1MACH
46      INTEGER J,KOUNT,L
47C
48      DIMENSION FV(19),W(9),X(3),Y(3),ZETA1(9),ZETA2(9),U(3),V(3)
49C
50C
51C            FIRST HOMOGENEOUS COORDINATES OF POINTS IN DEGREE-6
52C            AND DEGREE-8 FORMULA, TAKEN WITH MULTIPLICITY 3
53      DATA ZETA1(1),ZETA1(2),ZETA1(3),ZETA1(4),ZETA1(5),ZETA1(6),ZETA1(7
54     *  ),ZETA1(8),ZETA1(9)/0.5014265096581342E+00,
55     *  0.8738219710169965E+00,0.6365024991213939E+00,
56     *  0.5314504984483216E-01,0.8141482341455413E-01,
57     *  0.8989055433659379E+00,0.6588613844964797E+00,
58     *  0.8394777409957211E-02,0.7284923929554041E+00/
59C            SECOND HOMOGENEOUS COORDINATES OF POINTS IN DEGREE-6
60C            AND DEGREE-8 FORMULA, TAKEN WITH MUNLTIPLICITY 3
61      DATA ZETA2(1),ZETA2(2),ZETA2(3),ZETA2(4),ZETA2(5),ZETA2(6),ZETA2(7
62     *  ),ZETA2(8),ZETA2(9)/0.2492867451709329E+00,
63     *  0.6308901449150177E-01,0.5314504984483216E-01,
64     *  0.6365024991213939E+00,0.4592925882927229E+00,
65     *  0.5054722831703103E-01,0.1705693077517601E+00,
66     *  0.7284923929554041E+00,0.8394777409957211E-02/
67C           WEIGHTS OF MID-POINT OF TRIANGLE IN DEGREE-6
68C           RESP. DEGREE-8 FORMULAE
69      DATA W60/0.0E+00/
70      DATA W80/0.1443156076777862E+00/
71C           WEIGHTS IN DEGREE-6 AND DEGREE-8 RULE
72      DATA W(1),W(2),W(3),W(4),W(5),W(6),W(7),W(8),W(9)/
73     *  0.1167862757263407E+00,0.5084490637020547E-01,
74     *  0.8285107561839291E-01,0.8285107561839291E-01,
75     *  0.9509163426728497E-01,0.3245849762319813E-01,
76     *  0.1032173705347184E+00,0.2723031417443487E-01,
77     *  0.2723031417443487E-01/
78C
79C           LIST OF MAJOR VARIABLES
80C           ----------------------
81C           DJ      - AREA OF THE TRIANGLE
82C           DRESC   - APPROXIMATION TO INTEGRAL OF
83C                     ABS(F- IF/DJ)  OVER THE TRIANGLE
84C           RESAB6  - APPROXIMATION TO INTEGRAL OF
85C                     ABS(F) OVER THE TRIANGLE
86C           X       - CARTESIAN ABSCISSAE OF THE INTEGRATION
87C                     POINTS
88C           Y       - CARTESIAN ORDINATES OF THE INTEGRATION
89C                     POINTS
90C           FV      - FUNCTION VALUES
91C
92C           COMPUTE DEGREE-6 AND DEGREE-8 RESULTS FOR IF/DJ AND
93C           DEGREE-6 APPROXIMATION FOR ABS(F)
94C
95      EMACH = R1MACH(4)
96      UFLOW = R1MACH(1)
97      U1=U(1)
98      U2=U(2)
99      U3=U(3)
100      V1=V(1)
101      V2=V(2)
102      V3=V(3)
103      DJ = ABS(U1*V2-U2*V1-U1*V3+V1*U3+U2*V3-V2*U3)*0.5E+00
104      F0 = F((U1+U2+U3)/3.0E+00,(V1+V2+V3)/3.0E+00)
105      RES6 = F0*W60
106      RESAB6 = ABS(F0)*W60
107      FV(1) = F0
108      KOUNT = 1
109      RES8 = F0*W80
110      DO 50 J=1,9
111        Z1 = ZETA1(J)
112        Z2 = ZETA2(J)
113        Z3 = 1.0E+00-Z1-Z2
114        X(1) = Z1*U1+Z2*U2+Z3*U3
115        Y(1) = Z1*V1+Z2*V2+Z3*V3
116        X(2) = Z2*U1+Z3*U2+Z1*U3
117        Y(2) = Z2*V1+Z3*V2+Z1*V3
118        X(3) = Z3*U1+Z1*U2+Z2*U3
119        Y(3) = Z3*V1+Z1*V2+Z2*V3
120        IF(J.LE.4) THEN
121           F0 = 0.0E+00
122           DF0 = 0.0E+00
123           DO 10 L=1,3
124             KOUNT = KOUNT+1
125             FV(KOUNT) = F(X(L),Y(L))
126             F0 = F0+FV(KOUNT)
127             DF0 = DF0+ABS(FV(KOUNT))
128   10      CONTINUE
129           RES6 = RES6+F0*W(J)
130           RESAB6 = RESAB6+DF0*W(J)
131        ELSE
132           F0 = F(X(1),Y(1))+F(X(2),Y(2))+F(X(3),Y(3))
133           RES8 = RES8+F0*W(J)
134        ENDIF
135   50 CONTINUE
136C
137C           COMPUTE DEGREE-6 APPROXIMATION FOR THE INTEGRAL OF
138C           ABS(F-IF/DJ)
139C
140      DRESC = ABS(FV(1)-RES6)*W60
141      KOUNT = 2
142      DO 60 J=1,4
143        DRESC = DRESC+(ABS(FV(KOUNT)-RES6)+ABS(FV(KOUNT+1)-RES6)+ABS(
144     *  FV(KOUNT+2)-RES6))*W(J)
145        KOUNT = KOUNT+3
146   60 CONTINUE
147C
148C           COMPUTE DEGREE-6 AND DEGREE-8 APPROXIMATIONS FOR IF,
149C           AND ERROR ESTIMATE
150C
151      RES6 = RES6*DJ
152      RES8 = RES8*DJ
153      RESAB6 = RESAB6*DJ
154      DRESC = DRESC*DJ
155      EST = ABS(RES6-RES8)
156      IF(DRESC.NE.0.0E+00) EST = MAX(EST,DRESC*MIN(1.0E+00,(20.0E+00
157     *  *EST/DRESC)**1.5E+00))
158      IF(RESAB6.GT.UFLOW) EST = MAX(EMACH*RESAB6,EST)
159      RETURN
160      END
161