SUBROUTINE QUAD ( A, B, RESULT, K, EPSIL, NPTS, ICHECK, F ) c*********************************************************************72 c cc QUAD c C THIS SUBROUTINE ATTEMPTS TO CALCULATE THE INTEGRAL OF F(X) C OVER THE INVERVAL *A* TO *B* WITH RELATIVE ERROR NOT C EXCEEDING *EPSIL*. C THE RESULT IS OBTAINED USING A SEQUENCE OF 1, 3, 7, 15, 31, 63, C 127, AND 255 POINT INTERLACING FORMULAE (NO INTEGRAND C EVALUATIONS ARE WASTED) OF RESPECTIVE DEGREES 1, 5, 11, 23, C 47, 95, 191 AND 383. THE FORMULAE ARE BASED ON THE OPTIMAL C EXTENSIONS OF THE 3-POINT GAUSS FORMULA. DETAILS OF C THE FORMULAE ARE GIVEN IN 'THE OPTIMUM ADDITION OF POINTS C TO QUADRATURE FORMULAE' BY T.N.L. PATTERSON, MATHS. COMP. C VOL 22, 847-856, 1968. C *** INPUT *** C A LOWER LIMIT OF INTEGRATION. C B UPPER LIMIT OF INTEGRATION. C EPSIL RELATIVE ACCURACY REQUIRED. WHEN THE RELATIVE C DIFFERENCE OF TWO SUCCESSIVE FORMULAE DOES NOT C EXCEED *EPSIL* THE LAST FORMULA COMPUTED IS TAKEN C AS THE RESULT. C F F(X) IS THE INTEGRAND. C *** OUTPUT *** C RESULT THIS ARRAY, WHICH SHOULD BE DECLARED TO HAVE AT C LEAST 8 ELEMENTS, HOLDS THE RESULTS OBTAINED BY C THE 1, 3, 7, ETC. POINT FORMULAE. THE NUMBER OF C FORMULAE COMPUTED DEPENDS ON *EPSIL*. C K RESULT(K) HOLDS THE VALUE OF THE INTEGRAL TO THE C SPECIFIED RELATIVE ACCURACY. C NPTS NUMBER INTEGRAND EVALUATIONS. C ICHECK ON EXIT NORMALLY ICHECK = 0. HOWEVER, IF CONVERGENCE C TO THE ACCURACY REQUESTED IS NOT ACHIEVED, ICHECK = 1 C ON EXIT. C ABSCISSAE AND WEIGHTS OF QUADRATURE RULES ARE STACKED IN C ARRAY *P* IN THE ORDER IN WHICH THEY ARE NEEDED. c DIMENSION FUNCT(127), P(381), RESULT(8) DATA & P( 1),P( 2),P( 3),P( 4),P( 5),P( 6),P( 7), & P( 8),P( 9),P(10),P(11),P(12),P(13),P(14), & P(15),P(16),P(17),P(18),P(19),P(20),P(21), & P(22),P(23),P(24),P(25),P(26),P(27),P(28) / & 0.77459666924148337704E+00,0.55555555555555555556E+00, & 0.88888888888888888889E+00,0.26848808986833344073E+00, & 0.96049126870802028342E+00,0.10465622602646726519E+00, & 0.43424374934680255800E+00,0.40139741477596222291E+00, & 0.45091653865847414235E+00,0.13441525524378422036E+00, & 0.51603282997079739697E-01,0.20062852938698902103E+00, & 0.99383196321275502221E+00,0.17001719629940260339E-01, & 0.88845923287225699889E+00,0.92927195315124537686E-01, & 0.62110294673722640294E+00,0.17151190913639138079E+00, & 0.22338668642896688163E+00,0.21915685840158749640E+00, & 0.22551049979820668739E+00,0.67207754295990703540E-01, & 0.25807598096176653565E-01,0.10031427861179557877E+00, & 0.84345657393211062463E-02,0.46462893261757986541E-01, & 0.85755920049990351154E-01,0.10957842105592463824E+00 / DATA & P(29),P(30),P(31),P(32),P(33),P(34),P(35), & P(36),P(37),P(38),P(39),P(40),P(41),P(42), & P(43),P(44),P(45),P(46),P(47),P(48),P(49), & P(50),P(51),P(52),P(53),P(54),P(55),P(56) / & 0.99909812496766759766E+00,0.25447807915618744154E-02, & 0.98153114955374010687E+00,0.16446049854387810934E-01, & 0.92965485742974005667E+00,0.35957103307129322097E-01, & 0.83672593816886873550E+00,0.56979509494123357412E-01, & 0.70249620649152707861E+00,0.76879620499003531043E-01, & 0.53131974364437562397E+00,0.93627109981264473617E-01, & 0.33113539325797683309E+00,0.10566989358023480974E+00, & 0.11248894313318662575E+00,0.11195687302095345688E+00, & 0.11275525672076869161E+00,0.33603877148207730542E-01, & 0.12903800100351265626E-01,0.50157139305899537414E-01, & 0.42176304415588548391E-02,0.23231446639910269443E-01, & 0.42877960025007734493E-01,0.54789210527962865032E-01, & 0.12651565562300680114E-02,0.82230079572359296693E-02, & 0.17978551568128270333E-01,0.28489754745833548613E-01 / DATA & P(57),P(58),P(59),P(60),P(61),P(62),P(63), & P(64),P(65),P(66),P(67),P(68),P(69),P(70), & P(71),P(72),P(73),P(74),P(75),P(76),P(77), & P(78),P(79),P(80),P(81),P(82),P(83),P(84) / & 0.38439810249455532039E-01,0.46813554990628012403E-01, & 0.52834946790116519862E-01,0.55978436510476319408E-01, & 0.99987288812035761194E+00,0.36322148184553065969E-03, & 0.99720625937222195908E+00,0.25790497946856882724E-02, & 0.98868475754742947994E+00,0.61155068221172463397E-02, & 0.97218287474858179658E+00,0.10498246909621321898E-01, & 0.94634285837340290515E+00,0.15406750466559497802E-01, & 0.91037115695700429250E+00,0.20594233915912711149E-01, & 0.86390793819369047715E+00,0.25869679327214746911E-01, & 0.80694053195021761186E+00,0.31073551111687964880E-01, & 0.73975604435269475868E+00,0.36064432780782572640E-01, & 0.66290966002478059546E+00,0.40714410116944318934E-01, & 0.57719571005204581484E+00,0.44914531653632197414E-01, & 0.48361802694584102756E+00,0.48564330406673198716E-01 / DATA & P( 85),P( 86),P( 87),P( 88),P( 89),P( 90),P( 91), & P( 92),P( 93),P( 94),P( 95),P( 96),P( 97),P( 98), & P( 99),P(100),P(101),P(102),P(103),P(104),P(105), & P(106),P(107),P(108),P(109),P(110),P(111),P(112) / & 0.38335932419873034692E+00,0.51583253952048458777E-01, & 0.27774982202182431507E+00,0.53905499335266063927E-01, & 0.16823525155220746498E+00,0.55481404356559363988E-01, & 0.56344313046592789972E-01,0.56277699831254301273E-01, & 0.56377628360384717388E-01,0.16801938574103865271E-01, & 0.64519000501757369228E-02,0.25078569652949768707E-01, & 0.21088152457266328793E-02,0.11615723319955134727E-01, & 0.21438980012503867246E-01,0.27394605263981432516E-01, & 0.63260731936263354422E-03,0.41115039786546930472E-02, & 0.89892757840641357233E-02,0.14244877372916774306E-01, & 0.19219905124727766019E-01,0.23406777495314006201E-01, & 0.26417473395058259931E-01,0.27989218255238159704E-01, & 0.18073956444538835782E-03,0.12895240826104173921E-02, & 0.30577534101755311361E-02,0.52491234548088591251E-02 / DATA & P(113),P(114),P(115),P(116),P(117),P(118),P(119), & P(120),P(121),P(122),P(123),P(124),P(125),P(126), & P(127),P(128),P(129),P(130),P(131),P(132),P(133), & P(134),P(135),P(136),P(137),P(138),P(139),P(140) / & 0.77033752332797418482E-02,0.10297116957956355524E-01, & 0.12934839663607373455E-01,0.15536775555843982440E-01, & 0.18032216390391286320E-01,0.20357755058472159467E-01, & 0.22457265826816098707E-01,0.24282165203336599358E-01, & 0.25791626976024229388E-01,0.26952749667633031963E-01, & 0.27740702178279681994E-01,0.28138849915627150636E-01, & 0.99998243035489159858E+00,0.50536095207862517625E-04, & 0.99959879967191068325E+00,0.37774664632698466027E-03, & 0.99831663531840739253E+00,0.93836984854238150079E-03, & 0.99572410469840718851E+00,0.16811428654214699063E-02, & 0.99149572117810613240E+00,0.25687649437940203731E-02, & 0.98537149959852037111E+00,0.35728927835172996494E-02, & 0.97714151463970571416E+00,0.46710503721143217474E-02, & 0.96663785155841656709E+00,0.58434498758356395076E-02 / DATA & P(141),P(142),P(143),P(144),P(145),P(146),P(147), & P(148),P(149),P(150),P(151),P(152),P(153),P(154), & P(155),P(156),P(157),P(158),P(159),P(160),P(161), & P(162),P(163),P(164),P(165),P(166),P(167),P(168) / & 0.95373000642576113641E+00,0.70724899954335554680E-02, & 0.93832039777959288365E+00,0.83428387539681577056E-02, & 0.92034002547001242073E+00,0.96411777297025366953E-02, & 0.89974489977694003664E+00,0.10955733387837901648E-01, & 0.87651341448470526974E+00,0.12275830560082770087E-01, & 0.85064449476835027976E+00,0.13591571009765546790E-01, & 0.82215625436498040737E+00,0.14893641664815182035E-01, & 0.79108493379984836143E+00,0.16173218729577719942E-01, & 0.75748396638051363793E+00,0.17421930159464173747E-01, & 0.72142308537009891548E+00,0.18631848256138790186E-01, & 0.68298743109107922809E+00,0.19795495048097499488E-01, & 0.64227664250975951377E+00,0.20905851445812023852E-01, & 0.59940393024224289297E+00,0.21956366305317824939E-01, & 0.55449513263193254887E+00,0.22940964229387748761E-01 / DATA & P(169),P(170),P(171),P(172),P(173),P(174),P(175), & P(176),P(177),P(178),P(179),P(180),P(181),P(182), & P(183),P(184),P(185),P(186),P(187),P(188),P(189), & P(190),P(191),P(192),P(193),P(194),P(195),P(196) / & 0.50768775753371660215E+00,0.23854052106038540080E-01, & 0.45913001198983233287E+00,0.24690524744487676909E-01, & 0.40897982122988867241E+00,0.25445769965464765813E-01, & 0.35740383783153215238E+00,0.26115673376706097680E-01, & 0.30457644155671404334E+00,0.26696622927450359906E-01, & 0.25067873030348317661E+00,0.27185513229624791819E-01, & 0.19589750271110015392E+00,0.27579749566481873035E-01, & 0.14042423315256017459E+00,0.27877251476613701609E-01, & 0.84454040083710883710E-01,0.28076455793817246607E-01, & 0.28184648949745694339E-01,0.28176319033016602131E-01, & 0.28188814180192358694E-01,0.84009692870519326354E-02, & 0.32259500250878684614E-02,0.12539284826474884353E-01, & 0.10544076228633167722E-02,0.58078616599775673635E-02, & 0.10719490006251933623E-01,0.13697302631990716258E-01 / DATA & P(197),P(198),P(199),P(200),P(201),P(202),P(203), & P(204),P(205),P(206),P(207),P(208),P(209),P(210), & P(211),P(212),P(213),P(214),P(215),P(216),P(217), & P(218),P(219),P(220),P(221),P(222),P(223),P(224) / & 0.31630366082222647689E-03,0.20557519893273465236E-02, & 0.44946378920320678616E-02,0.71224386864583871532E-02, & 0.96099525623638830097E-02,0.11703388747657003101E-01, & 0.13208736697529129966E-01,0.13994609127619079852E-01, & 0.90372734658751149261E-04,0.64476204130572477933E-03, & 0.15288767050877655684E-02,0.26245617274044295626E-02, & 0.38516876166398709241E-02,0.51485584789781777618E-02, & 0.64674198318036867274E-02,0.77683877779219912200E-02, & 0.90161081951956431600E-02,0.10178877529236079733E-01, & 0.11228632913408049354E-01,0.12141082601668299679E-01, & 0.12895813488012114694E-01,0.13476374833816515982E-01, & 0.13870351089139840997E-01,0.14069424957813575318E-01, & 0.25157870384280661489E-04,0.18887326450650491366E-03, & 0.46918492424785040975E-03,0.84057143271072246365E-03 / DATA & P(225),P(226),P(227),P(228),P(229),P(230),P(231), & P(232),P(233),P(234),P(235),P(236),P(237),P(238), & P(239),P(240),P(241),P(242),P(243),P(244),P(245), & P(246),P(247),P(248),P(249),P(250),P(251),P(252) / & 0.12843824718970101768E-02,0.17864463917586498247E-02, & 0.23355251860571608737E-02,0.29217249379178197538E-02, & 0.35362449977167777340E-02,0.41714193769840788528E-02, & 0.48205888648512683476E-02,0.54778666939189508240E-02, & 0.61379152800413850435E-02,0.67957855048827733948E-02, & 0.74468208324075910174E-02,0.80866093647888599710E-02, & 0.87109650797320868736E-02,0.93159241280693950932E-02, & 0.98977475240487497440E-02,0.10452925722906011926E-01, & 0.10978183152658912470E-01,0.11470482114693874380E-01, & 0.11927026053019270040E-01,0.12345262372243838455E-01, & 0.12722884982732382906E-01,0.13057836688353048840E-01, & 0.13348311463725179953E-01,0.13592756614812395910E-01, & 0.13789874783240936517E-01,0.13938625738306850804E-01, & 0.14038227896908623303E-01,0.14088159516508301065E-01 / DATA & P(253),P(254),P(255),P(256),P(257),P(258),P(259), & P(260),P(261),P(262),P(263),P(264),P(265),P(266), & P(267),P(268),P(269),P(270),P(271),P(272),P(273), & P(274),P(275),P(276),P(277),P(278),P(279),P(280) / & 0.99999759637974846476E+00,0.69379364324103267170E-05, & 0.99994399620705437576E+00,0.53275293669780613125E-04, & 0.99976049092443204733E+00,0.13575491094922871973E-03, & 0.99938033802502358193E+00,0.24921240048299729402E-03, & 0.99874561446809511470E+00,0.38974528447328229322E-03, & 0.99780535449595727456E+00,0.55429531493037471492E-03, & 0.99651414591489027385E+00,0.74028280424450333046E-03, & 0.99483150280062100052E+00,0.94536151685852538246E-03, & 0.99272134428278861533E+00,0.11674841174299594077E-02, & 0.99015137040077015918E+00,0.14049079956551446427E-02, & 0.98709252795403406719E+00,0.16561127281544526052E-02, & 0.98351865757863272876E+00,0.19197129710138724125E-02, & 0.97940628167086268381E+00,0.21944069253638388388E-02, & 0.97473445975240266776E+00,0.24789582266575679307E-02 / DATA & P(281),P(282),P(283),P(284),P(285),P(286),P(287), & P(288),P(289),P(290),P(291),P(292),P(293),P(294), & P(295),P(296),P(297),P(298),P(299),P(300),P(301), & P(302),P(303),P(304),P(305),P(306),P(307),P(308) / & 0.96948465950245923177E+00,0.27721957645934509940E-02, & 0.96364062156981213252E+00,0.30730184347025783234E-02, & 0.95718821610986096274E+00,0.33803979910869203823E-02, & 0.95011529752129487656E+00,0.36933779170256508183E-02, & 0.94241156519108305981E+00,0.40110687240750233989E-02, & 0.93406843615772578800E+00,0.43326409680929828545E-02, & 0.92507893290707565236E+00,0.46573172997568547773E-02, & 0.91543758715576504064E+00,0.49843645647655386012E-02, & 0.90514035881326159519E+00,0.53130866051870565663E-02, & 0.89418456833555902286E+00,0.56428181013844441585E-02, & 0.88256884024734190684E+00,0.59729195655081658049E-02, & 0.87029305554811390585E+00,0.63027734490857587172E-02, & 0.85735831088623215653E+00,0.66317812429018878941E-02, & 0.84376688267270860104E+00,0.69593614093904229394E-02 / DATA & P(309),P(310),P(311),P(312),P(313),P(314),P(315), & P(316),P(317),P(318),P(319),P(320),P(321),P(322), & P(323),P(324),P(325),P(326),P(327),P(328),P(329), & P(330),P(331),P(332),P(333),P(334),P(335),P(336) / & 0.82952219463740140018E+00,0.72849479805538070639E-02, & 0.81462878765513741344E+00,0.76079896657190565832E-02, & 0.79909229096084140180E+00,0.79279493342948491103E-02, & 0.78291939411828301639E+00,0.82443037630328680306E-02, & 0.76611781930376009072E+00,0.85565435613076896192E-02, & 0.74869629361693660282E+00,0.88641732094824942641E-02, & 0.73066452124218126133E+00,0.91667111635607884067E-02, & 0.71203315536225203459E+00,0.94636899938300652943E-02, & 0.69281376977911470289E+00,0.97546565363174114611E-02, & 0.67301883023041847920E+00,0.10039172044056840798E-01, & 0.65266166541001749610E+00,0.10316812330947621682E-01, & 0.63175643771119423041E+00,0.10587167904885197931E-01, & 0.61031811371518640016E+00,0.10849844089337314099E-01, & 0.58836243444766254143E+00,0.11104461134006926537E-01 / DATA & P(337),P(338),P(339),P(340),P(341),P(342),P(343), & P(344),P(345),P(346),P(347),P(348),P(349),P(350), & P(351),P(352),P(353),P(354),P(355),P(356),P(357), & P(358),P(359),P(360),P(361),P(362),P(363),P(364) / & 0.56590588542365442262E+00,0.11350654315980596602E-01, & 0.54296566649831149049E+00,0.11588074033043952568E-01, & 0.51955966153745702199E+00,0.11816385890830235763E-01, & 0.49570640791876146017E+00,0.12035270785279562630E-01, & 0.47142506587165887693E+00,0.12244424981611985899E-01, & 0.44673539866202847374E+00,0.12443560190714035263E-01, & 0.42165768662616330006E+00,0.12632403643542078765E-01, & 0.39621280605761593918E+00,0.12810698163877361967E-01, & 0.37042208795007823014E+00,0.12978202239537399286E-01, & 0.34430734159943802278E+00,0.13134690091960152836E-01, & 0.31789081206847668318E+00,0.13279951743930530650E-01, & 0.29119514851824668196E+00,0.13413793085110098513E-01, & 0.26424337241092676194E+00,0.13536035934956213614E-01, & 0.23705884558982972721E+00,0.13646518102571291428E-01 / DATA & P(365),P(366),P(367),P(368),P(369),P(370),P(371), & P(372),P(373),P(374),P(375),P(376),P(377),P(378), & P(379),P(380),P(381) / & 0.20966523824318119477E+00,0.13745093443001896632E-01, & 0.18208649675925219825E+00,0.13831631909506428676E-01, & 0.15434681148137810869E+00,0.13906019601325461264E-01, & 0.12647058437230196685E+00,0.13968158806516938516E-01, & 0.98482396598119202090E-01,0.14017968039456608810E-01, & 0.70406976042855179063E-01,0.14055382072649964277E-01, & 0.42269164765363603212E-01,0.14080351962553661325E-01, & 0.14093886410782462614E-01,0.14092845069160408355E-01, & 0.14094407090096179347E-01 / ICHECK = 0 C CHECK FOR TRIVIAL CASE. IF ( A .EQ. B ) GO TO 70 C SCALE FACTORS. SUM = ( B + A ) / 2.0 DIFF = ( B - A ) / 2.0 C 1-POINT GAUSS. FZERO = F ( SUM ) RESULT(1) = 2.0 * FZERO * DIFF I = 0 IOLD = 0 INEW = 1 K = 2 ACUM = 0.0 GO TO 30 10 IF ( K .EQ. 8 ) GO TO 50 K = K + 1 ACUM = 0.0 C CONTRIBUTION FROM FUNCTION VALUES ALREADY COMPUTED. DO 20 J = 1, IOLD I = I + 1 ACUM = ACUM + P(I) * FUNCT(J) 20 CONTINUE C CONTRIBUTION FROM NEW FUNCTION VALUES. 30 IOLD = IOLD + INEW DO 40 J = INEW, IOLD I = I + 1 X = P(I) * DIFF FUNCT(J) = F ( SUM + X ) + F ( SUM - X ) I = I + 1 ACUM = ACUM + P(I) * FUNCT(J) 40 CONTINUE INEW = IOLD + 1 I = I + 1 RESULT(K) = ( ACUM + P(I) * FZERO ) * DIFF C CHECK FOR CONVERGENCE. IF ( ABS ( RESULT(K) - RESULT(K-1) ) - EPSIL * ABS ( RESULT(K) ) ) & 60, 60, 10 C CONVERGENCE NOT ACHIEVED. 50 ICHECK = 1 C NORMAL TERMINATION. 60 NPTS = INEW + IOLD RETURN C TRIVIAL CASE. 70 K = 2 RESULT(1) = 0.0 RESULT(2) = 0.0 NPTS = 0 RETURN END FUNCTION QSUB ( A, B, EPSIL, NPTS, ICHECK, RELERR, F ) c*********************************************************************72 c cc QSUB c C THIS FUNCTION ROUTINE PERFORMS AUTOMATIC INTEGRATION C OVER A FINITE INTERVAL USING THE BASIC INTEGRATION C ALGORITHM QUAD, TOGETHER WITH, IF NECESSARY, A NON- C ADAPTIVE SUBDIVISION PROCESS. C THE CALL TAKES THE FORM C QSUB(A,B,EPSIL,NPTS,ICHECK,RELERR,F) C AND CAUSES F(X) TO BE INTEGRATED OVER (A,B) WITH RELATIVE C ERROR HOPEFULLY NOT EXCEEDING EPSIL. SHOULD QUAD CONVERGE C (ICHECK=0) THEN QSUB WILL RETURN THE VALUE OBTAINED BY IT; C OTHERWISE SUBDIVISION WILL BE INVOKED AS A RESCUE C OPERATION IN A NON-ADAPTIVE MANNER. THE ARGUMENT RELERR C GIVES A CRUDE ESTIMATE OF THE ACTUAL RELATIVE ERROR C OBTAINED. C THE SUBDIVISION STRATEGY IS AS FOLLOWS: C LET THE INTERVAL (A,B) BE DIVIDED INTO 2**N PANELS AT STEP C N OF THE SUBDIVISION PROCESS. QUAD IS APPLIED FIRST TO C THE SUBDIVIDED INTERVAL ON WHICH QUAD LAST FAILED TO C CONVERGE, AND IF CONVERGENCE IS NOW ACHIEVED THE REMAINING C PANELS ARE INTEGRATED. SHOULD A CONVERGENCE FAILURE OCCUR C ON ANY PANEL, THE INTEGRATION AT THAT POINT IS TERMINATED C AND THE PROCEDURE REPEATED WITH N INCREASED BY 1. THE C STRATEGY ENSURES THAT POSSIBLY DELINQUENT INTERVALS ARE C EXAMINED BEFORE WORK, WHICH LATER MIGHT HAVE TO BE C DISCARDED, IS INVESTED IN WELL-BEHAVED PANELS. THE C PROCESS IS COMPLETE WHEN NO CONVERGENCE FAILURE OCCURS ON C ANY PANEL AND THE SUM OF THE RESULTS OBTAINED BY QUAD ON C EACH PANEL IS TAKEN AS THE VALUE OF THE INTEGRAL. C THE PROCESS IS VERY CAUTIOUS IN THAT THE SUBDIVISION OF C THE INTERVAL (A,B) IS UNIFORM, THE FINENESS OF WHICH IS C CONTROLLED BY THE SUCCESS OF QUAD. IN THIS WAY, IT IS C RATHER DIFFICULT FOR A SPURIOUS CONVERGENCE TO SLIP C THROUGH. C THE CONVERGENCE CRITERION OF QUAD IS SLIGHTLY RELAXED C IN THAT A PANEL IS DEEMED TO HAVE BEEN SUCCESSFULLY C INTEGRATED IF EITHER QUAD CONVERGES OR THE ESTIMATED C ABSOLUTE ERROR COMMITTED ON THIS PANEL DOES NOT EXCEED C EPSIL TIMES THE ESTIMATED ABSOLUTE VALUE OF THE INTEGRAL C OVER (A,B). THIS RELAXATION IS TO TRY TO TAKE ACCOUNT OF C A COMMON SITUATION WHERE ONE PARTICULAR PANEL CAUSES C SPECIAL DIFFICULTY, PERHAPS DUE TO SINGULARITY OF SOME C TYPE. IN THIS CASE, QUAD COULD OBTAIN NEARLY EXACT C ANSWERS ON ALL OTHER PANELS, AND SO THE RELATIVE ERROR FOR C THE TOTAL INTEGRATION WOULD BE ALMOST ENTIRELY DUE TO THE C DELINQUENT PANEL. WITHOUT THIS CONDITION, THE COMPUTATION C MIGHT CONTINUE DESPITE THE REQUESTED RELATIVE ERROR BEING C ACHIEVED. C THE OUTCOME OF THE INTEGRATION IS INDICATED BY ICHECK. C ICHECK = 0 - CONVERGENCE OBTAINED WITHOUT INVOKING C SUBDIVISION. THIS CORRESPONDS TO THE C DIRECT USE OF QUAD. C ICHECK = 1 - RESULT OBTAINED AFTER INVOKING SUBDIVISION. C ICHECK = 2 - AS FOR ICHECK=1 BUT AT SOME POINT THE C RELAXED CONVERGENCE CRITERION WAS USED. C THE RISK OF UNDERESTIMATING THE RELATIVE C ERROR WILL BE INCREASED. IF NECESSARY, C CONFIDENCE MAY BE RESTORED BY CHECKING C EPSIL AND RELERR FOR A SERIOUS DISCREPANCY. C ICHECK NEGATIVE C IF DURING THE SUBDIVISION PROCESS THE C ALLOWED UPPER LIMIT ON THE NUMBER OF PANELS C THAT MAY BE GENERATED (PRESENTLY 4096) IS C REACHED, A RESULT IS OBTAINED WHICH MAY BE C UNRELIABLE BY CONTINUING THE INTEGRATION C WITHOUT FURTHER SUBDIVISION IGNORING C CONVERGENCE FAILURES. THIS OCCURRENCE IS C FLAGGED BY RETURNING ICHECK WITH NEGATIVE C SIGN. C THE RELIABILITY OF THE ALGORITHM WILL DECREASE FOR LARGE C VALUES OF EPSIL. IT IS RECOMMENDED THAT EPSIL SHOULD C GENERALLY BE LESS THAN ABOUT 0.001. DIMENSION RESULT(8) INTEGER BAD, OUT LOGICAL RHS EXTERNAL F DATA NMAX / 4096 / CALL QUAD ( A, B, RESULT, K, EPSIL, NPTS, ICHECK, F ) QSUB = RESULT(K) RELERR = 0.0 IF ( QSUB .NE. 0.0 ) RELERR = & ABS ( ( RESULT(K) - RESULT(K-1) ) / QSUB ) C CHECK IF SUBDIVISION IS NEEDED. IF ( ICHECK .EQ. 0 ) RETURN C SUBDIVIDE. ESTIM = ABS ( QSUB * EPSIL ) IC = 1 RHS = .FALSE. N = 1 H = B - A BAD = 1 10 QSUB = 0.0 RELERR = 0.0 H = H * 0.5 N = N + N C INTERVAL (A,B) DIVIDED INTO N EQUAL SUBINTERVALS. C INTEGRATE OVER SUBINTERVALS BAD TO (BAD+1) WHERE TROUBLE C HAS OCCURRED. M1 = BAD M2 = BAD + 1 OUT = 1 GO TO 50 C INTEGRATE OVER SUBINTERVALS 1 TO (BAD-1). 20 M1 = 1 M2 = BAD - 1 RHS = .FALSE. OUT = 2 GO TO 50 C INTEGRATE OVER SUBINTERVALS (BAD+2) TO N. 30 M1 = BAD + 2 M2 = N OUT = 3 GO TO 50 C SUBDIVISION RESULT. 40 ICHECK = IC RELERR = RELERR / ABS ( QSUB ) RETURN C INTEGRATE OVER SUBINTERVALS M1 TO M2. 50 IF ( M1 .GT. M2 ) GO TO 90 DO 80 JJ = M1, M2 J = JJ C EXAMINE FIRST THE LEFT OR RIGHT HALF OF THE SUBDIVIDED C TROUBLESOME INTERVAL DEPENDING ON THE OBSERVED TREND. IF ( RHS ) J = M2 + M1 - JJ ALPHA = A + H * ( J - 1 ) BETA = ALPHA + H CALL QUAD ( ALPHA, BETA, RESULT, M, EPSIL, NF, ICHECK, F ) COMP = ABS ( RESULT(M) - RESULT(M-1) ) NPTS = NPTS + NF IF ( ICHECK .NE. 1 ) GO TO 70 IF ( COMP .LE. ESTIM ) GO TO 100 C SUBINTERVAL J HAS CAUSED TROUBLE. C CHECK IF FURTHER SUBDIVISION SHOULD BE CARRIED OUT. IF ( N .EQ. NMAX ) GO TO 60 BAD = 2 * J - 1 RHS = .FALSE. IF ( ( J - 2 * ( J / 2 ) ) .EQ. 0 ) RHS = .TRUE. GO TO 10 60 IC = -IABS ( IC ) 70 QSUB = QSUB + RESULT ( M ) 80 CONTINUE RELERR = RELERR + COMP 90 GO TO ( 20, 30, 40 ), OUT C RELAXED CONVERGENCE 100 IC = ISIGN ( 2, IC ) GO TO 70 END FUNCTION QSUBA ( A, B, EPSIL, NPTS, ICHECK, RELERR, F ) c*********************************************************************72 c cc QSUBA c C THIS FUNCTION ROUTINE PERFORMS AUTOMATIC INTEGRATION C OVER A FINITE INTERVAL USING THE BASIC INTEGRATION C ALGORITHM QUAD, TOGETHER WITH, IF NECESSARY, AN ADAPTIVE C SUBDIVISION PROCESS. IT IS GENERALLY MORE EFFICIENT THAN C THE NON-ADAPTIVE ALGORITHM QSUB BUT IS LIKELY TO BE LESS C RELIABLE (SEE COMP.J., 14, 189, 1971 ). C THE CALL TAKES THE FORM C QSUBA(A,B,EPSIL,NPTS,ICHECK,RELERR,F) C AND CAUSES F(X) TO BE INTEGRATED OVER (A,B) WITH RELATIVE C ERROR HOPEFULLY NOT EXCEEDING EPSIL. SHOULD QUAD CONVERGE C (ICHECK=0) THEN QSUBA WILL RETURN THE VALUE OBTAINED BY IT; C OTHERWISE SUBDIVISION WILL BE INVOKED AS A RESCUE C OPERATION IN AN ADAPTIVE MANNER. THE ARGUMENT RELERR GIVES C A CRUDE ESTIMATE OF THE ACTUAL RELATIVE ERROR OBTAINED. C THE SUBDIVISION STRATEGY IS AS FOLLOWS: C AT EACH STAGE OF THE PROCESS, AN INTERVAL IS PRESENTED FOR C SUBDIVISION. (INITIALLY THIS WILL BE THE WHOLE INTERVAL C (A,B)). THE INTERVAL IS HALVED AND QUAD APPLIED TO EACH C SUBINTERVAL. SHOULD QUAD FAIL ON THE FIRST SUBINTERVAL, C THE SUBINTERVAL IS STACKED FOR FUTURE SUBDIVISION AND THE C SECOND SUBINTERVAL IMMEDIATELY EXAMINED. SHOULD QUAD FAIL C ON THE SECOND SUBINTERVAL, THE SUBINTERVAL IS C IMMEDIATELY SUBDIVIDED AND THE WHOLE PROCESS REPEATED. C EACH TIME A CONVERGED RESULT IS OBTAINED, IT IS C ACCUMULATED AS THE PARTIAL VALUE OF THE INTEGRAL. WHEN C QUAD CONVERGES ON BOTH SUBINTERVALS, THE INTERVAL LAST C STACKED IS CHOSEN NEXT FOR SUBDIVISION AND THE PROCESS C REPEATED. A SUBINTERVAL IS NOT EXAMINED AGAIN ONCE A C CONVERGED RESULT IS OBTAINED FOR IT, SO THAT A SPURIOUS C CONVERGENCE IS MORE LIKELY TO SLIP THROUGH THAN FOR THE C NON-ADAPTIVE ALGORITHM QSUB. C THE CONVERGENCE CRITERION OF QUAD IS SLIGHTLY RELAXED C IN THAT A PANEL IS DEEMED TO HAVE BEEN SUCCESSFULLY C INTEGRATED IF EITHER QUAD CONVERGES OR THE ESTIMATED C ABSOLUTE ERROR COMMITTED ON THIS PANEL DOES NOT EXCEED C EPSIL TIMES THE ESTIMATED ABSOLUTE VALUE OF THE INTEGRAL C OVER (A,B). THIS RELAXATION IS TO TRY TO TAKE ACCOUNT OF C A COMMON SITUATION WHERE ONE PARTICULAR PANEL CAUSES C SPECIAL DIFFICULTY, PERHAPS DUE TO SINGULARITY OF SOME C TYPE. IN THIS CASE, QUAD COULD OBTAIN NEARLY EXACT C ANSWERS ON ALL OTHER PANELS, AND SO THE RELATIVE ERROR FOR C THE TOTAL INTEGRATION WOULD BE ALMOST ENTIRELY DUE TO THE C DELINQUENT PANEL. WITHOUT THIS CONDITION, THE COMPUTATION C MIGHT CONTINUE DESPITE THE REQUESTED RELATIVE ERROR BEING C ACHIEVED. C THE OUTCOME OF THE INTEGRATION IS INDICATED BY ICHECK. C ICHECK = 0 - CONVERGENCE OBTAINED WITHOUT INVOKING C SUBDIVISION. THIS WOULD CORRESPOND TO THE C DIRECT USE OF QUAD. C ICHECK = 1 - RESULT OBTAINED AFTER INVOKING SUBDIVISION. C ICHECK = 2 - AS FOR ICHECK=1 BUT AT SOME POINT THE C RELAXED CONVERGENCE CRITERION WAS USED. C THE RISK OF UNDERESTIMATING THE RELATIVE C ERROR WILL BE INCREASED. IF NECESSARY, C CONFIDENCE MAY BE RESTORED BY CHECKING C EPSIL AND RELERR FOR A SERIOUS DISCREPANCY. C ICHECK NEGATIVE C IF DURING THE SUBDIVISION PROCESS THE STACK C OF DELINQUENT INTERVALS BECOMES FULL (IT IS C PRESENTLY SET TO HOLD AT MOST 100 NUMBERS) C A RESULT IS OBTAINED BY CONTINUING THE C INTEGRATION IGNORING CONVERGENCE FAILURES C WHICH CANNOT BE ACCOMMODATED BY THE STACK. C THIS OCCURRENCE IS FLAGGED BY RETURNING C ICHECK WITH NEGATIVE SIGN. C THE RELIABILITY OF THE ALGORITHM WILL DECREASE FOR LARGE C VALUES OF EPSIL. IT IS RECOMMENDED THAT EPSIL SHOULD C GENERALLY BE LESS THAN ABOUT 0.001. DIMENSION RESULT(8), STACK(100) EXTERNAL F DATA ISMAX / 100 / CALL QUAD ( A, B, RESULT, K, EPSIL, NPTS, ICHECK, F ) QSUBA = RESULT(K) RELERR = 0.0 IF ( QSUBA .NE. 0.0 ) & RELERR = ABS ( ( RESULT(K) - RESULT(K-1) ) / QSUBA ) C CHECK IF SUBDIVISION IS NEEDED. IF ( ICHECK .EQ. 0 ) RETURN C SUBDIVIDE. ESTIM = ABS ( QSUBA * EPSIL ) RELERR = 0.0 QSUBA = 0.0 IS = 1 IC = 1 SUB1 = A SUB3 = B 10 SUB2 = ( SUB1 + SUB3 ) * 0.5 CALL QUAD ( SUB1, SUB2, RESULT, K, EPSIL, NF, ICHECK, F ) NPTS = NPTS + NF COMP = ABS ( RESULT(K) - RESULT(K-1) ) IF ( ICHECK .EQ. 0 ) GO TO 30 IF ( COMP .LE. ESTIM ) GO TO 70 IF ( IS .GE. ISMAX ) GO TO 20 C STACK SUBINTERVAL (SUB1,SUB2) FOR FUTURE EXAMINATION. STACK(IS) = SUB1 IS = IS + 1 STACK(IS) = SUB2 IS = IS + 1 GO TO 40 20 IC = -IABS ( IC ) 30 QSUBA = QSUBA + RESULT(K) RELERR = RELERR + COMP 40 CALL QUAD ( SUB2, SUB3, RESULT, K, EPSIL, NF, ICHECK, F ) NPTS = NPTS + NF COMP = ABS ( RESULT(K) - RESULT(K-1) ) IF ( ICHECK .EQ. 0 ) GO TO 50 IF ( COMP .LE. ESTIM ) GO TO 80 C SUBDIVIDE INTERVAL (SUB2,SUB3) SUB1 = SUB2 GO TO 10 50 QSUBA = QSUBA + RESULT(K) RELERR = RELERR + COMP IF ( IS .EQ. 1 ) GO TO 60 C SUBDIVIDE THE DELINQUENT INVERVAL LAST STACKED. IS = IS - 1 SUB3 = STACK(IS) IS = IS - 1 SUB1 = STACK(IS) GO TO 10 C SUBDIVISION RESULT. 60 ICHECK = IC RELERR = RELERR / ABS ( QSUBA ) RETURN C RELAXED CONVERGENCE. 70 IC = ISIGN ( 2, IC ) GO TO 30 80 IC = ISIGN ( 2, IC ) GO TO 50 END