12 #include "../base/base_uses.f90"
17 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'ps_wavelet_fft3d'
34 INTEGER,
INTENT(in) :: n
35 INTEGER,
INTENT(out) :: n_next
37 INTEGER,
PARAMETER :: ndata = 149, ndata1024 = 149
38 INTEGER,
DIMENSION(ndata),
PARAMETER :: idata = (/3, 4, 5, 6, 8, 9, 12, 15, 16, 18, 20, 24, &
39 25, 27, 30, 32, 36, 40, 45, 48, 54, 60, 64, 72, 75, 80, 81, 90, 96, 100, 108, 120, 125, &
40 128, 135, 144, 150, 160, 162, 180, 192, 200, 216, 225, 240, 243, 256, 270, 288, 300, 320, &
41 324, 360, 375, 384, 400, 405, 432, 450, 480, 486, 500, 512, 540, 576, 600, 625, 640, 648, &
42 675, 720, 729, 750, 768, 800, 810, 864, 900, 960, 972, 1000, 1024, 1080, 1125, 1152, 1200,&
43 1215, 1280, 1296, 1350, 1440, 1458, 1500, 1536, 1600, 1620, 1728, 1800, 1875, 1920, 1944, &
44 2000, 2025, 2048, 2160, 2250, 2304, 2400, 2430, 2500, 2560, 2592, 2700, 2880, 3000, 3072, &
45 3125, 3200, 3240, 3375, 3456, 3600, 3750, 3840, 3888, 4000, 4050, 4096, 4320, 4500, 4608, &
46 4800, 5000, 5120, 5184, 5400, 5625, 5760, 6000, 6144, 6400, 6480, 6750, 6912, 7200, 7500, &
53 loop_data:
DO i = 1, ndata1024
54 IF (n <= idata(i))
THEN
59 WRITE (unit=*, fmt=*)
"fourier_dim: ", n,
" is bigger than ", idata(ndata1024)
104 SUBROUTINE ctrig(n, trig, after, before, now, isign, ic)
114 REAL(kind=
dp) :: trig
115 INTEGER :: after, before, now, isign, ic
117 INTEGER :: i, itt, j, nh
118 REAL(kind=
dp) ::
angle, trigc, trigs, twopi
120 dimension now(7), after(7), before(7), trig(2,
ctrig_length)
121 INTEGER,
DIMENSION(7, 149) :: idata
123 DATA((idata(i, j), i=1, 7), j=1, 76)/ &
124 3, 3, 1, 1, 1, 1, 1, 4, 4, 1, 1, 1, 1, 1, &
125 5, 5, 1, 1, 1, 1, 1, 6, 6, 1, 1, 1, 1, 1, &
126 8, 8, 1, 1, 1, 1, 1, 9, 3, 3, 1, 1, 1, 1, &
127 12, 4, 3, 1, 1, 1, 1, 15, 5, 3, 1, 1, 1, 1, &
128 16, 4, 4, 1, 1, 1, 1, 18, 6, 3, 1, 1, 1, 1, &
129 20, 5, 4, 1, 1, 1, 1, 24, 8, 3, 1, 1, 1, 1, &
130 25, 5, 5, 1, 1, 1, 1, 27, 3, 3, 3, 1, 1, 1, &
131 30, 6, 5, 1, 1, 1, 1, 32, 8, 4, 1, 1, 1, 1, &
132 36, 4, 3, 3, 1, 1, 1, 40, 8, 5, 1, 1, 1, 1, &
133 45, 5, 3, 3, 1, 1, 1, 48, 4, 4, 3, 1, 1, 1, &
134 54, 6, 3, 3, 1, 1, 1, 60, 5, 4, 3, 1, 1, 1, &
135 64, 8, 8, 1, 1, 1, 1, 72, 8, 3, 3, 1, 1, 1, &
136 75, 5, 5, 3, 1, 1, 1, 80, 5, 4, 4, 1, 1, 1, &
137 81, 3, 3, 3, 3, 1, 1, 90, 6, 5, 3, 1, 1, 1, &
138 96, 8, 4, 3, 1, 1, 1, 100, 5, 5, 4, 1, 1, 1, &
139 108, 4, 3, 3, 3, 1, 1, 120, 8, 5, 3, 1, 1, 1, &
140 125, 5, 5, 5, 1, 1, 1, 128, 8, 4, 4, 1, 1, 1, &
141 135, 5, 3, 3, 3, 1, 1, 144, 6, 8, 3, 1, 1, 1, &
142 150, 6, 5, 5, 1, 1, 1, 160, 8, 5, 4, 1, 1, 1, &
143 162, 6, 3, 3, 3, 1, 1, 180, 5, 4, 3, 3, 1, 1, &
144 192, 6, 8, 4, 1, 1, 1, 200, 8, 5, 5, 1, 1, 1, &
145 216, 8, 3, 3, 3, 1, 1, 225, 5, 5, 3, 3, 1, 1, &
146 240, 6, 8, 5, 1, 1, 1, 243, 3, 3, 3, 3, 3, 1, &
147 256, 8, 8, 4, 1, 1, 1, 270, 6, 5, 3, 3, 1, 1, &
148 288, 8, 4, 3, 3, 1, 1, 300, 5, 5, 4, 3, 1, 1, &
149 320, 5, 4, 4, 4, 1, 1, 324, 4, 3, 3, 3, 3, 1, &
150 360, 8, 5, 3, 3, 1, 1, 375, 5, 5, 5, 3, 1, 1, &
151 384, 8, 4, 4, 3, 1, 1, 400, 5, 5, 4, 4, 1, 1, &
152 405, 5, 3, 3, 3, 3, 1, 432, 4, 4, 3, 3, 3, 1, &
153 450, 6, 5, 5, 3, 1, 1, 480, 8, 5, 4, 3, 1, 1, &
154 486, 6, 3, 3, 3, 3, 1, 500, 5, 5, 5, 4, 1, 1, &
155 512, 8, 8, 8, 1, 1, 1, 540, 5, 4, 3, 3, 3, 1, &
156 576, 4, 4, 4, 3, 3, 1, 600, 8, 5, 5, 3, 1, 1, &
157 625, 5, 5, 5, 5, 1, 1, 640, 8, 5, 4, 4, 1, 1, &
158 648, 8, 3, 3, 3, 3, 1, 675, 5, 5, 3, 3, 3, 1, &
159 720, 5, 4, 4, 3, 3, 1, 729, 3, 3, 3, 3, 3, 3, &
160 750, 6, 5, 5, 5, 1, 1, 768, 4, 4, 4, 4, 3, 1, &
161 800, 8, 5, 5, 4, 1, 1, 810, 6, 5, 3, 3, 3, 1/
162 DATA((idata(i, j), i=1, 7), j=77, 149)/ &
163 864, 8, 4, 3, 3, 3, 1, 900, 5, 5, 4, 3, 3, 1, &
164 960, 5, 4, 4, 4, 3, 1, 972, 4, 3, 3, 3, 3, 3, &
165 1000, 8, 5, 5, 5, 1, 1, 1024, 4, 4, 4, 4, 4, 1, &
166 1080, 6, 5, 4, 3, 3, 1, 1125, 5, 5, 5, 3, 3, 1, &
167 1152, 6, 4, 4, 4, 3, 1, 1200, 6, 8, 5, 5, 1, 1, &
168 1215, 5, 3, 3, 3, 3, 3, 1280, 8, 8, 5, 4, 1, 1, &
169 1296, 6, 8, 3, 3, 3, 1, 1350, 6, 5, 5, 3, 3, 1, &
170 1440, 6, 5, 4, 4, 3, 1, 1458, 6, 3, 3, 3, 3, 3, &
171 1500, 5, 5, 5, 4, 3, 1, 1536, 6, 8, 8, 4, 1, 1, &
172 1600, 8, 8, 5, 5, 1, 1, 1620, 5, 4, 3, 3, 3, 3, &
173 1728, 6, 8, 4, 3, 3, 1, 1800, 6, 5, 5, 4, 3, 1, &
174 1875, 5, 5, 5, 5, 3, 1, 1920, 6, 5, 4, 4, 4, 1, &
175 1944, 6, 4, 3, 3, 3, 3, 2000, 5, 5, 5, 4, 4, 1, &
176 2025, 5, 5, 3, 3, 3, 3, 2048, 8, 4, 4, 4, 4, 1, &
177 2160, 6, 8, 5, 3, 3, 1, 2250, 6, 5, 5, 5, 3, 1, &
178 2304, 6, 8, 4, 4, 3, 1, 2400, 6, 5, 5, 4, 4, 1, &
179 2430, 6, 5, 3, 3, 3, 3, 2500, 5, 5, 5, 5, 4, 1, &
180 2560, 8, 5, 4, 4, 4, 1, 2592, 6, 4, 4, 3, 3, 3, &
181 2700, 5, 5, 4, 3, 3, 3, 2880, 6, 8, 5, 4, 3, 1, &
182 3000, 6, 5, 5, 5, 4, 1, 3072, 6, 8, 4, 4, 4, 1, &
183 3125, 5, 5, 5, 5, 5, 1, 3200, 8, 5, 5, 4, 4, 1, &
184 3240, 6, 5, 4, 3, 3, 3, 3375, 5, 5, 5, 3, 3, 3, &
185 3456, 6, 4, 4, 4, 3, 3, 3600, 6, 8, 5, 5, 3, 1, &
186 3750, 6, 5, 5, 5, 5, 1, 3840, 6, 8, 5, 4, 4, 1, &
187 3888, 6, 8, 3, 3, 3, 3, 4000, 8, 5, 5, 5, 4, 1, &
188 4050, 6, 5, 5, 3, 3, 3, 4096, 8, 8, 4, 4, 4, 1, &
189 4320, 6, 5, 4, 4, 3, 3, 4500, 5, 5, 5, 4, 3, 3, &
190 4608, 6, 8, 8, 4, 3, 1, 4800, 6, 8, 5, 5, 4, 1, &
191 5000, 8, 5, 5, 5, 5, 1, 5120, 8, 8, 5, 4, 4, 1, &
192 5184, 6, 8, 4, 3, 3, 3, 5400, 6, 5, 5, 4, 3, 3, &
193 5625, 5, 5, 5, 5, 3, 3, 5760, 6, 8, 8, 5, 3, 1, &
194 6000, 6, 8, 5, 5, 5, 1, 6144, 6, 8, 8, 4, 4, 1, &
195 6400, 8, 8, 5, 5, 4, 1, 6480, 6, 8, 5, 3, 3, 3, &
196 6750, 6, 5, 5, 5, 3, 3, 6912, 6, 8, 4, 4, 3, 3, &
197 7200, 6, 5, 5, 4, 4, 3, 7500, 5, 5, 5, 5, 4, 3, &
198 7680, 6, 8, 8, 5, 4, 1, 8000, 8, 8, 5, 5, 5, 1, &
199 8192, 8, 8, 8, 4, 4, 1/
203 print *,
'VALUE OF', n,
'NOT ALLOWED FOR FFT, ALLOWED VALUES ARE:'
205 WRITE (*, 37) (idata(1, j), j=1, 149)
208 IF (n .EQ. idata(1, i))
THEN
211 itt = idata(1 + j, i)
214 now(j) = idata(1 + j, i)
226 after(i) = after(i - 1)*now(i - 1)
227 before(ic - i + 1) = before(ic - i + 2)*now(ic - i + 2)
230 twopi = 6.283185307179586_dp
231 angle = isign*twopi/n
232 IF (mod(n, 2) .EQ. 0)
THEN
236 trig(1, nh + 1) = -1._dp
237 trig(2, nh + 1) = 0._dp
241 trig(1, i + 1) = trigc
242 trig(2, i + 1) = trigs
243 trig(1, n - i + 1) = trigc
244 trig(2, n - i + 1) = -trigs
253 trig(1, i + 1) = trigc
254 trig(2, i + 1) = trigs
255 trig(1, n - i + 1) = trigc
256 trig(2, n - i + 1) = -trigs
279 SUBROUTINE fftstp(mm, nfft, m, nn, n, zin, zout, trig, after, now, before, isign)
285 INTEGER :: mm, nfft, m, nn, n
286 REAL(kind=
dp) :: zin, zout, trig
287 INTEGER :: after, now, before, isign
289 INTEGER :: atb, atn, ia, ias, ib, itrig, itt, j, &
290 nin1, nin2, nin3, nin4, nin5, nin6, &
291 nin7, nin8, nout1, nout2, nout3, &
292 nout4, nout5, nout6, nout7, nout8
293 REAL(kind=
dp) :: am, ap, bb, bm, bp, ci2, ci3, ci4, ci5, ci6, ci7, ci8, cm, cos2, cos4, cp, &
294 cr2, cr3, cr4, cr5, cr6, cr7, cr8, dm, dpp, r, r1, r2, r25, r3, r34, r4, r5, r6, r7, r8, &
295 rt2i, s, s1, s2, s25, s3, s34, s4, s5, s6, s7, s8, sin2, sin4, ui1, ui2, ui3, ur1, ur2, &
296 ur3, vi1, vi2, vi3, vr1, vr2, vr3
298 dimension trig(2,
ctrig_length), zin(2, mm, m), zout(2, nn, n)
303 rt2i = 0.7071067811865475_dp
312 nout2 = nout1 + after
318 zout(1, j, nout1) = r2 + r1
319 zout(2, j, nout1) = s2 + s1
320 zout(1, j, nout2) = r1 - r2
321 zout(2, j, nout2) = s1 - s2
323 DO 2000, ia = 2, after
325 IF (2*ias .EQ. after)
THEN
326 IF (isign .EQ. 1)
THEN
333 nout2 = nout1 + after
339 zout(1, j, nout1) = r1 - r2
340 zout(2, j, nout1) = s2 + s1
341 zout(1, j, nout2) = r2 + r1
342 zout(2, j, nout2) = s1 - s2
351 nout2 = nout1 + after
357 zout(1, j, nout1) = r2 + r1
358 zout(2, j, nout1) = s1 - s2
359 zout(1, j, nout2) = r1 - r2
360 zout(2, j, nout2) = s2 + s1
363 ELSE IF (4*ias .EQ. after)
THEN
364 IF (isign .EQ. 1)
THEN
371 nout2 = nout1 + after
379 zout(1, j, nout1) = r2 + r1
380 zout(2, j, nout1) = s2 + s1
381 zout(1, j, nout2) = r1 - r2
382 zout(2, j, nout2) = s1 - s2
391 nout2 = nout1 + after
399 zout(1, j, nout1) = r2 + r1
400 zout(2, j, nout1) = s2 + s1
401 zout(1, j, nout2) = r1 - r2
402 zout(2, j, nout2) = s1 - s2
405 ELSE IF (4*ias .EQ. 3*after)
THEN
406 IF (isign .EQ. 1)
THEN
413 nout2 = nout1 + after
421 zout(1, j, nout1) = r1 - r2
422 zout(2, j, nout1) = s2 + s1
423 zout(1, j, nout2) = r2 + r1
424 zout(2, j, nout2) = s1 - s2
433 nout2 = nout1 + after
441 zout(1, j, nout1) = r2 + r1
442 zout(2, j, nout1) = s1 - s2
443 zout(1, j, nout2) = r1 - r2
444 zout(2, j, nout2) = s2 + s1
448 itrig = ias*before + 1
457 nout2 = nout1 + after
465 zout(1, j, nout1) = r2 + r1
466 zout(2, j, nout1) = s2 + s1
467 zout(1, j, nout2) = r1 - r2
468 zout(2, j, nout2) = s1 - s2
472 ELSE IF (now .EQ. 4)
THEN
473 IF (isign .EQ. 1)
THEN
483 nout2 = nout1 + after
484 nout3 = nout2 + after
485 nout4 = nout3 + after
497 zout(1, j, nout1) = r + s
498 zout(1, j, nout3) = r - s
501 zout(1, j, nout2) = r - s
502 zout(1, j, nout4) = r + s
505 zout(2, j, nout1) = r + s
506 zout(2, j, nout3) = r - s
509 zout(2, j, nout2) = r + s
510 zout(2, j, nout4) = r - s
512 DO 4000, ia = 2, after
514 IF (2*ias .EQ. after)
THEN
523 nout2 = nout1 + after
524 nout3 = nout2 + after
525 nout4 = nout3 + after
541 zout(1, j, nout1) = r + s
542 zout(1, j, nout3) = r - s
545 zout(1, j, nout2) = r - s
546 zout(1, j, nout4) = r + s
549 zout(2, j, nout1) = r + s
550 zout(2, j, nout3) = r - s
553 zout(2, j, nout2) = r + s
554 zout(2, j, nout4) = r - s
575 nout2 = nout1 + after
576 nout3 = nout2 + after
577 nout4 = nout3 + after
595 zout(1, j, nout1) = r + s
596 zout(1, j, nout3) = r - s
599 zout(1, j, nout2) = r - s
600 zout(1, j, nout4) = r + s
603 zout(2, j, nout1) = r + s
604 zout(2, j, nout3) = r - s
607 zout(2, j, nout2) = r + s
608 zout(2, j, nout4) = r - s
622 nout2 = nout1 + after
623 nout3 = nout2 + after
624 nout4 = nout3 + after
636 zout(1, j, nout1) = r + s
637 zout(1, j, nout3) = r - s
640 zout(1, j, nout2) = r + s
641 zout(1, j, nout4) = r - s
644 zout(2, j, nout1) = r + s
645 zout(2, j, nout3) = r - s
648 zout(2, j, nout2) = r - s
649 zout(2, j, nout4) = r + s
651 DO 4100, ia = 2, after
653 IF (2*ias .EQ. after)
THEN
662 nout2 = nout1 + after
663 nout3 = nout2 + after
664 nout4 = nout3 + after
680 zout(1, j, nout1) = r + s
681 zout(1, j, nout3) = r - s
684 zout(1, j, nout2) = r + s
685 zout(1, j, nout4) = r - s
688 zout(2, j, nout1) = r + s
689 zout(2, j, nout3) = r - s
692 zout(2, j, nout2) = r - s
693 zout(2, j, nout4) = r + s
714 nout2 = nout1 + after
715 nout3 = nout2 + after
716 nout4 = nout3 + after
734 zout(1, j, nout1) = r + s
735 zout(1, j, nout3) = r - s
738 zout(1, j, nout2) = r + s
739 zout(1, j, nout4) = r - s
742 zout(2, j, nout1) = r + s
743 zout(2, j, nout3) = r - s
746 zout(2, j, nout2) = r - s
747 zout(2, j, nout4) = r + s
752 ELSE IF (now .EQ. 8)
THEN
753 IF (isign .EQ. -1)
THEN
767 nout2 = nout1 + after
768 nout3 = nout2 + after
769 nout4 = nout3 + after
770 nout5 = nout4 + after
771 nout6 = nout5 + after
772 nout7 = nout6 + after
773 nout8 = nout7 + after
807 zout(1, j, nout1) = ap + bp
808 zout(2, j, nout1) = cp + dpp
809 zout(1, j, nout5) = ap - bp
810 zout(2, j, nout5) = cp - dpp
811 zout(1, j, nout3) = am + dm
812 zout(2, j, nout3) = cm - bm
813 zout(1, j, nout7) = am - dm
814 zout(2, j, nout7) = cm + bm
834 dpp = (cm - dpp)*rt2i
835 zout(1, j, nout2) = ap + r
836 zout(2, j, nout2) = bm + s
837 zout(1, j, nout6) = ap - r
838 zout(2, j, nout6) = bm - s
839 zout(1, j, nout4) = am + cp
840 zout(2, j, nout4) = bp + dpp
841 zout(1, j, nout8) = am - cp
842 zout(2, j, nout8) = bp - dpp
844 DO 8000, ia = 2, after
880 nout2 = nout1 + after
881 nout3 = nout2 + after
882 nout4 = nout3 + after
883 nout5 = nout4 + after
884 nout6 = nout5 + after
885 nout7 = nout6 + after
886 nout8 = nout7 + after
934 zout(1, j, nout1) = ap + bp
935 zout(2, j, nout1) = cp + dpp
936 zout(1, j, nout5) = ap - bp
937 zout(2, j, nout5) = cp - dpp
938 zout(1, j, nout3) = am + dm
939 zout(2, j, nout3) = cm - bm
940 zout(1, j, nout7) = am - dm
941 zout(2, j, nout7) = cm + bm
961 dpp = (cm - dpp)*rt2i
962 zout(1, j, nout2) = ap + r
963 zout(2, j, nout2) = bm + s
964 zout(1, j, nout6) = ap - r
965 zout(2, j, nout6) = bm - s
966 zout(1, j, nout4) = am + cp
967 zout(2, j, nout4) = bp + dpp
968 zout(1, j, nout8) = am - cp
969 zout(2, j, nout8) = bp - dpp
987 nout2 = nout1 + after
988 nout3 = nout2 + after
989 nout4 = nout3 + after
990 nout5 = nout4 + after
991 nout6 = nout5 + after
992 nout7 = nout6 + after
993 nout8 = nout7 + after
1000 s3 = zin(2, j, nin3)
1001 r4 = zin(1, j, nin4)
1002 s4 = zin(2, j, nin4)
1003 r5 = zin(1, j, nin5)
1004 s5 = zin(2, j, nin5)
1005 r6 = zin(1, j, nin6)
1006 s6 = zin(2, j, nin6)
1007 r7 = zin(1, j, nin7)
1008 s7 = zin(2, j, nin7)
1009 r8 = zin(1, j, nin8)
1010 s8 = zin(2, j, nin8)
1027 zout(1, j, nout1) = ap + bp
1028 zout(2, j, nout1) = cp + dpp
1029 zout(1, j, nout5) = ap - bp
1030 zout(2, j, nout5) = cp - dpp
1031 zout(1, j, nout3) = am - dm
1032 zout(2, j, nout3) = cm + bm
1033 zout(1, j, nout7) = am + dm
1034 zout(2, j, nout7) = cm - bm
1053 cp = (cm + dpp)*rt2i
1054 dpp = (dpp - cm)*rt2i
1055 zout(1, j, nout2) = ap + r
1056 zout(2, j, nout2) = bm + s
1057 zout(1, j, nout6) = ap - r
1058 zout(2, j, nout6) = bm - s
1059 zout(1, j, nout4) = am + cp
1060 zout(2, j, nout4) = bp + dpp
1061 zout(1, j, nout8) = am - cp
1062 zout(2, j, nout8) = bp - dpp
1065 DO 8001, ia = 2, after
1069 cr2 = trig(1, itrig)
1070 ci2 = trig(2, itrig)
1072 cr3 = trig(1, itrig)
1073 ci3 = trig(2, itrig)
1075 cr4 = trig(1, itrig)
1076 ci4 = trig(2, itrig)
1078 cr5 = trig(1, itrig)
1079 ci5 = trig(2, itrig)
1081 cr6 = trig(1, itrig)
1082 ci6 = trig(2, itrig)
1084 cr7 = trig(1, itrig)
1085 ci7 = trig(2, itrig)
1087 cr8 = trig(1, itrig)
1088 ci8 = trig(2, itrig)
1101 nout2 = nout1 + after
1102 nout3 = nout2 + after
1103 nout4 = nout3 + after
1104 nout5 = nout4 + after
1105 nout6 = nout5 + after
1106 nout7 = nout6 + after
1107 nout8 = nout7 + after
1109 r1 = zin(1, j, nin1)
1110 s1 = zin(2, j, nin1)
1155 zout(1, j, nout1) = ap + bp
1156 zout(2, j, nout1) = cp + dpp
1157 zout(1, j, nout5) = ap - bp
1158 zout(2, j, nout5) = cp - dpp
1159 zout(1, j, nout3) = am - dm
1160 zout(2, j, nout3) = cm + bm
1161 zout(1, j, nout7) = am + dm
1162 zout(2, j, nout7) = cm - bm
1181 cp = (cm + dpp)*rt2i
1182 dpp = (dpp - cm)*rt2i
1183 zout(1, j, nout2) = ap + r
1184 zout(2, j, nout2) = bm + s
1185 zout(1, j, nout6) = ap - r
1186 zout(2, j, nout6) = bm - s
1187 zout(1, j, nout4) = am + cp
1188 zout(2, j, nout4) = bp + dpp
1189 zout(1, j, nout8) = am - cp
1190 zout(2, j, nout8) = bp - dpp
1195 ELSE IF (now .EQ. 3)
THEN
1197 bb = isign*0.8660254037844387_dp
1206 nout2 = nout1 + after
1207 nout3 = nout2 + after
1209 r1 = zin(1, j, nin1)
1210 s1 = zin(2, j, nin1)
1211 r2 = zin(1, j, nin2)
1212 s2 = zin(2, j, nin2)
1213 r3 = zin(1, j, nin3)
1214 s3 = zin(2, j, nin3)
1217 zout(1, j, nout1) = r + r1
1218 zout(2, j, nout1) = s + s1
1223 zout(1, j, nout2) = r1 - s2
1224 zout(2, j, nout2) = s1 + r2
1225 zout(1, j, nout3) = r1 + s2
1226 zout(2, j, nout3) = s1 - r2
1228 DO 3000, ia = 2, after
1230 IF (4*ias .EQ. 3*after)
THEN
1231 IF (isign .EQ. 1)
THEN
1239 nout2 = nout1 + after
1240 nout3 = nout2 + after
1242 r1 = zin(1, j, nin1)
1243 s1 = zin(2, j, nin1)
1244 r2 = zin(2, j, nin2)
1245 s2 = zin(1, j, nin2)
1246 r3 = zin(1, j, nin3)
1247 s3 = zin(2, j, nin3)
1250 zout(1, j, nout1) = r1 - r
1251 zout(2, j, nout1) = s + s1
1256 zout(1, j, nout2) = r1 - s2
1257 zout(2, j, nout2) = s1 - r2
1258 zout(1, j, nout3) = r1 + s2
1259 zout(2, j, nout3) = s1 + r2
1269 nout2 = nout1 + after
1270 nout3 = nout2 + after
1272 r1 = zin(1, j, nin1)
1273 s1 = zin(2, j, nin1)
1274 r2 = zin(2, j, nin2)
1275 s2 = zin(1, j, nin2)
1276 r3 = zin(1, j, nin3)
1277 s3 = zin(2, j, nin3)
1280 zout(1, j, nout1) = r + r1
1281 zout(2, j, nout1) = s1 - s
1286 zout(1, j, nout2) = r1 + s2
1287 zout(2, j, nout2) = s1 + r2
1288 zout(1, j, nout3) = r1 - s2
1289 zout(2, j, nout3) = s1 - r2
1292 ELSE IF (8*ias .EQ. 3*after)
THEN
1293 IF (isign .EQ. 1)
THEN
1301 nout2 = nout1 + after
1302 nout3 = nout2 + after
1304 r1 = zin(1, j, nin1)
1305 s1 = zin(2, j, nin1)
1310 r3 = zin(2, j, nin3)
1311 s3 = zin(1, j, nin3)
1314 zout(1, j, nout1) = r + r1
1315 zout(2, j, nout1) = s + s1
1320 zout(1, j, nout2) = r1 - s2
1321 zout(2, j, nout2) = s1 + r2
1322 zout(1, j, nout3) = r1 + s2
1323 zout(2, j, nout3) = s1 - r2
1333 nout2 = nout1 + after
1334 nout3 = nout2 + after
1336 r1 = zin(1, j, nin1)
1337 s1 = zin(2, j, nin1)
1342 r3 = zin(2, j, nin3)
1343 s3 = zin(1, j, nin3)
1346 zout(1, j, nout1) = r + r1
1347 zout(2, j, nout1) = s + s1
1352 zout(1, j, nout2) = r1 - s2
1353 zout(2, j, nout2) = s1 + r2
1354 zout(1, j, nout3) = r1 + s2
1355 zout(2, j, nout3) = s1 - r2
1361 cr2 = trig(1, itrig)
1362 ci2 = trig(2, itrig)
1364 cr3 = trig(1, itrig)
1365 ci3 = trig(2, itrig)
1373 nout2 = nout1 + after
1374 nout3 = nout2 + after
1376 r1 = zin(1, j, nin1)
1377 s1 = zin(2, j, nin1)
1388 zout(1, j, nout1) = r + r1
1389 zout(2, j, nout1) = s + s1
1394 zout(1, j, nout2) = r1 - s2
1395 zout(2, j, nout2) = s1 + r2
1396 zout(1, j, nout3) = r1 + s2
1397 zout(2, j, nout3) = s1 - r2
1401 ELSE IF (now .EQ. 5)
THEN
1403 cos2 = 0.3090169943749474_dp
1405 cos4 = -0.8090169943749474_dp
1407 sin2 = isign*0.9510565162951536_dp
1409 sin4 = isign*0.5877852522924731_dp
1420 nout2 = nout1 + after
1421 nout3 = nout2 + after
1422 nout4 = nout3 + after
1423 nout5 = nout4 + after
1425 r1 = zin(1, j, nin1)
1426 s1 = zin(2, j, nin1)
1427 r2 = zin(1, j, nin2)
1428 s2 = zin(2, j, nin2)
1429 r3 = zin(1, j, nin3)
1430 s3 = zin(2, j, nin3)
1431 r4 = zin(1, j, nin4)
1432 s4 = zin(2, j, nin4)
1433 r5 = zin(1, j, nin5)
1434 s5 = zin(2, j, nin5)
1439 zout(1, j, nout1) = r1 + r25 + r34
1440 r = r1 + cos2*r25 + cos4*r34
1441 s = sin2*s25 + sin4*s34
1442 zout(1, j, nout2) = r - s
1443 zout(1, j, nout5) = r + s
1444 r = r1 + cos4*r25 + cos2*r34
1445 s = sin4*s25 - sin2*s34
1446 zout(1, j, nout3) = r - s
1447 zout(1, j, nout4) = r + s
1452 zout(2, j, nout1) = s1 + s25 + s34
1453 r = s1 + cos2*s25 + cos4*s34
1454 s = sin2*r25 + sin4*r34
1455 zout(2, j, nout2) = r + s
1456 zout(2, j, nout5) = r - s
1457 r = s1 + cos4*s25 + cos2*s34
1458 s = sin4*r25 - sin2*r34
1459 zout(2, j, nout3) = r + s
1460 zout(2, j, nout4) = r - s
1462 DO 5000, ia = 2, after
1464 IF (8*ias .EQ. 5*after)
THEN
1465 IF (isign .EQ. 1)
THEN
1475 nout2 = nout1 + after
1476 nout3 = nout2 + after
1477 nout4 = nout3 + after
1478 nout5 = nout4 + after
1480 r1 = zin(1, j, nin1)
1481 s1 = zin(2, j, nin1)
1486 r3 = zin(2, j, nin3)
1487 s3 = zin(1, j, nin3)
1492 r5 = zin(1, j, nin5)
1493 s5 = zin(2, j, nin5)
1498 zout(1, j, nout1) = r1 + r25 - r34
1499 r = r1 + cos2*r25 - cos4*r34
1500 s = sin2*s25 + sin4*s34
1501 zout(1, j, nout2) = r - s
1502 zout(1, j, nout5) = r + s
1503 r = r1 + cos4*r25 - cos2*r34
1504 s = sin4*s25 - sin2*s34
1505 zout(1, j, nout3) = r - s
1506 zout(1, j, nout4) = r + s
1511 zout(2, j, nout1) = s1 + s25 + s34
1512 r = s1 + cos2*s25 + cos4*s34
1513 s = sin2*r25 + sin4*r34
1514 zout(2, j, nout2) = r + s
1515 zout(2, j, nout5) = r - s
1516 r = s1 + cos4*s25 + cos2*s34
1517 s = sin4*r25 - sin2*r34
1518 zout(2, j, nout3) = r + s
1519 zout(2, j, nout4) = r - s
1531 nout2 = nout1 + after
1532 nout3 = nout2 + after
1533 nout4 = nout3 + after
1534 nout5 = nout4 + after
1536 r1 = zin(1, j, nin1)
1537 s1 = zin(2, j, nin1)
1542 r3 = zin(2, j, nin3)
1543 s3 = zin(1, j, nin3)
1548 r5 = zin(1, j, nin5)
1549 s5 = zin(2, j, nin5)
1554 zout(1, j, nout1) = r1 + r25 + r34
1555 r = r1 + cos2*r25 + cos4*r34
1556 s = sin2*s25 + sin4*s34
1557 zout(1, j, nout2) = r - s
1558 zout(1, j, nout5) = r + s
1559 r = r1 + cos4*r25 + cos2*r34
1560 s = sin4*s25 - sin2*s34
1561 zout(1, j, nout3) = r - s
1562 zout(1, j, nout4) = r + s
1567 zout(2, j, nout1) = s1 + s25 - s34
1568 r = s1 + cos2*s25 - cos4*s34
1569 s = sin2*r25 + sin4*r34
1570 zout(2, j, nout2) = r + s
1571 zout(2, j, nout5) = r - s
1572 r = s1 + cos4*s25 - cos2*s34
1573 s = sin4*r25 - sin2*r34
1574 zout(2, j, nout3) = r + s
1575 zout(2, j, nout4) = r - s
1582 cr2 = trig(1, itrig)
1583 ci2 = trig(2, itrig)
1585 cr3 = trig(1, itrig)
1586 ci3 = trig(2, itrig)
1588 cr4 = trig(1, itrig)
1589 ci4 = trig(2, itrig)
1591 cr5 = trig(1, itrig)
1592 ci5 = trig(2, itrig)
1602 nout2 = nout1 + after
1603 nout3 = nout2 + after
1604 nout4 = nout3 + after
1605 nout5 = nout4 + after
1607 r1 = zin(1, j, nin1)
1608 s1 = zin(2, j, nin1)
1629 zout(1, j, nout1) = r1 + r25 + r34
1630 r = r1 + cos2*r25 + cos4*r34
1631 s = sin2*s25 + sin4*s34
1632 zout(1, j, nout2) = r - s
1633 zout(1, j, nout5) = r + s
1634 r = r1 + cos4*r25 + cos2*r34
1635 s = sin4*s25 - sin2*s34
1636 zout(1, j, nout3) = r - s
1637 zout(1, j, nout4) = r + s
1642 zout(2, j, nout1) = s1 + s25 + s34
1643 r = s1 + cos2*s25 + cos4*s34
1644 s = sin2*r25 + sin4*r34
1645 zout(2, j, nout2) = r + s
1646 zout(2, j, nout5) = r - s
1647 r = s1 + cos4*s25 + cos2*s34
1648 s = sin4*r25 - sin2*r34
1649 zout(2, j, nout3) = r + s
1650 zout(2, j, nout4) = r - s
1654 ELSE IF (now .EQ. 6)
THEN
1656 bb = isign*0.8660254037844387_dp
1669 nout2 = nout1 + after
1670 nout3 = nout2 + after
1671 nout4 = nout3 + after
1672 nout5 = nout4 + after
1673 nout6 = nout5 + after
1675 r2 = zin(1, j, nin3)
1676 s2 = zin(2, j, nin3)
1677 r3 = zin(1, j, nin5)
1678 s3 = zin(2, j, nin5)
1681 r1 = zin(1, j, nin1)
1682 s1 = zin(2, j, nin1)
1694 r2 = zin(1, j, nin6)
1695 s2 = zin(2, j, nin6)
1696 r3 = zin(1, j, nin2)
1697 s3 = zin(2, j, nin2)
1700 r1 = zin(1, j, nin4)
1701 s1 = zin(2, j, nin4)
1713 zout(1, j, nout1) = ur1 + vr1
1714 zout(2, j, nout1) = ui1 + vi1
1715 zout(1, j, nout5) = ur2 + vr2
1716 zout(2, j, nout5) = ui2 + vi2
1717 zout(1, j, nout3) = ur3 + vr3
1718 zout(2, j, nout3) = ui3 + vi3
1719 zout(1, j, nout4) = ur1 - vr1
1720 zout(2, j, nout4) = ui1 - vi1
1721 zout(1, j, nout2) = ur2 - vr2
1722 zout(2, j, nout2) = ui2 - vi2
1723 zout(1, j, nout6) = ur3 - vr3
1724 zout(2, j, nout6) = ui3 - vi3
1727 cpabort(
"error fftstp")
pure real(kind=dp) function angle(a, b)
Calculation of the angle between the vectors a and b. The angle is returned in radians.
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public ctrig_length
subroutine, public fftstp(mm, nfft, m, nn, n, zin, zout, trig, after, now, before, isign)
...
subroutine, public ctrig(n, trig, after, before, now, isign, ic)
...
subroutine, public fourier_dim(n, n_next)
Give a number n_next > n compatible for the FFT.