      SUBROUTINE SFAC
C
C THIS ROUTINE REVISED FEB. 10, 1972
C
C             STRUCTURE FACTOR GENERATION
C
      INCLUDE 'RFINE.PRM'
      INCLUDE 'A.INC'
      INCLUDE 'C.INC'
      INCLUDE 'D.INC'
      INCLUDE 'CNSTRS.INC'
      INCLUDE 'E.INC'
      INCLUDE 'F.INC'
      INCLUDE 'LS.INC'
      INCLUDE 'LSMAT.INC'
      INCLUDE 'TWIN.INC'
      DIMENSION V(3),VV(6),HSAVE(3),HSTANG(4),GSTANG(4),HTEMP(3)
      EQUIVALENCE (V,DAX),(VV,HHJ)
C
      DATA TWOPI/6.283185307/
C
      ISAVE=MFSEL
      IF(MFSEL.NE.0.AND.(IREJ.EQ.1.OR.IREJ.EQ.2).AND.MREJ.EQ.1)MFSEL=0
	ISAVE2=MFSEL
C SAVE INDICES
	DO 300 I=1,3
	  HSAVE(I)=H(I)
  300	CONTINUE
C BEGIN TWIN CALCULATION -- GENERATE HOST ANGLES IF NEEDED
      IF(NT.GT.0)THEN
	CALL GETANG(HSAVE,UBMAT,HSTANG,WAVE,PHI0,0)
	PHI0=HSTANG(4)
	SINTH=SIN(HSTANG(1)/114.59)
      END IF
C GENERATE TWIN INDICES
      DO 970 NTWN=NT,0,-1
	IF(NTWN.GT.0)THEN
	  ATWN(NTWN)=0.0
	  BTWN(NTWN)=0.0
	  MFSEL=0
	  DO 320 I=1,3
	    H(I)=0.0
	    DO 310 J=1,3
		H(I)=H(I)+TWINTM(I,J,NTWN)*HSAVE(J)
  310	    CONTINUE
	    H(I)=ANINT(H(I))
  320	  CONTINUE
C BACK TRANSFORM TWIN INDICES
	  DO 322 I=1,3
	    HTEMP(I)=0.0
	    DO 321 J=1,3
		HTEMP(I)=HTEMP(I)+TWINTM(J,I,NTWN)*H(J)
  321	    CONTINUE
  322	  CONTINUE
C GENERATE GUEST ANGLES
	  CALL GETANG(HTEMP,UBMAT,GSTANG,WAVE,PHI0,1)
C CHECK IF GUEST ANGLES OUTSIDE WINDOWS
	  DO 323 I=1,3
	    TEMP=ABS(GSTANG(I)-HSTANG(I))
C CONVERT DELTA CHI TO RADIANS AND CORRECT FOR SIN(THETA) EFFECT
	    IF(I.EQ.3)TEMP=TEMP*SINTH/28.648
	    IF(TEMP.GT.XWIN(I))GO TO 960
  323	  CONTINUE
	ELSE
	  DO 330 I=1,3
	    H(I)=HSAVE(I)
  330	  CONTINUE
	END IF
      ACALC=0.0
      BCALC=0.0
      AMINUS=0.0
      BMINUS=0.0
C CHECK IF REFLECTION ALLOWED BY LATTICE
      DO 338 L=1,JCELL
        SUM=0.0
        DO 336 K=1,3
          SUM=SUM+CELLP(K,L)*H(K)
336     CONTINUE
C IF SUM NOT AN INTEGER, REFLECTION NOT ALLOWED
        IF(ABS(ANINT(SUM)-SUM).GT.0.05)GO TO 952
338   CONTINUE
C             GENERATE TRANSFORMED INDICES
      DO 340 I=1,NSYM
        TPR(I)=0.0
        DO 340 J=1,3
          HP(J,I)=0.0
          TPR(I)=TPR(I)+H(J)*TRANS(J,I)
          DO 340 K=1,3
            HP(J,I)=HP(J,I)+H(K)*RMAT(K,J,I)
340   CONTINUE
C             START LOOP THROUGH ATOMS
      DO 950 L=1,NATOM
C             INITIALIZE
      TF=JCELL*(2-ICENT)
      EXPIJ=1.0
      AAIND=0.0
      BBIND=0.0
      IF(MFSEL)350,430,350
  350                        DO 360 M=1,3
      DAX(M)=0.0
  360 DBX(M)=0.0
      IF(ISOT(L))370,440,370
  370                        DO 380 N=1,6
      DAB(N)=0.0
  380 DBB(N)=0.0
      IF(ISOT(L)-1)450,450,390
  390 DO 400 N=1,10
      DAC(N)=0.0
  400 DBC(N)=0.0
      IF(ISOT(L)-2)450,450,410
  410 DO 420 N=1,15
      DAD(N)=0.0
  420 DBD(N)=0.0
      GO TO 450
  430 IF(ISOT(L))450,440,450
