      SUBROUTINE SETUP(NGUIDE,GCON)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MAXGUIDE=7,MAXPTS=65)
      DIMENSION GCON(MAXGUIDE)
      COMMON /PTS/ SPTS(MAXPTS),SVAL(MAXPTS),SINT(MAXPTS),SDER(MAXPTS)
      EXTERNAL POLY
C
C SET UP A BIASED AS RANDOM GUIDING FUNCTION AS EXPONENTIAL OF
C A SPLINE  -- USE BLI TO FIND WHERE TO EVALUATE IT --
C SPTS(I) CONTAINS THE X VALUES OF THE INTERPOLATING SPLINE
C SVAL(I) CONTAINS THE Y VALUES OF THE INTERPOLATING SPLINE
C SINT(I) CONTAINS THE INTEGRAL [0,X(I)] OF SVAL(I)
C SDER(I) CONTAINS THE DERIVATIVE OF SVAL(I)
C MAXPTS IS THE NUMBER OF POINTS IN THE INTERPOLATION
C
      A=0
      B=32
      FTEST=BLI(A,B,POLY,SPTS,SVAL,SINT,BLIER,MAXPTS)
      WRITE(6,'('' IN SETUP A='',G10.3,'' B='',G10.3,'' NGUIDE='',I5/
     2 '' GCON=''/(4G20.6))') A,B,NGUIDE,(GCON(I),I=1,NGUIDE)
C      DO 20 I=1,MAXPTS
C      WRITE(6,'(2G20.6)') SPTS(I),SVAL(I)
C20    CONTINUE
      WRITE(6,*)' FTEST=',FTEST,'+-',BLIER
      SINT(1)=0.0
      DO 30 I=2,MAXPTS
      IM=I-1
      SINT(I)=.5*(SPTS(I)-SPTS(IM))*(SVAL(I)+SVAL(IM))+SINT(IM)
      SDER(IM)=(SVAL(I)-SVAL(IM))/(SPTS(I)-SPTS(IM))
30    CONTINUE
      WRITE(6,*)' ST LINE INTEGRAL ',SINT(MAXPTS-1)
      RETURN
      END
C
      FUNCTION BLI(A,B,FUN,XI,FI,ERR,BLIER,NMAX)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION XI(NMAX),FI(NMAX),ERR(NMAX)
C
C *** USE HALF THE POINTS IN A MIDPOINT MESH
C
      NDO=NMAX/4
      H=(B-A)/(NDO-1)
      DO I=1,NDO
        XI(I)=A+H*(I-1)
        FI(I)=FUN(XI(I))
      ENDDO
      AIC=0
      DO I=2,NDO
        AIC=AIC+.5D0*H*(FI(I)+FI(I-1))
      ENDDO
      AIN=AIC
      H3=H**3
C *** ERR(1) IS THE ERROR IN THE INTERVAL A TO X(1), FPP USES PTS 1,2,3
      FPL=(FI(2)-FI(1))/(XI(2)-XI(1))
      FPU=(FI(3)-FI(2))/(XI(3)-XI(2))
      FPP=2*(FPU-FPL)/(XI(3)-XI(1))
C *** ERR(1) IS THE ERROR IN THE INTERVAL X(1) TO X(2), USES SAME FPP AS
      ERR(1)=ABS(8.333333D-2*H3*FPP)
C *** MID REGION USES I+1 TO I+2 FOR FPU AND I-1 TO I FOR FPL
      FPTEMP=FPU
      DO I=2,NDO-2
        FPU=(FI(I+2)-FI(I+1))/(XI(I+2)-XI(I+1))
        FPP=2*(FPU-FPL)/(XI(I+2)+XI(I+1)-XI(I)-XI(I-1))
        ERR(I)=ABS(8.333333D-2*H3*FPP)
        FPL=FPTEMP
        FPTEMP=FPU
      ENDDO
C *** INTERVAL FROM XI(N-1) TO XI(N) USES FPU AND FPTEMP TO ESTIMATE FPP
      FPP=2*(FPU-FPL)/(XI(NDO)-XI(NDO-2))
      ERR(NDO-1)=ABS(8.333333D-2*H3*FPP)
      ERR(NDO)=0
C *** INITIALIZING THE CURVE FIT PARTS
      SMA=0
      SMB=0
      SMD=0
      B0=0
      B1=0
      SSQ=0
C *** FINDING THE LARGEST ERROR AND THE TOTAL SUM OF THE ERRORS
50    ERRM=ERR(1)
      ILE=1
      ERRS=ERRM
      DO I=2,NDO-1
        ERRS=ERRS+ERR(I)
        IF(ERR(I).GT.ERRM)THEN
          ERRM=ERR(I)
          ILE=I
        ENDIF
      ENDDO
