	DIMENSION WG(8),FNA(10)
	DIMENSION WGP(8),WGSN(10)
	DIMENSION TITLE(20),ATME(10)
	COMMON JCDI,JTYP,JCDO,NDAT,ITP1
	COMMON KN1,KN2,KN3,KN4,KN5,KN6
	COMMON SCALE,SCALA,TEMP,TEMP1,TEMP2,CONST
	COMMON FNA,WG,ATME,TITLE,NTOT
	COMMON WGP,WGSN
C     DEFINE INPUT/OUTPUT UNIT
	JCDI=5
	JCDO=7
	JTYP=6
	ITP1=1
10	CALL NORME
	IF(KN5) 10,10,20
20	CALL ESORT
	GO TO 20
30	CALL EXIT
	END

	SUBROUTINE NORME
C     PREPARATION OF NORMALIZED STRUCTURE FACTORS
C
C                                   MAY 1964        R.SHIONO
C                            JULY 1966
C               FOR IBM 1130 COMPUTER WITH 1442 READ/PUNCH
C                       APRIL AND DECEMBER 1968
C ADAPDET FOR THE XDS SIGMA 5 8/72 BY R. SNYDER
	DIMENSION WG(8),FM(10),FNA(10),SNORM(8),SNRMS(8),SENMB(8)
	DIMENSION  TAS(8),TAN(8),JNMB(10)
	DIMENSION  WGP(8),WGSN(10)
	DIMENSION  NMB(8)
	DIMENSION  SNORE(8)
	DIMENSION  NMC(10),SINE(10),SINES(10),SINEM(10)
	DIMENSION TITLE(20), ATME(10)
	DIMENSION  DAMP(3)
	COMMON  JCDI,JTYP,JCDO,NDAT, ITP1
	COMMON  KN1,KN2,KN3,KN4,KN5,KN6
	COMMON  SCALE,SCALA,TEMP,TEMP1,TEMP2,CONST
	COMMON  FNA,WG,ATME,TITLE,NTOT
	COMMON WGP,WGSN
C     READ TITLE AND CONTROL CARD
155	READ ( JCDI,10,END=190 ) (TITLE(I),I=1,20)
	READ (JCDI,12) KN1,KN2,KN3,KN4,KN5,KN6,DAMP(1),DAMP(2),DAMP(3)
	DO  194  I=1,3
	IF  (DAMP(I) )  193,  193,  194
193	DAMP(I)=0.4
194	CONTINUE
	IF  ( KN3-1 )  189,  189,  192
189	READ (JCDI,100) SCALE,TEMP,CONST,SCALA,TEMP1,TEMP2
192	IF ( KN6 )  187,  187,  188
187	READ (JCDI,101) (FNA(I),I=1,10)
	READ (JCDI,101) (WG(I),I=1,7)
	READ  ( JCDI,101 ) (WGP(I),I=1,8)
	READ (JCDI,101) (WGSN(I),I=1,10)
188	WRITE (JTYP,11 ) (TITLE(I),I=1,20)
	WRITE ( JTYP,205) (FNA(I),I=1,10)
	EPS2=0.0
	EPS3=0.0
	FSUM=0.0
	WG(8)=1.
	DO  113 K=1,8
	SNORM(K)=0.
	SNORE(K)=0.0
	SNRMS(K)=0.
	SENMB(K)=0.
	TAS(K)=0.
	TAN(K)=0.
	IF ( WGP(K) )  115,  115,  116
115	WGP(K)=1.0
116	IF ( WG(K) )  114,  114,  113
114	WG(K)=1.
113	CONTINUE
	WRITE ( JTYP,13) (WG(I),I=1,7)
	DO  117  K=1,10
	NMC(K)=0
	SINE(K)=0.0
	SINES(K)=0.0
	SINEM(K)=0.0
	IF  ( WGSN(K) )  112,  112,  117
112	WGSN(K)=1.0
117	CONTINUE
	WRITE ( JTYP,14 ) (WGP(I),I=1,8 )
	WRITE (JTYP,210)
	WRITE (JTYP,15 ) (WGSN(I),I=1,10)
	DO  111  K=1,10