C             CALCULATE ISOTROPIC TEMPERATURE FACTOR
  440 TF=EXP(-BETA(1,L)*RHO)*TF
C             START LOOP THROUGH SYMMETRY OPERATIONS
  450                        DO 790 IA=1,NSYM
      ARG=TPR(IA)
                             DO 460 IB=1,3
      HPR(IB)=HP(IB,IA)
  460 ARG=ARG+HPR(IB)*XYZ(IB,L)
      CUM3=0.0
      IF(ISOT(L))470,570,470
C             GENERATE TRANSFORMED INDEX PRODUCTS FOR ANISO. ATOMS
  470                        DO 480 IC=1,3
  480 HHJ(IC)=HPR(IC)**2
      HHJ(4)=2.0*HPR(1)*HPR(2)
      HHJ(5)=2.0*HPR(1)*HPR(3)
      HHJ(6)=2.0*HPR(2)*HPR(3)
C             CALCULATE ANISOTROPIC TEMPERATURE FACTOR
      TEMP=0.0
                             DO 490 ID=1,6
  490 TEMP=TEMP-HHJ(ID)*BETA(ID,L)
      CUM4=0.0
      IF(ISOT(L)-1)560,560,500
C        CALCULATE THIRD CUMULANT INDEX PRODUCTS
  500 DO 510 ID=1,3
  510 HHHJ(ID)=HPR(ID)*HHJ(ID)
      HHHJ(4)=3.0*HHJ(1)*HPR(2)
      HHHJ(5)=3.0*HHJ(2)*HPR(1)
      HHHJ(6)=3.0*HHJ(1)*HPR(3)
      HHHJ(7)=3.0*HHJ(3)*HPR(1)
      HHHJ(8)=3.0*HHJ(2)*HPR(3)
      HHHJ(9)=3.0*HHJ(3)*HPR(2)
      HHHJ(10)=3.0*HHJ(4)*HPR(3)
      DO 520 ID=1,10
  520 CUM3=CUM3+HHHJ(ID)*CR(ID,L)
      IF(ISOT(L)-2)560,560,530
C        CALCULATE FOURTH CUMULANT INDEX PRODUCTS
  530 DO 540 ID=1,3
  540 HHHH(ID)=HHJ(ID)**2
      HHHH( 4)= 4.0*HHHJ(1)*HPR(2)
      HHHH( 5)= 4.0*HHHJ(1)*HPR(3)
      HHHH( 6)= 6.0*HHJ(1)*HHJ(2)
      HHHH(7)=12.0*HHJ(1)*HPR(2)*HPR(3)
      HHHH(8)=6.0*HHJ(1)*HHJ(3)
      HHHH(9)=4.0*HHHJ(2)*HPR(1)
      HHHH(10)=12.0*HHJ(2)*HPR(1)*HPR(3)
      HHHH(11)=12.0*HHJ(3)*HPR(1)*HPR(2)
      HHHH(12)=4.0*HHHJ(3)*HPR(1)
      HHHH(13)=4.0*HHHJ(2)*HPR(3)
      HHHH(14)=6.0*HHJ(2)*HHJ(3)
      HHHH(15)=4.0*HHHJ(3)*HPR(2)
      DO 550 ID=1,15
  550 CUM4=CUM4+HHHH(ID)*DR(ID,L)
C        CALCULATE TEMPERATURE FACTOR
  560 EXPIJ=EXP(TEMP+CUM4)
  570 T=TWOPI*ARG-CUM3
      COSIJ=COS (T)*EXPIJ
      SINIJ=SIN (T)*EXPIJ
      AAIND=AAIND+COSIJ
      IF(ICENT)580,590,580
  580 BBIND=BBIND+SINIJ
  590 IF(MFSEL)600,790,600
C             CALCULATE POSITIONAL DERIVATIVES
  600                        DO 610 IAA=1,3
  610 DAX(IAA)=DAX(IAA)-HPR(IAA)*SINIJ
      IF(ICENT)620,640,620
  620                        DO 630 IAB=1,3
  630 DBX(IAB)=DBX(IAB)+HPR(IAB)*COSIJ
  640 IF(ISOT(L))650,790,650
