C	WESTERN MICHIGAN UNIVERSITY
C	ITGTN.FOR (FILENAME ON LIBRARY DECTAPE)
C	ITGTN.FOR IS CALLED (BY RUNUUO FROM INTGR.FOR
C	EXTERNAL FUNCTIONS USED (GENERATED BY INTGR.FOR): SEE
C	 COMMENTS AT BEGINNING OF SOURCE LISTING OF INTGR.FOR
C	INTERNAL FUNCTIONS USED:  QMULT2, QMULT3,SQUANK
C	INTERNAL SUBR. USED:  INT2, INT3, GLO16
C	ABOVE COMMENTS AND RIGHT ADJUSTED COMMENTS PUT IN BY WG
C
C	THIS IS THE SECOND PART OF THE INTEGRATION PROGRAM IN WHICH
C	THE FOLLOWING  VARIABLES HAVE TO BE READ IN FROM THE TEMPORARY
C	FILE:
C
C	     IDIM---DIMENSION
C	     FF-----THE FUNCTION TO BE INTEGRATED
C	     FFL1---LOWER LIMIT OF THE INNER INTEGRAL
C	     FFU1---UPPER LIMIT OF THE INNER INTEGRAL
C	     FFL2---LOWER LIMIT OF THE MIDDLE INTEGRAL IN THE 3-
C	            DIMENSION INTEGRATION
C	     FFU2---UPPER LIMIT OF THE MIDDLE INTEGRAL
C
C
C	THE FOLLOWING ARE FUNCTIONS AND/OR LIMITS TO BE USED IN THE
C	INTEGRATION:
C
C	     F----FUNCTION TO BE INTEGRATED
C	     FL1--LOWER LIMIT OF THE INNER INTEGRAL
C	     FU1--UPPER LIMIT OF THE INNER INTEGRAL
C	     FL2--LOWER LIMIT OF THE MIDDLE INTEGRAL IN THE 3-DIMENSION
C	          INTEGRATION
C	     FU2--UPPER LIMIT OF THE MIDDLE INTEGRAL
C	     FL3--LOWER LIMIT OF THE OUTER INTEGRAL
C	     FU3--UPPER LIMIT OF THEOUTER INTEGRAL
C
C***********************************************************************
C
	DIMENSION FF(12),FFU1(12),FFL1(12),FFU2(12),FFL2(12)
	EXTERNAL F,FU1,FL1,FU2,FL2
C
C***********************************************************************
C	DEVICES USED:  SEE PART 1
C***********************************************************************
C
	IDLG=-1
	ICC=-4
	IOUT=30
	IDSK=1
	ITEMP=20
C
C***********************************************************************
C	READ IN INFORMATION THAT WERE TEMPORARILY STORED ON THE
C	DISK BY THE FIRST PART OF THE PROGRAM.
C***********************************************************************
C
	READ(ITEMP) IDIM,FF
	READ(ITEMP) FFL1,FFU1
	READ(ITEMP) FFL2,FFU2
	CALL RELEAS (ITEMP)
C
C***********************************************************************
C	ASK AND OBTAIN OUTER LIMITS--MUST BE REAL NUMBERS
C**********************************************************************
C
	IF (IDIM.EQ.1) WRITE(IDLG,61)
	IF (IDIM.GT.1) WRITE(IDLG,580)
580	FORMAT('-ENTER OUTER LIMITS:')
	WRITE(IDLG,101)
101	FORMAT(' LOWER: ',$)
	READ(ICC,53) FL3
	WRITE(IDLG,102)
102	FORMAT('+UPPER: ',$)
	READ(ICC,53) FU3
	IF (IDIM.GT.1) GO TO 540
C
C***********************************************************************
C	ONE-DIMENSION  INTEGRATION
C***********************************************************************
C
	WRITE(IDLG,541)
541	FORMAT(' ENTER VALUE FOR THE REQUIRED TOLERANCE: ',$)
	READ(ICC,53) ERROR
53	FORMAT(F)
C---------------SQUANK MAKES USE OF FUNCTION F GENERATED BY INTGR.FOR
	CALL SQUANK(FL3,FU3,ERROR,FIFTH,RUM,NO)
	S=SQUANK(FL3,FU3,ERROR,FIFTH,RUM,NO)
C
C***********************************************************************
C	START OF THE REPORT
C***********************************************************************
C
540	DO 5410 J=12,1,-1
	IF (FF(J).NE.' ') GO TO 542
5410	CONTINUE
544	WRITE(IDLG,543) FF
543	FORMAT('- ERROR IN FUNCTION, JOB TERMINATES'/1X,12A5)
	CALL EXIT
542	WRITE(IOUT,54) IDIM,(FF(JJ),JJ=1,J)
54	FORMAT('1',I1,'-DIMENSION INTEGRATION'//' FUNCTION: ',12A5)
	IF (IDIM.LT.2) GO TO 60
C
C***********************************************************************
C	TWO AND THREE-DIMENSION  INTEGRATION
C***********************************************************************
C
	WRITE(IOUT,30)
30	FORMAT('-LIMITS FOR INNER INTEGRAL:')
	DO 550 J=12,1,-1
	IF (FFL1(J).NE.' ') GO TO 551
550	CONTINUE
	GO TO 544
551	WRITE(IOUT,55) (FFL1(JJ),JJ=1,J)
55	FORMAT(' LOWER: ',12A5)
	DO 560 J=12,1,-1
	IF (FFU1(J).NE.' ') GO TO 561
560	CONTINUE
	GO TO 544
561	WRITE(IOUT,56) (FFU1(JJ),JJ=1,J)
56	FORMAT(' UPPER: ',12A5)
	IF (IDIM.LT.3) GO TO 57
	WRITE(IOUT,43)
43	FORMAT('-LIMITS FOR MIDDLE INTEGRAL:')
	DO 5500 J=12,1,-1
	IF (FFL2(J).NE.' ') GO TO 5501
5500	CONTINUE
	GO TO 544
5501	WRITE(IOUT,55) (FFL2(JJ),JJ=1,J)
	DO 5503 J=12,1,-1
	IF (FFU2(J).NE.' ') GO TO 5504
5503	CONTINUE
	GO TO 544
5504	WRITE(IOUT,56) (FFU2(JJ),JJ=1,J)
57	WRITE(IOUT,58)
58	FORMAT('-OUTER LIMITS:')
	WRITE(IOUT,5510) FL3
	WRITE(IOUT,5610) FU3
	IF (IDIM.EQ.2) CALL INT2(FU3,FL3,IOUT)
	IF (IDIM.EQ.3) CALL INT3(FU3,FL3,IOUT)
	GO TO 70
C
C***********************************************************************
C	CONTINUATION OF REPORT FOR ONE-DIMENSION INTEGRATION
C***********************************************************************
C
60	WRITE(IOUT,61)
61	FORMAT('-LIMITS OF INTEGRATION')
	WRITE(IOUT,5510) FL3
5510	FORMAT(' LOWER: ',G15.7)
	WRITE(IOUT,5610) FU3
5610	FORMAT(' UPPER: ', G15.7)
	WRITE(IOUT,65) S,FIFTH,NO,RUM
65	FORMAT('-ANSWER: ',G15.7/'-FIFTH-ORDER ADJUSTMENT TERM ',G15.7/
     1 ' NUMBER OF FUNCTION EVALUATION REQUIRED: ',G15.7/' THE CLAIMED
     1 TOLERANCE (ADJUSTED FOR ROUNDOFF ERROR): ',G15.7)
C
C***********************************************************************
C	END OF ONE INTEGRATION, EITHER CALL EXIT OR BRANCH TO FIRST
C	PART OF THE PROGRAM
C***********************************************************************
C
70	WRITE(IDLG,71)
71	FORMAT('-END OF INTEGRATION, TYPE 1 TO CONTINUE OR 0 TO
     1 TERMINATE'/)
	READ(ICC,14) LAST
14	FORMAT(I)
C---------------INTGR.EXE IS IN SYS:
	IF (LAST.GT.0) CALL RUNUUO('R INTGR')
	CALL EXIT
	END
C---------------ALL ARGS. ARE INPUT.
	SUBROUTINE INT2(FU3,FL3,IOUT)
C	THIS PROGRAM APPROXIMATES A 2-DIMENSIONAL ITERATED
C	INTEGRAL OF F(X,Y) AS DISCUSSED IN SECTION 1.6.
C	SUBROUTINE GLO16 PROVIDES A TABLE OF THE 16-POINT
C	GAUSS-LEGENDRE FORMULA.  THE LIMITS ON THE VARIABLES ARE
C	         AA .LE. X .LE. BB
C	     FL1(X) .LE. Y. LE. FU1(X)
C
	DOUBLE PRECISION DX(16),DA(16)
	DIMENSION X(40,2),A(40,2),MM(2)
	EXTERNAL F,FL1,FU1
	CALL GLO16(DX,DA,-1.D0,1.D0)
	DO 2 I = 1,16
	XX=DX(I)
	AB=DA(I)
	X(I,1)=XX
	X(I,2)=XX
	A(I,1)=AB
2	A(I,2)=AB
	MM(1)=16
	MM(2)=16
	Q = QMULT2(F,FL3,FU3,FU1,FL1,X,A,MM)
	WRITE(IOUT,4) Q
4	FORMAT('-ANSWER: ',G15.7)
	RETURN
	END
C---------------ALL ARGS. ARE INPUT
	FUNCTION QMULT2(FCN,AA,BB,FU,FL,X,A,MM)
	DIMENSION X(40,2),A(40,2),MM(2)
	H1=(BB-AA)/2.
	G1=(BB+AA)/2.
	Q1=0.
	M1=MM(1)
	M2=MM(2)
	DO 4 I=1,M1
	UI=H1*X(I,1)+G1
	AI=H1*A(I,1)
	D=FU(UI)
	C=FL(UI)
	H=(D-C)/2.
	G=(D+C)/2.
	Q=0.
	DO 2 J=1,M2
	VJ=H*X(J,2)+G
C---------------THERE IS NO FUNCTION FCN;  HOWEVER THE EXTERNAL F
C--------------- IN SUBR. INT2 AND THE F IN ARG. LIST OF QMULT2
C--------------- IN SUBR. INT2 ALLOWS FOR SUBST. OF F FOR FCN
2	Q=Q+A(J,2)*FCN(UI,VJ)
4	Q1=Q1+AI*H*Q
	QMULT2 = Q1
	RETURN
	END
C---------------ALL ARGS. ARE INPUT.
	SUBROUTINE INT3(FU3,FL3,IOUT)
C
C	THIS PROGRAM APPROXIMATES A 3-DIMENSIONAL ITERATED
C	INTEGRAL OF F(X,Y,Z) AS DISCUSSED IN SECTION 1.6.
C	SUBROUTINE GLO16 PROVIDES A TABLE OF THE 16-POINT
C	GAUSS-LEGENDRE FORMULA.  THE LIMITS ON THE VARIABLES ARE
C	            AA .LE. X .LE. BB
C	        FL2(X) .LE. Y .LE. FU2(X)
C	      FL1(X,Y) .LE. Z .LE. FU1(X,Y)
C
	DOUBLE PRECISION DX(16), DA(16)
	DIMENSION X(40,3),A(40,3), MM(3)
	EXTERNAL F,FL1,FU1,FL2,FU2
	CALL GLO16(DX,DA,-1.D0,1.D0)
	DO 2 I=1,16
	XX=DX(I)
	X(I,1)=XX
	X(I,2)=XX
	X(I,3)=XX
	AB=DA(I)
	A(I,1)=AB
	A(I,2)=AB
2	A(I,3)=AB
	MM(1)=16
	MM(2)=16
	MM(3)=16
	Q = QMULT3(F, FL3, FU3, FU2, FL2, FU1, FL1, X, A, MM)
	WRITE(IOUT,4) Q
4	FORMAT('-ANSWER: ',G15.7)
	RETURN
	END
C---------------ALL ARGS. ARE INPUT.
	FUNCTION QMULT3(FCN,AA,BB,FU1,FL1,FU2,FL2,X,A,MM)
	DIMENSION X(40,3), A(40,3), MM(3)
	H1= (BB-AA)/2.
	G1= (BB+AA)/2.
	Q1= 0.
	M1= MM(1)
	M2= MM(2)
	M3= MM(3)
	DO 6 I = 1,M1
	UI = H1*X(I,1) + G1
	AI = H1*A(I,1) 
	D1 = FU1(UI)
	C1 = FL1(UI)
	H2 = (D1-C1)/2.
	G2 = (D1+C1)/2.
	Q2 = 0.
	DO 4 J=1,M2
	VJ = H2*X(J,2) + G2
	AJ = H2*A(J,2)
	D = FU2(UI,VJ)
	C = FL2(UI,VJ)
	H = (D-C)/2.
	G = (D+C)/2.
	Q = 0.
	DO 2 K=1,M3
	WK = H*X(K,3)+G
2	Q = Q + A(K,3)*FCN(UI,VJ,WK)
4	Q2= Q2+ AJ*H*Q
6	Q1= Q1 + AI*Q2
	QMULT3=Q1
	RETURN
	END
C---------------A, BIG, ERROR ARE INPUT.  FIFTH, RUM, NO ARE RETURNED.
	FUNCTION SQUANK(A,BIG,ERROR,FIFTH,RUM,NO)
C
C	S*Q*U*A*N*K  STANDS FOR * SIMPSON QUADRATURE USED
C	       ADAPTIVELY.  NOISE KILLED.*
C
C
C
C
C
C		CALLING PROGRAM REQUIRES
C	EXTERNAL F
C	THIS IS FUNCTION TO BE INTEGRATED
C	A	THE LOWER LIMIT OF INTEGRATION
C	BIG	THE UPPER LIMIT
C	ERROR	THE REQUIRED TOLERANCE (ABSOLUTE ERROR)
C
C
C
C
C	OUTPUT
C
C	SQUANK	THE FIFTH ORDER RESULT = THIRD + FIFTH
C	FIFTH	THE FIFTH ORDER ADJUSTMENT TERM
C	RUM	THE CLAIMED TOLERANCE (ADJUSTED FOR ROUNDOFF ERROR)
C	NO	THE NUMBER OF FUNCTION EVALUATIONS REQUIRED
C
C
C
C
C	NOTES ON USE.  (1)  DISCONTINUOUS FUNCTIONS
C
C	THIS ROUTINE IS BASED ON DEGREE 3 AND DEGREE 5 LOCAL
C	POLYNOMIAL APPROXIMATION.  CONSEQUENTLY IT SHOULD NOT BE USED
C	WITH FUNCTIONS WHICH HAVE DISCONTINUITIES IN THE FOURTH OR LOWER
C	DERIVATIVES WITHIN THE INTERVAL OF INTEGRATION.  IF THERE ARE
C	SUCH DISCONTINUITIES, THIS ROUTINE WILL TAKE THIS TO BE EVIDENCE
C	OF ROUND OFF ERROR IN FUNCTION VALUES AND WILL ADJUST THE
C	TOLERANCE.
C
C	IF THE LOCATIONS OF SUCH DISCONTINUITIES ARE KNOWN, THE
C	ROUTINE MAY BE USED SEPARATELY FOR EACH INTERVAL BETWEEN
C	CONSECUTIVE DISCONTINUITIES.  WHILE, LIKE ALL SUCH ROUTINES, IT
C	DISLIKES DISCONTINUITIES, IT CAN HANDLE THEM IF THEY ARE LOCATED
C	AT THE END POINTS OF THE INTEGRATION INTERVAL.
C
C
C
C
C	NOTES ON USE.  (2)  FUNCTIONS WITH HIGH-FREQUENCY OSCILLATIONS.
C
C	THE ROUTINE WILL RETURN UNRELIABLE RESULTS FOR FUNCTIONS LIKE
C	G(X) TIMES COS(100*X).  IF THE HIGHEST PERIOD LIKELY TO BE
C	ENCOUNTERED IS KNOWN, THE INTERVAL SHOULD BE SUB-DIVIDED IN
C	SUCH A WAY THAT,
C
C	(ONE)  THERE ARE NOT MORE THAN THREE PERIODS PER INTERVAL, AND
C	(TWO)  THE PERIOD *P* DIVIDED BY SUB-INTERVAL, *(B-A)* IS NOT A
C	       SIMPLE FRACTION *N/M* WITH N OR M LESS THAN 9.
C
C
C
C
C	NOTES ON USE.  (3)  INTERVAL SUB-DIVISION
C
C	THE FAILURES DESCRIBED ABOVE ARE GENERALLY WORSE FOR *SQUANK*
C	THAN FOR OTHER ROUTINES BECAUSE *SQUANK* TAKES THE INCONVENIENT
C	BEHAVIOR AS AN INDICATION OF ROUND OFF ERROR.  IN GENERAL SUB-
C	DIVISION OF THE INTERVAL IS ADVOCATED.  ESSENTIALLY THE USER
C	CARRIES OUT, UNDER DRIVING PROGRAM CONTROL, A SEQUENCE OF
C	CALCULATIONS WHICH SHOULD HAVE BEEN CARRIED OUT IN THE
C	SUBROUTINE IN ANY CASE.  IN THIS WAY HE PREVENTS CHANCE LOW
C	ORDER FALSE CONVERGENCE AT VIRTUALLY NO ADDITIONAL COST.  NOTE
C	THAT THE SUM OF THE PARAMETERS *ERROR* FOR THE SUB-INTERVALS
C	SHOULD CORRESPOND TO THE VALUE REQUIRED FOR THE WHOLE INTERVAL.
C
C
C
C
C	NIM NUMBERING SYSTEM AND LOGIC
C
C	THE INTERVAL (A,B) IS DEFINED NIM = 1, LEVEL = 0.
C
C	THE INTERVAL NIM = N, LEVEL = L IS BISECTED, IF NECESSARY,
C		INTO TWO INTERVALS, NIM = 2*N AND NIM = 2*N+1,
C		BOTH AT LEVEL = L+1.
C
C	IF INTERVAL NIM = N, LEVEL = L DOES NOT CONVERGE, THE NEXT
C		INTERVAL CONSIDERED IS NIM = 2*N, LEVEL = L+1.
C
C	IF INTERVAL NIM=N, LEVEL = L DOES CONVERGE, THE NEXT INTERVAL
C		CONSIDERED IS NIM = M(R) + 1, LEVEL = L-R, WHERE M(R)
C		IS THE FIRST EVEN MEMBER OF THE SEQUENCE M(0) = N,
C		M(S+1) =(M(S)-1)/2.  IF THIS GIVES LEVEL = 0,
C		THE CALCULATION IS COMPLETE.
C
C
C
C	SCALING TO AVOIDING EXCESSIVE DIVISION BY TWO.
C
C
	THE INTERVAL (X1,X5) IS OF LENGTH H = X5-X1.  THE POINTS
C	X1,X2,X3,X4,X5 ARE THE POINTS OF QUARTERSECTION OF THIS INTERVAL
C	AND FX1,FX2,FX3,FX4,FX5 ARE THE CORRESPONDING FUNCTION VALUES.
C
C	EST IS APPROXIMATION TO (6.0/H)* INTEGRAL(X1,X5).
C
C	(EST1+EST2) IS APPROXIMATION TO (12.0/H)*INTEGRAL(X1,X5).
C
C	SUM IS APPROXIMATION TO (12.0) * INTEGRAL(A,X1).
C
C
C
C
C	STORAGE
C
C	X3ST(L) =0.5*(X5ST(L)+X1).  THUS X3ST(L) COULD BE RECALCULATED
C		AT EACH STAGE TO AVOID STORAGE.  ESTST(L) IS SAME
C		IN THIS RESPECT.  THE RESULTS OF ABOVE RECALCULATION
C		ARE IDENTICAL MACHINE NUMBERS.
C
C	X5ST(L) = X1 + (B-A)*(2**(-L)).  THIS COULD ALSO BE RECALCU-
C		LATED.  BUT IN THIS CASE CALCULATION IS EXCESSIVE
C		AND THERE IS A POSSIBILITY OF ROUND OFF ERROR ARISING
C		BECAUSE THE SAME POINT IS BEING CALCULATED IN TWO OR
C		MORE DIFFERENT WAYS.
C
C
C
C
C	AVOIDANCE OF ROUND OFF ERROR TROUBLE
C
C	IF INTERVAL DOES NOT CONVERGE, FOLLOWING INTERVAL SHOULD HAVE
C	ADIFF VALUE APPROXIMATELY EQUAL TO (1/16) TIMES PREVIOUS ADIFF
C	VALUES, CALLED ADIFF1 IN THE CODE.  THERE IS A THOREM WHICH
C	STATES THAT, UNLESS THE FOURTH DERIVATION OF F(X) VANISHES
C	IN THE PREVIOUS INTERVAL, ADIFF IS LESS THAN OR EQUAL TO ADIFF1.
C	IF THIS DOES NOT HAPPEN, IT IS TAKEN TO BE AN INDICATION OF
C	POSSIBLE ROUND OFF LEVEL.  IN THIS CASE, UNLESS LEV IS LESS
C	THAN FIVE, THE CURRENT TOLERANCE LEVEL, CEPS, IS APPROPRIATELY
C	ADJUSTED.  HOWEVER CEPS IS RESET AS AND WHEN  APPEARS THAT IT
C	SHOULD BE ADJUSTED EITHER UP OR DOWN.  IT IS REDUCED IF
C	CONVERGENCE OCCURS WITH A NON-ZERO ADIFF STRICTLY LESS THAN
C	0.25*CEPS.  AN INVOLVED  SECTION OF CODING GUARDS TO SOME
C	EXTENT AGAINST AN UNREALISTIC VALUE ARISING AS A RESULT OF 
C	A ZERO IN THE FOURTH DERIVATIVE.  A FACTOR EFACT IS CALCULATED
C	WHICH ADJUSTS THE CLAIMED TOLERANCE TO TAKE INTO ACCOUNT
C	THESE ALTERATIONS IN THE TOLERANCE LEVEL.
C
C	THE ROUTINE ENTERS THESE INVOLVED SECTIONS OF CODING ONLY IF
C	ROUND OFF ERROR APPEARS TO BE PRESENT.  IN A NORMAL (ROUND OFF
C	ERROR FREE)  RUN, THESE SECTIONS ARE SKIPPED AT A COST OF A
C	SINGLE COMPARISION PER ITERATION (TWO FUNCTION EVALUATIONS).
C
C
C	ARBITRARY CONSTANTS
C
C	THE FOLLOWING CONSTANTS HAVE BEEN ASSIGNED IN THE LIGHT OF
C	EXPERIENCE WITH NO THEORETICAL JUSTIFICATON.
C
C	(1)  NO CONVERGENCE IS ALLOWED AT LEVEL = **ZERO**.  THIS MEANS
C	     THAT THE ROUTINE IS CONSTRAINED TO BASE THE RESULT ON AT
C	     LEAST 9 FUNCTION VALUES.
C
C	(2)  NO UPWARD ADJUSTMENT OF THE TOLERANCE LEVEL IS CONSIDERED
C	     AT LEVELS LOWER THAN LEVEL = **FIVE**.  THE POINT SPACING
C	     IS THEN (BIG-A)/128.0.
C
C	(3)  PHYSICAL LIMIT.  HIGHEST LEVEL ALLOWED IS LEVEL =
C	     **THIRTY**.  HERE CONVERGENCE IS ASSIGNED WHETHER OR NOT
C	     THE INTERVAL HAS CONVERGED.  THE POINT SPACING IS THEN
C	     ABOUT (BIG-A)*2.0*10**-10.
C
C	(4)  UPWARD ADJUSTMENT OF TOLERANCE LEVEL IS LIMITED IN GENERAL
C	     TO A FACTOR **2.0** OR LESS.
C
C	(5)  DOWNWARD ADJUSTMENT OF TOLERANCE LEVEL IS INHIBITED IN
C	     GENERAL UNLESS BY A FACTOR GREATER THAN **4.0**.
C
C
C
C
C	SOME NOTATION
C
C	SUM AND SIM ARE RUNNING SUMS, INCREASED AT STAGE EIGHT.  THEY 
C	ARE RESPECTIVELY 12.0 * (THIRD ORDER APPROXIMATION TO THE
C	INTEGRAL) AND -180.0* (FIFTH ORDER ADJUSTMENT TO THE INTEGRAL).
C
C	CEPSF IS THE REQUIRED (SCALED) TOLERANCE.
C
C	CEPS IS THE RUNNING VALUE OF THE ADJUSTED TOLERANCE.
C
C	QCEPS = 0.25 * CEPS
C
C	LEVTAG = -1 OR 0,2,3 INDICATES WHETHER TOLERANCE IS NOT OR IS
C	CURRENTLY ADJUSTED.  (SEE COMMENT IN STAGE SEVEN).
C
C	EFACT IS RUNNING SUM CORRESPONDING TO 180.0 * RUM
C
C	FACERR = 15.0 OR 1.0 DEPENDING ON WHETHER TOLERANCE IS OR
C	IS NOT CURRENTLY ADJUSTED.  IF IT IS, THERE IS NO JUSTIFICATION
C	FOR THE DIFFERENCE OF APPROXIMATIONS.  FACERR = 15.0 REMOVES
C	THE BUILT IN 15.0 FACTOR FOR CALCULATION OF EFACT.
C
C	EPMACH  THE MACHINE ACCURACY PARAMETER.  THE ROUND OFF ERROR
C	GUARD DOES NOT REQUIRE THIS NUMBER.  IT IS MACHINE INDEPENDENT.
C	THIS IS ONLY USED TO HELP IN AN INITIAL GUESS IN STAGE TWO
C	IF THE VALUE OF ERROR HAPPENS TO BE ZERO.  ANY NON-ZERO NUMBER
C	MAY BE USED INSTEAD, WITH A VERY SMALL PENALTY IN NUMBER OF
C	FUNCTION EVALUATIONS IF A COMPLETELY UNREASONABLE NUMBER IS
C	 USED.
C
C
C
	DIMENSION FX3ST(30),X3ST(30),ESTST(30),FX5ST(30),X5ST(30)
	DIMENSION PREDIF(30)
C---------------EXTERNAL F IS APPARENTLY NOT NECESSARY
	EXTERNAL F
	DOUBLE PRECISION SUM, SIM
	EPMACH = 0.0000000000075
C
C	     ****  STAGE ONE   ****
C	INITIALISE ALL QUANTITIES REQUIRED FOR CENTRAL CALCULATION
C	(STAGE 3).
C
	SUM=0.
	SIM=0.
	CEPSF=180.0*ERROR/(BIG-A)
	CEPS=CEPSF
	ADIFF = 0.0
	LEVTAG = -1
	FACERR = 1.0
	XZERO = A
	EFACT= 0.0
	NIM=1
	LEV=0
C
C	FIRST INTERVAL
C
	X1=A
	X5=BIG
	X3=0.5*(A+BIG)
	FX1=F(X1)
	FX3=F(X3)
	FX5=F(X5)
	NO=3
	EST= FX1 + FX5 + 4.0*FX3
C
C	     ****  STAGE TWO  ****
C	SET A STARTING VALUE FOR TOLERANCE IN CASE THAT CEPSF = 0.0
C
	IF (CEPSF) 295,205,295
205	LEVTAG = 0
	FACERR = 15.0
	CEPS = EPMACH*ABS(FX1)
	IF (FX1) 295,210,295
210	CEPS = EPMACH*ABS(FX3)
	LEVTAG = 3
	IF (FX3) 295,215,295
215	CEPS =EPMACH*ABS(FX5)
	IF (FX5) 295,220,295
220	CEPS = EPMACH
295	QCEPS = 0.25*CEPS
C
C	INITIALISING COMPLETE
C
C	     ****  STAGE THREE  ****
C	CENTRAL CALCULATION.
C		REQUIRES X1,X3,X5,FX1,FX3,FX5,EST,ADIFF.
C
300	CONTINUE
	X2=0.5*(X1+X3)
	X4=0.5*(X3+X5)
	FX2=F(X2)
	FX4=F(X4)
	NO=NO+2
	EST1 = FX1 + 4.0*FX2 + FX3
	EST2 = FX3 + 4.0*FX4 + FX5
	ADIFF1 = ADIFF
	DIFF = EST + EST - EST1 - EST2
	IF (LEV-30) 305,800,800
305	ADIFF = ABS(DIFF)
	CRIT = ADIFF - CEPS
	IF (CRIT) 700,700,400
C
C	END OF CENTRAL LOOP
C	     NEXT STAGE IS STAGE FOUR IN CASE OF NO NATURAL CONVERGENCE
C	     NEXT STAGE IS STAGE SEVEN IN CASE OF NATURAL CONVERGENCE
C
C	     ****  STAGE FOUR  ****
C	NO NATURAL CONVERGENCE.  A COMPLEX SEQUENCE OF INSTRUCTIONS
C	FOLLOWS WHICH ASSIGNS CONVERGENCE AND/OR ALTERS TOLERANCE
C	LEVEL IN UPWARD DIRECTION IF THERE ARE INDICATIONS
C	OF ROUND OFF ERROR.
C
400	CONTINUE
	IF (ADIFF1-ADIFF) 410,410,500
C
C	IN A NORMAL RUN WITH NO ROUND OFF ERROR PROBLEM, ADIFF1 IS
C	GREATER THAN ADIFF AND THE REST OF STAGE FOUR IS OMITTED.
C
410	IF (LEV-5) 500,415,415
415	EFACT =EFACT + CEPS *(X1-XZERO)*FACERR
	XZERO = X1
	FACERR = 15.0
C
C	THE REST OF STAGE FOUR DEALS WITH UPWARD ADJUSTMENT OF TOLERANCE
C	     (CEPS) BECAUSE OF SUSPECTED ROUND OFF ERROR TROUBLE.
C
	IF (ADIFF-2.0*CEPS) 420,420,425
C
C	SMALL JUMP IN CEPS.  ASSIGN CONVERGENCE.
C
420	CEPS = ADIFF
	LEVTAG = 0
	GO TO 780
425	IF (ADIFF1-ADIFF) 435,430,435
C
C	LARGE JUMP IN CEPS
C
430	CEPS = ADIFF
	GO TO 445
C
C	FACTOR TWO JUMP IN CEPS
C
435	CEPS = 2.0 * CEPS
	IF (LEVTAG-3) 440,445,445
440	LEVTAG = 2
445	QCEPS = 0.25*CEPS
C
C	     ****  STAGE FIVE  ****
C	NO ACTUAL CONVERGENCE.
C	STORE RIGHT HAND ELEMENTS
C
500	CONTINUE
	NIM= 2* NIM
	LEV = LEV + 1
	ESTST(LEV) = EST2
	X3ST(LEV) = X4
	X5ST(LEV) = X5
	FX3ST(LEV) = FX4
	FX5ST(LEV) = FX5
	PREDIF(LEV) = ADIFF
C
C	     ****  STAGE SIX  ****
C	SET UP QUANTITIES FOR CENTRAL CALCULATION.
C
C	READY TO GO AHEAD AT LEVEL LOWER WITH LEFT HAND ELEMENTS
C	X1 AND FX1 ARE THE SAME AS BEFORE
C
	X5=X3
	X3=X2
	FX5=FX3
	FX3=FX2
	EST = EST1
	GO TO 300
C
C	     ****  STAGE SEVEN  ****
C	NATURAL CONVERGENCE IN PREVIOUS INTERVAL.  THE FOLLOWING
C	COMPLEX SEQUENCE CHECKS PRIMARILY THAT TOLERANCE LEVEL IS
C	NOT TOO HIGH.  UNDER CERTAIN CIRCUMSTANCES NON CONVERGENCE
C	IS ASSIGNED AND/OR TOLERANCE LEVEL IS RE-SET.
C
700	CONTINUE
C
C	CHECK THAT IT WAS NOT LEVEL ZERO INTERVAL.  IF SO ASSIGN
C	NON CONVERGENCE
C
	IF (LEV) 400,400,705
C
C	LEVTAG = -1  CEPS = CEPSF, ITS ORIGINAL VALUE.
C	LEVTAG =  0  CEPS IS GREATER THAN CEPSF.  REGULAR SITUATION.
C	LEVTAG =  2  CEPS IS GREATER THAN CEPSF.  CEPS PREVIOUSLY ASKED
C	             FOR A BIG JUMP, BUT DID NOT GET ONE.
C	LEVTAG =  3  CEPS IS GREATER THAN CEPSF.  CEPS PREVIOUSLY HAD
C	             A BIG JUMP.
C
705	IF (LEVTAG) 800,710,710
C
C	IN A NORMAL RUN WITH NO ROUND OFF ERROR PROBLEM, LEVTAG = -1
C		AND THE REST OF STAGE SEVEN IS OMITTED.
C
710	CEPST = 15.0*CEPS
C
C	CEPST HERE IS FACERR*CURRENT VALUE OF CEPS
C
	IF (CRIT) 715,800,800
715	IF (LEVTAG-2) 720,740,750
C
C	LEVTAG = 0
C
720	IF (ADIFF) 800,800,725
725	IF (ADIFF-QCEPS) 730,800,800
730	IF (ADIFF-CEPSF) 770,770,735
735	LEVTAG = 0
	CEPS = ADIFF
	EFACT = EFACT + CEPST*(X1-XZERO)
	XZERO=X1
	GO TO 445
C
C	LEVTAG = 2
C
740	LEVTAG = 0
	IF (ADIFF) 765,765,725
C
C	LEVTAG = 3
C
750	LEVTAG = 0
	IF (ADIFF) 775,775,730
765	CEPS = ADIFF1
	GO TO 775
770	LEVTAG = -1
	FACERR = 1.0
	CEPS = CEPSF
775	EFACT = EFACT + CEPST*(X1-XZERO)
	XZERO = X1
780	CONTINUE
	QCEPS = 0.25*CEPS
C
C	     ****  STAGE EIGHT  ****
C	ACTUAL CONVERGENCE IN PREVIOUS INTERVAL.  INCREMENTS ADDED INTO
C	RUNNING SUMS
C
C	ADD INTO SUM AND SIM
C
800	CONTINUE
	SUM= SUM+ (EST1+EST2)*(X5-X1)
	IF (LEVTAG) 805,810,810
C
C	WE ADD INTO SIM ONLY IF WE ARE CLEAR OF ROUND OFF LEVEL
C
805	SIM=SIM + DIFF*(X5-X1)
810	CONTINUE
C
C	     ****  STAGE NINE  ****
C	SORT OUT WHICH LEVEL TO GO TO.  THIS INVOLVES NIM NUMBERING
C	SYSTEM DESCRIBED BEFORE STAGE ONE.
C
905	NUM = NIM/2
	NOM= NIM - 2*NUM
	IF (NOM) 910,915,910
910	NIM = NUM
	LEV = LEV - 1
	GO TO 905
915	NIM = NIM + 1
C
C	NEW LEVEL IS SET.  IF LEV = 0 WE HAVE FINISHED
C
	IF (LEV) 1100,1100,1000
C
C	     ****  STAGE TEN  ****
C	SET UP QUANTITIES FOR CENTRAL CALCULATION.
C
1000	CONTINUE
	X1=X5
	FX1=FX5
	X3=X3ST(LEV)
	X5=X5ST(LEV)
	FX3=FX3ST(LEV)
	FX5=FX5ST(LEV)
	EST=ESTST(LEV)
	ADIFF = PREDIF(LEV)
	GO TO 300
C
C	     ****  STAGE ELEVENT  ****
C
1100	CONTINUE
	EFACT = EFACT + CEPS * (BIG-XZERO)*FACERR
	RUM = EFACT/180.0
	THIRD = SUM/12.0
	FIFTH = -SIM/180.0
	SQUANK = THIRD + FIFTH
	RETURN
C
C	END OF SQUANK
C
	END
C---------------X, A ARE RETURNED.  C, D ARE INPUT.
	SUBROUTINE GLO16(X,A,C,D)
	DOUBLE PRECISION DMC,DPC,C,D,X(16),A(16),XX(8),AA(8)
	DATA (XX(I),AA(I),I=1,8)/
     1.9894009349916499325961541734D0 ,.2715245941175409485178057245D-1,
     1.9445750230732325760779884155D0 ,.6225352393864789286284383699D-1,
     1.8656312023878317438804678977D0 ,.9515851168249278480992510760D-1,
     1.7554044083550030338951011948D0 ,.1246289712555338720524762821D0 ,
     1.6178762444026437484466717640D0 ,.1495959888165767320815017305D0 ,
     1.4580167776572273863424194429D0 ,.1691565193950025381893120790D0 ,
     1.2816035507792589132304605014D0 ,.1826034150449235888667636679D0 ,
     1.9501250983763744018531933542D-1,.1894506104550684962853967232D0/
	DMC = .5D0*(D-C)
	DPC = .5D0*(D+C)
	DO 2 I=1,8
	NI = 17 - I
	X(I) = -DMC*XX(I) + DPC
	X(NI)=  DMC*XX(I) + DPC
	A(I) =  DMC*AA(I)
2	A(NI)=  DMC*AA(I)
	RETURN
	END

