      	PROGRAM LUZATI
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 (MNP= 128)
     	DIMENSION NSHL(MNP), FOBS(MNP), SOBS(MNP), ROBS(MNP)
        DIMENSION RCAL(MNP,4), DELR(4), PR(5), KP(5)
     	CHARACTER FNAME*32, DNAME*32, PNAME*32, TEXT*80, KEY*2
     	DIMENSION RDELS(25), RP(25), RPP(25), S(25), R(25)
     	INTEGER   HKL(3)
        LOGICAL   YES
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)
	SHELL = 0.01	!SHELL WIDTH IN 1/2D
      	ICENT = 0	!NON-CENTROSYMMETRIC SPACE GROUP
      	NREFS = 0	!NUMBER OF REFLECTIONS FROM INPUT FILE
      	MREFS = 0	!NUMBER OF REFLECTIONS AFTER CUTOFFS
     	CUTL  = 0.0	!REFLECTION CUTOFF IN 1/2D-SPACING
     	CUTH  = 0.0	!REFLECTION CUTOFF IN 1/2D-SPACING
        BDELR = 0.0	!BEST FIT FOR CURRENT DATA SET
C
C  LAYOUT THE DIALOGUE MENU...
C
     	CALL IFINIT_SESSION
     	CALL IFANSI_CLEAR(0)
        TEXT(1:43) = 'LUZATI: A program to prepare a Luzzati plot'
        TEXT(44: ) = ' for GRAFIT''S XY-Plot'
      	CALL IFANSI_TITLE( TEXT )
      	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).'/
     +' DF=               ',
     +'    ==> Output data file to be plotted by GRAFIT.'/
     +'                   ',
     +'    ==> Output plot parameter/setup file.'/
     +' PR=               ',
     +'    ==> Plot resolution (shell width in 1/2d).'/
     +' DS=               ',
     +'    ==> Reflection cutoff window in 1/2d space.'/
     +' SG=               ',
     +'    ==> 0/1 for non-centro/centrosymmetric space group.'/
     +' 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. Plot resolution and data cutoff should be set',
     +' before'/,5X,'entering the data filename.     Enter GO',
     +' to start computing.')
     	CALL IFANSI_VIDEO(90)
     	WRITE(6,1002)
     	CALL IFANSI_VIDEO(47)
     	WRITE(6,1005)
1005	FORMAT(' oops:')
C
C  PRESENT THE DEFAULT RUN PARAMETERS...
C
2000	CALL IFANSI_VIDEO(30)				!REVERSE VIDEO ON
      	TEXT = ' '
     	CALL IFANSI_TEXT( FNAME, 19, 5, 4 )
     	CALL IFANSI_IFL( NREFS, 1, 'I5', 5, 53 )
     	CALL IFANSI_IFL( MREFS, 1, 'I5', 5, 59 )
     	CALL IFANSI_TEXT( DNAME, 19, 6, 4 )
     	CALL IFANSI_TEXT( PNAME, 19, 7, 4 )
     	CALL IFANSI_FFL( SHELL, 1, 'F5.3', 8, 4 )
     	CALL IFANSI_TEXT( TEXT,  14, 8, 9 )
     	CALL IFANSI_FFL( CUTL, 1, 'F5.3', 9, 4 )
     	CALL IFANSI_FFL( CUTH, 1, 'F5.3', 9, 10)
     	CALL IFANSI_TEXT( TEXT,  8, 9, 15 )
     	CALL IFANSI_IFL( ICENT, 1, 'I2', 10, 4 )
     	CALL IFANSI_TEXT( TEXT,  17, 10, 6 )
      	DO I=1, 4
     	   CALL IFANSI_FFL( DELR(I), 1, 'F4.2', 11, 4+(I-1)*5 )
      	ENDDO
110	IF( BDELR.GT.0.0 ) CALL IFANSI_VIDEO(80)
      	CALL IFANSI_FFL( BDELR, 1, 'F4.2', 11, 72 )
C
C  START THE DIALOGUE PROMPT...
C
111	CALL IFANSI_VIDEO(47)
        CALL IFANSI_TEXT('INPUT> ', 7, 22, 1 )
     	CALL IFANSI_CLEAR(4)
        CALL IFANSI_VIDEO(30)
      	READ(5,222,ERR=999,END=999) TEXT
