      	PROGRAM LUZZATI
C
C  Program to generate a Luzzati plot, a R-value vs. mean |Fo| plot,
C  and determine an estimate of the mean coordinate error (DELTA-R).
C  bcf   6-84
C
C  MODIFICATIONS BY R.S. ZHOU TO PREPARE A DATA FILE FOR GRAFIT (8/88).
C
     	PARAMETER MMD = 50000, MNP = 512
     	DIMENSION SOBS(MNP), ROBS(MNP), RCAL(MNP,4), DELR(4)
     	DIMENSION RDELS(25), RP(25), RPP(25), S(25), R(25)
     	CHARACTER FNAME*32, DNAME*32, PNAME*32, TEXT*64, KEY*2
     	DIMENSION PR(5), KP(5)
     	DIMENSION FBUF(504), IBUF(504)
     	INTEGER   HKL(3)
     	EQUIVALENCE (FBUF,IBUF)
     	COMMON/DATA/ NREFS, NBINS, CUTLOS, CUTHIS, CUTLOF, CUTHIF,
     +               FOS(MMD), FCS(MMD), SLS(MMD), SIG(MMD)
C
C  FOLLOWING DATA FROM LUZZATI, ACTA. CRYST. VOL.5,P807 (TABLEAU 2)
C
	DATA RDELS / 0.00, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07,
     +		     0.08, 0.09, 0.10, 0.12, 0.14, 0.16, 0.18, 0.20,
     +	       0.22, 0.24, 0.26, 0.28, 0.30, 0.35, 0.40, 0.45, 0.50 /
	DATA RP    /0.000,0.039,0.078,0.115,0.152,0.188,0.223,0.256,
     +		    0.288,0.320,0.350,0.410,0.462,0.510,0.554,0.595,
     +	      0.631,0.662,0.689,0.713,0.735,0.776,0.802,0.817,0.823 /
	DATA RPP   /0.000,0.025,0.050,0.074,0.098,0.122,0.145,0.168,
     +		    0.191,0.214,0.237,0.281,0.319,0.353,0.385,0.414,
     +	      0.440,0.463,0.483,0.502,0.518,0.548,0.564,0.574,0.580 /
C
C  INITIALIZE RUN PARAMETERS...
C
      	FNAME = ' '	!INPUT REFLECTION FILENAME      (*.BIN)
      	DNAME = ' '	!OUTPUT DATA FILENAME           (*.PLT)
      	PNAME = ' '	!OUTPUT PLOT PARAMETER FILENAME (*.PAR)
	NBINS = 1000	!NUMBER OF REFLECTIONS IN A SHELL (BIN SIZE)
      	ICENT = 0	!NON-CENTROSYMMETRIC SPACE GROUP
      	IGENT = 0	!DO NOT GENERATE A R-VALUE VS |Fo| CURVE
      	NDELR = 1	!NUMBER OF THEORETICAL DELR CURVES TO MAKE
      	BDELR = 0.0	!BEST FIT FOR CURRENT DATA SET
      	NREFS = 0	!NUMBER OF REFLECTIONS FROM INPUT FILE
      	MREFS = 0	!NUMBER OF REFLECTIONS AFTER CUTOFFS
     	IFORM = 1	!FILE FORMAT: DEFAULT TO PROLSQ FORMAT
     	CUTLOS = 0.0	!REFLECTION CUTOFF IN D-SPACING
     	CUTHIS = 0.0	!REFLECTION CUTOFF IN D-SPACING
     	CUTLOF = 0.0	!REFLECTION CUTOFF IN (FO-FC)/SIGF SPACE
     	CUTHIF = 0.0	!REFLECTION CUTOFF IN (FO-FC)/SIGF SPACE