C *** UPDATING THE FITTING PARTS
      ERRS=DMAX1(ERRS,1.D-14*AIN)
      IF(ERRS.GT.0.D0)THEN
        ERRI=1/ERRS
        ERRI2=ERRI*ERRI
        SMA=SMA+ERRI2
        SMB=SMB+ERRI
        SMD=SMD+1
        AINMAC=AIN-AIC
        B0=B0+AINMAC*ERRI2
        B1=B1+AINMAC*ERRI
        SSQ=SSQ+AINMAC*AINMAC*ERRI2
      ENDIF
      IF(NDO.EQ.NMAX)THEN
        IF(SMD.EQ.0.D0)THEN
          BLI=0
          BLIER=0
          RETURN
        ENDIF
        DEN=SMA*SMD-SMB*SMB
        IF(DEN.LT.1.D-14*SMA*SMD)THEN
          BLI=B0/SMA+AIC
          BLIER=1.D-12*BLI
          RETURN
        ENDIF
        C0=(SMD*B0-SMB*B1)/DEN
        C1=(-SMB*B0+SMA*B1)/DEN
        CHIS=DMAX1(1.D-14*SSQ,(SSQ-C0*B0-C1*B1))
        BLI=C0+AIC
        BLIER=DSQRT(CHIS*SMD/DEN)
        RETURN
      ENDIF
C *** POINT TO BE INSERTED WILL BE ILE+1, THUS CURRENT XI(ILE+1) BECOMES
C *** XI(ILE+2)
      DO I=NDO,ILE+1,-1
        XI(I+1)=XI(I)
        FI(I+1)=FI(I)
        ERR(I+1)=ERR(I)
      ENDDO
      NDO=NDO+1
      IADD=ILE+1
      XI(IADD)=.5D0*(XI(ILE)+XI(ILE+2))
      FI(IADD)=FUN(XI(IADD))
      AIN=AIN-.5D0*(XI(ILE+2)-XI(ILE))*(FI(ILE)+FI(ILE+2))
     2   +.5D0*(XI(ILE+1)-XI(ILE))*(FI(ILE)+FI(ILE+1))
     3   +.5D0*(XI(ILE+2)-XI(ILE+1))*(FI(ILE+1)+FI(ILE+2))
      IF(ILE.EQ.1)THEN
        FPL=(FI(2)-FI(1))/(XI(2)-XI(1))
        FPU=(FI(3)-FI(2))/(XI(3)-XI(2))
        FPP=2*(FPU-FPL)/(XI(3)-XI(1))
C *** ERR(1) IS THE ERROR IN THE INTERVAL X(1) TO X(2)
        ERR(1)=ABS(8.333333D-2*(XI(2)-XI(1))**3*FPP)
C *** ERR3 IS BASED ON 12 FOR FPL AND 34 FOR FPU
        FTEMP=FPU
        FPU=(FI(4)-FI(3))/(XI(4)-XI(3))
        FPP=2*(FPU-FPL)/(XI(4)+XI(3)-XI(2)-XI(1))
        ERR(2)=ABS(8.333333D-2*(XI(3)-XI(2))**3*FPP)
        IF(NDO.GE.5)THEN
          FPL=FTEMP
          FPU=(FI(5)-FI(4))/(XI(5)-XI(4))
          FPP=2*(FPU-FPL)/(XI(4)+XI(4)-XI(3)-XI(2))
        ENDIF
        ERR(3)=ABS(8.333333D-2*(XI(4)-XI(3))**3*FPP)
        GOTO 50
      ENDIF
      IF(ILE.EQ.NDO-1)THEN
        FPU=(FI(NDO)-FI(NDO-1))/(XI(NDO)-XI(NDO-1))
        FPL=(FI(NDO-1)-FI(NDO-2))/(XI(NDO-1)-XI(NDO-2))
        FTEMP=FPL
        FPP=2*(FPU-FPL)/(XI(NDO)-XI(NDO-2))
        ERR(NDO-1)=ABS(8.333333D-2*FPP*(XI(NDO)-XI(NDO-1))**3)
        FPL=(FI(NDO-2)-FI(NDO-3))/(XI(NDO-2)-XI(NDO-3))
        FPP=2*(FPU-FPL)/(XI(NDO)+XI(NDO-1)-XI(NDO-2)-XI(NDO-3))
        ERR(NDO-2)=ABS(8.333333D-2*(XI(NDO-1)-XI(NDO-2))**3*FPP)
        IF(NDO-4.GT.0)THEN
          FPL=(FI(NDO-3)-FI(NDO-4))/(XI(NDO-3)-XI(NDO-4))
          FPU=FTEMP
          FPP=2*(FPU-FPL)/(XI(NDO-1)+XI(NDO-2)-XI(NDO-3)-XI(NDO-4))
        ENDIF
        ERR(NDO-3)=ABS(8.333333D-2*(XI(NDO-2)-XI(NDO-3))**3*FPP)
        GOTO 50
      ENDIF
