*
* $Id: gpindp.F,v 1.1.1.1 1996/04/01 15:03:12 mclareni Exp $
*
* $Log: gpindp.F,v $
* Revision 1.1.1.1  1996/04/01 15:03:12  mclareni
* Mathlib gen
*
*
#include "sys/CERNLIB_machine.h"
#include "_gen/pilot.h"
      DOUBLE PRECISION FUNCTION GPINDP(A,B,EPSIN,EPSOUT,FUNC,IOP)
C
C     PARAMETERS
C
C     A       = LOWER BOUNDARY
C     B       = UPPER BOUNDARY
C     EPSIN   = ACCURACY REQUIRED FOR THE APPROXINATION
C     EPSOUT  = IMPROVED ERROR ESTIMATE FOR THE APPROXIMATION
C     FUNC    = FUNCTION ROUTINE FOR THE FUNCTION FUNC(X).TO BE DE-
C               CLARED EXTERNAL IN THE CALLING ROUTINE
C     IOP     = OPTION PARAMETER , IOP=1 , MODIFIED ROMBERG ALGORITHM,
C                                          ORDINARY CASE
C                                  IOP=2 , MODIFIED ROMBERG ALGORITHM,
C                                          COSINE TRANSFORMED CASE
C                                  IOP=3 , MODIFIED CLENSHAW-CURTIS AL
C                                          GORITHM
C
C     PARAMETERS IN COMMON BLOCK / GPINT /
C
C     TEND    = UPPER BOUND FOR VALUE OF INTEGRAL
C     UMID    = LOWER BOUND FOR VALUE OF INTEGRALC
C     N       = THE NUMBER OF INTEGRAND VALUES USED IN THE CALCULATION
C     LINE    = LINE NO IN ROMBERG TABLE (RELATED TO N THROUGH
C               N-1=2**(LINE-1) , APPLICABLE ONLY FOR IOP=1 OR 2)
C     IOUT    = ELEMENT NO IN LINE (APPLICABLE ONLY FOR IOP=1 OR 2)
C     JOP     = OPTION PARAMETER , JOP=0 , NO PRINTING OF INTERMEDIATE
C                                          CALCULATIONS
C                                  JOP=1 , PRINT INTERMEDIATE CALCULA-
C                                          TIONS
C     KOP     = OPTION PARAMETER , KOP=0 , NO TIME ESTIMATE
C                                  KOP=1 , ESTIMATE TIME
C     T       = TIME USED FOR CALCULATION IN MSEC.
C
C     INTEGRATION PARAMETERS
C
C     NUPPER  = 9 , CORRESPONDS TO 1024 SUB-INTERVALS FOR THE UNFOLDED
C               INTEGRAL.THE MAX.NO OF FUNCTION EVALUATIONS THUS BEEING
C               1025.THE HIGHEST END-POINT APPROXIMATION IS THUS USING
C               1024 INTERVALS WHILE THE HIGHEST MID-POINT APPROXIMA-
C               TION IS USING 512 INTERVALS.
C
C     INPUT/OUTPUT PARAMETERS
C
      EXTERNAL FUNC
      DOUBLE PRECISION A,B,EPSIN,EPSOUT,BOUND,FUNC
C
C     INTERNAL ARRAYS
C
      DOUBLE PRECISION ACOF(513),BCOF(513),CCOF(1025)
C
C     CONSTANTS IN DATA STATEMENTS
C
C*NS  DOUBLE PRECISION ZERO,FOURTH,HALF,ONE,TWO,THREE,FOUR,FAC1,FAC2,PI,
C*NS 1RANDER
      DOUBLE PRECISION ZERO,FOURTH,HALF,ONE,TWO,      FOUR,FAC1,FAC2,PI
C
C     VARIABLES DEPENDING ON STEPSIZE
C
      DOUBLE PRECISION ALF,BET,RN,HNSTEP,TEND,UMID,WMEAN,DELN,TNEW,AR
