/[MITgcm]/MITgcm/pkg/zonal_filt/fftpack.F
ViewVC logotype

Contents of /MITgcm/pkg/zonal_filt/fftpack.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph

Revision 1.5 - (show annotations) (download)
Wed Jun 6 15:47:54 2001 UTC (22 years, 11 months ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint46n_post, checkpoint47e_post, checkpoint44e_post, checkpoint46l_post, checkpoint46g_pre, checkpoint47c_post, release1_p13_pre, checkpoint50c_post, checkpoint46f_post, checkpoint48e_post, checkpoint50c_pre, checkpoint44f_post, checkpoint46b_post, checkpoint43a-release1mods, ecco_c50_e32, ecco_c50_e33, ecco_c50_e30, ecco_c50_e31, release1_p13, checkpoint40pre3, checkpoint40pre1, checkpoint40pre7, checkpoint40pre6, checkpoint48i_post, checkpoint40pre9, checkpoint40pre8, checkpoint46l_pre, chkpt44d_post, checkpoint51, checkpoint50, release1_p8, release1_p9, checkpoint50d_post, release1_p1, release1_p2, release1_p3, release1_p4, release1_p5, release1_p6, release1_p7, checkpoint50b_pre, checkpoint44e_pre, checkpoint51f_post, release1_b1, ecco_c51_e34d, ecco_c51_e34e, ecco_c51_e34f, ecco_c51_e34g, ecco_c51_e34a, ecco_c51_e34b, ecco_c51_e34c, checkpoint48b_post, checkpoint43, checkpoint51d_post, checkpoint48c_pre, checkpoint47d_pre, release1_chkpt44d_post, checkpoint47a_post, checkpoint48d_pre, checkpoint51j_post, checkpoint47i_post, checkpoint47d_post, icebear5, icebear4, icebear3, icebear2, checkpoint46d_pre, checkpoint40pre2, checkpoint48d_post, release1-branch_tutorials, checkpoint48f_post, checkpoint45d_post, checkpoint46j_pre, chkpt44a_post, checkpoint44h_pre, checkpoint48h_post, checkpoint40pre4, ecco_c50_e29, checkpoint51b_pre, checkpoint46a_post, checkpoint47g_post, checkpoint46j_post, checkpoint51h_pre, checkpoint46k_post, ecco_c50_e28, chkpt44c_pre, checkpoint48a_post, checkpoint45a_post, checkpoint50f_post, checkpoint50a_post, checkpoint50f_pre, ecco_c44_e19, ecco_c44_e18, ecco_c44_e17, ecco_c44_e16, release1_p12, release1_p10, release1_p11, release1_p16, release1_p17, release1_p14, release1_p15, checkpoint47j_post, ecco_c50_e33a, branch-exfmods-tag, checkpoint44g_post, branchpoint-genmake2, checkpoint46e_pre, checkpoint48c_post, checkpoint45b_post, checkpoint46b_pre, release1-branch-end, release1_final_v1, checkpoint51b_post, checkpoint51c_post, checkpoint46c_pre, checkpoint46, checkpoint47b_post, checkpoint44b_post, ecco_c51_e34, checkpoint46h_pre, checkpoint46m_post, checkpoint46a_pre, checkpoint50g_post, checkpoint45c_post, ecco_ice2, ecco_ice1, checkpoint44h_post, checkpoint46g_post, release1_p12_pre, ecco_c44_e22, checkpoint50h_post, checkpoint50e_pre, checkpoint50i_post, ecco_c44_e25, checkpoint51i_pre, checkpoint40pre5, checkpoint47f_post, checkpoint50e_post, chkpt44a_pre, checkpoint46i_post, ecco_c44_e23, ecco_c44_e20, ecco_c44_e21, ecco_c44_e26, ecco_c44_e27, ecco_c44_e24, checkpoint46c_post, ecco-branch-mod1, ecco-branch-mod2, ecco-branch-mod3, ecco-branch-mod4, ecco-branch-mod5, checkpoint50d_pre, checkpoint46e_post, release1_beta1, checkpoint51e_post, checkpoint44b_pre, checkpoint42, checkpoint40, checkpoint41, checkpoint47, checkpoint44, checkpoint45, checkpoint48, checkpoint49, checkpoint46h_post, checkpoint51f_pre, chkpt44c_post, checkpoint48g_post, checkpoint47h_post, checkpoint44f_pre, checkpoint51g_post, checkpoint46d_post, checkpoint50b_post, release1-branch_branchpoint, checkpoint51a_post
Branch point for: c24_e25_ice, branch-exfmods-curt, release1_final, release1-branch, branch-genmake2, release1, ecco-branch, release1_50yr, icebear, release1_coupled
Changes since 1.4: +50 -50 lines 
 o Added some missing "IMPLICIT REAL*8" (they were "IMPLICIT REAL"!).
 o Converted constants to D0 notation to avoid real*4 conversions.
1 C $Header: /u/gcmpack/models/MITgcmUV/pkg/zonal_filt/fftpack.F,v 1.4 2001/05/29 14:01:40 adcroft Exp $
2 C $Name: $
3 C SUBROUTINE RFFTB (N,R,WSAVE)
4 C DIMENSION R(1) ,WSAVE(1)
5 C IF (N .EQ. 1) RETURN
6 C CALL RFFTB1 (N,R,WSAVE,WSAVE(N+1),WSAVE(2*N+1))
7 C RETURN
8 C END
9 SUBROUTINE RADB2 (IDO,L1,CC,CH,WA1)
10
11 IMPLICIT REAL*8 (A-H,O-Z)
12 IMPLICIT INTEGER (I-N)
13
14 DIMENSION CC(IDO,2,L1) ,CH(IDO,L1,2) ,
15 1 WA1(1)
16 DO 101 K=1,L1
17 CH(1,K,1) = CC(1,1,K)+CC(IDO,2,K)
18 CH(1,K,2) = CC(1,1,K)-CC(IDO,2,K)
19 101 CONTINUE
20 IF (IDO-2) 107,105,102
21 102 IDP2 = IDO+2
22 DO 104 K=1,L1
23 DO 103 I=3,IDO,2
24 IC = IDP2-I
25 CH(I-1,K,1) = CC(I-1,1,K)+CC(IC-1,2,K)
26 TR2 = CC(I-1,1,K)-CC(IC-1,2,K)
27 CH(I,K,1) = CC(I,1,K)-CC(IC,2,K)
28 TI2 = CC(I,1,K)+CC(IC,2,K)
29 CH(I-1,K,2) = WA1(I-2)*TR2-WA1(I-1)*TI2
30 CH(I,K,2) = WA1(I-2)*TI2+WA1(I-1)*TR2
31 103 CONTINUE
32 104 CONTINUE
33 IF (MOD(IDO,2) .EQ. 1) RETURN
34 105 DO 106 K=1,L1
35 CH(IDO,K,1) = CC(IDO,1,K)+CC(IDO,1,K)
36 CH(IDO,K,2) = -(CC(1,2,K)+CC(1,2,K))
37 106 CONTINUE
38 107 RETURN
39 END
40 SUBROUTINE RADB3 (IDO,L1,CC,CH,WA1,WA2)
41
42 IMPLICIT REAL*8 (A-H,O-Z)
43 IMPLICIT INTEGER (I-N)
44
45 DIMENSION CC(IDO,3,L1) ,CH(IDO,L1,3) ,
46 1 WA1(1) ,WA2(1)
47 DATA TAUR,TAUI /-.5D0,.866025403784439D0/
48 DO 101 K=1,L1
49 TR2 = CC(IDO,2,K)+CC(IDO,2,K)
50 CR2 = CC(1,1,K)+TAUR*TR2
51 CH(1,K,1) = CC(1,1,K)+TR2
52 CI3 = TAUI*(CC(1,3,K)+CC(1,3,K))
53 CH(1,K,2) = CR2-CI3
54 CH(1,K,3) = CR2+CI3
55 101 CONTINUE
56 IF (IDO .EQ. 1) RETURN
57 IDP2 = IDO+2
58 DO 103 K=1,L1
59 DO 102 I=3,IDO,2
60 IC = IDP2-I
61 TR2 = CC(I-1,3,K)+CC(IC-1,2,K)
62 CR2 = CC(I-1,1,K)+TAUR*TR2
63 CH(I-1,K,1) = CC(I-1,1,K)+TR2
64 TI2 = CC(I,3,K)-CC(IC,2,K)
65 CI2 = CC(I,1,K)+TAUR*TI2
66 CH(I,K,1) = CC(I,1,K)+TI2
67 CR3 = TAUI*(CC(I-1,3,K)-CC(IC-1,2,K))
68 CI3 = TAUI*(CC(I,3,K)+CC(IC,2,K))
69 DR2 = CR2-CI3
70 DR3 = CR2+CI3
71 DI2 = CI2+CR3
72 DI3 = CI2-CR3
73 CH(I-1,K,2) = WA1(I-2)*DR2-WA1(I-1)*DI2
74 CH(I,K,2) = WA1(I-2)*DI2+WA1(I-1)*DR2
75 CH(I-1,K,3) = WA2(I-2)*DR3-WA2(I-1)*DI3
76 CH(I,K,3) = WA2(I-2)*DI3+WA2(I-1)*DR3
77 102 CONTINUE
78 103 CONTINUE
79 RETURN
80 END
81 SUBROUTINE RADB4 (IDO,L1,CC,CH,WA1,WA2,WA3)
82
83 IMPLICIT REAL*8 (A-H,O-Z)
84 IMPLICIT INTEGER (I-N)
85
86 DIMENSION CC(IDO,4,L1) ,CH(IDO,L1,4) ,
87 1 WA1(1) ,WA2(1) ,WA3(1)
88 DATA SQRT2 /1.414213562373095D0/
89 DO 101 K=1,L1
90 TR1 = CC(1,1,K)-CC(IDO,4,K)
91 TR2 = CC(1,1,K)+CC(IDO,4,K)
92 TR3 = CC(IDO,2,K)+CC(IDO,2,K)
93 TR4 = CC(1,3,K)+CC(1,3,K)
94 CH(1,K,1) = TR2+TR3
95 CH(1,K,2) = TR1-TR4
96 CH(1,K,3) = TR2-TR3
97 CH(1,K,4) = TR1+TR4
98 101 CONTINUE
99 IF (IDO-2) 107,105,102
100 102 IDP2 = IDO+2
101 DO 104 K=1,L1
102 DO 103 I=3,IDO,2
103 IC = IDP2-I
104 TI1 = CC(I,1,K)+CC(IC,4,K)
105 TI2 = CC(I,1,K)-CC(IC,4,K)
106 TI3 = CC(I,3,K)-CC(IC,2,K)
107 TR4 = CC(I,3,K)+CC(IC,2,K)
108 TR1 = CC(I-1,1,K)-CC(IC-1,4,K)
109 TR2 = CC(I-1,1,K)+CC(IC-1,4,K)
110 TI4 = CC(I-1,3,K)-CC(IC-1,2,K)
111 TR3 = CC(I-1,3,K)+CC(IC-1,2,K)
112 CH(I-1,K,1) = TR2+TR3
113 CR3 = TR2-TR3
114 CH(I,K,1) = TI2+TI3
115 CI3 = TI2-TI3
116 CR2 = TR1-TR4
117 CR4 = TR1+TR4
118 CI2 = TI1+TI4
119 CI4 = TI1-TI4
120 CH(I-1,K,2) = WA1(I-2)*CR2-WA1(I-1)*CI2
121 CH(I,K,2) = WA1(I-2)*CI2+WA1(I-1)*CR2
122 CH(I-1,K,3) = WA2(I-2)*CR3-WA2(I-1)*CI3
123 CH(I,K,3) = WA2(I-2)*CI3+WA2(I-1)*CR3
124 CH(I-1,K,4) = WA3(I-2)*CR4-WA3(I-1)*CI4
125 CH(I,K,4) = WA3(I-2)*CI4+WA3(I-1)*CR4
126 103 CONTINUE
127 104 CONTINUE
128 IF (MOD(IDO,2) .EQ. 1) RETURN
129 105 CONTINUE
130 DO 106 K=1,L1
131 TI1 = CC(1,2,K)+CC(1,4,K)
132 TI2 = CC(1,4,K)-CC(1,2,K)
133 TR1 = CC(IDO,1,K)-CC(IDO,3,K)
134 TR2 = CC(IDO,1,K)+CC(IDO,3,K)
135 CH(IDO,K,1) = TR2+TR2
136 CH(IDO,K,2) = SQRT2*(TR1-TI1)
137 CH(IDO,K,3) = TI2+TI2
138 CH(IDO,K,4) = -SQRT2*(TR1+TI1)
139 106 CONTINUE
140 107 RETURN
141 END
142 SUBROUTINE RADB5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
143
144 IMPLICIT REAL*8 (A-H,O-Z)
145 IMPLICIT INTEGER (I-N)
146
147 DIMENSION CC(IDO,5,L1) ,CH(IDO,L1,5) ,
148 1 WA1(1) ,WA2(1) ,WA3(1) ,WA4(1)
149 DATA TR11,TI11,TR12,TI12 /.309016994374947D0,.951056516295154D0,
150 1-.809016994374947D0,.587785252292473D0/
151 DO 101 K=1,L1
152 TI5 = CC(1,3,K)+CC(1,3,K)
153 TI4 = CC(1,5,K)+CC(1,5,K)
154 TR2 = CC(IDO,2,K)+CC(IDO,2,K)
155 TR3 = CC(IDO,4,K)+CC(IDO,4,K)
156 CH(1,K,1) = CC(1,1,K)+TR2+TR3
157 CR2 = CC(1,1,K)+TR11*TR2+TR12*TR3
158 CR3 = CC(1,1,K)+TR12*TR2+TR11*TR3
159 CI5 = TI11*TI5+TI12*TI4
160 CI4 = TI12*TI5-TI11*TI4
161 CH(1,K,2) = CR2-CI5
162 CH(1,K,3) = CR3-CI4
163 CH(1,K,4) = CR3+CI4
164 CH(1,K,5) = CR2+CI5
165 101 CONTINUE
166 IF (IDO .EQ. 1) RETURN
167 IDP2 = IDO+2
168 DO 103 K=1,L1
169 DO 102 I=3,IDO,2
170 IC = IDP2-I
171 TI5 = CC(I,3,K)+CC(IC,2,K)
172 TI2 = CC(I,3,K)-CC(IC,2,K)
173 TI4 = CC(I,5,K)+CC(IC,4,K)
174 TI3 = CC(I,5,K)-CC(IC,4,K)
175 TR5 = CC(I-1,3,K)-CC(IC-1,2,K)
176 TR2 = CC(I-1,3,K)+CC(IC-1,2,K)
177 TR4 = CC(I-1,5,K)-CC(IC-1,4,K)
178 TR3 = CC(I-1,5,K)+CC(IC-1,4,K)
179 CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3
180 CH(I,K,1) = CC(I,1,K)+TI2+TI3
181 CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3
182 CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3
183 CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3
184 CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3
185 CR5 = TI11*TR5+TI12*TR4
186 CI5 = TI11*TI5+TI12*TI4
187 CR4 = TI12*TR5-TI11*TR4
188 CI4 = TI12*TI5-TI11*TI4
189 DR3 = CR3-CI4
190 DR4 = CR3+CI4
191 DI3 = CI3+CR4
192 DI4 = CI3-CR4
193 DR5 = CR2+CI5
194 DR2 = CR2-CI5
195 DI5 = CI2-CR5
196 DI2 = CI2+CR5
197 CH(I-1,K,2) = WA1(I-2)*DR2-WA1(I-1)*DI2
198 CH(I,K,2) = WA1(I-2)*DI2+WA1(I-1)*DR2
199 CH(I-1,K,3) = WA2(I-2)*DR3-WA2(I-1)*DI3
200 CH(I,K,3) = WA2(I-2)*DI3+WA2(I-1)*DR3
201 CH(I-1,K,4) = WA3(I-2)*DR4-WA3(I-1)*DI4
202 CH(I,K,4) = WA3(I-2)*DI4+WA3(I-1)*DR4
203 CH(I-1,K,5) = WA4(I-2)*DR5-WA4(I-1)*DI5
204 CH(I,K,5) = WA4(I-2)*DI5+WA4(I-1)*DR5
205 102 CONTINUE
206 103 CONTINUE
207 RETURN
208 END
209 SUBROUTINE RADBG (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
210
211 IMPLICIT REAL*8 (A-H,O-Z)
212 IMPLICIT INTEGER (I-N)
213
214 DIMENSION CH(IDO,L1,IP) ,CC(IDO,IP,L1) ,
215 1 C1(IDO,L1,IP) ,C2(IDL1,IP),
216 2 CH2(IDL1,IP) ,WA(1)
217 DATA TPI/6.28318530717959D0/
218 ARG = TPI/FLOAT(IP)
219 DCP = COS(ARG)
220 DSP = SIN(ARG)
221 IDP2 = IDO+2
222 NBD = (IDO-1)/2
223 IPP2 = IP+2
224 IPPH = (IP+1)/2
225 IF (IDO .LT. L1) GO TO 103
226 DO 102 K=1,L1
227 DO 101 I=1,IDO
228 CH(I,K,1) = CC(I,1,K)
229 101 CONTINUE
230 102 CONTINUE
231 GO TO 106
232 103 DO 105 I=1,IDO
233 DO 104 K=1,L1
234 CH(I,K,1) = CC(I,1,K)
235 104 CONTINUE
236 105 CONTINUE
237 106 DO 108 J=2,IPPH
238 JC = IPP2-J
239 J2 = J+J
240 DO 107 K=1,L1
241 CH(1,K,J) = CC(IDO,J2-2,K)+CC(IDO,J2-2,K)
242 CH(1,K,JC) = CC(1,J2-1,K)+CC(1,J2-1,K)
243 107 CONTINUE
244 108 CONTINUE
245 IF (IDO .EQ. 1) GO TO 116
246 IF (NBD .LT. L1) GO TO 112
247 DO 111 J=2,IPPH
248 JC = IPP2-J
249 DO 110 K=1,L1
250 DO 109 I=3,IDO,2
251 IC = IDP2-I
252 CH(I-1,K,J) = CC(I-1,2*J-1,K)+CC(IC-1,2*J-2,K)
253 CH(I-1,K,JC) = CC(I-1,2*J-1,K)-CC(IC-1,2*J-2,K)
254 CH(I,K,J) = CC(I,2*J-1,K)-CC(IC,2*J-2,K)
255 CH(I,K,JC) = CC(I,2*J-1,K)+CC(IC,2*J-2,K)
256 109 CONTINUE
257 110 CONTINUE
258 111 CONTINUE
259 GO TO 116
260 112 DO 115 J=2,IPPH
261 JC = IPP2-J
262 DO 114 I=3,IDO,2
263 IC = IDP2-I
264 DO 113 K=1,L1
265 CH(I-1,K,J) = CC(I-1,2*J-1,K)+CC(IC-1,2*J-2,K)
266 CH(I-1,K,JC) = CC(I-1,2*J-1,K)-CC(IC-1,2*J-2,K)
267 CH(I,K,J) = CC(I,2*J-1,K)-CC(IC,2*J-2,K)
268 CH(I,K,JC) = CC(I,2*J-1,K)+CC(IC,2*J-2,K)
269 113 CONTINUE
270 114 CONTINUE
271 115 CONTINUE
272 116 AR1 = 1.D0
273 AI1 = 0.D0
274 DO 120 L=2,IPPH
275 LC = IPP2-L
276 AR1H = DCP*AR1-DSP*AI1
277 AI1 = DCP*AI1+DSP*AR1
278 AR1 = AR1H
279 DO 117 IK=1,IDL1
280 C2(IK,L) = CH2(IK,1)+AR1*CH2(IK,2)
281 C2(IK,LC) = AI1*CH2(IK,IP)
282 117 CONTINUE
283 DC2 = AR1
284 DS2 = AI1
285 AR2 = AR1
286 AI2 = AI1
287 DO 119 J=3,IPPH
288 JC = IPP2-J
289 AR2H = DC2*AR2-DS2*AI2
290 AI2 = DC2*AI2+DS2*AR2
291 AR2 = AR2H
292 DO 118 IK=1,IDL1
293 C2(IK,L) = C2(IK,L)+AR2*CH2(IK,J)
294 C2(IK,LC) = C2(IK,LC)+AI2*CH2(IK,JC)
295 118 CONTINUE
296 119 CONTINUE
297 120 CONTINUE
298 DO 122 J=2,IPPH
299 DO 121 IK=1,IDL1
300 CH2(IK,1) = CH2(IK,1)+CH2(IK,J)
301 121 CONTINUE
302 122 CONTINUE
303 DO 124 J=2,IPPH
304 JC = IPP2-J
305 DO 123 K=1,L1
306 CH(1,K,J) = C1(1,K,J)-C1(1,K,JC)
307 CH(1,K,JC) = C1(1,K,J)+C1(1,K,JC)
308 123 CONTINUE
309 124 CONTINUE
310 IF (IDO .EQ. 1) GO TO 132
311 IF (NBD .LT. L1) GO TO 128
312 DO 127 J=2,IPPH
313 JC = IPP2-J
314 DO 126 K=1,L1
315 DO 125 I=3,IDO,2
316 CH(I-1,K,J) = C1(I-1,K,J)-C1(I,K,JC)
317 CH(I-1,K,JC) = C1(I-1,K,J)+C1(I,K,JC)
318 CH(I,K,J) = C1(I,K,J)+C1(I-1,K,JC)
319 CH(I,K,JC) = C1(I,K,J)-C1(I-1,K,JC)
320 125 CONTINUE
321 126 CONTINUE
322 127 CONTINUE
323 GO TO 132
324 128 DO 131 J=2,IPPH
325 JC = IPP2-J
326 DO 130 I=3,IDO,2
327 DO 129 K=1,L1
328 CH(I-1,K,J) = C1(I-1,K,J)-C1(I,K,JC)
329 CH(I-1,K,JC) = C1(I-1,K,J)+C1(I,K,JC)
330 CH(I,K,J) = C1(I,K,J)+C1(I-1,K,JC)
331 CH(I,K,JC) = C1(I,K,J)-C1(I-1,K,JC)
332 129 CONTINUE
333 130 CONTINUE
334 131 CONTINUE
335 132 CONTINUE
336 IF (IDO .EQ. 1) RETURN
337 DO 133 IK=1,IDL1
338 C2(IK,1) = CH2(IK,1)
339 133 CONTINUE
340 DO 135 J=2,IP
341 DO 134 K=1,L1
342 C1(1,K,J) = CH(1,K,J)
343 134 CONTINUE
344 135 CONTINUE
345 IF (NBD .GT. L1) GO TO 139
346 IS = -IDO
347 DO 138 J=2,IP
348 IS = IS+IDO
349 IDIJ = IS
350 DO 137 I=3,IDO,2
351 IDIJ = IDIJ+2
352 DO 136 K=1,L1
353 C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J)
354 C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J)
355 136 CONTINUE
356 137 CONTINUE
357 138 CONTINUE
358 GO TO 143
359 139 IS = -IDO
360 DO 142 J=2,IP
361 IS = IS+IDO
362 DO 141 K=1,L1
363 IDIJ = IS
364 DO 140 I=3,IDO,2
365 IDIJ = IDIJ+2
366 C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J)
367 C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J)
368 140 CONTINUE
369 141 CONTINUE
370 142 CONTINUE
371 143 RETURN
372 END
373 SUBROUTINE RFFTB1 (N,C,CH,WA,IFAC)
374
375 IMPLICIT REAL*8 (A-H,O-Z)
376 IMPLICIT INTEGER (I-N)
377
378 DIMENSION CH(1) ,C(1) ,WA(1) ,IFAC(1)
379 NF = IFAC(2)
380 NA = 0
381 L1 = 1
382 IW = 1
383 DO 116 K1=1,NF
384 IP = IFAC(K1+2)
385 L2 = IP*L1
386 IDO = N/L2
387 IDL1 = IDO*L1
388 IF (IP .NE. 4) GO TO 103
389 IX2 = IW+IDO
390 IX3 = IX2+IDO
391 IF (NA .NE. 0) GO TO 101
392 CALL RADB4 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3))
393 GO TO 102
394 101 CALL RADB4 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3))
395 102 NA = 1-NA
396 GO TO 115
397 103 IF (IP .NE. 2) GO TO 106
398 IF (NA .NE. 0) GO TO 104
399 CALL RADB2 (IDO,L1,C,CH,WA(IW))
400 GO TO 105
401 104 CALL RADB2 (IDO,L1,CH,C,WA(IW))
402 105 NA = 1-NA
403 GO TO 115
404 106 IF (IP .NE. 3) GO TO 109
405 IX2 = IW+IDO
406 IF (NA .NE. 0) GO TO 107
407 CALL RADB3 (IDO,L1,C,CH,WA(IW),WA(IX2))
408 GO TO 108
409 107 CALL RADB3 (IDO,L1,CH,C,WA(IW),WA(IX2))
410 108 NA = 1-NA
411 GO TO 115
412 109 IF (IP .NE. 5) GO TO 112
413 IX2 = IW+IDO
414 IX3 = IX2+IDO
415 IX4 = IX3+IDO
416 IF (NA .NE. 0) GO TO 110
417 CALL RADB5 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4))
418 GO TO 111
419 110 CALL RADB5 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4))
420 111 NA = 1-NA
421 GO TO 115
422 112 IF (NA .NE. 0) GO TO 113
423 CALL RADBG (IDO,IP,L1,IDL1,C,C,C,CH,CH,WA(IW))
424 GO TO 114
425 113 CALL RADBG (IDO,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW))
426 114 IF (IDO .EQ. 1) NA = 1-NA
427 115 L1 = L2
428 IW = IW+(IP-1)*IDO
429 116 CONTINUE
430 IF (NA .EQ. 0) RETURN
431 DO 117 I=1,N
432 C(I) = CH(I)
433 117 CONTINUE
434 RETURN
435 END
436 C SUBROUTINE RFFTF (N,R,WSAVE)
437
438 C IMPLICIT REAL*8 (A-H,O-Z)
439
440 C DIMENSION R(1) ,WSAVE(1)
441 C IF (N .EQ. 1) RETURN
442 C CALL RFFTF1 (N,R,WSAVE,WSAVE(N+1),WSAVE(2*N+1))
443 C RETURN
444 C END
445 SUBROUTINE RADF2 (IDO,L1,CC,CH,WA1)
446
447 IMPLICIT REAL*8 (A-H,O-Z)
448 IMPLICIT INTEGER (I-N)
449
450 DIMENSION CH(IDO,2,L1) ,CC(IDO,L1,2) ,
451 1 WA1(1)
452 DO 101 K=1,L1
453 CH(1,1,K) = CC(1,K,1)+CC(1,K,2)
454 CH(IDO,2,K) = CC(1,K,1)-CC(1,K,2)
455 101 CONTINUE
456 IF (IDO-2) 107,105,102
457 102 IDP2 = IDO+2
458 DO 104 K=1,L1
459 DO 103 I=3,IDO,2
460 IC = IDP2-I
461 TR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
462 TI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
463 CH(I,1,K) = CC(I,K,1)+TI2
464 CH(IC,2,K) = TI2-CC(I,K,1)
465 CH(I-1,1,K) = CC(I-1,K,1)+TR2
466 CH(IC-1,2,K) = CC(I-1,K,1)-TR2
467 103 CONTINUE
468 104 CONTINUE
469 IF (MOD(IDO,2) .EQ. 1) RETURN
470 105 DO 106 K=1,L1
471 CH(1,2,K) = -CC(IDO,K,2)
472 CH(IDO,1,K) = CC(IDO,K,1)
473 106 CONTINUE
474 107 RETURN
475 END
476 SUBROUTINE RADF3 (IDO,L1,CC,CH,WA1,WA2)
477
478 IMPLICIT REAL*8 (A-H,O-Z)
479 IMPLICIT INTEGER (I-N)
480
481 DIMENSION CH(IDO,3,L1) ,CC(IDO,L1,3) ,
482 1 WA1(1) ,WA2(1)
483 DATA TAUR,TAUI /-.5D0,.866025403784439D0/
484 DO 101 K=1,L1
485 CR2 = CC(1,K,2)+CC(1,K,3)
486 CH(1,1,K) = CC(1,K,1)+CR2
487 CH(1,3,K) = TAUI*(CC(1,K,3)-CC(1,K,2))
488 CH(IDO,2,K) = CC(1,K,1)+TAUR*CR2
489 101 CONTINUE
490 IF (IDO .EQ. 1) RETURN
491 IDP2 = IDO+2
492 DO 103 K=1,L1
493 DO 102 I=3,IDO,2
494 IC = IDP2-I
495 DR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
496 DI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
497 DR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3)
498 DI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3)
499 CR2 = DR2+DR3
500 CI2 = DI2+DI3
501 CH(I-1,1,K) = CC(I-1,K,1)+CR2
502 CH(I,1,K) = CC(I,K,1)+CI2
503 TR2 = CC(I-1,K,1)+TAUR*CR2
504 TI2 = CC(I,K,1)+TAUR*CI2
505 TR3 = TAUI*(DI2-DI3)
506 TI3 = TAUI*(DR3-DR2)
507 CH(I-1,3,K) = TR2+TR3
508 CH(IC-1,2,K) = TR2-TR3
509 CH(I,3,K) = TI2+TI3
510 CH(IC,2,K) = TI3-TI2
511 102 CONTINUE
512 103 CONTINUE
513 RETURN
514 END
515 SUBROUTINE RADF4 (IDO,L1,CC,CH,WA1,WA2,WA3)
516
517 IMPLICIT REAL*8 (A-H,O-Z)
518 IMPLICIT INTEGER (I-N)
519
520 DIMENSION CC(IDO,L1,4) ,CH(IDO,4,L1) ,
521 1 WA1(1) ,WA2(1) ,WA3(1)
522 DATA HSQT2 /.7071067811865475D0/
523 DO 101 K=1,L1
524 TR1 = CC(1,K,2)+CC(1,K,4)
525 TR2 = CC(1,K,1)+CC(1,K,3)
526 CH(1,1,K) = TR1+TR2
527 CH(IDO,4,K) = TR2-TR1
528 CH(IDO,2,K) = CC(1,K,1)-CC(1,K,3)
529 CH(1,3,K) = CC(1,K,4)-CC(1,K,2)
530 101 CONTINUE
531 IF (IDO-2) 107,105,102
532 102 IDP2 = IDO+2
533 DO 104 K=1,L1
534 DO 103 I=3,IDO,2
535 IC = IDP2-I
536 CR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
537 CI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
538 CR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3)
539 CI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3)
540 CR4 = WA3(I-2)*CC(I-1,K,4)+WA3(I-1)*CC(I,K,4)
541 CI4 = WA3(I-2)*CC(I,K,4)-WA3(I-1)*CC(I-1,K,4)
542 TR1 = CR2+CR4
543 TR4 = CR4-CR2
544 TI1 = CI2+CI4
545 TI4 = CI2-CI4
546 TI2 = CC(I,K,1)+CI3
547 TI3 = CC(I,K,1)-CI3
548 TR2 = CC(I-1,K,1)+CR3
549 TR3 = CC(I-1,K,1)-CR3
550 CH(I-1,1,K) = TR1+TR2
551 CH(IC-1,4,K) = TR2-TR1
552 CH(I,1,K) = TI1+TI2
553 CH(IC,4,K) = TI1-TI2
554 CH(I-1,3,K) = TI4+TR3
555 CH(IC-1,2,K) = TR3-TI4
556 CH(I,3,K) = TR4+TI3
557 CH(IC,2,K) = TR4-TI3
558 103 CONTINUE
559 104 CONTINUE
560 IF (MOD(IDO,2) .EQ. 1) RETURN
561 105 CONTINUE
562 DO 106 K=1,L1
563 TI1 = -HSQT2*(CC(IDO,K,2)+CC(IDO,K,4))
564 TR1 = HSQT2*(CC(IDO,K,2)-CC(IDO,K,4))
565 CH(IDO,1,K) = TR1+CC(IDO,K,1)
566 CH(IDO,3,K) = CC(IDO,K,1)-TR1
567 CH(1,2,K) = TI1-CC(IDO,K,3)
568 CH(1,4,K) = TI1+CC(IDO,K,3)
569 106 CONTINUE
570 107 RETURN
571 END
572 SUBROUTINE RADF5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
573
574 IMPLICIT REAL*8 (A-H,O-Z)
575 IMPLICIT INTEGER (I-N)
576
577 DIMENSION CC(IDO,L1,5) ,CH(IDO,5,L1) ,
578 1 WA1(1) ,WA2(1) ,WA3(1) ,WA4(1)
579 DATA TR11,TI11,TR12,TI12 /.309016994374947D0,.951056516295154D0,
580 1-.809016994374947D0,.587785252292473D0/
581 DO 101 K=1,L1
582 CR2 = CC(1,K,5)+CC(1,K,2)
583 CI5 = CC(1,K,5)-CC(1,K,2)
584 CR3 = CC(1,K,4)+CC(1,K,3)
585 CI4 = CC(1,K,4)-CC(1,K,3)
586 CH(1,1,K) = CC(1,K,1)+CR2+CR3
587 CH(IDO,2,K) = CC(1,K,1)+TR11*CR2+TR12*CR3
588 CH(1,3,K) = TI11*CI5+TI12*CI4
589 CH(IDO,4,K) = CC(1,K,1)+TR12*CR2+TR11*CR3
590 CH(1,5,K) = TI12*CI5-TI11*CI4
591 101 CONTINUE
592 IF (IDO .EQ. 1) RETURN
593 IDP2 = IDO+2
594 DO 103 K=1,L1
595 DO 102 I=3,IDO,2
596 IC = IDP2-I
597 DR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
598 DI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
599 DR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3)
600 DI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3)
601 DR4 = WA3(I-2)*CC(I-1,K,4)+WA3(I-1)*CC(I,K,4)
602 DI4 = WA3(I-2)*CC(I,K,4)-WA3(I-1)*CC(I-1,K,4)
603 DR5 = WA4(I-2)*CC(I-1,K,5)+WA4(I-1)*CC(I,K,5)
604 DI5 = WA4(I-2)*CC(I,K,5)-WA4(I-1)*CC(I-1,K,5)
605 CR2 = DR2+DR5
606 CI5 = DR5-DR2
607 CR5 = DI2-DI5
608 CI2 = DI2+DI5
609 CR3 = DR3+DR4
610 CI4 = DR4-DR3
611 CR4 = DI3-DI4
612 CI3 = DI3+DI4
613 CH(I-1,1,K) = CC(I-1,K,1)+CR2+CR3
614 CH(I,1,K) = CC(I,K,1)+CI2+CI3
615 TR2 = CC(I-1,K,1)+TR11*CR2+TR12*CR3
616 TI2 = CC(I,K,1)+TR11*CI2+TR12*CI3
617 TR3 = CC(I-1,K,1)+TR12*CR2+TR11*CR3
618 TI3 = CC(I,K,1)+TR12*CI2+TR11*CI3
619 TR5 = TI11*CR5+TI12*CR4
620 TI5 = TI11*CI5+TI12*CI4
621 TR4 = TI12*CR5-TI11*CR4
622 TI4 = TI12*CI5-TI11*CI4
623 CH(I-1,3,K) = TR2+TR5
624 CH(IC-1,2,K) = TR2-TR5
625 CH(I,3,K) = TI2+TI5
626 CH(IC,2,K) = TI5-TI2
627 CH(I-1,5,K) = TR3+TR4
628 CH(IC-1,4,K) = TR3-TR4
629 CH(I,5,K) = TI3+TI4
630 CH(IC,4,K) = TI4-TI3
631 102 CONTINUE
632 103 CONTINUE
633 RETURN
634 END
635 SUBROUTINE RADFG (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
636
637 IMPLICIT REAL*8 (A-H,O-Z)
638 IMPLICIT INTEGER (I-N)
639
640 DIMENSION CH(IDO,L1,IP) ,CC(IDO,IP,L1) ,
641 1 C1(IDO,L1,IP) ,C2(IDL1,IP),
642 2 CH2(IDL1,IP) ,WA(1)
643 DATA TPI/6.28318530717959D0/
644 ARG = TPI/FLOAT(IP)
645 DCP = COS(ARG)
646 DSP = SIN(ARG)
647 IPPH = (IP+1)/2
648 IPP2 = IP+2
649 IDP2 = IDO+2
650 NBD = (IDO-1)/2
651 IF (IDO .EQ. 1) GO TO 119
652 DO 101 IK=1,IDL1
653 CH2(IK,1) = C2(IK,1)
654 101 CONTINUE
655 DO 103 J=2,IP
656 DO 102 K=1,L1
657 CH(1,K,J) = C1(1,K,J)
658 102 CONTINUE
659 103 CONTINUE
660 IF (NBD .GT. L1) GO TO 107
661 IS = -IDO
662 DO 106 J=2,IP
663 IS = IS+IDO
664 IDIJ = IS
665 DO 105 I=3,IDO,2
666 IDIJ = IDIJ+2
667 DO 104 K=1,L1
668 CH(I-1,K,J) = WA(IDIJ-1)*C1(I-1,K,J)+WA(IDIJ)*C1(I,K,J)
669 CH(I,K,J) = WA(IDIJ-1)*C1(I,K,J)-WA(IDIJ)*C1(I-1,K,J)
670 104 CONTINUE
671 105 CONTINUE
672 106 CONTINUE
673 GO TO 111
674 107 IS = -IDO
675 DO 110 J=2,IP
676 IS = IS+IDO
677 DO 109 K=1,L1
678 IDIJ = IS
679 DO 108 I=3,IDO,2
680 IDIJ = IDIJ+2
681 CH(I-1,K,J) = WA(IDIJ-1)*C1(I-1,K,J)+WA(IDIJ)*C1(I,K,J)
682 CH(I,K,J) = WA(IDIJ-1)*C1(I,K,J)-WA(IDIJ)*C1(I-1,K,J)
683 108 CONTINUE
684 109 CONTINUE
685 110 CONTINUE
686 111 IF (NBD .LT. L1) GO TO 115
687 DO 114 J=2,IPPH
688 JC = IPP2-J
689 DO 113 K=1,L1
690 DO 112 I=3,IDO,2
691 C1(I-1,K,J) = CH(I-1,K,J)+CH(I-1,K,JC)
692 C1(I-1,K,JC) = CH(I,K,J)-CH(I,K,JC)
693 C1(I,K,J) = CH(I,K,J)+CH(I,K,JC)
694 C1(I,K,JC) = CH(I-1,K,JC)-CH(I-1,K,J)
695 112 CONTINUE
696 113 CONTINUE
697 114 CONTINUE
698 GO TO 121
699 115 DO 118 J=2,IPPH
700 JC = IPP2-J
701 DO 117 I=3,IDO,2
702 DO 116 K=1,L1
703 C1(I-1,K,J) = CH(I-1,K,J)+CH(I-1,K,JC)
704 C1(I-1,K,JC) = CH(I,K,J)-CH(I,K,JC)
705 C1(I,K,J) = CH(I,K,J)+CH(I,K,JC)
706 C1(I,K,JC) = CH(I-1,K,JC)-CH(I-1,K,J)
707 116 CONTINUE
708 117 CONTINUE
709 118 CONTINUE
710 GO TO 121
711 119 DO 120 IK=1,IDL1
712 C2(IK,1) = CH2(IK,1)
713 120 CONTINUE
714 121 DO 123 J=2,IPPH
715 JC = IPP2-J
716 DO 122 K=1,L1
717 C1(1,K,J) = CH(1,K,J)+CH(1,K,JC)
718 C1(1,K,JC) = CH(1,K,JC)-CH(1,K,J)
719 122 CONTINUE
720 123 CONTINUE
721 C
722 AR1 = 1.D0
723 AI1 = 0.D0
724 DO 127 L=2,IPPH
725 LC = IPP2-L
726 AR1H = DCP*AR1-DSP*AI1
727 AI1 = DCP*AI1+DSP*AR1
728 AR1 = AR1H
729 DO 124 IK=1,IDL1
730 CH2(IK,L) = C2(IK,1)+AR1*C2(IK,2)
731 CH2(IK,LC) = AI1*C2(IK,IP)
732 124 CONTINUE
733 DC2 = AR1
734 DS2 = AI1
735 AR2 = AR1
736 AI2 = AI1
737 DO 126 J=3,IPPH
738 JC = IPP2-J
739 AR2H = DC2*AR2-DS2*AI2
740 AI2 = DC2*AI2+DS2*AR2
741 AR2 = AR2H
742 DO 125 IK=1,IDL1
743 CH2(IK,L) = CH2(IK,L)+AR2*C2(IK,J)
744 CH2(IK,LC) = CH2(IK,LC)+AI2*C2(IK,JC)
745 125 CONTINUE
746 126 CONTINUE
747 127 CONTINUE
748 DO 129 J=2,IPPH
749 DO 128 IK=1,IDL1
750 CH2(IK,1) = CH2(IK,1)+C2(IK,J)
751 128 CONTINUE
752 129 CONTINUE
753 C
754 IF (IDO .LT. L1) GO TO 132
755 DO 131 K=1,L1
756 DO 130 I=1,IDO
757 CC(I,1,K) = CH(I,K,1)
758 130 CONTINUE
759 131 CONTINUE
760 GO TO 135
761 132 DO 134 I=1,IDO
762 DO 133 K=1,L1
763 CC(I,1,K) = CH(I,K,1)
764 133 CONTINUE
765 134 CONTINUE
766 135 DO 137 J=2,IPPH
767 JC = IPP2-J
768 J2 = J+J
769 DO 136 K=1,L1
770 CC(IDO,J2-2,K) = CH(1,K,J)
771 CC(1,J2-1,K) = CH(1,K,JC)
772 136 CONTINUE
773 137 CONTINUE
774 IF (IDO .EQ. 1) RETURN
775 IF (NBD .LT. L1) GO TO 141
776 DO 140 J=2,IPPH
777 JC = IPP2-J
778 J2 = J+J
779 DO 139 K=1,L1
780 DO 138 I=3,IDO,2
781 IC = IDP2-I
782 CC(I-1,J2-1,K) = CH(I-1,K,J)+CH(I-1,K,JC)
783 CC(IC-1,J2-2,K) = CH(I-1,K,J)-CH(I-1,K,JC)
784 CC(I,J2-1,K) = CH(I,K,J)+CH(I,K,JC)
785 CC(IC,J2-2,K) = CH(I,K,JC)-CH(I,K,J)
786 138 CONTINUE
787 139 CONTINUE
788 140 CONTINUE
789 RETURN
790 141 DO 144 J=2,IPPH
791 JC = IPP2-J
792 J2 = J+J
793 DO 143 I=3,IDO,2
794 IC = IDP2-I
795 DO 142 K=1,L1
796 CC(I-1,J2-1,K) = CH(I-1,K,J)+CH(I-1,K,JC)
797 CC(IC-1,J2-2,K) = CH(I-1,K,J)-CH(I-1,K,JC)
798 CC(I,J2-1,K) = CH(I,K,J)+CH(I,K,JC)
799 CC(IC,J2-2,K) = CH(I,K,JC)-CH(I,K,J)
800 142 CONTINUE
801 143 CONTINUE
802 144 CONTINUE
803 RETURN
804 END
805 SUBROUTINE RFFTF1 (N,C,CH,WA,IFAC)
806
807 IMPLICIT REAL*8 (A-H,O-Z)
808 IMPLICIT INTEGER (I-N)
809
810 DIMENSION CH(1) ,C(1) ,WA(1) ,IFAC(1)
811 NF = IFAC(2)
812 NA = 1
813 L2 = N
814 IW = N
815 DO 111 K1=1,NF
816 KH = NF-K1
817 IP = IFAC(KH+3)
818 L1 = L2/IP
819 IDO = N/L2
820 IDL1 = IDO*L1
821 IW = IW-(IP-1)*IDO
822 NA = 1-NA
823 IF (IP .NE. 4) GO TO 102
824 IX2 = IW+IDO
825 IX3 = IX2+IDO
826 IF (NA .NE. 0) GO TO 101
827 CALL RADF4 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3))
828 GO TO 110
829 101 CALL RADF4 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3))
830 GO TO 110
831 102 IF (IP .NE. 2) GO TO 104
832 IF (NA .NE. 0) GO TO 103
833 CALL RADF2 (IDO,L1,C,CH,WA(IW))
834 GO TO 110
835 103 CALL RADF2 (IDO,L1,CH,C,WA(IW))
836 GO TO 110
837 104 IF (IP .NE. 3) GO TO 106
838 IX2 = IW+IDO
839 IF (NA .NE. 0) GO TO 105
840 CALL RADF3 (IDO,L1,C,CH,WA(IW),WA(IX2))
841 GO TO 110
842 105 CALL RADF3 (IDO,L1,CH,C,WA(IW),WA(IX2))
843 GO TO 110
844 106 IF (IP .NE. 5) GO TO 108
845 IX2 = IW+IDO
846 IX3 = IX2+IDO
847 IX4 = IX3+IDO
848 IF (NA .NE. 0) GO TO 107
849 CALL RADF5 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4))
850 GO TO 110
851 107 CALL RADF5 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4))
852 GO TO 110
853 108 IF (IDO .EQ. 1) NA = 1-NA
854 IF (NA .NE. 0) GO TO 109
855 CALL RADFG (IDO,IP,L1,IDL1,C,C,C,CH,CH,WA(IW))
856 NA = 1
857 GO TO 110
858 109 CALL RADFG (IDO,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW))
859 NA = 0
860 110 L2 = L1
861 111 CONTINUE
862 IF (NA .EQ. 1) RETURN
863 DO 112 I=1,N
864 C(I) = CH(I)
865 112 CONTINUE
866 RETURN
867 END
868 C SUBROUTINE RFFTI (N,WSAVE)
869
870 C IMPLICIT REAL*8 (A-H,O-Z)
871 C IMPLICIT INTEGER (I-N)
872
873 C DIMENSION WSAVE(1)
874 C IF (N .EQ. 1) RETURN
875 C CALL RFFTI1 (N,WSAVE(N+1),WSAVE(2*N+1))
876 C RETURN
877 C END
878 SUBROUTINE RFFTI1 (N,WA,IFAC)
879
880 IMPLICIT REAL*8 (A-H,O-Z)
881 IMPLICIT INTEGER (I-N)
882
883 DIMENSION WA(1) ,IFAC(1) ,NTRYH(4)
884 DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/4,2,3,5/
885 NL = N
886 NF = 0
887 J = 0
888 101 J = J+1
889 IF (J-4) 102,102,103
890 102 NTRY = NTRYH(J)
891 GO TO 104
892 103 NTRY = NTRY+2
893 104 NQ = NL/NTRY
894 NR = NL-NTRY*NQ
895 IF (NR) 101,105,101
896 105 NF = NF+1
897 IFAC(NF+2) = NTRY
898 NL = NQ
899 IF (NTRY .NE. 2) GO TO 107
900 IF (NF .EQ. 1) GO TO 107
901 DO 106 I=2,NF
902 IB = NF-I+2
903 IFAC(IB+2) = IFAC(IB+1)
904 106 CONTINUE
905 IFAC(3) = 2
906 107 IF (NL .NE. 1) GO TO 104
907 IFAC(1) = N
908 IFAC(2) = NF
909 TPI = 6.28318530717959D0
910 ARGH = TPI/FLOAT(N)
911 IS = 0
912 NFM1 = NF-1
913 L1 = 1
914 IF (NFM1 .EQ. 0) RETURN
915 DO 110 K1=1,NFM1
916 IP = IFAC(K1+2)
917 LD = 0
918 L2 = L1*IP
919 IDO = N/L2
920 IPM = IP-1
921 DO 109 J=1,IPM
922 LD = LD+L1
923 I = IS
924 ARGLD = FLOAT(LD)*ARGH
925 FI = 0.D0
926 DO 108 II=3,IDO,2
927 I = I+2
928 FI = FI+1.D0
929 ARG = FI*ARGLD
930 WA(I-1) = COS(ARG)
931 WA(I) = SIN(ARG)
932 108 CONTINUE
933 IS = IS+IDO
934 109 CONTINUE
935 L1 = L2
936 110 CONTINUE
937 RETURN
938 END
939 C SUBROUTINE R8FFTB (N,R,WSAVE)
940 C implicit real*8 (A-H,O-Z)
941 C DIMENSION R(1) ,WSAVE(1)
942 C IF (N .EQ. 1) RETURN
943 C CALL R8FFTB1 (N,R,WSAVE,WSAVE(N+1),WSAVE(2*N+1))
944 C RETURN
945 C END
946 SUBROUTINE R8ADB2 (IDO,L1,CC,CH,WA1)
947 implicit real*8 (A-H,O-Z)
948 IMPLICIT INTEGER (I-N)
949 DIMENSION CC(IDO,2,L1) ,CH(IDO,L1,2) ,
950 1 WA1(1)
951 DO 101 K=1,L1
952 CH(1,K,1) = CC(1,1,K)+CC(IDO,2,K)
953 CH(1,K,2) = CC(1,1,K)-CC(IDO,2,K)
954 101 CONTINUE
955 IF (IDO-2) 107,105,102
956 102 IDP2 = IDO+2
957 DO 104 K=1,L1
958 DO 103 I=3,IDO,2
959 IC = IDP2-I
960 CH(I-1,K,1) = CC(I-1,1,K)+CC(IC-1,2,K)
961 TR2 = CC(I-1,1,K)-CC(IC-1,2,K)
962 CH(I,K,1) = CC(I,1,K)-CC(IC,2,K)
963 TI2 = CC(I,1,K)+CC(IC,2,K)
964 CH(I-1,K,2) = WA1(I-2)*TR2-WA1(I-1)*TI2
965 CH(I,K,2) = WA1(I-2)*TI2+WA1(I-1)*TR2
966 103 CONTINUE
967 104 CONTINUE
968 IF (MOD(IDO,2) .EQ. 1) RETURN
969 105 DO 106 K=1,L1
970 CH(IDO,K,1) = CC(IDO,1,K)+CC(IDO,1,K)
971 CH(IDO,K,2) = -(CC(1,2,K)+CC(1,2,K))
972 106 CONTINUE
973 107 RETURN
974 END
975 SUBROUTINE R8ADB3 (IDO,L1,CC,CH,WA1,WA2)
976 implicit real*8 (A-H,O-Z)
977 IMPLICIT INTEGER (I-N)
978 DIMENSION CC(IDO,3,L1) ,CH(IDO,L1,3) ,
979 1 WA1(1) ,WA2(1)
980 DATA TAUR,TAUI /-.5D0,.866025403784439D0/
981 DO 101 K=1,L1
982 TR2 = CC(IDO,2,K)+CC(IDO,2,K)
983 CR2 = CC(1,1,K)+TAUR*TR2
984 CH(1,K,1) = CC(1,1,K)+TR2
985 CI3 = TAUI*(CC(1,3,K)+CC(1,3,K))
986 CH(1,K,2) = CR2-CI3
987 CH(1,K,3) = CR2+CI3
988 101 CONTINUE
989 IF (IDO .EQ. 1) RETURN
990 IDP2 = IDO+2
991 DO 103 K=1,L1
992 DO 102 I=3,IDO,2
993 IC = IDP2-I
994 TR2 = CC(I-1,3,K)+CC(IC-1,2,K)
995 CR2 = CC(I-1,1,K)+TAUR*TR2
996 CH(I-1,K,1) = CC(I-1,1,K)+TR2
997 TI2 = CC(I,3,K)-CC(IC,2,K)
998 CI2 = CC(I,1,K)+TAUR*TI2
999 CH(I,K,1) = CC(I,1,K)+TI2
1000 CR3 = TAUI*(CC(I-1,3,K)-CC(IC-1,2,K))
1001 CI3 = TAUI*(CC(I,3,K)+CC(IC,2,K))
1002 DR2 = CR2-CI3
1003 DR3 = CR2+CI3
1004 DI2 = CI2+CR3
1005 DI3 = CI2-CR3
1006 CH(I-1,K,2) = WA1(I-2)*DR2-WA1(I-1)*DI2
1007 CH(I,K,2) = WA1(I-2)*DI2+WA1(I-1)*DR2
1008 CH(I-1,K,3) = WA2(I-2)*DR3-WA2(I-1)*DI3
1009 CH(I,K,3) = WA2(I-2)*DI3+WA2(I-1)*DR3
1010 102 CONTINUE
1011 103 CONTINUE
1012 RETURN
1013 END
1014 SUBROUTINE R8ADB4 (IDO,L1,CC,CH,WA1,WA2,WA3)
1015 implicit real*8 (A-H,O-Z)
1016 IMPLICIT INTEGER (I-N)
1017 DIMENSION CC(IDO,4,L1) ,CH(IDO,L1,4) ,
1018 1 WA1(1) ,WA2(1) ,WA3(1)
1019 DATA SQRT2 /1.414213562373095D0/
1020 DO 101 K=1,L1
1021 TR1 = CC(1,1,K)-CC(IDO,4,K)
1022 TR2 = CC(1,1,K)+CC(IDO,4,K)
1023 TR3 = CC(IDO,2,K)+CC(IDO,2,K)
1024 TR4 = CC(1,3,K)+CC(1,3,K)
1025 CH(1,K,1) = TR2+TR3
1026 CH(1,K,2) = TR1-TR4
1027 CH(1,K,3) = TR2-TR3
1028 CH(1,K,4) = TR1+TR4
1029 101 CONTINUE
1030 IF (IDO-2) 107,105,102
1031 102 IDP2 = IDO+2
1032 DO 104 K=1,L1
1033 DO 103 I=3,IDO,2
1034 IC = IDP2-I
1035 TI1 = CC(I,1,K)+CC(IC,4,K)
1036 TI2 = CC(I,1,K)-CC(IC,4,K)
1037 TI3 = CC(I,3,K)-CC(IC,2,K)
1038 TR4 = CC(I,3,K)+CC(IC,2,K)
1039 TR1 = CC(I-1,1,K)-CC(IC-1,4,K)
1040 TR2 = CC(I-1,1,K)+CC(IC-1,4,K)
1041 TI4 = CC(I-1,3,K)-CC(IC-1,2,K)
1042 TR3 = CC(I-1,3,K)+CC(IC-1,2,K)
1043 CH(I-1,K,1) = TR2+TR3
1044 CR3 = TR2-TR3
1045 CH(I,K,1) = TI2+TI3
1046 CI3 = TI2-TI3
1047 CR2 = TR1-TR4
1048 CR4 = TR1+TR4
1049 CI2 = TI1+TI4
1050 CI4 = TI1-TI4
1051 CH(I-1,K,2) = WA1(I-2)*CR2-WA1(I-1)*CI2
1052 CH(I,K,2) = WA1(I-2)*CI2+WA1(I-1)*CR2
1053 CH(I-1,K,3) = WA2(I-2)*CR3-WA2(I-1)*CI3
1054 CH(I,K,3) = WA2(I-2)*CI3+WA2(I-1)*CR3
1055 CH(I-1,K,4) = WA3(I-2)*CR4-WA3(I-1)*CI4
1056 CH(I,K,4) = WA3(I-2)*CI4+WA3(I-1)*CR4
1057 103 CONTINUE
1058 104 CONTINUE
1059 IF (MOD(IDO,2) .EQ. 1) RETURN
1060 105 CONTINUE
1061 DO 106 K=1,L1
1062 TI1 = CC(1,2,K)+CC(1,4,K)
1063 TI2 = CC(1,4,K)-CC(1,2,K)
1064 TR1 = CC(IDO,1,K)-CC(IDO,3,K)
1065 TR2 = CC(IDO,1,K)+CC(IDO,3,K)
1066 CH(IDO,K,1) = TR2+TR2
1067 CH(IDO,K,2) = SQRT2*(TR1-TI1)
1068 CH(IDO,K,3) = TI2+TI2
1069 CH(IDO,K,4) = -SQRT2*(TR1+TI1)
1070 106 CONTINUE
1071 107 RETURN
1072 END
1073 SUBROUTINE R8ADB5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
1074 implicit real*8 (A-H,O-Z)
1075 IMPLICIT INTEGER (I-N)
1076 DIMENSION CC(IDO,5,L1) ,CH(IDO,L1,5) ,
1077 1 WA1(1) ,WA2(1) ,WA3(1) ,WA4(1)
1078 DATA TR11,TI11,TR12,TI12 /.309016994374947D0,.951056516295154D0,
1079 1-.809016994374947D0,.587785252292473D0/
1080 DO 101 K=1,L1
1081 TI5 = CC(1,3,K)+CC(1,3,K)
1082 TI4 = CC(1,5,K)+CC(1,5,K)
1083 TR2 = CC(IDO,2,K)+CC(IDO,2,K)
1084 TR3 = CC(IDO,4,K)+CC(IDO,4,K)
1085 CH(1,K,1) = CC(1,1,K)+TR2+TR3
1086 CR2 = CC(1,1,K)+TR11*TR2+TR12*TR3
1087 CR3 = CC(1,1,K)+TR12*TR2+TR11*TR3
1088 CI5 = TI11*TI5+TI12*TI4
1089 CI4 = TI12*TI5-TI11*TI4
1090 CH(1,K,2) = CR2-CI5
1091 CH(1,K,3) = CR3-CI4
1092 CH(1,K,4) = CR3+CI4
1093 CH(1,K,5) = CR2+CI5
1094 101 CONTINUE
1095 IF (IDO .EQ. 1) RETURN
1096 IDP2 = IDO+2
1097 DO 103 K=1,L1
1098 DO 102 I=3,IDO,2
1099 IC = IDP2-I
1100 TI5 = CC(I,3,K)+CC(IC,2,K)
1101 TI2 = CC(I,3,K)-CC(IC,2,K)
1102 TI4 = CC(I,5,K)+CC(IC,4,K)
1103 TI3 = CC(I,5,K)-CC(IC,4,K)
1104 TR5 = CC(I-1,3,K)-CC(IC-1,2,K)
1105 TR2 = CC(I-1,3,K)+CC(IC-1,2,K)
1106 TR4 = CC(I-1,5,K)-CC(IC-1,4,K)
1107 TR3 = CC(I-1,5,K)+CC(IC-1,4,K)
1108 CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3
1109 CH(I,K,1) = CC(I,1,K)+TI2+TI3
1110 CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3
1111 CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3
1112 CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3
1113 CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3
1114 CR5 = TI11*TR5+TI12*TR4
1115 CI5 = TI11*TI5+TI12*TI4
1116 CR4 = TI12*TR5-TI11*TR4
1117 CI4 = TI12*TI5-TI11*TI4
1118 DR3 = CR3-CI4
1119 DR4 = CR3+CI4
1120 DI3 = CI3+CR4
1121 DI4 = CI3-CR4
1122 DR5 = CR2+CI5
1123 DR2 = CR2-CI5
1124 DI5 = CI2-CR5
1125 DI2 = CI2+CR5
1126 CH(I-1,K,2) = WA1(I-2)*DR2-WA1(I-1)*DI2
1127 CH(I,K,2) = WA1(I-2)*DI2+WA1(I-1)*DR2
1128 CH(I-1,K,3) = WA2(I-2)*DR3-WA2(I-1)*DI3
1129 CH(I,K,3) = WA2(I-2)*DI3+WA2(I-1)*DR3
1130 CH(I-1,K,4) = WA3(I-2)*DR4-WA3(I-1)*DI4
1131 CH(I,K,4) = WA3(I-2)*DI4+WA3(I-1)*DR4
1132 CH(I-1,K,5) = WA4(I-2)*DR5-WA4(I-1)*DI5
1133 CH(I,K,5) = WA4(I-2)*DI5+WA4(I-1)*DR5
1134 102 CONTINUE
1135 103 CONTINUE
1136 RETURN
1137 END
1138 SUBROUTINE R8ADBG (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
1139 implicit real*8 (A-H,O-Z)
1140 IMPLICIT INTEGER (I-N)
1141 DIMENSION CH(IDO,L1,IP) ,CC(IDO,IP,L1) ,
1142 1 C1(IDO,L1,IP) ,C2(IDL1,IP),
1143 2 CH2(IDL1,IP) ,WA(1)
1144 DATA TPI/6.28318530717959D0/
1145 ARG = TPI/FLOAT(IP)
1146 DCP = COS(ARG)
1147 DSP = SIN(ARG)
1148 IDP2 = IDO+2
1149 NBD = (IDO-1)/2
1150 IPP2 = IP+2
1151 IPPH = (IP+1)/2
1152 IF (IDO .LT. L1) GO TO 103
1153 DO 102 K=1,L1
1154 DO 101 I=1,IDO
1155 CH(I,K,1) = CC(I,1,K)
1156 101 CONTINUE
1157 102 CONTINUE
1158 GO TO 106
1159 103 DO 105 I=1,IDO
1160 DO 104 K=1,L1
1161 CH(I,K,1) = CC(I,1,K)
1162 104 CONTINUE
1163 105 CONTINUE
1164 106 DO 108 J=2,IPPH
1165 JC = IPP2-J
1166 J2 = J+J
1167 DO 107 K=1,L1
1168 CH(1,K,J) = CC(IDO,J2-2,K)+CC(IDO,J2-2,K)
1169 CH(1,K,JC) = CC(1,J2-1,K)+CC(1,J2-1,K)
1170 107 CONTINUE
1171 108 CONTINUE
1172 IF (IDO .EQ. 1) GO TO 116
1173 IF (NBD .LT. L1) GO TO 112
1174 DO 111 J=2,IPPH
1175 JC = IPP2-J
1176 DO 110 K=1,L1
1177 DO 109 I=3,IDO,2
1178 IC = IDP2-I
1179 CH(I-1,K,J) = CC(I-1,2*J-1,K)+CC(IC-1,2*J-2,K)
1180 CH(I-1,K,JC) = CC(I-1,2*J-1,K)-CC(IC-1,2*J-2,K)
1181 CH(I,K,J) = CC(I,2*J-1,K)-CC(IC,2*J-2,K)
1182 CH(I,K,JC) = CC(I,2*J-1,K)+CC(IC,2*J-2,K)
1183 109 CONTINUE
1184 110 CONTINUE
1185 111 CONTINUE
1186 GO TO 116
1187 112 DO 115 J=2,IPPH
1188 JC = IPP2-J
1189 DO 114 I=3,IDO,2
1190 IC = IDP2-I
1191 DO 113 K=1,L1
1192 CH(I-1,K,J) = CC(I-1,2*J-1,K)+CC(IC-1,2*J-2,K)
1193 CH(I-1,K,JC) = CC(I-1,2*J-1,K)-CC(IC-1,2*J-2,K)
1194 CH(I,K,J) = CC(I,2*J-1,K)-CC(IC,2*J-2,K)
1195 CH(I,K,JC) = CC(I,2*J-1,K)+CC(IC,2*J-2,K)
1196 113 CONTINUE
1197 114 CONTINUE
1198 115 CONTINUE
1199 116 AR1 = 1.D0
1200 AI1 = 0.D0
1201 DO 120 L=2,IPPH
1202 LC = IPP2-L
1203 AR1H = DCP*AR1-DSP*AI1
1204 AI1 = DCP*AI1+DSP*AR1
1205 AR1 = AR1H
1206 DO 117 IK=1,IDL1
1207 C2(IK,L) = CH2(IK,1)+AR1*CH2(IK,2)
1208 C2(IK,LC) = AI1*CH2(IK,IP)
1209 117 CONTINUE
1210 DC2 = AR1
1211 DS2 = AI1
1212 AR2 = AR1
1213 AI2 = AI1
1214 DO 119 J=3,IPPH
1215 JC = IPP2-J
1216 AR2H = DC2*AR2-DS2*AI2
1217 AI2 = DC2*AI2+DS2*AR2
1218 AR2 = AR2H
1219 DO 118 IK=1,IDL1
1220 C2(IK,L) = C2(IK,L)+AR2*CH2(IK,J)
1221 C2(IK,LC) = C2(IK,LC)+AI2*CH2(IK,JC)
1222 118 CONTINUE
1223 119 CONTINUE
1224 120 CONTINUE
1225 DO 122 J=2,IPPH
1226 DO 121 IK=1,IDL1
1227 CH2(IK,1) = CH2(IK,1)+CH2(IK,J)
1228 121 CONTINUE
1229 122 CONTINUE
1230 DO 124 J=2,IPPH
1231 JC = IPP2-J
1232 DO 123 K=1,L1
1233 CH(1,K,J) = C1(1,K,J)-C1(1,K,JC)
1234 CH(1,K,JC) = C1(1,K,J)+C1(1,K,JC)
1235 123 CONTINUE
1236 124 CONTINUE
1237 IF (IDO .EQ. 1) GO TO 132
1238 IF (NBD .LT. L1) GO TO 128
1239 DO 127 J=2,IPPH
1240 JC = IPP2-J
1241 DO 126 K=1,L1
1242 DO 125 I=3,IDO,2
1243 CH(I-1,K,J) = C1(I-1,K,J)-C1(I,K,JC)
1244 CH(I-1,K,JC) = C1(I-1,K,J)+C1(I,K,JC)
1245 CH(I,K,J) = C1(I,K,J)+C1(I-1,K,JC)
1246 CH(I,K,JC) = C1(I,K,J)-C1(I-1,K,JC)
1247 125 CONTINUE
1248 126 CONTINUE
1249 127 CONTINUE
1250 GO TO 132
1251 128 DO 131 J=2,IPPH
1252 JC = IPP2-J
1253 DO 130 I=3,IDO,2
1254 DO 129 K=1,L1
1255 CH(I-1,K,J) = C1(I-1,K,J)-C1(I,K,JC)
1256 CH(I-1,K,JC) = C1(I-1,K,J)+C1(I,K,JC)
1257 CH(I,K,J) = C1(I,K,J)+C1(I-1,K,JC)
1258 CH(I,K,JC) = C1(I,K,J)-C1(I-1,K,JC)
1259 129 CONTINUE
1260 130 CONTINUE
1261 131 CONTINUE
1262 132 CONTINUE
1263 IF (IDO .EQ. 1) RETURN
1264 DO 133 IK=1,IDL1
1265 C2(IK,1) = CH2(IK,1)
1266 133 CONTINUE
1267 DO 135 J=2,IP
1268 DO 134 K=1,L1
1269 C1(1,K,J) = CH(1,K,J)
1270 134 CONTINUE
1271 135 CONTINUE
1272 IF (NBD .GT. L1) GO TO 139
1273 IS = -IDO
1274 DO 138 J=2,IP
1275 IS = IS+IDO
1276 IDIJ = IS
1277 DO 137 I=3,IDO,2
1278 IDIJ = IDIJ+2
1279 DO 136 K=1,L1
1280 C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J)
1281 C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J)
1282 136 CONTINUE
1283 137 CONTINUE
1284 138 CONTINUE
1285 GO TO 143
1286 139 IS = -IDO
1287 DO 142 J=2,IP
1288 IS = IS+IDO
1289 DO 141 K=1,L1
1290 IDIJ = IS
1291 DO 140 I=3,IDO,2
1292 IDIJ = IDIJ+2
1293 C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J)
1294 C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J)
1295 140 CONTINUE
1296 141 CONTINUE
1297 142 CONTINUE
1298 143 RETURN
1299 END
1300 SUBROUTINE R8FFTB1 (N,C,CH,WA,IFAC)
1301 implicit real*8 (A-H,O-Z)
1302 IMPLICIT INTEGER (I-N)
1303 DIMENSION CH(1) ,C(1) ,WA(1) ,IFAC(1)
1304 NF = IFAC(2)
1305 NA = 0
1306 L1 = 1
1307 IW = 1
1308 DO 116 K1=1,NF
1309 IP = IFAC(K1+2)
1310 L2 = IP*L1
1311 IDO = N/L2
1312 IDL1 = IDO*L1
1313 IF (IP .NE. 4) GO TO 103
1314 IX2 = IW+IDO
1315 IX3 = IX2+IDO
1316 IF (NA .NE. 0) GO TO 101
1317 CALL R8ADB4 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3))
1318 GO TO 102
1319 101 CALL R8ADB4 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3))
1320 102 NA = 1-NA
1321 GO TO 115
1322 103 IF (IP .NE. 2) GO TO 106
1323 IF (NA .NE. 0) GO TO 104
1324 CALL R8ADB2 (IDO,L1,C,CH,WA(IW))
1325 GO TO 105
1326 104 CALL R8ADB2 (IDO,L1,CH,C,WA(IW))
1327 105 NA = 1-NA
1328 GO TO 115
1329 106 IF (IP .NE. 3) GO TO 109
1330 IX2 = IW+IDO
1331 IF (NA .NE. 0) GO TO 107
1332 CALL R8ADB3 (IDO,L1,C,CH,WA(IW),WA(IX2))
1333 GO TO 108
1334 107 CALL R8ADB3 (IDO,L1,CH,C,WA(IW),WA(IX2))
1335 108 NA = 1-NA
1336 GO TO 115
1337 109 IF (IP .NE. 5) GO TO 112
1338 IX2 = IW+IDO
1339 IX3 = IX2+IDO
1340 IX4 = IX3+IDO
1341 IF (NA .NE. 0) GO TO 110
1342 CALL R8ADB5 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4))
1343 GO TO 111
1344 110 CALL R8ADB5 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4))
1345 111 NA = 1-NA
1346 GO TO 115
1347 112 IF (NA .NE. 0) GO TO 113
1348 CALL R8ADBG (IDO,IP,L1,IDL1,C,C,C,CH,CH,WA(IW))
1349 GO TO 114
1350 113 CALL R8ADBG (IDO,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW))
1351 114 IF (IDO .EQ. 1) NA = 1-NA
1352 115 L1 = L2
1353 IW = IW+(IP-1)*IDO
1354 116 CONTINUE
1355 IF (NA .EQ. 0) RETURN
1356 DO 117 I=1,N
1357 C(I) = CH(I)
1358 117 CONTINUE
1359 RETURN
1360 END
1361 C SUBROUTINE R8FFTF (N,R,WSAVE)
1362 C implicit real*8 (A-H,O-Z)
1363 C DIMENSION R(1) ,WSAVE(1)
1364 C IF (N .EQ. 1) RETURN
1365 C CALL R8FFTF1 (N,R,WSAVE,WSAVE(N+1),WSAVE(2*N+1))
1366 C RETURN
1367 C END
1368 SUBROUTINE R8ADF2 (IDO,L1,CC,CH,WA1)
1369 implicit real*8 (A-H,O-Z)
1370 IMPLICIT INTEGER (I-N)
1371 DIMENSION CH(IDO,2,L1) ,CC(IDO,L1,2) ,
1372 1 WA1(1)
1373 DO 101 K=1,L1
1374 CH(1,1,K) = CC(1,K,1)+CC(1,K,2)
1375 CH(IDO,2,K) = CC(1,K,1)-CC(1,K,2)
1376 101 CONTINUE
1377 IF (IDO-2) 107,105,102
1378 102 IDP2 = IDO+2
1379 DO 104 K=1,L1
1380 DO 103 I=3,IDO,2
1381 IC = IDP2-I
1382 TR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
1383 TI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
1384 CH(I,1,K) = CC(I,K,1)+TI2
1385 CH(IC,2,K) = TI2-CC(I,K,1)
1386 CH(I-1,1,K) = CC(I-1,K,1)+TR2
1387 CH(IC-1,2,K) = CC(I-1,K,1)-TR2
1388 103 CONTINUE
1389 104 CONTINUE
1390 IF (MOD(IDO,2) .EQ. 1) RETURN
1391 105 DO 106 K=1,L1
1392 CH(1,2,K) = -CC(IDO,K,2)
1393 CH(IDO,1,K) = CC(IDO,K,1)
1394 106 CONTINUE
1395 107 RETURN
1396 END
1397 SUBROUTINE R8ADF3 (IDO,L1,CC,CH,WA1,WA2)
1398 implicit real*8 (A-H,O-Z)
1399 IMPLICIT INTEGER (I-N)
1400 DIMENSION CH(IDO,3,L1) ,CC(IDO,L1,3) ,
1401 1 WA1(1) ,WA2(1)
1402 DATA TAUR,TAUI /-.5D0,.866025403784439D0/
1403 DO 101 K=1,L1
1404 CR2 = CC(1,K,2)+CC(1,K,3)
1405 CH(1,1,K) = CC(1,K,1)+CR2
1406 CH(1,3,K) = TAUI*(CC(1,K,3)-CC(1,K,2))
1407 CH(IDO,2,K) = CC(1,K,1)+TAUR*CR2
1408 101 CONTINUE
1409 IF (IDO .EQ. 1) RETURN
1410 IDP2 = IDO+2
1411 DO 103 K=1,L1
1412 DO 102 I=3,IDO,2
1413 IC = IDP2-I
1414 DR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
1415 DI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
1416 DR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3)
1417 DI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3)
1418 CR2 = DR2+DR3
1419 CI2 = DI2+DI3
1420 CH(I-1,1,K) = CC(I-1,K,1)+CR2
1421 CH(I,1,K) = CC(I,K,1)+CI2
1422 TR2 = CC(I-1,K,1)+TAUR*CR2
1423 TI2 = CC(I,K,1)+TAUR*CI2
1424 TR3 = TAUI*(DI2-DI3)
1425 TI3 = TAUI*(DR3-DR2)
1426 CH(I-1,3,K) = TR2+TR3
1427 CH(IC-1,2,K) = TR2-TR3
1428 CH(I,3,K) = TI2+TI3
1429 CH(IC,2,K) = TI3-TI2
1430 102 CONTINUE
1431 103 CONTINUE
1432 RETURN
1433 END
1434 SUBROUTINE R8ADF4 (IDO,L1,CC,CH,WA1,WA2,WA3)
1435 implicit real*8 (A-H,O-Z)
1436 IMPLICIT INTEGER (I-N)
1437 DIMENSION CC(IDO,L1,4) ,CH(IDO,4,L1) ,
1438 1 WA1(1) ,WA2(1) ,WA3(1)
1439 DATA HSQT2 /.7071067811865475D0/
1440 DO 101 K=1,L1
1441 TR1 = CC(1,K,2)+CC(1,K,4)
1442 TR2 = CC(1,K,1)+CC(1,K,3)
1443 CH(1,1,K) = TR1+TR2
1444 CH(IDO,4,K) = TR2-TR1
1445 CH(IDO,2,K) = CC(1,K,1)-CC(1,K,3)
1446 CH(1,3,K) = CC(1,K,4)-CC(1,K,2)
1447 101 CONTINUE
1448 IF (IDO-2) 107,105,102
1449 102 IDP2 = IDO+2
1450 DO 104 K=1,L1
1451 DO 103 I=3,IDO,2
1452 IC = IDP2-I
1453 CR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
1454 CI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
1455 CR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3)
1456 CI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3)
1457 CR4 = WA3(I-2)*CC(I-1,K,4)+WA3(I-1)*CC(I,K,4)
1458 CI4 = WA3(I-2)*CC(I,K,4)-WA3(I-1)*CC(I-1,K,4)
1459 TR1 = CR2+CR4
1460 TR4 = CR4-CR2
1461 TI1 = CI2+CI4
1462 TI4 = CI2-CI4
1463 TI2 = CC(I,K,1)+CI3
1464 TI3 = CC(I,K,1)-CI3
1465 TR2 = CC(I-1,K,1)+CR3
1466 TR3 = CC(I-1,K,1)-CR3
1467 CH(I-1,1,K) = TR1+TR2
1468 CH(IC-1,4,K) = TR2-TR1
1469 CH(I,1,K) = TI1+TI2
1470 CH(IC,4,K) = TI1-TI2
1471 CH(I-1,3,K) = TI4+TR3
1472 CH(IC-1,2,K) = TR3-TI4
1473 CH(I,3,K) = TR4+TI3
1474 CH(IC,2,K) = TR4-TI3
1475 103 CONTINUE
1476 104 CONTINUE
1477 IF (MOD(IDO,2) .EQ. 1) RETURN
1478 105 CONTINUE
1479 DO 106 K=1,L1
1480 TI1 = -HSQT2*(CC(IDO,K,2)+CC(IDO,K,4))
1481 TR1 = HSQT2*(CC(IDO,K,2)-CC(IDO,K,4))
1482 CH(IDO,1,K) = TR1+CC(IDO,K,1)
1483 CH(IDO,3,K) = CC(IDO,K,1)-TR1
1484 CH(1,2,K) = TI1-CC(IDO,K,3)
1485 CH(1,4,K) = TI1+CC(IDO,K,3)
1486 106 CONTINUE
1487 107 RETURN
1488 END
1489 SUBROUTINE R8ADF5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
1490 implicit real*8 (A-H,O-Z)
1491 IMPLICIT INTEGER (I-N)
1492 DIMENSION CC(IDO,L1,5) ,CH(IDO,5,L1) ,
1493 1 WA1(1) ,WA2(1) ,WA3(1) ,WA4(1)
1494 DATA TR11,TI11,TR12,TI12 /.309016994374947D0,.951056516295154D0,
1495 1-.809016994374947D0,.587785252292473D0/
1496 DO 101 K=1,L1
1497 CR2 = CC(1,K,5)+CC(1,K,2)
1498 CI5 = CC(1,K,5)-CC(1,K,2)
1499 CR3 = CC(1,K,4)+CC(1,K,3)
1500 CI4 = CC(1,K,4)-CC(1,K,3)
1501 CH(1,1,K) = CC(1,K,1)+CR2+CR3
1502 CH(IDO,2,K) = CC(1,K,1)+TR11*CR2+TR12*CR3
1503 CH(1,3,K) = TI11*CI5+TI12*CI4
1504 CH(IDO,4,K) = CC(1,K,1)+TR12*CR2+TR11*CR3
1505 CH(1,5,K) = TI12*CI5-TI11*CI4
1506 101 CONTINUE
1507 IF (IDO .EQ. 1) RETURN
1508 IDP2 = IDO+2
1509 DO 103 K=1,L1
1510 DO 102 I=3,IDO,2
1511 IC = IDP2-I
1512 DR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
1513 DI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
1514 DR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3)
1515 DI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3)
1516 DR4 = WA3(I-2)*CC(I-1,K,4)+WA3(I-1)*CC(I,K,4)
1517 DI4 = WA3(I-2)*CC(I,K,4)-WA3(I-1)*CC(I-1,K,4)
1518 DR5 = WA4(I-2)*CC(I-1,K,5)+WA4(I-1)*CC(I,K,5)
1519 DI5 = WA4(I-2)*CC(I,K,5)-WA4(I-1)*CC(I-1,K,5)
1520 CR2 = DR2+DR5
1521 CI5 = DR5-DR2
1522 CR5 = DI2-DI5
1523 CI2 = DI2+DI5
1524 CR3 = DR3+DR4
1525 CI4 = DR4-DR3
1526 CR4 = DI3-DI4
1527 CI3 = DI3+DI4
1528 CH(I-1,1,K) = CC(I-1,K,1)+CR2+CR3
1529 CH(I,1,K) = CC(I,K,1)+CI2+CI3
1530 TR2 = CC(I-1,K,1)+TR11*CR2+TR12*CR3
1531 TI2 = CC(I,K,1)+TR11*CI2+TR12*CI3
1532 TR3 = CC(I-1,K,1)+TR12*CR2+TR11*CR3
1533 TI3 = CC(I,K,1)+TR12*CI2+TR11*CI3
1534 TR5 = TI11*CR5+TI12*CR4
1535 TI5 = TI11*CI5+TI12*CI4
1536 TR4 = TI12*CR5-TI11*CR4
1537 TI4 = TI12*CI5-TI11*CI4
1538 CH(I-1,3,K) = TR2+TR5
1539 CH(IC-1,2,K) = TR2-TR5
1540 CH(I,3,K) = TI2+TI5
1541 CH(IC,2,K) = TI5-TI2
1542 CH(I-1,5,K) = TR3+TR4
1543 CH(IC-1,4,K) = TR3-TR4
1544 CH(I,5,K) = TI3+TI4
1545 CH(IC,4,K) = TI4-TI3
1546 102 CONTINUE
1547 103 CONTINUE
1548 RETURN
1549 END
1550 SUBROUTINE R8ADFG (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
1551 implicit real*8 (A-H,O-Z)
1552 IMPLICIT INTEGER (I-N)
1553 DIMENSION CH(IDO,L1,IP) ,CC(IDO,IP,L1) ,
1554 1 C1(IDO,L1,IP) ,C2(IDL1,IP),
1555 2 CH2(IDL1,IP) ,WA(1)
1556 DATA TPI/6.28318530717959D0/
1557 ARG = TPI/FLOAT(IP)
1558 DCP = COS(ARG)
1559 DSP = SIN(ARG)
1560 IPPH = (IP+1)/2
1561 IPP2 = IP+2
1562 IDP2 = IDO+2
1563 NBD = (IDO-1)/2
1564 IF (IDO .EQ. 1) GO TO 119
1565 DO 101 IK=1,IDL1
1566 CH2(IK,1) = C2(IK,1)
1567 101 CONTINUE
1568 DO 103 J=2,IP
1569 DO 102 K=1,L1
1570 CH(1,K,J) = C1(1,K,J)
1571 102 CONTINUE
1572 103 CONTINUE
1573 IF (NBD .GT. L1) GO TO 107
1574 IS = -IDO
1575 DO 106 J=2,IP
1576 IS = IS+IDO
1577 IDIJ = IS
1578 DO 105 I=3,IDO,2
1579 IDIJ = IDIJ+2
1580 DO 104 K=1,L1
1581 CH(I-1,K,J) = WA(IDIJ-1)*C1(I-1,K,J)+WA(IDIJ)*C1(I,K,J)
1582 CH(I,K,J) = WA(IDIJ-1)*C1(I,K,J)-WA(IDIJ)*C1(I-1,K,J)
1583 104 CONTINUE
1584 105 CONTINUE
1585 106 CONTINUE
1586 GO TO 111
1587 107 IS = -IDO
1588 DO 110 J=2,IP
1589 IS = IS+IDO
1590 DO 109 K=1,L1
1591 IDIJ = IS
1592 DO 108 I=3,IDO,2
1593 IDIJ = IDIJ+2
1594 CH(I-1,K,J) = WA(IDIJ-1)*C1(I-1,K,J)+WA(IDIJ)*C1(I,K,J)
1595 CH(I,K,J) = WA(IDIJ-1)*C1(I,K,J)-WA(IDIJ)*C1(I-1,K,J)
1596 108 CONTINUE
1597 109 CONTINUE
1598 110 CONTINUE
1599 111 IF (NBD .LT. L1) GO TO 115
1600 DO 114 J=2,IPPH
1601 JC = IPP2-J
1602 DO 113 K=1,L1
1603 DO 112 I=3,IDO,2
1604 C1(I-1,K,J) = CH(I-1,K,J)+CH(I-1,K,JC)
1605 C1(I-1,K,JC) = CH(I,K,J)-CH(I,K,JC)
1606 C1(I,K,J) = CH(I,K,J)+CH(I,K,JC)
1607 C1(I,K,JC) = CH(I-1,K,JC)-CH(I-1,K,J)
1608 112 CONTINUE
1609 113 CONTINUE
1610 114 CONTINUE
1611 GO TO 121
1612 115 DO 118 J=2,IPPH
1613 JC = IPP2-J
1614 DO 117 I=3,IDO,2
1615 DO 116 K=1,L1
1616 C1(I-1,K,J) = CH(I-1,K,J)+CH(I-1,K,JC)
1617 C1(I-1,K,JC) = CH(I,K,J)-CH(I,K,JC)
1618 C1(I,K,J) = CH(I,K,J)+CH(I,K,JC)
1619 C1(I,K,JC) = CH(I-1,K,JC)-CH(I-1,K,J)
1620 116 CONTINUE
1621 117 CONTINUE
1622 118 CONTINUE
1623 GO TO 121
1624 119 DO 120 IK=1,IDL1
1625 C2(IK,1) = CH2(IK,1)
1626 120 CONTINUE
1627 121 DO 123 J=2,IPPH
1628 JC = IPP2-J
1629 DO 122 K=1,L1
1630 C1(1,K,J) = CH(1,K,J)+CH(1,K,JC)
1631 C1(1,K,JC) = CH(1,K,JC)-CH(1,K,J)
1632 122 CONTINUE
1633 123 CONTINUE
1634 C
1635 AR1 = 1.D0
1636 AI1 = 0.D0
1637 DO 127 L=2,IPPH
1638 LC = IPP2-L
1639 AR1H = DCP*AR1-DSP*AI1
1640 AI1 = DCP*AI1+DSP*AR1
1641 AR1 = AR1H
1642 DO 124 IK=1,IDL1
1643 CH2(IK,L) = C2(IK,1)+AR1*C2(IK,2)
1644 CH2(IK,LC) = AI1*C2(IK,IP)
1645 124 CONTINUE
1646 DC2 = AR1
1647 DS2 = AI1
1648 AR2 = AR1
1649 AI2 = AI1
1650 DO 126 J=3,IPPH
1651 JC = IPP2-J
1652 AR2H = DC2*AR2-DS2*AI2
1653 AI2 = DC2*AI2+DS2*AR2
1654 AR2 = AR2H
1655 DO 125 IK=1,IDL1
1656 CH2(IK,L) = CH2(IK,L)+AR2*C2(IK,J)
1657 CH2(IK,LC) = CH2(IK,LC)+AI2*C2(IK,JC)
1658 125 CONTINUE
1659 126 CONTINUE
1660 127 CONTINUE
1661 DO 129 J=2,IPPH
1662 DO 128 IK=1,IDL1
1663 CH2(IK,1) = CH2(IK,1)+C2(IK,J)
1664 128 CONTINUE
1665 129 CONTINUE
1666 C
1667 IF (IDO .LT. L1) GO TO 132
1668 DO 131 K=1,L1
1669 DO 130 I=1,IDO
1670 CC(I,1,K) = CH(I,K,1)
1671 130 CONTINUE
1672 131 CONTINUE
1673 GO TO 135
1674 132 DO 134 I=1,IDO
1675 DO 133 K=1,L1
1676 CC(I,1,K) = CH(I,K,1)
1677 133 CONTINUE
1678 134 CONTINUE
1679 135 DO 137 J=2,IPPH
1680 JC = IPP2-J
1681 J2 = J+J
1682 DO 136 K=1,L1
1683 CC(IDO,J2-2,K) = CH(1,K,J)
1684 CC(1,J2-1,K) = CH(1,K,JC)
1685 136 CONTINUE
1686 137 CONTINUE
1687 IF (IDO .EQ. 1) RETURN
1688 IF (NBD .LT. L1) GO TO 141
1689 DO 140 J=2,IPPH
1690 JC = IPP2-J
1691 J2 = J+J
1692 DO 139 K=1,L1
1693 DO 138 I=3,IDO,2
1694 IC = IDP2-I
1695 CC(I-1,J2-1,K) = CH(I-1,K,J)+CH(I-1,K,JC)
1696 CC(IC-1,J2-2,K) = CH(I-1,K,J)-CH(I-1,K,JC)
1697 CC(I,J2-1,K) = CH(I,K,J)+CH(I,K,JC)
1698 CC(IC,J2-2,K) = CH(I,K,JC)-CH(I,K,J)
1699 138 CONTINUE
1700 139 CONTINUE
1701 140 CONTINUE
1702 RETURN
1703 141 DO 144 J=2,IPPH
1704 JC = IPP2-J
1705 J2 = J+J
1706 DO 143 I=3,IDO,2
1707 IC = IDP2-I
1708 DO 142 K=1,L1
1709 CC(I-1,J2-1,K) = CH(I-1,K,J)+CH(I-1,K,JC)
1710 CC(IC-1,J2-2,K) = CH(I-1,K,J)-CH(I-1,K,JC)
1711 CC(I,J2-1,K) = CH(I,K,J)+CH(I,K,JC)
1712 CC(IC,J2-2,K) = CH(I,K,JC)-CH(I,K,J)
1713 142 CONTINUE
1714 143 CONTINUE
1715 144 CONTINUE
1716 RETURN
1717 END
1718 SUBROUTINE R8FFTF1 (N,C,CH,WA,IFAC)
1719 implicit real*8 (A-H,O-Z)
1720 IMPLICIT INTEGER (I-N)
1721 DIMENSION CH(1) ,C(1) ,WA(1) ,IFAC(1)
1722 NF = IFAC(2)
1723 NA = 1
1724 L2 = N
1725 IW = N
1726 DO 111 K1=1,NF
1727 KH = NF-K1
1728 IP = IFAC(KH+3)
1729 L1 = L2/IP
1730 IDO = N/L2
1731 IDL1 = IDO*L1
1732 IW = IW-(IP-1)*IDO
1733 NA = 1-NA
1734 IF (IP .NE. 4) GO TO 102
1735 IX2 = IW+IDO
1736 IX3 = IX2+IDO
1737 IF (NA .NE. 0) GO TO 101
1738 CALL R8ADF4 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3))
1739 GO TO 110
1740 101 CALL R8ADF4 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3))
1741 GO TO 110
1742 102 IF (IP .NE. 2) GO TO 104
1743 IF (NA .NE. 0) GO TO 103
1744 CALL R8ADF2 (IDO,L1,C,CH,WA(IW))
1745 GO TO 110
1746 103 CALL R8ADF2 (IDO,L1,CH,C,WA(IW))
1747 GO TO 110
1748 104 IF (IP .NE. 3) GO TO 106
1749 IX2 = IW+IDO
1750 IF (NA .NE. 0) GO TO 105
1751 CALL R8ADF3 (IDO,L1,C,CH,WA(IW),WA(IX2))
1752 GO TO 110
1753 105 CALL R8ADF3 (IDO,L1,CH,C,WA(IW),WA(IX2))
1754 GO TO 110
1755 106 IF (IP .NE. 5) GO TO 108
1756 IX2 = IW+IDO
1757 IX3 = IX2+IDO
1758 IX4 = IX3+IDO
1759 IF (NA .NE. 0) GO TO 107
1760 CALL R8ADF5 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4))
1761 GO TO 110
1762 107 CALL R8ADF5 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4))
1763 GO TO 110
1764 108 IF (IDO .EQ. 1) NA = 1-NA
1765 IF (NA .NE. 0) GO TO 109
1766 CALL R8ADFG (IDO,IP,L1,IDL1,C,C,C,CH,CH,WA(IW))
1767 NA = 1
1768 GO TO 110
1769 109 CALL R8ADFG (IDO,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW))
1770 NA = 0
1771 110 L2 = L1
1772 111 CONTINUE
1773 IF (NA .EQ. 1) RETURN
1774 DO 112 I=1,N
1775 C(I) = CH(I)
1776 112 CONTINUE
1777 RETURN
1778 END
1779 C SUBROUTINE R8FFTI (N,WSAVE)
1780 C implicit real*8 (A-H,O-Z)
1781 C DIMENSION WSAVE(1)
1782 C IF (N .EQ. 1) RETURN
1783 C CALL R8FFTI1 (N,WSAVE(N+1),WSAVE(2*N+1))
1784 C RETURN
1785 C END
1786 SUBROUTINE R8FFTI1 (N,WA,IFAC)
1787 implicit real*8 (A-H,O-Z)
1788 IMPLICIT INTEGER (I-N)
1789 DIMENSION WA(1) ,IFAC(1) ,NTRYH(4)
1790 DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/4,2,3,5/
1791 NL = N
1792 NF = 0
1793 J = 0
1794 101 J = J+1
1795 IF (J-4) 102,102,103
1796 102 NTRY = NTRYH(J)
1797 GO TO 104
1798 103 NTRY = NTRY+2
1799 104 NQ = NL/NTRY
1800 NR = NL-NTRY*NQ
1801 IF (NR) 101,105,101
1802 105 NF = NF+1
1803 IFAC(NF+2) = NTRY
1804 NL = NQ
1805 IF (NTRY .NE. 2) GO TO 107
1806 IF (NF .EQ. 1) GO TO 107
1807 DO 106 I=2,NF
1808 IB = NF-I+2
1809 IFAC(IB+2) = IFAC(IB+1)
1810 106 CONTINUE
1811 IFAC(3) = 2
1812 107 IF (NL .NE. 1) GO TO 104
1813 IFAC(1) = N
1814 IFAC(2) = NF
1815 TPI = 6.28318530717959D0
1816 ARGH = TPI/FLOAT(N)
1817 IS = 0
1818 NFM1 = NF-1
1819 L1 = 1
1820 IF (NFM1 .EQ. 0) RETURN
1821 DO 110 K1=1,NFM1
1822 IP = IFAC(K1+2)
1823 LD = 0
1824 L2 = L1*IP
1825 IDO = N/L2
1826 IPM = IP-1
1827 DO 109 J=1,IPM
1828 LD = LD+L1
1829 I = IS
1830 ARGLD = FLOAT(LD)*ARGH
1831 FI = 0.D0
1832 DO 108 II=3,IDO,2
1833 I = I+2
1834 FI = FI+1.D0
1835 ARG = FI*ARGLD
1836 WA(I-1) = COS(ARG)
1837 WA(I) = SIN(ARG)
1838 108 CONTINUE
1839 IS = IS+IDO
1840 109 CONTINUE
1841 L1 = L2
1842 110 CONTINUE
1843 RETURN
1844 END
  ViewVC Help
Powered by ViewVC 1.1.22