111	JNMB(K)=0
C     IF DATA ARE ON DISK, BYPASS SCAT
	IF ( KN1-2 )  148,  149,  148
C     GO TO PICK UP ATOMIC SCATTERING FACTORS
148	CALL SCAT
	DO  152  K=1,10
152	FSUM=FSUM+ATME(K)*FNA(K)
	DO  153  K=1,10
	FFSUM=ATME(K)/FSUM
	EPS2=EPS2+FNA(K)*FFSUM*FFSUM
153	EPS3=EPS3+FNA(K)*FFSUM*FFSUM*FFSUM
	CONST=EPS3/SQRT(EPS2*EPS2*EPS2)
	WRITE ( JTYP,103) EPS2,EPS3,CONST
149	IF  ( KN3 )  150,  150,  151
150	CALL  WILSN
151	REWIND ITP1
	IF  ( KN2 )  140,  140,  157
157	WRITE (JTYP, 212 )
140	READ (ITP1)   IH,IK,IL,FOBS,FCAL,AO,AC,BO,BC,SIG,Q,ID,RHO,SINT,
	1    (FM(I),I=1,10)
	BO=0.0
	BC=0.0
	FH=IH
	FK=IK
	FL=IL
	IF  ( IH-999 )  145,  141,  145
145	IF  ( FH )  120,  121,  120
121	IF  (FK)  130,  122,  130
122	IF  (FL)  123,  141,  123
123	IREF=7
	GO TO 135
120	IF  (FK)  125,  124,  125
124	IF  (FL)  127,  126,  127
126	IREF=5
	GO TO 135
127	IREF=4
	GO TO 135
125	IF  (FL)  129,  128,  129
128	IREF=2
	GO TO 135
129	IREF=1
	GO TO 135
130	IF  (FL)  132,  131,  132
131	IREF=6
	GO TO 135
132	IREF=3
	GO TO 135
135	SUMF=0.
	SUMFU=0.0
	FIDH=IH/2
	FIDK=IK/2
	FIDL=IL/2
	IF  (.5*FH-FIDH)   161,   160,   161
160	IF  (.5*FK-FIDK)   165,  162,  165
162	IF   (.5*FL-FIDL)   164,   163,   164
163	KLAS=1
	GO TO 175
164	KLAS=2
	GO TO 175
165	IF  (.5*FL-FIDL)   167,  166,  167
166	KLAS=3
	GO TO 175
167	KLAS=4
	GO TO 175
161	IF  (.5*FK-FIDK)  169,  168,  169
168	IF  (.5*FL-FIDL)   171,  170,  171
170	KLAS=5
	GO TO 175
171	KLAS=6
	GO TO 175
169	IF  (.5*FL-FIDL)   173,  172,  173
172	KLAS=7
	GO TO 175
173	KLAS=8
175	DO  136  K=1,10
	SUMFU=SUMFU+FM(K)*FNA(K)
136	SUMF=SUMF+FM(K)*FM(K)*FNA(K)
	SINH=SINT
	DO  185  KSIN=1,10
	SINH=SINH-0.1
	IF  ( SINH )  186,  186,  185
185	CONTINUE
	KSIN=10
186	IF ( KN2 )  137,  137,  138
137	SFOE=SCALE*FOBS*EXP(TEMP*RHO)
	GO TO 139
138	SFOE=SCALA*FOBS*EXP((TEMP+TEMP1*RHO)*RHO)
139	ESF=SFOE*WGP(KLAS)*WGSN(KSIN)/SQRT(WG(IREF)*SUMF)
	UNIF=SFOE/SUMFU
	ESS=ESF*ESF
	ESM=ESS-1.
	IF  ( KN4 )  183,  183,  184
