File indexing completed on 2018-03-02 18:44:40 UTC
view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
6d54cf9ca1 Ed H*0001 #include "ZONAL_FILT_OPTIONS.h"
0002
42c525bfb4 Alis*0003
0004
0005
0006
0007
0008
0009 SUBROUTINE RADB2 (IDO,L1,CC,CH,WA1)
aea29c8517 Alis*0010
9bf05f3f81 Alis*0011 IMPLICIT REAL*8 (A-H,O-Z)
aea29c8517 Alis*0012 IMPLICIT INTEGER (I-N)
0013
6d54cf9ca1 Ed H*0014 #ifdef ALLOW_ZONAL_FILT
0015
42c525bfb4 Alis*0016 DIMENSION CC(IDO,2,L1) ,CH(IDO,L1,2) ,
3efde2d9d6 Jean*0017 1 WA1(*)
42c525bfb4 Alis*0018 DO 101 K=1,L1
0019 CH(1,K,1) = CC(1,1,K)+CC(IDO,2,K)
0020 CH(1,K,2) = CC(1,1,K)-CC(IDO,2,K)
0021 101 CONTINUE
0022 IF (IDO-2) 107,105,102
0023 102 IDP2 = IDO+2
0024 DO 104 K=1,L1
0025 DO 103 I=3,IDO,2
0026 IC = IDP2-I
0027 CH(I-1,K,1) = CC(I-1,1,K)+CC(IC-1,2,K)
0028 TR2 = CC(I-1,1,K)-CC(IC-1,2,K)
0029 CH(I,K,1) = CC(I,1,K)-CC(IC,2,K)
0030 TI2 = CC(I,1,K)+CC(IC,2,K)
0031 CH(I-1,K,2) = WA1(I-2)*TR2-WA1(I-1)*TI2
0032 CH(I,K,2) = WA1(I-2)*TI2+WA1(I-1)*TR2
0033 103 CONTINUE
0034 104 CONTINUE
0035 IF (MOD(IDO,2) .EQ. 1) RETURN
0036 105 DO 106 K=1,L1
0037 CH(IDO,K,1) = CC(IDO,1,K)+CC(IDO,1,K)
0038 CH(IDO,K,2) = -(CC(1,2,K)+CC(1,2,K))
0039 106 CONTINUE
0040 107 RETURN
0041 END
0042 SUBROUTINE RADB3 (IDO,L1,CC,CH,WA1,WA2)
aea29c8517 Alis*0043
9bf05f3f81 Alis*0044 IMPLICIT REAL*8 (A-H,O-Z)
aea29c8517 Alis*0045 IMPLICIT INTEGER (I-N)
0046
42c525bfb4 Alis*0047 DIMENSION CC(IDO,3,L1) ,CH(IDO,L1,3) ,
3efde2d9d6 Jean*0048 1 WA1(*) ,WA2(*)
9bf05f3f81 Alis*0049 DATA TAUR,TAUI /-.5D0,.866025403784439D0/
42c525bfb4 Alis*0050 DO 101 K=1,L1
0051 TR2 = CC(IDO,2,K)+CC(IDO,2,K)
0052 CR2 = CC(1,1,K)+TAUR*TR2
0053 CH(1,K,1) = CC(1,1,K)+TR2
0054 CI3 = TAUI*(CC(1,3,K)+CC(1,3,K))
0055 CH(1,K,2) = CR2-CI3
0056 CH(1,K,3) = CR2+CI3
0057 101 CONTINUE
0058 IF (IDO .EQ. 1) RETURN
0059 IDP2 = IDO+2
0060 DO 103 K=1,L1
0061 DO 102 I=3,IDO,2
0062 IC = IDP2-I
0063 TR2 = CC(I-1,3,K)+CC(IC-1,2,K)
0064 CR2 = CC(I-1,1,K)+TAUR*TR2
0065 CH(I-1,K,1) = CC(I-1,1,K)+TR2
0066 TI2 = CC(I,3,K)-CC(IC,2,K)
0067 CI2 = CC(I,1,K)+TAUR*TI2
0068 CH(I,K,1) = CC(I,1,K)+TI2
0069 CR3 = TAUI*(CC(I-1,3,K)-CC(IC-1,2,K))
0070 CI3 = TAUI*(CC(I,3,K)+CC(IC,2,K))
0071 DR2 = CR2-CI3
0072 DR3 = CR2+CI3
0073 DI2 = CI2+CR3
0074 DI3 = CI2-CR3
0075 CH(I-1,K,2) = WA1(I-2)*DR2-WA1(I-1)*DI2
0076 CH(I,K,2) = WA1(I-2)*DI2+WA1(I-1)*DR2
0077 CH(I-1,K,3) = WA2(I-2)*DR3-WA2(I-1)*DI3
0078 CH(I,K,3) = WA2(I-2)*DI3+WA2(I-1)*DR3
0079 102 CONTINUE
0080 103 CONTINUE
0081 RETURN
0082 END
0083 SUBROUTINE RADB4 (IDO,L1,CC,CH,WA1,WA2,WA3)
aea29c8517 Alis*0084
9bf05f3f81 Alis*0085 IMPLICIT REAL*8 (A-H,O-Z)
aea29c8517 Alis*0086 IMPLICIT INTEGER (I-N)
0087
42c525bfb4 Alis*0088 DIMENSION CC(IDO,4,L1) ,CH(IDO,L1,4) ,
3efde2d9d6 Jean*0089 1 WA1(*) ,WA2(*) ,WA3(*)
9bf05f3f81 Alis*0090 DATA SQRT2 /1.414213562373095D0/
42c525bfb4 Alis*0091 DO 101 K=1,L1
0092 TR1 = CC(1,1,K)-CC(IDO,4,K)
0093 TR2 = CC(1,1,K)+CC(IDO,4,K)
0094 TR3 = CC(IDO,2,K)+CC(IDO,2,K)
0095 TR4 = CC(1,3,K)+CC(1,3,K)
0096 CH(1,K,1) = TR2+TR3
0097 CH(1,K,2) = TR1-TR4
0098 CH(1,K,3) = TR2-TR3
0099 CH(1,K,4) = TR1+TR4
0100 101 CONTINUE
0101 IF (IDO-2) 107,105,102
0102 102 IDP2 = IDO+2
0103 DO 104 K=1,L1
0104 DO 103 I=3,IDO,2
0105 IC = IDP2-I
0106 TI1 = CC(I,1,K)+CC(IC,4,K)
0107 TI2 = CC(I,1,K)-CC(IC,4,K)
0108 TI3 = CC(I,3,K)-CC(IC,2,K)
0109 TR4 = CC(I,3,K)+CC(IC,2,K)
0110 TR1 = CC(I-1,1,K)-CC(IC-1,4,K)
0111 TR2 = CC(I-1,1,K)+CC(IC-1,4,K)
0112 TI4 = CC(I-1,3,K)-CC(IC-1,2,K)
0113 TR3 = CC(I-1,3,K)+CC(IC-1,2,K)
0114 CH(I-1,K,1) = TR2+TR3
0115 CR3 = TR2-TR3
0116 CH(I,K,1) = TI2+TI3
0117 CI3 = TI2-TI3
0118 CR2 = TR1-TR4
0119 CR4 = TR1+TR4
0120 CI2 = TI1+TI4
0121 CI4 = TI1-TI4
0122 CH(I-1,K,2) = WA1(I-2)*CR2-WA1(I-1)*CI2
0123 CH(I,K,2) = WA1(I-2)*CI2+WA1(I-1)*CR2
0124 CH(I-1,K,3) = WA2(I-2)*CR3-WA2(I-1)*CI3
0125 CH(I,K,3) = WA2(I-2)*CI3+WA2(I-1)*CR3
0126 CH(I-1,K,4) = WA3(I-2)*CR4-WA3(I-1)*CI4
0127 CH(I,K,4) = WA3(I-2)*CI4+WA3(I-1)*CR4
0128 103 CONTINUE
0129 104 CONTINUE
0130 IF (MOD(IDO,2) .EQ. 1) RETURN
0131 105 CONTINUE
0132 DO 106 K=1,L1
0133 TI1 = CC(1,2,K)+CC(1,4,K)
0134 TI2 = CC(1,4,K)-CC(1,2,K)
0135 TR1 = CC(IDO,1,K)-CC(IDO,3,K)
0136 TR2 = CC(IDO,1,K)+CC(IDO,3,K)
0137 CH(IDO,K,1) = TR2+TR2
0138 CH(IDO,K,2) = SQRT2*(TR1-TI1)
0139 CH(IDO,K,3) = TI2+TI2
0140 CH(IDO,K,4) = -SQRT2*(TR1+TI1)
0141 106 CONTINUE
0142 107 RETURN
0143 END
0144 SUBROUTINE RADB5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
aea29c8517 Alis*0145
9bf05f3f81 Alis*0146 IMPLICIT REAL*8 (A-H,O-Z)
aea29c8517 Alis*0147 IMPLICIT INTEGER (I-N)
0148
42c525bfb4 Alis*0149 DIMENSION CC(IDO,5,L1) ,CH(IDO,L1,5) ,
3efde2d9d6 Jean*0150 1 WA1(*) ,WA2(*) ,WA3(*) ,WA4(*)
9bf05f3f81 Alis*0151 DATA TR11,TI11,TR12,TI12 /.309016994374947D0,.951056516295154D0,
0152 1-.809016994374947D0,.587785252292473D0/
42c525bfb4 Alis*0153 DO 101 K=1,L1
0154 TI5 = CC(1,3,K)+CC(1,3,K)
0155 TI4 = CC(1,5,K)+CC(1,5,K)
0156 TR2 = CC(IDO,2,K)+CC(IDO,2,K)
0157 TR3 = CC(IDO,4,K)+CC(IDO,4,K)
0158 CH(1,K,1) = CC(1,1,K)+TR2+TR3
0159 CR2 = CC(1,1,K)+TR11*TR2+TR12*TR3
0160 CR3 = CC(1,1,K)+TR12*TR2+TR11*TR3
0161 CI5 = TI11*TI5+TI12*TI4
0162 CI4 = TI12*TI5-TI11*TI4
0163 CH(1,K,2) = CR2-CI5
0164 CH(1,K,3) = CR3-CI4
0165 CH(1,K,4) = CR3+CI4
0166 CH(1,K,5) = CR2+CI5
0167 101 CONTINUE
0168 IF (IDO .EQ. 1) RETURN
0169 IDP2 = IDO+2
0170 DO 103 K=1,L1
0171 DO 102 I=3,IDO,2
0172 IC = IDP2-I
0173 TI5 = CC(I,3,K)+CC(IC,2,K)
0174 TI2 = CC(I,3,K)-CC(IC,2,K)
0175 TI4 = CC(I,5,K)+CC(IC,4,K)
0176 TI3 = CC(I,5,K)-CC(IC,4,K)
0177 TR5 = CC(I-1,3,K)-CC(IC-1,2,K)
0178 TR2 = CC(I-1,3,K)+CC(IC-1,2,K)
0179 TR4 = CC(I-1,5,K)-CC(IC-1,4,K)
0180 TR3 = CC(I-1,5,K)+CC(IC-1,4,K)
0181 CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3
0182 CH(I,K,1) = CC(I,1,K)+TI2+TI3
0183 CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3
0184 CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3
0185 CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3
0186 CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3
0187 CR5 = TI11*TR5+TI12*TR4
0188 CI5 = TI11*TI5+TI12*TI4
0189 CR4 = TI12*TR5-TI11*TR4
0190 CI4 = TI12*TI5-TI11*TI4
0191 DR3 = CR3-CI4
0192 DR4 = CR3+CI4
0193 DI3 = CI3+CR4
0194 DI4 = CI3-CR4
0195 DR5 = CR2+CI5
0196 DR2 = CR2-CI5
0197 DI5 = CI2-CR5
0198 DI2 = CI2+CR5
0199 CH(I-1,K,2) = WA1(I-2)*DR2-WA1(I-1)*DI2
0200 CH(I,K,2) = WA1(I-2)*DI2+WA1(I-1)*DR2
0201 CH(I-1,K,3) = WA2(I-2)*DR3-WA2(I-1)*DI3
0202 CH(I,K,3) = WA2(I-2)*DI3+WA2(I-1)*DR3
0203 CH(I-1,K,4) = WA3(I-2)*DR4-WA3(I-1)*DI4
0204 CH(I,K,4) = WA3(I-2)*DI4+WA3(I-1)*DR4
0205 CH(I-1,K,5) = WA4(I-2)*DR5-WA4(I-1)*DI5
0206 CH(I,K,5) = WA4(I-2)*DI5+WA4(I-1)*DR5
0207 102 CONTINUE
0208 103 CONTINUE
0209 RETURN
0210 END
0211 SUBROUTINE RADBG (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
aea29c8517 Alis*0212
9bf05f3f81 Alis*0213 IMPLICIT REAL*8 (A-H,O-Z)
aea29c8517 Alis*0214 IMPLICIT INTEGER (I-N)
0215
42c525bfb4 Alis*0216 DIMENSION CH(IDO,L1,IP) ,CC(IDO,IP,L1) ,
0217 1 C1(IDO,L1,IP) ,C2(IDL1,IP),
3efde2d9d6 Jean*0218 2 CH2(IDL1,IP) ,WA(*)
9bf05f3f81 Alis*0219 DATA TPI/6.28318530717959D0/
42c525bfb4 Alis*0220 ARG = TPI/FLOAT(IP)
0221 DCP = COS(ARG)
0222 DSP = SIN(ARG)
0223 IDP2 = IDO+2
0224 NBD = (IDO-1)/2
0225 IPP2 = IP+2
0226 IPPH = (IP+1)/2
0227 IF (IDO .LT. L1) GO TO 103
0228 DO 102 K=1,L1
0229 DO 101 I=1,IDO
0230 CH(I,K,1) = CC(I,1,K)
0231 101 CONTINUE
0232 102 CONTINUE
0233 GO TO 106
0234 103 DO 105 I=1,IDO
0235 DO 104 K=1,L1
0236 CH(I,K,1) = CC(I,1,K)
0237 104 CONTINUE
0238 105 CONTINUE
0239 106 DO 108 J=2,IPPH
0240 JC = IPP2-J
0241 J2 = J+J
0242 DO 107 K=1,L1
0243 CH(1,K,J) = CC(IDO,J2-2,K)+CC(IDO,J2-2,K)
0244 CH(1,K,JC) = CC(1,J2-1,K)+CC(1,J2-1,K)
0245 107 CONTINUE
0246 108 CONTINUE
0247 IF (IDO .EQ. 1) GO TO 116
0248 IF (NBD .LT. L1) GO TO 112
0249 DO 111 J=2,IPPH
0250 JC = IPP2-J
0251 DO 110 K=1,L1
0252 DO 109 I=3,IDO,2
0253 IC = IDP2-I
0254 CH(I-1,K,J) = CC(I-1,2*J-1,K)+CC(IC-1,2*J-2,K)
0255 CH(I-1,K,JC) = CC(I-1,2*J-1,K)-CC(IC-1,2*J-2,K)
0256 CH(I,K,J) = CC(I,2*J-1,K)-CC(IC,2*J-2,K)
0257 CH(I,K,JC) = CC(I,2*J-1,K)+CC(IC,2*J-2,K)
0258 109 CONTINUE
0259 110 CONTINUE
0260 111 CONTINUE
0261 GO TO 116
0262 112 DO 115 J=2,IPPH
0263 JC = IPP2-J
0264 DO 114 I=3,IDO,2
0265 IC = IDP2-I
0266 DO 113 K=1,L1
0267 CH(I-1,K,J) = CC(I-1,2*J-1,K)+CC(IC-1,2*J-2,K)
0268 CH(I-1,K,JC) = CC(I-1,2*J-1,K)-CC(IC-1,2*J-2,K)
0269 CH(I,K,J) = CC(I,2*J-1,K)-CC(IC,2*J-2,K)
0270 CH(I,K,JC) = CC(I,2*J-1,K)+CC(IC,2*J-2,K)
0271 113 CONTINUE
0272 114 CONTINUE
0273 115 CONTINUE
9bf05f3f81 Alis*0274 116 AR1 = 1.D0
0275 AI1 = 0.D0
42c525bfb4 Alis*0276 DO 120 L=2,IPPH
0277 LC = IPP2-L
0278 AR1H = DCP*AR1-DSP*AI1
0279 AI1 = DCP*AI1+DSP*AR1
0280 AR1 = AR1H
0281 DO 117 IK=1,IDL1
0282 C2(IK,L) = CH2(IK,1)+AR1*CH2(IK,2)
0283 C2(IK,LC) = AI1*CH2(IK,IP)
0284 117 CONTINUE
0285 DC2 = AR1
0286 DS2 = AI1
0287 AR2 = AR1
0288 AI2 = AI1
0289 DO 119 J=3,IPPH
0290 JC = IPP2-J
0291 AR2H = DC2*AR2-DS2*AI2
0292 AI2 = DC2*AI2+DS2*AR2
0293 AR2 = AR2H
0294 DO 118 IK=1,IDL1
0295 C2(IK,L) = C2(IK,L)+AR2*CH2(IK,J)
0296 C2(IK,LC) = C2(IK,LC)+AI2*CH2(IK,JC)
0297 118 CONTINUE
0298 119 CONTINUE
0299 120 CONTINUE
0300 DO 122 J=2,IPPH
0301 DO 121 IK=1,IDL1
0302 CH2(IK,1) = CH2(IK,1)+CH2(IK,J)
0303 121 CONTINUE
0304 122 CONTINUE
0305 DO 124 J=2,IPPH
0306 JC = IPP2-J
0307 DO 123 K=1,L1
0308 CH(1,K,J) = C1(1,K,J)-C1(1,K,JC)
0309 CH(1,K,JC) = C1(1,K,J)+C1(1,K,JC)
0310 123 CONTINUE
0311 124 CONTINUE
0312 IF (IDO .EQ. 1) GO TO 132
0313 IF (NBD .LT. L1) GO TO 128
0314 DO 127 J=2,IPPH
0315 JC = IPP2-J
0316 DO 126 K=1,L1
0317 DO 125 I=3,IDO,2
0318 CH(I-1,K,J) = C1(I-1,K,J)-C1(I,K,JC)
0319 CH(I-1,K,JC) = C1(I-1,K,J)+C1(I,K,JC)
0320 CH(I,K,J) = C1(I,K,J)+C1(I-1,K,JC)
0321 CH(I,K,JC) = C1(I,K,J)-C1(I-1,K,JC)
0322 125 CONTINUE
0323 126 CONTINUE
0324 127 CONTINUE
0325 GO TO 132
0326 128 DO 131 J=2,IPPH
0327 JC = IPP2-J
0328 DO 130 I=3,IDO,2
0329 DO 129 K=1,L1
0330 CH(I-1,K,J) = C1(I-1,K,J)-C1(I,K,JC)
0331 CH(I-1,K,JC) = C1(I-1,K,J)+C1(I,K,JC)
0332 CH(I,K,J) = C1(I,K,J)+C1(I-1,K,JC)
0333 CH(I,K,JC) = C1(I,K,J)-C1(I-1,K,JC)
0334 129 CONTINUE
0335 130 CONTINUE
0336 131 CONTINUE
0337 132 CONTINUE
0338 IF (IDO .EQ. 1) RETURN
0339 DO 133 IK=1,IDL1
0340 C2(IK,1) = CH2(IK,1)
0341 133 CONTINUE
0342 DO 135 J=2,IP
0343 DO 134 K=1,L1
0344 C1(1,K,J) = CH(1,K,J)
0345 134 CONTINUE
0346 135 CONTINUE
0347 IF (NBD .GT. L1) GO TO 139
0348 IS = -IDO
0349 DO 138 J=2,IP
0350 IS = IS+IDO
0351 IDIJ = IS
0352 DO 137 I=3,IDO,2
0353 IDIJ = IDIJ+2
0354 DO 136 K=1,L1
0355 C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J)
0356 C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J)
0357 136 CONTINUE
0358 137 CONTINUE
0359 138 CONTINUE
0360 GO TO 143
0361 139 IS = -IDO
0362 DO 142 J=2,IP
0363 IS = IS+IDO
0364 DO 141 K=1,L1
0365 IDIJ = IS
0366 DO 140 I=3,IDO,2
0367 IDIJ = IDIJ+2
0368 C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J)
0369 C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J)
0370 140 CONTINUE
0371 141 CONTINUE
0372 142 CONTINUE
0373 143 RETURN
0374 END
0375 SUBROUTINE RFFTB1 (N,C,CH,WA,IFAC)
aea29c8517 Alis*0376
9bf05f3f81 Alis*0377 IMPLICIT REAL*8 (A-H,O-Z)
aea29c8517 Alis*0378 IMPLICIT INTEGER (I-N)
0379
3efde2d9d6 Jean*0380 DIMENSION CH(*) ,C(*) ,WA(*) ,IFAC(*)
42c525bfb4 Alis*0381 NF = IFAC(2)
0382 NA = 0
0383 L1 = 1
0384 IW = 1
0385 DO 116 K1=1,NF
0386 IP = IFAC(K1+2)
0387 L2 = IP*L1
0388 IDO = N/L2
0389 IDL1 = IDO*L1
0390 IF (IP .NE. 4) GO TO 103
0391 IX2 = IW+IDO
0392 IX3 = IX2+IDO
0393 IF (NA .NE. 0) GO TO 101
0394 CALL RADB4 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3))
0395 GO TO 102
0396 101 CALL RADB4 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3))
0397 102 NA = 1-NA
0398 GO TO 115
0399 103 IF (IP .NE. 2) GO TO 106
0400 IF (NA .NE. 0) GO TO 104
0401 CALL RADB2 (IDO,L1,C,CH,WA(IW))
0402 GO TO 105
0403 104 CALL RADB2 (IDO,L1,CH,C,WA(IW))
0404 105 NA = 1-NA
0405 GO TO 115
0406 106 IF (IP .NE. 3) GO TO 109
0407 IX2 = IW+IDO
0408 IF (NA .NE. 0) GO TO 107
0409 CALL RADB3 (IDO,L1,C,CH,WA(IW),WA(IX2))
0410 GO TO 108
0411 107 CALL RADB3 (IDO,L1,CH,C,WA(IW),WA(IX2))
0412 108 NA = 1-NA
0413 GO TO 115
0414 109 IF (IP .NE. 5) GO TO 112
0415 IX2 = IW+IDO
0416 IX3 = IX2+IDO
0417 IX4 = IX3+IDO
0418 IF (NA .NE. 0) GO TO 110
0419 CALL RADB5 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4))
0420 GO TO 111
0421 110 CALL RADB5 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4))
0422 111 NA = 1-NA
0423 GO TO 115
0424 112 IF (NA .NE. 0) GO TO 113
0425 CALL RADBG (IDO,IP,L1,IDL1,C,C,C,CH,CH,WA(IW))
0426 GO TO 114
0427 113 CALL RADBG (IDO,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW))
0428 114 IF (IDO .EQ. 1) NA = 1-NA
0429 115 L1 = L2
0430 IW = IW+(IP-1)*IDO
0431 116 CONTINUE
0432 IF (NA .EQ. 0) RETURN
0433 DO 117 I=1,N
0434 C(I) = CH(I)
0435 117 CONTINUE
0436 RETURN
0437 END
0438
aea29c8517 Alis*0439
9bf05f3f81 Alis*0440
aea29c8517 Alis*0441
42c525bfb4 Alis*0442
0443
0444
0445
0446
0447 SUBROUTINE RADF2 (IDO,L1,CC,CH,WA1)
aea29c8517 Alis*0448
9bf05f3f81 Alis*0449 IMPLICIT REAL*8 (A-H,O-Z)
aea29c8517 Alis*0450 IMPLICIT INTEGER (I-N)
0451
42c525bfb4 Alis*0452 DIMENSION CH(IDO,2,L1) ,CC(IDO,L1,2) ,
3efde2d9d6 Jean*0453 1 WA1(*)
42c525bfb4 Alis*0454 DO 101 K=1,L1
0455 CH(1,1,K) = CC(1,K,1)+CC(1,K,2)
0456 CH(IDO,2,K) = CC(1,K,1)-CC(1,K,2)
0457 101 CONTINUE
0458 IF (IDO-2) 107,105,102
0459 102 IDP2 = IDO+2
0460 DO 104 K=1,L1
0461 DO 103 I=3,IDO,2
0462 IC = IDP2-I
0463 TR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
0464 TI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
0465 CH(I,1,K) = CC(I,K,1)+TI2
0466 CH(IC,2,K) = TI2-CC(I,K,1)
0467 CH(I-1,1,K) = CC(I-1,K,1)+TR2
0468 CH(IC-1,2,K) = CC(I-1,K,1)-TR2
0469 103 CONTINUE
0470 104 CONTINUE
0471 IF (MOD(IDO,2) .EQ. 1) RETURN
0472 105 DO 106 K=1,L1
0473 CH(1,2,K) = -CC(IDO,K,2)
0474 CH(IDO,1,K) = CC(IDO,K,1)
0475 106 CONTINUE
0476 107 RETURN
0477 END
0478 SUBROUTINE RADF3 (IDO,L1,CC,CH,WA1,WA2)
aea29c8517 Alis*0479
9bf05f3f81 Alis*0480 IMPLICIT REAL*8 (A-H,O-Z)
aea29c8517 Alis*0481 IMPLICIT INTEGER (I-N)
0482
42c525bfb4 Alis*0483 DIMENSION CH(IDO,3,L1) ,CC(IDO,L1,3) ,
3efde2d9d6 Jean*0484 1 WA1(*) ,WA2(*)
9bf05f3f81 Alis*0485 DATA TAUR,TAUI /-.5D0,.866025403784439D0/
42c525bfb4 Alis*0486 DO 101 K=1,L1
0487 CR2 = CC(1,K,2)+CC(1,K,3)
0488 CH(1,1,K) = CC(1,K,1)+CR2
0489 CH(1,3,K) = TAUI*(CC(1,K,3)-CC(1,K,2))
0490 CH(IDO,2,K) = CC(1,K,1)+TAUR*CR2
0491 101 CONTINUE
0492 IF (IDO .EQ. 1) RETURN
0493 IDP2 = IDO+2
0494 DO 103 K=1,L1
0495 DO 102 I=3,IDO,2
0496 IC = IDP2-I
0497 DR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
0498 DI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
0499 DR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3)
0500 DI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3)
0501 CR2 = DR2+DR3
0502 CI2 = DI2+DI3
0503 CH(I-1,1,K) = CC(I-1,K,1)+CR2
0504 CH(I,1,K) = CC(I,K,1)+CI2
0505 TR2 = CC(I-1,K,1)+TAUR*CR2
0506 TI2 = CC(I,K,1)+TAUR*CI2
0507 TR3 = TAUI*(DI2-DI3)
0508 TI3 = TAUI*(DR3-DR2)
0509 CH(I-1,3,K) = TR2+TR3
0510 CH(IC-1,2,K) = TR2-TR3
0511 CH(I,3,K) = TI2+TI3
0512 CH(IC,2,K) = TI3-TI2
0513 102 CONTINUE
0514 103 CONTINUE
0515 RETURN
0516 END
0517 SUBROUTINE RADF4 (IDO,L1,CC,CH,WA1,WA2,WA3)
aea29c8517 Alis*0518
9bf05f3f81 Alis*0519 IMPLICIT REAL*8 (A-H,O-Z)
aea29c8517 Alis*0520 IMPLICIT INTEGER (I-N)
0521
42c525bfb4 Alis*0522 DIMENSION CC(IDO,L1,4) ,CH(IDO,4,L1) ,
3efde2d9d6 Jean*0523 1 WA1(*) ,WA2(*) ,WA3(*)
9bf05f3f81 Alis*0524 DATA HSQT2 /.7071067811865475D0/
42c525bfb4 Alis*0525 DO 101 K=1,L1
0526 TR1 = CC(1,K,2)+CC(1,K,4)
0527 TR2 = CC(1,K,1)+CC(1,K,3)
0528 CH(1,1,K) = TR1+TR2
0529 CH(IDO,4,K) = TR2-TR1
0530 CH(IDO,2,K) = CC(1,K,1)-CC(1,K,3)
0531 CH(1,3,K) = CC(1,K,4)-CC(1,K,2)
0532 101 CONTINUE
0533 IF (IDO-2) 107,105,102
0534 102 IDP2 = IDO+2
0535 DO 104 K=1,L1
0536 DO 103 I=3,IDO,2
0537 IC = IDP2-I
0538 CR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
0539 CI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
0540 CR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3)
0541 CI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3)
0542 CR4 = WA3(I-2)*CC(I-1,K,4)+WA3(I-1)*CC(I,K,4)
0543 CI4 = WA3(I-2)*CC(I,K,4)-WA3(I-1)*CC(I-1,K,4)
0544 TR1 = CR2+CR4
0545 TR4 = CR4-CR2
0546 TI1 = CI2+CI4
0547 TI4 = CI2-CI4
0548 TI2 = CC(I,K,1)+CI3
0549 TI3 = CC(I,K,1)-CI3
0550 TR2 = CC(I-1,K,1)+CR3
0551 TR3 = CC(I-1,K,1)-CR3
0552 CH(I-1,1,K) = TR1+TR2
0553 CH(IC-1,4,K) = TR2-TR1
0554 CH(I,1,K) = TI1+TI2
0555 CH(IC,4,K) = TI1-TI2
0556 CH(I-1,3,K) = TI4+TR3
0557 CH(IC-1,2,K) = TR3-TI4
0558 CH(I,3,K) = TR4+TI3
0559 CH(IC,2,K) = TR4-TI3
0560 103 CONTINUE
0561 104 CONTINUE
0562 IF (MOD(IDO,2) .EQ. 1) RETURN
0563 105 CONTINUE
0564 DO 106 K=1,L1
0565 TI1 = -HSQT2*(CC(IDO,K,2)+CC(IDO,K,4))
0566 TR1 = HSQT2*(CC(IDO,K,2)-CC(IDO,K,4))
0567 CH(IDO,1,K) = TR1+CC(IDO,K,1)
0568 CH(IDO,3,K) = CC(IDO,K,1)-TR1
0569 CH(1,2,K) = TI1-CC(IDO,K,3)
0570 CH(1,4,K) = TI1+CC(IDO,K,3)
0571 106 CONTINUE
0572 107 RETURN
0573 END
0574 SUBROUTINE RADF5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
aea29c8517 Alis*0575
9bf05f3f81 Alis*0576 IMPLICIT REAL*8 (A-H,O-Z)
aea29c8517 Alis*0577 IMPLICIT INTEGER (I-N)
0578
42c525bfb4 Alis*0579 DIMENSION CC(IDO,L1,5) ,CH(IDO,5,L1) ,
3efde2d9d6 Jean*0580 1 WA1(*) ,WA2(*) ,WA3(*) ,WA4(*)
9bf05f3f81 Alis*0581 DATA TR11,TI11,TR12,TI12 /.309016994374947D0,.951056516295154D0,
0582 1-.809016994374947D0,.587785252292473D0/
42c525bfb4 Alis*0583 DO 101 K=1,L1
0584 CR2 = CC(1,K,5)+CC(1,K,2)
0585 CI5 = CC(1,K,5)-CC(1,K,2)
0586 CR3 = CC(1,K,4)+CC(1,K,3)
0587 CI4 = CC(1,K,4)-CC(1,K,3)
0588 CH(1,1,K) = CC(1,K,1)+CR2+CR3
0589 CH(IDO,2,K) = CC(1,K,1)+TR11*CR2+TR12*CR3
0590 CH(1,3,K) = TI11*CI5+TI12*CI4
0591 CH(IDO,4,K) = CC(1,K,1)+TR12*CR2+TR11*CR3
0592 CH(1,5,K) = TI12*CI5-TI11*CI4
0593 101 CONTINUE
0594 IF (IDO .EQ. 1) RETURN
0595 IDP2 = IDO+2
0596 DO 103 K=1,L1
0597 DO 102 I=3,IDO,2
0598 IC = IDP2-I
0599 DR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
0600 DI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
0601 DR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3)
0602 DI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3)
0603 DR4 = WA3(I-2)*CC(I-1,K,4)+WA3(I-1)*CC(I,K,4)
0604 DI4 = WA3(I-2)*CC(I,K,4)-WA3(I-1)*CC(I-1,K,4)
0605 DR5 = WA4(I-2)*CC(I-1,K,5)+WA4(I-1)*CC(I,K,5)
0606 DI5 = WA4(I-2)*CC(I,K,5)-WA4(I-1)*CC(I-1,K,5)
0607 CR2 = DR2+DR5
0608 CI5 = DR5-DR2
0609 CR5 = DI2-DI5
0610 CI2 = DI2+DI5
0611 CR3 = DR3+DR4
0612 CI4 = DR4-DR3
0613 CR4 = DI3-DI4
0614 CI3 = DI3+DI4
0615 CH(I-1,1,K) = CC(I-1,K,1)+CR2+CR3
0616 CH(I,1,K) = CC(I,K,1)+CI2+CI3
0617 TR2 = CC(I-1,K,1)+TR11*CR2+TR12*CR3
0618 TI2 = CC(I,K,1)+TR11*CI2+TR12*CI3
0619 TR3 = CC(I-1,K,1)+TR12*CR2+TR11*CR3
0620 TI3 = CC(I,K,1)+TR12*CI2+TR11*CI3
0621 TR5 = TI11*CR5+TI12*CR4
0622 TI5 = TI11*CI5+TI12*CI4
0623 TR4 = TI12*CR5-TI11*CR4
0624 TI4 = TI12*CI5-TI11*CI4
0625 CH(I-1,3,K) = TR2+TR5
0626 CH(IC-1,2,K) = TR2-TR5
0627 CH(I,3,K) = TI2+TI5
0628 CH(IC,2,K) = TI5-TI2
0629 CH(I-1,5,K) = TR3+TR4
0630 CH(IC-1,4,K) = TR3-TR4
0631 CH(I,5,K) = TI3+TI4
0632 CH(IC,4,K) = TI4-TI3
0633 102 CONTINUE
0634 103 CONTINUE
0635 RETURN
0636 END
0637 SUBROUTINE RADFG (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
aea29c8517 Alis*0638
9bf05f3f81 Alis*0639 IMPLICIT REAL*8 (A-H,O-Z)
aea29c8517 Alis*0640 IMPLICIT INTEGER (I-N)
0641
42c525bfb4 Alis*0642 DIMENSION CH(IDO,L1,IP) ,CC(IDO,IP,L1) ,
0643 1 C1(IDO,L1,IP) ,C2(IDL1,IP),
3efde2d9d6 Jean*0644 2 CH2(IDL1,IP) ,WA(*)
9bf05f3f81 Alis*0645 DATA TPI/6.28318530717959D0/
42c525bfb4 Alis*0646 ARG = TPI/FLOAT(IP)
0647 DCP = COS(ARG)
0648 DSP = SIN(ARG)
0649 IPPH = (IP+1)/2
0650 IPP2 = IP+2
0651 IDP2 = IDO+2
0652 NBD = (IDO-1)/2
0653 IF (IDO .EQ. 1) GO TO 119
0654 DO 101 IK=1,IDL1
0655 CH2(IK,1) = C2(IK,1)
0656 101 CONTINUE
0657 DO 103 J=2,IP
0658 DO 102 K=1,L1
0659 CH(1,K,J) = C1(1,K,J)
0660 102 CONTINUE
0661 103 CONTINUE
0662 IF (NBD .GT. L1) GO TO 107
0663 IS = -IDO
0664 DO 106 J=2,IP
0665 IS = IS+IDO
0666 IDIJ = IS
0667 DO 105 I=3,IDO,2
0668 IDIJ = IDIJ+2
0669 DO 104 K=1,L1
0670 CH(I-1,K,J) = WA(IDIJ-1)*C1(I-1,K,J)+WA(IDIJ)*C1(I,K,J)
0671 CH(I,K,J) = WA(IDIJ-1)*C1(I,K,J)-WA(IDIJ)*C1(I-1,K,J)
0672 104 CONTINUE
0673 105 CONTINUE
0674 106 CONTINUE
0675 GO TO 111
0676 107 IS = -IDO
0677 DO 110 J=2,IP
0678 IS = IS+IDO
0679 DO 109 K=1,L1
0680 IDIJ = IS
0681 DO 108 I=3,IDO,2
0682 IDIJ = IDIJ+2
0683 CH(I-1,K,J) = WA(IDIJ-1)*C1(I-1,K,J)+WA(IDIJ)*C1(I,K,J)
0684 CH(I,K,J) = WA(IDIJ-1)*C1(I,K,J)-WA(IDIJ)*C1(I-1,K,J)
0685 108 CONTINUE
0686 109 CONTINUE
0687 110 CONTINUE
0688 111 IF (NBD .LT. L1) GO TO 115
0689 DO 114 J=2,IPPH
0690 JC = IPP2-J
0691 DO 113 K=1,L1
0692 DO 112 I=3,IDO,2
0693 C1(I-1,K,J) = CH(I-1,K,J)+CH(I-1,K,JC)
0694 C1(I-1,K,JC) = CH(I,K,J)-CH(I,K,JC)
0695 C1(I,K,J) = CH(I,K,J)+CH(I,K,JC)
0696 C1(I,K,JC) = CH(I-1,K,JC)-CH(I-1,K,J)
0697 112 CONTINUE
0698 113 CONTINUE
0699 114 CONTINUE
0700 GO TO 121
0701 115 DO 118 J=2,IPPH
0702 JC = IPP2-J
0703 DO 117 I=3,IDO,2
0704 DO 116 K=1,L1
0705 C1(I-1,K,J) = CH(I-1,K,J)+CH(I-1,K,JC)
0706 C1(I-1,K,JC) = CH(I,K,J)-CH(I,K,JC)
0707 C1(I,K,J) = CH(I,K,J)+CH(I,K,JC)
0708 C1(I,K,JC) = CH(I-1,K,JC)-CH(I-1,K,J)
0709 116 CONTINUE
0710 117 CONTINUE
0711 118 CONTINUE
0712 GO TO 121
0713 119 DO 120 IK=1,IDL1
0714 C2(IK,1) = CH2(IK,1)
0715 120 CONTINUE
0716 121 DO 123 J=2,IPPH
0717 JC = IPP2-J
0718 DO 122 K=1,L1
0719 C1(1,K,J) = CH(1,K,J)+CH(1,K,JC)
0720 C1(1,K,JC) = CH(1,K,JC)-CH(1,K,J)
0721 122 CONTINUE
0722 123 CONTINUE
0723
9bf05f3f81 Alis*0724 AR1 = 1.D0
0725 AI1 = 0.D0
42c525bfb4 Alis*0726 DO 127 L=2,IPPH
0727 LC = IPP2-L
0728 AR1H = DCP*AR1-DSP*AI1
0729 AI1 = DCP*AI1+DSP*AR1
0730 AR1 = AR1H
0731 DO 124 IK=1,IDL1
0732 CH2(IK,L) = C2(IK,1)+AR1*C2(IK,2)
0733 CH2(IK,LC) = AI1*C2(IK,IP)
0734 124 CONTINUE
0735 DC2 = AR1
0736 DS2 = AI1
0737 AR2 = AR1
0738 AI2 = AI1
0739 DO 126 J=3,IPPH
0740 JC = IPP2-J
0741 AR2H = DC2*AR2-DS2*AI2
0742 AI2 = DC2*AI2+DS2*AR2
0743 AR2 = AR2H
0744 DO 125 IK=1,IDL1
0745 CH2(IK,L) = CH2(IK,L)+AR2*C2(IK,J)
0746 CH2(IK,LC) = CH2(IK,LC)+AI2*C2(IK,JC)
0747 125 CONTINUE
0748 126 CONTINUE
0749 127 CONTINUE
0750 DO 129 J=2,IPPH
0751 DO 128 IK=1,IDL1
0752 CH2(IK,1) = CH2(IK,1)+C2(IK,J)
0753 128 CONTINUE
0754 129 CONTINUE
0755
0756 IF (IDO .LT. L1) GO TO 132
0757 DO 131 K=1,L1
0758 DO 130 I=1,IDO
0759 CC(I,1,K) = CH(I,K,1)
0760 130 CONTINUE
0761 131 CONTINUE
0762 GO TO 135
0763 132 DO 134 I=1,IDO
0764 DO 133 K=1,L1
0765 CC(I,1,K) = CH(I,K,1)
0766 133 CONTINUE
0767 134 CONTINUE
0768 135 DO 137 J=2,IPPH
0769 JC = IPP2-J
0770 J2 = J+J
0771 DO 136 K=1,L1
0772 CC(IDO,J2-2,K) = CH(1,K,J)
0773 CC(1,J2-1,K) = CH(1,K,JC)
0774 136 CONTINUE
0775 137 CONTINUE
0776 IF (IDO .EQ. 1) RETURN
0777 IF (NBD .LT. L1) GO TO 141
0778 DO 140 J=2,IPPH
0779 JC = IPP2-J
0780 J2 = J+J
0781 DO 139 K=1,L1
0782 DO 138 I=3,IDO,2
0783 IC = IDP2-I
0784 CC(I-1,J2-1,K) = CH(I-1,K,J)+CH(I-1,K,JC)
0785 CC(IC-1,J2-2,K) = CH(I-1,K,J)-CH(I-1,K,JC)
0786 CC(I,J2-1,K) = CH(I,K,J)+CH(I,K,JC)
0787 CC(IC,J2-2,K) = CH(I,K,JC)-CH(I,K,J)
0788 138 CONTINUE
0789 139 CONTINUE
0790 140 CONTINUE
0791 RETURN
0792 141 DO 144 J=2,IPPH
0793 JC = IPP2-J
0794 J2 = J+J
0795 DO 143 I=3,IDO,2
0796 IC = IDP2-I
0797 DO 142 K=1,L1
0798 CC(I-1,J2-1,K) = CH(I-1,K,J)+CH(I-1,K,JC)
0799 CC(IC-1,J2-2,K) = CH(I-1,K,J)-CH(I-1,K,JC)
0800 CC(I,J2-1,K) = CH(I,K,J)+CH(I,K,JC)
0801 CC(IC,J2-2,K) = CH(I,K,JC)-CH(I,K,J)
0802 142 CONTINUE
0803 143 CONTINUE
0804 144 CONTINUE
0805 RETURN
0806 END
0807 SUBROUTINE RFFTF1 (N,C,CH,WA,IFAC)
aea29c8517 Alis*0808
9bf05f3f81 Alis*0809 IMPLICIT REAL*8 (A-H,O-Z)
aea29c8517 Alis*0810 IMPLICIT INTEGER (I-N)
0811
3efde2d9d6 Jean*0812 DIMENSION CH(*) ,C(*) ,WA(*) ,IFAC(*)
42c525bfb4 Alis*0813 NF = IFAC(2)
0814 NA = 1
0815 L2 = N
0816 IW = N
0817 DO 111 K1=1,NF
0818 KH = NF-K1
0819 IP = IFAC(KH+3)
0820 L1 = L2/IP
0821 IDO = N/L2
0822 IDL1 = IDO*L1
0823 IW = IW-(IP-1)*IDO
0824 NA = 1-NA
0825 IF (IP .NE. 4) GO TO 102
0826 IX2 = IW+IDO
0827 IX3 = IX2+IDO
0828 IF (NA .NE. 0) GO TO 101
0829 CALL RADF4 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3))
0830 GO TO 110
0831 101 CALL RADF4 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3))
0832 GO TO 110
0833 102 IF (IP .NE. 2) GO TO 104
0834 IF (NA .NE. 0) GO TO 103
0835 CALL RADF2 (IDO,L1,C,CH,WA(IW))
0836 GO TO 110
0837 103 CALL RADF2 (IDO,L1,CH,C,WA(IW))
0838 GO TO 110
0839 104 IF (IP .NE. 3) GO TO 106
0840 IX2 = IW+IDO
0841 IF (NA .NE. 0) GO TO 105
0842 CALL RADF3 (IDO,L1,C,CH,WA(IW),WA(IX2))
0843 GO TO 110
0844 105 CALL RADF3 (IDO,L1,CH,C,WA(IW),WA(IX2))
0845 GO TO 110
0846 106 IF (IP .NE. 5) GO TO 108
0847 IX2 = IW+IDO
0848 IX3 = IX2+IDO
0849 IX4 = IX3+IDO
0850 IF (NA .NE. 0) GO TO 107
0851 CALL RADF5 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4))
0852 GO TO 110
0853 107 CALL RADF5 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4))
0854 GO TO 110
0855 108 IF (IDO .EQ. 1) NA = 1-NA
0856 IF (NA .NE. 0) GO TO 109
0857 CALL RADFG (IDO,IP,L1,IDL1,C,C,C,CH,CH,WA(IW))
0858 NA = 1
0859 GO TO 110
0860 109 CALL RADFG (IDO,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW))
0861 NA = 0
0862 110 L2 = L1
0863 111 CONTINUE
0864 IF (NA .EQ. 1) RETURN
0865 DO 112 I=1,N
0866 C(I) = CH(I)
0867 112 CONTINUE
0868 RETURN
0869 END
0870
aea29c8517 Alis*0871
9bf05f3f81 Alis*0872
aea29c8517 Alis*0873
0874
42c525bfb4 Alis*0875
0876
0877
0878
0879
0880 SUBROUTINE RFFTI1 (N,WA,IFAC)
aea29c8517 Alis*0881
9bf05f3f81 Alis*0882 IMPLICIT REAL*8 (A-H,O-Z)
aea29c8517 Alis*0883 IMPLICIT INTEGER (I-N)
0884
3efde2d9d6 Jean*0885 DIMENSION WA(*) ,IFAC(*) ,NTRYH(4)
42c525bfb4 Alis*0886 DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/4,2,3,5/
0887 NL = N
0888 NF = 0
0889 J = 0
0890 101 J = J+1
0891 IF (J-4) 102,102,103
0892 102 NTRY = NTRYH(J)
0893 GO TO 104
0894 103 NTRY = NTRY+2
0895 104 NQ = NL/NTRY
0896 NR = NL-NTRY*NQ
0897 IF (NR) 101,105,101
0898 105 NF = NF+1
0899 IFAC(NF+2) = NTRY
0900 NL = NQ
0901 IF (NTRY .NE. 2) GO TO 107
0902 IF (NF .EQ. 1) GO TO 107
0903 DO 106 I=2,NF
0904 IB = NF-I+2
0905 IFAC(IB+2) = IFAC(IB+1)
0906 106 CONTINUE
0907 IFAC(3) = 2
0908 107 IF (NL .NE. 1) GO TO 104
0909 IFAC(1) = N
0910 IFAC(2) = NF
9bf05f3f81 Alis*0911 TPI = 6.28318530717959D0
42c525bfb4 Alis*0912 ARGH = TPI/FLOAT(N)
0913 IS = 0
0914 NFM1 = NF-1
0915 L1 = 1
0916 IF (NFM1 .EQ. 0) RETURN
0917 DO 110 K1=1,NFM1
0918 IP = IFAC(K1+2)
0919 LD = 0
0920 L2 = L1*IP
0921 IDO = N/L2
0922 IPM = IP-1
0923 DO 109 J=1,IPM
0924 LD = LD+L1
0925 I = IS
0926 ARGLD = FLOAT(LD)*ARGH
9bf05f3f81 Alis*0927 FI = 0.D0
42c525bfb4 Alis*0928 DO 108 II=3,IDO,2
0929 I = I+2
9bf05f3f81 Alis*0930 FI = FI+1.D0
42c525bfb4 Alis*0931 ARG = FI*ARGLD
0932 WA(I-1) = COS(ARG)
0933 WA(I) = SIN(ARG)
0934 108 CONTINUE
0935 IS = IS+IDO
0936 109 CONTINUE
0937 L1 = L2
0938 110 CONTINUE
0939 RETURN
0940 END
0941
0942
0943
0944
0945
0946
0947
0948 SUBROUTINE R8ADB2 (IDO,L1,CC,CH,WA1)
0949 implicit real*8 (A-H,O-Z)
aea29c8517 Alis*0950 IMPLICIT INTEGER (I-N)
42c525bfb4 Alis*0951 DIMENSION CC(IDO,2,L1) ,CH(IDO,L1,2) ,
3efde2d9d6 Jean*0952 1 WA1(*)
42c525bfb4 Alis*0953 DO 101 K=1,L1
0954 CH(1,K,1) = CC(1,1,K)+CC(IDO,2,K)
0955 CH(1,K,2) = CC(1,1,K)-CC(IDO,2,K)
0956 101 CONTINUE
0957 IF (IDO-2) 107,105,102
0958 102 IDP2 = IDO+2
0959 DO 104 K=1,L1
0960 DO 103 I=3,IDO,2
0961 IC = IDP2-I
0962 CH(I-1,K,1) = CC(I-1,1,K)+CC(IC-1,2,K)
0963 TR2 = CC(I-1,1,K)-CC(IC-1,2,K)
0964 CH(I,K,1) = CC(I,1,K)-CC(IC,2,K)
0965 TI2 = CC(I,1,K)+CC(IC,2,K)
0966 CH(I-1,K,2) = WA1(I-2)*TR2-WA1(I-1)*TI2
0967 CH(I,K,2) = WA1(I-2)*TI2+WA1(I-1)*TR2
0968 103 CONTINUE
0969 104 CONTINUE
0970 IF (MOD(IDO,2) .EQ. 1) RETURN
0971 105 DO 106 K=1,L1
0972 CH(IDO,K,1) = CC(IDO,1,K)+CC(IDO,1,K)
0973 CH(IDO,K,2) = -(CC(1,2,K)+CC(1,2,K))
0974 106 CONTINUE
0975 107 RETURN
0976 END
0977 SUBROUTINE R8ADB3 (IDO,L1,CC,CH,WA1,WA2)
0978 implicit real*8 (A-H,O-Z)
aea29c8517 Alis*0979 IMPLICIT INTEGER (I-N)
42c525bfb4 Alis*0980 DIMENSION CC(IDO,3,L1) ,CH(IDO,L1,3) ,
3efde2d9d6 Jean*0981 1 WA1(*) ,WA2(*)
9bf05f3f81 Alis*0982 DATA TAUR,TAUI /-.5D0,.866025403784439D0/
42c525bfb4 Alis*0983 DO 101 K=1,L1
0984 TR2 = CC(IDO,2,K)+CC(IDO,2,K)
0985 CR2 = CC(1,1,K)+TAUR*TR2
0986 CH(1,K,1) = CC(1,1,K)+TR2
0987 CI3 = TAUI*(CC(1,3,K)+CC(1,3,K))
0988 CH(1,K,2) = CR2-CI3
0989 CH(1,K,3) = CR2+CI3
0990 101 CONTINUE
0991 IF (IDO .EQ. 1) RETURN
0992 IDP2 = IDO+2
0993 DO 103 K=1,L1
0994 DO 102 I=3,IDO,2
0995 IC = IDP2-I
0996 TR2 = CC(I-1,3,K)+CC(IC-1,2,K)
0997 CR2 = CC(I-1,1,K)+TAUR*TR2
0998 CH(I-1,K,1) = CC(I-1,1,K)+TR2
0999 TI2 = CC(I,3,K)-CC(IC,2,K)
1000 CI2 = CC(I,1,K)+TAUR*TI2
1001 CH(I,K,1) = CC(I,1,K)+TI2
1002 CR3 = TAUI*(CC(I-1,3,K)-CC(IC-1,2,K))
1003 CI3 = TAUI*(CC(I,3,K)+CC(IC,2,K))
1004 DR2 = CR2-CI3
1005 DR3 = CR2+CI3
1006 DI2 = CI2+CR3
1007 DI3 = CI2-CR3
1008 CH(I-1,K,2) = WA1(I-2)*DR2-WA1(I-1)*DI2
1009 CH(I,K,2) = WA1(I-2)*DI2+WA1(I-1)*DR2
1010 CH(I-1,K,3) = WA2(I-2)*DR3-WA2(I-1)*DI3
1011 CH(I,K,3) = WA2(I-2)*DI3+WA2(I-1)*DR3
1012 102 CONTINUE
1013 103 CONTINUE
1014 RETURN
1015 END
1016 SUBROUTINE R8ADB4 (IDO,L1,CC,CH,WA1,WA2,WA3)
1017 implicit real*8 (A-H,O-Z)
aea29c8517 Alis*1018 IMPLICIT INTEGER (I-N)
42c525bfb4 Alis*1019 DIMENSION CC(IDO,4,L1) ,CH(IDO,L1,4) ,
3efde2d9d6 Jean*1020 1 WA1(*) ,WA2(*) ,WA3(*)
9bf05f3f81 Alis*1021 DATA SQRT2 /1.414213562373095D0/
42c525bfb4 Alis*1022 DO 101 K=1,L1
1023 TR1 = CC(1,1,K)-CC(IDO,4,K)
1024 TR2 = CC(1,1,K)+CC(IDO,4,K)
1025 TR3 = CC(IDO,2,K)+CC(IDO,2,K)
1026 TR4 = CC(1,3,K)+CC(1,3,K)
1027 CH(1,K,1) = TR2+TR3
1028 CH(1,K,2) = TR1-TR4
1029 CH(1,K,3) = TR2-TR3
1030 CH(1,K,4) = TR1+TR4
1031 101 CONTINUE
1032 IF (IDO-2) 107,105,102
1033 102 IDP2 = IDO+2
1034 DO 104 K=1,L1
1035 DO 103 I=3,IDO,2
1036 IC = IDP2-I
1037 TI1 = CC(I,1,K)+CC(IC,4,K)
1038 TI2 = CC(I,1,K)-CC(IC,4,K)
1039 TI3 = CC(I,3,K)-CC(IC,2,K)
1040 TR4 = CC(I,3,K)+CC(IC,2,K)
1041 TR1 = CC(I-1,1,K)-CC(IC-1,4,K)
1042 TR2 = CC(I-1,1,K)+CC(IC-1,4,K)
1043 TI4 = CC(I-1,3,K)-CC(IC-1,2,K)
1044 TR3 = CC(I-1,3,K)+CC(IC-1,2,K)
1045 CH(I-1,K,1) = TR2+TR3
1046 CR3 = TR2-TR3
1047 CH(I,K,1) = TI2+TI3
1048 CI3 = TI2-TI3
1049 CR2 = TR1-TR4
1050 CR4 = TR1+TR4
1051 CI2 = TI1+TI4
1052 CI4 = TI1-TI4
1053 CH(I-1,K,2) = WA1(I-2)*CR2-WA1(I-1)*CI2
1054 CH(I,K,2) = WA1(I-2)*CI2+WA1(I-1)*CR2
1055 CH(I-1,K,3) = WA2(I-2)*CR3-WA2(I-1)*CI3
1056 CH(I,K,3) = WA2(I-2)*CI3+WA2(I-1)*CR3
1057 CH(I-1,K,4) = WA3(I-2)*CR4-WA3(I-1)*CI4
1058 CH(I,K,4) = WA3(I-2)*CI4+WA3(I-1)*CR4
1059 103 CONTINUE
1060 104 CONTINUE
1061 IF (MOD(IDO,2) .EQ. 1) RETURN
1062 105 CONTINUE
1063 DO 106 K=1,L1
1064 TI1 = CC(1,2,K)+CC(1,4,K)
1065 TI2 = CC(1,4,K)-CC(1,2,K)
1066 TR1 = CC(IDO,1,K)-CC(IDO,3,K)
1067 TR2 = CC(IDO,1,K)+CC(IDO,3,K)
1068 CH(IDO,K,1) = TR2+TR2
1069 CH(IDO,K,2) = SQRT2*(TR1-TI1)
1070 CH(IDO,K,3) = TI2+TI2
1071 CH(IDO,K,4) = -SQRT2*(TR1+TI1)
1072 106 CONTINUE
1073 107 RETURN
1074 END
1075 SUBROUTINE R8ADB5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
1076 implicit real*8 (A-H,O-Z)
aea29c8517 Alis*1077 IMPLICIT INTEGER (I-N)
42c525bfb4 Alis*1078 DIMENSION CC(IDO,5,L1) ,CH(IDO,L1,5) ,
3efde2d9d6 Jean*1079 1 WA1(*) ,WA2(*) ,WA3(*) ,WA4(*)
9bf05f3f81 Alis*1080 DATA TR11,TI11,TR12,TI12 /.309016994374947D0,.951056516295154D0,
1081 1-.809016994374947D0,.587785252292473D0/
42c525bfb4 Alis*1082 DO 101 K=1,L1
1083 TI5 = CC(1,3,K)+CC(1,3,K)
1084 TI4 = CC(1,5,K)+CC(1,5,K)
1085 TR2 = CC(IDO,2,K)+CC(IDO,2,K)
1086 TR3 = CC(IDO,4,K)+CC(IDO,4,K)
1087 CH(1,K,1) = CC(1,1,K)+TR2+TR3
1088 CR2 = CC(1,1,K)+TR11*TR2+TR12*TR3
1089 CR3 = CC(1,1,K)+TR12*TR2+TR11*TR3
1090 CI5 = TI11*TI5+TI12*TI4
1091 CI4 = TI12*TI5-TI11*TI4
1092 CH(1,K,2) = CR2-CI5
1093 CH(1,K,3) = CR3-CI4
1094 CH(1,K,4) = CR3+CI4
1095 CH(1,K,5) = CR2+CI5
1096 101 CONTINUE
1097 IF (IDO .EQ. 1) RETURN
1098 IDP2 = IDO+2
1099 DO 103 K=1,L1
1100 DO 102 I=3,IDO,2
1101 IC = IDP2-I
1102 TI5 = CC(I,3,K)+CC(IC,2,K)
1103 TI2 = CC(I,3,K)-CC(IC,2,K)
1104 TI4 = CC(I,5,K)+CC(IC,4,K)
1105 TI3 = CC(I,5,K)-CC(IC,4,K)
1106 TR5 = CC(I-1,3,K)-CC(IC-1,2,K)
1107 TR2 = CC(I-1,3,K)+CC(IC-1,2,K)
1108 TR4 = CC(I-1,5,K)-CC(IC-1,4,K)
1109 TR3 = CC(I-1,5,K)+CC(IC-1,4,K)
1110 CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3
1111 CH(I,K,1) = CC(I,1,K)+TI2+TI3
1112 CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3
1113 CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3
1114 CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3
1115 CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3
1116 CR5 = TI11*TR5+TI12*TR4
1117 CI5 = TI11*TI5+TI12*TI4
1118 CR4 = TI12*TR5-TI11*TR4
1119 CI4 = TI12*TI5-TI11*TI4
1120 DR3 = CR3-CI4
1121 DR4 = CR3+CI4
1122 DI3 = CI3+CR4
1123 DI4 = CI3-CR4
1124 DR5 = CR2+CI5
1125 DR2 = CR2-CI5
1126 DI5 = CI2-CR5
1127 DI2 = CI2+CR5
1128 CH(I-1,K,2) = WA1(I-2)*DR2-WA1(I-1)*DI2
1129 CH(I,K,2) = WA1(I-2)*DI2+WA1(I-1)*DR2
1130 CH(I-1,K,3) = WA2(I-2)*DR3-WA2(I-1)*DI3
1131 CH(I,K,3) = WA2(I-2)*DI3+WA2(I-1)*DR3
1132 CH(I-1,K,4) = WA3(I-2)*DR4-WA3(I-1)*DI4
1133 CH(I,K,4) = WA3(I-2)*DI4+WA3(I-1)*DR4
1134 CH(I-1,K,5) = WA4(I-2)*DR5-WA4(I-1)*DI5
1135 CH(I,K,5) = WA4(I-2)*DI5+WA4(I-1)*DR5
1136 102 CONTINUE
1137 103 CONTINUE
1138 RETURN
1139 END
1140 SUBROUTINE R8ADBG (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
1141 implicit real*8 (A-H,O-Z)
aea29c8517 Alis*1142 IMPLICIT INTEGER (I-N)
42c525bfb4 Alis*1143 DIMENSION CH(IDO,L1,IP) ,CC(IDO,IP,L1) ,
1144 1 C1(IDO,L1,IP) ,C2(IDL1,IP),
3efde2d9d6 Jean*1145 2 CH2(IDL1,IP) ,WA(*)
9bf05f3f81 Alis*1146 DATA TPI/6.28318530717959D0/
42c525bfb4 Alis*1147 ARG = TPI/FLOAT(IP)
1148 DCP = COS(ARG)
1149 DSP = SIN(ARG)
1150 IDP2 = IDO+2
1151 NBD = (IDO-1)/2
1152 IPP2 = IP+2
1153 IPPH = (IP+1)/2
1154 IF (IDO .LT. L1) GO TO 103
1155 DO 102 K=1,L1
1156 DO 101 I=1,IDO
1157 CH(I,K,1) = CC(I,1,K)
1158 101 CONTINUE
1159 102 CONTINUE
1160 GO TO 106
1161 103 DO 105 I=1,IDO
1162 DO 104 K=1,L1
1163 CH(I,K,1) = CC(I,1,K)
1164 104 CONTINUE
1165 105 CONTINUE
1166 106 DO 108 J=2,IPPH
1167 JC = IPP2-J
1168 J2 = J+J
1169 DO 107 K=1,L1
1170 CH(1,K,J) = CC(IDO,J2-2,K)+CC(IDO,J2-2,K)
1171 CH(1,K,JC) = CC(1,J2-1,K)+CC(1,J2-1,K)
1172 107 CONTINUE
1173 108 CONTINUE
1174 IF (IDO .EQ. 1) GO TO 116
1175 IF (NBD .LT. L1) GO TO 112
1176 DO 111 J=2,IPPH
1177 JC = IPP2-J
1178 DO 110 K=1,L1
1179 DO 109 I=3,IDO,2
1180 IC = IDP2-I
1181 CH(I-1,K,J) = CC(I-1,2*J-1,K)+CC(IC-1,2*J-2,K)
1182 CH(I-1,K,JC) = CC(I-1,2*J-1,K)-CC(IC-1,2*J-2,K)
1183 CH(I,K,J) = CC(I,2*J-1,K)-CC(IC,2*J-2,K)
1184 CH(I,K,JC) = CC(I,2*J-1,K)+CC(IC,2*J-2,K)
1185 109 CONTINUE
1186 110 CONTINUE
1187 111 CONTINUE
1188 GO TO 116
1189 112 DO 115 J=2,IPPH
1190 JC = IPP2-J
1191 DO 114 I=3,IDO,2
1192 IC = IDP2-I
1193 DO 113 K=1,L1
1194 CH(I-1,K,J) = CC(I-1,2*J-1,K)+CC(IC-1,2*J-2,K)
1195 CH(I-1,K,JC) = CC(I-1,2*J-1,K)-CC(IC-1,2*J-2,K)
1196 CH(I,K,J) = CC(I,2*J-1,K)-CC(IC,2*J-2,K)
1197 CH(I,K,JC) = CC(I,2*J-1,K)+CC(IC,2*J-2,K)
1198 113 CONTINUE
1199 114 CONTINUE
1200 115 CONTINUE
9bf05f3f81 Alis*1201 116 AR1 = 1.D0
1202 AI1 = 0.D0
42c525bfb4 Alis*1203 DO 120 L=2,IPPH
1204 LC = IPP2-L
1205 AR1H = DCP*AR1-DSP*AI1
1206 AI1 = DCP*AI1+DSP*AR1
1207 AR1 = AR1H
1208 DO 117 IK=1,IDL1
1209 C2(IK,L) = CH2(IK,1)+AR1*CH2(IK,2)
1210 C2(IK,LC) = AI1*CH2(IK,IP)
1211 117 CONTINUE
1212 DC2 = AR1
1213 DS2 = AI1
1214 AR2 = AR1
1215 AI2 = AI1
1216 DO 119 J=3,IPPH
1217 JC = IPP2-J
1218 AR2H = DC2*AR2-DS2*AI2
1219 AI2 = DC2*AI2+DS2*AR2
1220 AR2 = AR2H
1221 DO 118 IK=1,IDL1
1222 C2(IK,L) = C2(IK,L)+AR2*CH2(IK,J)
1223 C2(IK,LC) = C2(IK,LC)+AI2*CH2(IK,JC)
1224 118 CONTINUE
1225 119 CONTINUE
1226 120 CONTINUE
1227 DO 122 J=2,IPPH
1228 DO 121 IK=1,IDL1
1229 CH2(IK,1) = CH2(IK,1)+CH2(IK,J)
1230 121 CONTINUE
1231 122 CONTINUE
1232 DO 124 J=2,IPPH
1233 JC = IPP2-J
1234 DO 123 K=1,L1
1235 CH(1,K,J) = C1(1,K,J)-C1(1,K,JC)
1236 CH(1,K,JC) = C1(1,K,J)+C1(1,K,JC)
1237 123 CONTINUE
1238 124 CONTINUE
1239 IF (IDO .EQ. 1) GO TO 132
1240 IF (NBD .LT. L1) GO TO 128
1241 DO 127 J=2,IPPH
1242 JC = IPP2-J
1243 DO 126 K=1,L1
1244 DO 125 I=3,IDO,2
1245 CH(I-1,K,J) = C1(I-1,K,J)-C1(I,K,JC)
1246 CH(I-1,K,JC) = C1(I-1,K,J)+C1(I,K,JC)
1247 CH(I,K,J) = C1(I,K,J)+C1(I-1,K,JC)
1248 CH(I,K,JC) = C1(I,K,J)-C1(I-1,K,JC)
1249 125 CONTINUE
1250 126 CONTINUE
1251 127 CONTINUE
1252 GO TO 132
1253 128 DO 131 J=2,IPPH
1254 JC = IPP2-J
1255 DO 130 I=3,IDO,2
1256 DO 129 K=1,L1
1257 CH(I-1,K,J) = C1(I-1,K,J)-C1(I,K,JC)
1258 CH(I-1,K,JC) = C1(I-1,K,J)+C1(I,K,JC)
1259 CH(I,K,J) = C1(I,K,J)+C1(I-1,K,JC)
1260 CH(I,K,JC) = C1(I,K,J)-C1(I-1,K,JC)
1261 129 CONTINUE
1262 130 CONTINUE
1263 131 CONTINUE
1264 132 CONTINUE
1265 IF (IDO .EQ. 1) RETURN
1266 DO 133 IK=1,IDL1
1267 C2(IK,1) = CH2(IK,1)
1268 133 CONTINUE
1269 DO 135 J=2,IP
1270 DO 134 K=1,L1
1271 C1(1,K,J) = CH(1,K,J)
1272 134 CONTINUE
1273 135 CONTINUE
1274 IF (NBD .GT. L1) GO TO 139
1275 IS = -IDO
1276 DO 138 J=2,IP
1277 IS = IS+IDO
1278 IDIJ = IS
1279 DO 137 I=3,IDO,2
1280 IDIJ = IDIJ+2
1281 DO 136 K=1,L1
1282 C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J)
1283 C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J)
1284 136 CONTINUE
1285 137 CONTINUE
1286 138 CONTINUE
1287 GO TO 143
1288 139 IS = -IDO
1289 DO 142 J=2,IP
1290 IS = IS+IDO
1291 DO 141 K=1,L1
1292 IDIJ = IS
1293 DO 140 I=3,IDO,2
1294 IDIJ = IDIJ+2
1295 C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J)
1296 C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J)
1297 140 CONTINUE
1298 141 CONTINUE
1299 142 CONTINUE
1300 143 RETURN
1301 END
1302 SUBROUTINE R8FFTB1 (N,C,CH,WA,IFAC)
1303 implicit real*8 (A-H,O-Z)
aea29c8517 Alis*1304 IMPLICIT INTEGER (I-N)
3efde2d9d6 Jean*1305 DIMENSION CH(*) ,C(*) ,WA(*) ,IFAC(*)
42c525bfb4 Alis*1306 NF = IFAC(2)
1307 NA = 0
1308 L1 = 1
1309 IW = 1
1310 DO 116 K1=1,NF
1311 IP = IFAC(K1+2)
1312 L2 = IP*L1
1313 IDO = N/L2
1314 IDL1 = IDO*L1
1315 IF (IP .NE. 4) GO TO 103
1316 IX2 = IW+IDO
1317 IX3 = IX2+IDO
1318 IF (NA .NE. 0) GO TO 101
1319 CALL R8ADB4 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3))
1320 GO TO 102
1321 101 CALL R8ADB4 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3))
1322 102 NA = 1-NA
1323 GO TO 115
1324 103 IF (IP .NE. 2) GO TO 106
1325 IF (NA .NE. 0) GO TO 104
1326 CALL R8ADB2 (IDO,L1,C,CH,WA(IW))
1327 GO TO 105
1328 104 CALL R8ADB2 (IDO,L1,CH,C,WA(IW))
1329 105 NA = 1-NA
1330 GO TO 115
1331 106 IF (IP .NE. 3) GO TO 109
1332 IX2 = IW+IDO
1333 IF (NA .NE. 0) GO TO 107
1334 CALL R8ADB3 (IDO,L1,C,CH,WA(IW),WA(IX2))
1335 GO TO 108
1336 107 CALL R8ADB3 (IDO,L1,CH,C,WA(IW),WA(IX2))
1337 108 NA = 1-NA
1338 GO TO 115
1339 109 IF (IP .NE. 5) GO TO 112
1340 IX2 = IW+IDO
1341 IX3 = IX2+IDO
1342 IX4 = IX3+IDO
1343 IF (NA .NE. 0) GO TO 110
1344 CALL R8ADB5 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4))
1345 GO TO 111
1346 110 CALL R8ADB5 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4))
1347 111 NA = 1-NA
1348 GO TO 115
1349 112 IF (NA .NE. 0) GO TO 113
1350 CALL R8ADBG (IDO,IP,L1,IDL1,C,C,C,CH,CH,WA(IW))
1351 GO TO 114
1352 113 CALL R8ADBG (IDO,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW))
1353 114 IF (IDO .EQ. 1) NA = 1-NA
1354 115 L1 = L2
1355 IW = IW+(IP-1)*IDO
1356 116 CONTINUE
1357 IF (NA .EQ. 0) RETURN
1358 DO 117 I=1,N
1359 C(I) = CH(I)
1360 117 CONTINUE
1361 RETURN
1362 END
1363
1364
1365
1366
1367
1368
1369
1370 SUBROUTINE R8ADF2 (IDO,L1,CC,CH,WA1)
1371 implicit real*8 (A-H,O-Z)
aea29c8517 Alis*1372 IMPLICIT INTEGER (I-N)
42c525bfb4 Alis*1373 DIMENSION CH(IDO,2,L1) ,CC(IDO,L1,2) ,
3efde2d9d6 Jean*1374 1 WA1(*)
42c525bfb4 Alis*1375 DO 101 K=1,L1
1376 CH(1,1,K) = CC(1,K,1)+CC(1,K,2)
1377 CH(IDO,2,K) = CC(1,K,1)-CC(1,K,2)
1378 101 CONTINUE
1379 IF (IDO-2) 107,105,102
1380 102 IDP2 = IDO+2
1381 DO 104 K=1,L1
1382 DO 103 I=3,IDO,2
1383 IC = IDP2-I
1384 TR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
1385 TI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
1386 CH(I,1,K) = CC(I,K,1)+TI2
1387 CH(IC,2,K) = TI2-CC(I,K,1)
1388 CH(I-1,1,K) = CC(I-1,K,1)+TR2
1389 CH(IC-1,2,K) = CC(I-1,K,1)-TR2
1390 103 CONTINUE
1391 104 CONTINUE
1392 IF (MOD(IDO,2) .EQ. 1) RETURN
1393 105 DO 106 K=1,L1
1394 CH(1,2,K) = -CC(IDO,K,2)
1395 CH(IDO,1,K) = CC(IDO,K,1)
1396 106 CONTINUE
1397 107 RETURN
1398 END
1399 SUBROUTINE R8ADF3 (IDO,L1,CC,CH,WA1,WA2)
1400 implicit real*8 (A-H,O-Z)
aea29c8517 Alis*1401 IMPLICIT INTEGER (I-N)
42c525bfb4 Alis*1402 DIMENSION CH(IDO,3,L1) ,CC(IDO,L1,3) ,
3efde2d9d6 Jean*1403 1 WA1(*) ,WA2(*)
9bf05f3f81 Alis*1404 DATA TAUR,TAUI /-.5D0,.866025403784439D0/
42c525bfb4 Alis*1405 DO 101 K=1,L1
1406 CR2 = CC(1,K,2)+CC(1,K,3)
1407 CH(1,1,K) = CC(1,K,1)+CR2
1408 CH(1,3,K) = TAUI*(CC(1,K,3)-CC(1,K,2))
1409 CH(IDO,2,K) = CC(1,K,1)+TAUR*CR2
1410 101 CONTINUE
1411 IF (IDO .EQ. 1) RETURN
1412 IDP2 = IDO+2
1413 DO 103 K=1,L1
1414 DO 102 I=3,IDO,2
1415 IC = IDP2-I
1416 DR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
1417 DI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
1418 DR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3)
1419 DI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3)
1420 CR2 = DR2+DR3
1421 CI2 = DI2+DI3
1422 CH(I-1,1,K) = CC(I-1,K,1)+CR2
1423 CH(I,1,K) = CC(I,K,1)+CI2
1424 TR2 = CC(I-1,K,1)+TAUR*CR2
1425 TI2 = CC(I,K,1)+TAUR*CI2
1426 TR3 = TAUI*(DI2-DI3)
1427 TI3 = TAUI*(DR3-DR2)
1428 CH(I-1,3,K) = TR2+TR3
1429 CH(IC-1,2,K) = TR2-TR3
1430 CH(I,3,K) = TI2+TI3
1431 CH(IC,2,K) = TI3-TI2
1432 102 CONTINUE
1433 103 CONTINUE
1434 RETURN
1435 END
1436 SUBROUTINE R8ADF4 (IDO,L1,CC,CH,WA1,WA2,WA3)
1437 implicit real*8 (A-H,O-Z)
aea29c8517 Alis*1438 IMPLICIT INTEGER (I-N)
42c525bfb4 Alis*1439 DIMENSION CC(IDO,L1,4) ,CH(IDO,4,L1) ,
3efde2d9d6 Jean*1440 1 WA1(*) ,WA2(*) ,WA3(*)
9bf05f3f81 Alis*1441 DATA HSQT2 /.7071067811865475D0/
42c525bfb4 Alis*1442 DO 101 K=1,L1
1443 TR1 = CC(1,K,2)+CC(1,K,4)
1444 TR2 = CC(1,K,1)+CC(1,K,3)
1445 CH(1,1,K) = TR1+TR2
1446 CH(IDO,4,K) = TR2-TR1
1447 CH(IDO,2,K) = CC(1,K,1)-CC(1,K,3)
1448 CH(1,3,K) = CC(1,K,4)-CC(1,K,2)
1449 101 CONTINUE
1450 IF (IDO-2) 107,105,102
1451 102 IDP2 = IDO+2
1452 DO 104 K=1,L1
1453 DO 103 I=3,IDO,2
1454 IC = IDP2-I
1455 CR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
1456 CI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
1457 CR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3)
1458 CI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3)
1459 CR4 = WA3(I-2)*CC(I-1,K,4)+WA3(I-1)*CC(I,K,4)
1460 CI4 = WA3(I-2)*CC(I,K,4)-WA3(I-1)*CC(I-1,K,4)
1461 TR1 = CR2+CR4
1462 TR4 = CR4-CR2
1463 TI1 = CI2+CI4
1464 TI4 = CI2-CI4
1465 TI2 = CC(I,K,1)+CI3
1466 TI3 = CC(I,K,1)-CI3
1467 TR2 = CC(I-1,K,1)+CR3
1468 TR3 = CC(I-1,K,1)-CR3
1469 CH(I-1,1,K) = TR1+TR2
1470 CH(IC-1,4,K) = TR2-TR1
1471 CH(I,1,K) = TI1+TI2
1472 CH(IC,4,K) = TI1-TI2
1473 CH(I-1,3,K) = TI4+TR3
1474 CH(IC-1,2,K) = TR3-TI4
1475 CH(I,3,K) = TR4+TI3
1476 CH(IC,2,K) = TR4-TI3
1477 103 CONTINUE
1478 104 CONTINUE
1479 IF (MOD(IDO,2) .EQ. 1) RETURN
1480 105 CONTINUE
1481 DO 106 K=1,L1
1482 TI1 = -HSQT2*(CC(IDO,K,2)+CC(IDO,K,4))
1483 TR1 = HSQT2*(CC(IDO,K,2)-CC(IDO,K,4))
1484 CH(IDO,1,K) = TR1+CC(IDO,K,1)
1485 CH(IDO,3,K) = CC(IDO,K,1)-TR1
1486 CH(1,2,K) = TI1-CC(IDO,K,3)
1487 CH(1,4,K) = TI1+CC(IDO,K,3)
1488 106 CONTINUE
1489 107 RETURN
1490 END
1491 SUBROUTINE R8ADF5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
1492 implicit real*8 (A-H,O-Z)
aea29c8517 Alis*1493 IMPLICIT INTEGER (I-N)
42c525bfb4 Alis*1494 DIMENSION CC(IDO,L1,5) ,CH(IDO,5,L1) ,
3efde2d9d6 Jean*1495 1 WA1(*) ,WA2(*) ,WA3(*) ,WA4(*)
9bf05f3f81 Alis*1496 DATA TR11,TI11,TR12,TI12 /.309016994374947D0,.951056516295154D0,
1497 1-.809016994374947D0,.587785252292473D0/
42c525bfb4 Alis*1498 DO 101 K=1,L1
1499 CR2 = CC(1,K,5)+CC(1,K,2)
1500 CI5 = CC(1,K,5)-CC(1,K,2)
1501 CR3 = CC(1,K,4)+CC(1,K,3)
1502 CI4 = CC(1,K,4)-CC(1,K,3)
1503 CH(1,1,K) = CC(1,K,1)+CR2+CR3
1504 CH(IDO,2,K) = CC(1,K,1)+TR11*CR2+TR12*CR3
1505 CH(1,3,K) = TI11*CI5+TI12*CI4
1506 CH(IDO,4,K) = CC(1,K,1)+TR12*CR2+TR11*CR3
1507 CH(1,5,K) = TI12*CI5-TI11*CI4
1508 101 CONTINUE
1509 IF (IDO .EQ. 1) RETURN
1510 IDP2 = IDO+2
1511 DO 103 K=1,L1
1512 DO 102 I=3,IDO,2
1513 IC = IDP2-I
1514 DR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
1515 DI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
1516 DR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3)
1517 DI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3)
1518 DR4 = WA3(I-2)*CC(I-1,K,4)+WA3(I-1)*CC(I,K,4)
1519 DI4 = WA3(I-2)*CC(I,K,4)-WA3(I-1)*CC(I-1,K,4)
1520 DR5 = WA4(I-2)*CC(I-1,K,5)+WA4(I-1)*CC(I,K,5)
1521 DI5 = WA4(I-2)*CC(I,K,5)-WA4(I-1)*CC(I-1,K,5)
1522 CR2 = DR2+DR5
1523 CI5 = DR5-DR2
1524 CR5 = DI2-DI5
1525 CI2 = DI2+DI5
1526 CR3 = DR3+DR4
1527 CI4 = DR4-DR3
1528 CR4 = DI3-DI4
1529 CI3 = DI3+DI4
1530 CH(I-1,1,K) = CC(I-1,K,1)+CR2+CR3
1531 CH(I,1,K) = CC(I,K,1)+CI2+CI3
1532 TR2 = CC(I-1,K,1)+TR11*CR2+TR12*CR3
1533 TI2 = CC(I,K,1)+TR11*CI2+TR12*CI3
1534 TR3 = CC(I-1,K,1)+TR12*CR2+TR11*CR3
1535 TI3 = CC(I,K,1)+TR12*CI2+TR11*CI3
1536 TR5 = TI11*CR5+TI12*CR4
1537 TI5 = TI11*CI5+TI12*CI4
1538 TR4 = TI12*CR5-TI11*CR4
1539 TI4 = TI12*CI5-TI11*CI4
1540 CH(I-1,3,K) = TR2+TR5
1541 CH(IC-1,2,K) = TR2-TR5
1542 CH(I,3,K) = TI2+TI5
1543 CH(IC,2,K) = TI5-TI2
1544 CH(I-1,5,K) = TR3+TR4
1545 CH(IC-1,4,K) = TR3-TR4
1546 CH(I,5,K) = TI3+TI4
1547 CH(IC,4,K) = TI4-TI3
1548 102 CONTINUE
1549 103 CONTINUE
1550 RETURN
1551 END
1552 SUBROUTINE R8ADFG (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
1553 implicit real*8 (A-H,O-Z)
aea29c8517 Alis*1554 IMPLICIT INTEGER (I-N)
42c525bfb4 Alis*1555 DIMENSION CH(IDO,L1,IP) ,CC(IDO,IP,L1) ,
1556 1 C1(IDO,L1,IP) ,C2(IDL1,IP),
3efde2d9d6 Jean*1557 2 CH2(IDL1,IP) ,WA(*)
9bf05f3f81 Alis*1558 DATA TPI/6.28318530717959D0/
42c525bfb4 Alis*1559 ARG = TPI/FLOAT(IP)
1560 DCP = COS(ARG)
1561 DSP = SIN(ARG)
1562 IPPH = (IP+1)/2
1563 IPP2 = IP+2
1564 IDP2 = IDO+2
1565 NBD = (IDO-1)/2
1566 IF (IDO .EQ. 1) GO TO 119
1567 DO 101 IK=1,IDL1
1568 CH2(IK,1) = C2(IK,1)
1569 101 CONTINUE
1570 DO 103 J=2,IP
1571 DO 102 K=1,L1
1572 CH(1,K,J) = C1(1,K,J)
1573 102 CONTINUE
1574 103 CONTINUE
1575 IF (NBD .GT. L1) GO TO 107
1576 IS = -IDO
1577 DO 106 J=2,IP
1578 IS = IS+IDO
1579 IDIJ = IS
1580 DO 105 I=3,IDO,2
1581 IDIJ = IDIJ+2
1582 DO 104 K=1,L1
1583 CH(I-1,K,J) = WA(IDIJ-1)*C1(I-1,K,J)+WA(IDIJ)*C1(I,K,J)
1584 CH(I,K,J) = WA(IDIJ-1)*C1(I,K,J)-WA(IDIJ)*C1(I-1,K,J)
1585 104 CONTINUE
1586 105 CONTINUE
1587 106 CONTINUE
1588 GO TO 111
1589 107 IS = -IDO
1590 DO 110 J=2,IP
1591 IS = IS+IDO
1592 DO 109 K=1,L1
1593 IDIJ = IS
1594 DO 108 I=3,IDO,2
1595 IDIJ = IDIJ+2
1596 CH(I-1,K,J) = WA(IDIJ-1)*C1(I-1,K,J)+WA(IDIJ)*C1(I,K,J)
1597 CH(I,K,J) = WA(IDIJ-1)*C1(I,K,J)-WA(IDIJ)*C1(I-1,K,J)
1598 108 CONTINUE
1599 109 CONTINUE
1600 110 CONTINUE
1601 111 IF (NBD .LT. L1) GO TO 115
1602 DO 114 J=2,IPPH
1603 JC = IPP2-J
1604 DO 113 K=1,L1
1605 DO 112 I=3,IDO,2
1606 C1(I-1,K,J) = CH(I-1,K,J)+CH(I-1,K,JC)
1607 C1(I-1,K,JC) = CH(I,K,J)-CH(I,K,JC)
1608 C1(I,K,J) = CH(I,K,J)+CH(I,K,JC)
1609 C1(I,K,JC) = CH(I-1,K,JC)-CH(I-1,K,J)
1610 112 CONTINUE
1611 113 CONTINUE
1612 114 CONTINUE
1613 GO TO 121
1614 115 DO 118 J=2,IPPH
1615 JC = IPP2-J
1616 DO 117 I=3,IDO,2
1617 DO 116 K=1,L1
1618 C1(I-1,K,J) = CH(I-1,K,J)+CH(I-1,K,JC)
1619 C1(I-1,K,JC) = CH(I,K,J)-CH(I,K,JC)
1620 C1(I,K,J) = CH(I,K,J)+CH(I,K,JC)
1621 C1(I,K,JC) = CH(I-1,K,JC)-CH(I-1,K,J)
1622 116 CONTINUE
1623 117 CONTINUE
1624 118 CONTINUE
1625 GO TO 121
1626 119 DO 120 IK=1,IDL1
1627 C2(IK,1) = CH2(IK,1)
1628 120 CONTINUE
1629 121 DO 123 J=2,IPPH
1630 JC = IPP2-J
1631 DO 122 K=1,L1
1632 C1(1,K,J) = CH(1,K,J)+CH(1,K,JC)
1633 C1(1,K,JC) = CH(1,K,JC)-CH(1,K,J)
1634 122 CONTINUE
1635 123 CONTINUE
1636
9bf05f3f81 Alis*1637 AR1 = 1.D0
1638 AI1 = 0.D0
42c525bfb4 Alis*1639 DO 127 L=2,IPPH
1640 LC = IPP2-L
1641 AR1H = DCP*AR1-DSP*AI1
1642 AI1 = DCP*AI1+DSP*AR1
1643 AR1 = AR1H
1644 DO 124 IK=1,IDL1
1645 CH2(IK,L) = C2(IK,1)+AR1*C2(IK,2)
1646 CH2(IK,LC) = AI1*C2(IK,IP)
1647 124 CONTINUE
1648 DC2 = AR1
1649 DS2 = AI1
1650 AR2 = AR1
1651 AI2 = AI1
1652 DO 126 J=3,IPPH
1653 JC = IPP2-J
1654 AR2H = DC2*AR2-DS2*AI2
1655 AI2 = DC2*AI2+DS2*AR2
1656 AR2 = AR2H
1657 DO 125 IK=1,IDL1
1658 CH2(IK,L) = CH2(IK,L)+AR2*C2(IK,J)
1659 CH2(IK,LC) = CH2(IK,LC)+AI2*C2(IK,JC)
1660 125 CONTINUE
1661 126 CONTINUE
1662 127 CONTINUE
1663 DO 129 J=2,IPPH
1664 DO 128 IK=1,IDL1
1665 CH2(IK,1) = CH2(IK,1)+C2(IK,J)
1666 128 CONTINUE
1667 129 CONTINUE
1668
1669 IF (IDO .LT. L1) GO TO 132
1670 DO 131 K=1,L1
1671 DO 130 I=1,IDO
1672 CC(I,1,K) = CH(I,K,1)
1673 130 CONTINUE
1674 131 CONTINUE
1675 GO TO 135
1676 132 DO 134 I=1,IDO
1677 DO 133 K=1,L1
1678 CC(I,1,K) = CH(I,K,1)
1679 133 CONTINUE
1680 134 CONTINUE
1681 135 DO 137 J=2,IPPH
1682 JC = IPP2-J
1683 J2 = J+J
1684 DO 136 K=1,L1
1685 CC(IDO,J2-2,K) = CH(1,K,J)
1686 CC(1,J2-1,K) = CH(1,K,JC)
1687 136 CONTINUE
1688 137 CONTINUE
1689 IF (IDO .EQ. 1) RETURN
1690 IF (NBD .LT. L1) GO TO 141
1691 DO 140 J=2,IPPH
1692 JC = IPP2-J
1693 J2 = J+J
1694 DO 139 K=1,L1
1695 DO 138 I=3,IDO,2
1696 IC = IDP2-I
1697 CC(I-1,J2-1,K) = CH(I-1,K,J)+CH(I-1,K,JC)
1698 CC(IC-1,J2-2,K) = CH(I-1,K,J)-CH(I-1,K,JC)
1699 CC(I,J2-1,K) = CH(I,K,J)+CH(I,K,JC)
1700 CC(IC,J2-2,K) = CH(I,K,JC)-CH(I,K,J)
1701 138 CONTINUE
1702 139 CONTINUE
1703 140 CONTINUE
1704 RETURN
1705 141 DO 144 J=2,IPPH
1706 JC = IPP2-J
1707 J2 = J+J
1708 DO 143 I=3,IDO,2
1709 IC = IDP2-I
1710 DO 142 K=1,L1
1711 CC(I-1,J2-1,K) = CH(I-1,K,J)+CH(I-1,K,JC)
1712 CC(IC-1,J2-2,K) = CH(I-1,K,J)-CH(I-1,K,JC)
1713 CC(I,J2-1,K) = CH(I,K,J)+CH(I,K,JC)
1714 CC(IC,J2-2,K) = CH(I,K,JC)-CH(I,K,J)
1715 142 CONTINUE
1716 143 CONTINUE
1717 144 CONTINUE
1718 RETURN
1719 END
1720 SUBROUTINE R8FFTF1 (N,C,CH,WA,IFAC)
1721 implicit real*8 (A-H,O-Z)
aea29c8517 Alis*1722 IMPLICIT INTEGER (I-N)
3efde2d9d6 Jean*1723 DIMENSION CH(*) ,C(*) ,WA(*) ,IFAC(*)
42c525bfb4 Alis*1724 NF = IFAC(2)
1725 NA = 1
1726 L2 = N
1727 IW = N
1728 DO 111 K1=1,NF
1729 KH = NF-K1
1730 IP = IFAC(KH+3)
1731 L1 = L2/IP
1732 IDO = N/L2
1733 IDL1 = IDO*L1
1734 IW = IW-(IP-1)*IDO
1735 NA = 1-NA
1736 IF (IP .NE. 4) GO TO 102
1737 IX2 = IW+IDO
1738 IX3 = IX2+IDO
1739 IF (NA .NE. 0) GO TO 101
1740 CALL R8ADF4 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3))
1741 GO TO 110
1742 101 CALL R8ADF4 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3))
1743 GO TO 110
1744 102 IF (IP .NE. 2) GO TO 104
1745 IF (NA .NE. 0) GO TO 103
1746 CALL R8ADF2 (IDO,L1,C,CH,WA(IW))
1747 GO TO 110
1748 103 CALL R8ADF2 (IDO,L1,CH,C,WA(IW))
1749 GO TO 110
1750 104 IF (IP .NE. 3) GO TO 106
1751 IX2 = IW+IDO
1752 IF (NA .NE. 0) GO TO 105
1753 CALL R8ADF3 (IDO,L1,C,CH,WA(IW),WA(IX2))
1754 GO TO 110
1755 105 CALL R8ADF3 (IDO,L1,CH,C,WA(IW),WA(IX2))
1756 GO TO 110
1757 106 IF (IP .NE. 5) GO TO 108
1758 IX2 = IW+IDO
1759 IX3 = IX2+IDO
1760 IX4 = IX3+IDO
1761 IF (NA .NE. 0) GO TO 107
1762 CALL R8ADF5 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4))
1763 GO TO 110
1764 107 CALL R8ADF5 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4))
1765 GO TO 110
1766 108 IF (IDO .EQ. 1) NA = 1-NA
1767 IF (NA .NE. 0) GO TO 109
1768 CALL R8ADFG (IDO,IP,L1,IDL1,C,C,C,CH,CH,WA(IW))
1769 NA = 1
1770 GO TO 110
1771 109 CALL R8ADFG (IDO,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW))
1772 NA = 0
1773 110 L2 = L1
1774 111 CONTINUE
1775 IF (NA .EQ. 1) RETURN
1776 DO 112 I=1,N
1777 C(I) = CH(I)
1778 112 CONTINUE
1779 RETURN
1780 END
1781
1782
1783
1784
1785
1786
1787
1788 SUBROUTINE R8FFTI1 (N,WA,IFAC)
1789 implicit real*8 (A-H,O-Z)
aea29c8517 Alis*1790 IMPLICIT INTEGER (I-N)
3efde2d9d6 Jean*1791 DIMENSION WA(*) ,IFAC(*) ,NTRYH(4)
42c525bfb4 Alis*1792 DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/4,2,3,5/
1793 NL = N
1794 NF = 0
1795 J = 0
1796 101 J = J+1
1797 IF (J-4) 102,102,103
1798 102 NTRY = NTRYH(J)
1799 GO TO 104
1800 103 NTRY = NTRY+2
1801 104 NQ = NL/NTRY
1802 NR = NL-NTRY*NQ
1803 IF (NR) 101,105,101
1804 105 NF = NF+1
1805 IFAC(NF+2) = NTRY
1806 NL = NQ
1807 IF (NTRY .NE. 2) GO TO 107
1808 IF (NF .EQ. 1) GO TO 107
1809 DO 106 I=2,NF
1810 IB = NF-I+2
1811 IFAC(IB+2) = IFAC(IB+1)
1812 106 CONTINUE
1813 IFAC(3) = 2
1814 107 IF (NL .NE. 1) GO TO 104
1815 IFAC(1) = N
1816 IFAC(2) = NF
9bf05f3f81 Alis*1817 TPI = 6.28318530717959D0
42c525bfb4 Alis*1818 ARGH = TPI/FLOAT(N)
1819 IS = 0
1820 NFM1 = NF-1
1821 L1 = 1
1822 IF (NFM1 .EQ. 0) RETURN
1823 DO 110 K1=1,NFM1
1824 IP = IFAC(K1+2)
1825 LD = 0
1826 L2 = L1*IP
1827 IDO = N/L2
1828 IPM = IP-1
1829 DO 109 J=1,IPM
1830 LD = LD+L1
1831 I = IS
1832 ARGLD = FLOAT(LD)*ARGH
9bf05f3f81 Alis*1833 FI = 0.D0
42c525bfb4 Alis*1834 DO 108 II=3,IDO,2
1835 I = I+2
9bf05f3f81 Alis*1836 FI = FI+1.D0
42c525bfb4 Alis*1837 ARG = FI*ARGLD
1838 WA(I-1) = COS(ARG)
1839 WA(I) = SIN(ARG)
1840 108 CONTINUE
1841 IS = IS+IDO
1842 109 CONTINUE
1843 L1 = L2
1844 110 CONTINUE
1845 RETURN
6d54cf9ca1 Ed H*1846
1847 #endif /* ALLOW_ZONAL_FILT */
1848
42c525bfb4 Alis*1849 END
