*     HISCOMP file

      SUBROUTINE HISCPR(H,H2,N,IB,NP)
*
*     compress histogram
*                    - -
*        CALL HISPRR(H,N,B,NP)
*                        - --
*     where H(1:N)    = histogram content
*           H2(1:N)   = error **2 content
*           B(1:NP+1) = bin compression data (bin limits)
*
*           Ith bin: FIX(B(I)+1.5) ... FIX(B(I+1)+0.5)
*
*     Criteria: neff = (Y/dY)**2
*               NP as a function of Neff
*
*
*
      REAL H(*),H2(*)
      REAL YY(256),Y1(256),Y2(256),YE(256),FD(256),Z(256),T(256)
      INTEGER IB(*)
*     ...
C 111 FORMAT(1X,A/(1X,10F10.3))
C     WRITE(*,*) 'HCOMPX executing ......................'
      NP=0
*     determine indices IA...IB of nonzero histogram bins
      DO I=1,N
       IF(H(I).NE.0.0) GOTO 10
      END DO
      GOTO 100
   10 IA=I-1
      IB(1)=IA
      DO I=N,1,-1
       IF(H(I).NE.0.0) GOTO 20
      END DO
   20 IZ=I
*     copy nonzero part of histogram... IA+1 ... IA+NB=IZ
      NB=IZ-IA
      SUMI=0.0
      SUME=0.0
      DO I=1,NB
       YY(I)=H(IA+I)
       SUMI=SUMI+ABS(H(IA+I))
       IF(H2(IA+I).GT.0.0) THEN
          YEFF=H(IA+I)**2/H2(IA+I)
          YE(I)=YEFF
          SUME=SUME+YEFF
       END IF
      END DO
C     WRITE(*,*) 'SUMI,SUME',SUMI,SUME,SUME/SUMI
*     mean scaling factor
      FAC=SQRT(SUME/SUMI)
C     WRITE(*,*) IA,IZ,NB
C     WRITE(*,111) 'Array YY', (YY(L),L=1,NB)
*     ... and smooth the histogram data
      CALL SPLFT(YY,NB)
C     WRITE(*,111) 'Array YY after smooth', (YY(L),L=1,NB)
*     transformation to fixed bin variance
      DO I=1,NB
       YY(I)=FAC*SQRT(ABS(YY(I)))
      END DO
C     WRITE(*,111) 'Array YY after sqrt', (YY(L),L=1,NB)
*     determine first derivative ...
      DO I=1,NB-1
       Y1(I)=YY(I+1)-YY(I)
      END DO
      Y1(NB)=0.0
*     ... and second derivative ...
      DO I=2,NB-1
       Y2(I)=Y1(I)-Y1(I-1)
      END DO
      Y2(1)=Y2(2)
      Y2(NB)=Y2(NB-1)
      SUM1=0.0
      SUM2=0.0
      DO I=1,NB
        Y1(I)=Y1(I)**2
        Y2(I)=Y2(I)**2
        SUM1=SUM1+Y1(I)
        SUM2=SUM2+Y2(I)
      END DO
C     WRITE(*,111) 'Array Y1 squared     ', (Y1(L),L=1,NB)
C     WRITE(*,111) 'Array Y2 squared     ', (Y2(L),L=1,NB)
C     WRITE(*,*) 'SUM1 SUM2',SUM1,SUM2
*     effective events per bin
      EVBIN=SUME/FLOAT(NB)
      RATL=2.00
*     next two parameters are for "structure" of histogram
      RAT1=0.15
      RAT2=0.05
      DO I=1,NB
       FD(I)=1.0+ALOG(AMAX1(1.0,YE(I)))/RATL
     +          +SQRT(Y1(I)/RAT1+Y2(I)/RAT2)
      END DO
C     WRITE(*,111) 'Array FD             ', (FD(L),L=1,NB)
      DO ITER=1,4
*      hanning (first and last)
       FD1=0.25*(3.0*FD(1)+FD(  2))
       FDN=0.25*(3.0*FD(N)+FD(NB-1))
       ZL    =FD1
*      hanning
       DO I=1,NB-2
        ZN    =0.25*(FD(I)+FD(I+1)+FD(I+1)+FD(I+2))
        FD(I) =ZL
        ZL    =ZN
       END DO
       FD(NB-1)=ZL
       FD(NB  )=FDN
C      WRITE(*,111) 'Array FD             ', (FD(L),L=1,NB)
      END DO
*     initial state: 1 bin each
      DO I=1,NB
       IB(I+1)=1
       FD(NB+2-I)=FD(NB+1-I)
      END DO
      FD(1)=0.0