C     UNITARY STRUCTURE FACTOR REPLACES FOBS
184	FOBS=UNIF
183	SENMB(8)=SENMB(8)+1.0
	SNORM(8)=SNORM(8)+ESS
	SNRMS(8)=SNRMS(8)+ABS(ESM)
	SNORE(8)=SNORE(8)+ESF
	SNORM(IREF)=SNORM(IREF)+ESS
	SNORE(IREF)=SNORE(IREF)+ESF
	SNRMS(IREF)=SNRMS(IREF)+ABS(ESM)
	SENMB(IREF)=SENMB(IREF)+1.
	TAS(KLAS)=TAS(KLAS)+ESS
	TAN(KLAS)=TAN(KLAS)+1.
	NMC(KSIN)=NMC(KSIN)+1
	SINE(KSIN)=SINE(KSIN)+ESF
	SINES(KSIN)=SINES(KSIN)+ESS
	SINEM(KSIN)=SINEM(KSIN)+ABS(ESM)
	ESFA=ESF
	DO  180  K=1,10
	IF  (ESFA)   181,  181,  180
180	ESFA=ESFA-0.3
	K=10
181	JNMB(K)=JNMB(K)+1
	IF  ( KN5-1 )  140,  133,  133
133	BACKSPACE ITP1
	WRITE (ITP1)   IH,IK,IL,FOBS,ESF,KLAS,K,AC,BO,BC,SIG,Q,ID,RHO,
	1   SINT,(FM(I),I=1,10)
	GO TO 140
C      END OF REFLEXION DATA
141	NTOT=SENMB(8)
C     WRITE LAST RECOED
	IH=999
	WRITE (ITP1)   IH,IK,IL,FOBS,ESF,KLAS,K,AC,BO,BC,SIG,Q,ID,RHO,
	1  SINT,(FM(I),I=1,10)
	WRITE ( JTYP,102)
	DO  142  K=1,8
	IF  ( SENMB(K))  144,  143,  144
143	SNORM(K)=0.0
	SNRMS(K)=0.0
	GO TO 142
144	SNORM(K)=SNORM(K)/SENMB(K)
	SNRMS(K)=SNRMS(K)/SENMB(K)
	SNORE(K)=SNORE(K)/SENMB(K)
	WG(K)=WG(K)*(SNORM(K)-DAMP(1)*(SNORM(K)-1.0))
142	NMB(K)=SENMB(K)
	WRITE ( JTYP,106 ) (NMB(K),K=1,8 )
	WRITE ( JTYP,109 ) (SNORE(K),K=1,8)
	WRITE (JTYP,107 ) (SNORM(K),K=1,8)
	WRITE (JTYP,108 ) (SNRMS(K),K=1,8)
	DO  147  K=1,10
	IF  ( NMC(K) )  147,  147,  156
156	ANM=NMC(K)
	SINE(K)=SINE(K)/ANM
	SINES(K)=SINES(K)/ANM
	SINEM(K)=SINEM(K)/ANM
	WGSN(K)=WGSN(K)/SQRT(SINES(K)-DAMP(3)*(SINES(K)-1.0))
147	CONTINUE
	WRITE (JTYP,210 )
	WRITE (JTYP,211) (NMC(K),K=1,10)
	WRITE (JTYP,207) (SINE(K),K=1,10)
	WRITE (JTYP,208 ) (SINES(K),K=1,10)
	WRITE (JTYP,209 ) (SINEM(K),K=1,10)
	WRITE ( JTYP,201 )
	DO  182  K=1,8
	TAS(K)=TAS(K)/TAN(K)
	WGP(K)=WGP(K)/SQRT(TAS(K)-DAMP(2)*(TAS(K)-1.0))
182	NMB(K)=TAN(K)
	WRITE ( JTYP,202 ) (NMB(K),K=1,8 )
	WRITE ( JTYP,203 ) ( TAS(K),K=1,8 )
	WRITE ( JTYP,204 )
	WRITE ( JTYP,202 )  (JNMB(K),K=1,10)
