C   07/06/96 606071908  MEMBER NAME  SMHIST   (S4)          FORTG1
      SUBROUTINE SMHIST(YI,Y,DY,N,IDY)
      REAL YI(1),Y(1),DY(1)
      REAL A(64,5),YL(128),YS(128),X(128),C(129),W(128)
C
C     ------------------------------------------------------------------
C
C
C     PURPOSE - SMOOTH HISTOGRAM DATA
C
C     USAGE                 --   -- - ---
C               CALL SMHIST(YI,Y,DY,N,IDY)
C                              -
C               YI( ) = INPUT HISTOGRAM OF N BINS
C               Y ( ) = SMOOTHED HISTOGRAM
C               DY( ) = ERRORS OF INPUT HISTOGRAM DATA (FOR IDY=1 ONLY)
C               IDY   = 0   HISTOGRAM WITH ERROR OF BIN CONTENT EQUAL
C                           TO SQRT(CONTENT), POISSON DISTR., OR
C                     = 1   WEIGTHED HISTOGRAM WITH ERRORS GIVEN.
C
C               YI( ) AND Y( ) MAY BE THE SAME ARRAY. THERE IS NO
C               LIMITATION ON THE NUMBER OF BINS.
C
C     METHOD  - IT IS ASSUMED, THAT THE HISTOGRAM DATA CORRESPOND TO
C               BINS OF EQUAL SIZE. THE METHOD CONSISTS OF SEVERAL
C     PARTS. PEAKS AND VALLEYS OF ONLY ONE BIN WILL COMPLETELY VANISH,
C     THUS THE BIN SIZE SHOULD BE CHOOSEN SMALL ENOUGH, THAT EXPECTED
C     PEAKS WILL BE AT LEAST 3 BINS WIDE.
C
C     1. PART      IF THREE ADJACENT BINS CONTAIN LESS THAN 8 (EQUIV.)
C                  ENTRIES, THEY WILL BE AVERAGED.
C     2. PART *    TRANSFORMATION OF THE POISSON DISTRIBUTED DATA TO
C                  DATA OF VARIANCE 1.0
C     3. PART      353C ALGORITHM TO REDUCE LOCAL CHI SQUARE
C                  (SEE SUBROUTINE SMOPER)
C     4. PART      SPLINE FIT (SEE BELOW)
C     5. PART *    BACK TRANSFORMATION
C
C     *  FOR IDY = 0 ONLY.
C
C     SPLINE FIT - THE SUBROUTINE SPLESQ IS USED TO PERFORM A SPLINE
C     FIT, WHICH CONSISTS IN A LEAST SQUARE OF A CERTAIN NUMBER OF
C     REFERENCE POINTS, CONNECTED BY A SPLINE FUNCTION. THE BAND MATRIX
C     OCCURING IN THE FIT IS SOLVED BY A FAST RTR DECOMPOSITION METHOD
C     IN SUBROUTINE BANDSL. THE SPLINE FIT IS PERFORMED ITERATIVELY,
C     REFINING THE SETTING OF REFERENCE POINTS ACCORDING TO THE LOCAL
C     CHI SQUARE. IN THE FIRST ITERATION, THE REFERENCE POINTS ARE
C     CHOOSEN EQUIDISTANT WITH ABOUT 8 BINS BETWEEN 2 REFERENCE POINTS.
C     IN THE FOLLOWING ITERATIONS, USING THE LOCAL CHI SQUARE, THE
C     DISTANCE BETWEEN 2 REFERENCE POINTS IS REDUCED.
C
C     THE RESULTING SMOOTHED HISTOGRAM CONTENT WILL BE SMOOTH WITH
C     LITTLE SHORT LENGTH OSCILLATIONS, WHILE SIGNIFICANT STRUCTURES
C     WILL BE PRESENT.
C
C     SUBROUTINES CALLED - SMOPER
C                          SPLESQ (CALLS BANDSL)
C
C
C     ------------------------------------------------------------------
C
C
C     TRANSFER TO Y( ) AND GET FIRST AND LAST BIN WITH NONZERO CONTENT
C
      IA=0
      IB=0
      WL=0.0
      IL=0
      DO 2 I=1,N
      Y(I)=YI(I)
      IF(IDY.EQ.0.AND.Y(I).EQ.0.0) GOTO 2
      IF(IDY.NE.0.AND.DY(I).EQ.0.0) GOTO 2
      IF(IDY.NE.0.AND.IA.EQ.0) WL=DY(I)
      IF(IA.EQ.0) IA=I
      IB=I
      IF(IDY.EQ.0.OR.I.EQ.1.OR.I.EQ.N) GOTO 2
      IF(DY(I-1).NE.0.0.OR.DY(I+1).NE.0.0) GOTO 2
      WL=WL+DY(I)**2
      IL=IL+1
    2 CONTINUE
      IF(IDY.NE.0.AND.IL.NE.0.AND.WL.EQ.0.0) GOTO 100
      IF(IDY.NE.0) WL=FLOAT(IL)/WL


      NC=IB-IA+1
      IF(NC.LE.3) GOTO 100
