C     A THREE DIMENSIONAL FOURIER PROGRAM
C                         FOR IBM 1130 WITH 16K CORE STORAGE
C                    1442  CARD READ/PUNCH, NO ON-LINE PRINTER
C              PUNCH OUT THE RESULT AT THE END, ALL SCALED.
C                  JUNE 1968 AND JANUARY 1969
C                                                R.SHIONO
C  REVISD FOR XDS SIGMA 5 BY R. SNYDER 8/1/72
	DIMENSION ITEE(32),ITEO(32),ITOE(32),ITOO(32),IEXP(32)
	DIMENSION  FMULT(6)
	DIMENSION COSK(1581),SINK(1581)
	DIMENSION  IRHO(150)
	DIMENSION  GRS(3),GRE(3),GRI(3),ITP1(3),ITP2(3),ITP3(3),ITP4(3)
	DIMENSION LDR(3)
	DIMENSION  RHO(150)
	DATA LDR/'X','Y','Z'/
100	FORMAT (  5I4, 6F5.3  )
101	FORMAT ( 9F4.2   )
102	FORMAT (F7.5,F7.2,E10.3,2I2)
103	FORMAT (32I2)
104	FORMAT ( / 9HHORIZ.   ,A2,2X3F6.3/9HVERT.    ,A2,2X3F6.3 /
	1  9HSECTION     A2,2X 3F6.3   )
105	FORMAT (3I9,F9.4,23X,I1)
106	FORMAT(/10H(KMAX*X1)=I6,16H IS LARGER THAN I5,
	112H, CALL EXIT/)
