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