C
C     FOLD BINS OF VERY LOW STATISTIC
C
      NM=IB-2
      ZL=Y(IA)
      DO 6 I=IA,NM
      ZN=Y(I+1)
      SUM=Y(I)+Y(I+1)+Y(I+2)
      IF(IDY.EQ.0) GOTO 3
      SUM2=DY(I)**2+DY(I+1)**2+DY(I+2)**2
      IF(SUM2.EQ.0.0) GOTO 4
      IF(SUM**2/SUM2.GT.8.0) GOTO 4
      ZN=0.25*(ZN+Y(I+1))
      GOTO 4
    3 IF(SUM.GT.8.0) GOTO 4
      ZN=ZN+(SUM-3.0*ZN)/(SUM+3.0)
    4 Y(I)=ZL
    6 ZL=ZN
      Y(IB-1)=ZL
C
C     TRANSFORM POISSON DISTRIBUTED DATA TO DATA OF VARIANCES 1.0
C
      IF(IDY.NE.0) GOTO 12
      DO 10 I=IA,IB
      IF(Y(I).GT.0.0) GOTO 8
      Y(I)=1.0
      GOTO 10
    8 Y(I)=SQRT(Y(I))+SQRT(Y(I)+1.0)
   10 CONTINUE
C
C     APPLY SIMPLE 353 ALGORITHM TO REDUCE CHI SQUARES
C
   12 CALL SMOPER(Y(IA),NC,321)
C
C     SMOOTH BY REPEATED SPLINE FITS
C
      IF(IDY.EQ.0) IEQ=-1
      IF(IDY.NE.0) IEQ=+1
      KL=(NC+101)/112
      IF(KL.LE.0) KL=1
      LK=(NC+12*KL-11)/KL
      JA=IA
      DO 80 LL=1,KL
      JB=MIN0(IB,JA+LK-1)
      I=0
      DO 20 J=JA,JB
      I=I+1
      IF(LL.EQ.1) YL(I)=Y(J)
      YS(I)=YL(I)
      IF(IDY.EQ.0) GOTO 20
      W(I)=WL
      IF(DY(J).NE.0.0) W(I)=1.0/DY(J)**2
   20 X(I)=0.25
      X(1)=0.0
      M=I
      NSP=3+MAX0(M/8,1)
      CALL SPLESQ(YS,X,W,M,A,NSP,IEQ,CHDF)
      CLIM=1.0
C     DO 70 K=1,2
      DO 70 K=1,4
      NBEL=0
      NVLG=0
      C(1)=0.0
      DO 30 I=1,M
      IF(IDY.EQ.0) C(I+1)=C(I)+(YS(I)-YL(I))**2
      IF(IDY.NE.0) C(I+1)=C(I)+W(I)*(YS(I)-YL(I))**2
      IF(C(I+1)-C(I).LT.CLIM) NBEL=NBEL+1
   30 YS(I)=YL(I)
      DC=C(M+1)/FLOAT(M-NSP)
      WRITE(6,101) K,DC,DVMAX,M,NSP,IDY,NBEL
  101 FORMAT(1X,I2,2G15.5,4I15)
      DV=C(4)/3.0
      CFAC=1.0
      IF(NBEL+NBEL.GT.M) CFAC=1.1665
      FAC=CFAC