C     WRITE(*,111) 'Array FD             ', (FD(L),L=1,NB+1)
*
      NP=NB
      FDMAX=0.0
*     find maximum
      DO I=2,NP+1
       IF(FD(I).GT.FDMAX) FDMAX=FD(I)
      END DO
      HALF=0.5*FDMAX
C     WRITE(*,*) 'HALF=0.5*FDMAX',HALF,FDMAX
*     find largest region less than half max
   30 NLAST=NP
      HALF=SQRT(4.0 *HALF)
C     WRITE(*,*) 'HALF=         ',HALF
      IA=0
      JA=0
      DO I=2,NP+2
       IF(I.GT.NP+1.OR.FD(I).GT.HALF) THEN
*         more than half - do not modify
          IF(IA.NE.0) THEN
C            WRITE(*,*) IA,IZ,' IS IA,IZ'
*            end of region IA...IB reached
             IF(JA.EQ.0.OR.IZ-IA.GT.JZ-JA) THEN
                JA=IA
                JZ=IZ
             END IF
             IA=0
          END IF
       ELSE
*         less than half
          IF(IA.EQ.0) THEN
*            start region (one earlier)
             IA=MAX(2,I-1)
             IZ=I
          ELSE
*            save bin index (one later)
             IZ=MIN(NP,I)+1
          END IF
       END IF
      END DO
C     WRITE(*,*) JA,JZ,JZ-JA+1,' bins'
      IF(JA.EQ.0) GOTO 100
      IF(JZ-JA.LE.1) GOTO 100
C     WRITE(*,*) JA,JZ,JZ-JA+1,' bins'
      IF(MOD(JZ-JA,2).EQ.0) THEN
*        increase by one
         IF(FD(JA).LT.FD(JZ)) THEN
            IF(JA.NE.   2) JA=JA-1
         ELSE
            IF(JZ.NE.NP+1) JZ=JZ+1
         END IF
      END IF
C     WRITE(*,*) JA,JZ,JZ-JA+1,' bins'
      IF(JZ-JA.LE.1) GOTO 100
      IF(MOD(JZ-JA,2).EQ.0) THEN
*        decrease by one
         IF(FD(JZ).LT.FD(JA)) THEN
            IF(JA.NE.   2) JA=JA+1
         ELSE
            IF(JZ.NE.NP+1) JZ=JZ-1
         END IF
      END IF
C     WRITE(*,*) JA,JZ,JZ-JA+1,' bins'
      IF(JZ-JA.LE.1) GOTO 100
      IF(MOD(JZ-JA,2).EQ.0) THEN
*        decrease by one
         IF(FD(JZ).LT.FD(JA)) THEN
            JA=JA+1
         ELSE
            JZ=JZ-1
         END IF
      END IF
C     WRITE(*,*) JA,JZ,JZ-JA+1,' bins'
*     combine bins
      DO J=JA,JZ,2
         FD(J)=SQRT(FD(J)**2+FD(J+1)**2)
         FD(J+1)=0.0
         IB(J)=IB(J)+IB(J+1)
         IB(J+1)=0
      END DO
C     WRITE(*,112) 'Array IB             ', (IB(L),L=1,NP+1)
*     compress arrays
      I=1
      DO J=2,NP+1
         IF(IB(J).NE.0) THEN
            I=I+1
            IB(I)=IB(J)
            FD(I)=FD(J)
         END IF
      END DO
      NP=I-1
C     WRITE(*,111) 'Array FD <<<<<<<<<<  ', (FD(L),L=1,NP+1)
C     WRITE(*,112) 'Array IB             ', (IB(L),L=1,NP+1)
C     HALF=HALF*(1.0+FLOAT(JZ-JA)/FLOAT(NLAST)*0.25)
      GOTO 30
  112 FORMAT(1X,A/(1X,20I5))
  100 RETURN
      END
      SUBROUTINE HCOMPX(H,N,B,NP)
*
*     compress histogram
*                    - -
*        CALL HCOMPX(H,N,B,NP)
*                        - --
*     where H(1:N)  = histogram content
*           B(1:NP) = bin compression
*
*
*
      REAL H(N)
      REAL YY(256),Y1(256),Y2(256),Y3(256),Y4(256)
      REAL B(*)
*     results from smoothing is used via common
      COMMON/CBSP4/NKNOT,T(64),AMP(2,64),ENSA,LNEG,SNOISE,A1,A2,
     +             FD(256),Z(256),Q(320)
*     ...
      NP=0
*     determine indices IA...IZ of nonzero histogram bins
      DO I=1,N
       IF(H(I).NE.0.0) GOTO 10
      END DO
      GOTO 100
   10 IA=I-1
      DO I=N,1,-1
       IF(H(I).NE.0.0) GOTO 20
      END DO
   20 IB=I
