C  PROGRAM  STRUFAC     PROGRAMMED BY BOB SNYDER   5/1/67
C  PROGRAMMED FOR THE I.B.M. 1620
C  RUv UNDER KINGSTRAN THIS PROGRAM DOES ABOUT 50 REFLECTIONS PER MINUTE
C  STRUFAC IS A GENERAL PROGRAM FOR OBTAINING STRUCTOR FACTORS FROM INTE
C  VALUES FOR THE EQUI-INCLINATION WEISSENBERG TECHNIQUE.  OPTIONS ALLOW
C  SCALING OF LAYERS TO A COMMON SCALE AND THE CALCULATION OF WEIGHTS BA
C  ON A MODIFIED CRUCHSHANK SCHEME.
C***********************************************************************
C  INPUT  SET  UP
C***********************************************************************
C  CARD 1  DIRECT CELL CONSTANTS (6F10.6)
C        COL 1-10  A IN ANGSTROMS
C        COL11-20  B *    *
C        COL21-30  C *    *
C        COL31-40  COS ALPHA
C        COL41-50  COS BETA
C        COL51-60  COS GAMMA
C  CARD 2   CONTROL CARD  (8I1)
C        COL 1  I1  =1 IF A IS OSCILLATION AXIS    =0 IFNOT
C        COL 2  I2  =1 IF B IS OSCILLATION AXIS    =0 IFNOT
C        COL 3  I3  =1 IF C IS OSCILLATION AXIS    =0 IFNOT
C        COL 4  LP  =0 IF LP CORRECTIONS ARE TO BE APPLIED    =1 IF NOT
C        COL 5  ISPOT  =0 IF SPOTSHAPE CORRECTION IS TO BE APPLIED  =1 I
C        COL 6  IWT  =0 IF WEIGHTS ARE TO BE CALCULATED    =1 IF NOT
C        COL 7  =0  PROGRAM GENERATES BOTH PUNCHED AND PRINTED OUTPUT
C               =1 PRINTED OUTPUT ONLY
C        COL 8  ISCALE  =1 SCALE FACTORS WILL BE APPLIED
C                       =0  ALL SCALE FACTORS WILL BE SET TO 1.0
C***********************************************************************
C  IF ISCALE =1 THE NEXT TWO CARDS SHOULD BE TREATED AS BATCH CARDS AND
C  PRECEED EACH LAYER OF DATA FOR WHICH THE SCALE FACTOR DIFFERS
C  IF ISCALE=0  THEY NEED ONLY BE PUT IN ONCE
C***********************************************************************
C  CARD 3  BATCH CARD 1  TITLE  (80A1)
C        COL 1-80  ANY ALPHAMERIC INFORMATION
C  CARD 4  BATCH CARD 2  (4F10.6)
C        COL 1-10  WAVELENGTH OF RADIATION USED IN ANGSTROMS
C        COL11-20  C1 =THE ERROR OF INTENSE REFLECTIONS IN PERCENT(USUAL
C        COL21-30  C2 =THE ERROR IN READING WEAK REFLECTIONS IN RELATIVE
C                         INTENSITY UNITS (USUALLY 2. OR 3.)
C        COL31-40  SF  = THE SCALE FACTOR BY WHICH ALL REFLECTIONS IN TH
C                             WILL BE MULTIPLIED
C  NOTE* IFIWT=1  C1 AND C2  MAY BE LEFT BLANK*
C      * IF ISCALE=0 SF MAY BE LEFT BLANK*
C  CARD 5 (AND ON)  DATA CARDS  (3I3,F10.0)
C        COL1-3  H
C        COL4-6  K
C        COL7-9  L
C        COL10-19  INTENSITY
C    FOLLOW THE DATA WITH A CARD CONTAINING A PUNCH ONLY IN COL 72  NAME
C  LAST CARD   (71XI1)
C        COL 72  =1 MEANS READ TWO NEW BATCH CARDS AND A NEW SET OF REFL
C        COL 72  =2  MEANS END OF DATA CALL EXIT
C***********************************************************************
C***********************************************************************
	DIMENSION T(16)
990	PRINT32
32	FORMAT(1H1)
	READ3,A,B,C,CA,CB,CG
3	FORMAT(6F10.6)
	PRINT4,A,B,C,CA,CB,CG
4	FORMAT(15H0CELL CONSTANTS//3H A=F10.6,3H B=F10.6,3H C=F10.6,
	111H COS ALPHA=F10.6, 10H COS BETA=F10.6, 11H COS GAMMA=F10.6)
	READ5,I1,I2,I3,LP,ISPOT,IWT,IPCH,ISCALE
5	FORMAT(8I1)
	IF(I1)7,7,6
6	ROT=A
	PRINT8
8	FORMAT(35H0THE OSCILLATION AXIS IS THE A AXIS)
	GO TO 13
7	IF(I2)10,10,9
9	ROT=B
	PRINT11
11	FORMAT(35H0THE OSCILLATION AXIS IS THE B AXIS)
	GO TO 13
10	ROT=C
	PRINT12
12	FORMAT(35H0THE OSCILLATION AXIS IS THE C AXIS)
13	IF(LP)15,14,15
14	PRINT16
16	FORMAT(31H0LP CORRECTIONS WILL BE APPLIED)
15	IF(ISPOT)18,17,18
17	PRINT19
19	FORMAT(72H0EQUI-INCLINATION WEISSENBERG SPOT EXTENSION CORRECTIONS
	1 WILL BE APPLIED)
18	IF(IWT)21,20,21
20	PRINT22
22	FORMAT(27H0WEIGHTS WILL BE CALCULATED)
21	IF(IPCH)63,23,63
23	PRINT25
25	FORMAT(22H0CARDS WILL BE PUNCHED)
63	IF(ISCALE)997,997,61
61	PRINT62
62	FORMAT(38H0SCALE FACTORS WILL BE APPLIED TO DATA)
997	READ1,(T(J),J=1,16)
1	FORMAT(16A4)
	PRINT2,(T(J),J=1,16)
2	FORMAT(1H116A4)
24	READ26,ELAMDA,C1,C2,SF
26	FORMAT(4F10.6)
	PRINT27,ELAMDA
27	FORMAT(16H0THE WAVELENGTH=F10.6)
	IF(IWT)30,28,30
28	PRINT29,C1
29	FORMAT(39H0THE ERROR IN INTENSE REFLECTIONS(C1) =F5.2,8H PERCENT)
	PRINT31,C2
31	FORMAT(46H0THE ERROR IN READING WEAK REFLECTIONS(C2) IS ,F5.2,
	125H RELATIVE INTENSITY UNITS)
	C1=C1/100.
30	IF(ISCALE)64,64,65
65	PRINT66,SF
66	FORMAT(51H0THE SCALE FACTOR FOR THIS BATCH OF REFLECTIONS IS ,
	1F10.6)
	GO TO 67
64	SF=1.0
67	NUM=0
	NTAG=14
	PRINT32
60	READ33,J,K,L,AI,M
33	FORMAT(3I3,F10.0,52XI1)
	IF(M)998,34,998
34	NUM=NUM+1
	NTAG=NTAG+1
	IF(AI)101,100,101
100	AI=-C2
101	CONTINUE
C THE FOLLOWING CALCULATES SIN THETA AND SIN2*THETA
	AJ=J
	AK=K
	AL=L
C  SA# SIN OF ALPHA
	SA=SQRTF(1.0-CA**2)
	SB=SQRTF(1.0-CB**2)
	SG=SQRTF(1.0-CG**2)
C  V# THE VOLUME OF THE UNIT CELL
	V=A*B*C*SQRTF(1.0+2.0*CA*CB*CG-CA**2-CB**2-CG**2)
C  AS# A STAR
	AS=B*C*SA/V
	BS=A*C*SB/V
	CS=A*B*SG/V
C  CAS= COS OF ALPHA STAR
	CAS=(CB*CG-CA)/(SB*SG)
	CBS=(CG*CA-CB)/(SG*SA)
	CGS=(CA*CB-CG)/(SA*SB)
C  STOL=SIN THETA ON LAMDA   ST=SIN THETA   CT=  COS THETA
	STOL=0.5*SQRTF((AJ*AS)**2+(AK*BS)**2+(AL*CS)**2+2.0*(AJ*AS)*(AK*BS
	1)*CGS+2.0*(AJ*AS)*AL*CS *CBS+2.0*(AK*BS)*(AL*CS)*CAS)
	ST=STOL*ELAMDA
	STTEST=ST**2
	IF(STTEST.LE.1.0)GO TO 84
	PRINT85,J,K,L,AI
85	FORMAT(4H0THE 3I3,F7.2,82H REFLECTION HAS SIN**2 THETA GREATER THA
	1N ONE AND HAS BEEN OMITTED FROM PROCESSING,/)
	IF(M-1)60,990,999
84	CT=SQRTF(1.0-STTEST)
C  S2T=  SIN TWO THETA      C2T=  COS TWO THETA
	S2T=2.0*ST*CT
	C2T=1.0-2.0*ST**2
C  THE FOLLOWING CALCULATES THE LORENTZ POLORIZATION CORRECTION
	IF(LP)36,35,36
35	IF(I1)38,37,38
38	AROT=ABSF(J)
	GO TO 75
37	IF(I2)40,39,40
40	AROT=ABSF(K)
	GO TO 75
39	AROT=ABSF(L)
75	SINMU=AROT*ELAMDA/(2.0*ROT)
	 PL=(4.0*CT*SQRTF(ST**2-SINMU**2))/(1.0+(1.0-2.0*ST**2)**2)
	GO TO 70
36	PL=1.0
C  THE FOLLOWING CALCULATES THE SPOT EXTENSION CORRECTION AFTER
C  D.C. PHILLIPS ACTA CRYST. 7,746(1954)
70	IF(ISPOT)41,42,41
42	ZETA=AROT*ELAMDA/ROT
	XI=4.0*ST**2/ZETA**2
	PI=3.1415927
	SPOT=  180./(4.0*PI)*ZETA*SQRTF(SQRTF(4.0-ZETA**2)/(XI-1.0))
	SPOT=1.0+  SPOT/(28.7/SQRTF(1.0-ZETA**2/4.0)+  75.0)
	GO TO 71
41	SPOT=1.0
71	FSQ=SPOT*PL*SF*AI
	IF(FSQ)80,80,81
80	FSQ=ABS(FSQ)
	F=SQRT(FSQ)
	F=-F
	FSQ=-FSQ
	GO TO 82
81	F=SQRTF(FSQ)
82	IF(IWT)44,43,44
43	WTPR=4000. /(FSQ+C2**2*FSQ/(C1**2*AI**2))
	WTPU=FSQ*WTPR/1000.
	GO TO 72
44	WTPR=1.0
	WTPU=1.0
72	IF(NTAG-15)46,45,45
45	PRINT47
47	FORMAT( 55H0NUM  H  K  L    IOBS SPOTSHAPE  LP      F**2       F
	148H     WT.    WT*FSQ   SINT/L     SINT      **SF**)
55	FORMAT(1H0)
	NTAG=0
46	IF(IPCH)48,49,48
49	PUNCH50,J,K,L,FSQ,WTPU
50	FORMAT(3I9  ,4F9.2)
48	PRINT51,NUM,J,K,L,AI,SPOT,PL,FSQ,F,WTPR,WTPU,STOL,ST,SF
	WRITE(1,50)J,K,L,FSQ,WTPU,SF,AI
51	FORMAT(I5,3I3,F8.0,2F8.4,F10.2,3F9.3,2F10.6,F10.4)
998	IF(M-1)60,990,999
999	CALL EXIT
	END