222	FORMAT(A)
     	IF( TEXT.EQ.' ' ) GOTO 111
     	NCH = IFMISC_STRIP(TEXT)
     	CALL  IFMISC_UPPER(NCH,TEXT)
     	KEY = TEXT(1:2)
     	IF( KEY.EQ.'GO' ) GOTO 777
     	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:)
     	NCH  = NCH - 3
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 data file (CTRL/C to abort)...')
            NREFS = 0
     	    MREFS = 0
            YES = .FALSE.
            IF( CUTH.GT.CUTL ) YES = .TRUE.
            DO IS=1,128
               NSHL(IS) = 0
               FOBS(IS) = 0
               SOBS(IS) = 0
               ROBS(IS) = 0
            ENDDO
            CALL IFQAST_CTRLC     !TRAP A CTRL/C INTERRUPT
 1234       READ(1,ERR=2222,END=555) HKL, A, B, FO, SL
            IF( IFCTRL_C() ) THEN
                NREFS = 0
                MREFS = 0
                CALL IFOOPS_CLEAR
                GOTO 111
            ENDIF
            MREFS = MREFS + 1
            IF( YES .AND. (SL.LT.CUTL.OR.SL.GT.CUTH) ) GOTO 1234
            IS = INT(SL/SHELL) + 1
            IF( IS.GT.128 ) GOTO 5555    !RESOLUTION TOO HIGH
            FC = SQRT(A*A+B*B)
            NSHL(IS) = NSHL(IS) + 1
            ROBS(IS) = ROBS(IS) + ABS(FO-FC)
            SOBS(IS) = SOBS(IS) + SL
            FOBS(IS) = FOBS(IS) + FO
            NREFS = NREFS + 1
            GOTO 1234
 555	    CALL IFQAST_CLEAR            !CLEAR THE CTRL/C TRAP
            IF( NREFS.EQ.0 ) GOTO 1112
            DO IS=1, 128
               NCH = NSHL(IS)
               IF( NCH.GT.0 ) THEN
                   ROBS(IS) = ROBS(IS)/FOBS(IS)
                   SOBS(IS) = SOBS(IS)/NCH
               ENDIF
            ENDDO
     	    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 THE OUTPUT DATA FILENAME...
C
     	IF( KEY.EQ.'DF' ) THEN
     	    DNAME = TEXT
     	    PNAME = DNAME
      	    CALL IFMISC_NAMEXT( PNAME, '.PAR' )
     	    CALL IFANSI_TEXT( DNAME, 19, 6, 4 )
     	    CALL IFANSI_TEXT( PNAME, 19, 7, 4 )
      	    GOTO 111
     	ENDIF
C
C  GET PLOT RESOLUTION...
C
     	IF( KEY.EQ.'PR' ) THEN		!GET BIN SIZE
      	    NCH = IFMISC_DECODE( TEXT, NCH, 1, PR, KP )
      	    IF( NCH.EQ.0 ) GOTO 3333
      	    IF( PR(1).LT.0.001 .OR. PR(1).GT.0.1 ) GOTO 3333
      	    SHELL = PR(1)
     	    CALL IFANSI_FFL( SHELL, 1, 'F5.3', 8, 4 )
      	    GOTO 111
      	ENDIF
C
C  GET D-SPACING CUTOFFS...
C
     	IF( KEY.EQ.'DS' ) THEN
     	    NCH = IFMISC_DECODE( TEXT, NCH, 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.9.0 ) GOTO 3333
     		CUTL = PR(1)
     	    ENDIF
     	    IF( KP(2).NE.0.0 ) THEN
     		IF( PR(2).LT.0.0 .OR. PR(2).GT.9.0 ) GOTO 3333
     		CUTH = PR(2)
     	    ENDIF
            IF( CUTL.GT.CUTH ) THEN
                 PR(1)= CUTL
                 CUTL = CUTH
                 CUTH = PR(1)
            ENDIF
            CALL IFANSI_FFL( CUTL, 1, 'F5.3', 9, 4 )
            CALL IFANSI_FFL( CUTH, 1, 'F5.3', 9, 10 )
     	    GOTO 111
     	ENDIF
C
C  GET CENTROID FLAG...
C
      	IF( KEY.EQ.'SG' ) THEN
      	    NCH = IFMISC_DECODE( TEXT, NCH, 1, PR, KP )
      	    IF( NCH.EQ.0 ) GOTO 3333
      	    ICENT = PR(1)
      	    IF( ICENT.NE.0 ) ICENT = 1
      	    CALL IFANSI_IFL( ICENT, 1, 'I2', 10, 4 )
     	    BDELR = 0.0
      	    GOTO 111
      	ENDIF
C
C  GET DELTA_R VALUES...
C
      	IF( KEY.EQ.'DR' ) THEN
      	    NCH = IFMISC_DECODE( TEXT, NCH, 4, PR, KP )
      	    IF( NCH.EQ.0 ) GOTO 3333
      	    DO I=1, 4
      	       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,'F4.2',11,4+(I-1)*5 )
      	       ENDIF
      	    ENDDO
     	    GOTO 111
      	ENDIF
     	GOTO 3333
