      SUBROUTINE READFX
      CHARACTER*8 AN
      COMMON NC,NV,IRT,IPR,ISYM,NF,NT,NA,NS,NP,IC,IWARN,I10,
     1 WVL(2),SC(2),BO(2),SINT,DIDBO(2),SINANG,NDT,KDT,
     2 PRCP(6),DR(9,2),DI(9,2),BI(9),DIDB(9),SF(9,32,2),
     3 TS(3,48),IS(2,3,48),KI(253),P(253),DA(253),
     4 AN(60),ISF(60),ITF(60),AI(60),DIDA(60),
     5 X(3,60),DIDX(3,60)
      DIMENSION  DCP(6),RCP(6)
      DO 10  K=1,NDT
      DO 9   I=1,NF
      READ(60,1010)(SF(I,J,K),J=1,32)
 1010 FORMAT(8F9.3)
      IF(IRT.EQ.0)  GO TO 9
      WRITE(61,1020)I,K,(SF(I,J,K),J=1,32)
 1020 FORMAT('0SCATTERING TABLE'I3,',  SET'I3,/ (1X,8F12.3))
    9 CONTINUE
   10 CONTINUE
      IF(ISYM.LT.3)  GO TO 30
      IF(IRT.GT.0)  WRITE(61,1030)
 1030 FORMAT(/'0ANOMALOUS DISPERSIONS')
      DO 20  K=1,NDT
      DO 19  I=1,NF
      READ(60,1040)DR(I,K),DI(I,K)
 1040 FORMAT(2F9.3)
      IF(IRT.EQ.0)  GO TO 19
      WRITE(61,1050)I,K,DR(I,K),DI(I,K)
 1050 FORMAT('0TYPE'I3,'  SET'I3,4X,'REAL='F12.3,4X,'IMAGINARY='F12.3)
   19 CONTINUE
   20 CONTINUE
   30 CONTINUE
      READ(60,1120)  ((TS(J,I),(IS(K,J,I),K=1,2),J=1,3),I=1,NS)
 1120 FORMAT(3(F11.6,2I2))
C***  READ AND TRANSFORM CELL PARAMETERS
      READ(60,1060)(DCP(J),J=1,6)
 1060 FORMAT(6F9.6)
      V=SQRT(1.0+2.0*DCP(4)*DCP(5)*DCP(6)-DCP(4)**2-DCP(5)**2-
     1  DCP(6)**2)
      DO 40  J=1,3
      RCP(J)=SQRT(1.0-DCP(J+3)**2)/(DCP(J)*V)
   40 CONTINUE
      RCP(4)=(DCP(5)*DCP(6)-DCP(4))/SQRT((1.0-DCP(5)**2)*(1.0-DCP(6)**2
     1 ))
      RCP(5)=(DCP(6)*DCP(4)-DCP(5))/SQRT((1.0-DCP(6)**2)*(1.0-DCP(4)**2
     1 ))
      RCP(6)=(DCP(4)*DCP(5)-DCP(6))/SQRT((1.0-DCP(4)**2)*(1.0-DCP(5)**2
     1 ))
      IF(DCP(1).GT.1.0)  GO TO 70
      DO 60  J=1,6
      CP=DCP(J)  
      DCP(J)=RCP(J)
      RCP(J)=CP
   60 CONTINUE
   70 CONTINUE
      PRCP(1)=0.25*RCP(1)**2
      PRCP(2)=0.25*RCP(2)**2
      PRCP(3)=0.25*RCP(3)**2
      PRCP(4)=0.5*RCP(2)*RCP(3)*RCP(4)
      PRCP(5)=0.5*RCP(1)*RCP(3)*RCP(5)
      PRCP(6)=0.5*RCP(1)*RCP(2)*RCP(6)
      READ(60,1090)WVL(1),WVL(2)
 1090 FORMAT(2F9.6)
      IF(WVL(1).EQ.0.0)  WVL(1)=1.54178
      IF(IRT.EQ.0)  GO TO 71
      WRITE(61,1130) ((TS(J,I),(IS(K,J,I),K=1,2),J=1,3),I=1,NS)
 1130 FORMAT(/'0SYMMETRY TRANSFORMATIONS'/(F13.6,2I2,2(F20.6,2I2)))
      WRITE(61,1070)(DCP(J),J=1,6)
 1070 FORMAT(/'0REAL CELL PARAMETERS'3F12.4,3F12.6)
      WRITE(61,1080)(RCP(J),J=1,6)
 1080 FORMAT(/'0RECIPROCAL PARAMETERS'3F12.7,3F12.6)
      WRITE(61,1100)WVL(1)
 1100 FORMAT(/'0WAVELENGTH FOR SET 1 = 'F9.6)
      IF(NDT.GT.1)WRITE(61,1101)WVL(2)
 1101 FORMAT('                    2 = 'F9.6)
   71 NP=4*NA+NT+4
      DO 72  I=1,NP
   72 KI(I)=0
      IF(NC.EQ.1)  GO TO 80
      READ(60,1110)(KI(I),I=1,NP)
 1110 FORMAT(80I1)
      ICV=0  
      IWARN=0
      DO 75  I=1,NP
      IF(KI(I).EQ.1)  ICV=ICV+1
   75 CONTINUE
      IF(ICV.EQ.NV)  GO TO 80
      IWARN=1
      WRITE(61,1140)
 1140 FORMAT(/'0INCONSISTENCY IN PARAMETERS TO BE VARIED')
   80 CONTINUE
      RETURN
      END
      SUBROUTINE READPR(LUC)