C
C     CONSTANTS RELATED TO CALCULATION OF TRIGONOMETRIC FUNCTIONS
C
      DOUBLE PRECISION TRIARG,ALFN0,BETN0,GAMMAN,DELTAN,ALFNJ,BETNJ,ETAN
     1K,KSINK
C
C     OTHER VARIABLES USED
C
C*NS  DOUBLE PRECISION CONST1,CONST2,XPLUS,XMIN,ERROR,RK,A0,A1,A2,COF,FA
C*NS 1CTOR,ENDPTS
      DOUBLE PRECISION CONST1,CONST2,XPLUS,XMIN,      RK,A0,A1,A2,COF,FA
     1CTOR,ENDPTS
C
      COMMON /GPINT/ TEND, UMID, N, LINE, IOUT, JOP, KOP, T
      DATA ZERO,FOURTH,HALF,ONE,TWO,FOUR/0.D0,.25D0,.5D0,1.D0,2.D0,4.D0/
      DATA PI,FAC1,FAC2/3.141592653589793238462643383279D0,.411233516712
     1056609118103791649D0,.822467033424113218236207583298D0/
C
      DATA NUPPER/9/
C
C     TIMEX(T) IS A LIBRARY SUBROUTINE GIVING THE ELAPSED CP TIME
C
      IF(KOP .NE. 0)  CALL TIMEX ( T1)
C
C     INITIAL CALCULATIONS
C
      ALF=HALF*(B-A)
      BET=HALF*(B+A)
      CONST1=FUNC(A)+FUNC(B)
      CONST2=FUNC(BET)
      HNSTEP=TWO
      IF(IOP.EQ.2) HNSTEP=PI
C
      IF(IOP.GT.1) GOTO 10
C
C     MODIFIED ROMBERG ALGORITHM,ORDINARY CASE
C
      BCOF(1)=HNSTEP*CONST2
      ACOF(1)=HALF*(CONST1+BCOF(1))
      FACTOR=ONE
      ACOF(2)=ACOF(1)-(ACOF(1)-BCOF(1))/(FOUR*FACTOR-ONE)
      GOTO 30
C
10    IF(IOP.GT.2) GOTO 20
C
C     MODIFIED ROMBERG ALGORITHM,COSINE TRANSFORMED CASE
C
      AR=FAC1
      ENDPTS=CONST1
      ACOF(1)=FAC2*CONST1
      BCOF(1)=HNSTEP*CONST2-AR*CONST1
      FACTOR=FOUR
      ACOF(1)=HALF*(ACOF(1)+BCOF(1))
      ACOF(2)=ACOF(1)-(ACOF(1)-BCOF(1))/(FOUR*FACTOR-ONE)
      AR=FOURTH*AR
      GOTO 30
C
20    CONST1=HALF*CONST1
      ACOF(1)=HALF*(CONST1+CONST2)
      ACOF(2)=HALF*(CONST1-CONST2)
      BCOF(2)=ACOF(2)
      TEND=TWO*(ACOF(1)-ACOF(2)/(ONE+TWO))
C
C     MODIFIED CLENSHAW-CURTIS ALGORITHM
C
30    HNSTEP=HALF*HNSTEP
      NHALF=1
      N=2
      RN=TWO
C
      IF(IOP.NE.1) THEN
C
C     INITIAL PARAMETERS SPECIAL FOR THE MODIFIED ROMBERG ALGORITHM,
C     COSINE TRANSFORMED CASE AND THE MODIFIED CLENSHAW-CURTIS ALGORITHM
C
        TRIARG=FOURTH*PI
        ALFN0=-ONE
      ENDIF
C
C     END OF INITIAL CALCULATIONS
C
C     START ACTUAL CALCULATION
C
C---  Transform this DO-loop into a GOTO to avoid illegal jumps into it
C
C     DO 350 I=1,NUPPER
      I=0
41    I=I+1
      IF(I.GT.NUPPER) GOTO 350
      LINE=I+2
C
      IF(IOP.GT.1) GOTO 60