1111	CALL IFOOPS(' Can not open/find the data file.')
      	CLOSE(UNIT=1)
      	GOTO 111
1112	CALL IFOOPS(' No data found in the input file.')
     	CLOSE(UNIT=1)
     	GOTO 111
2222	CALL IFOOPS(' Data 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
5555    CALL IFOOPS(' Plot resolution too high.')
        MREFS = 0
        NREFS = 0
      	GOTO 111
777	IF( MREFS.EQ.0 ) GOTO 4444
C
C  COMPUTE LUZZATI CURVES...
C
      	CALL IFOOPS(' Working... Please wait...')
        IF( BDELR.GT.0.0 ) GOTO 4321
C
C  FIND A VALUE OF THE ERROR THAT MOST CLOSELY RESEMBLES 
C  A SERIES OF POINTS REPRESENTING A LUZZATI CURVE.
C
      	SUMD = 1.0E+30
      	SUMF = 0.01
3100	SUMF = SUMF + 0.01
	SUMS = 0.0
	DO IS=1, 25
      	   S(IS) = RDELS(IS) / SUMF
      	   IF( ICENT.EQ.0 ) THEN
               R(IS) = RPP(IS)
           ELSE
      	       R(IS) = RP(IS)
           ENDIF
      	ENDDO
	DO IS=1, 128
           IF( NSHL(IS).GT.0 ) THEN
      	       DO K=1, 25
      	          IF( S(K).GT.SOBS(IS) ) GO TO 3200
      	       ENDDO
3200	       RCAL(IS,2)=FINTERP(SOBS(IS),S(K-1),R(K-1),S(K),R(K))
      	       SUMS = SUMS + (RCAL(IS,2)-ROBS(IS))**2
           ENDIF
      	ENDDO
	IF( SUMS.LT.SUMD ) THEN
      	    BDELR= SUMF
      	    SUMD = SUMS
      	    DO IS=1, 128
      	       RCAL(IS,1) = RCAL(IS,2)
      	    ENDDO
	ENDIF
	IF( SUMF.LT.2.0 ) GO TO 3100
      	CALL IFOOPS_CLEAR
      	GOTO 110
C
C  NOW GENERATE THEORETICAL CURVES FOR DIFFERENT DELTA-R VALUES.
C
4321	DO J=1, 4
           IF( DELR(J).GT.0.0 ) GOTO 4322
        ENDDO
        DELR(1) = BDELR
        GOTO 4000
4322    DO J=1, 4
           IF( DELR(J).GT.0.0 ) THEN
      	       DO IS=1, 25
	          S(IS) = RDELS(IS) / DELR(J)
      	          IF( ICENT.EQ.0 ) THEN
                      R(IS) = RPP(IS)
      	          ELSE
                      R(IS) = RP(IS)
                  ENDIF
      	       ENDDO
      	       DO IS=1, 128
                  IF( NSHL(IS).GT.0 ) THEN
      	              DO K=1, 25
      	                 IF( S(K).GT.SOBS(IS) ) GO TO 3000
      	              ENDDO
3000	           RCAL(IS,J)=FINTERP(SOBS(IS),S(K-1),R(K-1),S(K),R(K))
                  ENDIF
      	       ENDDO
           ENDIF
      	ENDDO
C
C  WRITE THE PLOT DATA FILE...
C
4000	OPEN(UNIT=2,FILE=DNAME,STATUS='NEW',FORM='FORMATTED',
     +  CARRIAGECONTROL='LIST' )
        J = 1
        DO IS=1, 4
           IF( DELR(IS).GT.0 ) J = J + 1
        ENDDO
      	WRITE(2,4001) -J
4001    FORMAT(I3)
        DO IS = 1, 128
           IF( NSHL(IS).GT.0 ) WRITE(2,4002)
     +     SOBS(IS), ROBS(IS), (RCAL(IS,I),I=1,J-1)
        ENDDO
4002	FORMAT(<J+1>G13.5)
        CLOSE(UNIT=2)
C
C  WRITE THE 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'/
     +'SETS= 2 1 1 1 1 1 0 0.08 0.0 0.0 1.0 1.0')
     	DO I=1, J-1
     	   WRITE(2,5001)
5001	   FORMAT('SETS= 4 0 2 1 1 1 0 0.08 0.0 0.0 1.0 1.0')
     	ENDDO
      	CLOSE(UNIT=2)
      	CALL IFOOPS_CLEAR
        GOTO 111
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