C             CALCULATE THERMAL DERIVATIVES FOR ANISOTROPIC ATOMS
  650                        DO 660 IAC=1,6
  660 DAB(IAC)=DAB(IAC)-HHJ(IAC)*COSIJ
      IF(ICENT)670,690,670
  670 DO 680 ID=1,6
  680 DBB(ID)=DBB(ID)-HHJ(ID)*SINIJ
  690 IF(ISOT(L)-1)790,790,700
C        CALCULATE THIRD CUM. DERIVS.
  700 DO 710 ID=1,10
  710 DAC(ID)=DAC(ID)+HHHJ(ID)*SINIJ
      IF(ICENT)720,740,720
  720 DO 730 ID=1,10
  730 DBC(ID)=DBC(ID)-HHHJ(ID)*COSIJ
  740 IF(ISOT(L)-2)790,790,750
C        CALCULATE FOURTH CUMULANT DERIVS.
  750 DO 760 ID=1,15
  760 DAD(ID)=DAD(ID)+HHHH(ID)*COSIJ
      IF(ICENT)770,790,770
  770 DO 780 ID=1,15
  780 DBD(ID)=DBD(ID)+HHHH(ID)*SINIJ
  790 CONTINUE
  800 ITY1=ISCAT(1,L)
      ITY2=ISCAT(2,L)
      TEMP4=TF*SITE(L)
      IF(ITY2)810,820,810
C             EVALUATE SCATTERING FACTOR FOR MULTIPLE OCCUPANCY
  810 OCCB=OCCUP(L)-OCCA(L)
      REAL1=SCAT(ITY1)
      REAL2=SCAT(ITY2)
      TEMP=OCCA(L)*REAL1+OCCB*REAL2
      TEMP1=OCCA(L)*ANOM(2,ITY1)+OCCB*ANOM(2,ITY2)
      TEMP2=REAL1-REAL2
      TEMP3=ANOM(2,ITY1)-ANOM(2,ITY2)
      GO TO 830
C             EVALUATE SCATTERING FACTOR FOR SINGLE OCCUPANCY
  820 REAL1=SCAT(ITY1)
      TEMP=OCCUP(L)*REAL1
      TEMP3=ANOM(2,ITY1)
      TEMP1=OCCUP(L)*TEMP3
      TEMP2=REAL1
  830 BBIND=FISGN*BBIND
      IF(MFSEL)840,850,840
C             EVALUATE DERIVATIVES WITH RESPECT TO OCCUPANCY
  840 DADO(L)=TEMP4*(TEMP2*AAIND-TEMP3*BBIND)
      DBDO(L)=TEMP4*(TEMP2*BBIND+TEMP3*AAIND)
      IF(IFSGN.EQ.0)GO TO 850
C        EVALUATE STRUCTURE FACTOR FOR ENANTIOMORPHIC ORIENTATION
      AMINUS =TEMP4*(TEMP*AAIND+TEMP1*BBIND)+AMINUS
      BMINUS=TEMP4*(-TEMP*BBIND+TEMP1*AAIND)+BMINUS
C             EVALUATE INDIVIDUAL ATOM CONTRIBUTIONS
  850 AIND(L)=TEMP4*(TEMP*AAIND-TEMP1*BBIND)
      BIND(L)=TEMP4*(TEMP*BBIND+TEMP1*AAIND)
      ACALC=ACALC+AIND(L)
      BCALC=BCALC+BIND(L)
      IF(MFSEL)860,950,860
C             MODIFY DERIVATIVES
  860 TEMP5=TWOPI*TEMP4
                             DO 870 IAG=1,3
      DADX(IAG,L)=TEMP5*(TEMP*DAX(IAG)-TEMP1*DBX(IAG))
  870 DBDX(IAG,L)=TEMP5*(TEMP*DBX(IAG)+TEMP1*DAX(IAG))
      IF(ISOT(L))890,880,890
C             DERIVATIVE WITH RESPECT TO ISOTROPIC TEMPERATURE FACTOR
  880 DADB(1,L)=-RHO*AIND(L)
      DBDB(1,L)=-RHO*BIND(L)
      GO TO 950
C             DERIVATIVES WITH RESPECT TO ANISOTROPIC TEMPERATURE FACTOR
  890                        DO 900 IAH=1,6
      DADB(IAH,L)=TEMP4*(TEMP*DAB(IAH)-TEMP1*DBB(IAH))
  900 DBDB(IAH,L)=TEMP4*(TEMP*DBB(IAH)+TEMP1*DAB(IAH))
      IF(ISOT(L)-1)950,950,910
