PROGRAM SFACTR
* Need to define the data types before defining parameters.
* Makes the string STRUCTURE equal to 3 letters.
	CHARACTER*(3) STRUCTURE
	INTEGER H, K, L
	REAL PI, SREAL, SIMG
	COMPLEX S
* Parameter statement gives the value of pi the label PI.
	PARAMETER (PI = 3.1415926)
* Ask user to tell the program which structure they want the structure factor of.
	WRITE(UNIT=*, FMT=*) 'What structure is the lattice (BCC, FCC,
     $ SC, HCP)?'
	READ(UNIT=*, FMT=*) STRUCTURE
* Ask user to define the Miller indicies of the diffracted beam.
	WRITE(UNIT=*, FMT=*) 'Enter the Miller indicies of the 
     $ diffracted beam (h, k, l) '
	READ(UNIT=*, FMT=*) H, K, L
* Call the required subroutine.
	IF(STRUCTURE .EQ. 'BCC' .OR. STRUCTURE .EQ. 'bcc') THEN
		CALL BCCS(PI,H,K,L,S)
	ELSE IF(STRUCTURE .EQ. 'FCC' .OR. STRUCTURE .EQ. 'fcc') THEN
		CALL FCCS(PI,H,K,L,S)
	ELSE IF(STRUCTURE .EQ. 'SC' .OR. STRUCTURE .EQ. 'sc') THEN
		CALL SCS(PI,H,K,L,S)
	ELSE IF(STRUCTURE .EQ. 'HCP' .OR. STRUCTURE .EQ. 'hcp') THEN
		CALL HCPS(PI,H,K,L,S)
	ELSE
		STOP '! Incorrect structure specified !'
	END IF
* Find the real and imaginary terms in the structure factor.
	SREAL = REAL(S)
	SIMG = AIMAG(S)
* Return the structure factor.
	WRITE(UNIT=*, FMT=*) STRUCTURE
	WRITE(UNIT=*, FMT=* ) 'Structure factor = ( ', SREAL,
     $  ' + ', SIMG, 'i ) xf'
	WRITE(UNIT=*, FMT=*) 'Number format is a complex number'
	WRITE(UNIT=*, FMT=*) '(A + Bi) where A is the real term,'
	WRITE(UNIT=*, FMT=*) ' and B is the imaginary term' 
	END

	SUBROUTINE BCCS(PI,H,K,L,S)
* Define data types
	INTEGER H, K, L, MILLERS
	COMPLEX EXPONENT, S
	REAL PI, FACTOR, X, Y
* Sum Miller indicies and find exponential of pi*i*sum
	MILLERS = H + K + L
	FACTOR = PI * MILLERS
	X = COS(FACTOR)
	Y = SIN(FACTOR)
	EXPONENT = CMPLX(X, Y)
* Find the structure factor
	S = (1, 0) + EXPONENT
* Diffraction occurs for h+k+l is even
	END

	SUBROUTINE FCCS(PI,H,K,L,S)
* Define data types
	INTEGER H, K, L, HPLUSK, LPLUSK, HPLUSL
	REAL FACTHK, FACTLK, FACTHL
	COMPLEX EXPHK, EXPLK, EXPHL, S
* Sum h + k and find the exponential of pi*i*(h+k)
	HPLUSK = H + K
	FACTHK = PI * HPLUSK
	XHK = COS(FACTHK)
	YHK = SIN(FACTHK)
	EXPHK = CMPLX(XHK, YHK)
* Sum l + k and find the exponential of pi*i*(l+k)
	LPLUSK = L + K
	FACTLK = PI * LPLUSK
	XLK = COS(FACTLK)
	YLK = SIN(FACTLK)
	EXPLK = CMPLX(XLK, YLK)
* Sum h + l and find the exponential of pi*i*(h+l)
	HPLUSL = H + L
	FACTHL = PI * HPLUSL
	XHL = COS(FACTHL)
	YHL = SIN(FACTHL)
	EXPHL = CMPLX(XHL, YHL)
* Find the structure factor by adding the exponentials to 1
	S = (1,0) + EXPHK + EXPLK + EXPHL
* Diffraction occurs for hkl are all even or all even
	END

	SUBROUTINE SCS(PI,H,K,L,S)
* Define data types
	COMPLEX S
* Simple cubic structures allow diffraction from all h,k,l
	S = (1,0)
	END

	SUBROUTINE HCPS(PI,H,K,L,S)
* Define data types
	INTEGER H, K, L
	REAL PI, MILLER, FACTOR, X, Y
	COMPLEX EXPONENT, S
* Calculate relavent values
	MILLER = (2./3.)*H + (4./3.)*K + L
	FACTOR = PI * MILLER
	X = COS(FACTOR)
	Y = SIN(FACTOR)
	EXPONENT = CMPLX(X, Y)
* Find structure factor
	S = (1,0) + EXPONENT
* Diffraction is forbidden when l is odd and h+2k is divisible by 3
	END