      PROGRAM STRUCFAC
C  THIS CALCULATES STRUCTURE FACTORS FOR ALL THE 17 PLANE GROUPS.
C  IT WILL DO SO BOTH FOR NORMAL AND FOR ANOMALOUS SCATTERING.
C  THE INITIAL PART OF THE PROGRAMME REQUIRES THE INPUT OF THE 
C  FOLLOWING INFORMATION:
C  X-RAY WAVELENGTH - XLAM
C  CELL DIMENSIONS AND ANGLE - A, B, GAMMA
C    NOTE THAT OTHER PROJECTIONS CAN BE USED IN WHICH CASE 
C    A, B, GAMMA MAY, FOR EXAMPLE, BE REPLACED BY A, C, BETA SO 
C    THE INSTRUCTIONS GIVEN BY THE PROGRAMME WILL HAVE TO BE 
C    INTERPRETED ACCORDINGLY.    
C  PLANE GROUP - PG
C  CONTENTS OF THE CELL IN THE FORM OF:
C    TYPES OF ATOMS - ATYPE(I)
C    NUMBER OF ATOMS OF EACH TYPE IN THE ASYMMETRIC UNIT - NATOM(I)
C    THEIR COORDINATES - X, Y
C    WHETHER OR NOT THIS TYPE OF ATOM SCATTERS ANOMALOUSLY
C    IF SO THE VALUES OF F' AND F".
C  THE MAXIMUM NUMBER OF ATOMS IN THE UNIT CELL IS LIMITED TO 100 BY
C  THE DIMENSION OF THE ARRAYS X &amp; Y.   THE NUMBER OF DIFFERENT TYPES
C  OF ATOM IS LIMITED TO 20 BY THE DIMENSIONS OF IBEG, IEND, FP AND 
C  FPP AND THE MAXIMUM NUMBER OF REFLECTIONS FOR WHICH STRUCTURE 
C  FACTORS ARE CALCULATED IS LIMITED TO 1000 BY THE DIMENSIONS OF
C  THE ARRAYS HT, KT, F &amp; PHI.
      DIMENSION EQP(17,72),X(100),Y(100),F(1000),PHI(1000),FP(20),
     +FPP(20),NAU(17),MTYPE(13),NATOM(20),IBEG(20),IEND(20),
     +HT(1000),KT(1000)
      INTEGER HMAX,H,HT
      CHARACTER PGG*(*),PG*4,ATOM*(*),ASYMB*3,Q*4,FNAME*10
      PARAMETER (PGG='p1  p2  pm  pg  cm  pmm pmg pgg cmm p4  p4m p4g p3
     +  p3m1p31mp6  p6m ') 
      PARAMETER (ATOM='H  LI C  N  O  F  NA AL P  S  CL FE BR HG END')
C  THE FOLLOWING DATA STATEMENTS GIVE FOR ALL THE 17 PLANE GROUPS
C  THE NUMBER OF ASYMMETRIC UNITS IN THE CELL AND CONSTANTS FOR
C  CALCULATING THE POSITIONS OF ATOMS EQUIVALENT TO (X, Y) IN THE
C  FORM (C*X+D*Y+E, F*X+G*Y+H).
      DATA NAU(1),(EQP(1,I),I=1,6)/1,1,0,0,0,1,0/
      DATA NAU(2),(EQP(2,I),I=1,12)/2,1,0,0,0,1,0,-1,0,0,0,-1,0/
      DATA NAU(3),(EQP(3,I),I=1,12)/2,1,0,0,0,1,0,-1,0,0,0,1,0/
      DATA NAU(4),(EQP(4,I),I=1,12)/2,1,0,0,0,1,0,-1,0,0,0,1,0.5/
      DATA NAU(5),(EQP(5,I),I=1,24)/4,1,0,0,0,1,0,-1,0,0,0,1,0,
     +1,0,0.5,0,1,0.5,-1,0,0.5,0,1,0.5/
      DATA NAU(6),(EQP(6,I),I=1,24)/4,1,0,0,0,1,0,-1,0,0,0,1,0,
     +-1,0,0,0,-1,0,1,0,0,0,-1,0/
      DATA NAU(7),(EQP(7,I),I=1,24)/4,1,0,0,0,1,0,-1,0,0,0,-1,0,
     +1,0,0.5,0,-1,0,-1,0,0.5,0,1,0/
      DATA NAU(8),(EQP(8,I),I=1,24)/4,1,0,0,0,1,0,-1,0,0,0,-1,0,
     +1,0,0.5,0,-1,0.5,-1,0,0.5,0,1,0.5/
      DATA NAU(9),(EQP(9,I),I=1,48)/8,1,0,0,0,1,0,-1,0,0,0,1,0,
     +-1,0,0,0,-1,0,1,0,0,0,-1,0,1,0,0.5,0,1,0.5,-1,0,0.5,0,1,0.5,
     +-1,0,0.5,0,-1,0.5,1,0,0.5,0,-1,0.5/
      DATA NAU(10),(EQP(10,I),I=1,24)/4,1,0,0,0,1,0,-1,0,0,0,-1,0,
     +0,1,0,-1,0,0,0,-1,0,1,0,0/
      DATA NAU(11),(EQP(11,I),I=1,48)/8,1,0,0,0,1,0,-1,0,0,0,-1,0,
     +0,1,0,-1,0,0,0,-1,0,1,0,0,-1,0,0,0,1,0,1,0,0,0,-1,0,
     +0,-1,0,-1,0,0,0,1,0,1,0,0/
      DATA NAU(12),(EQP(12,I),I=1,48)/8,1,0,0,0,1,0,0,1,0,-1,0,0,
     +-1,0,0.5,0,1,0.5,0,-1,0.5,-1,0,0.5,-1,0,0,0,-1,0,0,-1,0,1,0,0,
     +1,0,0.5,0,-1,0.5,0,1,0.5,1,0,0.5/
      DATA NAU(13),(EQP(13,I),I=1,18)/3,1,0,0,0,1,0,0,-1,0,1,-1,0,
     +-1,1,0,-1,0,0/
      DATA NAU(14),(EQP(14,I),I=1,36)/6,1,0,0,0,1,0,0,-1,0,1,-1,0,
     +-1,1,0,-1,0,0,0,-1,0,-1,0,0,1,0,0,1,-1,0,-1,1,0,0,1,0/
      DATA NAU(15),(EQP(15,I),I=1,36)/6,1,0,0,0,1,0,0,-1,0,1,-1,0,
     +-1,1,0,-1,0,0,0,1,0,1,0,0,-1,0,0,-1,1,0,1,-1,0,0,-1,0/
      DATA NAU(16),(EQP(16,I),I=1,36)/6,1,0,0,0,1,0,0,-1,0,1,-1,0,
     +-1,1,0,-1,0,0,-1,0,0,0,-1,0,0,1,0,-1,1,0,1,-1,0,1,0,0/
      DATA NAU(17),(EQP(17,I),I=1,72)/12,1,0,0,0,1,0,0,-1,0,1,-1,0,
     +-1,1,0,-1,0,0,0,1,0,1,0,0,-1,0,0,-1,1,0,1,-1,0,0,-1,0,
     +-1,0,0,0,-1,0,0,1,0,-1,1,0,1,-1,0,1,0,0,0,-1,0,-1,0,0,
     +1,0,0,1,-1,0,-1,1,0,0,1,0/
C
      OPEN (UNIT=9, FILE='LPT1')
C  CHANGE NAME OF OUTPUT FILE IF REQUIRED
      FNAME='SF1.DAT'
   68 WRITE(6,'('' THE OUTPUT FILE, WHICH CONTAINS H K F AND PHI '')')
      WRITE(6,'('' IS NAMED "SF1.DAT". DO YOU WANT TO CHANGE IT Y/N'')')
      READ(5,50)Q
      IF(Q.EQ.'Y'.OR.Q.EQ.'y')GOTO 66
      IF(Q.EQ.'N')GOTO 67
      IF(Q.EQ.'n')GOTO 67
      GOTO 68
   66 WRITE(6,'('' READ IN THE FILENAME <= 10 CHARACTERS] '')')
      READ(5,58)FNAME
   58 FORMAT(A10)
   67 OPEN (UNIT=10, FILE=FNAME)
      WRITE(6,'('' INPUT X-RAY WAVELENGTH IN ANGSTROMS '')')
      READ(5,*)XLAM
      WRITE(6,'('' INPUT CELL DIMENSIONS A,B [IN ANGSTROMS '')')
      WRITE(6,'('' AND GAMMA IN DEGREES '')')
      READ(5,*)A,B,GAMMA
C  READ IN PLANE GROUP SYMBOL AND TRANSLATE INTO IDENTIFYING INTEGER
    1 WRITE(6,'('' INPUT PLANE GROUP - ONE OF THE FOLLOWING '')')
      WRITE(6,'('' p1, p2, pm, pg, pgg, cm, pmm, pmg, cmm, p4,'')')
      WRITE(6,'('' p4m, p4g, p3, p3m1, p31m, p6, p6m '')')
      WRITE(6,'('' Note: use lower case symbols as shown '')')
      READ(5,50)PG
   50 FORMAT(A4)
C  FORM IDENTIFYING INTEGER 1 TO 17 FOR PLANE GROUP
      NX=INDEX(PGG,PG)+3
      NPG=NX/4      
      IF(NPG*4.NE.NX)GOTO 1
C  NOW INITIALISE THE FOLLOWING QUANTITIES: NTOTAT -TOTAL NUMBER OF 
C  ATOMS IN THE CELL; NOTYPE - THE NUMBER OF DIFFERENT TYPE OF ATOMS
C  IBEG AND IEND ARE ARRAYS WHICH FLAG THE FIRST AND LAST ATOMS OF
C  THE VARIOUS TYPES IN THE COORDINATE LIST
      NTOTAT=0
      NOTYPE=0
      IBEG(1)=1
    2 WRITE(6,'('' INPUT ATOMIC TYPE-ONE OF H, LI, C, N, O, F, NA,'')')
      WRITE(6,'('' AL, P, S, CL, FE, BR, HG - OR INPUT END IF ALL '')')
      WRITE(6,'('' COORDINATES HAVE BEEN READ IN '')')
      READ(5,60)ASYMB
   60 FORMAT(A3)
C  GENERATE A NUMBER FROM 1 TO 14 TO IDENTIFY THE TYPE OF ATOM WHERE
C  1 = H, 2 = LI, ......... 14 = HG, OR 15 TO END INPUT.
      NX=INDEX(ATOM,ASYMB)+2
      NASYMB=NX/3      
      IF(NASYMB*3.NE.NX)GOTO 2
      IF(NASYMB.EQ.15)GOTO 33      
C  A NEW TYPE OF ATOM HAS BEEN IDENTIFIED. 
      NOTYPE=NOTYPE+1
      MTYPE(NOTYPE)=NASYMB
      WRITE(6,'('' INPUT THE NUMBER OF ATOMS OF THIS TYPE IN ONE '')')
      WRITE(6,'('' ASYMMETRIC UNIT '')')
      READ(5,*)NATOM(NOTYPE)
C  INPUT COORDINATES
      WRITE(6,'('' INPUT COORDINATES X, Y - ONE ATOM PER LINE '')')
      WRITE(6,'('' IF YOU ENTER A COORDINATE INCORRECTLY THERE '')')
      WRITE(6,'('' IS AN OPPORTUNITY LATER TO CORRECT IT '')')
      DO 5 I=1,NATOM(NOTYPE)
      READ(5,*)XX,YY
C  GENERATE COORDINATES OF ALL SYMMETRY RELATED ATOMS
      DO 5 J=1,NAU(NPG)
      NTOTAT=NTOTAT+1
      X(NTOTAT)=EQP(NPG,6*J-5)*XX+EQP(NPG,6*J-4)*YY+EQP(NPG,6*J-3)
      Y(NTOTAT)=EQP(NPG,6*J-2)*XX+EQP(NPG,6*J-1)*YY+EQP(NPG,6*J)
    5 CONTINUE
      IEND(NOTYPE)=NTOTAT
      IBEG(NOTYPE+1)=NTOTAT+1
C  IF THIS ATOMIC TYPE IS AN ANOMALOUS SCATTERER THEN THE REAL AND
C  IMAGINARY PARTS OF THE ANOMALOUS COMPONENT OF THE SCATTERING 
C  FACTOR MUST BE GIVEN
   16 WRITE(6,'('' IS THIS TYPE AN ANOMALOUS SCATTERER? Y/N'')')
      READ(5,50)Q
      IF(Q.EQ.'Y   '.OR.Q.EQ.'y   ')GOTO 6
      IF(Q.EQ.'N   ')GOTO 2
      IF(Q.EQ.'n   ')GOTO 2
      GOTO 16
    6 WRITE(6,'('' INPUT REAL AND IMAGINARY PART OF A.S FACTOR'')')
      READ(5,*)FP(NOTYPE),FPP(NOTYPE)
      GOTO 2
C  AT THIS STAGE ALL ATOMIC COORDINATES HAVE BEEN INPUT
C  IN CASE OF ERROR THEY MAY NOW BE DISPAYED WITH A POSSIBILITY
C  OF CHANGING ANY INCORRECT ONE
   33 WRITE(6,'('' DO YOU WANT TO INSPECT THE COORDINATES? Y/N '')')
      READ(5,50)Q
      IF(Q.EQ.'Y'.OR.Q.EQ.'y')GOTO 80
      IF(Q.EQ.'N')GOTO 3
      IF(Q.EQ.'n')GOTO 3
      GOTO 33
   80 DO 90 I=1,NTOTAT,NAU(NPG)
      WRITE(6,200)X(I),Y(I)
  200 FORMAT(2F8.4)
   70 WRITE(6,'('' DO YOU WANT TO CHANGE THIS? Y/N'')')
      READ(5,50)Q
      IF(Q.EQ.'Y'.OR.Q.EQ.'y')GOTO 71
      IF(Q.EQ.'N')GOTO 90
      IF(Q.EQ.'n')GOTO 90
      GOTO 70
   71 WRITE(6,'('' INPUT REPLACEMENT X,Y '')')
      READ(5,*)XX,YY
      DO 72 J=1,NAU(NPG)
      X(I+J-1)=EQP(NPG,6*J-5)*XX+EQP(NPG,6*J-4)*YY+EQP(NPG,6*J-3)
   72 Y(I+J-1)=EQP(NPG,6*J-2)*XX+EQP(NPG,6*J-1)*YY+EQP(NPG,6*J)
   90 CONTINUE     
    3 WRITE(6,'('' INPUT TEMPERATURE FACTOR.   IF ONE IS NOT '')')
      WRITE(6,'('' REQUIRED THEN INPUT 0 '')')
      READ(5,*)BTEMP
C  PROVIDE A TEST QUANTITY FOR THE PRESENCE OF ANOMALOUS SCATTERING
      TESTAS=0
      DO 7 I=1,NOTYPE
    7 TESTAS=TESTAS+FP(I)**2+FPP(I)**2
C  ALL INPUT IS NOW COMPLETE.   NOW STRUCTURE FACTORS ARE CALCULATED
C  TO THE LIMIT OF RESOLUTION GIVEN BY THE X-RAY WAVELENGTH.  
C  NORMALLY A COMPLETE SEMICIRCLE OF STRUCTURE FACTORS WILL BE 
C  CALCULATED ALTHOUGH THE ASYMMETRIC UNIT OF RECIPROCAL SPACE 
C  MAY BE SMALLER THAN THIS.   CALCULATIONS ARE FOR ALL INDEPENDENT 
C  REFLECTIONS  WITH k &gt;= 0.   IF ANOMALOUS SCATTERING IS PRESENT THEN
C  THE COMPLETE CIRCLE IS CALCULATED.  FIRST CALCULATE MAXIMUM INDICES.
      HMAX=INT(2.0*A/XLAM)
      KMAX=INT(2.0*B/XLAM)
C  CO1, CO2 AND CO3 ARE USED FOR CALCULATING SIN(THETA)/LAMBDA
      CO1=0.25/A/A
      CO2=0.25/B/B
      CO3=0.5/A/B*COS(GAMMA*0.0174533)
C  SLIM IS THE LIMITING VALUE OF (SIN(THETA)/LAMBDA)**2 AND NOREF 
C  RECORDS THE TOTAL NUMBER OF REFLECTIONS.
      SLIM=1/XLAM/XLAM
      NOREF=0
      DO 10 H=-HMAX,HMAX
      DO 10 K=0,KMAX
      IF(K.EQ.0.AND.H.LT.0)GOTO 10
C  NOW CALCULATE (SIN(THETA)/LAMBDA)**2 AND CHECK THAT &lt; SLIM
      SL2=CO1*H*H+CO2*K*K-CO3*H*K
      IF(SL2.GT.SLIM)GOTO 10
      SL=SQRT(SL2)
C  NOW THE STRUCTURE FACTOR IS CALCULATED.   IF THERE IS ANOMALOUS 
C  SCATTERING THEN BOTH F(h k l) AND F(-h -k -l) ARE CALCULATED.    
      SUMA=0
      SUMB=0
      SMAASP=0
      SMAASM=0
      SMBASP=0
      SMBASM=0
      DO 15 I=1,NOTYPE              
C  CALL SUBROUTINE TO DERIVE SCATTERING FACTOR [SCF] FOR ATOM 
C  OF TYPE I FOR SIN(THETA)/LAMBDA = SL
      CALL SCAT(MTYPE(I),SL,SCF)
      PSUMA=0
      PSUMB=0
C  CALCULATE CONTRIBUTIONS TO REAL AND IMAGINARY PARTS OF THE
C  STRUCTURE FACTOR OF ATOMS OF TYPE I
      DO 26 J=IBEG(I),IEND(I)
      Z=6.2831853*(H*X(J)+K*Y(J))
      PSUMA=PSUMA+COS(Z)
   26 PSUMB=PSUMB+SIN(Z)
      SUMA=SUMA+SCF*PSUMA
      SUMB=SUMB+SCF*PSUMB
C  TEST FOR PRESENCE OF ANOMALOUS SCATTERING
      IF(TESTAS.LT.1.0E-12)GOTO 15
C  TEST FOR ATOM TYPE I AS ANOMALOUS SCATTERER
      IF(FP(I)**2+FPP(I)**2.LT.1.0E-12)GOTO 15
      SMAASP=SMAASP+PSUMA*FP(I)-PSUMB*FPP(I)
      SMAASM=SMAASM+PSUMA*FP(I)+PSUMB*FPP(I)
      SMBASP=SMBASP+PSUMB*FP(I)+PSUMA*FPP(I)
      SMBASM=SMBASM-PSUMB*FP(I)+PSUMA*FPP(I)
   15 CONTINUE
      NOREF=NOREF+1
      HT(NOREF)=H
      KT(NOREF)=K
      BTFAC=1.0
      IF(BTEMP.GT.1.0E-12)BTFAC=EXP(-BTEMP*SL2)
      F(NOREF)=SQRT((SUMA+SMAASP)**2+(SUMB+SMBASP)**2)*BTFAC
      PHI(NOREF)=ATAN2(SUMB+SMBASP,SUMA+SMAASP)*57.2958
      IF(F(NOREF).LT.1.0E-4)PHI(NOREF)=0
      IF(TESTAS.LT.1.0E-12)GOTO 10
      NOREF=NOREF+1
      HT(NOREF)=-H
      KT(NOREF)=-K
      F(NOREF)=SQRT((SUMA+SMAASM)**2+(-SUMB+SMBASM)**2)*BTFAC
      PHI(NOREF)=ATAN2(-SUMB+SMBASM,SUMA+SMAASM)*57.2958
      IF(F(NOREF).LT.1.0E-4)PHI(NOREF)=0
   10 CONTINUE           
      WRITE(10,*)NOREF,(HT(I),KT(I),F(I),PHI(I),I=1,NOREF)
      WRITE(9,700)
  700 FORMAT(3(24H    h    k     F    PHI ))
      WRITE(9,100)(HT(I),KT(I),F(I),PHI(I),I=1,NOREF)
  100 FORMAT(3(I6,I4,F8.2,F6.0))
      STOP
      END


      SUBROUTINE SCAT(I,SL,SCF)
C  THERE NOW FOLLOW THE SCATTERING FACTORS FROM INTERNATIONAL TABLES
C  FOR ATOMS H, LI, C, N, O, NA, AL, P, S, CL, FE, BR AND HG AT 
C  INTERVALS OF 0.1 FOR SIN(THETA)/LAMBDA. PARABOLIC INTERPOLATION IS
C  USED TO OBTAIN THE SCATTERING FACTOR FOR SIN(THETA)/LAMBDA = SL.
      DIMENSION SF(14,0:10)
      DATA (SF(1,I),I=0,10)/1.000,0.811,0.481,0.251,0.130,0.071,
     +0.040,0.024,0.015,0.010,0.007/
      DATA (SF(2,I),I=0,10)/3.000,2.215,1.741,1.521,1.269,1.032,
     +0.823,0.650,0.513,0.404,0.320/
      DATA (SF(3,I),I=0,10)/6.000,5.126,3.581,2.502,1.950,1.685,
     +1.536,1.426,1.322,1.218,1.114/
      DATA (SF(4,I),I=0,10)/7.000,6.203,4.600,3.241,2.397,1.944,
     +1.698,1.550,1.444,1.350,1.263/
      DATA (SF(5,I),I=0,10)/8.000,7.250,5.634,4.094,3.010,2.338,
     +1.944,1.714,1.566,1.462,1.374/
      DATA (SF(6,I),I=0,10)/9.000,8.293,6.691,5.044,3.760,2.878,
     +2.312,1.958,1.735,1.587,1.481/
      DATA (SF(7,I),I=0,10)/11.00,9.76,8.34,6.89,5.47,4.29,3.40,
     +2.76,2.31,2.00,1.78/
      DATA (SF(8,I),I=0,10)/13.00,11.23,9.16,7.88,6.72,5.69,4.71,
     +3.88,3.21,2.71,2.32/
      DATA (SF(9,I),I=0,10)/15.00,13.17,10.34,8.59,7.54,6.67,5.83,
     +5.02,4.28,3.64,3.11/
      DATA (SF(10,I),I=0,10)/16.00,14.33,11.21,8.99,7.83,7.05,6.31,
     +5.56,4.82,4.15,3.56/
      DATA (SF(11,I),I=0,10)/17.00,15.33,12.00,9.44,8.07,7.29,6.64,
     +5.96,5.27,4.60,4.00/
      DATA (SF(12,I),I=0,10)/26.00,23.68,20.09,16.77,13.84,11.47,
     +9.71,8.47,7.60,6.99,6.51/
      DATA (SF(13,I),I=0,10)/35.00,32.43,27.70,23.82,20.84,18.27,
     +15.91,13.78,11.93,10.41,9.19/
      DATA (SF(14,I),I=0,10)/80.00,75.48,67.14,59.31,52.65,47.04,
     +42.31,38.22,34.64,31.43,28.59/
      J=INT(SL/0.1)
      IF(J.EQ.10)J=9
      IF(J.EQ.0)J=1
      A=(SF(I,J+1)+SF(I,J-1)-2*SF(I,J))/0.02
      B=(SF(I,J+1)-SF(I,J-1))/0.2
      SCF=A*(SL-0.1*J)**2+B*(SL-0.1*J)+SF(I,J)
      RETURN
      END
      