C        DERIVATIVES WRT THE THIRD CUMULANT
  910 DO 920 ID=1,10
      DADC(ID,L)=TEMP4*(TEMP*DAC(ID)-TEMP1*DBC(ID))
  920 DBDC(ID,L)=TEMP4*(TEMP*DBC(ID)+TEMP1*DAC(ID))
      IF(ISOT(L)-2)950,950,930
C        DERIVATIVES WRT THE FOURTH CUMULANT
  930 DO 940 ID=1,15
      DADD(ID,L)=TEMP4*(TEMP*DAD(ID)-TEMP1*DBD(ID))
  940 DBDD(ID,L)=TEMP4*(TEMP*DBD(ID)+TEMP1*DAD(ID))
  950 CONTINUE
  952 IF(NTWN.GT.0)THEN
        ATWN(NTWN)=ACALC
        BTWN(NTWN)=BCALC
      END IF
  960 MFSEL=ISAVE2
  970 CONTINUE
  980 MFSEL=ISAVE
      RETURN
      END
      SUBROUTINE GETANG(H,UB,ANG,WAVE,PHI0,IFLG)
C
C ROUTINE TO CALCULATE EULERIAN ANGLES FOR REFLECTION WITH INDICES 'H'
C  AND MATRIX 'UB'.  IF 'IFLG' IS ZERO, CALCULATE IN BISECTING MODE,
C  OTHERWISE CALCULATE ANGLES WITH PHI = 'PHI0'.
C  THIS VERSION ASSUMES THE BUSING AND LEVY ANGLE CONVENTIONS.
C
      DIMENSION H(3),UB(3,3),ANG(4),HPHI(3)
      DATA RAD/57.29578/
C
	SUM=0.0
	DO 110 I=1,3
	  HPHI(I)=0.0
	  DO 100 J=1,3
	    HPHI(I)=HPHI(I)+UB(I,J)*H(J)
100	  CONTINUE
	  SUM=SUM+HPHI(I)**2
110	CONTINUE
	SUM=SQRT(SUM)
	TEMP=SUM*WAVE/2.0
	IF(ABS(TEMP).GT.1.0)TEMP=1.0
	ANG(1)=2.0*RAD*ASIN(TEMP)
	ANG(2)=0.0
	ANG(3)=RAD*ATAN2(HPHI(3),SQRT(HPHI(1)**2+HPHI(2)**2))
	if(abs(hphi(2)).gt.1.0e-15)then
	    ANG(4)=RAD*ATAN2(HPHI(2),HPHI(1))
	else
	    ang(4)=0.0
	endif
	IF(IFLG.EQ.0)RETURN
C SET PHI TO PHI0 AND SOLVE FOR DELTA CHI AND DELTA OMEGA
	ANG(4)=PHI0
	NTRY=0
120	R1=HPHI(1)/SUM-COS(ANG(2)/RAD)*COS(ANG(3)/RAD)*COS(ANG(4)/RAD)
     1  +SIN(ANG(2)/RAD)*SIN(ANG(4)/RAD)
	R2=HPHI(2)/SUM-COS(ANG(2)/RAD)*COS(ANG(3)/RAD)*SIN(ANG(4)/RAD)
     1  -SIN(ANG(2)/RAD)*COS(ANG(4)/RAD)
	A1 = -SIN(ANG(2)/RAD)*COS(ANG(3)/RAD)*COS(ANG(4)/RAD)
     1  -COS(ANG(2)/RAD)*SIN(ANG(4)/RAD)
	A2 = -SIN(ANG(2)/RAD)*COS(ANG(3)/RAD)*SIN(ANG(4)/RAD)
     1  +COS(ANG(2)/RAD)*COS(ANG(4)/RAD)
	B1 = -COS(ANG(2)/RAD)*SIN(ANG(3)/RAD)*COS(ANG(4)/RAD)
	B2 = -COS(ANG(2)/RAD)*SIN(ANG(3)/RAD)*SIN(ANG(4)/RAD)
	DEN = A1*B2-A2*B1
	IF(ABS(DEN).GT.1.0E-3)THEN
	  Y = (A2*R1-A1*R2)/DEN
	  X = (B2*R1-B1*R2)/DEN
	  ANG(2)=ANG(2)+X*RAD/2.0
	  ANG(3)=ANG(3)+Y*RAD/2.0
	  NTRY=NTRY+1
	  IF(AMAX1(ABS(X),ABS(Y)).GE.0.004.AND.NTRY.LT.5)GO TO 120
	END IF
	RETURN
	END