*     copy nonzero part of histogram... IA+1 ... IA+NB=IB
      NB=IB-IA
      DO I=1,NB
       YY(I)=H(IA+I)
      END DO
*     ... and smooth the histogram data
      CALL SPLFT(YY,NB)
*     determine first derivative ...
      DO I=1,NB-1
       Y1(I)=YY(I+1)-YY(I)
      END DO
      Y1(NB)=0.0
*     ... and second derivative ...
      DO I=2,NB-1
       Y2(I)=Y1(I)-Y1(I-1)
      END DO
      Y2(1)=Y2(2)
      Y2(NB)=Y2(NB-1)
*     ... and third derivative ...
      DO I=2,NB-2
       Y3(I)=Y2(I+1)-Y2(I)
      END DO
      Y3(1)=0.0
      Y3(NB-1)=0.0
      Y3(NB)=0.0
*     ... and fourth derivative
      DO I=3,NB-2
       Y4(I)=ABS(Y3(I)-Y3(I-1))
      END DO
      Y4(1)=0.5*(Y4(3)+Y4(4))
      Y4(2)=0.5*(Y4(3)+Y4(4))
      Y4(NB-1)=0.5*(Y4(NB-2)+Y4(NB-3))
      Y4(NB  )=0.5*(Y4(NB-2)+Y4(NB-3))
*     smoothing of fourth derivative by folding
      AV=0.25*(Y4(1)+Y4(2)+Y4(2)+Y4(3))
      DO I=2,NB-1
       VA=AV
       AV=0.25*(Y4(I-1)+Y4(I)+Y4(I)+Y4(I+1))
       Y4(I-1)=VA
      END DO
      Y4(NB-1)=AV
      Y4(NB  )=AV
*     second derivative
      DO I=1,NB
       Y2(I)=ABS(Y2(I))
      END DO
*     smoothing by folding
      AV=0.25*(Y2(1)+Y2(2)+Y2(2)+Y2(3))
      DO I=2,NB-1
       VA=AV
       AV=0.25*(Y2(I-1)+Y2(I)+Y2(I)+Y2(I+1))
       Y2(I-1)=VA
      END DO
      Y2(NB-1)=AV
      Y2(NB  )=AV
*     FD from SPLFT
      DO I=NB,2,-1
       FD(I)=FD(I)-FD(I-1)
      END DO
*
      SA=AMAX1(1.0,ENSA)
      NP=SQRT(2.0*FLOAT(NB)*SA)+0.5
      IF(MOD(NB-NP,2).EQ.1) NP=NP+1
*     smooth FD = fourth derivative
      AV=0.25*(FD(1)+FD(2)+FD(2)+FD(3))
      DO I=2,NB-1
       VA=AV
       AV=0.25*(FD(I-1)+FD(I)+FD(I)+FD(I+1))
       FD(I-1)=VA
      END DO
      FD(NB-1)=AV
      FD(NB  )=AV
*     integrate FD ...
      DO I=2,NB
       FD(I)=FD(I)+FD(I-1)
      END DO
*     ... adding constant to each value
      CSTS=0.0
      DO I=1,NB
       CSTS=CSTS+20.0
       Z(I)=FD(I)+CSTS
      END DO
*     define "constant step" subdivision T(.)
      DO I=0,NP
       YF=FLOAT(I)/FLOAT(NP)*Z(NB)
*      interpolate - find index for YF in array Z(1)...Z(NB)
       L=LEFTF(YF,Z,NB)
       IF(L.EQ.0) THEN
          YL=0.0
       ELSE
          YL=Z(L)
       END IF
*      determine XF by interpolation
       XF=FLOAT(L)+(YF-YL)/(Z(L+1)-YL)
       T(I+1)=XF
      END DO
*
      ZS=0.0
      DO I=1,NP
*      choose 1 or 3 or 5 or 7 bins
       DZ=1.0
       DO J=1,3
        DT=T(I+1)-ZS-DZ
        IF(DT.LT.1.0) GOTO 26
        DZ=DZ+2.0
       END DO
*      use DZ bins
   26  ZS=ZS+DZ
       Z(I)=DZ
      END DO
*     determine final result by sum
      B(1)=IA
      DO I=1,NP
       B(I+1)=B(I)+Z(I)
      END DO
  100 RETURN
      END
      SUBROUTINE HPRT(H,N,A,E)