C***  READS INPUT PARAMETERS FROM CARDS OR LU1
      CHARACTER*8 AN
      COMMON NC,NV,IRT,IPR,ISYM,NF,NT,NA,NS,NP,IC,IWARN,I10,
     1 WVL(2),SC(2),BO(2),SINT,DIDBO(2),SINANG,NDT,KDT,
     2 PRCP(6),DR(9,2),DI(9,2),BI(9),DIDB(9),SF(9,32,2),
     3 TS(3,48),IS(2,3,48),KI(253),P(253),DA(253),
     4 AN(60),ISF(60),ITF(60),AI(60),DIDA(60),
     5 X(3,60),DIDX(3,60)
      WRITE(61,1000)
 1000 FORMAT(/21H0**** INPUT DATA ****)
      READ(LUC,1015)SC(1),SC(2)
      READ(LUC,1010)BO(1),BO(2)
 1015 FORMAT(E12.6)
      READ(LUC,1010)(BI(I),I=1,NT)
 1010 FORMAT(F9.6)
      WRITE(61,1020)SC(1),KI(1),SC(2),KI(2),BO(1),KI(3),BO(2),KI(4),
     +(I,BI(I),KI(I+4),I=1,NT)
 1020 FORMAT('0SCALE FACTOR 1='F12.6,I4,/'              2='F12.6,I4,/
     +'0OVERALL TEMPERATURE FACTOR 1='F12.6,I4,/' ',27X,'2='F12.6,I4,/
     +'0INDIVIDUAL TEMPERATURE FACTORS'(/1X,I3,F14.6,I4))
      WRITE(61,1030)
 1030 FORMAT(/'0  ATOM     SF   TF  MULTIPLICITY'9X,'X'14X,'Y'14X,'Z')
      DO 10  I=1,NA
      READ(LUC,1040)AN(I),ISF(I),ITF(I),AI(I),(X(J,I),J=1,3)
 1040 FORMAT(A8,2X,I1,2X,I1,4X,4F9.6)
   10 WRITE(61,1050)AN(I),ISF(I),ITF(I),AI(I),KI(4*I+NT+1),
     1 (X(J,I),KI(4*I+NT+J+1),J=1,3)
 1050 FORMAT(1H ,A8,2I5,4(F12.6,I3))
      WRITE(61,1060)NP
 1060 FORMAT(/1H0,I3,'  PARAMETERS READ')
      CALL COPTOP(P,SC,BO,BI,AI,X,NA,NT)
      RETURN
      END
      SUBROUTINE CALCIS(H,M,PL,RHO,IO)
      CHARACTER*8 AN
      INTEGER H
      COMMON NC,NV,IRT,IPR,ISYM,NF,NT,NA,NS,NP,IC,IWARN,I10,
     1 WVL(2),SC(2),BO(2),SINT,DIDBO(2),SINANG,NDT,KDT,
     2 PRCP(6),DR(9,2),DI(9,2),BI(9),DIDB(9),SF(9,32,2),
     3 TS(3,48),IS(2,3,48),KI(253),P(253),DA(253),
     4 AN(60),ISF(60),ITF(60),AI(60),DIDA(60),
     5 X(3,60),DIDX(3,60)
      DIMENSION DADA(60),DBDA(60),DADX(3,60),DBDX(3,60),
     1 DADAO(60),DBDAO(60),DADXO(3,60),DBDXO(3,60)
      DIMENSION H(3),HT(3),HTO(3),DADB(9),DBDB(9),SR(9),TF(9),SI(9)
      ISKIP=(IC-NC)*IO
      A=0.0
      B=0.0
      RSA=SQRT(RHO)/0.05
      ISA=RSA
      FSA=RSA-ISA