C
C     MODIFIED ROMBERG ALGORITHM,ORDINARY CASE
C
C     COMPUTE FIRST ELEMENT IN MID-POINT FORMULA FOR ORDINARY CASE
C
      UMID=ZERO
      ALFNJ=HALF*HNSTEP
      DO 50 J=1,NHALF
      XPLUS=ALF*ALFNJ+BET
      XMIN=-ALF*ALFNJ+BET
      UMID=UMID+FUNC(XPLUS)+FUNC(XMIN)
      ALFNJ=ALFNJ+HNSTEP
50    CONTINUE
      UMID=HNSTEP*UMID
      GOTO 100
C
C     COMPUTE FUNCTION VALUES FOR MODIFIED ROMBERG ALGORITHM,COSINE
C     TRANSFORMED CASE AND MODIFIED CLENSHAW-CURTIS ALGORITHM
C
60    CONST1=-SIN(TRIARG)
      CONST2=HALF*ALFN0/CONST1
      IF(IOP.EQ.2) ETANK=CONST2
      ALFN0=CONST1
      BETN0=CONST2
      GAMMAN=ONE-TWO*ALFN0**2
      DELTAN=-TWO*ALFN0*BETN0
C
      DO 70 J=1,NHALF
      ALFNJ=GAMMAN*CONST1+DELTAN*CONST2
      BETNJ=GAMMAN*CONST2-DELTAN*CONST1
      XPLUS=ALF*ALFNJ+BET
      XMIN=-ALF*ALFNJ+BET
      CCOF(J)=FUNC(XPLUS)+FUNC(XMIN)
      CONST1=ALFNJ
      CONST2=BETNJ
70    CONTINUE
C
      IF(IOP.EQ.3) GOTO 190
C
C     COMPUTE FIRST ELEMENT IN MID-POINT FORMULA FOR COSINE TRANSFORMED
C     ROMBERG ALGORITHM
C
      NCOF=NHALF-1
      COF=TWO*(TWO*ETANK**2-ONE)
      A2=ZERO
      A1=ZERO
      A0=CCOF(NHALF)
      IF(NCOF.EQ.0) GOTO 90
      DO 80 J=1,NCOF
      A2=A1
      A1=A0
      INDEX=NHALF-J
      A0=CCOF(INDEX)+COF*A1-A2
80    CONTINUE
90    UMID=HNSTEP*(A0-A1)*ETANK-AR*ENDPTS
      AR=FOURTH*AR
C
C     MODIFIED ROMBERG ALGORITHM,CALCULATE (I+1)-TH ROW IN U-TABLE
C
100   CONST1=FOUR*FACTOR
      INDEX=I+1
      DO 110 J=2,INDEX
      TEND=UMID+(UMID-BCOF(J-1))/(CONST1-ONE)
      BCOF(J-1)=UMID
      UMID=TEND
      CONST1=FOUR*CONST1
110   CONTINUE
      BCOF(INDEX)=TEND
      XPLUS=CONST1
C
C     CALCULATION OF (I+1)-TH ROW IN U-TABLE FINISHED
C
C     PRINT INTERMEDIATE RESULTS IF WANTED
C
      IF(JOP.EQ.0) GOTO 120
C
      ICHECK=0
      ASSIGN 120 TO JUMP
      GOTO 360
C
C     TEST IF REQUIRED ACCURACY IS OBTAINED
C
120   EPSOUT=ONE
      IOUT=1
      DO 140 J=1,INDEX
      CONST1=HALF*(ACOF(J)+BCOF(J))
      CONST2=HALF*ABS((ACOF(J)-BCOF(J))/CONST1)
      IF(CONST2.GT.EPSOUT) GOTO 130
      EPSOUT=CONST2
      IOUT=J
130   ACOF(J)=CONST1
140   CONTINUE
C
C     TESTING ON ACCURACY FINISHED
C
      IF(IOUT.EQ.INDEX) IOUT=IOUT+1
      ACOF(INDEX+1)=ACOF(INDEX)-(ACOF(INDEX)-BCOF(INDEX))/(XPLUS-ONE)
C
      IF(EPSOUT.GT.EPSIN) GOTO 340
