C   07/06/96 606071908  MEMBER NAME  SMHIST1  (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(128),W(128),DX(128)
      INTEGER NR/0/
      NR=NR+1
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,5321)
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
      CLIM=0.20
      DO 22 I=1,M
      C(I)=YL(I)**2
      IF(IDY.NE.0) C(I)=C(I)*W(I)
   22 CONTINUE
      CALL SMOPER(C,M,7)
      SUM=0.0
      DO 24 I=1,M
   24 SUM=SUM+C(I)
      SUM=SUM/FLOAT(M)
      IF(IDY.EQ.0) SUM=0.25*SUM
      DXA=2.1/FLOAT(M)
      DXB=0.125*SQRT(0.02*SUM+0.1)
      DXB=AMIN1(0.125,DXB)
      DXB=AMAX1(DXA,DXB)
      DXQ=0.5*(DXA+DXB)
      DO 26 I=1,M
      X(I)=FLOAT(I-1)*DXQ*2.0
      FAC=SQRT(C(I)/SUM)
      FAC=AMAX1(0.3,FAC)
      FAC=AMIN1(2.0,FAC)
      DX(I)=FAC*DXQ
      DX(I)=AMAX1(DXA,DX(I))
   26 DX(I)=AMIN1(0.5,DX(I))
      DO 70 K=1,4
   30 NSP=3+X(M)*0.5
      CALL SPLESQ(DX,X,W,M,A,NSP,2*IEQ,CHDF)
      X(1)=0.0
      YS(1)=YL(1)
      DO 32 I=2,M
      YS(I)=YL(I)
   32 X(I)=X(I-1)+DX(I)+DX(I-1)
      NSP=3+X(M)*0.5
      CALL SPLESQ(YS,X,W,M,A,NSP,2*IEQ,CHDF)
      FACGEN=1.0
      IF(CHDF.LT.0.75*CLIM) FACGEN=0.8
      IF(CHDF.GT.1.50*CLIM) FACGEN=1.2
      DO 34 I=1,M
      C(I)=(YS(I)-YL(I))**2
      IF(IDY.NE.0) C(I)=C(I)*W(I)
   34 CONTINUE
      CALL SMOPER(C,M,7)
      SUM=0.0
      DO 40 I=1,M
      FAC=1.0
      SUM=SUM+C(I)
      IF(C(I).GT.1.1*CLIM) GOTO 36
      IF(C(I).GT.0.9*CLIM) GOTO 38
      FAC=0.8
      IF(C(I).LT.0.5*CLIM) FAC=FAC*FAC
      GOTO 38
   36 FAC=1.2
      IF(C(I).GT.2.0*CLIM) FAC=FAC*FAC
   38 DX(I)=DX(I)*FAC*FACGEN
   39 DX(I)=AMAX1(DXA,DX(I))
      DX(I)=AMIN1(0.5,DX(I))
   40 CONTINUE
      SUM=SUM/FLOAT(M)
   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