*
*     print histogram  =  content  of  array  h(1)...h(n)  with  several
*
*     SUBROUTINE CALLED - PVERT
*
      CHARACTER*128 PX
      COMMON/CCONVT/PX(10)
      DOUBLE PRECISION GSUM,SUM
      REAL H(1),G(120),Q(16),SP(13)
      CHARACTER*10 RPR(13),STR*12
*     ...
      IC=ICHAR('0')-1
*     eventually combine bins together
      MM=1+(N-1)/120
      M =N/MM
      LM=N/MM
*
*     sum, integrate content, find limits etc
*
*     COPY h TO g
      K=0
      DO I=1,N,MM
       JEND=MIN0(N,I+MM-1)
       SUM=0.0
       DO J=I,JEND
        SUM=SUM+H(J)
       END DO
       K=K+1
       G(K)=SUM
      END DO
*     find LM, IMIN, IMAX, NPOS, NNEG ...
      LM=0
      IMIN=1
      IMAX=1
      NPOS=0
      NNEG=0
      GSUM=0.0
      DO I=1,M
       IF(G(I).NE.0.0) THEN
          GSUM=GSUM+G(I)
          LM=I
          IF(G(I).GT.0.0) THEN
             NPOS=NPOS+1
             IF(G(I).GT.G(IMAX)) IMAX=I
          ELSE
             NNEG=NNEG+1
             IF(G(I).LT.G(IMIN)) IMIN=I
          END IF
       END IF
      END DO
      WRITE(*,102)
      WRITE(*,101) GSUM,N,A,E
      LASTM=LM
*
      IF(LM.EQ.0) GOTO 100
*     print up to 10 lines of X's
      DO J=1,10
       PX(J)=' '
      END DO
      IF(NNEG*10.LT.LM) THEN
*        PLOT ONLY POSITIVE CONTENTS
         FAC=10.0/G(IMAX)
         DO I=1,LM
          K=G(I)*FAC+0.5
          IF(K.GT.0) THEN
             JA=11-K
             DO J=JA,10
              PX(J)(I+8:I+8)='X'
             END DO
          END IF
         END DO
      ELSE
*        print positive and negative contents
         GMIN=AMIN1(G(IMIN),0.0)
         GMAX=AMAX1(G(IMAX),0.0)
         FAC=9.0/(GMAX-GMIN)
   37    LMIN=-GMIN*FAC+0.5
         LMAX=+GMAX*FAC+0.5
         IF(LMAX+LMIN.GT.9) THEN
            FAC=FAC*0.98
            GOTO 37
         END IF
         LMIN=-LMIN
         DO I=1,LM
         IF(G(I).NE.0.0) THEN
             IF(G(I).GT.0.0) THEN
                L=+G(I)*FAC+0.5
                LA=LMAX+1-L
                LB=LMAX
                LL=LMAX+1
             ELSE
                L=-G(I)*FAC+0.5
                L=-L
                LA=LMAX+2
                LB=LMAX+1-L
                LL=LMAX+1
             END IF
             DO L=LA,LB
              PX(L)(I+8:I+8)='X'
             END DO
             PX(LL)(I+8:I+8)='0'
          END IF
         END DO
      END IF
      WRITE(*,102)
      DO J=1,10
       WRITE(*,102)  PX(J)
      END DO
*     print bin content vertical
      CALL PVERT(G(1),LM,-5)
*     print numbered line
      ML=((LM+9)/10)*10
      DO I=1,2
       PX(I)=' '
      END DO
      DO J=1,ML
       PX(1)(8+J:8+J)='-'
       IF(MOD(J,2).NE.1) THEN
          I=MOD(J,10)+1
          PX(2)(8+J:8+J)=CHAR(I+IC)
          IF(I.EQ.1) THEN
             I=MOD(J/10,10)+1
             PX(1)(8+J:8+J)=CHAR(I+IC)
          END IF
       END IF
      END DO
      DO J=1,2
       WRITE(*,102)  PX(J)
      END DO
*     print scale
      ML=((LM+9)/10)*10
      LN=ML/10+1
      ST=(E-A)/FLOAT(M)
      IF(LN.GT.12) LN=12
      DO I=1,LN
       SP(I)=A+ST*10.0*FLOAT(I-1)
       IF(ABS(SP(I)).LT.5.0E-4*ST) SP(I)=0.0
       CALL PNVG(SP(I),STR,JS)
*      COPY STRING
       RPR(I)=' '
       IF(JS.LE.10) THEN
          JT=(10-JS)/2
          RPR(I)(JT+1:JT+JS)=STR(1:JS)
       END IF
      END DO
      WRITE(*,103)    (RPR(I),I=1,LN)
  100 RETURN
  101 FORMAT(' Histogram content',G12.4,' in',I4,' bins between',
     1   G12.4,' and',G12.4)
  102 FORMAT(3X,A128)
  103 FORMAT(6X,12A10)
      END


      SUBROUTINE CNVF(FLT,STR,JS)