C***  CLEAR ARRAYS NEEDED FOR SUMS
      DO 5  I=1,NT
      DADB(I)=0.0
    5 DBDB(I)=0.0
      DO 20  I=1,NA
      DADAO(I)=0.0
      DBDAO(I)=0.0
      DBDA(I)=0.0
      DO 10  J=1,3
      DADXO(J,I)=0.0
      DBDXO(J,I)=0.0
   10 DBDX(J,I)=0.0
   20 CONTINUE
C***  START LOOP THROUGH NS SYMMETRY CARDS
      DO 70  K=1,NS
C***  CALCULATE TRANSFORMED INDICES
      TT=0.0
      DO 32  J=1,3
      TT=TT+H(J)*TS(J,K)
      HT(J)=0.0
      DO 31  JJ=1,3
      DO 30  L=1,2
      IM=IS(L,JJ,K)
      IF(IM.EQ.J)  HT(J)=HT(J)+H(JJ)
      IF(IM.EQ.-J) HT(J)=HT(J)-H(JJ)
   30 CONTINUE
   31 CONTINUE
   32 CONTINUE
C***  END INDEX TRANSFORMATION
C***     RETAIN TRANSFORMED INDICES AGAINST ANY PATCHING
      TTO=TT
      DO 35 J=1,3
   35 HTO(J)=HT(J)
C***  CALCULATE TRIGONOMETRIC TERMS
C***  START LOOP THROUGHT NA ATOMS
      DO 65  I=1,NA
C***  CALL PATCH
      CALL PATCH(TT,HT,I)
      HX=TT
      DO 40  J=1,3
   40 HX=HX+HT(J)*X(J,I)
      COSIK=COS(6.28318531*HX)
      SINIK=SIN(6.28318531*HX)
      DADAO(I)=DADAO(I)+COSIK
      IF(ISKIP.EQ.0)  GO TO 50
      DO 48  J=1,3
   48 DADXO(J,I)=DADXO(J,I)+HT(J)*SINIK
   50 CONTINUE
      IF(ISYM.EQ.1.OR.ISYM.EQ.3)  GO TO 60
      DBDAO(I)=DBDAO(I)+SINIK
      IF(ISKIP.EQ.0)  GO TO 60
      DO 58  J=1,3
   58 DBDXO(J,I)=DBDXO(J,I)+HT(J)*COSIK
   60 CONTINUE
C***  RESET TRANSFORMED INDICES
      TT=TTO
      DO 61 J=1,3
   61 HT(J)=HTO(J)
   65 CONTINUE
C***  END LOOP THROUGH NA ATOMS
   70 CONTINUE