C *** NEED TO RECALCULATE ERR(ILE-2) TO ERR(ILE+1)
      IF(ILE.GT.4)THEN
        FPL=(FI(ILE-2)-FI(ILE-3))/(XI(ILE-2)-XI(ILE-3))
        FPU=(FI(ILE)-FI(ILE-1))/(XI(ILE)-XI(ILE-1))
        FPP=2*(FPU-FPL)/(XI(ILE)+XI(ILE-1)-XI(ILE-2)-XI(ILE-3))
        FPL=(FI(ILE-1)-FI(ILE-2))/(XI(ILE-1)-XI(ILE-2))
      ELSE
        ILE=MAX0(3,ILE)
        FPL=(FI(ILE-1)-FI(ILE-2))/(XI(ILE-1)-XI(ILE-2))
        FPU=(FI(ILE)-FI(ILE-1))/(XI(ILE)-XI(ILE-1))
        FPP=2*(FPU-FPL)/(XI(ILE)-XI(ILE-2))
      ENDIF
      ERR(ILE-2)=ABS(8.333333D-2*(XI(ILE-1)-XI(ILE-2))**3*FPP)
      FTEMP=FPU
      FPU=(FI(ILE+1)-FI(ILE))/(XI(ILE+1)-XI(ILE))
      FPP=2*(FPU-FPL)/(XI(ILE+1)+XI(ILE)-XI(ILE-1)-XI(ILE-2))
      ERR(ILE-1)=ABS(8.333333D-2*(XI(ILE)-XI(ILE-1))**3*FPP)
      FPL=FTEMP
      FTEMP=FPU
      FPU=(FI(ILE+2)-FI(ILE+1))/(XI(ILE+2)-XI(ILE+1))
      FPP=2*(FPU-FPL)/(XI(ILE+2)+XI(ILE+1)-XI(ILE)-XI(ILE-1))
      ERR(ILE)=ABS(8.333333D-2*(XI(ILE+1)-XI(ILE))**3*FPP)
      FPL=FTEMP
      IF(ILE+2.EQ.NDO)THEN
        FPP=2*(FPU-FPL)/(XI(ILE+2)-XI(ILE))
      ELSE
        FPU=(FI(ILE+3)-FI(ILE+2))/(XI(ILE+3)-XI(ILE+2))
        FPP=2*(FPU-FPL)/(XI(ILE+3)+XI(ILE+2)-XI(ILE+1)-XI(ILE))
      ENDIF
      ERR(ILE+1)=ABS(8.333333D-2*(XI(ILE+2)-XI(ILE+1))**3*FPP)
      GOTO 50
      END
C
      FUNCTION POLY(X)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (MAXGUIDE=7)
      COMMON /GUIDE/ GCON(MAXGUIDE),NGUIDE
C *** The first coefficient is an exponential c(i)*x
C *** The rest are in the form c(i)*exp((x-c(i+1))/c(i+2))**2
      POLY=GCON(1)*X
      NEXP=(NGUIDE-1)/3
      IOFF=1
      DO I=1,NEXP
        POLY=POLY+GCON(IOFF+1)*EXP(-((X-GCON(IOFF+2))/GCON(IOFF+3))**2)
        IOFF=IOFF+3
      ENDDO
      POLY=EXP(POLY)
      RETURN
      END
C
      SUBROUTINE LOCATE(ARR,M,X,J)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION ARR(1)
C
C FIND THE LAST ELEMENT IN ARR(J) S.T. ARR(J)<X
C J IS NEEDED FOR INTERPOLATIONS BUT J+1 IS NEEDED FOR POINT SELECTION
C ITS ASSUMED THAT ARR(1) <= ARR(2) <= ... ARR(M)
C REFERENCE: NUMERICAL RECIPIES PG. 90
C
      JL=0
      JU=M+1
10    IF(JU-JL.GT.1) THEN
        JM=(JU+JL)/2
        IF(X.GT.ARR(JM)) THEN
          JL=JM
        ELSE
          JU=JM
        ENDIF
        GO TO 10
      ENDIF
      J=JL
      RETURN
      END
C
