SUBROUTINE SINTA (ANSWER,WORK,IFLAG) c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. c ALL RIGHTS RESERVED. c Based on Government Sponsored Research NAS7-03001. c>> 2009-11/03 SINTA Krogh Set initial value for a few variables. c>> 2009-10-17 SINTA Krogh AACUM now set with ACUM on coalescing nodes. c>> 2008-02-23 SINTA Krogh/Snyder Misc. Changes c>> 2001-09-11 SINTA Krogh Fixed so really small intervals not picked. c>> 1996-04-27 SINTA Krogh Changes to use .C. and C%%. C>> 1996-03-31 SINTA Krogh Removed unused variable in common. c>> 1996-03-30 SINTA Krogh Change specific intrinsics to generics. c>> 1995-11-28 SINTA Krogh Converted from SFTRAN to Fortran 77. C>> 1995-10-24 SINTA Krogh Removed blanks in numbers for C conversion. C>> 1994-11-14 SINTA Krogh Declared all vars. c>> 1994-10-19 SINTA Krogh Changes to use M77CON c>> 1994-09-01 SINTA Snyder don't give round-off msg if err << eps C>> 1994-08-22 SINTA Krogh -- Modified data for C conversion. c>> 1994-08-19 SINTA Snyder fix divide by 0 if ERR=0 and don't believe. c>> 1994-08-15 SINTA Snyder convert MAX(A,B,C) to MAX(A,MAX(B,C)) c>> 1994-07-07 SINTA Snyder set up for CHGTYP. C>> 1994-06-29 SINTA Snyder suppress discontinuity msg if one jump. C>> 1994-06-22 SINTA Snyder don't use RE when reducing stepsize. C>> 1993-05-26 SINTA Krogh -- Added data for C conversion. Needs edit. C>> 1993-05-18 SINTA Krogh -- Changed "END" to "END PROGRAM" C>> 1992-03-03 SINTA Added call to SINTO for an error message. C>> 1991-09-20 SINTA Krogh converted '(1)' dimensioning to '(*)'. C>> 1989-11-08 SINTA Snyder Revise coalesced abscissa test for CRAY C>> 1988-05-17 SINTA Snyder Initial code. c--S replaces "?": ?INT, ?INTA, ?INTC, ?INTDL, ?INTDU, ?INTEC c--& ?INTF, ?INTM, ?INTMA, ?INTNC, ?INTNS, ?INTO, ?INTSM, ?INT1 C C THIS IS THE INTEGRATION PROGRAM FOR SINT. C THE RESULT IS OBTAINED USING QUADRATURE FORMULAE DUE TO C T. N. L. PATTERSON, MATHEMATICS OF COMPUTATION, VOLUME 22, C PAGES 847-856, 1968. C C ***** WARNING *********************************************** C C THE RELIABILITY AND EFFICIENCY OF THIS PROGRAM ARE STRONGLY C INFLUENCED BY DISCONTINUITIES IN THE FUNCTION OR IT'S C DERIVATIVES, INTEGRABLE SINGULARITIES ON THE REAL AXIS, C NON-INTEGRABLE SINGULARITIES NEAR THE REGION, AND POLES C NOT ON THE REAL AXIS. THE EFFICIENCY AND RELIABILITY C OF INTEGRATING SUCH FUNCTIONS MAY BE GREATLY IMPROVED BY MANUALLY C SUBDIVIDING THE INTERVAL OF INTEGRATION AT THE DISCONTINUITY, C SINGULARITY OR REAL PART OF THE POLE, AND SUMMING THE ANSWERS. C A CHANGE OF VARIABLE TO ELIMINATE OR REDUCE THE STRENGTH OF THE C SINGULARITY WILL SIGNIFICANTLY IMPROVE PERFORMANCE. C C ***** FORMAL ARGUMENTS ************************************* C C ANSWER IS THE COMPUTED INTEGRAL IF IFLAG .LT. 0. IF IFLAG = 6, C ANSWER IS THE APPROXIMATE LOCATION OF A SINGULARITY WHICH C APPEARS NOT TO BE INTEGRABLE. IF IFLAG = 0, ANSWER IS C USED TO PASS A FUNCTION VALUE FROM THE USER PROGRAM TO SINTA. REAL ANSWER C WORK CONTAINS THE ESTIMATED ABSOLUTE ERROR IF IFLAG .LT. 0. C IF IFLAG = 0, WORK IS USED TO PASS AN ABSCISSA TO THE C USER PROGRAM (FOR REVERSE COMMUNICATION). REAL WORK(*) C IFLAG IS USED TO CONTROL REVERSE COMMUNICATION AND TO RETURN C STATUS INFORMATION TO THE USER. VALUES ARE C 0 - INTEGRAND VALUE NEEDED. COMPUTE ANSWER=INTEGRAND C EVALUATED AT WORK(1), ..., WORK(NDIM). C -1 - NORMAL TERMINATION, EITHER THE ABSOLUTE OR THE RELATIVE C ERROR CRITERIA ARE SATISFIED. C -2 - NORMAL TERMINATION, NEITHER THE ABSOLUTE NOR RELATIVE C ERROR CRITERIA ARE SATISFIED, BUT THE RELATIVE TO THE C OBTAINABLE ERROR CRITERION IS SATISFIED. C -3 - NORMAL TERMINATION, BUT NONE OF THE ERROR CRITERIA ARE C SATISFIED. C 4 - BAD IOPT VALUE (SEE SINT1 OR SINTM). C 5 - TOO MANY FUNCTION VALUES NEEDED. C KDIM+5 - APPARENT NON-INTEGRABLE SINGULARITY NEAR ANSWER. C KDIM IS CURRENT DIMENSION - 1 IN 1 DIMENSIONAL CASE. INTEGER IFLAG(*) C C ***** EXTERNAL REFERENCES *********************************** C C SINTDL TO FORM AND ANALYZE DIFFERENCE LINES. C C SINTDU TO UPDATE DIFFERENCE LINES. C C SINTF PROVIDE THE INTEGRAND. SINTF IS EXPLAINED IN C EXPLAINED IN SINT1 AND SINTM. C C SINTNS TO INCREASE AND REDUCE NSUB, THE DEGREE OF THE INDEPENDENT C VARIABLE TRANSFORMATION. C C SINTO USED FOR PRINTING OUTPUT. C C SINTSM USED TO CALCULATE THE MINIMUM STEPSIZE. EXTERNAL SINTSM REAL SINTSM INTEGER NSRB,NSRA,NSIB,NSIA,MAXK,NSUBMX,KONVRG,IC1,IC2,JUMPS C C ***** PARAMETERS ***************************************** C C NSIA AS FOR NSRA, BUT NSUB IS TO BE INCREASED. C NSIB AS FOR NSRB, BUT NSUB IS TO BE INCREASED. C NSRA IS AN ARGUMENT OF SINTNS THAT MEANS REDUCE NSUB, BUT ONLY C ALOCAL NEEDS TO BE RECALCULATED. C NSRB AS FOR NSRA, BUT BOTH ALOCAL AND BLOCAL MUST BE RECALCULATED. C PARAMETER (NSRB=1, NSRA=NSRB+1, NSIB=NSRA+1, NSIA=NSIB+1) C C ***** INTERNAL AND COMMON VARIABLES ********************** C C AACUM IS THE INTEGRAL OF THE ABSOLUTE VALUE OF THE INTEGRAND OVER C CURRENT PANEL. C ABSCIS IS THE ABSCISSA TO BE TRANSFORMED TO USER COORDINATES. C SEE NSUB. C ABSDIF IS THE ABSOLUTE VALUE OF ONE HALF OF THE LENGTH OF THE PANEL C (SMALLEST SUBDIVISION) CURRENTLY BEING INTEGRATED. C ACUM IS THE ESTIMATE OF THE INTEGRAL OVER THE CURRENT PANEL. C AIMFOR IS AN INDEX OF THE FORMULA TO AIM FOR AFTER A PANEL HAS BEEN C SUCCESSFULLY INTEGRATED. THE CURRENT VALUE OF AIMFOR C SPECIFIES A 63-POINT FORMULA. INTEGER AIMFOR C AINIT IS THE INITIAL LOWER REGION BOUNDARY. C ALOCAL IS THE 'LEFT' BOUNDARY OF THE PANEL BEING INTEGRATED. REAL ALOCAL C ALPHA IS USED DURING DETERMINATION OF THE NEXT STEP SIZE AS AN C INDICATOR OF THE PERFORMANCE OF THE INTEGRATION OF THE C PREVIOUS STEP, WITH RESPECT TO THE COMMITTED AND PERMITTED C ERROR. REAL ALPHA C BETA IS USED DURING COMPUTATION OF ABSCISSAE USED BY THE C QUADRATURE FORMULA AND DURING COMPUTATION OF THE NEXT STEP. REAL BETA C BINIT IS THE INITIAL UPPER REGION BOUNDARY. C BLOCAL IS THE 'RIGHT' BOUNDARY OF THE PANEL CURRENTLY BEING C INTEGRATED. REAL BLOCAL C COUNT IS USED DURING A SEARCH TO DISTINGUISH INTEGRABLE FROM C NON-INTEGRABLE SINGULARITIES. C DELMIN IS THE MINIMUM STEPSIZE ALLOWED, LEST THE ABSCISSAE COALESCE. C DELTA IS THE ABSOLUTE VALUE OF THE CURRENT STEP SIZE. C DID1 IS INITIALLY FALSE, AND SET TRUE WHEN ONE PANEL HAS C BEEN SUCCESSFULLY INTEGRATED. C DIFF IS ONE HALF OF THE SIGNED CURRENT STEP SIZE. C DISCF IS NON-ZERO WHEN A DISCONTINUITY IS DETECTED. IF DISCF IS C POSITIVE AND IPRINT IS NOT ZERO AFTER SUCCESSFULLY C INTEGRATING A PANEL IN PART = 1, THE BOUNDARIES FOR THE C PANEL ARE PRINTED. IN ADDITION, IF ABS(DISCF) IS 2 AFTER A C PANEL IS INTEGRATED, THE REST OF THE PART IS USED FOR THE C NEXT PANEL. C DISCHK IS SET TO MAX(DISCHK,1) WHEN A JUMP IS SEEN IN THE DIFFERENCE C LINE, BUT THE NEXT FORMULA IS TRIED, OR WHEN A DISCONTINUITY C IS DETECTED. IF THERE ARE JUMPS IN THE DIFFERENCE LINE WHICH C SUBSEQUENTLY VANISH, DISCHK IS SET TO 2. WHEN A PANEL IS C ACCEPTED, DISCHK IS DECREMENTED AS FAR AS ZERO. DISCHK IS C SET TO -1 WHEN NSUB IS CHANGED. DISCHK IS SET TO -2 WHEN A C CAUTION POINT IS ENCOUNTERED. IF DISCHK IS POSITIVE, THE C DIFFERENCE LINES MUST BE FORMED EVERY QUADRATURE STEP, AND C THE END POINTS MUST BE ADDED BEFORE THE ANSWER CAN BE C ACCEPTED. IF DISCHK IS -1, THE DIFFERENCE LINE MUST BE C FORMED BEFORE THE ANSWER CAN BE ACCEPTED. IF DISCHK IS -2 C THE DIFFERENCE LINES MUST BE FORMED AND THE END POINTS ADDED C BEFORE THE ANSWER CAN BE ACCEPTED. C DISCX IS THE APPROXIMATE LOCATION OF A DISCONTINUITY. C EDUE2A, EDUE2B ARE THE ERRORS IN THE INTEGRAL DUE TO IMPRECISE LIMITS C A AND B RESPECTIVELY. C END STORES THE 'RIGHT' BOUNDARIES OF MAJOR SUBDIVISIONS. C ENDPHI CONTAINS THE END OF THE DIFFERENCE LINE ON WHICH THERE CAN BE C A JUMP. ENDPHI(1) = 1 BECAUSE THERE CAN ONLY BE A JUMP IN PHI C AT 1. ENDPHI(2) = LENDT BECAUSE THERE CAN ONLY BE A JUMP IN C PHIT AT LENDT. ENDPHI IS USED ONLY WHEN CHECKING WHETHER THE C ERROR IN THE FUNCTION AT THE END OF A PANEL IS SO LARGE THAT C A JUMP IS NOT NECESSARILY BELIEVABLE. INTEGER ENDPHI(2) C ENDPTS IS USED TO CONTROL ADDITION OF THE END POINTS TO THE C DIFFERENCE LINES: C 1. FUNCTION NEED NOT BE EVALUATED AT A. C 2. FUNCTION NEED NOT BE EVALUATED AT B. C 3. FUNCTION NEED NOT BE EVALUATED AT EITHER A AND B. C EP IS THE PREVIOUS VALUE OF ERRC. C EPNOIZ IS THE MAXIMUM ERROR THRESHOLD FOR NOISE DETECTION. C EPS IS THE ERROR ALLOWED ON THE CURRENT PANEL. C EPSMAX IS THE MOST STRINGENT TOLERANCE ALLOWED ON THE CURRENT PANEL. C SEE XEPS. C EPSMIN IS AN ESTIMATE OF THE SMALLEST ERROR POSSIBLE ON THE CURRENT C PANEL. C EPSO IS THE ORIGINAL ERROR TOLERANCE. C EPSR IS THE REMAINING ERROR TOLERANCE. C EPSS IS A STORAGE CELL FOR THE MAXIMUM OF EPS AND EPSMAX. EPS IS C CHANGED AFTER EPSS IS COMPUTED. C ERR IS THE EXTRAPOLATED ERROR COMMITTED DURING THE CURRENT C INTEGRATION STEP. C ERRAT ERRAT(I) CONTAINS THE ERROR IN FAT(I) IF FATAS(I), I=1 OR 2. C ERRC IS THE ERROR COMMITTED DURING THE CURRENT INTEGRATION STEP, C FOR PURPOSES OF DETERMINING THE NEXT STEPSIZE. C ERRCF IS USED TO COMPUTE A MINIMUM VALUE FOR THE ERROR C WHEN K.LE.4, USING THE DIFFERENCE LINES. REAL ERRCF(3) C ERRF IS THE ESTIMATED ERROR IN THE INTEGRAND. WHEN KDIM=1, ERRF C IS OBTAINED FROM WORK(FEA) (IF FEA .NE. 0). WHEN KDIM.GT.1 C ERRF IS THE ESTIMATED ERROR IN THE INNER INTEGRAL. C ERRI IS THE ERROR ESTIMATED BY THE DIFFERENCE BETWEEN TWO C SUCCESSIVE INTEGRATION FORMULAE APPLIED TO THE CURRENT PANEL. C ERRINA, ERRINB ARE THE ERRORS THE USER DECLARES BY WAY OF AN OPTION C TO BE PRESENT IN THE LIMITS A AND B RESPECTIVELY. C ERRT IS THE TOTAL ERROR COMMITTED ON A MAJOR SUBDIVISION. C ESOLD IS THE PREVIOUS ERROR ESTIMATED ACROSS A JUMP BY THE SEARCH. C EXTRA IS ANY EXTRA AMOUNT ADDED TO THE ERROR BECAUSE OF OUR C SUSPICIOUS NATURE. EXTRA IS BASED UPON HIGH-ORDER C DIFFERENCES IF THERE ARE NO JUMPS. C FAIL INDICATES TO THE STEPSIZE SELECTION ROUTINE THAT INTEGRATION C WAS NOT SUCCESSFUL. FAIL IS ALSO SET IF A JUMP IN THE C DIFFERENCE LINE VANISHES WHEN ACCURATE LOCATION IS ATTEMPTED. C FAT IS A VECTOR EQUIVALENT TO FATA,FATB. C FATA SAVES THE FUNCTION VALUE AT ALOCAL. REAL FATA C FATAS IS TRUE WHEN FATA AND ERRAT(1) CONTAIN VALID DATA. LOGICAL FATAS C FATB SAVES THE FUNCTION VALUE AT BLOCAL. REAL FATB C FATBS IS TRUE WHEN FATB AND ERRAT(2) CONTAIN VALID DATA. LOGICAL FATBS C FATS IS A VECTOR EQUIVALENT TO FATAS,FATBS. C FEA IS THE INDEX IN WORK OF THE ABSOLUTE ERROR COMMITTED C COMPUTING F. WORK(FEA) MAY BE CHANGED DURING THE INTEGRATION. C FEA APPLIES ONLY TO THE INNERMOST INTEGRATION OF A MULTIPLE C QUADRATURE. C FER IS THE RELATIVE ERROR COMMITTED COMPUTING F. C FNCVAL IS THE COMPUTED FUNCTION VALUE. C FSAVE IS USED TO STORE THE FUNCTION VALUE AT THE ABSCISSA WHERE C THE INTERVAL IS SUBDIVIDED. C FSAVED IS SET IF FSAVE CONTAINS USEFUL DATA. C FSTORE IS USED TO DETERMINE WHICH FUNCTION VALUES ARE TO BE STORED, C AND WHERE THEY ARE TO BE STORED IN THE FUNCTION TABLE. C FSTORE IS READ-ONLY. THE CONTENT OF EACH ELEMENT OF FSTORE C IS A TWO DIGIT BASE-100 NUMBER, WHERE THE FIRST DIGIT C INDICATES WHERE THE FIRST FUNCTION VALUE COMPUTED DURING THE C CURRENT INTEGRATION STEP IS TO BE STORED, AND THE SECOND DIGIT C INDICATES THE LAST CELL OF THE FUNCTION TABLE TO BE USED FOR C THE CURRENT INTEGRATION STEP. INTEGER FSTORE(7) C FT USED TO HOLD, TEMPORARILY, THE FUNCTION VALUES USED. C THESE FUNCTION VALUES ARE USED TO COMPUTE THE DIFFERENCE C LINE. SEE XT. C FUDGE IS USED DURING THE CONVERGENCE TEST TO RELATE THE ESTIMATED C ERROR TO THE ROUND-OFF LEVEL. FUDGE(K-KAIMT) GIVES THE C RELATIVE PRECISION THAT MUST BE OBTAINED IF THE ERROR C ESTIMATE IS NOT TO BE INCREASED. REAL FUDGE(5) C FUNCT IS USED TO STORE FUNCTION VALUES WHICH MAY BE RE-USED BY C THE PATTERSON FORMULAE, AND ADDITIONAL FUNCTION VALUES TO C COMPUTE THE DIFFERENCE LINE. C FUNCT(1-17) IS STORAGE NEEDED FOR THE PATTERSON METHOD C FUNCT(18-24) IS STORAGE TO RECOVER F(X1) THROUGH F(X15). C F1, F2 CONTAIN NEWLY COMPUTED FUNCTION VALUES DURING THE QUADRATURE C STEP. THEY ARE ALSO USED ELSEWHERE FOR TEMPORARY STORAGE. C GAMMA IS USED DURING COMPUTATION OF THE NEXT STEPSIZE. REAL GAMMA(7) C HAVDIF IS SET TO INDICATE THAT THE DIFFERENCE LINE HAS BEEN COMPUTED C USING 15 FUNCTION VALUES. C I IS USED FREELY AS AN INDEX. C IC1 IC2 ARE THE INDICES OF CONVERGENCE IN THE FORWARD AND BACKWARD C DIFFERENCE LINES. IC1 AND IC2 ARE EQUIVALENT TO ISTOP(1,1) C AND ISTOP(1,2). (SEE BELOW). C IEND IS SET WHEN THE END OF A MAJOR SUBDIVISION IS IN THE CURRENT C PANEL. C IH ARE INDICES USED TO LOCATE ABSCISSAE IN THE TABLE OF C QUADRATURE RULES. IH IS USED DURING COMPUTATION OF THE C DIFFERENCE LINE. INTEGER IH(7) C INC INC2 ARE USED DURING ANALYSIS OF THE DIFFERENCE LINE AS INDEX C INCREMENTORS. THEY ARE ALSO USED FREELY AS INDICES. C INEW IS USED DURING THE QUADRATURE STEP TO CONTROL A LOOP. INEW C MAY NOT BE CHANGED BETWEEN QUADRATURE STEPS. C INIT IS SET TO INDICATE THAT THE ORIGINAL INTERVAL IS BEING C INTEGRATED. C INSTOP IS USED TO CONTROL A LOOP DURING ANALYSIS OF THE DIFFERENCE C LINE. C IOLD IS USED DURING THE QUADRATURE STEP TO CONTROL A LOOP. IOLD C MAY NOT BE CHANGED BETWEEN QUADRATURE STEPS. C IP IS USED DURING THE QUADRATURE STEP TO INDEX THE TABLE OF C ABSCISSAE AND WEIGHTS. IP MAY NOT BE CHANGED BETWEEN C QUADRATURE STEPS. C IPRINT CONTROLS THE LEVEL OF DIAGNOSTIC PRINTING. C 0 - NO PRINTING C 1 - MINIMUM PRINTING - ERROR MESSAGES (DEFAULT) C 2 - PANEL BOUNDARIES AND ANSWERS C 3 - ERROR ESTIMATES FOR EACH QUADRATURE FORMULA C 4 - DETAILED OUTPUT (DIFFERENCE LINES, ETC). C ISTOP HOLDS RESULTS OF ANALYSIS OF DIFFERENCE LINE. C ENTRIES IN ISTOP MEAN C ISTOP(1,1) - INDEX OF XT TO WHICH CONVERGENCE EXTENDS FROM C THE ALOCAL END OF THE INTERVAL, C ISTOP(1,2) - INDEX OF XT TO WHICH CONVERGENCE EXTENDS FROM C THE BLOCAL END OF THE INTERVAL C ISTOP(2,1) AND ISTOP(2,2) MAY INDICATE INDICES OF XT WHERE C THE INTERVAL SHOULD BE SUBDIVIDED. C IXKDIM INDEXES THE CELL IN IFLAG IN WHICH THE CURRENT DIMENSION OF A C MULTIDIMENSIONAL INTEGRATION IS STORED. THE DEFAULT IS 1, C CORRESPONDING WITH THE USUAL CELL USED FOR COMMUNICATION WITH C SINTF. THE USER MAY NEED TO CHANGE THIS BY USING OPTION 12 C IF NON-STANDARD DIMENSION CHANGES ARE NECESSARY, AND A CHANGE C OF VARIABLE IS SIMULTANEOUSLY NECESSARY. SEE SINTM AND C SINTMA. C J IS FREELY USED AS AN INDEX. C JUMPS STORES THE QUALITATIVE RESULT OF THE DIFFERENCE LINE ANALYSIS. C 1 - NO JUMPS C 2 - ONE JUMP, ON THE END OF THE LINE C 3 - ONE JUMP, NOT ON THE END C 4 - TWO JUMPS. C J1 J2 ARE USED DURING THE QUADRATURE STEP TO CONTROL STORING INTO C THE FUNCT ARRAY. THEY ARE RECOMPUTED AT EACH STEP. THEY C ARE ALSO USED DURING THE ANALYSIS OF THE DIFFERENCE LINE TO C STORE 'JUMP' INDICES. C J1OLD ARE USED TO REMEMBER J1 AND J2 DURING A SEARCH C J2OLD FOR A JUMP. C K IS THE INDEX FOR THE QUADRATURE FORMULA LAST USED. THE C NUMBER OF POINTS USED IN THE KTH FORMULA IS 2**K - 1. C KAIMT IS THE INDEX OF THE QUADRATURE FORMULA EXPECTED TO CONVERGE C ON THE NEXT STEP. C KDIM IS THE CURRENT DIMENSION - ALWAYS 1 IN 1 DIMENSIONAL CASE. C KK IS A TEMPORARY INTEGER STORAGE. C KMAX IS THE MAXIMUM VALUE ALLOWED FOR K ON THE CURRENT INTEGRATION C STEP. C KMAXF IS THE VALUE TO USE FOR KMAX INITIALLY AND AFTER A FAILURE. C KMIN IS THE MINIMUM QUADRATURE FORMULA INDEX TO ACCEPT IF THERE C ARE NO JUMPS. KMIN IS USUALLY ZERO, BUT IS SET TO K IF THE C STEP SIZE INCREASES. C KONVRG KONVRG IS THE TOTAL NUMBER OF NODES OF CONVERGENCE OBSERVED C IN THE DIFFERENCE LINES. KONVRG=IC1+LENDT-IC2. C KORECT IS USED TO DETERMINE WHICH OF THE STORED FUNCTION VALUES ARE C TO BE USED TO CORRECT ONE-HALF OF THE PREVIOUSLY COMPUTED C INTEGRAL ON THE CURRENT INTERVAL. THE INDICES ARE STORED IN C KORECT AS BASE-20 DIGITS. KORECT IS READ-ONLY, INDEXED BY K. INTEGER KORECT(6) C L IS USED DURING ANALYSIS OF THE DIFFERENCE LINE. C LENDT IS THE CURRENT LENGTH OF THE DIFFERENCE LINE. C LOCAL IS A VECTOR EQUIVALENT TO ALOCAL,BLOCAL. C LOCAL(3),LOCAL(4) ARE ALOCAL,BLOCAL IN USER COORDINATES. C MAXK IS THE MAXIMUM VALUE ALLOWED FOR K. C NEEDH IS SET IF THE HEADER FOR DIAGNOSTIC INFORMATION IS NEEDED. C NEWEPS MEANS EPS MUST BE RECOMPUTED TO ACCOUNT FOR ERROR IN LIMITS LOGICAL NEWEPS C NFEVAL IS THE CURRENT NUMBER OF FUNCTION EVALUATIONS. C NFJUMP IS THE VALUE OF NFEVAL WHEN A JUMP IS FIRST NOTICED. C NFMAX IS THE INDEX OF THE ENTRY IN THE OPTION VECTOR (IOPT) TO C USE TO CONTROL THE MAXIMUM NUMBER OF FUNCTION EVALUATIONS. C NOMOUT IS THE NOMINAL OUTPUT UNIT. C NSUB IS THE POWER OF THE TRANSFORMATION OF THE FORM C X=TA+(T-TA)*((T-TA)/(TB-TA))**(NSUB-1) C WHERE X IS THE USER'S INDEPENDENT VARIABLE AND T IS THE C INDEPENDENT VARIABLE TO USE FOR INTEGRATION. NSUB MAY C ONLY BE A POWER OF 2. NSUB=0 REALLY MEANS NSUB=1, BUT C IS SIMPLER TO USE. C NSUBMX IS THE MAXIMUM VALUE FOR NSUB. C NSUBSV IS THE VALUE TO USE FOR NSUB WHEN RESTARTING C AFTER A SUBDIVISION. C NXKDIM IS THE VALUE OF KDIM BEFORE AN INNER INTEGRAL WAS REQUESTED. C THE INNER INTEGRAL IS NOT NECESSARILY OF DIMENSION KDIM-1. C SEE IXKDIM ABOVE, SINTM AND SINTMA. C P IS THE ARRAY OF ABSCISSAE AND WEIGHTS. DATA FOR P IS INCLUDED C AFTER THE DESCRIPTION OF INTERNAL VARIABLES. REAL P(305) C PAACUM IS THE PREVIOUS VALUE OF AACUM. C PACUM IS THE PREVIOUS VALUE OF ACUM. ABS(ACUM-PACUM) IS USED AS AN C INITIAL ERROR ESTIMATE. C PART IS THE INDEX OF THE CURRENT MAJOR SUBDIVISION. START, END, C STEP, ERRT AND RESULT ARE SUBSCRIPTED BY PART. PART BEGINS C SET TO 1, IS CHANGED TO 2 WHEN THE INTERVAL IS SUBDIVIDED, C AND SET BACK TO 1 WHEN PART 2 IS INTEGRATED. WHEN PART 1 C IS FINISHED, THE ENTIRE PROBLEM IS FINISHED. C PEPSMN IS THE PREVIOUS VALUE OF EPSMIN. PEPSMN MAY NOT BE CHANGED C BETWEEN QUADRATURE STEPS. C PERR IS THE ERROR COMMITTED ON THE LAST SUCCESSFULLY INTEGRATED C PANEL. PERR IS USED DURING THE NOISE TEST. C PF1 PF2 ARE THE PREVIOUS VALUES OF F1 AND F2. PF1 AND PF2 ARE USED TO C COMPUTE EPSMIN. C PHI CONTAINS THE DIFFERENCE LINES. PHI IS DIMENSIONED AND C EQUIVALENCED WITH PHIT AFTER THE DESCRIPTION OF INTERNAL C VARIABLES. C PHIS IS A TWO-DIMENSIONAL ARRAY EQUIVALENCED TO PHI AND PHIT. IT C IS USED TO ACCESS THE DIFFERENCE LINES USING AN INDEX. REAL PHIS(17,2) C PHISUM IS THE SUM OF 3 MIDDLE ENTRIES OF PHI. C PHIT IS THE BACKWARD DIFFERENCE LINE. SEE PHI. REAL PHIT(17) C PHTSUM IS THE SUM OF 3 MIDDLE ENTRIES OF PHIT. C PX IS THE PREVIOUS VALUE OF THE OFFSET FROM THE ENDS OF THE C CURRENT INTERVAL. PX IS USED TO COMPUTE EPSMIN. C RE IS THE RATIO OF THE ERROR COMMITTED WITH THE CURRENT C QUADRATURE FORMULA TO THAT COMMITTED WITH THE PREVIOUS C QUADRATURE FORMULA. RE IS USED TO DETECT CONVERGENCE AND C TO AVOID FALSE CONVERGENCE, AND TO DETECT NOISE. C RECORR IS USED TO CORRECT RE WHEN NSUB IS NOT ZERO. ONLY THE FIRST C AND THIRD ELEMENTS OF RECORR ARE USED EXPLICITLY. REAL RECORR(3) C RELEPS IS THE ABSOLUTE ERROR TOLERANCE, BASED UPON THE RELATIVE C RELATIVE TOLERANCE REQUESTED AND THE BEST AVAILABLE ESTIMATE C OF THE ANSWER. C RELOBT IS THE PRECISION RELATIVE TO THE OBTAINABLE PRECISION. C THE ESTIMATED ERROR IS COMPARED TO C (EPSMIN/AACUM)**RELOBT*AACUM. C RELOOK IS USED TO FORCE EXAMINATION OF THE DIFFERENCE LINES, IF C THE METHOD IS TO ACCEPT AN ANSWER BUT RE*REP .GT. RELOOK. REAL RELOOK C RELTOL IS THE INDEX INTO WORK OF THE PARAMETERS FOR ABSOLUTE, LOCAL C RELATIVE, AND GLOBAL RELATIVE TOLERANCE PARAMETERS. C REPROD IS THE PRODUCT OF ALL VALUES OF RE FOR THE CURRENT PANEL. C REPROD IS NOT ALLOWED TO BE GREATER THAN 1.0. C REP IS THE PREVIOUS VALUE OF RE. IT IS USED SIMILARLY TO RE. C RESULT IS THE CUMULATIVE INTEGRAL FOR THE CURRENT MAJOR SUBDIVISION. C REVERS IS THE INDEX OF THE ENTRY IN THE OPTION VECTOR (IOPT) TO C USE TO CONTROL REVERSE COMMUNICATION. C RNDC IS THE NOISE OR ROUND-OFF LEVEL, INITIALLY EQUAL TO FER. C RNDC IS CHANGED IF NOISE IS DETECTED. C ROUNDF IS SET TO INDICATE THAT ROUND-OFF IS LIMITING THE ERROR, OR C THAT THE INTEGRAND IS BEING TREATED AS PURE NOISE WITH RESPECT C TO THE REQUESTED ERROR TOLERANCE. C S IS A SCRATCH VARIABLE. C SEARCH IS USED TO INDICATE WHETHER POINTS NOT USED BY THE INTEGRATION C FORMULAE ARE BEING ADDED TO THE DIFFERENCE LINES. C 1 - POINTS NOT BEING ADDED C 2 - SEARCHING FOR INITIAL INTERVAL SIZE C 3 - MAKING SURE AN END-POINT JUMP IS NOT INSIDE THE INTERVAL C 4 - SEARCHING FOR A DISCONTINUITY C 5 - CHANGING FROM SEARCH=2 OR 3 TO SEARCH=4 C 6 - NSUB WAS JUST REDUCED BECAUSE WE SAW JUMPS. C NSUB CANNOT BE INCREASED WITHOUT DOING A TYPE-2 SEARCH. C NEGATIVE - INITIALIZE THE SEARCH (TRANSIENT STATE). C SMIN IS THE MINIMUM STEPSIZE TO USE FOR A PANEL BEGINNING AT SX. REAL SMIN C SPACE RESERVES SOME SPACE IN /SINTC/ IN CASE IT IS NEEDED. THIS C IS DONE SO USERS' CODE WON'T NEED TO CHANGE. C START STORES THE 'LEFT' BOUNDARY OF A PART. C STEP STORES THE SIZE OF THE FIRST STEP TO USE ON A MAJOR C SUBDIVISION. THE SIGN OF STEP IS THE SIGN TO USE FOR C ACCUMULATION OF THE ANSWER INTO RESULT(PART). C SUM IS THE CENTER OF THE CURRENT INTERVAL. C T IS THE TRANSFORMED INDEPENDENT VARIABLE. SEE NSUB. C TA TB ARE THE INTERVAL ON WHICH A VARIABLE TRANSFORMATION C IS DEFINED. SEE NSUB. THE SIGN OF TB IS ALSO USED TO C INDICATE THE DIRECTION OF STEP ACCUMULATION. C TALOC IF NON-ZERO IS THE INDEX IN WORK OF THE LOCATION OF A C SINGULARITY. C TASAVE IS THE VALUE TO USE FOR TA WHEN RESTARTING PART 1 C AFTER SUBDIVIDING. C TDECR IS AN ARITHMETIC STATEMENT FUNCTION DEFINED BELOW. REAL TDECR C TEND IS THE POINT AT WHICH USAGE OF A T**N SUBSTITUTION SHOULD BE C DISCONTINUED. C THROW IS USED DURING THE SEARCH FOR A DISCONTINUITY TO INDICATE C WHICH END OF THE DIFFERENCE LINES SHOULD BE DISCARDED. INTEGER THROW C TINCR IS AN ARITHMETIC STATEMENT FUNCTION DEFINED BELOW. REAL TINCR C TLEN IS THE TOTAL LENGTH OF THE UN-INTEGRATED PART. TLEN IS USED C TO PARTION ERROR ALLOTTMENTS. C TP TP1 ARE TEMPORARY VARIABLES. REAL TP1 C TPS IS A TEMPORARY VARIABLE. C WHERE IS USED AS A COMPUTED GO TO INDEX. C WHERE2 IS USED AS A COMPUTED GO TO INDEX. C WORRY IS USED TO STORE THE CAUTION POINTS FOR EACH PART. C X DURING THE QUADRATURE STEP, X IS THE DISTANCE FROM THE ENDS OF C THE INTERVAL, AT WHICH ABSCISSAE THE INTEGRAND IS EVALUATED. C X IS ALSO USED ELSEWHERE FOR TEMPORARY STORAGE. C XEPS IS THE MINIMUM PART OF THE REMAINING ERROR TOLERANCE TO C REQUIRE ON THE CURRENT INTERVAL. REAL XEPS C XJ IS THE MAGNITUDE OF THE LARGEST JUMP IN THE BACKWARD C DIFFERENCE LINE. C XJP IS THE MAGNITUDE OF THE LARGEST JUMP IN THE FORWARD C DIFFERENCE LINE. C XJUMP IS A THRESHOLD FOR JUMPS IN THE DIFFERENCE LINE. C XJUMPS IS THE VALUE OF XJUMP TO USE DURING A SEARCH. REAL XJUMPS C XSTEP IS USED TO INCREASE THE CURRENT EFFECTIVE STEPSIZE AFTER C SUCCESSFULLY INTEGRATING AN INTERVAL. THE CURRENT EFFECTIVE C STEPSIZE IS USED TO COMPUTE THE NEXT STEPSIZE. XSTEP IS C READ-ONLY, INDEXED BY MAX(AIMFOR-K+3,1). REAL XSTEP(7) C XT THE ABSCISSAE OF THE QUADRATURE FORMULA. XT ARE USED TO C COMPUTE THE DIFFERENCE LINE. C X1 X2 ARE TEMPORARY STORAGE. C Z ARE ADJUSTABLE PARAMETERS USED FOR DEVELOPMENT. REAL Z(4) C ZL1 IS THE ARGUMENT OF ARITHMETIC STATEMENT FUNCTIONS. REAL ZL1 C C ***** COMMON STORAGE ****************************************** C C COMMON /SINTNC/ CONTAINS VARIABLES NOT SEPARATELY SAVED FOR C EACH DIMENSION OF A MULTIPLE QUADRATURE. COMMON /SINTC/ C CONTAINS VARIABLES THAT MUST BE SAVED FOR EACH DIMENSION OF THE C QUADRATURE. THE VARIABLES IN EACH COMMON BLOCK ARE STORED IN THE C ORDER - ALWAYS DOUBLE, DOUBLE IF DOUBLE PRECISION PROGRAM, DOUBLE C IF DOUBLE PRECISION PROGRAM AND EXPONENT RANGE OF DOUBLE AND C SINGLE VERY DIFFERENT, SINGLE, INTEGER, LOGICAL. A PAD OF LOGICAL C VARIABLES IS INCLUDED AT THE END OF /SINTC/. THE DIMENSION OF C THE PAD MAY NEED TO BE VARIED SO THAT NO VARIABLES BEYOND THE END C OF THE COMMON BLOCK ARE ALTERED. C C DECLARATIONS OF COMMON /SINTNC/ VARIABLES. C REAL AINIT, BINIT, FNCVAL, S, TP REAL FER, FER1, RELOBT, TPS, XJ, XJP INTEGER FEA, FEA1, INC, INC2, IPRINT, 1 ISTOP(2,2),JPRINT, KDIM, KK, KMAXF, NDIM, 2 NFINDX, NFMAX, NFMAXM, RELTOL, REVERM, REVERS, 3 WHEREM LOGICAL NEEDH C C DECLARATIONS OF COMMON /SINTC/ VARIABLES. C c--D Next line special: S => D, X => Q, D => D, P => D DOUBLE PRECISION ACUM, PACUM, RESULT(2) C 139 $.TYPE.$ VARIABLES REAL 1 AACUM, ABSCIS, DELMIN, DELTA, DIFF, DISCX(2), 2 END(2), ERRINA, ERRINB, FAT(2), FSAVE, 3 FUNCT(24), F1, F2, LOCAL(4), PAACUM, PF1, 4 PF2, PHISUM, PHTSUM, PX, SPACE(6), 5 STEP(2), START(2), SUM, T, TA, TASAVE, 6 TB, TEND, WORRY(2), X, X1, 7 X2, XT(17), FT(17), PHI(34) c Note XT, FT, and PHI above are last, because they must be in adjacent c locations in SINTC. C 30 $DSTYP$ VARIABLES REAL 1 ABSDIF, COUNT, EDUE2A, EDUE2B, EP, EPNOIZ, 2 EPS, EPSMAX, EPSMIN, EPSO, EPSR, EPSS, 3 ERR, ERRAT(2), ERRC, ERRF, ERRI, ERRT(2), 4 ESOLD, EXTRA, PEPSMN, RE, RELEPS, REP, 5 REPROD, RNDC, TLEN, XJUMP C 29 INTEGER VARIABLES INTEGER DISCF, DISCHK, ENDPTS, I, INEW, 1 IOLD, IP, IXKDIM, J, J1, J1OLD, 2 J2, J2OLD, K, KAIMT, KMAX, KMIN, 3 L, LENDT, NFEVAL, NFJUMP, NSUB, NSUBSV, 4 NXKDIM, PART, SEARCH, TALOC, WHERE, WHERE2 C 11 TO 18 LOGICALS (7 ARE PADDING). LOGICAL DID1, FAIL, FATS(2), FSAVED, HAVDIF, 1 IEND, INIT, ROUNDF, XCDOBT(2), PAD(7) C C THE COMMON BLOCKS. C COMMON /SINTNC/ c 1 2 3 4 5 6 7 8 W AINIT, BINIT, FNCVAL, S, TP, FER, FER1, RELOBT, c 9 10 11 12 13 1 2 3 X TPS, XJ, XJP, FEA, FEA1, KDIM, INC, INC2, c 4 (2,2) 8 9 10 11 12 13 14 Y ISTOP, JPRINT, IPRINT, KK, KMAXF, NDIM, NFINDX, NFMAX, c 15 16 17 18 19 20 Z NFMAXM, RELTOL, REVERM, REVERS, WHEREM, NEEDH COMMON /SINTC/ 1 ACUM, PACUM, RESULT COMMON /SINTC/ c 1 2 (4) 6 7 8 9 10 11 (2) 1 AACUM, LOCAL, ABSCIS, TA, DELTA, DELMIN, DIFF, DISCX, c 13 (2) 15 16 17 (2) 19 20 (24) 44 2 END, ERRINA, ERRINB, FAT, FSAVE, FUNCT, F2, c 45 46 47 48 49 50 51 (6) 3 PAACUM, PF1, PF2, PHISUM, PHTSUM, PX, SPACE, c 57 (2) 59 (2) 61 62 63 64 65 4 STEP, START, SUM, T, TASAVE, TB, TEND, c 66 (2) 68 69 70 71 72 5 WORRY, X1, X2, X, F1, COUNT, c 73 (17) 90 (17) 107 (34) 6 XT, FT, PHI COMMON /SINTC/ c 141 142 143 144 145 146 1 ABSDIF, EDUE2A, EDUE2B, EP, EPNOIZ, EPSMAX, c 147 148 149 150 (2) 152 153 2 EPSO, EPSR, EPSS, ERRAT, ERRC, ERRF, c 154 (2) 156 157 158 159 160 3 ERRT, ESOLD, EXTRA, PEPSMN, RELEPS, REP, c 161 162 163 4 RNDC, TLEN, XJUMP, c 164 165 166 167 168 169 5 ERRI, ERR, EPSMIN, EPS, RE, REPROD COMMON /SINTC/ c 170 171 172 1 DISCF, DISCHK, ENDPTS, INEW, IOLD, IP, IXKDIM, 2 J, J1, J1OLD, J2, J2OLD, KMAX, KMIN, 3 L, LENDT, NFJUMP, NSUBSV, NXKDIM, TALOC, WHERE2, c 1 2 3 4 5 6 7 8 4 I, K, KAIMT, NSUB, PART, SEARCH, WHERE, NFEVAL COMMON /SINTC/ 1 DID1, FAIL, FATS, FSAVED, HAVDIF, IEND, INIT, ROUNDF, 2 XCDOBT, PAD SAVE /SINTNC/, /SINTC/ C C THE VARIABLES HERE DEFINE THE MACHINE ENVIRONMENT. ALL ARE SET C IN DINTOP. THE MEANING ATTACHED TO THESE VARIABLES CAN BE C FOUND BY LOOKING AT THE DEFINITIONS IN DINTOP. REAL 1 EMEPS, EEPSM8, EDELM2, EDELM3, ESQEPS, ERSQEP, ERSQE6, EMINF, 2 ESMALL, ENZER, EDELM1, ENINF COMMON /SINTEC/ 1 EMEPS, EEPSM8, EDELM2, EDELM3, ESQEPS, ERSQEP, ERSQE6, EMINF, 2 ESMALL, ENZER, EDELM1, ENINF SAVE /SINTEC/ C C ***** EQUIVALENCE STATEMENTS ******************************* C EQUIVALENCE (PHI(18),PHIT),(PHI(1),PHIS(1,1)) EQUIVALENCE (FAT(1),FATA), (FAT(2),FATB) EQUIVALENCE (FATS(1),FATAS), (FATS(2),FATBS) EQUIVALENCE (ISTOP(1,1),IC1), (ISTOP(1,2),IC2) EQUIVALENCE (LOCAL(1),ALOCAL), (LOCAL(2),BLOCAL) C C ***** DATA STATEMENTS *************************************** C DATA AIMFOR /6/ DATA ENDPHI(1) /1/ DATA ERRCF /1.579E0,5.233E-2,8.501E-3/ DATA FSTORE /202,304,508,916,1217,1417,100/ DATA FUDGE /1.0E-8,1.0E-6,1.0E-5,1.0E-4,3.0E-4/ C K-KAIMT = -4 -3 -2 -1 0 DATA GAMMA /1.42e0,2.28e0,4.62e0,19.00e0,75.00e0,75.00e0,75.00e0/ DATA IH /14,6,16,2,18,8,20/ DATA KORECT /41,81,161,321,55890063,56395705/ DATA MAXK /8/ DATA NSUBMX /4/ DATA RECORR /1.0E-3, 0.0E0, 1.0E-1/ DATA RELOOK /5.0E-3/ DATA XEPS / 0.02E0 / DATA XJUMPS /0.5E0/ DATA XSTEP /0.5e0,0.875e0,1.0e0,2.25e0,5.0e0,15.0e0,50.0e0/ DATA Z /1.0e0, 1.0e0, 1.75e0, 6.0e0/ C MULTIPLY DELTA BY Z(1) IN INITIAL INTERVAL SEARCH C DELTA.LE.Z(2)*ABSDIF IN INITIAL INTERVAL SEARCH MEANS DONT GO ON C DELTA.LT.Z(3)*ABSDIF IF K.LT.6 AS FOR Z(2) C DELTA.GT.Z(4)*ABSDIF MEANS KMAX=5 AFTER INITIAL INTERVAL SEARCH C C IN THE COMMENTS BELOW, F(K,I) REFERS TO THE FUNCTION VALUE C COMPUTED FOR THE I'TH NODE OF THE K'TH FORMULA. C P Points C Indices Use (C = Correction Coefficient N&W = Nodes and Weights) C 1 3 C for F(1,1) C 2- 3 N&W for F(2,1) C 4- 5 7 C for F(1,1) and F(2,1) C 6- 9 N&W for F(3, 1-2) C 10- 13 15 C for F(1,1), F(2,1), F(3,1-2) C 14- 21 N&W for F(4, 1-4) C 22- 29 31 C for F(1,1), F(2,1), F(3,1-2), F(4,1-4). C 30- 45 N&W FOR F(5,1-8). C 46- 61 63 C for F(1,1), F(2,1), F(3,1-2), F(4,1-4), F(5,1-8). C 62- 93 N&W FOR F(6,1-16). C 94-105 127 C for F(3,1), F(4,1-2), F(5,1-3), F(6,1-6). C 106-169 N&W FOR F(7,1-32). C 170-177 256 C for F(4,1), F(5,1), F(6,1-2), F(7,1-4). C 178:305 N&W FOR F(8,1-64). C C CORRECTIONS, NODES, AND WEIGHTS FOR THE 3-POINT FORMULA. C c++ Save data by elements if ~.C. DATA P(1) / -.11111111111111111111E+00 / DATA P(2) / +.22540333075851662296E+00 / DATA P(3) / +.55555555555555555556E+00 / DATA P(4) / +.647209421402969791E-02 / DATA P(5) / -.928968790944433705E-02 / DATA P(6) / +.39508731291979716579E-01 / DATA P(7) / +.10465622602646726519E+00 / DATA P(8) / +.56575625065319744200E+00 / DATA P(9) / +.40139741477596222291E+00 / DATA P(10) / +.5223046896961622E-04 / DATA P(11) / +.17121030961750000E-03 / DATA P(12) / -.724830016153892898E-03 / DATA P(13) / -.7017801099209042E-04 / DATA P(14) / +.61680367872449777899E-02 / DATA P(15) / +.17001719629940260339E-01 / DATA P(16) / +.11154076712774300110E+00 / DATA P(17) / +.92927195315124537686E-01 / DATA P(18) / +.37889705326277359705E+00 / DATA P(19) / +.17151190913639138079E+00 / DATA P(20) / +.77661331357103311837E+00 / DATA P(21) / +.21915685840158749640E+00 / DATA P(22) / +.682166534792E-08 / DATA P(23) / +.12667409859336E-06 / DATA P(24) / +.59565976367837165E-05 / DATA P(25) / +.1392330106826E-07 / DATA P(26) / -.6629407564902392E-04 / DATA P(27) / -.704395804282302E-06 / DATA P(28) / -.34518205339241E-07 / DATA P(29) / -.814486910996E-08 / DATA P(30) / +.90187503233240234038E-03 / DATA P(31) / +.25447807915618744154E-02 / DATA P(32) / +.18468850446259893130E-01 / DATA P(33) / +.16446049854387810934E-01 / DATA P(34) / +.70345142570259943330E-01 / DATA P(35) / +.35957103307129322097E-01 / DATA P(36) / +.16327406183113126449E+00 / DATA P(37) / +.56979509494123357412E-01 / DATA P(38) / +.29750379350847292139E+00 / DATA P(39) / +.76879620499003531043E-01 / DATA P(40) / +.46868025635562437602E+00 / DATA P(41) / +.93627109981264473617E-01 / DATA P(42) / +.66886460674202316691E+00 / DATA P(43) / +.10566989358023480974E+00 / DATA P(44) / +.88751105686681337425E+00 / DATA P(45) / +.11195687302095345688E+00 / DATA P(46) / +.371583E-15 / DATA P(47) / +.21237877E-12 / DATA P(48) / +.10522629388435E-08 / DATA P(49) / +.1748029E-14 / DATA P(50) / +.3475718983017160E-06 / DATA P(51) / +.90312761725E-11 / DATA P(52) / +.12558916E-13 / DATA P(53) / +.54591E-15 / DATA P(54) / -.72338395508691963E-05 / DATA P(55) / -.169699579757977E-07 / DATA P(56) / -.854363907155E-10 / DATA P(57) / -.12281300930E-11 / DATA P(58) / -.462334825E-13 / DATA P(59) / -.42244055E-14 / DATA P(60) / -.88501E-15 / DATA P(61) / -.40904E-15 / DATA P(62) / +.12711187964238806027E-03 / DATA P(63) / +.36322148184553065969E-03 / DATA P(64) / +.27937406277780409196E-02 / DATA P(65) / +.25790497946856882724E-02 / DATA P(66) / +.11315242452570520059E-01 / DATA P(67) / +.61155068221172463397E-02 / DATA P(68) / +.27817125251418203419E-01 / DATA P(69) / +.10498246909621321898E-01 / DATA P(70) / +.53657141626597094849E-01 / DATA P(71) / +.15406750466559497802E-01 / DATA P(72) / +.89628843042995707499E-01 / DATA P(73) / +.20594233915912711149E-01 / DATA P(74) / +.13609206180630952284E+00 / DATA P(75) / +.25869679327214746911E-01 / DATA P(76) / +.19305946804978238813E+00 / DATA P(77) / +.31073551111687964880E-01 / DATA P(78) / +.26024395564730524132E+00 / DATA P(79) / +.36064432780782572640E-01 / DATA P(80) / +.33709033997521940454E+00 / DATA P(81) / +.40715510116944318934E-01 / DATA P(82) / +.42280428994795418516E+00 / DATA P(83) / +.44914531653632197414E-01 / DATA P(84) / +.51638197305415897244E+00 / DATA P(85) / +.48564330406673198716E-01 / DATA P(86) / +.61664067580126965307E+00 / DATA P(87) / +.51583253952048458777E-01 / DATA P(88) / +.72225017797817568492E+00 / DATA P(89) / +.53905499335266063927E-01 / DATA P(90) / +.83176474844779253501E+00 / DATA P(91) / +.55481404356559363988E-01 / DATA P(92) / +.94365568695340721002E+00 / DATA P(93) / +.56277699831254301273E-01 / DATA P(94) / +.1041098E-15 / DATA P(95) / +.249472054598E-10 / DATA P(96) / +.55E-20 / DATA P(97) / +.290412475995385E-07 / DATA P(98) / +.367282126E-13 / DATA P(99) / +.5568E-18 / DATA P(100) / -.871176477376972025E-06 / DATA P(101) / -.8147324267441E-09 / DATA P(102) / -.8830920337E-12 / DATA P(103) / -.18018239E-14 / DATA P(104) / -.70528E-17 / DATA P(105) / -.506E-19 / DATA P(106) / +.17569645108401419961E-04 / DATA P(107) / +.50536095207862517625E-04 / DATA P(108) / +.40120032808931675009E-03 / DATA P(109) / +.37774664632698466027E-03 / DATA P(110) / +.16833646815926074696E-02 / DATA P(111) / +.93836984854238150079E-03 / DATA P(112) / +.42758953015928114900E-02 / DATA P(113) / +.16811428654214699063E-02 / DATA P(114) / +.85042788218938676006E-02 / DATA P(115) / +.25687649437940203731E-02 / DATA P(116) / +.14628500401479628890E-01 / DATA P(117) / +.35728927835172996494E-02 / DATA P(118) / +.22858485360294285840E-01 / DATA P(119) / +.46710503721143217474E-02 / DATA P(120) / +.33362148441583432910E-01 / DATA P(121) / +.58434498758356395076E-02 / DATA P(122) / +.46269993574238863589E-01 / DATA P(123) / +.70724899954335554680E-02 / DATA P(124) / +.61679602220407116350E-01 / DATA P(125) / +.83428387539681577056E-02 / DATA P(126) / +.79659974529987579270E-01 / DATA P(127) / +.96411777297025366953E-02 / DATA P(128) / +.10025510022305996335E+00 / DATA P(129) / +.10955733387837901648E-01 / DATA P(130) / +.12348658551529473026E+00 / DATA P(131) / +.12275830560082770087E-01 / DATA P(132) / +.14935550523164972024E+00 / DATA P(133) / +.13591571009765546790E-01 / DATA P(134) / +.17784374563501959262E+00 / DATA P(135) / +.14893641664815182035E-01 / DATA P(136) / +.20891506620015163857E+00 / DATA P(137) / +.16173218729577719942E-01 / DATA P(138) / +.24251603361948636206E+00 / DATA P(139) / +.17421930159464173747E-01 / DATA P(140) / +.27857691462990108452E+00 / DATA P(141) / +.18631848256138790186E-01 / DATA P(142) / +.31701256890892077191E+00 / DATA P(143) / +.19795495048097499488E-01 / DATA P(144) / +.35772335749024048622E+00 / DATA P(145) / +.20905851445812023852E-01 / DATA P(146) / +.40059606975775710702E+00 / DATA P(147) / +.21956366305317824939E-01 / DATA P(148) / +.44550486736806745112E+00 / DATA P(149) / +.22940964229387748761E-01 / DATA P(150) / +.49231224246628339785E+00 / DATA P(151) / +.23854052106038540080E-01 / DATA P(152) / +.54086998801016766712E+00 / DATA P(153) / +.24690524744487676909E-01 / DATA P(154) / +.59102017877011132759E+00 / DATA P(155) / +.25445769965464765813E-01 / DATA P(156) / +.64259616216846784762E+00 / DATA P(157) / +.26115673376706097680E-01 / DATA P(158) / +.69542355844328595666E+00 / DATA P(159) / +.26696622927450359906E-01 / DATA P(160) / +.74932126969651682339E+00 / DATA P(161) / +.27185513229624791819E-01 / DATA P(162) / +.80410249728889984607E+00 / DATA P(163) / +.27579749566481873035E-01 / DATA P(164) / +.85957576684743982540E+00 / DATA P(165) / +.27877251476613701609E-01 / DATA P(166) / +.91554595991628911629E+00 / DATA P(167) / +.28076455793817246607E-01 / DATA P(168) / +.97181535105025430566E+00 / DATA P(169) / +.28176319033016602131E-01 / DATA P(170) / +.3326E-18 / DATA P(171) / +.114094770478E-11 / DATA P(172) / +.2952436056970351E-08 / DATA P(173) / +.51608328E-15 / DATA P(174) / -.110177219650597323E-06 / DATA P(175) / -.58656987416475E-10 / DATA P(176) / -.23340340645E-13 / DATA P(177) / -.1248950E-16 / DATA P(178) / +.24036202515353807630E-05 / DATA P(179) / +.69379364324108267170E-05 / DATA P(180) / +.56003792945624240417E-04 / DATA P(181) / +.53275293669780613125E-04 / DATA P(182) / +.23950907556795267013E-03 / DATA P(183) / +.13575491094922871973E-03 / DATA P(184) / +.61966197497641806982E-03 / DATA P(185) / +.24921240048299729402E-03 / DATA P(186) / +.12543855319048853002E-02 / DATA P(187) / +.38974528447328229322E-03 / DATA P(188) / +.21946455040427254399E-02 / DATA P(189) / +.55429531493037471492E-03 / DATA P(190) / +.34858540851097261500E-02 / DATA P(191) / +.74028280424450333046E-03 / DATA P(192) / +.51684971993789994803E-02 / DATA P(193) / +.94536151685852538246E-03 / DATA P(194) / +.72786557172113846706E-02 / DATA P(195) / +.11674841174299594077E-02 / DATA P(196) / +.98486295992298408193E-02 / DATA P(197) / +.14049079956551446427E-02 / DATA P(198) / +.12907472045965932809E-01 / DATA P(199) / +.16561127281544526052E-02 / DATA P(200) / +.16481342421367271240E-01 / DATA P(201) / +.19197129710138724125E-02 / DATA P(202) / +.20593718329137316189E-01 / DATA P(203) / +.21944069253638388388E-02 / DATA P(204) / +.25265540247597332240E-01 / DATA P(205) / +.24789582266575679307E-02 / DATA P(206) / +.30515340497540768229E-01 / DATA P(207) / +.27721957645934509940E-02 / DATA P(208) / +.36359378430187867480E-01 / DATA P(209) / +.30730184347025783234E-02 / DATA P(210) / +.42811783890139037259E-01 / DATA P(211) / +.33803979910869203823E-02 / DATA P(212) / +.49884702478705123440E-01 / DATA P(213) / +.36933779170256508183E-02 / DATA P(214) / +.57588434808916940190E-01 / DATA P(215) / +.40110687240750233989E-02 / DATA P(216) / +.65931563842274211999E-01 / DATA P(217) / +.43326409680929828545E-02 / DATA P(218) / +.74921067092924347640E-01 / DATA P(219) / +.46573172997568547773E-02 / DATA P(220) / +.84562412844234959360E-01 / DATA P(221) / +.49843645647655386012E-02 / DATA P(222) / +.94859641186738404810E-01 / DATA P(223) / +.53130866051870565663E-02 / DATA P(224) / +.10581543166444097714E+00 / DATA P(225) / +.56428181013844441585E-02 / DATA P(226) / +.11743115975265809315E+00 / DATA P(227) / +.59729195655081658049E-02 / DATA P(228) / +.12970694445188609414E+00 / DATA P(229) / +.63027734490857587172E-02 / DATA P(230) / +.14264168911376784347E+00 / DATA P(231) / +.66317812429018878941E-02 / DATA P(232) / +.15623311732729139895E+00 / DATA P(233) / +.69593614093904229394E-02 / DATA P(234) / +.17047780536259859981E+00 / DATA P(235) / +.72849479805538070639E-02 / DATA P(236) / +.18537121234486258656E+00 / DATA P(237) / +.76079896657190565832E-02 / DATA P(238) / +.20090770903915859819E+00 / DATA P(239) / +.79279493342948491103E-02 / DATA P(240) / +.21708060588171698360E+00 / DATA P(241) / +.82443037630328680306E-02 / DATA P(242) / +.23388218069623990928E+00 / DATA P(243) / +.85565435613076896192E-02 / DATA P(244) / +.25130370638306339718E+00 / DATA P(245) / +.88641732094824942641E-02 / DATA P(246) / +.26933547875781873867E+00 / DATA P(247) / +.91667111635607884067E-02 / DATA P(248) / +.28796684463774796540E+00 / DATA P(249) / +.94636899938300652943E-02 / DATA P(250) / +.30718623022088529711E+00 / DATA P(251) / +.97546565363174114611E-02 / DATA P(252) / +.32698116976958152079E+00 / DATA P(253) / +.10039172044056840798E-01 / DATA P(254) / +.34733833458998250389E+00 / DATA P(255) / +.10316812330947621682E-01 / DATA P(256) / +.36824356228880576959E+00 / DATA P(257) / +.10587167904885197931E-01 / DATA P(258) / +.38968188628481359983E+00 / DATA P(259) / +.10849844089337314099E-01 / DATA P(260) / +.41163756555233745857E+00 / DATA P(261) / +.11104461134006926537E-01 / DATA P(262) / +.43409411457634557737E+00 / DATA P(263) / +.11350654315980596602E-01 / DATA P(264) / +.45703433350168850951E+00 / DATA P(265) / +.11588074033043952568E-01 / DATA P(266) / +.48044033846254297801E+00 / DATA P(267) / +.11816385890830235763E-01 / DATA P(268) / +.50429359208123853983E+00 / DATA P(269) / +.12035270785279562630E-01 / DATA P(270) / +.52857493412834112307E+00 / DATA P(271) / +.12244424981611985899E-01 / DATA P(272) / +.55326461233797152625E+00 / DATA P(273) / +.12443560190714035263E-01 / DATA P(274) / +.57834231337383669993E+00 / DATA P(275) / +.12632403643542078765E-01 / DATA P(276) / +.60378719394238406082E+00 / DATA P(277) / +.12810698163877361967E-01 / DATA P(278) / +.62957791204992176986E+00 / DATA P(279) / +.12978202239537399286E-01 / DATA P(280) / +.65569265840056197721E+00 / DATA P(281) / +.13134690091960152836E-01 / DATA P(282) / +.68210918793152331682E+00 / DATA P(283) / +.13279951743930530650E-01 / DATA P(284) / +.70880485148175331803E+00 / DATA P(285) / +.13413793085110098513E-01 / DATA P(286) / +.73575662758907323806E+00 / DATA P(287) / +.13536035934956213614E-01 / DATA P(288) / +.76294115441017027278E+00 / DATA P(289) / +.13646518102571291428E-01 / DATA P(290) / +.79033476175681880523E+00 / DATA P(291) / +.13745093443001896632E-01 / DATA P(292) / +.81791350324074780175E+00 / DATA P(293) / +.13831631909506428676E-01 / DATA P(294) / +.84565318851862189130E+00 / DATA P(295) / +.13906019601325461264E-01 / DATA P(296) / +.87352941562769803314E+00 / DATA P(297) / +.13968158806516938516E-01 / DATA P(298) / +.90151760340188079791E+00 / DATA P(299) / +.14017968039456608810E-01 / DATA P(300) / +.92959302395714482093E+00 / DATA P(301) / +.14055382072649964277E-01 / DATA P(302) / +.95773083523463639678E+00 / DATA P(303) / +.14080351962553661325E-01 / DATA P(304) / +.98590611358921753738E+00 / DATA P(305) / +.14092845069160408355E-01 / C C ***** STATEMENT FUNCTIONS ********************************* C C TDECR IS USED TO TRANSFORM AN ABSCISSA FROM THE CURRENT COORDINATE C SYSTEM TO ONE IN WHICH NSUB IS DECREMENTED BY A FACTOR OF 2. C TDECR(ZL1)=TA+(ZL1-TA)*((ZL1-TA)/TB) TDECR(ZL1)=TA*(1.0e0+TA/TB)+ZL1*((ZL1-TA)/TB-TA/TB) C TINCR IS USED TO TRANSFORM AN ABSCISSA FROM THE CURRENT COORDINATE C SYSTEM TO ONE IN WHICH NSUB IS INCREMENTED BY A FACTOR OF 2. TINCR(ZL1)=TA+SIGN(SQRT(ABS(TB*(ZL1-TA))),TB) C C ***** PROCEDURES ****************************************** C IF (WHERE.NE.0) GO TO 2900 C C CHECK FOR TRIVIAL CASE C PART=1 RELEPS=0.0e0 10 ERRT(1)=0.0e0 c--D Next line special: S => D, X => Q RESULT(1)=0.0D0 XCDOBT(1)=.FALSE. EDUE2A=0.0e0 EDUE2B=0.0e0 IF (AINIT.EQ.BINIT) GO TO 2320 C C SETUP FOR COMPLETE INTERVAL. C SET INITIAL VALUES. C IFLAG(1)=0 XJ = 0.E0 XJP = 0.E0 FAT(1) = 0.E0 FAT(2) = 0.E0 ERRC = 0.0e0 ERRT(2)=0.0e0 START(1)=AINIT END(1)=BINIT EPSR=MAX(EPSO,ABS(RELEPS)) EPSMAX=XEPS*EPSR NSUB=0 TA=START(1) TLEN=ABS(BINIT-AINIT) STEP(1)=1.0e0 WORRY(1)=END(1) INIT=.TRUE. DID1=.FALSE. HAVDIF=.FALSE. EPS=EPSR IF (TALOC.NE.0) THEN I=ABS(TALOC) TA=WORK(I) NSUB=2 IF (TALOC.LT.0) NSUB=4 IF (SIGN(1.0e0,TA-START(1))*(END(1)-TA).GT.0.0e0) THEN IF (TA.NE.START(1)) THEN C TA IS INSIDE THE INTERVAL. SUBDIVIDE IMMEDIATELY. PART=2 STEP(2)=-1.0e0 END(2)=START(1) WORRY(2)=END(2) c--D Next line special: S => D, X => Q RESULT(2)=0.0d0 XCDOBT(2)=.FALSE. START(1)=TA START(2)=TA TB=END(2)-TA STEP(1)=SINTSM(TA) NSUBSV=NSUB TASAVE=TA EPS=EPSR*ABS(END(2)-START(2))/TLEN END IF ELSE C TA IS OUTSIDE THE INTERVAL. IF (ABS(TA-START(1)).GE.ABS(TA-END(1))) THEN C TA IS NEARER BINIT. TURN THE INTERVAL AROUND. X=END(1) END(1)=START(1) START(1)=X STEP(1)=-STEP(1) END IF IF (TA.NE.START(1)) THEN TB=END(1)-TA START(1)=TINCR(START(1)) IF (NSUB.EQ.4) START(1)=TINCR(START(1)) END IF END IF END IF 40 RNDC=FER EPNOIZ=ENINF FATAS=.FALSE. DISCF=0 DISCHK=0 FSAVED=.FALSE. KMAX=KMAXF KAIMT=MIN(MAX(5,KMAXF),MAXK-2) ALOCAL=START(PART) TB=END(PART)-TA C GO COMPUTE MINIMUM STEPSIZE 50 DELMIN=SINTSM(ALOCAL) 60 BLOCAL=END(PART) TEND=BLOCAL IEND=.TRUE. 70 FATBS=.FALSE. KMIN=0 GO TO 80 75 EPS=EPSR*DELTA/TLEN INIT=.FALSE. C C INITIALIZE VALUES NECESSARY TO INTEGRATE ONE PANEL. C 80 LOCAL(3)=ALOCAL LOCAL(4)=BLOCAL IF (NSUB.NE.0) THEN LOCAL(3)=TDECR(LOCAL(3)) LOCAL(4)=TDECR(LOCAL(4)) IF (NSUB.NE.2) THEN LOCAL(3)=TDECR(LOCAL(3)) LOCAL(4)=TDECR(LOCAL(4)) END IF END IF DELTA=BLOCAL-ALOCAL DIFF=0.5e0*DELTA ABSDIF=ABS(DIFF) PX=DIFF DELTA=ABS(DELTA) SUM=ALOCAL+0.5e0*(BLOCAL-ALOCAL) EPSS=MAX(EPS,EPSMAX) DID1=START(PART).NE.ALOCAL ROUNDF=.FALSE. RE=0.2e0 REPROD=1.0e0 IP=1 IOLD=0 K=2 IF (IPRINT.GT.1) CALL SINTO (1,WORK) C C APPLY 1-POINT GAUSS FORMULA (MIDPOINT RULE). C ABSCIS=SUM WHERE=1 GO TO 2850 90 FUNCT(1)=FNCVAL PF1=FUNCT(1) PF2=FUNCT(1) PACUM=FUNCT(1)*DELTA PAACUM=ABS(PACUM) ERRI=PACUM ERR=ERRI PEPSMN=ABS(PACUM*RNDC) IF (FEA.NE.0) PEPSMN=PEPSMN+ABS(DELTA*T*WORK(FEA)) C C GO ON TO NEXT FORMULA. C ACUM=FUNCT(1)*P(1) AACUM=ABS(ACUM) EPSMIN=0.0e0 GO TO 130 95 DISCHK=MAX(DISCHK,1) 100 K=K+1 PACUM=ACUM PEPSMN=EPSMIN EPSMIN=ENZER c--D Next line special: S => D, X => Q ACUM=0.0D0 AACUM=0.0E0 C C COMPUTE CONTRIBUTION FROM PREVIOUSLY CALCULATED FUNCTION VALUES. C KK=KORECT(K-2) 110 J=MOD(KK,400) KK=KK/400 J1=MOD(J,20) J2=J/20 DO 120 J=J1,J2 IP=IP+1 ACUM=ACUM+P(IP)*FUNCT(J) 120 CONTINUE IF (KK.NE.0) GO TO 110 C C COMPUTE CONTRIBUTION FROM NEW FUNCTION VALUES. C 130 INEW=IOLD+1 IOLD=IOLD+INEW J2=MOD(FSTORE(K-1),100) J1=FSTORE(K-1)/100 J=INEW 140 X=P(IP+1)*DIFF X1=BLOCAL-X X2=ALOCAL+X ABSCIS=X1 ERRF=0.0e0 WHERE=2 IF (J.EQ.INEW) GO TO 2840 C THE FIRST FUNCTION VALUE REQUESTED FOR EACH FORMULA IS THE C ONE NEAREST THE EDGE OF THE PANEL. GO TO 2850 150 F1=FNCVAL IF (FEA.NE.0) ERRF=ABS(T*WORK(FEA)) ABSCIS=X2 WHERE=3 IF (J.EQ.INEW) GO TO 2840 GO TO 2850 160 F2=FNCVAL IF (FEA.NE.0) ERRF=ERRF+ABS(T*WORK(FEA)) C UPDATE TOLERANCE IF A FUNCTION VALUE WAS REQUESTED AT THE EDGE OF C THE INITIAL INTERVAL. IF (J.EQ.INEW) THEN NEWEPS=.FALSE. IF (ALOCAL.EQ.AINIT) THEN EDUE2A=ABS(F1*ERRINA) NEWEPS=.TRUE. ELSE IF (ALOCAL.EQ.BINIT) THEN EDUE2B=ABS(F1*ERRINB) NEWEPS=.TRUE. END IF IF (BLOCAL.EQ.AINIT) THEN EDUE2A=ABS(F2*ERRINA) NEWEPS=.TRUE. ELSE IF (BLOCAL.EQ.BINIT) THEN EDUE2B=ABS(F2*ERRINB) NEWEPS=.TRUE. END IF IF (NEWEPS) THEN EPS=EDUE2A+EDUE2B EPS=MAX(EPSR-EPS,0.1e0*(ERRT(1)+ERRT(2)+EPS))*DELTA/TLEN END IF END IF C COMPUTE EPSMIN X1=(ABS(X1*(F1-PF1))+ABS(X2*(F2-PF2)))*EMEPS X1=X1/ABS(X-PX) PX=X PF1=F1 PF2=F2 IP=IP+2 S=ABS(F2) IF (NSUB.NE.0) S=MAX(S,ABS(ANSWER*(TA/(X2-TA)))) EPSMIN=EPSMIN+P(IP)*(X1+(ABS(F1)+S)*RNDC+ERRF) AACUM=AACUM+P(IP)*(ABS(F1)+ABS(F2)) F2=F2+F1 ACUM=ACUM+P(IP)*F2 IF (J1.LE.J2) THEN FUNCT(J1)=F2 J1=J1+1 END IF IF (J.LE.7) FUNCT(17+J)=F1 J=J+1 IF (J.LE.IOLD) GO TO 140 c--D Next line special: S => D, X => Q ACUM=ABSDIF*ACUM+0.5D0*PACUM AACUM=MAX(EDELM1,ABSDIF*AACUM+0.5E0*PAACUM) EPSMIN=MAX(EDELM1,ABSDIF*EPSMIN+0.5e0*PEPSMN) EPSR=MAX(EPSR,EPSMIN) C C CHECK FOR CONVERGENCE. C EP=MAX(ABS(ERRI),0.125e0*EPSMIN) ERRI=ACUM-PACUM ERR=ABS(ERRI) EXTRA=0.0e0 REP=RE IF (INIT) THEN IF (RELTOL.NE.0.0e0) THEN IF (RELEPS.LE.0.0e0) THEN RELEPS=ABS(ACUM)-ERR RELEPS=-ABS(WORK(RELTOL+2))*MAX(RELEPS,0.0e0) EPSR=MAX(EPS,ABS(RELEPS)) EPSMAX=XEPS*EPSR EPSS=MAX(EPSMAX,EPSR) END IF END IF END IF EPS=MAX(EPSS,(EPSMIN/AACUM)**RELOBT*AACUM) RE=ERR/EP TPS=RE+RE IF (NSUB.NE.0) THEN TPS=0.25e0-1000.0e0*ABS((TA-LOCAL(3))/(LOCAL(4)-LOCAL(3)))**NSUB TPS=2.0e0*MAX(RE,RECORR(NSUB-1)*TPS) END IF REPROD=MIN(RE*REPROD,1.0e0) FAIL=.TRUE. SEARCH=1 WHERE=2 HAVDIF=HAVDIF.AND.K.GT.4 J=MIN(0,K-KAIMT) IF (J.LT.0) TPS=MAX(TPS,2.0e0*MIN(5.0e0*TPS,MIN(REP,0.125e0))) TP=EPSMIN*(FUDGE(J+5)/RNDC) IF (ERR.GT.TP) THEN C IF ERROR IS LARGE RELATIVE TO THE OBTAINABLE PRECISION, INCREASE C THE ERROR ESTIMATE. (THIS TEST IS CONSERVATIVE WHEN K IS SMALL) TP1=0.25e0/MAX(1.0e0,25.0e0*REP*REP) 190 ERR=4.0e0*ERR TP=32.0e0*TP TPS=TPS*(1.0e0+TPS)/(TP1+TPS) IF (ERR.GT.TP) GO TO 190 C END OF INCREASING ERROR ESTIMATE WHEN RELATIVE ERROR IS LARGE. IF ((ERR*RNDC).GE.EPSMIN) THEN IF (ERR.GT.EPS) IF (DISCHK) 220,220,230 C GOT ACCURACY REQUESTED, INTEGRAND ASSUMED TO BE ALL NOISE. C DO NOT BLINDLY ACCEPT RESULT IF ON INITIAL INTERVAL. IF (K-2) 220,220,230 END IF END IF ERR=ERR*TPS C DO NOT ACCEPT RESULT WITHOUT LOOKING AT DIFFERENCE LINE IF C WORKING ON THE INITIAL INTERVAL OR THE LEFT BOUNDARY WAS A C STOP POINT OR AN END-POINT SINGULARITY, OR IF RE*REP .GT. RELOOK. IF (DISCHK.GT.0) GO TO 230 IF (ERR.GT.EPS) GO TO 220 IF (K.NE.KAIMT) GO TO 230 IF (DISCHK.LT.0) GO TO 230 IF (RE*REP.GT.RELOOK) GO TO 230 IF (INIT) GO TO 230 IF (RE.GE.1.0e0) GO TO 220 WHERE=3 GO TO 230 220 WHERE=1 230 IF (IPRINT.GT.2) CALL SINTO (2,WORK) EP=ERRC ERRC=ERR ERRI=ABS(ERRI) GO TO (240,250,2250), WHERE 240 IF (K.LT.MIN(KMAX,KAIMT)) GO TO 100 C C THE ERROR IS NOT SMALL ENOUGH, OR IT SEEMS SMALL ENOUGH BUT WE C HAVE A REASON TO BE SUSPICIOUS. C 250 IF (HAVDIF) GO TO 500 C COMPUTE ABSCISSAE AND UNSCRAMBLE FUNCTION VALUES INTO C USABLE ORDER. FT(9)=FUNCT(1) XT(9)=SUM J=1 L=8 KK=MIN(K,4) DO 270 J1=2,KK INC=L L=L/2 DO 260 I=L,7,INC J=J+1 J2=IH(I) BETA=DIFF*P(J2) XT(I+1)=ALOCAL+BETA XT(17-I)=BLOCAL-BETA FT(I+1)=FUNCT(J)-FUNCT(16+J) FT(17-I)=FUNCT(16+J) 260 CONTINUE 270 CONTINUE ENDPTS=1 IF (K.NE.2) GO TO 300 ENDPTS=3 C COMPUTE F AT (OR NEAR) THE ENDS OF THE PANEL. WHERE=4 WHERE2=1 GO TO 2810 280 FT(1)=FATA XT(1)=X WHERE2=2 GO TO 2810 290 FT(17)=FATB XT(17)=X 300 WHERE=0 C WHERE=0 TELLS SINTDL TO FORM AND ANALYZE DIFFERENCE LINES. GO TO 370 C A NEW FUNCTION VALUE HAS BEEN COMPUTED. 310 WHERE=3 C NSUB HAS BEEN CHANGED OR THE DIFFERENCE LINES HAVE BEEN OTHERWISE C MUTILATED. UPDATE THE DIFFERENCE LINES. 350 IF (KDIM.EQ.1) ERRF=0.0e0 IF (FEA.NE.0) ERRF=ABS(T*WORK(FEA)) C WHERE TELLS SINTDU WHETHER TO ADD A FUNCTION VALUE IN THE C MIDDLE, ON THE ALOCAL END OR THE BLOCAL END, OR TO RECOMPUTE THE C DIFFERENCE LINES. CALL SINTDU 360 WHERE=1 C WHERE .NE. 0 TELLS SINTDL TO ANALYZE DIFFERENCE LINES THAT C HAVE ALREADY BEEN FORMED. 370 CALL SINTDL (WORK) 500 KONVRG=IC1+LENDT-IC2 C C EXAMINE THE ISTOP ARRAY. C 510 JUMPS=1 IF (ISTOP(2,2)+18-ISTOP(2,1).EQ.0) GO TO 530 JUMPS=4 IF (ISTOP(2,2)*(18-ISTOP(2,1)).NE.0) GO TO 530 JUMPS=3 IF ((LENDT-ISTOP(2,1))*(1-ISTOP(2,2)).EQ.0) JUMPS=2 GO TO 530 C IF WE WERE NOT DOING A JUMP SEARCH AND FOUND JUMPS, START ONE. 520 ESOLD=-1.0e0 COUNT=0.0e0 XJUMP=XJUMPS SEARCH=5 C SET A CAUTION POINT. WORRY(PART)=XT(LENDT) IF (PART.EQ.2.AND.STEP(1)*STEP(2).LT.0.0e0) WORRY(1)=XT(1) 530 CONTINUE C DO BLOCK J1=ISTOP(2,1) J2=ISTOP(2,2) ERR=ERRC GO TO (540,2150,770,860,800,535), SEARCH c One hopes SEARCH==6 and JUMPS==4 can't happen 535 GO TO (2160,620,970,520), JUMPS 540 COUNT=0.0e0 IF (JUMPS.EQ.1) GO TO 1920 IF (K.LE.2) GO TO 95 IF (DELTA.LT.DELMIN+DELMIN) GO TO 2170 IF (DISCHK.EQ.0) DISCHK=-1 IF (JUMPS.EQ.3) GO TO 1920 IF (JUMPS.NE.2) GO TO 790 C 88888888888 C END POINT JUMP. C IF ((TA-ALOCAL)*SIGN(1.0e0,BLOCAL-TA).GT.0.0e0) * IF (ENDPTS-2) 2750,2760,560 560 CONTINUE IF (NSUB.NE.0) THEN C NSUB CAN ONLY BE NON-ZERO WHEN .NOT.DID1 IF THE JUMP WAS EXPECTED. IF (.NOT.DID1) XJ=ABS(XJ) END IF IF (XJ.GE.0.0e0) THEN IF (XJP.GT.0.0e0) THEN IF (J2.EQ.0) GO TO 630 IF (RE.GT.1.0e0) GO TO 670 ERR=ERR+EXTRA IF (ERR.LE.EPS) GO TO 750 c?? Do we need to ask (re .lt. rep) and compute two ways? c## S=RE*(RE/REP)*ERR c## IF (S.GT.EPS) GO TO 670 IF (RE*RE*ERR.GT.REP*EPS) GO TO 670 C IF KAIMT IS EVER PERMITTED .GT. MAXK-2, MAXK-K MUST BE TESTED. IF (K-KAIMT-1) 2740,2740,670 END IF END IF IF (ENDPTS-2) 2750,2760,600 600 CONTINUE IF (SEARCH.EQ.1) THEN C IF THE ERROR IN THE FUNCTION AT THE BOUNDARY OF THE PANEL IS C LARGER THAN THE DIFFERENCE BETWEEN THE FUNCTION AND THE C VALUE PREDICTED BY EXTRAPOLATING THE DIFFERENCE LINE, C REPLACE THE FUNCTION BY THE EXTRAPOLATED VALUE, AND DECREASE C THE ERROR BY THE CORRECTION. THIS DOESN'T AFFECT THE ERROR C IN THE INTEGRAL BECAUSE FUNCTION VALUES AT THE ENDS DON'T C CONTRIBUTE TO THE INTEGRAL. ENDPHI(2)=LENDT DO 610 WHERE2 = 1, 2 IF (FATS(WHERE2)) THEN IF (ERRAT(WHERE2).GT.ABS(PHIS(ENDPHI(3-WHERE2),WHERE2))) THEN FAT(WHERE2)=FAT(WHERE2)-PHIS(ENDPHI(3-WHERE2),WHERE2) FT(ENDPHI(WHERE2))=FAT(WHERE2) ERRAT(WHERE2)=ERRAT(WHERE2)-ABS(PHIS(ENDPHI(3-WHERE2),WHERE2)) GO TO 310 END IF END IF 610 CONTINUE END IF DISCHK=MAX(DISCHK,1) IF (RE.LE.1.0e0) THEN IF (J2.EQ.0) GO TO 630 IF (ERR.LE.EPS) GO TO 760 END IF IF (K.LT.4) GO TO 100 C END BLOCK 620 IF (J2.NE.0) GO TO 670 C C IF WE SEE A JUMP NEAR BLOCAL AND WE ARE AT THE END OF A PART, C TURN THE INTERVAL AROUND. IF WE ARE NOT AT THE END OF A C PART, TAKE HALF THE CURRENT INTERVAL AND START OVER. C IF (SEARCH.NE.1) GO TO 1440 630 CONTINUE C DO BLOCK IF (IEND) THEN IF (NSUB.EQ.0) GO TO 640 IF (ABS(ALOCAL-TA).GE.0.125e0*ABS(TB)) THEN IF (.NOT.DID1.OR.XJP.LE.0.1e0) GO TO 640 ISTOP(2,1)=18 C IF THE JUMP IS WEAK, PRETEND ITS NOT THERE. GO TO 510 END IF END IF DELTA=0.5e0*DELTA GO TO 2610 C END BLOCK 640 CONTINUE if (NSUB.NE.0) then CALL SINTNS (NSRB) GO TO 640 END IF START(PART)=BLOCAL BLOCAL=ALOCAL END(PART)=ALOCAL ALOCAL=START(PART) TA=ALOCAL TB=END(PART)-TA IF (DISCHK.LE.0) DISCHK=-2 CALL SINTNS (NSIA) C TAKE HALF OF THE INTERVAL TO AVOID GETTING IN A LOOP. DELTA=0.5e0*DELTA TEND=ALOCAL+SIGN(DELTA,TB) STEP(PART)=-STEP(PART) GO TO 2610 C C IF WE SEE AN EXPECTED AND BELIEVABLE JUMP NEAR ALOCAL, WE C SHOULD INCREASE NSUB. OTHERWISE, GO INTO THE INTERVAL SEARCH. C 670 CONTINUE IF (XJ.GT.0.0e0.AND.NSUB.LT.NSUBMX.AND.(INIT.OR.LENDT-IC2.GE.2)) *THEN IF (RE.GE.1.0e0.AND.DID1.AND.DISCHK.NE.-2) THEN IF (K.LT.KAIMT) GO TO 100 END IF IF ((TA-ALOCAL)*(BLOCAL-TA).GT.0.0e0) THEN TA=ALOCAL TB=END(PART)-TA END IF KAIMT=MIN(MAX(5,KMAXF),MAXK-2) KMAX=KMAXF IF (SEARCH.NE.1) THEN BLOCAL=XT(LENDT) IEND=.FALSE. END IF CALL SINTNS (NSIB) IF (.NOT.IEND.AND.(PART.NE.2.OR.DID1)) THEN C SET A CAUTION POINT. IF (ABS(BLOCAL-START(PART)).LT.ABS(TEND-START(PART))) TEND=BLOCAL END IF IF (SEARCH.EQ.1) THEN IF (IC2.LT.LENDT-2.OR.K.GT.3) THEN DELMIN=SINTSM(ALOCAL) EPS=EPSS GO TO 70 END IF BLOCAL=ALOCAL+0.5e0*(BLOCAL-ALOCAL) DELTA=0.5e0*DELTA END IF IEND=.FALSE. DELMIN=SINTSM(ALOCAL) GO TO 75 END IF C WE CAN NOT INCREASE NSUB. START A SEARCH FOR THE CORRECT C INTERVAL LENGTH. IF (ENDPTS-2) 2750,2760,2155 C C ADD POINTS AT THE END NEAR THE JUMP UNTIL WE ARE SURE THE C JUMP IS NOT INSIDE THE INTERVAL. THIS PROCESS CONTINUES UNTIL C THE JUMP VANISHES, MOVES INSIDE THE INTERVAL, OR THE ERROR C ESTIMATED ACROSS THE JUMP IS SUFFICIENTLY SMALL. IF ALOCAL = TA C (WHICH MEANS WE ARE FAIRLY CERTAIN THE SINGULARITY IS NOT IN THE C INTERVAL), NSUB IS NOT ZERO, AND THE JUMP IS EXPECTED, ACCEPT THE C ANSWER WITHOUT DOING A SEARCH. C 750 IF (ENDPTS-2) 2750,2760,760 760 SEARCH=-3 770 GO TO (1190,780,980,870), JUMPS 780 IF (NSUB.EQ.0) GO TO 980 IF (XJ.LE.0.0e0) GO TO 980 IF (ALOCAL-TA) 980,2240,980 C C JUMPS WERE NOTICED IN THE ORIGINAL DIFFERENCE LINE ON THIS PANEL. C 790 IF (K+JUMPS.LT.6) GO TO 100 IF (MAX(J1,J2).LE.NSUB+2) GO TO 2130 EPSS=EPS IF (NSUB.NE.0) GO TO 810 SEARCH=-4 XJUMP=XJUMPS GO TO 940 800 IF (NSUB.EQ.0) GO TO 850 C IF WE HAVE NOT DONE ANYTHING, REDUCE DELTA INSTEAD OF SEARCHING. C THERE WAS A REASON FOR SETTING NSUB NON-ZERO. 810 IF (J2.LT.2) J2=J1 J=MIN(J1,J2) IF (ABS(2.0e0*(XT(1)-TA)).LE.ABS(XT(J)-TA)) THEN WORRY(PART)=XT(J) DELTA=0.75e0*ABS(XT(J)-ALOCAL) DISCHK=2 GO TO 2610 END IF C REDUCE NSUB TO ZERO BEFORE CONTINUING THE SEARCH. X1=REAL(NSUB) J=1 DO 820 I=1,LENDT T=(XT(I)-TA)/TB IF (NSUB.NE.2) T=T**(NSUB-1) FT(I)=FT(I)/(X1*T) XT(I)=TA+T*(XT(I)-TA) IF (I.NE.1) THEN IF (XT(I).NE.XT(I-1)) THEN J=J+1 XT(J)=XT(I) FT(J)=FT(I) END IF END IF 820 CONTINUE LENDT=J SEARCH=6 C FINISH DECREASING NSUB 840 IF (NSUB.NE.0) THEN CALL SINTNS (NSRB) GO TO 840 END IF GO TO 310 C ARE JUMPS STILL WITHIN RANGE OF PREVIOUS JUMPS? 850 SEARCH=4 860 IF (JUMPS.LT.4) GO TO 870 IF (J1.LT.MIN(J1OLD,J2OLD)) GO TO 870 IF (J2.LE.MAX(J1OLD,J2OLD)) GO TO 940 870 IF (NFEVAL.LE.NFJUMP+3) GO TO 2690 GO TO 310 C C THE DIFFERENCE LINE HAS BEEN BALANCED. CONTINUE THE JUMP SEARCH C OR COMPUTE THE INITIAL INTERVAL. C 910 CONTINUE IF (SEARCH.NE.4) THEN IF (JUMPS-3) 920,930,915 915 IF (J2-(NSUB+2)) 930,520,520 920 I=6-NSUB/2 J=I IF (IC1.GT.J) J=MIN(I+1,6) KK=4 KMAX=MAX(KMAX,5) DELTA=ABS(XT(J+1)-ALOCAL) IF (NSUB.EQ.0) 1 DELTA=MIN(DELTA,ABS(XT(6)-ALOCAL)+2.0e0*ABS(XT(6)-XT(5))) DELTA=Z(1)*DELTA TPS=EMEPS*EPSMIN/RNDC C TPS DOESN'T DEPEND ON TOLERANCE; INIT. INTERVAL DOESN'T EITHER. IF (NSUB.NE.0) TPS=MIN(EPSMIN,TPS) ALPHA=(DELTA*ERRCF(2)*(ABS(PHI(7))+ABS(PHI(8))))/TPS IF (ALPHA.LT.1.0e0) KK=5 GO TO 2480 930 CONTINUE END IF IF (NFEVAL.GT.NFJUMP+1) GO TO 310 IF (JUMPS.EQ.1) GO TO 1190 C C JUMPS HAVE NOT ESCAPED. C 940 IF (J1.GE.J2) GO TO 970 945 CONTINUE IF (J2-J1.LE.1) THEN I=J1 J1=J2 J2=I ELSE C JUMPS DO NOT OVERLAP, CHOOSE ONE ON WHICH TO CONVERGE. IF (PART.NE.1.AND.STEP(1)*STEP(2).LE.0.0e0) THEN J1=18 ELSE J2=0 END IF END IF C JUMPS OVERLAP NOW. 970 J1OLD=J1 IF (J1OLD.EQ.LENDT) J1OLD=LENDT-1 980 J2OLD=J2 IF (J2OLD.EQ.1) J2OLD=2 C WE ONLY GET HERE IF WE HAVE AT LEAST 1 JUMP. IF (J1.EQ.18) J1=J2+1 IF (J2.EQ.0) J2=J1-1 IF (MIN(ISTOP(2,1),LENDT+1-ISTOP(2,2)).LT.(LENDT+2)/3) GO TO 1350 IF (ABS(SEARCH).NE.4) GO TO 990 c IF (JUMPS.NE.4) GO TO 1350 c The above test prevented getting out of the search early on c rapidly decaying integrands. IF (ISTOP(2,1).LT.ISTOP(2,2)) GO TO 1350 C C ESTIMATE ERROR ACROSS JUMP. C 990 WHERE=0 1000 CONTINUE IF (ABS(SEARCH).NE.4.AND.ABS(SEARCH).NE.3) THEN IF (JUMPS.GE.3) THEN c IF (WHERE.EQ.0) GO TO 1350 c The above test prevented estimating error if JUMPS == 3 INC=LENDT INC2=1 END IF ELSE if (jumps .eq. 3) then inc = min(lendt,j1old+1) inc2 = max(1,j2old-1) else INC=MIN(LENDT,J1+1) INC2=MAX(1,J2-1) end if END IF TPS=0.0e0 F1=-1.0e0 DO 1010 I=INC2,INC-1 F2=MAX(ABS(XT(I+1)-XT(I)),ABS(EEPSM8*XT(I))) F2=ABS((FT(I+1)-FT(I))*F2) C F1 IS THE ERROR COMMITTED ON THE TRAPEZOID WITH THE LARGEST ERR F1=MAX(F1,F2) TPS=TPS+F2 1010 CONTINUE TPS=0.5e0*TPS IF (MAX(XJ,XJP).GT.0.05e0) TPS=4.0e0*TPS X=ABS(0.5e0*(XT(INC)-XT(INC2))) IF (INC2.GE.INC-2) THEN IF (INC2.GT.1) TPS=TPS+ABS(X*(FT(INC2)-FT(INC2-1))) IF (INC.LT.LENDT) TPS=TPS+ABS(X*(FT(INC+1)-FT(INC))) END IF C COMPUTE THE ALLOWABLE ERROR. C DO BLOCK IF (WHERE.EQ.0) THEN C THE ABSCISSAE DID NOT COALESCE DURING INSERTION OF A POINT C IN THE DIFFERENCE LINES S=MAX(0.25e0*EPS,EPSMAX) C TRY TO GET AT LEAST 3 DIGITS IF THE JUMPS LOOK REAL: if (jumps .ne. 3) S=MIN(S, 1.0e-3*(EPSMIN/RNDC)) IF (SEARCH.NE.6) THEN IF (ABS(SEARCH)-3) 1060,1080,1100 1060 CONTINUE END IF C ABS(SEARCH)=2 OR 6 HERE TPS=10.0e0*TPS IF (IPRINT.GT.3) CALL SINTO (6,WORK) IF (TPS.LE.0.01e0*S) IF (ERR-S) 2240, 2240, 1150 if (jumps .ge. 3) GO TO 1350 IF (NSUB.NE.0) IF (ABS(XT(4)-ALOCAL)-DELMIN) 1150,1150,1440 IF (ABS(XT(9)-ALOCAL).GE.8.0e0*DELMIN) GO TO 1440 C SET STEPSIZE TO DELMIN. DELTA=0.0e0 GO TO 2620 C ABS(SEARCH)=3 HERE 1080 CONTINUE IF (TPS.LE.EPS) THEN IF (XJ.GT.0.0e0) THEN IF (XJP.GT.0.0e0) THEN IF (RE.LE.1.0e0) THEN IF (NSUB.NE.0) THEN C AUGMENT EPSMIN DUE TO NOISE INTRODUCED BY C T**2N TRANSFORMATION. X1=ALOCAL-TA IF (X1.EQ.0.0e0) X1=XT(2)-TA X1=ABS(TB/X1) IF (NSUB.NE.2) X1=X1**(NSUB-1) EPSMIN=EPSMIN+ABS(RNDC*FT(1)*TA*X1) END IF GO TO 2240 END IF END IF END IF IF (ISTOP(2,2).NE.0.AND.TPS.LE.S) GO TO 1150 END IF 1100 CONTINUE IF (SEARCH.LE.0) THEN IF (SEARCH.EQ.-3) GO TO 1120 C THE JUMP SEARCH IS JUST STARTING. IF THE SUM OF THE C ERROR ESTIMATED BY THE QUADRATURE STEP, AND THE ERROR C ACROSS THE JUMP IS LESS THAN THE TOLERANCE FOR THE PANEL, C ACCEPT THE ANSWER. ERR=TPS+ERRC IF (ERR.GT.EPS) GO TO 1120 DISCF=1 S=EPS GO TO 1180 END IF END IF IF (ESOLD.GT.0.0e0) COUNT=COUNT+MIN((F1-ESOLD)/ESOLD,1.0e0) C END BLOCK 1120 CONTINUE ESOLD=F1 IF (IPRINT.GT.3) CALL SINTO (6,WORK) SEARCH=ABS(SEARCH) IF (WHERE.EQ.0) THEN C THE ABSCISSAE DID NOT COALESCE DURING INSERTION OF A POINT IN C THE DIFFERENCE LINES. IF (SEARCH.EQ.3) GO TO 1350 IF (TPS-S) 2670,2670,1350 END IF IF (SEARCH.EQ.3) IF (COUNT) 2237,2237,2190 IF (COUNT.LE.0.0e0) GO TO 1150 LOCAL(3)=XT(INC) GO TO 2190 C C THE ESTIMATED ERROR ACROSS THE JUMP IS SMALL ENOUGH TO ACCEPT. C ESTIMATE THE INTEGRAL USING THE COMPOSITE TRAPEZOID RULE. C THIS IS SATISFACTORY SINCE THE PRESENCE OF THE JUMP NEGATES C ANY POSSIBILITY OF HIGH-ORDER CONVERGENCE. C 1150 EPS=S 1160 ERR=TPS ERRI=ERR c--D Next line special: S => D, X => Q ACUM=0.0D0 X1=0.0e0 X2=0.0e0 J=INC-1 C DO BLOCK 1170 S=ABS(XT(J+1)-XT(J)) C COMPUTE EPSMIN ALSO F1=ABS(FT(J)) F2=ABS(FT(J+1)) X1=X1+ABS((XT(J+1)+XT(J))*(F2-F1)) X2=X2+S*(F2+F1) ACUM=ACUM+(FT(J+1)+FT(J))*S J=J-1 IF (J.GE.INC2) GO TO 1170 C END BLOCK EPSMIN=0.5e0*(X1*EMEPS+X2*RNDC) c--D Next line special: S => D, X => Q ACUM=0.5D0*ACUM IEND=.FALSE. HAVDIF=.FALSE. C SET K SMALL FOR STEPSIZE SELECTION PURPOSES. C K=2 1180 FAIL=.FALSE. I=INC J=INC2 IF (SEARCH-4) 1240,1230,1240 C C IF BOTH JUMPS VANISH, AND THE ORIGINAL ERROR WAS SUFFICIENTLY C SMALL, AND ENOUGH OF THE INTERVAL IS BEING EXAMINED, ACCEPT C THE ANSWER. C 1190 IF (ERRC.LE.EPS.AND.10.0e0*ABS(XT(LENDT)-XT(1)).GT.DELTA) 1 GO TO 2237 C C COMPUTE A NEW PANEL TO INTEGRATE. C DISCHK=2 IF (PART.NE.1) THEN IF (STEP(1)*STEP(2).LE.0.0e0) THEN IF (J1.EQ.18) J1=0 J=MAX(2,MAX(J2,J1-2)) I=LENDT IF (JUMPS.NE.1) I=MIN(LENDT,J+3) GO TO 1240 END IF END IF IF (J2.EQ.0) J2=18 I=MIN(LENDT-1,MIN(J1,J2+2)) J=1 IF (JUMPS.NE.1) J=MAX(1,I-3) GO TO 1250 C C GET THE BOUNDARIES OF THE SUBDIVISION. C 1230 IF (PART.EQ.1.OR.STEP(1)*STEP(2).GT.0.0e0) GO TO 1250 1240 L=I I=J J=L X=XJ XJ=XJP XJP=X 1250 X1=XT(I) X2=XT(J) FATA=FT(I) FATB=FT(J) FATAS=.TRUE. FATBS=.TRUE. WHERE=5 IF (.NOT.FAIL) THEN IF (NSUB.EQ.0) THEN DISCX(1)=X1 DISCX(2)=X2 ELSE DISCX(1)=TA+(X1-TA)*((X1-TA)/TB)**(NSUB-1) DISCX(2)=TA+(X2-TA)*((X2-TA)/TB)**(NSUB-1) END IF IF (SEARCH.EQ.-4) GO TO 2240 DELTA=ABS(X1-X2) IF (IPRINT.GT.1) CALL SINTO (7,WORK) IF (SEARCH.NE.4) THEN BLOCAL=X2 ABSDIF=0.5e0*DELTA IF (MIN(ABS(XJ),ABS(XJP)).LE.1.0E-3) DISCF=1 IF (SEARCH.EQ.2.AND.NSUB.NE.0) CALL SINTNS (NSRB) GO TO 2240 END IF END IF C DECREASE NSUB TO ZERO. 1260 IF (NSUB.NE.0) THEN CALL SINTNS (NSRB) go to 1260 END IF IF (FAIL) THEN C C TRY NOT TO CHOOSE A PANEL TO BE INTEGRATED THAT WILL LEAVE C A RIDICULOUSLY SMALL UN-INTEGRATED PIECE. C C DO BLOCK IF (ABS(X2-ALOCAL).GE.1.5e0*DELMIN) THEN IF (ABS(X2-X1).LE.8.0e0*ABS(X1-ALOCAL)) THEN SMIN=SINTSM(X1) IF (ABS(X1-BLOCAL).GE.1.5e0*SMIN) THEN IF (ABS(X2-X1).LE.8.0e0*ABS(X1-BLOCAL)) GO TO 1280 END IF IF (PART.EQ.2) GO TO 1270 IF (.NOT.IEND) GO TO 1270 X1=X2 END IF END IF DELTA=ABS(X1-ALOCAL) FATAS=.FALSE. GO TO 2610 C END BLOCK 1270 CONTINUE X1=BLOCAL FATAS=.FALSE. END IF 1280 CONTINUE FATBS=FATAS C C SPLIT THE PROBLEM INTO TWO PARTS AND INTEGRATE AWAY FROM C THE MIDDLE IN BOTH DIRECTIONS. C C DO BLOCK IF (PART.NE.1) THEN I=2 c## EPSR=EPSR+ERRT(2) IF (STEP(1)*STEP(2).LE.0.0e0) GO TO 1290 END IF END(2)=ALOCAL I=1 C END BLOCK 1290 CONTINUE TLEN=ABS(END(2)-END(1)) FSAVE=FATA FSAVED=FATAS S=START(1) START(1)=X1 IF (FAIL) THEN C Jump disappeared. X1..X2 not yet integrated. START(2)=X1 ! Why was the following ever done ??? ! TA=X1+0.5e0*(X2-X1) TA=X1 TASAVE=TA ELSE C Jump didn't disappear. X1..X2 was integrated using trapezoid START(2)=X2 TA=START(1) TASAVE=START(2) END IF TB=END(2)-TA F1=WORRY(1) WORRY(1)=LOCAL(3-I) IF (ABS(F1-S).GE.ABS(X1-S)) THEN IF (ABS(F1-S).LT.ABS(WORRY(1)-S)) WORRY(1)=F1 END IF WORRY(2)=LOCAL(I) DELTA=ABS(X1-X2) STEP(2)=SIGN(1.0e0,SIGN(1.0e0,TB)*(BINIT-AINIT)) ALOCAL=X1 BLOCAL=X2 PART=2 DID1=.FALSE. INIT=.FALSE. NSUBSV=0 IF (.NOT.FAIL) THEN IF (XJ.GT.0.0e0) NSUBSV=2 IF (XJP.GT.0.0e0) CALL SINTNS (NSIB) END IF ABSDIF=0.5e0*DELTA ERRT(2)=0.0e0 RESULT(2)=0.0e0 XCDOBT(2)=.FALSE. DISCHK=MAX(DISCHK,1) IF (.NOT.FAIL) THEN SMIN=SINTSM(BLOCAL) IEND=ABS(END(2)-BLOCAL).LE.SMIN C SOMEDAY, WE SHOULD USE A TRAPEZOIDAL ESTIMATE ACROSS ANY C TINY PIECE LEFT UN-INTEGRATED DUE TO THE ABOVE STATEMENT. IF (JUMPS .EQ. 3) then DISCF = -2 else DISCF=2 end if S=WORRY(1) IF (XJ.LE.0.0e0) S=END(1) C SET THE RESTARTING STEPSIZE FOR PART 1. SMIN=SINTSM(START(1)) STEP(1)=SIGN(MAX(ABS(S-START(1)),SMIN),STEP(1)) GO TO 2250 END IF C SET THE RESTARTING STEPSIZE FOR PART 1 TO THE STEPSIZE FOR PART 2. SMIN=SINTSM(X1) STEP(1)=SIGN(MAX(DELTA,SMIN),STEP(1)) KMIN=0 GO TO 2615 C C THE ESTIMATED ERROR ACROSS THE INTERVAL IS TOO LARGE. REDUCE C THE INTERVAL LENGTH AND ADD A FUNCTION VALUE IN THE NEIGHBORHOOD C OF THE JUMP. THEN RE-FORM AND RE-EXAMINE THE DIFFERENCES. C 1350 I=J1 J=J2 C DO FOREVER 1360 CONTINUE C DO BLOCK IF (I.NE.LENDT) THEN IF (J.NE.1) THEN IF (JUMPS.GE.3) THEN X1=ABS(XT(I)-XT(I+1)) X2=ABS(XT(J)-XT(J-1)) IF (0.3e0*X2.GE.MIN(X1,ABS(XT(J+1)-XT(J)))) THEN I=J-1 GO TO 1400 END IF IF (0.3e0*X1.GE.MIN(X2,ABS(XT(I)-XT(I-1)))) THEN J=I I=I+1 GO TO 1400 END IF END IF END IF END IF IF (I-J-1) 1390,1450,1410 1390 I=J+1 IF (ABS(XT(J)-XT(J+1)).LT.ABS(XT(J)-XT(J-1))) I=J-1 C END BLOCK 1400 CONTINUE ABSCIS=XT(J)+0.25e0*(XT(I)-XT(J)) GO TO 1490 1410 CONTINUE IF (ABS(XJ).LE.ABS(XJP)) THEN XJ=5.0e0*XJ I=I-1 ELSE XJP=5.0e0*XJP J=J+1 END IF GO TO 1360 C END FOREVER 1440 J=1 1450 I=J+1 ABSCIS=XT(I)+0.5e0*(XT(J)-XT(I)) 1490 L=MAX(I,J) IF (SEARCH.LE.0) THEN SEARCH=-SEARCH C IF (SEARCH.EQ.4) THEN C DETERMINE WHETHER NODES HAVE COALESCED. WE ASSUME THAT NODES C WILL NOT COALESCE WHEN SEARCH = 2. C END IF END IF 1500 IF (ABSCIS.EQ.XT(I).OR.ABSCIS.EQ.XT(J)) THEN X1=XT(I)+0.5e0*(XT(J)-XT(I)) IF (X1.EQ.ABSCIS) GO TO 2680 ABSCIS=X1 go to 1500 END IF 1520 CONTINUE IF (LENDT.GT.7) THEN C DECIDE HOW MUCH OF THE DIFFERENCE LINE TO THROW AWAY. C THROW INDICATES WHICH END TO TRIM - 1=ALOCAL, 2=BLOCAL, 3=BOTH. THROW=0 IF (JUMPS.GE.3) IF (ABS(SEARCH)-4) 1555,1530,1555 IF (SEARCH.EQ.2.OR.SEARCH.EQ.6) IF (LENDT-9) 1720,1680,1630 1530 CONTINUE IF (ISTOP(2,2).GT.ISTOP(2,1)+1) THEN C IF JUMPS DO NOT OVERLAP, THROW AWAY A LOT. IF (PART.NE.1) THEN IF (STEP(1)*STEP(2).LE.0.0e0) THEN J1=ISTOP(2,2)-4 J2=0 IF (J1) 1610,1610,1640 END IF END IF J2=LENDT-ISTOP(2,1)-3 J1=0 IF (J2) 1610,1610,1640 END IF IF (JUMPS.EQ.1) GO TO 1610 1555 J=J2 J2=MAX(MIN(LENDT-L,(IC2-J1-1)/2),0) J1=MAX(MIN(L-1,(J-IC1-1)/2),0) IF (SEARCH.NE.4) THEN IF (LENDT.LT.17) J1=0 END IF C DO BLOCK IF (J2-1) 1580,1570,1640 1570 THROW=2 1580 IF (J1-1) 1600,1590,1640 1590 THROW=THROW+1 1600 IF (THROW.NE.0) IF (THROW-2) 1670,1680,1670 C END BLOCK 1610 CONTINUE IF (LENDT.LT.15) GO TO 1740 IF (JUMPS.NE.1) THEN IF (ISTOP(2,1).EQ.18) GO TO 1680 IF (ISTOP(2,2).EQ.0) GO TO 1670 END IF IF (LENDT+1-J1OLD-J2OLD) 1670,1660,1680 1630 J1=0 J2=LENDT-9 C THROW OUT J1 NODES FROM THE HEAD AND J2 NODES FROM THE TAIL C OF THE DIFFERENCE LINES. 1640 LENDT=LENDT-J1-J2 C FORCE RE-FORMATION OF DIFFERENCES FROM XT AND FT. NFJUMP=NFEVAL-10 IF (J1.EQ.0) GO TO 1740 IF (SEARCH.NE.4) GO TO 520 IF (J1OLD.NE.18) J1OLD=J1OLD-J1 J2OLD=J2OLD-J1 L=L-J1 DO 1050 I=1,LENDT XT(I)=XT(I+J1) FT(I)=FT(I+J1) 1050 CONTINUE GO TO 1740 1660 THROW=3 1670 IF (SEARCH.NE.4) GO TO 520 C THROW OUT X(1) AND F(1), RE-FORM THE DIFFERENCE LINES. KK=1 INC=1 J=1 J1=1 J2=LENDT FATAS=.FALSE. L=L-1 IF (J1OLD.NE.0) J1OLD=J1OLD-1 IF (J2OLD.NE.0) J2OLD=J2OLD-1 GO TO 1690 C THROW OUT X(LENDT) AND F(LENDT), RE-FORM THE DIFFERENCE LINES. 1680 KK=2 INC=-1 J=17+LENDT J1=LENDT J2=1 FATBS=.FALSE. 1690 BETA=1.0e0 I=J1 PHI(J)=PHI(J)-PHI(J+INC) J=J+INC I=I+INC C DO BLOCK 1700 BETA=BETA*((XT(J1+INC)-XT(I+INC))/(XT(J1)-XT(I))) PHI(J)=BETA*(PHI(J)-PHI(J+INC)) J=J+INC I=I+INC IF (I.NE.J2) GO TO 1700 C END BLOCK LENDT=LENDT-1 DO 1710 I=1,LENDT PHIT(I)=PHIT(I+1) IF (INC.GE.0) THEN XT(I)=XT(I+1) FT(I)=FT(I+1) END IF 1710 CONTINUE THROW=THROW-KK IF (THROW.GT.0) GO TO 1680 C THE NEXT STATEMENT OCCASSIONALLY CAUSES A DIFFERENCE LINE TO C BE SHORTENED, BUT NO NEW POINT ADDED. IF (KK+L.EQ.2) GO TO 360 1720 CONTINUE END IF 1740 CONTINUE IF (L.GT.LENDT) GO TO 310 C ADD A POINT IN THE INTERIOR OF THE INTERVAL. X=ABSCIS WHERE=6 GO TO 2850 C C NO JUMPS, OR ONLY ONE INTERIOR JUMP, IN THE ORIGINAL C DIFFERENCE LINE FOR THIS PANEL. C 1920 CONTINUE C DO BLOCK IF (RE.LE.1.0e0) THEN IF (K.LE.5.AND.KAIMT.GT.5.AND.KONVRG.LT.4) 1 ERR=MAX(ERR,ERRI*(ERRI/EPSMIN)) IF (ERR.GT.EPS) THEN IF (K.LT.4) GO TO 100 IF (K.GT.5) GO TO 2060 IF (KONVRG-6) 2060,2060,2130 END IF C C ERROR CRITERIA HAVE APPARENTLY BEEN SATISFIED. LOOK AT THE C DIFFERENCE LINES TO AVOID BEING FOOLED BY ANSWERS FROM C SUCCESSIVE FORMULAE BEING NEARLY EQUAL BUT BOTH WRONG. C IF(KONVRG.LE.LENDT.AND.K.LT.KMIN.AND.PHISUM.LE.5.0e0*PHTSUM) 1 GO TO 100 IF (.NOT.INIT.AND.NSUB.EQ.0) THEN IF (ERR.GT.EPSMIN) THEN I=KAIMT-2 IF (KONVRG.GT.LENDT+LENDT-5) I=KAIMT-3 IF (K-I) 100,1950,1960 1950 IF (ENDPTS-2) 2750,2760,1960 1960 CONTINUE END IF END IF IF (K.LE.4) THEN EXTRA=ERRCF(K-1) IF (K.NE.2) THEN IF (KONVRG.LT.K+K-1+ENDPTS/2) EXTRA=10.0e0*EXTRA END IF EXTRA=EXTRA*DELTA*(ABS(PHI(LENDT-1))+ABS(PHIT(2))) J=(K-2)*(KONVRG-7) IF (J.GT.20) THEN EXTRA=0.0e0 ELSE IF (J.GT.0) THEN EXTRA=EXTRA*(0.1e0**J) END IF IF (JUMPS.GE.3) GO TO 2140 IF (EXTRA.GE.32.0e0*EPSMIN) THEN S=ERR+EXTRA IF (S.GT.EPS) IF (K-KAIMT) 2730,2050,2050 ERR=S END IF END IF IF (JUMPS.EQ.3) GO TO 750 C DO BLOCK IF (DISCHK.LE.0) THEN IF (DISCHK.GE.-1) THEN IF (.NOT.INIT) GO TO 2235 IF (K.GE.4) GO TO 2235 KMAX=MAXK GO TO 2030 END IF END IF C COMPUTE F AT (OR NEAR) THE ENDS OF THE PANEL, UNLESS C THERE IS A JUMP AT THE OTHER END OF THE PANEL. KMAX=KAIMT+1 C END BLOCK 2030 CONTINUE IF (ENDPTS-2) 2750,2760,2235 END IF C C CHECK FOR NOISE IN THE INTEGRAND. C 2050 IF (K.LT.4) GO TO 100 2060 CONTINUE IF (RE.GE.0.01e0) THEN C DO BLOCK IF (RE.LT.1.0e0) THEN IF (NSUB.NE.0.AND. 1 ABS((LOCAL(3)-TA)/(LOCAL(4)-LOCAL(3))) 2 .LE.1.0e0) GO TO 2130 ELSE IF (ERRI.LT.100.0e0*EPSMIN) THEN IF (4.0e0*ERR.GT.EPSMIN) GO TO 2100 IF (.NOT.INIT.AND.DISCHK.LE.0) GO TO 2250 IF (ENDPTS-2) 2750,2760,2250 END IF IF (REPROD.GE.1.0E-4) GO TO 2130 IF (K.LT.KAIMT.AND.ERR.GT.100.0e0*EPSMIN) GO TO 2130 IF (ERRI*REPROD.GT.0.01e0*EPNOIZ) GO TO 2130 IF (RE*ERR.LT.EPS) GO TO 2130 C END BLOCK 2100 CONTINUE IF (ENDPTS-2) 2750,2760,2120 C NOISE IN THE INTEGRAND - RESET THE ROUND-OFF CONSTANT. 2120 RNDC=4.0e0*RNDC*(ERR/EPSMIN) IF (IPRINT.GT.0.and.100.0e0*err.gt.eps) CALL SINTO (8,WORK) EPS=ERR ROUNDF=.TRUE. DELTA=DELTA+DELTA GO TO 2250 END IF C END BLOCK C NO NOISE DETECTED. 2130 IF (DELTA.LT.DELMIN+DELMIN) GO TO 2250 KMAX=MAXK IF (DISCHK.GT.0.OR.DISCHK.LT.-1) KMAX=KAIMT+1 C IF WE HAVE A LOT OF CONVERGENCE, WE SHOULD MAKE IT. C IF WE DONT MAKE IT, SOMETHING FISHY IS HAPPENING. C ENDPTS IS INCLUDED BECAUSE CONVERGENCE IS EASIER WHEN THE END C POINTS HAVE BEEN ADDED TO THE DIFFERENCE LINES. IF (KONVRG.GE.ENDPTS+7.AND.RE.GT.1.0E-3) KMAX=5 IF (10.0e0*RE*RE*ERR.LT.REP*EPS) KMAX=MIN(MAXK,MAX(KMAX,K+1)) IF (JUMPS.LT.3) THEN IF (K.LT.KMAX) THEN IF (DISCHK.GE.-1) THEN IF (ABS(WORRY(PART)-START(PART))-DELMIN.GT.ABS(BLOCAL-START(PART)) 1.OR.DID1) GO TO 2730 IF (K.LT.KAIMT) THEN S=RE*ERR IF (K.LT.KAIMT-1) S=S*RE/REP IF (S.LE.EPS) GO TO 2730 END IF END IF END IF END IF C C USE THE SEARCH MECHANISM TO FIND THE PROPER INTERVAL LENGTH. C 2140 IF (ENDPTS.LT.2) GO TO 2750 ENDPTS=3 SEARCH=2 2150 GO TO (2160,620,970,2165), JUMPS 2155 SEARCH=2 2160 INC=4 INC2=1 IF (JUMPS.NE.1) GO TO 990 IF (IC1+IC1+NSUB-12) 990,2690,2690 2165 IF (J1.LT.J2) GO TO 945 XJUMP=XJUMPS GO TO 870 C C ABSCISSAE HAVE COALESCED. C 2170 IF (K.LT.4) GO TO 100 2180 DELTA=100.0e0*DELTA IF (RE.LT.0.875e0) GO TO 2200 IF (ERRI.LT.MAX(EPSMIN,PEPSMN)) GO TO 2200 IF (ERR.LT.MAX(EPSMIN,PEPSMN)) GO TO 2200 IF (REPROD.LE.1.0E-4) GO TO 2200 IF (HAVDIF.AND.ABS(JUMPS-3).NE.1) GO TO 2200 2190 CONTINUE CALL SINTO (9,WORK) IFLAG(1)=5+KDIM WORK(1)=EMINF ANSWER=ALOCAL GO TO 2340 2200 CONTINUE IF (WHERE.LT.0) THEN C WHERE .LT. 0 MEANS ABSCISSAE COALESCED DURING THE QUADRATURE C STEP. ACUM IS GARBAGE DURING THE QUADRATURE STEP SO USE PACUM. ACUM=PACUM AACUM = abs(ACUM) EPSMIN=PEPSMN ELSE TPS=ACUM-PACUM ACUM=PACUM TPS=(RE*TPS)/(1.0e0-RE) IF (ABS(RE-REP).GE.0.125e0*RE) THEN C DO NOT TRUST THE ERROR ESTIMATE. ERR=ERR+ABS(TPS) ELSE C TRUST THE ERROR ESTIMATE, APPLY AITKEN ACCELERATION. ACUM=ACUM-TPS END IF END IF IF (IPRINT.GT.1) CALL SINTO (10,WORK) GO TO 2250 C C INTEGRATION WAS SUCCESSFUL ON LAST INTERVAL. CHECK IF THE C COMPLETE INTERVAL HAS BEEN INTEGRATED. C 2235 IF ((ERR*RNDC).GT.EPSMIN) DELTA=1000.0e0*DELTA GO TO 2240 2237 ERR=ERRC 2240 DISCHK=MAX(DISCHK-1,0) 2250 SEARCH=1 I=PART IF (PART.NE.1) THEN IF (STEP(1)*STEP(2).GT.0.0e0) I=1 END IF ERRT(I)=ERRT(I)+MAX(ERR,EPSMIN) c At one time, we had max(err,epsmin,edmeps*result(i)), where edmeps c is 0.5*r1mach(4), to account for round-off introduced during c accumulation. This effect should be negligible. XCDOBT(I)=XCDOBT(I).OR.ERR.GT.(EPSMIN/AACUM)**RELOBT*AACUM TPS=MAX(EPSO,ABS(RELEPS)) epsr=errt(1)+errt(2)+edue2a+edue2b EPSR=MAX(TPS-epsr,MAX(0.1e0*epsr,0.0e0)) FAIL=.FALSE. KMIN=0 RESULT(I)=RESULT(I)+ACUM IF (IPRINT.GT.1) CALL SINTO (11,WORK) IF (DISCF.gt.0) THEN IF (PART.EQ.1) THEN DISCF=0 ELSE IF (.NOT.IEND) GO TO 2350 END IF IF (IPRINT.GT.0) CALL SINTO (13,WORK) END IF IF (.NOT.IEND) GO TO 2350 IF (PART.NE.1) GO TO 2612 2320 ANSWER=RESULT(1) IF (BINIT.LT.AINIT) ANSWER=-ANSWER WORK(1)=ERRT(1)+EDUE2A+EDUE2B C DO BLOCK IF (WORK(1).GT.EPSO) THEN IFLAG(1)=-2 IF (XCDOBT(1)) IFLAG(1)=-3 IF (RELTOL.EQ.0) GO TO 2340 TPS=ABS(RELEPS) RELEPS=ABS(ANSWER)-WORK(1) RELEPS=MAX(ABS(WORK(RELTOL+2))*RELEPS,ENZER) IF (WORK(1).GT.RELEPS) THEN IF (WORK(1).GT.TPS) GO TO 2340 IF (IFLAG(1).GE.-2) GO TO 2340 GO TO 10 END IF END IF IFLAG(1)=-1 C END BLOCK 2340 WHERE=0 2345 IF (NFINDX.NE.0.AND.KDIM.EQ.1) IFLAG(NFINDX)=NFEVAL RETURN C C THE INTEGRATION IS COMPLETE THROUGH THE PANEL ENDING AT BLOCAL. C 2350 FATA=FATB FATAS=FATBS TLEN=TLEN-ABSDIF-ABSDIF IF (abs(DISCF).GE.2) THEN DISCF=sign(1,discf) ALOCAL=BLOCAL EPS=EPSR*(ABS(END(PART)-ALOCAL)/TLEN) IF (PART-1) 40,40,50 END IF KMAX=MAXK KK=K IF (DID1) THEN IF (K.LT.KAIMT) THEN IF (ERR.GT.ERSQE6*EPSMIN) KK=KK+1 END IF ELSE IF (K.LT.6) THEN KK=MIN(6,K+NSUB/2) END IF IF (NSUB.NE.0) THEN C DO BLOCK IF (ABS(BLOCAL-START(PART)).LT.ABS(TEND-START(PART))-DELMIN) 1 THEN IF (.NOT.HAVDIF) GO TO 2390 IF (PHISUM.LT.REAL(10*(K-KAIMT))*PHTSUM) GO TO 2380 IF (PHISUM.GE.16.0e0*PHTSUM) GO TO 2390 IF (PHISUM.GE.PHTSUM*0.5e0**(KAIMT-K+1)) GO TO 2390 GO TO 2380 END IF TEND=END(PART) KMAX=KMAXF C END BLOCK 2380 CONTINUE HAVDIF=.FALSE. CALL SINTNS (NSRB) END IF 2390 CONTINUE ALOCAL=BLOCAL ERRC=MAX(ERR,ERRI) EPNOIZ=MIN(ERRC,4.0e0*EPNOIZ) IF (.NOT.ROUNDF) THEN IF (REP.GT.0.3e0) KMIN=K-1 S=10.0e0*ERRI IF (S.LE.EPSMIN.AND.RNDC.GT.FER) RNDC=MAX(FER,S*RNDC/EPSMIN) END IF S=MAX(REP,1.0E-5) TP1=1.0e0 C DO BLOCK 2420 TP1=TP1*1.125e0 S=S*6.0e0 IF (S.LT.0.4e0) GO TO 2420 C END BLOCK C DO BLOCK IF (K.GE.MAX(KAIMT-1,3).AND.K.LE.AIMFOR.AND.ERRI.LE.EPSMIN.AND. 1 REP.LE.0.1e0) THEN I=AIMFOR-K+4 KK=K-1 ELSE I=AIMFOR-KK+3 C 1 .LE. I .LE. 7 IF (I.GE.3) THEN IF (I.EQ.3.AND.ERRI.GE.EPSMIN.AND.REP.GT.0.3e0)GO TO 2450 IF (ERR.EQ.0.0e0) GO TO 2450 IF (ERR.GT.EPSMIN) GO TO 2450 IF (ERRI.GT.10.0e0*EPSMIN) GO TO 2450 C ADJUST STEP SIZE WHEN APPROACHING NOISE OR ROUND-OFF C LIMIT. DELTA=DELTA*TP1 GO TO 2450 END IF END IF ERRC=MAX(ERRC,ABS(EP)) ERR=ERRC RE=REP IF (REP.GT.1.0e0) DELTA=0.875e0*DELTA C END BLOCK 2450 CONTINUE TPS=EPS IF (K.GE.KK) TPS=MIN(ESQEPS*EPSMIN/RNDC,EPS) ALPHA=0.0e0 IF (TPS.NE.0.0e0) ALPHA=MAX(ERRC*MIN(1.0e0,SQRT(RE)),EPSMIN)/TPS 2480 DELTA=DELTA*XSTEP(I) TP=ABS(ALOCAL-START(PART)) KAIMT=AIMFOR C DO BLOCK IF (SEARCH.NE.1.OR.TP.LT.ABS(WORRY(PART)-START(PART))-DELMIN) THEN C C CHOOSE THE NEXT STEPSIZE BASED UPON THE ORDER OF THE CURRENT C FORMULA (SOME FUNCTION OF K), AND THE RATIO OF THE ERROR C COMMITTED TO THE ERROR REQUESTED ON THE CURRENT PANEL. C THE ERROR BEHAVES AS H**ORDER(K-1) WHERE H IS THE STEPSIZE. C THEN SET GAMMA = 1.125**ORDER(K-1), AND CHANGE THE ERROR C RATIO BY GAMMA AND THE STEPSIZE BY 1.125 UNTIL THE ERROR C RATIO IS GREATER THAN 2.0. C IF (ALPHA.NE.1.0e0) THEN IF (ALPHA.LT.1.0e0) THEN IF (ALPHA.EQ.0.0e0) GO TO 2500 BETA=1.125e0 ALPHA=ALPHA*GAMMA(KK-1) ELSE ALPHA=1.0e0/ALPHA BETA=1.0e0/1.125e0 DELTA=DELTA*BETA END IF 2490 CONTINUE ALPHA=ALPHA*GAMMA(KK-1) IF (ALPHA.GT.2.0e0) GO TO 2500 DELTA=DELTA*BETA GO TO 2490 END IF 2500 CONTINUE IF (SEARCH.NE.1) THEN DISCHK=MAX(DISCHK,0) IF (IPRINT.GT.3) CALL SINTO (12,WORK) IF (SEARCH.EQ.6) GO TO 2560 KMAX=MIN(MAXK-1,KMAX) IF (DELTA.LE.Z(2)*ABSDIF 1 .OR.(K.LT.6.AND.DELTA.LT.Z(3)*ABSDIF)) GO TO 2610 IF (DELTA.GT.Z(4)*ABSDIF) KMAX=MIN(KMAX,5) KAIMT=MIN(KMAX,KAIMT) IF (K.GE.KMAX) GO TO 2560 DELTA=ABSDIF+ABSDIF ERR=ERRC+EXTRA IF (K.GE.4) DISCHK=-1 IF (ERR.GT.EPS) GO TO 100 IF (RE.GT.1.0e0) GO TO 100 IF (KONVRG-6) 100,100,2240 END IF IF (KK.LT.K) DELTA=MAX(DELTA,2.0e0*TP1*XSTEP(AIMFOR-K+3)*ABSDIF) IF (DELTA.GT.ABSDIF+ABSDIF) KMIN=K IF (TP+DELTA.LT.ABS(END(PART)-START(PART))) THEN IF (MIN(ABS(TB),TP+DELTA).GE.ABS(WORRY(PART)-START(PART))) THEN DISCHK=-1 KMAX=KMAXF IF (TP+1.25e0*DELTA.GE.ABS(TB)) GO TO 2540 END IF END IF C DO BLOCK IF (TP.LT.ABS(TEND-START(PART))-DELMIN) THEN IF (K.GE.KAIMT) GO TO 2560 IF (.NOT.HAVDIF) GO TO 2560 IF (PHISUM.LE.10.0e0*PHTSUM) GO TO 2560 IF (NSUB.GE.NSUBMX) GO TO 2560 BLOCAL=ALOCAL+SIGN(DELTA,TB) CALL SINTNS (NSIB) GO TO 2620 END IF TEND=ALOCAL+SIGN(DELTA,TB) KMAX=KMAXF C END BLOCK IF (DID1.OR.PART.EQ.1.OR.FAIL) GO TO 2610 C AFTER THE FIRST SUCCESS IN PART 2, TREAT THE NEXT INTERVAL C LIKE A RESTART. ALSO, INCREASE THE STARTING STEPSIZE FOR C RESTARTING PART 1 IF POSSIBLE. KAIMT=3 KMAX=KMAXF X1=ALOCAL+SIGN(DELTA,TB) X2=ALOCAL IF (NSUB.NE.0) THEN X1=TDECR(X1) X2=TDECR(X2) IF (NSUB.NE.2) THEN X1=TDECR(X1) X2=TDECR(X2) END IF END IF STEP(1)=SIGN(MAX(ABS(STEP(1)),ABS(X1-X2)),STEP(1)) GO TO 2610 END IF 2540 continue DISCHK=-1 KMAX=KMAXF DELTA=MAX(DELTA,ABS(END(PART)-ALOCAL)) IF (DELTA.GT.3.0e0*ABSDIF) DISCHK=-2 WORRY(PART)=END(PART) GO TO 2610 C END BLOCK 2560 CONTINUE DELTA=MIN(ABSDIF,DELTA) WORRY(PART)=BLOCAL KMAX=KMAXF KAIMT=MIN(MAX(5,KMAXF),AIMFOR) 2610 FATBS=.FALSE. GO TO 2620 2612 IEND=.FALSE. PART=1 NSUB=0 DISCF=0 DISCHK=MAX(DISCHK,1) RESULT(1)=RESULT(1)+RESULT(2) ERRT(1)=ERRT(1)+ERRT(2) ERRT(2)=0.0e0 XCDOBT(1)=XCDOBT(1).OR.XCDOBT(2) ALOCAL=START(1) TA=TASAVE TB=END(1)-TA TLEN=ABS(END(1)-START(1)) DELMIN=SINTSM(ALOCAL) 2614 IF (NSUB.LT.NSUBSV) THEN CALL SINTNS (NSIA) GO TO 2614 END IF IF (DID1) DELTA=0.0e0 DELTA=MAX(ABS(STEP(1)),MAX(10.0e0*DELMIN,DELTA)) STEP(1)=SIGN(1.0e0,SIGN(1.0e0,TB)*(BINIT-AINIT)) FATBS=.FALSE. FATAS=FSAVED FATA=FSAVE FSAVED=.FALSE. DID1=.FALSE. 2615 KAIMT=3 KMAX=KMAXF 2620 X2=ABS(END(PART)-ALOCAL) C X2 IS THE LENGTH OF THE REST OF THE PART. INIT=.FALSE. EPS=EPSR*X2/TLEN DELMIN=SINTSM(ALOCAL) IF (DELTA.LT.X2) THEN IEND=.FALSE. IF (DELTA.LT.DELMIN) THEN DELTA=DELMIN FATBS=.FALSE. END IF C DO NOT LEAVE A PIECE TOO SMALL. IF (X2-DELTA.LE.DELMIN) GO TO 60 BLOCAL=ALOCAL+SIGN(DELTA,TB) GO TO 75 END IF C REDUCE KAIMT IF DELTA IS VERY MUCH GREATER THAN THE REMAINING C PART OF THE MAJOR SUBDIVISION. C DO FOREVER 2640 CONTINUE X2=3.0e0*X2 IF (X2.GT.DELTA) GO TO 60 KAIMT=KAIMT-1 IF (KAIMT.LT.3) GO TO 60 GO TO 2640 C END FOREVER C C ***** INTERNAL SUBROUTINES ******************************* C C COMPUTE MINIMUM RATIO OF STEPSIZE CHANGES. IF THE STEPSIZE C CHANGES MORE RAPIDLY THAN 0.3, ADD A POINT IN THE LARGEST C INTERVAL. C 2670 WHERE=3 KK=MAX(INC2+1,3) J=MAX(INC,3) EPS=S GO TO 2710 2680 WHERE=2 C ABSCISSAE HAVE COALESCED. J=LENDT GO TO 2700 2690 WHERE=1 J=LENDT IF (JUMPS.NE.1) GO TO 2700 IF (SEARCH.NE.2.AND.SEARCH.NE.6) 1 IF (KONVRG+KONVRG-LENDT) 2700,2700,1190 KK=4 J=6 GO TO 2710 2700 KK=3 2710 TP=0.3e0 TP1=ABS(XT(KK-1)-XT(KK-2)) DO 2720 I=KK,J S=ABS(XT(I)-XT(I-1)) IF (TP*S.GE.TP1) THEN TP=TP1/S L=I ELSE IF (TP*TP1.GT.S) THEN TP=S/TP1 L=I-1 END IF TP1=S 2720 CONTINUE IF (TP.EQ.0.3e0) GO TO (910,1000,1160), WHERE ABSCIS=XT(L-1)+0.5e0*(XT(L)-XT(L-1)) IF (WHERE.NE.2) THEN IF (ABSCIS.EQ.XT(L-1)) GO TO 1000 IF (ABSCIS.EQ.XT(L)) GO TO 1000 END IF IF (JUMPS.EQ.1) GO TO 1520 IF (J1.EQ.18) J1=J2+1 IF (J2.EQ.0) J2=J1-1 GO TO 1520 C C REQUEST A FUNCTION VALUE AT THE EDGE OF A PANEL. C 2730 IF (DISCHK.LE.0) GO TO 100 2740 IF (ENDPTS-2) 2750,2760,100 2750 IF (J1.EQ.LENDT) GO TO 2760 WHERE=5 WHERE2=1 ENDPTS=2 IF (J2.EQ.1) ENDPTS=3 GO TO 2810 2760 WHERE=5 WHERE2=2 ENDPTS=3 IF (J2.EQ.1) GO TO 530 2810 ABSCIS=LOCAL(WHERE2) X=ABSCIS IF (FATS(WHERE2)) GO TO 2830 GO TO 2850 2820 FAT(WHERE2)=FNCVAL FATS(WHERE2)=.TRUE. ERRAT(WHERE2)=FER*ABS(FNCVAL) IF (FEA.NE.0) ERRAT(WHERE2)=ERRAT(WHERE2)+ABS(T*WORK(FEA)) IF (LOCAL(WHERE2).EQ.AINIT) THEN EDUE2A=FNCVAL*ERRINA EPS=MAX(EPSR-(EDUE2A+EDUE2B),EDUE2A+EDUE2B)*DELTA/TLEN ELSE IF (LOCAL(WHERE2).EQ.BINIT) THEN EDUE2B=FNCVAL*ERRINB EPS=MAX(EPSR-(EDUE2A+EDUE2B),EDUE2A+EDUE2B)*DELTA/TLEN END IF EPS=MAX(EPS,(EPSMIN/AACUM)**RELOBT*AACUM) 2830 J=WHERE2+2*(WHERE-4) GO TO (280,290,350,350), J C C COMPUTE F USING FORWARD OR REVERSE COMMUNICATION. C 2840 IF (ABS(ABSCIS-LOCAL(4-WHERE)).GT.ESMALL) GO TO 2850 WHERE=-WHERE GO TO 2180 2850 IF (NFMAX.LE.0) GO TO 2860 IF (NFEVAL.LT.NFMAX) GO TO 2860 WHERE=7 IFLAG(1)=5 GO TO 2345 2860 NFEVAL=NFEVAL+1 T=1.0e0 IF (NSUB.NE.0) THEN c Put ABSCIS in user coordinates X1=TA/TB T=(ABSCIS-TA)/TB c? ABSCIS=TA*(1.0e0+X1)+ABSCIS*(T-X1) ABSCIS=TA+ABSCIS*T+X1*(TA-ABSCIS) IF (NSUB.NE.2) THEN F2=(ABSCIS-TA)/TB c? ABSCIS=TA*(1.0e0+X1)+ABSCIS*(F2-X1) ABSCIS=TA+ABSCIS*F2+X1*(TA-ABSCIS) T=T*F2 T=T+T END IF T=T+T END IF IF (ABSCIS.EQ.ALOCAL) THEN c## F2=SINTSM(ABSCIS)/EDELM2 c## ABSCIS=ABSCIS+SIGN(F2,TB) abscis=abscis+sign(max(enzer,emeps*abs(abscis)),tb) X=ABSCIS ELSE IF (ABSCIS.EQ.BLOCAL) THEN c## F2=SINTSM(ABSCIS)/EDELM2 c## ABSCIS=ABSCIS-SIGN(F2,TB) abscis=abscis-sign(max(enzer,emeps*abs(abscis)),tb) X=ABSCIS END IF WORK(KDIM)=ABSCIS IF (REVERS.NE.0) RETURN CALL SINTF (ANSWER,WORK,IFLAG(1)) 2900 FNCVAL=T*ANSWER GO TO (90,150,160,2820,2820,350,2860), WHERE C END .