C
C  LAYOUT THE DIALOGUE MENU...
C
     	CALL IFINIT(1)
     	CALL IFANSI_CLEAR(0)
      	CALL IFANSI_TITLE( 'LUZZATI: A program to prepare a
     + Luzzati plot for GRAFIT' )
      	CALL IFANSI_VIDEO(47)
      	CALL IFANSI_VIDEO(80)
      	CALL IFANSI_VIDEO(90)
      	WRITE(6,1000)
1000	FORMAT(/
     +,19X,'C U R R E N T   R U N   P A R A M E T E R S',18X,' '/)
     	CALL IFANSI_VIDEO(47)
     	WRITE(6,1001)
1001	FORMAT(
     +' FF=               ',
     +'    ==> Reflection file (contains      /      reflections).'/
     +' FM=               ',
     +'    ==> File format: 1=PROLSQ, 2=FMLS.'/
     +' SC=               ',
     +'    ==> Reflection cutoff window in d-spacing.'/
     +' FC=               ',
     +'    ==> Reflection cutoff window in (Fo-Fc)/Sigma(Fo) space.'/
     +' DF=               ',
     +'    ==> Output data file to be plotted by GRAFIT.'/
     +'                   ',
     +'    ==> Output plot parameter/setup file.'/
     +' BS=               ',
     +'    ==> Bin size (number of reflections in a shell).'/
     +' SG=               ',
     +'    ==> 0/1 for non-centro/centrosymmetric space group.'/
     +' RF=               ',
     +'    ==> 0/1 for not generating/generating a R-value vs |Fo|.'/
     +' NL=               ',
     +'    ==> Number of ideal Luzzati curves to generate (max=4).'/
     +' DR=0.00 0.00 0.00 0.00',
     +'==> Delta-R values for the ideal curves(best fit:    ).')

     	CALL IFANSI_VIDEO(90)
     	WRITE(6,1002)
1002	FORMAT(80X,' ')
     	CALL IFANSI_VIDEO(47)
      	WRITE(6,1003)
1003	FORMAT(5x'To edit the run parameters, enter AC=p1,p2...',
     +' where AC is the 2 letter'/5x'Access-Code designated for',
     +' the parameter and p1,p2... are the desired'/5x'parameter',
     +' values. To continue with the run, simply hit the RETURN key.')
     	CALL IFANSI_VIDEO(90)
     	WRITE(6,1002)
     	CALL IFANSI_VIDEO(47)
     	WRITE(6,1005)
1005	FORMAT(' oops:'/' INPUT>')
C
C  PRESENT THE DEFAULT RUN PARAMETERS...
C
2000	CALL IFANSI_VIDEO(30)				!REVERSE VIDEO ON
     	CALL IFANSI_TEXT( FNAME, 19, 5, 4 )
     	CALL IFANSI_IFL( NREFS, 1, 5, 5, 53 )
     	CALL IFANSI_IFL( MREFS, 1, 5, 5, 59 )
      	TEXT = ' '
     	CALL IFANSI_IFL( IFORM, 1, 2, 6, 4 )
     	CALL IFANSI_TEXT( TEXT,  17,  6, 6 )
     	CALL IFANSI_FFL( CUTLOS, 1, 9.4, 7, 4 )
     	CALL IFANSI_FFL( CUTHIS, 1, 9.4, 7, 14)
     	CALL IFANSI_FFL( CUTLOF, 1, 9.4, 8, 4 )
     	CALL IFANSI_FFL( CUTHIF, 1, 9.4, 8, 14)
     	CALL IFANSI_TEXT( DNAME, 19, 9, 4 )
     	CALL IFANSI_TEXT( PNAME, 19, 10, 4 )
     	CALL IFANSI_IFL( NBINS, 1, 4, 11, 4 )
     	CALL IFANSI_TEXT( TEXT,  15, 11, 8 )
     	CALL IFANSI_IFL( ICENT, 1, 2, 12, 4 )
     	CALL IFANSI_TEXT( TEXT,  17, 12, 6 )
     	CALL IFANSI_IFL( IGENT, 1, 2, 13, 4 )
     	CALL IFANSI_TEXT( TEXT,  17, 13, 6 )
     	CALL IFANSI_IFL( NDELR, 1, 2, 14, 4 )
     	CALL IFANSI_TEXT( TEXT,  17, 14, 6 )
      	DO I=1, 4
     	   CALL IFANSI_FFL( DELR(I), 1, 4.2, 15, 4+(I-1)*5 )
      	ENDDO
110	IF( BDELR.GT.0.0 ) CALL IFANSI_VIDEO(80)
      	CALL IFANSI_FFL( BDELR, 1, 4.2, 15, 72 )
C
C  START THE LOOP PROMPT...
C
111	CALL IFANSI_VIDEO(30)
     	CALL IFANSI_CURSOR(22,8)
     	CALL IFANSI_CLEAR(4)
      	READ(5,222,ERR=999,END=999) TEXT
222	FORMAT(A)
     	IF( TEXT.EQ.' ' ) GOTO 777
     	NCH = IFMISC_STRIP(TEXT)
     	CALL  IFMISC_UPPER(TEXT)
     	KEY = TEXT(1:2)
     	IF( KEY.EQ.'EX' .OR. KEY.EQ.'QU' ) GOTO 999
C
C  GET THE RUN PARAMETER NUMBER...
C
     	IF( NCH.LT.4 .OR. TEXT(3:3).NE.'=' ) GOTO 3333
     	TEXT = TEXT(4:)
C
C  READ REFLECTION FILE...
C
     	IF( KEY.EQ.'FF' ) THEN
     	    FNAME = TEXT
      	    CALL IFMISC_DEFEXT( FNAME, '.BIN' )
      	    OPEN(1,FILE=FNAME,STATUS='OLD',FORM='UNFORMATTED',ERR=1111)
     	    CALL IFOOPS(' Reading the reflection file...')

     	    IF( IFORM.EQ.1 ) THEN
     		DO  I=1, MMD
     		    READ(1,ERR=2222,END=555) HKL,A,B,FOS(I),SLS(I)
     		    FCS(I) = SQRT(A*A+B*B)
     		    SIG(I) = 1.0
     		ENDDO
  555		NREFS = I-1
     	    ELSE
     		NREFS = 0
     		READ(1,ERR=2222,END=444) I
  666		READ(1,ERR=2222,END=444) FBUF
     		DO  I = 1, 500, 10
     		IF( IBUF(I).EQ.99.0 ) GOTO 444
     		NREFS = NREFS + 1
     		FOS(NREFS) = FBUF(I+3)
     		FCS(NREFS) = FBUF(I+4)
     		SIG(NREFS) = FBUF(I+8)
     		SLS(NREFS) = FBUF(I+9)
     		ENDDO
     		GOTO 666
     	    ENDIF
  444	    IF( NREFS.EQ.0 ) GOTO 1112
     	    IF( NBINS.GT.NREFS/10 ) NBINS = NREFS/10
     	    IF( NBINS.EQ.0 ) NBINS = 1
     	    MREFS = NREFS
     	    DNAME = FNAME
     	    PNAME = FNAME
     	    CALL IFMISC_NAMEXT( DNAME, '.PLT' )	!GET DATA FILENAME
     	    CALL IFMISC_NAMEXT( PNAME, '.PAR' )	!GET SETUP FILENAME
     	    CLOSE( UNIT=1 )
     	    CALL IFOOPS_CLEAR
     	    BDELR = 0.0
     	    GOTO 2000
     	ENDIF
C
C  GET REFLECTION FILE FORMAT...
C
     	IF( KEY.EQ.'FM' ) THEN
     	    NCH = IFMISC_DECODE( TEXT, 1, PR, KP )
      	    IF( NCH.EQ.0 ) GOTO 3333
     	    IF( PR(1).LT.1.0 .OR. PR(1).GT.2.0 ) GOTO 3333
     	    IFORM = PR(1)
     	    CALL IFANSI_IFL( IFORM, 1, 2, 6, 4 )
      	    GOTO 111
     	ENDIF

      	IF( MREFS.EQ.0 ) GOTO 4444

     	IF( KEY.EQ.'SC' ) THEN		!GET CUTOFFS IN D-SPACING
     	    NCH = IFMISC_DECODE( TEXT, 2, PR, KP )
     	    IF( NCH.EQ.0 ) GOTO 3333
     	    IF( KP(1).NE.0.0 ) THEN
     		IF( PR(1).LT.0.0 .OR. PR(1).GT.999.0 ) GOTO 3333
     		CUTLOS = PR(1)
     	        CALL IFANSI_FFL( CUTLOS, 1, 9.4, 7, 4 )
     	    ENDIF
     	    IF( KP(2).NE.0.0 ) THEN
     		IF( PR(2).LT.0.0 .OR. PR(2).GT.999.0 ) GOTO 3333
     		CUTHIS = PR(2)
     	        CALL IFANSI_FFL( CUTHIS, 1, 9.4, 7, 14 )
     	    ENDIF
     	    MREFS = NCUTOFF( 0 )
     	    BDELR = 0.0
     	    GOTO 111
     	ENDIF

     	IF( KEY.EQ.'FC' ) THEN		!GET CUTOFFS IN F-SPACING
     	    NCH = IFMISC_DECODE( TEXT, 2, PR, KP )
     	    IF( NCH.EQ.0 ) GOTO 3333
     	    IF( KP(1).NE.0.0 ) THEN
     		IF( ABS(PR(1)).GT.999.0 ) GOTO 3333
     		CUTLOF = PR(1)
     	        CALL IFANSI_FFL( CUTLOF, 1, 9.4, 8, 4 )
     	    ENDIF
     	    IF( KP(2).NE.0.0 ) THEN
     		IF( ABS(PR(2)).GT.999.0 ) GOTO 3333
     		CUTHIF = PR(2)
     	        CALL IFANSI_FFL( CUTHIF, 1, 9.4, 8, 14 )
     	    ENDIF
     	    MREFS = NCUTOFF( 0 )
     	    BDELR = 0.0
     	    GOTO 111
     	ENDIF
C
C  GET THE OUTPUT DATA FILENAME...
C
     	IF( KEY.EQ.'DF' ) THEN
     	    DNAME = TEXT
     	    PNAME = DNAME
      	    CALL IFMISC_NAMEXT( PNAME, '.PAR' )
     	    CALL IFANSI_TEXT( DNAME, 19, 9, 4 )
     	    CALL IFANSI_TEXT( PNAME, 19,10, 4 )
      	    GOTO 111
     	ENDIF

     	IF( KEY.EQ.'BS' ) THEN		!GET BIN SIZE
      	    NCH = IFMISC_DECODE( TEXT, 1, PR, KP )
      	    IF( NCH.EQ.0 ) GOTO 3333
      	    IDUM = MREFS / INT(PR(1))
      	    IF( IDUM.LT.2 .OR. IDUM.GT.MNP ) GOTO 3333
      	    NBINS = PR(1)
     	    CALL IFANSI_IFL( NBINS, 1, 4, 11, 4 )
     	    BDELR = 0.0
      	    GOTO 111
      	ENDIF

      	IF( KEY.EQ.'SG' ) THEN		!GET CENTROID FLAG
      	    NCH = IFMISC_DECODE( TEXT, 1, PR, KP )
      	    IF( NCH.EQ.0 ) GOTO 3333
      	    ICENT = PR(1)
      	    IF( ICENT.NE.0 ) ICENT = 1
      	    CALL IFANSI_IFL( ICENT, 1, 2, 12, 4 )
     	    BDELR = 0.0
      	    GOTO 111
      	ENDIF

      	IF( KEY.EQ.'RF' ) THEN		!GET IGENT FLAG
      	    NCH = IFMISC_DECODE( TEXT, 1, PR, KP )
      	    IF( NCH.EQ.0 ) GOTO 3333
      	    IGENT = PR(1)
      	    IF( IGENT.NE.0 ) IGENT = 1
      	    CALL IFANSI_IFL( IGENT, 1, 2, 13, 4 )
     	    BDELR = 0.0
      	    GOTO 111
      	ENDIF

      	IF( KEY.EQ.'NL' ) THEN	!GET NUMBER OF LUZZATI CURVES
      	    NCH = IFMISC_DECODE( TEXT, 4, PR, KP )
      	    IF( NCH.EQ.0 ) GOTO 3333
      	    IF( PR(1).LT.0.0 .OR. PR(1).GT.4.0 ) GOTO 3333
      	    NDELR = PR(1)
     	    CALL IFANSI_IFL( NDELR, 1, 2, 14, 4 )
     	    GOTO 111
      	ENDIF

      	IF( KEY.EQ.'DR' ) THEN		!GET DELTA_R VALUES
      	    NCH = IFMISC_DECODE( TEXT, 4, PR, KP )
      	    IF( NCH.EQ.0 ) GOTO 3333
      	    DO I=1, NDELR
      	       IF( KP(I).EQ.1 ) THEN
      	           IF( PR(I).LT.0.01 .OR. PR(I).GT.2.0 ) GOTO 3333
      	           DELR(I) = PR(I)
     	           CALL IFANSI_FFL( DELR(I), 1, 4.2, 15, 4+(I-1)*5 )
      	       ENDIF
      	    ENDDO
     	    GOTO 111
      	ENDIF
     	GOTO 3333
1111	CALL IFOOPS(' Can not open/find the reflection file.')
      	CLOSE(UNIT=1)
      	GOTO 111
1112	CALL IFOOPS(' No reflection found in the input file.')
     	CLOSE(UNIT=1)
     	GOTO 111
2222	CALL IFOOPS(' Reflection file format not recognized.')
      	CLOSE(UNIT=1)
      	GOTO 111
3333	CALL IFOOPS(' Invalid entry or parameter out of range.')
      	GOTO 111
4444	CALL IFOOPS(' No data yet to work on.')
      	GOTO 111
C
C  CONTINUE WITH THE RUN...
C
777	IF( MREFS.EQ.0 ) GOTO 999
      	CALL IFOOPS(' Data reduction in process... please wait...')
      	IF( NDELR.EQ.1 .AND. BDELR.GT.0.0 ) GOTO 4000
      	IF( BDELR.GT.0.0 ) GOTO 2001
C
C  SORT THE REFLECTION ACCORDING TO SLS ARRAY FIRST...
C
     	MREFS = NCUTOFF( 1 )
	CALL IFMISC_SORT_TWO( MREFS, SLS, FOS, FCS )
      	NPTS = 0
      	DO I=1, MREFS, NBINS
      	   NPS  = 0
      	   SUMS = 0.0
      	   SUMD = 0.0
      	   SUMF = 0.0
      	   DO J= I, MIN(MREFS,I+NBINS-1)
      	      NPS  = NPS + 1
      	      SUMS = SUMS + SLS(J)
      	      SUMD = SUMD + ABS(FOS(J)-FCS(J))
      	      SUMF = SUMF + FOS(J)
      	   ENDDO
	   NPTS = NPTS + 1
      	   ROBS(NPTS) = SUMD/SUMF
      	   SOBS(NPTS) = SUMS/NPS
      	ENDDO
     	IF( NDELR.EQ.0 ) GOTO 4000
C
C  NOW GENERATE THEORETICAL CURVES FOR DIFFERENT DELTA-R VALUES.
C
2001	IDUM = NDELR
      	DO I=1, NDELR
      	   IF( DELR(I).EQ.0.0 ) THEN
      	       DO J=I, NDELR-1
      	          DELR(J) = DELR(J+1)
      	       ENDDO
      	       IDUM = IDUM - 1
      	   ENDIF
      	ENDDO
      	IF( NDELR.LT.1 ) NDELR = 1
      	IF( NDELR.EQ.1 .AND. DELR(1).EQ.0.0 ) GOTO 3001
      	DO J=1, NDELR
      	   DO I=1, 25
	      S(I) = RDELS(I) / DELR(J)
      	      IF( ICENT.EQ.0 ) R(I) = RPP(I)
      	      IF( ICENT.EQ.1 ) R(I) = RP(I)
      	   ENDDO
      	   DO I=1, NPTS
      	      DO K=1, 25
      	         IF( S(K).GT.SOBS(I) ) GO TO 3000
      	      ENDDO
3000	      RCAL(I,J)=FINTERP(SOBS(I),S(K-1),R(K-1),S(K),R(K))
      	   ENDDO
      	ENDDO
      	GOTO 4000
C
C  NOW FIND A VALUE OF THE ERROR THAT MOST CLOSELY RESEMBLES 
C  A SERIES OF POINTS REPRESENTING A LUZZATI CURVE.
C
3001	IF( BDELR.GT.0.0 ) GOTO 4000
      	SUMD = 1.0E+30
      	SUMF = 0.01
3100	SUMF = SUMF + 0.01
	SUMS = 0.0
	DO I=1, 25
      	   S(I) = RDELS(I) / SUMF
      	   IF( ICENT.EQ.0 ) R(I) = RPP(I)
      	   IF( ICENT.EQ.1 ) R(I) = RP(I)
      	ENDDO
	DO I=1, NPTS
      	   DO K=1, 25
      	      IF( S(K).GT.SOBS(I) ) GO TO 3200
      	   ENDDO
3200	   RCAL(I,2) = FINTERP(SOBS(I),S(K-1),R(K-1),S(K),R(K))
      	   SUMS = SUMS + (RCAL(I,2)-ROBS(I))**2
      	ENDDO
	IF( SUMS.LT.SUMD ) THEN
      	    BDELR= SUMF
      	    SUMD = SUMS
      	    DO I=1, NPTS
      	       RCAL(I,1) = RCAL(I,2)
      	    ENDDO
	ENDIF
	IF( SUMF.LT.2.0 ) GO TO 3100
      	CALL IFOOPS_CLEAR
      	GOTO 110
C
C  WRITE THE LUZZATI DATA FILE...
C
4000	OPEN(UNIT=2,FILE=DNAME,STATUS='NEW',FORM='FORMATTED',
     +  CARRIAGECONTROL='LIST' )
      	WRITE(2,4001) NPTS
4001	FORMAT(I5)
      	WRITE(2,4002) (SOBS(I),ROBS(I),I=1,NPTS)
4002	FORMAT(2G13.5)
      	DO J = 1, NDELR
      	   WRITE(2,4001) NPTS
      	   WRITE(2,4002) (SOBS(I),RCAL(I,J),I=1,NPTS)
      	ENDDO
      	IF( IGENT.EQ.0 ) GOTO 888
C
C  GENERATE THE R-VALUE VS MEAN |Fo| CURVE NEXT...
C
      	CALL IFMISC_SORT_TWO( MREFS, FOS, FCS, SLS )

      	NPTS = 0
      	DO I=1, MREFS, NBINS
      	   NPS  = 0
      	   SUMD = 0.0
      	   SUMF = 0.0
      	   DO J= I, MIN(MREFS,I+NBINS-1)
      	      NPS  = NPS + 1
      	      SUMD = SUMD + ABS(FOS(J)-FCS(J))
      	      SUMF = SUMF + FOS(J)
      	   ENDDO
	   NPTS = NPTS + 1
      	   ROBS(NPTS) = SUMD/SUMF
      	   SOBS(NPTS) = SUMF/NPS
      	ENDDO
C
C  WRITE THE R-VALUE VS MEAN |Fo| CURVE...
C
      	WRITE(2,4001) NPTS
      	WRITE(2,4002) (SOBS(I),ROBS(I),I=1,NPTS)
888	CLOSE( UNIT=2 )
C
C  WRITE THE PARAMETER OR SETUP FILE FOR GRAFIT...
C
      	OPEN( UNIT=2, FILE=PNAME, STATUS='NEW',
     +  FORM='FORMATTED', CARRIAGECONTROL='LIST' )

      	WRITE(2,5000) 
5000	FORMAT(
     +'TITL=Luzzati plot'/'XLB1=1/2d'/'YLB1=R-value'/
     +'XLB2=Mean |Fo|'/'YLB2=R-value'/
     +'SETS= 2, 1, 1, 1, 1, 1  0  0.08  0.0 0.0 1.0 1.0')
     	DO I=1, NDELR
     	WRITE(2,5001)
5001	FORMAT(
     +'SETS= 4, 0, 2, 1, 1, 1  0  0.08  0.0 0.0 1.0 1.0')
     	ENDDO
     	WRITE(2,5002)
5002	FORMAT(
     +'SETS= 3, 2, 1, 1, 2, 1  0  0.08  0.0 0.0 1.0 1.0')
      	CLOSE( UNIT=2 )

      	CALL IFOOPS_CLEAR
999	CALL IFANSI_VIDEO(47)
      	CALL EXIT
      	END

     	FUNCTION FINTERP (X, X1, Y1, X2, Y2)
     	FM = (Y2 - Y1) / (X2 - X1)
     	FINTERP = FM * (X - X1) + Y1
     	RETURN
     	END

     	SUBROUTINE IFOOPS( TEXT )
     	CHARACTER TEXT*(*)
     	CALL IFOOPS_CLEAR
     	NCH = IFMISC_LENG( TEXT )
     	CALL IFANSI_TEXT( TEXT, NCH, 21, 6 )
     	RETURN
     	END

     	SUBROUTINE IFOOPS_CLEAR
     	CALL IFANSI_CURSOR(21,6)
     	CALL IFANSI_CLEAR(4)
     	RETURN
     	END

     	FUNCTION NCUTOFF( ICODE )
     	PARAMETER MMD = 50000
     	COMMON/DATA/ NREFS, NBINS, CUTLOS, CUTHIS, CUTLOF, CUTHIF,
     +               FOS(MMD), FCS(MMD), SLS(MMD), SIG(MMD)

     	
     	NRS = NREFS
     	IF( CUTLOS.EQ.CUHIS .AND. CUTLOF.EQ.CUTHIF ) GOTO 999
     	ICD = ICODE
     	NRS = 0
     	DO  1111 I=1, NREFS
     	    IF( CUTLOS.NE.CUTHIS ) THEN
     		DSP = 1.0/(2.0*SLS(I))
     		IF( DSP.LT.CUTLOS .OR. DSP.GT.CUTHIS ) GOTO 1111
     	    ENDIF
     	    IF( CUTLOF.NE.CUTHIF ) THEN
     		DSP = (FOS(I)-FCS(I))/SIG(I)
     		IF( DSP.LT.CUTLOF .OR. DSP.GT.CUTHIF ) GOTO 1111
     	    ENDIF
     	    NRS = NRS + 1
     	    IF( ICD.EQ.1 ) THEN
     		FOS(NRS) = FOS(I)
     		FCS(NRS) = FCS(I)
     		SIG(NRS) = SIG(I)
     		SLS(NRS) = SLS(I)
     	    ENDIF
1111	CONTINUE
     	IF( NBINS.GT.NRS/10 ) NBINS = NRS/10
     	IF( NBINS.EQ.0 ) NBINS = 1

999	NCUTOFF = NRS
     	CALL IFANSI_IFL(   NRS, 1, 5, 5, 59 )
     	CALL IFANSI_IFL( NBINS, 1, 4, 11, 4 )

     	RETURN
     	END