C
C     CALCULATION FOR MODIFIED ROMBERG ALGORITHM FINISHED
C
150   N=2*N
C
C     PRINT INTERMEDIATE RESULTS IF WANTED
C
      IF(JOP.EQ.0) GOTO 170
C
      ICHECK=1
      INDEX=INDEX+1
      ASSIGN 160 TO JUMP
      GOTO 360
C
160   INDEX=INDEX-1
C
170   N=N+1
      J=IOUT
      IF((J-1).LT.INDEX) GOTO 180
      J=INDEX
180   TEND=ALF*(TWO*ACOF(J)-BCOF(J))
      UMID=ALF*BCOF(J)
      GPINDP=ALF*ACOF(IOUT)
C
      GOTO 310
C
C     START CALCULATION FOR MODIFIED CLENSHAW-CURTIS ALGORITHM
C
190   BCOF(1)=ZERO
      DO 200 J=1,NHALF
      BCOF(1)=BCOF(1)+CCOF(J)
200   CONTINUE
      BCOF(1)=HALF*HNSTEP*BCOF(1)
C
C     CALCULATION OF FIRST B-COEFFICIENT FINISHED.COMPUTE THE HIGHER
C     COEFFICIENTS IF NHALF GREATER THAN ONE
C
      IF(NHALF.EQ.1) GOTO 230
C
      CONST1=ONE
      CONST2=ZERO
      NCOF=NHALF-1
      KSIGN=-1
      DO 220 K=1,NCOF
C
C     COMPUTE TRIGONOMETRIC SUM FOR B-COEFFICIENT
C
      ETANK=GAMMAN*CONST1-DELTAN*CONST2
      KSINK=GAMMAN*CONST2+DELTAN*CONST1
      COF=TWO*(TWO*ETANK**2-ONE)
      A2=ZERO
      A1=ZERO
      A0=CCOF(NHALF)
      DO 210 J=1,NCOF
      A2=A1
      A1=A0
      INDEX=NHALF-J
      A0=CCOF(INDEX)+COF*A1-A2
210   CONTINUE
C
      BCOF(K+1)=HNSTEP*(A0-A1)*ETANK
      IF(KSIGN.EQ.-1) BCOF(K+1)=-BCOF(K+1)
      KSIGN=-KSIGN
C
      CONST1=ETANK
      CONST2=KSINK
C
220   CONTINUE
C
C     CALCULATION OF B-COEFFICIENTS FINISHED
C
C     COMPUTE NEW MODIFIED MID-POINT APPROXIMATION WHEN THE INTERVAL
C     OF INTEGRATION IS DIVIDED IN N EQUAL SUB INTERVALS
C
230   UMID=ZERO
      RK=RN
      NN=NHALF+1
      DO 240 K=1,NN
      INDEX=NN+1-K
      UMID=UMID+BCOF(INDEX)/(RK**2-ONE)
      RK=RK-TWO
240   CONTINUE
      UMID=-TWO*UMID
C
C     COMPUTE NEW C-COEFFICIENTS FOR END-POINT APPROXIMATION
C
      NN=N+2
      DO 250 J=1,NHALF
      INDEX=NN-J
      CCOF(J)=HALF*(ACOF(J)+BCOF(J))
      CCOF(INDEX)=HALF*(ACOF(J)-BCOF(J))
250   CONTINUE
      INDEX=NHALF+1
      CCOF(INDEX)=ACOF(INDEX)
C
C     CALCULATION OF NEW COEFFICIENTS FINISHED
C
C     COMPUTE NEW END-POINT APPROXIMATION WHEN THE INTERVAL OF INTEGRA-
C     TION IS DIVIDED IN 2N EQUAL SUB INTERVALS
C
      WMEAN=HALF*(TEND+UMID)
      BOUND=HALF*(TEND-UMID)
C
      DELN=ZERO
      RK=TWO*RN
      DO 260 J=1,NHALF
      INDEX=N+2-J
      DELN=DELN+CCOF(INDEX)/(RK**2-ONE)
      RK=RK-TWO
260   CONTINUE
      DELN=-TWO*DELN