C     IF(DV+DV.GT.CLIM) FAC=1.4142
C     IF(DV+DV.GT.CLIM) FAC=1.2599
C     IF(   DV.GT.CLIM) FAC=2.0000
C     IF(   DV.GT.CLIM) FAC=1.5874
      IF(DV.LT.0.2*DC) FAC=FAC/1.1665
      IF(DV.GT.2.0*DC) FAC=FAC*1.1665
      IF(DV.GT.3.0*DC) FAC=FAC*1.1665
      IF(DV.GT.3.0*DC) NVLG=NVLG+1
      X(2)=FAC*X(2)
      X(3)=FAC*X(3)
      DV=(C(M+1)-C(M-1))/3.0
      FAC=CFAC
      IF(DV.LT.0.2*DC) FAC=FAC/1.1665
C     IF(DV+DV.GT.CLIM) FAC=1.4142
C     IF(DV+DV.GT.CLIM) FAC=1.2599
C     IF(   DV.GT.CLIM) FAC=2.0000
C     IF(   DV.GT.CLIM) FAC=1.5874
      IF(DV.GT.2.0*DC) FAC=FAC*1.1665
      IF(DV.GT.3.0*DC) FAC=FAC*1.1665
      IF(DV.GT.3.0*DC) NVLG=NVLG+1
      X(M)=FAC*X(M)
      X(M-1)=FAC*X(M-1)
      N3=M-3
      IF(N3.LT.4) GOTO 45
      DVMAX=0.0
      DO 40 I=4,N3
      DV=0.125*(C(I+3)+C(I+1)-C(I-1)-C(I-3))
      DVMAX=AMAX1(DV,DVMAX)
      X(I)=X(I)*CFAC
C     IF(DV+DV.GT.CLIM) X(I)=X(I)*1.4142
C     IF(DV+DV.GT.CLIM) X(I)=X(I)*1.2599
C     IF(   DV.GT.CLIM) X(I)=X(I)*1.4142
C     IF(   DV.GT.CLIM) X(I)=X(I)*1.2599
      IF(DV.LT.0.2*DC) X(I)=X(I)/1.1665
      IF(DV.GT.2.0*DC) X(I)=X(I)*1.1665
      IF(DV.GT.3.0*DC) X(I)=X(I)*1.1665
      IF(DV.GT.3.0*DC) NVLG=NVLG+1
   40 CONTINUE
   45 CALL SMOPER(X(2),M-1,5321)
      DO 50 I=2,M
   50 X(I)=X(I-1)+X(I)
      NSP=3+X(M)*0.5
      NSP=MIN0(NSP,64)
      CALL SPLESQ(YS,X,W,M,A,NSP,2*IEQ,CHDF)
      IF(CFAC.EQ.1.0.AND.NVLG.EQ.0) GOTO 71
      XL=X(1)
      DO 60 I=2,M
      XN=X(I)
      X(I)=XN-XL
   60 XL=XN
C     CLIM=0.25*CLIM
C     CLIM=0.50*CLIM
   70 CONTINUE
   71 IF(LL.EQ.KL) GOTO 74
      KA=JB-10
      KB=MIN0(IB,KA+LK-1)
      I=0
      DO 72 K=KA,KB
      I=I+1
   72 YL(I)=Y(K)
   74 I=0
      DO 78 J=JA,JB
      I=I+1
      IF(I.GT.11.OR.JA.EQ.IA) GOTO 76
      Y(J)=0.1*(FLOAT(I-1)*YS(I)+FLOAT(11-I)*Y(J))
      GOTO 78
   76 Y(J)=YS(I)
   78 CONTINUE
   80 JA=JB-10
C
C     TRANSFORM BACK TO POISSON DISTRIBUTED DATA
C
      IF(IDY.NE.0) GOTO 100
      DO 90 I=1,N
      IF(Y(I).GT.1.0) GOTO 85
      Y(I)=0.0
      GOTO 90
   85 SY=(0.5*(Y(I)-1.0/Y(I)))**2
      NY=100.0*SY+0.5
      Y(I)=0.01*FLOAT(NY)
   90 CONTINUE
C
C
  100 RETURN
      END