C***  END LOOP THROUGH THE NS SYMMETRY TRANSFORMATIONS
      DO 75  IT=1,NT
   75 TF(IT)=EXP(-RHO*BI(IT))
      DO 80  IA=1,NF
      SR(IA)=SF(IA,ISA+1,KDT)+FSA*(SF(IA,ISA+2,KDT)-SF(IA,ISA+1,KDT))
      IF(ISA.EQ.0) Q=(-SF(IA,1,KDT)+2.0*SF(IA,2,KDT)-SF(IA,3,KDT))/2.0
      IF(ISA.NE.0) Q=(-SF(IA,ISA,KDT)+SF(IA,ISA+1,KDT)+SF(IA,ISA+2,KDT)-
     1 SF(IA,ISA+3,KDT))/4.0
      SR(IA)=SR(IA)+Q*FSA*(1.0-FSA)
      SI(IA)=0.0
      IF(ISYM.LT.3)  GO TO 80
      SR(IA)=SR(IA)+DR(IA,KDT)
      SI(IA)=DI(IA,KDT)
   80 CONTINUE
      DO 110  I=1,NA
      IA=ISF(I)   
      IT=ITF(I)
      TFSR=SR(IA)*TF(IT)
      TFSI=SI(IA)*TF(IT)
      DADA(I)=TFSR*DADAO(I)-TFSI*DBDAO(I)
      DADB(IT)=DADB(IT)+AI(I)*DADA(I)
      IF(ISKIP.EQ.0)  GO TO 95
      DO 90  J=1,3
   90 DADX(J,I)=AI(I)*(TFSR*DADXO(J,I)+TFSI*DBDXO(J,I))
   95 CONTINUE
      IF(ISYM.EQ.1)  GO TO 110
      DBDA(I)=TFSR*DBDAO(I)+TFSI*DADAO(I)
      DBDB(IT)=DBDB(IT)+AI(I)*DBDA(I)
      IF(ISKIP.EQ.0)  GO TO 110
      DO 100  J=1,3
  100 DBDX(J,I)=AI(I)*(TFSR*DBDXO(J,I)-TFSI*DADXO(J,I))
  110 CONTINUE
C***  END LOOP THROUGH THE NA ATOMS
      DO 142  IT=1,NT
      A=A+DADB(IT)
      IF(ISYM.EQ.1)  GO TO 142
      B=B+DBDB(IT)
  142 CONTINUE
      FMULT=SC(KDT)*M*PL/1000000
      IF(BO(KDT).NE.0.0)  FMULT=FMULT*EXP(-2.0*BO(KDT)*RHO)
      IF(ISYM.EQ.1.OR.ISYM.EQ.3)  FMULT=4.0*FMULT
      RINT=FMULT*(A**2+B**2)
      FMULT=2.0*FMULT
      IF(IC.LT.NC)  GO TO 146
      FCALC=EXP(-BO(KDT)*RHO)*SQRT(A**2+B**2)
      IF(ISYM.EQ.1.OR.ISYM.EQ.3) FCALC=2.0*FCALC
      WRITE(11,1100) (H(J),J=1,3),FCALC,A,B
 1100 FORMAT(1X,3I4,3F10.2)
      IF(IPR.NE.0) WRITE(61,1000)(H(J),J=1,3),M,PL,SINANG,RINT
      IF(I10.NE.0) WRITE(10,1000)(H(J),J=1,3),M,PL,SINANG,RINT
 1000 FORMAT(1X,3I4,1I6,1F11.3,1F11.6,1F11.2,64X)
  146 CONTINUE
C***  ACCUMULATE INTENSITY AND DERIVATIVE SUMS
      SINT=SINT+RINT
C***  BYPASS DERIVATIVES ON FINAL CYCLE
      IF(ISKIP.EQ.0)  GO TO 160
      DIDBO(KDT)=DIDBO(KDT)-2.0*RHO*RINT
      DO 147  I=1,NT
  147 DIDB(I)=DIDB(I)-RHO*FMULT*(A*DADB(I)+B*DBDB(I))
      DO 150  I=1,NA
      DIDA(I)=DIDA(I)+FMULT*(A*DADA(I)+B*DBDA(I))
      DO 148  J=1,3
  148 DIDX(J,I)=DIDX(J,I)-6.28318531*FMULT*(A*DADX(J,I)-B*DBDX(J,I))
  150 CONTINUE
  160 CONTINUE
      RETURN
      END
      SUBROUTINE SMI(AT,N,INDIC,PROW,PCOL)
