C***********************************************************************
C*                                                                     *
C*                             STATRIFL                                *
C*                                                                     *
C*                           by M.NARDELLI                             *
C* Dipartimento di Chimica Generale ed Inorganica, Chimica Analitica,  *
C* Chimica Fisica della Universita` degli Studi di Parma, Centro di    *
C* Studio per la Strutturistica Diffrattometrica del C.N.R., Viale     *
C* delle Scienze, I-43100 Parama, Italy.                               *
C*                                                                     *
C*    This routine considers the distribution of the unobserved        *
C*    reflections in reciprocal space, i.e. in shells of sin(theta)/   *
C*    lambda                                                           *
C*                                                                     *
C*    INPUT: the *.fcf  or *.hkl files for SHELXL-93                   *
C*                                                                     *
C*    SUBROUTINES:                                                     *
C*                 INPOUT                                              *
C*                 RCS                                                 *
C*                 SINTL                                               *
C***********************************************************************
      PROGRAM STATRIFL
      CHARACTER*1 TITL(80)
      CHARACTER*14 FILIN,FILOUT
      DIMENSION C(3),A(3),CR(3),CSR(3),N(50),IH(3),NP(50),SPO(50),
     1 SPU(50)
      COMMON/INPOU/IN,IO
C-----INPUT & OUTPUT files
      IN=13
      IO=15
      CALL INPOUT
      WRITE(*,102)
102   FORMAT(' Key in  0  if data are from the  .fcf  file'/
     1       ' key in  1  if data are from the  .hkl  file'/
     2 ' ?>')
      READ(*,*)KDC
C-----Read data
      WRITE(*,'('' Key in the name of the compound''/'' ?>'')')
      READ(*,'(80A1)') TITL
      WRITE(IO,'(1X,80A1/)') TITL
      WRITE(*,'('' Key in the lattice constants and wavelength ''
     1 ''(Free Format)''/'' ?>'')')
      READ(*,*) C,A,AL
      CALL RCS(C,A,CR,CSR)
      STLMX=1.0/AL
      IF(KDC.EQ.0)THEN
      DO 7 I=1,17
7     READ(IN,'(A80)')TIT
      ENDIF
      AA=0.05
      NT=0
      NO=0
      NN=0
      TSP=0.0
      OSP=0.0
      USP=0.0
      DO 1 I=1,50
      N(I)=0
      NP(I)=0
      SPO(I)=0.0
1     SPU(I)=0.0
2     CONTINUE
      IF(KDC.EQ.0)READ(IN,103,END=5,ERR=5)IH,F2,SF2
      IF(KDC.EQ.1)READ(IN,104,END=5,ERR=5)IH,F2,SF2
103   FORMAT(3I4,12X,F12.2,F10.2)
104   FORMAT(3I4,2F8.2)
      IF(IH(1).EQ.0.AND.IH(2).EQ.0.AND.IH(3).EQ.0)GO TO 5
      NT=NT+1
      TSP=TSP+F2
      IF(F2.GE.(2.0*SF2))THEN
      NO=NO+1
      OSP=OSP+F2
      CALL SINTL(IH,CR,CSR,STL)
      DO 3 I=1,50
      IF(STL.GT.(STLMX+0.05)) GO TO 5
      IF(STL.GT.(AA*(I-1)).AND.STL.LT.(AA*I)) THEN
      NP(I)=NP(I)+1
      SPO(I)=SPO(I)+F2
      ENDIF
3     CONTINUE
      GO TO 2
      ELSE
      NN=NN+1
      USP=USP+F2
      ENDIF
      CALL SINTL(IH,CR,CSR,STL)
      DO 4 I=1,50
      IF(STL.GT.(STLMX+0.05)) GO TO 5
      IF(STL.GT.(AA*(I-1)).AND.STL.LT.(AA*I)) THEN
      N(I)=N(I)+1
      SPU(I)=SPU(I)+F2
      ENDIF
4     CONTINUE
      GO TO 2
5     WRITE(IO,105)NT,TSP,NO,OSP,NN,USP
105   FORMAT(' No. Total reflections      =',I7,3X,
     1 'Total scattering power        =', F10.1/
     2       ' No. Observed reflections   =',I7,3X,
     3 '"Observed" scattering power   =', F10.1/
     4       ' No. Unobserved reflections =',I7,3X,
     5 '"Unobserved" scattering power =', F10.1)
      WRITE(IO,106)
106   FORMAT(/' Numbers of reflections and scattering power in the ',
     1 'sin(theta)/lambda ranges:'/' sin(theta)/lambda',2X,'No.',3X,
     2 '%',5X,'No.',4X,'%',7X,'scat. pow.',2X,'%',5X,'scat. pow.',
     3 1X,'%'/20X,'observed',4X,'unobserved',7X,'observed',9X,
     4 'unobserved')
      NK=STLMX/0.05+1
      DO 6 I=1,NK
      RN=N(I)
      RNN=NN
      PC=RN/RNN*100
      PN=NP(I)
      PNO=NO
      PP=PN/PNO*100
      PSPO=SPO(I)/OSP*100
      PSPU=SPU(I)/USP*100
      A1=AA*(I-1)
      A2=AA*I
6     WRITE(IO,107) A1,A2,NP(I),PP,N(I),PC,SPO(I),PSPO,SPU(I),PSPU
107   FORMAT(2X,F4.2,' - ',F4.2,I9,F6.1,I7,F6.1,F13.1,F6.1,F10.1,F6.1)
      STOP
      END
 
      SUBROUTINE INPOUT
C-----Assignes the input and output files
      CHARACTER*10 STR1
      CHARACTER*8,STR2
      CHARACTER*14 FILIN
      COMMON/INPOU/IN,IO
      STR1='Key in the'
      STR2='put file'
      WRITE(*,'(1X,A10,'' in'',A8/'' ?>'')')STR1,STR2
      READ(*,'(A14)')FILIN
      OPEN(UNIT=IN,FILE=FILIN,FORM='FORMATTED')
      REWIND(UNIT=IN)
      WRITE(*,'(1X,A10,'' out'',A8/'' ?>'')')STR1,STR2
      READ(*,'(A14)')FILIN
      OPEN(UNIT=IO,FILE=FILIN,FORM='FORMATTED')
      REWIND(UNIT=IO)
      RETURN
      END
 
      SUBROUTINE RCS(C,A,CR,CSR)
C-----Calculates reciprocal cell constants
      DIMENSION C(3),A(3),SN(3),CS(3),CR(3),CSR(3)
      RD=0.01745329
      S=0
      DO 1 I=1,3
      A(I)=A(I)*RD
      SN(I)=SIN(A(I))
      CS(I)=COS(A(I))
1     S=S+A(I)
      S=S/2.0
      V=2.0*C(1)*C(2)*C(3)*SQRT(SIN(S)*SIN(S-A(1))*SIN(S-A(2))*
     1 SIN(S-A(3)))
      CR(1)=C(2)*C(3)*SN(1)/V
      CR(2)=C(3)*C(1)*SN(2)/V
      CR(3)=C(1)*C(2)*SN(3)/V
      CSR(1)=(CS(2)*CS(3)-CS(1))/(SN(2)*SN(3))
      CSR(2)=(CS(3)*CS(1)-CS(2))/(SN(3)*SN(1))
      CSR(3)=(CS(1)*CS(2)-CS(3))/(SN(1)*SN(2))
      RETURN
      END
 
      SUBROUTINE SINTL(IH,CR,CSR,STL)
C-----Calculates sin(theta)/lambda
      DIMENSION IH(3),CR(3),CSR(3)
      IO=15
      STL=0
      DO 1 I=1,3
1     STL=STL+(IH(I)*CR(I))**2
      STL=0.5*SQRT(STL+2.0*(IH(2)*IH(3)*CR(2)*CR(3)*CSR(1)+
     1 IH(3)*IH(1)*CR(3)*CR(1)*CSR(2)+
     1 IH(1)*IH(2)*CR(1)*CR(2)*CSR(3)))
      RETURN
      END
