/[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.6 - (show annotations) (download)
Thu Oct 9 04:19:20 2003 UTC (20 years, 6 months ago) by edhill
Branch: MAIN
CVS Tags: checkpoint51k_post, checkpoint52d_pre, checkpoint51o_pre, checkpoint51l_post, checkpoint52, checkpoint51t_post, checkpoint51n_post, checkpoint51s_post, checkpoint51n_pre, checkpoint52b_pre, checkpoint51l_pre, checkpoint51q_post, checkpoint52b_post, checkpoint52c_post, checkpoint51r_post, checkpoint51i_post, checkpoint52d_post, checkpoint52a_pre, branch-netcdf, checkpoint51o_post, checkpoint52a_post, ecco_c52_e35, checkpoint51m_post, checkpoint51p_post, checkpoint51u_post
Branch point for: branch-nonh, tg2-branch, netcdf-sm0, checkpoint51n_branch
Changes since 1.5: +9 -1 lines 
 o first check-in for the "branch-genmake2" merge
 o verification suite as run on shelley (gcc 3.2.2):

Wed Oct  8 23:42:29 EDT 2003
                T           S           U           V
G D M    c        m  s        m  s        m  s        m  s
E p a R  g  m  m  e  .  m  m  e  .  m  m  e  .  m  m  e  .
N n k u  2  i  a  a  d  i  a  a  d  i  a  a  d  i  a  a  d
2 d e n  d  n  x  n  .  n  x  n  .  n  x  n  .  n  x  n  .

OPTFILE=NONE

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