*
      CHARACTER*(*) STR
*
************************************************************************
*
*     Convert floating point number FLT into string STR with JS
*     non-blank characters. JS = 2 ... 12
*
*     Examples:
*     CALL CNVF( -1.200,STR,JS)   -> STR='-1.2' ; JS = 4
*     CALL CNVF(1.0/3.0,STR,JS)   -> STR='0.333333' ; JS = 8
*     CALL CNVF(-123.00,STR,JS)   -> STR='-123.' ; JS = 5
*     STR will contain up to 6 digits, trailing zeros are suppressed.
*     STR will contain one decimal point.
*     Range of floating point numbers is 10E-36 ... 10E+36 (due to
*     limitations on certain machines.
*
*
************************************************************************
*
      CHARACTER*1 DIG(0:9)
      INTEGER ND(12),KD(6)
      DOUBLE PRECISION FEX(-3:5)
*     ...
      DATA DIG/'0','1','2','3','4','5','6','7','8','9'/
      DATA FEX/1.0D8,1.0D7,1.0D6,1.0D5,1.0D4,1.0D3,1.0D2,1.0D1,1.0D0/
      DATA KD/100000,10000,1000,100,10,1/
      DATA IJ0/0/
*
      IF(IJ0.EQ.0) THEN
         IJ0=ICHAR('0')
         IJ9=ICHAR('9')
      END IF
      JS=0
      X=FLT
      IF(X.NE.0.0) THEN
*        limited range on certain machines
         IF(ABS(X).LT.1.0E-36) X=SIGN(1.0E-36,FLT)
         IF(ABS(X).GT.1.0E+36) X=SIGN(1.0E+36,FLT)
      ELSE
         GOTO 22
      END IF
      N=IFIX(ALOG10(ABS(X))+100.0)-100
   20 IF(N.LE.5.AND.N+3.GE.0) THEN
         K=DBLE(ABS(X))*FEX(N)+0.5D0
      ELSE
         K=DBLE(ABS(X))*10.0D0**(5-N)+0.5D0
      END IF
   21 IF(K.GT.999999) THEN
         N=N+1
         GOTO 20
      END IF
*     rounding
      IF(MOD(K,100).EQ.99) THEN
         K=K+1
         GOTO 21
      ELSE IF(MOD(K,1000).EQ.998) THEN
         K=K+2
         GOTO 21
      ELSE IF(MOD(K,100).EQ.01) THEN
         K=K-1
      ELSE IF(MOD(K,1000).EQ.002) THEN
         K=K-2
      END IF
*     determine number of digits
      KDUP=K
      DO 40 I=1,6
      IF(MOD(KDUP,10).NE.0) GOTO 50
      KDUP=KDUP/10
   40 CONTINUE
      I=6
   50 NDIG=7-I
*
   22 IF(X.EQ.0.0) THEN
         GOTO 23
      ELSE IF(N.GT.5) THEN
*        N greater than +5
         NFL=(N/3)*3
         M=MOD(N,3)+1
      ELSE IF(N+3.GE.0) THEN
*        N between +5 and -3
         NFL=0
         M=MAX0(0,N+1)
      ELSE
*        N less than -3
         NFL=-((2-N)/3)*3
         M=3-MOD(2-N,3)
*        calculate length, if exponent modified by 3
         NLG=MAX0(M-3,NDIG)+4
         IF(X.LT.0.0)     NLG=NLG+1
         IF(M.LT.3)       NLG=NLG+1
         IF(M.LT.2)       NLG=NLG+1
         IF(NFL.LT.(-12)) NLG=NLG+1
*        if not more than 12 bytes, add 3 to exponent
         IF(NLG.LE.12) THEN
            NFL=NFL+3
            M  =M-3
            N  =N+3
         END IF
      END IF
*
*     K = 6-digit number
*     N = exponent (or +3)
*     NFL = printed exponent
*     M   = position of point
*     NDIG = number of nonzero digits
*
   23 STR=' '
*
      IF(X.EQ.0.0) THEN
         JS=JS+1
         STR(JS:JS+1)='0.'
         JS=JS+1
         GOTO 100
      ELSE IF(X.LT.0) THEN
         JS=JS+1
         STR(JS:JS)='-'
      END IF
*
      IF(NFL.GE.0.AND.(N+1.LE.0.AND.N+3.GE.0)) THEN
*         special case for N= -1 or -2 or -3 (0. or 0.0 or 0.00)
         JS=JS+1
         STR(JS:JS-N+1)='0.00'
         JS=JS-N
*        added code to avois some problems
      ELSE IF(NFL.LT.0.AND.M.LE.0) THEN
         NFL=NFL-3
         M=M+3
      END IF
*
      DO 24 I=1,6
      J=K/KD(I)
      K=K-KD(I)*J
      JS=JS+1
      STR(JS:JS)=DIG(J)
      IF(I.EQ.M) THEN
         JS=JS+1
         STR(JS:JS)='.'
      END IF
      IF(I.GE.M.AND.K.EQ.0) GOTO 26
   24 CONTINUE
*
   26 IF(NFL.EQ.0) GOTO 100
      JS=JS+1
      STR(JS:JS)='E'
*
      K=IABS(NFL)
      DO 30 I=1,2
      ND(I)=MOD(K,10)
      K    =K/10
      IF(K.EQ.0) GOTO 32
   30 CONTINUE
      I=2
   32 IF(NFL.LT.0) THEN
         JS=JS+1
         STR(JS:JS)='-'
      END IF
      DO 34 J=I,1,-1
      JS=JS+1
   34 STR(JS:JS)=DIG(ND(J))
*
  100 RETURN
      END

      SUBROUTINE CNVG(FLT,STR,JS)
*
      CHARACTER*(*) STR
*
************************************************************************
*
*     The same as CNVF, but for graphics.
*
************************************************************************
*
      CHARACTER*12 TTR
      IF(FLT.EQ.0.0) THEN
         STR='0'
         JS=1
         GOTO 100
      END IF
      CALL PNVF(FLT,STR,JS)
      IE=INDEX(STR(1:JS),'E')
      IF(IE.EQ.0) THEN
         IF(STR(JS:JS).EQ.'.') THEN
            STR(JS:JS)=' '
            JS=JS-1
         END IF
      ELSE
         IF(STR(IE-1:IE-1).EQ.'.') THEN
            NE=1+JS-IE
            TTR(1:NE)=STR(IE:JS)
            STR(IE-1:IE-2+NE)=TTR(1:NE)
            JS=JS-1
         END IF
      END IF
*
  100 RETURN
      END

      SUBROUTINE AHISTS(A,N,XA,XB)
      REAL A(*)
*
*     prepare array A(1)...A(N+10) for histogram with N bins,
*     region XA ... XB
*
*     A(1)    dimension of array
*     A(2)    XA
*     A(3)    XB
*     A(4)    NB (floating point)
*     A(5)    total nr of entries
*     A(6)    nr of entries low
*     A(7)    nr of entries high
*     A(8)    total sum of weights
*     A(9)    outside low
*     A(10)   outside high
*     A(11)...A(N+10) histogram
*
      NB=N
      DO I=1,N+10
       A(I)=0.0
      END DO
      A(1)=N+10
      A(2)=XA
      A(3)=XB
      A(4)=NB
      RETURN
      END

      SUBROUTINE AHISTU(A,X)
*     entry into histogram
      REAL A(*)
*     ...
      I=(X-A(2))/(A(3)-A(2))*A(4)+1.0
      NB=A(4)
      IF(I.LE. 0) I=-1
      IF(I.GT.NB) I= 0
      A(5)=A(5)+1.0
      A(8)=A(8)+1.0
      IF(I.LE.0) A(7+I)=A(7+I)+1.0
      A(10+I)=A(10+I)+1.0
      RETURN
      END

      SUBROUTINE AHISTW(A,X,W)
*     entry into histogram
      REAL A(*)
*     ...
      I=(X-A(2))/(A(3)-A(2))*A(4)+1.0
      NB=A(4)
      IF(I.LE. 0) I=-1
      IF(I.GT.NB) I= 0
      A(5)=A(5)+1.0
      A(8)=A(8)+W
      IF(I.LE.0) A(7+I)=A(7+I)+W
      A(10+I)=A(10+I)+W
      RETURN
      END

      SUBROUTINE AHISTP(A)
*     print histogram
      REAL A(*)
*     ...
      NB=A(4)
      WRITE(6,101) (A(I),I=5,10)
      CALL VPRINT(A(11),NB,A(2),A(3),' in AHISTP')  
C     CALL HPRT(A(11),NB,A(2),A(3))
  100 RETURN
  101 FORMAT(' A-Histogram  entries total/low/high=',3F9.0,
     +                 ' W-sum total/low/high=',3G12.4)
      END

      SUBROUTINE ROUNDB(XA,XB,NB,A,B)
*
*     Determine rounded limits A and B of a region with NB bins, which
*     includes the values XA and XB.
*
*                 -- -- --
*     CALL ROUNDB(XA,XB,NB,A,B)
*                          - -
*
*     Example: CALL ROUNDB(0.0172,2.5726,60,A,B)
*              resulting A=0.0  B=3.0
*
      REAL BINS(10)
      DATA BINS/0.1,0.15,0.2,0.25,0.4,0.5,0.6,0.8,1.0,1.5/
*     ...
      ZA=XA
      ZB=XB
      BIN=(ZB-ZA)/FLOAT(NB)
      IF(BIN.LE.0.0) BIN=1.0E-4
      NEX=INT(ALOG10(BIN)+100.0)-99
      FEX=10.0**NEX
      BIN=BIN/FEX
      DO I=1,9
       IF(BINS(I).GE.BIN) GOTO 10
      END DO
      I=9
   10 BIN=BINS(I)*FEX
      A=10.0*BIN*AINT(ZA/(10.0*BIN))
      IF(A.GT.ZA) A=A-10.0*BIN
      B=A+FLOAT(NB)*BIN
      IF(B.LT.ZB.AND.I.NE.10) THEN
         I=I+1
         GOTO 10
      END IF
      END

      SUBROUTINE BGTEXT(TEXT)

*     PRINTS UP TO 18 CHARACTERS WITH LARGE (7 LINES) CHARACTERS
*                 ----
*     CALL BGTEXT(TEXT)
*
*     TEXT = UP TO 18 CHARACTERS (> 1 CHARACTER)
*
*     IF TEXT IS ONLY ONE CHARACTER, THIS CHARACTER IS USED TO PRINT
*     THE LARGE CHARACTERS (DEFAULT IS THE CHARCATER O) IN FOLLOWING
*     CALLS WITH A ARGUMENT TEXT, LONGER THAN 1 CHARACTER
*
*
*     TO INSERT (UP TO 5) INTEGER NUMBERS I1, I2 ... INTO THE TEXT,
*     CALL BGINT  WITH THE INTEGER NUMBERS AS ARGUMENTS, AND THEN
*     BGTEXT WITH THE TEXT AS ARGUMENT, WHERE TEXT HAS GROUPS OF &
*     AT THE POSITIONS OF THE INTEGERS, E.G.
*
*     CALL BGINT(15)
*     CALL BGINT( 3)
*     CALL BGINT( 5)
*     CALL BGTEXT('&&& IS && TIMES &')
*
*     SPECIAL CALL
*
*     CALL BGTEXT('PAGE')   STARTS NEW PAGE
*
*
*
      PARAMETER (LUNPR=6)
      CHARACTER*(*) TEXT
      CHARACTER*18  TLINE
      CHARACTER*1   CHDIS/'O'/
      CHARACTER*126 PR
      CHARACTER*55 CH
     1  /'ABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789+-''/()$=,.#*<>%?!:0'/
      INTEGER ND/0/,IN(5),LUN/-1/
      INTEGER LINE(126),POWER(7)/64,32,16,8,4,2,1/

*     INTEGER CODE FOR BIG CHARACTERS

      INTEGER ICH(5,55)
     *      /63,68,68,68,63,127,73,73,73,54,62,65,65,65,34,127,65,65,65,
     *62,127,73,73,73,65,127,72,72,72,64,62,65,73,73,46,127,8,8,8,127,0,
     *0,127,0,0,2,1,1,1,126,127,8,20,34,65,127,1,1,1,1,127,32,16,32,127,
     *127,16,8,4,127,62,65,65,65,62,127,72,72,72,48,62,65,69,67,63,127,7
     *2,76,74,49,50,73,73,73,38,64,64,127,64,64,126,1,1,1,126,124,2,1,2,
     *124,127,2,4,2,127,99,20,8,20,99,96,16,15,16,96,67,69,73,81,97,62,6
     *5,65,65,62,0,16,32,127,0,33,67,69,73,49,34,65,73,73,54,6,10,23,34,
     *66,114,81,81,81,78,62,73,73,73,38,65,66,68,72,112,54,73,73,73,54,5
     *0,73,73,73,62,0,8,28,8,0,0,8,8,8,0,0,0,112,0,0,2,4,8,16,32,0,0,62,
     *65,0,0,65,62,0,0,48,73,127,73,6,0,20,20,20,0,0,0,7,0,0,0,0,1,0,0,
     *20,62,20,62,20,42,28,62,28,42,0,8,20,34,0,0,34,20,8,0,
     *114,84,127,21,39,32,64,77,72,48,0,0,125,0,0,0,0,20,0,0,
     *0,0,0,0,0/
      COMMON/BCS/IW(30)
*     ...
      IF(LUN.LT.0) THEN
         LUN=LUNPR
         IF(IW(30).EQ.12345) LUN=IW(6)
      END IF
      IF(LUN.LE.0) GOTO 100
*     SKIP TO NEW PAGE
      IF(TEXT.EQ.'PAGE') THEN
         WRITE(LUN,101)
         GOTO 100
      END IF
*     DEFINE DISPLAY CHARACTER
      IF(LEN(TEXT).EQ.1) THEN
         CHDIS=TEXT(1:1)
         GOTO 100
      END IF

*     INSERT INTEGERS INTO THE TEXT

      TLINE=TEXT
      IC=ICHAR('0')
      DO 60 I=1,ND

      DO 10 J=1,18
      IF(TLINE(J:J).EQ.'&') GOTO 20
   10 CONTINUE
      GOTO 70
   20 LA=J
      DO 30 J=LA,18
      IF(TLINE(J:J).NE.'&') GOTO 40
   30 CONTINUE
      J=19
   40 LB=J-1
      TLINE(LA:LB)=' '
      IJK=IABS(IN(I))
      DO 50 L=LA,LB
      IF(IJK.NE.0.OR.L.EQ.LA) THEN
          LC=MOD(IJK,10)
          IJK=IJK/10
          TLINE(LA+LB-L:LA+LB-L)=CHAR(IC+LC)
      END IF
   50 CONTINUE
      IF(IN(I).LT.0) TLINE(LA:LA)='-'
   60 CONTINUE
   70 ND=0

*     PRINT IN SEVEN LINES

      WRITE(LUN,102)

*     TRANSLATE CHARACTERS TO INTEGER CODE

      L=1
      DO 80  J=1,18
      CH(55:55)=TLINE(J:J)
      I=INDEX(CH,TLINE(J:J))

*     TWO BLANK COLUMNS

      LINE(L  )=0
      LINE(L+1)=0
      L=L+2

*     INSERT INTEGER CODE FOR CHARACTER IN 5 COLUMNS

      DO 75 K=1,5
      LINE(L)=ICH(K,I)
   75 L=L+1
   80 CONTINUE

*     PREPARE LINE AFTER LINE AND PRINT

      DO 90 K=1,7
      PR=' '
      DO 85 L=1,126
      IF(IAND(POWER(K),LINE(L)).NE.0) THEN
         PR(L:L)=CHDIS
      END IF
   85 CONTINUE
      WRITE(LUN,102) PR
   90 CONTINUE
      WRITE(LUN,102)
      GOTO 100


      ENTRY BGINT(NUMBER)

*     STORE INTEGER FOR NEXT TEXT PRINTOUT

      IF(ND.GE.5) GOTO 100
      ND=ND+1
      IN(ND)=NUMBER

  100 RETURN
  101 FORMAT(1H1)
  102 FORMAT(1X,A)
      END


      FUNCTION LEFTF(X,T,N)
*
*     Search and determine index.
*
*     function LEFTF returns the value L, which satisfies
*                   T(L)  LE  X  LT  T(L+1)
*     except in the case X=T(N), where
*                   T(L)  LT  X  LE  T(L+1)
*     holds.
*     T(1),T(2)...T(n) is a non-decreasing  sequence.  In  the  case  of
*     multiple T-values (e.g. T(2)=T(3)) L satisfies
*                       T(L)  LT  T(L+1).
*     L = 0  for X  LT  T(1)
*     L = N  for X  GT  T(N)
*
*               - - -
*        L=LEFTF(X,T,N)
*
*     Example:  X(1)=0.2  X(2)=1.7  X(3)=3.4  X(4)=7.9
*
*               L=LEFTF( 1.98,X,4)  =>   L=2
*               L=LEFTF(12.86,X,4)  =>   L=4
*
      REAL T(*)
*     handle x = outside cases
      IF(X.LT.T(1)) THEN
*        outside low
         LEFTF=0
         GOTO 100
      END IF
      IF(X.GE.T(N)) THEN
*        outside high
         LEFTF=N
   10    IF(X.GT.T(LEFTF)) GOTO 100
*        X on upper limit
         LEFTF=LEFTF-1
         GOTO 10
      END IF
*     binary search
      IL=1
      IR=N
   20 LEFTF=(IL+IR)/2
      IF(X-T(LEFTF))      45,70,50
   45 IR=LEFTF   
      GOTO 20
   50 IF(T(LEFTF+1)-X)    55,60,70
   55 IL=LEFTF+1
      GOTO 20
*     final equality checks
   60 LEFTF=LEFTF+1
   70 IF(X.EQ.T(LEFTF+1)) GOTO 60
  100 RETURN
      END