C     **********SYMMETRIC  MATRIX  INVERSION**********
      REAL*8 AT,AMAX,SINV,PROW,PCOL
      LOGICAL  ELIM
      DIMENSION ELIM(99)
      DIMENSION AT(1),PROW(1),PCOL(1)
      NMI = N-1
      DO 1  I = 1,N
    1 ELIM(I)=.FALSE.
C     **********MAIN  LOOP**********
      DO 7  I=1,N
C     **********PIVOT  SEARCH**********
      AMAX = 0.0
      LOC = 1
         DO 2  NROW=1,N
         IF(AMAX.GE.ABS(AT(LOC)).OR.ELIM(NROW))  GO TO 2
         NPIV = NROW
         LPIV = LOC
         AMAX = ABS(AT(LOC))
    2    LOC = LOC+NROW+1
C     **********EXIT  FOR  SINGULAR  MATRIX**********
      IF(AMAX.EQ.0.0)  GO TO 8
      ELIM(NPIV) = .TRUE.
      SINV = 1.0/AT(LPIV)
C     **********PREPARATION  FOR  EXCHANGE  STEP**********
      IF(NPIV.EQ.N)  GO TO 4
      LOC = LPIV+NPIV
         DO 3  NROW=NPIV,NMI
         PROW(NROW+1) = AT(LOC)
         PCOL(NROW+1) = -AT(LOC)*SINV
         IF(ELIM(NROW+1))  PCOL(NROW+1) = -PCOL(NROW+1)
         AT(LOC) = 0.0
    3    LOC = LOC+NROW+1
    4 IF(NPIV.EQ.1)  GO TO 6
      LOC = LPIV-NPIV+1
      NPIVM = NPIV-1
         DO 5  NROW=1,NPIVM
         PROW(NROW) = AT(LOC)
         IF(ELIM(NROW))  PROW(NROW) = -PROW(NROW)
         PCOL(NROW) = -AT(LOC)*SINV
         AT(LOC) = 0.0
    5    LOC = LOC+1
    6 AT(LPIV) = 0.0
      PROW(NPIV) = 1.0
      PCOL(NPIV) = SINV
C     **********EXCHANGE  STEP**********
      LOC = 1
         DO 7  NROW=1,N
            DO 7  NCOL =1,NROW
            AT(LOC) = AT(LOC)+PROW(NROW)*PCOL(NCOL)
    7    LOC = LOC+1
C     **********END  MAIN  LOOP**********
      INDIC = 0
      RETURN
C     **********RETURN  FOR  SINGULAR  MATRIX**********
    8 INDIC = 1
      RETURN
      END
      SUBROUTINE COPEXP(R,SC,BO,BI,AJ,Z,NA,NT)
      DIMENSION BI(9),R(2),AJ(1),Z(3,1),SC(2),BO(2)
      SC(1)=R(1)
      SC(2)=R(2)
      BO(1)=R(3)
      BO(2)=R(4)
      DO 10  I=1,NT
   10 BI(I)=R(I+4)
      DO 30  I=1,NA
      AJ(I)=R(4*I+NT+1)
      DO 20  J=1,3
   20 Z(J,I)=R(4*I+NT+J+1)
   30 CONTINUE
      RETURN
      END
      SUBROUTINE COPTOP(R,SC,BO,BI,AJ,Z,NA,NT)
      DIMENSION BI(9),R(2),AJ(1),Z(3,1),SC(2),BO(2)
      R(1)=SC(1)
      R(2)=SC(2)
      R(3)=BO(1)
      R(4)=BO(2)
      DO 10  I=1,NT
   10 R(I+4)=BI(I)
      DO 30  I=1,NA
      R(4*I+NT+1)=AJ(I)
      DO 20  J=1,3
   20 R(4*I+NT+J+1)=Z(J,I)
   30 CONTINUE
      RETURN
      END