C
C     PRINT INTERMEDIATE RESULTS IF WANTED
C
      IF(JOP.EQ.0) GOTO 270
C
      GOTO 400
C
C     PRINTING OF INTERMEDIATE RESULTS FINISHED
C
270   TNEW=WMEAN+DELN
      EPSOUT=ABS(BOUND/TNEW)
C
      IF(EPSOUT.GT.EPSIN) GOTO 320
C
C     REQUIRED ACCURACY OBTAINED
C
280   N=2*N+1
C
C*UL 290   TEND=ALF*(TEND+DELN)
      TEND=ALF*(TEND+DELN)
      UMID=ALF*(UMID+DELN)
C
C*UL 300   GPINDP=ALF*TNEW
      GPINDP=ALF*TNEW
C
310   IF(KOP.EQ.0) GOTO 315
      CALL TIMEX ( T)
      T=1000.*(T - T1)
C
 315  RETURN
C
320   DO 330 J=1,N
      ACOF(J)=CCOF(J)
330   CONTINUE
      ACOF(N+1)=CCOF(N+1)
      BCOF(N+1)=CCOF(N+1)
      TEND=TNEW
C
340   NHALF=N
      N=2*N
      RN=TWO*RN
      HNSTEP=HALF*HNSTEP
      IF(IOP.GT.1) TRIARG=HALF*TRIARG
C
      GOTO 41
350   CONTINUE
C
C     REQUIRED ACCURACY OF INTEGRAL NOT OBTAINED
C
      N=NHALF
      RN=HALF*RN
C
      IF(IOP.LT.3) GOTO 150
C
      TEND=TWO*(TNEW-DELN)-UMID
C
      GOTO 280
C     PRINT INTERMEDIATE RESULTS FOR THE MODIFIED ROMBERG ALGORITHM
C
360   IF((N.NE.2).AND.(N.NE.256)) GOTO 370
      IF(N.EQ.256) WRITE(6,460)
      WRITE(6,420)
370   DO 390 J=1,INDEX
      CONST1=ALF*ACOF(J)
      IF(ICHECK.EQ.1) GOTO 380
      CONST2=ALF*BCOF(J)
      IF(J.EQ.1) WRITE(6,430) N,J,CONST1,CONST2
      IF(J.GT.1) WRITE(6,440) J,CONST1,CONST2
      GOTO 390
C
380   IF(J.EQ.1) WRITE(6,430) N,J,CONST1
      IF(J.GT.1) WRITE(6,440) J,CONST1
390   CONTINUE
      GOTO JUMP,(120,160)
C
C     PRINTING FINISHED FOR THE MODIFIED ROMBERG ALGORITHM
C
C     PRINT INTERMEDIATE RESULTS FOR THE MODIFIED CLENSHAW-CURTIS AL-
C     GORITHM
C
400   A0=ALF*TEND
      A1=ALF*WMEAN
      A2=ALF*UMID
      CONST1=ALF*BOUND
      CONST2=ALF*DELN
C
      IF(N.GT.2) GOTO 410
C
      WRITE(6,470)
410   WRITE(6,480) N,A0
      WRITE(6,490) A1,CONST2,CONST1
      WRITE(6,490) A2
      GOTO 270
C
C     PRINTING FINISHED FOR THE MODIFIED CLENSHAW-CURTIS ALGORITHM
C
420   FORMAT(/,8X,'N',3X,'J',19X,'TEND(J)',34X,'UMID(J)',/)
430   FORMAT(5X,I4,2X,I2,4X,D36.29,5X,D36.29)
440   FORMAT(11X,I2,4X,D36.29,5X,D36.29)
450   FORMAT(/)
460   FORMAT('1'////)
470   FORMAT(6X,'N',17X,'TEND  ',/,20X,'(TEND+UMID)/2',28X,'DELN',34X,
     1 '(TEND-UMID)/2'/24X,'UMID')
480   FORMAT(/,4X,I3,2X,D36.29)
490   FORMAT(9X,D36.29,3X,D36.29,3X,D36.29)
      END