190	CALL EXIT
191	RETURN
10	FORMAT (20A4 )
11	FORMAT (//1H ,20A4 )
12	FORMAT (6I2, 3F6.2 )
13	FORMAT (     20HWEIGHT OF REFLEXION, 7X3HHKL,7X3HHK0,7X3H0KL,
	1 7X3HH0L,7X3HH00,7X3H0K0,7X3H00L///20X 7F10.2    // )
14	FORMAT (27X,3HEEE,7X3HEEO,7X3HEOE,7X3HEOO,7X3HOEE,7X3HOEO,7X3HOOE,
	1  7X3HOOO,//20X8F10.2 //)
15	FORMAT (1H ,10X 10F10.2 //)
16	FORMAT (//1H ,43H*****  END OF E, PRESS START TO RERUN ***** //)
100	FORMAT (8F8.3 )
101	FORMAT (10F5.1 )
102	FORMAT(/ 1H ,18X,3HHKL,7X,3HHK0,7X,3H0KL,7X,3HH0L,7X,3HH00,7X,3H0K
	10,7X,3H00L,7X,3HALL,4X,9H(CENTRIC),2X10H(NON-CENT)     // )
103	FORMAT (/1H , 5HEPS2=F8.5,8H,  EPS3=,F9.6,18H,  E(3)/E(2)**3/2= ,
	1    F8.5  /// )
106	FORMAT ( 1H , 11HNO. REFL.   , 8I10   // 1H     )
107	FORMAT ( 1H ,11HAVRG. E**2 ,8F10.3,6X5H1.000,7X5H1.000  //1H  )
108	FORMAT ( 1H ,11HAVRG. ESQ-1 ,8F10.3,6X5H0.968,7X5H0.736  ///1H  )
109	FORMAT (1H ,11HAVRG. E    ,8F10.3,6X5H0.798,7X5H0.886  //1H   )
201	FORMAT(/ 1H ,41HNUMBER AND AVERAGE OF E**2 FOR EACH CLASS //18X,
	13HEEE,7X3HEEO,7X3HEOE,7X3HEOO,7X3HOEE,7X3HOEO,7X3HOOE,7X3HOOO// )
202	FORMAT ( 1H , 10X, 10I10  )
203	FORMAT (/1H , 10X8F10.3 // 1H  )
204	FORMAT (1H ,45HDISTRIBUTION ON MAGNITUDE OF E (.3 INTERVAL) ///19X
	194HE=0.   0.0-0.3   0.3-0.6   0.6-0.9   0.9-1.2   1.2-1.5   1.5-1.
	28   1.8-2.1   2.1-2.4   2.4-     // )
205	FORMAT (/1H ,41HNUMBER OF EACH KIND OF ATOMS IN UNIT CELL   ///
	1   10X, 10F8.0  ///)
207	FORMAT (1H ,11HAVRG. E    , 10F10.3 // )
208	FORMAT (1H ,11HAVRG. E**2 ,10F10.3 // )
209	FORMAT (1H ,11HAVRG. ESQ-1 ,10F10.3 // )
210	FORMAT (1H ,16HSINE THETA UP TO,1X2H.1,8X2H.2,8X2H.3,8X2H.4,
	1  8X2H.5,8X2H.6,8X2H.7,8X2H.8,8X2H.9,7X3H1.0 // )
211	FORMAT (1H ,11HNO. REFL.   ,10I10 // )
212	FORMAT (1H ,32HTEMP. FACTORS B1 AND B2 ARE USED //)
	END
C     SUBROUTINE TO READ THE ATOMIC SCATTERING FACTORS
C            AND INTERPOLATE AND WRITE A DATA TAPE
C                                           R.SHIONO
C                  FOR NORMALIZED SF PROGRAM
C                            JULY 1966
C              FOR IBM 1130 SYATEM                   APRIL 1968
	SUBROUTINE SCAT
C
	DIMENSION   RCP(6),PRCP(6),AFM(11,51),FM(15)
	DIMENSION  TITLE(20), ATME(10)
	DIMENSION  WG(8), FNA(10)
	COMMON JCDI,JTYP,JCDO,NDAT,ITP1
	COMMON  KN1,KN2,KN3,KN4,KN5,KN6
	COMMON  SCALE,SCALA,TEMP,TEMP1,TEMP2,CONST
	COMMON  FNA,WG,ATME,TITLE
1	FORMAT ( 3F7.3, 4F7.4    )
2	FORMAT   ( 7F10.4  )
3	FORMAT ( 1H , 30HATOMIC SCATTERING FACTORS F(0)   //11X,10F8.1// )
4	FORMAT ( 3I9,2F9.2,F9.0,I5,2F6.3 )
5	FORMAT( 1H ,3H A=,F8.3,4H, B=,F8.3,4H, C=,F8.3,7H, COSA=,F7.4,
	1   7H, COSB=,F7.4,7H, COSC=,F7.4,9H, LAMBDA=,F7.4  /   )
6	FORMAT ( 3I4,6F7.2,F5.3,I3,F7.2 )
	RA=0.0
	RB=0.0
C     READ CELL CONSTANT CARD
	READ ( JCDI,1) A,B,C,CSAA,CSAB,CSAC,WAVE
	WRITE ( JTYP,5 ) A,B,C,CSAA,CSAB,CSAC,WAVE
	D=SQRT(1.+2.*CSAA*CSAB*CSAC-CSAA*CSAA-CSAB*CSAB-CSAC*CSAC)
	D=1./(D*A*B*C)
	SNAA=SQRT(1.-CSAA*CSAA)
	SNAB=SQRT(1.-CSAB*CSAB)
	SNAC=SQRT(1.-CSAC*CSAC)
	RCP(1)=D*B*C*SNAA
	RCP(2)=D*C*A*SNAB
	RCP(3)=D*A*B*SNAC
	RCP(4)=(CSAB*CSAC-CSAA)/(SNAB*SNAC)
	RCP(5)=(CSAC*CSAA-CSAB)/(SNAC*SNAA)
	RCP(6)=(CSAA*CSAB-CSAC)/(SNAA*SNAB)
	PRCP(1)=RCP(1)*RCP(1)*.25
	PRCP(2)=RCP(2)*RCP(2)*.25
	PRCP(3)=RCP(3)*RCP(3)*.25
	PRCP(4)=RCP(1)*RCP(2)*RCP(6)*.25
	PRCP(5)=RCP(1)*RCP(3)*RCP(5)*.25
	PRCP(6)=RCP(2)*RCP(3)*RCP(4)*.25
C     WRITE THE RECIPROCAL AND THE DIRECT CELL CONSTANTS AS FIRST RECORD
	FCAL=0.
	AO=0.
	BO=0.
	AC=0.
	BC=0.
	DO  31  I=1,11
	DO  31  J=1,51
31	AFM(I,J)=0.
C     READ ATOMIC SCATTERING FACTOR TABLES
	I=0
20	I=I+1
	READ ( JCDI,2 ) (AFM(I,J),J=1,7 )
	IF  ( AFM(I,1) )  21, 21, 22
22	READ ( JCDI,2 ) (AFM(I,J),J=8,51 )
	NATSC=I
	GO TO 20
C     END OF F TABLE,START READING REFLEXION DATA
21	JBB=KN1+1
C     STORE F(0),S FOR LATER USE
	DO  32  I=1,10
32	ATME(I)=AFM(I,1)
	WRITE ( JTYP,3 ) (ATME(I),I=1,NATSC )
	NDAT=1
61	IF ( KN1 )  60,  60,  65
60	READ ( JCDI,4 ) IH,IK,IL,FOBS,SIG,Q,ID,RA,RB
62	FH=IH
	FK=IK
	FL=IL
	IF  (IH)  25,  23,  25
23	IF  (IK)  25,  24,  25
24	IF  (IL)  25,  50,  25
C     STRUCTURE FACTOR OUTPUT-FOURIER INPUT FORMAT
65	READ (JCDI,6 ) IH,IK,IL,FOBS,FCAL,AO,AC,BO,BC,SINT,ID,SIG
	GO TO 62
25	IF  ( Q )  64,  64,  63
64	Q=1.0
63	RHO=FH*FH*PRCP(1)+FK*FK*PRCP(2)+FL*FL*PRCP(3)+2.*FH*FK*PRCP(4)+
	1       2.*FK*FL*PRCP(6)+2.*FL*FH*PRCP(5)
	SINT=SQRT(RHO)*WAVE
	IF  ( SINT-1.0 )  26,  61,  61
26	BRAG=SINT
	DO  27  I=2,51
	BRAG=BRAG-0.02
	IF  (BRAG)  28,  28,  27
27	CONTINUE
28	I=I-1
	FRAC=(BRAG+0.02)*50.
	DO  29  K=1,10
29	FM(K)=0.
	DO  30  K=1,NATSC
30	FM(K)=AFM(K,I)+(AFM(K,I+1)-AFM(K,I))*FRAC
53	WRITE (ITP1)   IH,IK,IL,FOBS,FCAL,AO,AC,BO,BC,SIG,Q,ID,RHO,
	1   SINT,(FM(K),K=1,10)
	GO TO 61
50	IH=999
	WRITE (ITP1)   IH,IK,IL,FOBS,FCAL,AO,AC,BO,BC,SIG,Q,ID,RHO,
	1    SINT,(FM(K),K=1,10)
	REWIND ITP1
	RETURN
	END
C     WILSON-S PLOT FOR SCALE AND TEMPERATURE FACTOR
C             SINE THETA MIN. AND MAX. MAY BE SPECIFIED TO DISCRIMINATE
C             NUMBER OF SINE THETA REGIONS           1  TO  30
C        CALCULATES LINEAR AS WELL AS THS SECOND ORDER POLYNOMIAL FIT
C                            APRIL AND JUNE 1964           R.SHIONO
C                            JULY 1966
C
	SUBROUTINE WILSN
	DIMENSION  SOL(30),RI(30),RN(30),NA(30),FNA(10),FM(10)
	DIMENSION  SOLS(30)
	DIMENSION  TITLE(14),ATME(10)
	DIMENSION  WG(8)
	COMMON JCDI,JTYP,JCDO,NDAT, ITP1
	COMMON  KN1,KN2,KN3,KN4,KN5,KN6
	COMMON  SCALE,SCALA,TEMP,TEMP1,TEMP2,CONST
	COMMON   FNA,WG,ATME,TITLE
C     READ CONTROL CARD
	READ  (JCDI,90 ) NR,IFA,IFB,AMIN,AMAX,FACT
	IFA=IFA+1
	IFB=IFB+1
	IF  ( AMIN )  80,  80,  81
80	AMIN=0.
81	IF  (AMAX)   82,  82,  83
82	AMAX=1.0
83	FNR=NR
	DTH=(AMAX-AMIN)/FNR
	NAF=10
	DO  60  I=1,NR
	NA(I)=0
	SOL(I)=0.
	SOLS(I)=0.
	RI(I)=0.
60	RN(I)=0.
	NDAT=1
66	READ (ITP1)   IH,IK,IL,FOBS,FCAL,AO,AC,BO,BC,SIG,Q,ID,RHO,SINT,
	1      (FM(I),I=1,10)
	FH=IH
	FK=IK
	FL=IL
	IF  ( FH-999. )  64,  63,  63
64	IF  ( SINT )  66,  66,  67
67	GO TO ( 69, 61, 65 ), IFA
61	IF  ( FOBS )  69,  66,  69
65	FOBS=FOBS*FACT
69	AFSQ=0.0
	IF  (SINT-AMIN)  66,  84,  84
84	IF  (SINT-AMAX)  85,  85,  66
85	IF  ( FH )  50,  51,  50
51	IF  ( FK )  53,  52,  53
52	IF  ( FL )  54,  63,  54
54	IREF=7
	GO TO 44
50	IF  ( FK )  57,  56,  57
56	IF  ( FL )  59,  58,  59
58	IREF=5
	GO TO 44
59	IREF=4
	GO TO 44
57	IF  ( FL )  41,  40,  41
40	IREF=2
	GO TO 44
41	IREF=1
	GO TO 44
53	IF  ( FL )  43,  42,  43
42	IREF=6
	GO TO 44
43	IREF=3
44	GO TO ( 62, 78), IFB
62	DO  68  K=1,NAF
68	AFSQ=AFSQ+FM(K)*FM(K)*FNA(K)*WG(IREF)
	GO TO 79
78	AFSQ=FCAL*FCAL
79	THETA=DTH+AMIN
	DO  70  K=1,NR
	IF  (SINT-THETA)  71,  71,  70
70	THETA=THETA+DTH
	K=NR
71	SOL(K)=SOL(K)+RHO
	SOLS(K)=SOLS(K)+RHO*RHO
	RI(K)=RI(K)+AFSQ
	RN(K)=RN(K)+FOBS*FOBS
	NA(K)=NA(K)+1
	GO TO 66
C     END OF DATA
63	NZERO=0
	REWIND ITP1
	DO  73  K=1,NR
	IF  ( NA(K) )  76,  76,  77
76	SOL(K)=0.
	NZERO=NZERO+1
	GO TO 73
77	RI(K)=ALOG(RI(K)/RN(K))
	FL=NA(K)
	SOL(K)=SOL(K)/FL
	SOLS(K)=SOLS(K)/FL
73	CONTINUE
	WRITE (JTYP,91)
	WRITE ( JTYP,93)
	AGST=AMIN
	DO  75  K=1,NR
	AGFN=AGST+DTH
	WRITE ( JTYP,92) AGST,AGFN,NA(K),SOL(K),RI(K)
75	AGST=AGFN
	S0=NR-NZERO
	S1=0.
	S2=0.
	S3=0.
	S4=0.
	V0=0.
	V1=0.
	V2=0.
	DO  74  K=1,NR
	S1=S1+SOL(K)
	SLS=SOL(K)*SOL(K)
	S2=S2+SLS
	S3=S3+SOL(K)*SLS
	S4=S4+SLS*SLS
	V0=V0+RI(K)
	V2=V2+RI(K)*SLS
74	V1=V1+RI(K)*SOL(K)
C     CALCULATE A LINEAR FIT
	DD=S0*S2-S1*S1
	SCALE=(V0*S2-V1*S1)/DD
	TEMP=0.5*(S0*V1-S1*V0)/DD
	SCALE=EXP(0.5*SCALE)
	WRITE ( JTYP,94)  SCALE,TEMP
C     CALCULATE A SECOND ORDER POLYNOMIAL FIT
	DD=S4*DD+S2*(S1*S3-S2*S2)-S3*(S0*S3-S2*S1)
	SCALA=(V0*S2*S4+S2*V1*S3+V2*S1*S3-V0*S3*S3-V1*S1*S4-V2*S2*S2)/DD
	SCALA=EXP(0.5*SCALA)
	DD=DD+DD
	TEMP1=(S0*V1*S4+S2*S1*V2+S2*V0*S3-S0*S3*V2-S1*V0*S4-S2*V1*S2)/DD
	TEMPA=(S0*S2*V2+V0*S1*S3+S2*S1*V1-S0*V1*S3-S1*S1*V2-S2*S2*V0)/DD
	WRITE ( JTYP,96) SCALA,TEMP1,TEMPA
C     SET TEMPERATURE FACTOR TO ZERO IF NEGATIVE
	IF  ( TEMP )  45,  46,  46
45	TEMP=0.0
46	RETURN
90	FORMAT ( 3I5, 3F5.0  )
91	FORMAT (46H SCALE AND TEMPERATURE FACTOR BY WILSON METHOD ///   )
92	FORMAT ( 1H  F4.2,2H -F5.2, I15, F15.3,F15.3  )
93	FORMAT ( 1H 4X3HSIN9X44HNO. REFLEXION (SIN/LAMBDA)**2    ELOG(RATI
	1O)  //     )
94	FORMAT(//14H SCALE FACTOR= F 8.3, 16H,  TEMP. FACTOR= F 8.3//)
96	FORMAT(/14H SCALE FACTOR=F 8.3,18H, TEMP. FACTOR B1=F 8.3,5H, B2=
	1  F 8.3  //)
	END
C     E VALUE SORTING PROGRAM
C              TO BE CALLED AFTER THE E CALCULATION PROGRAM
C                  FOR IBM 1130 MODEL 2C WITH 16 K CORE.
C
C                                        JULY  1968
C                                         R.SHIONO
	SUBROUTINE ESORT
	DIMENSION IHK(5,1000),FOB(1000),ES(1000),SINE(1000)
	DIMENSION  FNA(10),WG(8),ATME(10),TITLE(20),WGP(8),WGSN(10)
	COMMON  ICDI,ITYP,ICDO,NDAT,ITAP
	COMMON KN1,KN2,KN3,KN4,KN5,KN6
	COMMON  SCALE,SCALA,TEMP,TEMP1,TEMP2,CONST
	COMMON  FNA,WG,ATME,TITLE,NTOT
	COMMON  WGP,WGSN
	REWIND ITAP
12	FORMAT ( 3I4,4F7.2,14X F5.3,I3,I3  )
	SCALB=SCALE
	IF  ( KN2 )  54,  54,  55
55	SCALB=SCALA
54	IF (KN4 )  56,  56,  57
57	SCALB=1.0
56	IF  ( KN5 )  51,  51,  31
31	GO TO ( 33,  32), KN5
C     SORT AND PUNCH E
C     STARTING GROUP NUMBER
32	NK=10
30	NREF=0
	REWIND ITAP
40	NREF=NREF+1
41	READ (ITAP)   (IHK(J,NREF),J=1,3),FOB(NREF),ES(NREF),IHK(4,NREF),
	1  K,AC,BO,BC,SIG,Q,IHK(5,NREF),RHO,SINE(NREF)
	IF  ( IHK(1,NREF)-999 )  42,  43,  43
42	IF  ( NK-K )  41,  40,  41
C     START SORTING ON THE MAGNITUDE OF E,S
43	NDATA=NREF-1
	IF  ( NDATA )  53,  53,  44
44	M1=NDATA
151	M1=M1/2
	IF  ( M1 )  152,  162,  152
152	M2=NDATA-M1
	M3=1
153	M4=M3
154	M5=M4+M1
	IF  ( ES(M4)-ES(M5) )  155,  157,  157
155	DO  156  IJ=1,5
	IHOLD=IHK(IJ,M4)
	IHK(IJ,M4)=IHK(IJ,M5)
156	IHK(IJ,M5)=IHOLD
	AHOLD=FOB(M4)
	FOB(M4)=FOB(M5)
	FOB(M5)=AHOLD
	AHOLD=ES(M4)
	ES(M4)=ES(M5)
	ES(M5)=AHOLD
	AHOLD=SINE(M4)
	SINE(M4)=SINE(M5)
	SINE(M5)=AHOLD
	M4=M4-M1
	IF  ( M4-1 )  157,  154,  154
157	M3=M3+1
	IF  ( M3-M2 )  153,  153,  151
C     END OF SORTING
162	DO  50  I=1,NDATA
	FOB(I)=FOB(I)*SCALB
	ESS=ES(I)*ES(I)
	ESM=ESS-1.0
50	WRITE (ICDO,12)(IHK(IJ,I),IJ=1,3),FOB(I),ES(I),ESS,ESM,SINE(I),
	1  IHK(5,I),IHK(4,I)
53	NK=NK-1
	IF  ( NK )  51,  51,  52
C     CONSOLE ENTRY SWITCH 0 ON TO TERMINATE PUNCHING AND RETURN
52	GO TO 30
C     PUNCH E CARDS WITHOU SORTING
33	REWIND ITAP
35	READ (ITAP,END=51) IH,IK,IL,FOBS,ESF,KLAS,K,AC,BO,BC,SIG,Q,ID,RHO,
	1SINT
	IF  ( IH-999 )  34,  51,  51
34	ESS=ESF*ESF
	ESM=ESS-1.0
	FOBS=FOBS*SCALB
	WRITE (ICDO,12) IH,IK,IL,FOBS,ESF,ESS,ESM,SINT,ID,KLAS
	GO TO 35
C     END OF ALL REFLEXIONS
51	RETURN
C     GO BACK AND REPEAT THE PROGRAM WITH A NEW PARAMETERS
	END