107	FORMAT (/39HTOO MANY SECTIONS SPECIFIED FOR ONE RUN / )
108	FORMAT (1H )
109	FORMAT (/1H ,A2,2XA2,1H=,18I4 )
110	FORMAT (A2,2H =F6.3 )
111	FORMAT (F5.3,I2,1X18I4 )
112	FORMAT (//1H ,9HRHO(MAX)= E10.3,10H,  FACTOR=F11.5 //)
113	FORMAT  ( // 3HEND  //   )
114	FORMAT ( 3I4, 6F7.2, F5.3, I3  )
115	FORMAT (/17HNO. OF REFLEXIONS I6,10H, H2(MAX)=I3/ )
C     DEFINE INPUT OUTPUT DEVICE NUMBERS
	JTYP=108
	JCDI=105
	JCDO=108
C     THE DIMENSION OF THE INTERMEDIATE STORAGE
	JDIM=1581
	JRHO=150
C     THE MAXIMUM NUMBER OF VALUES TO BE PUNCHED ON A CARD, WILL BE
C             DIFFERENT FOR A PRINTER
	JJMX=18
C     READ CONTROL CARDS
208	READ (JCDI,100) ITYPE,ISORT,KMAX,ITES1,ITES2,FMULT
	READ ( JCDI,101 ) (GRS(I),GRE(I),GRI(I),I=1,3  )
	ITES1=ITES1+1
	ITES2=ITES2+1
	READ (JCDI,102 )SCALI,CONST,DMAX,LSW0,LSW1
	IF  ( SCALI )  246,  245,  246
245	SCALI=1.0
246	READ (JCDI,103 ) ITEE,ITEO,ITOE,ITOO
	GO TO ( 1, 2, 3, 4, 5, 6 ), ISORT
1	I=1
	J=2
	K=3
	GO TO 7
2	I=1
	J=3
	K=2
	GO TO 7
3	I=2
	J=1
	K=3
	GO TO 7
4	I=2
	J=3
	K=1
	GO TO 7
5	I=3
	J=1
	K=2
	GO TO 7
6	I=3
	J=2
	K=1
7	XMIN=GRS(I)
	XMAX=GRE(I)
	XINC=GRI(I)
	YMIN=GRS(J)
	YMAX=GRE(J)
	YINC=GRI(J)
	ZMIN=GRS(K)
	ZMAX=GRE(K)
	ZINC=GRI(K)
	JSEC=K
	JVER=J
	JHOR=I
	KMAXC=0
	INCX=XINC*1000.0+0.5
	MINX=XMIN*1000.0+0.5
	NH1=(XMAX-XMIN)/XINC+1.3
	IS=(YMAX-YMIN)/YINC+1.3
	IT=(ZMAX-ZMIN)/ZINC+1.3
	IR=NH1
	NIR=1
227	IF ( IR-JJMX*NIR )  229,  229,  228
228	NIR=NIR+1
	GO TO 227
229	Z=ZMIN
	WRITE ( JTYP,104 ) LDR(I),XMIN,XMAX,XINC,LDR(J),YMIN,YMAX,YINC,
	1   LDR(K),ZMIN,ZMAX,ZINC
C     REARRANGE THE TRIG. COMBINATIONS
	DO  210  N1=1,29,4
	DO  211  N2=1,3
	N3=N1+N2
	ITP1(N2)=ITEE(N3)
	ITP2(N2)=ITEO(N3)
	ITP3(N2)=ITOE(N3)
211	ITP4(N2)=ITOO(N3)
	ITEE(N1+1)=ITP1(I)
	ITEE(N1+2)=ITP1(J)
	ITEE(N1+3)=ITP1(K)
	ITEO(N1+1)=ITP2(I)
	ITEO(N1+2)=ITP2(J)
	ITEO(N1+3)=ITP2(K)
	ITOE(N1+1)=ITP3(I)
	ITOE(N1+2)=ITP3(J)
	ITOE(N1+3)=ITP3(K)
	ITOO(N1+1)=ITP4(I)
	ITOO(N1+2)=ITP4(J)
210	ITOO(N1+3)=ITP4(K)
	DO  226  I=1,32
	ITEE(I)=ITEE(I)-1
	ITEO(I)=ITEO(I)-1
	ITOE(I)=ITOE(I)-1
226	ITOO(I)=ITOO(I)-1
	LSW1 = LSW1 + 1
C     IF THE REFLEXION DATA ARE ON FILE
	GO TO ( 206, 205 ), LSW1
C     READ IN ALL DATA
205	NDAT = 0
	LSW0 = LSW0 + 1
221	GO TO  ( 220, 200 ), LSW0
C     FOURIER INPUT FORMAT
200	READ ( JCDI,114 ) IH,IK,IL,FOBS,FCAL,AOBS,ACAL,BOBS,BCAL,SINT
222	IF  ( IH )  201,  202,  201
202	IF  ( IK )  201,  203,  201
203	IF  ( IL )  201,  204,  201
C     SET UP ALL THE NECESSARY INFORMATION FOR SUMMATION AND
C             WRITE THE RESULTS ON DISK
C     SET UP THE AXIAL ARRANGEMENT
201	GO TO ( 11, 12, 13, 14, 15, 16), ISORT
11	H1=IH
	H2=IK
	H3=IL
	GO TO 17
12	H1=IH
	H2=IL
	H3=IK
	GO TO 17
13	H1=IK
	H2=IH
	H3=IL
	GO TO 17
14	H1=IK
	H2=IL
	H3=IH
	GO TO 17
15	H1=IL
	H2=IH
	H3=IK
	GO TO 17
16	H1=IL
	H2=IK
	H3=IH
17	KI=ABS(H2)
	KINH1=KI*NH1
C     PREPARE COEFFICIENTS
	GO TO  ( 19, 22, 23, 24 ), ITYPE
19	COEFA=AOBS
	COEFB=BOBS
	GO TO 26
22	COEFA=ACAL
	COEFB=BCAL
	GO TO 26
23	COEFA=AOBS-ACAL
	COEFB=BOBS-BCAL
	GO TO 26
24	COEFA=FOBS*FOBS
	COEFB=0.0
26	M=1
C     TEST INDICES FOR TRIG. COMBINATIONS
	GO TO ( 38, 28, 29, 30, 31, 32, 33, 34), ITES1
28	NO=IH
	GO TO (35,40,43),M
29	NO=IK
	GO TO (35,40,43),M
30	NO=IL
	GO TO (35,40,43),M
31	NO=IH+IK
	GO TO (35,40,43),M
32	NO=IH+IL
	GO TO (35,40,43),M
33	NO=IK+IL
	GO TO (35,40,43),M
34	NO=IH+IK+IL
	GO TO (35,40,43),M
35	IF ( NO-2*(NO/2) )  39,  38,  39
38	M=2
	GO TO ( 230, 28, 29, 30, 31, 32, 33, 34), ITES2
39	M=3
	GO TO ( 231, 28, 29, 30, 31, 32, 33, 34), ITES2
40	IF  ( NO-2*(NO/2) )  231,  230,  231
43	IF  ( NO-2*(NO/2) )  233,  232,  233
230	M=1
	GO TO 54
231	M=2
	GO TO 54
232	M=3
	GO TO 54
233	M=4
C     DETERMINE THE MULTIPLICITY
54	IF ( IH )  57,  55,  57
55	IF  ( IK )  60,  56,  60
56	FMUL=FMULT(6)
	GO TO 67
57	IF ( IL )  64,  58,  64
58	IF  ( IK )  63,  59,  63
59	FMUL=FMULT(4)
	GO TO 67
60	IF  ( IL )  62,  61,  62
61	FMUL=FMULT(5)
	GO TO 67
62	FMUL=FMULT(1)
	GO TO 67
63	FMUL=FMULT(3)
	GO TO 67
64	IF  ( IK )  66,  65,  66
65	FMUL=FMULT(2)
	GO TO 67
66	FMUL = 1.0
67	COEFA=COEFA*FMUL*SCALI
	COEFB=COEFB*FMUL*SCALI
	WRITE (1)      H1,H2,H3,COEFA,COEFB,M,KINH1
	NDAT = NDAT + 1
	GO TO ( 215, 216, 217 ), JVER
215	KCRT=IABS(IH)
	GO TO 218
216	KCRT=IABS(IK)
	GO TO 218
217	KCRT=IABS(IL)
218	IF  ( KMAXC-KCRT) 219,  221,  221
219	KMAXC=KCRT
	GO TO 221
C     STRUCTURE FACTOR INPUT FORMAT
C     ORFLS INPUT FOR PATTERSON ONLY
220	READ(JCDI,105) IH,IK,IL,FOBS,IUN
	IF(IUN.EQ.1) GO TO 220
	GO TO 222
C     END OF DATA
C     LAST RECORD INDICATION
204	H1=999.0
	WRITE ( 1      ) H1,H2,H3,COEFA,COEFB,M,KINH1
	KMAX=KMAXC
	NDAT=NDAT-1
	WRITE (JTYP,115 ) NDAT,KMAX
	IF ( NDAT )  252,  252,  206
206	KMAX=KMAX+1
	ITC=0
	NFSM=1
	JMAX=KMAX*NH1
	IF  ( JMAX-JDIM )  209,  209,  207
209	IF  ( IR-JRHO )  250,  250,  207
250	GO TO 97
207	WRITE(JTYP,106) JMAX,JDIM; CALL EXIT
252	CONTINUE
	GO TO 208
C
C     UP TO THIS PART MAY BE SET INTO A SEPARATE LINK,IF NECESSARY
C     ACTUAL CALCULATION STARTS FROM HERE
97	REWIND 1
C     CLEAR THE SUM AREA
	DO  8  I=1,JMAX
	COSK(I)=0.0
8	SINK(I)=0.0
	ITC=ITC+1
	ZPR=Z+0.0005
	WRITE (JTYP,110) LDR(JSEC),ZPR
C     START THE FIRST AND THE SECOND SUMMATION
9	READ (1      ) H1,H2,H3,COEFA,COEFB,M,KINH1
	IF  ( H1-999.0 )  10,  86,  86
10	COEF=COEFA
	THET=6.28319*H3*Z
	COSLZ=COS(THET)
	SINLZ=SIN(THET)
	GO TO ( 46, 48, 50, 52), M
46	DO 47 I=1,32
47	IEXP(I)=ITEE(I)
	GO TO 248
48	DO 49 I=1,32
49	IEXP(I)=ITEO(I)
	GO TO 248
50	DO 51 I=1,32
51	IEXP(I)=ITOE(I)
	GO TO 248
52	DO 53 I=1,32
53	IEXP(I)=ITOO(I)
248	DTH=6.28319*H1*XINC
	THETA=6.28319*H1*XMIN-DTH
	DO  85  I=1,17,16
C     BTPASS IF THE AMPLITUDE IS ZERO
	IF  ( COEF )  69,  75,  69
69	N=I+12
	DO  84  J=I,N,4
	IF  (IEXP(J+3) )  75,  70,  71
70	SUMA=COEF*COSLZ
	GO TO 72
71	SUMA=COEF*SINLZ
72	IF  ( IEXP(J) )  87,  74,  73
73	SUMA=-SUMA
74	THET=THETA
	DO  83  M=1,IR
	THET=THET+DTH
	MN=KINH1+M
	IF  ( IEXP(J+1) )  87,  76,  77
76	SUM=SUMA*COS(THET)
	GO TO 78
77	SUM=SUMA*SIN(THET)
78	IF  ( IEXP(J+2) )  87,  79,  80
79	COSK(MN)=COSK(MN)+SUM
	GO TO 83
80	IF  ( H2 )  81,  82,  82
81	SUM=-SUM
82	SINK(MN)=SINK(MN)+SUM
83	CONTINUE
84	CONTINUE
75	IF  ( COEFB )  85,  9,  85
85	COEF=COEFB
	GO TO 9
C     END OF TWO SUMMATIONS
87	CONTINUE
86	Y=YMIN
	WRITE (2)      ZPR
	DO  91  J=1,IS
	DTH=6.28319*Y
	DO  93  M=1,IR
	SUM=0.0
	THET=-DTH
	MN=M-NH1
	DO  92  KI=1,KMAX
	THET=THET+DTH
	MN=MN+NH1
92	SUM=SUM+COSK(MN)*COS(THET)+SINK(MN)*SIN(THET)
	SUM=SUM+CONST
	IF  ( SUM-DMAX )  93,  93,  224
224	DMAX=SUM
93	RHO(M)=SUM
	YO=Y+0.0005
	WRITE (2)      YO,(RHO(I),I=1,IR)
91	Y=Y+YINC
	GO TO 96
96	Z=Z+ZINC
	GO TO 98
244	IF (Z-ZMAX )  97,  97,  98
C     END OF CALCULATION, START PUNCHING THE RESULTS
C     THE PART AFTER THIS MAY BE MADE INTO A SEPARATE LINK
98	REWIND 2
	DFAC=999.5/DMAX
	RDFAC=1.0/DFAC
	WRITE (JTYP,112) DMAX,RDFAC
C     SECTIONS UP TO ITC (LAST SECTION CALCULATED)
	DO  243 J=1,ITC
	READ (2)      Z
	WRITE ( JCDO,110 ) LDR(JSEC),Z
	NFSMC=NFSM
C     SET UP THE HEADING  FOR HORIZONTAL COORDINATES
	IRHO(1)=MINX
	IRS=IR-1
	DO  247  I=1,IRS
247	IRHO(I+1)=IRHO(I)+INCX
C     PUNCH JJMX VALUES PER CARD, GO TO NEXT LINE
C     HORIZONTAL DIRECTION, IR GRID POINTS
	JJ=1
	DO  212  M=1,NIR
	JM=JJ+JJMX-1
	IF  ( JM-IR )  214,  214,  213
213	JM=IR
214	NFSM=NFSMC
	WRITE (JCDO,109 ) LDR(JVER),LDR(JHOR),(IRHO(I),I=JJ,JM)
	WRITE (JCDO,108 )
C     VERTICAL DIRECTION, IS LINES
	DO  240  KI=1,IS
	READ (2)      YO,(RHO(I),I=1,IR)
C     NORMALIZE TO THE MAXIMUM VALUE OF 999
	DO  241  I=JJ,JM
	IRHO(I)=RHO(I)*DFAC
	IF  ( IRHO(I)+999 )  242,  241,  241
242	IRHO(I)=-999
241	CONTINUE
	WRITE (JCDO,111) YO,M,(IRHO(I),I=JJ,JM)
240	CONTINUE
	WRITE (JCDO,108)
212	JJ=JM+1
	WRITE (JCDO,108)
243	CONTINUE
	WRITE (JTYP,113)
	CALL EXIT
	END
