36 #define SWAP(a,b) tempr=(a); (a)=(b); (b)=tempr
37 #define TWOPI 6.28318530717959
51 unsigned int dimCount = nn[idim];
52 if (!dimCount || (dimCount & (dimCount - 1)))
56 "fft::transform(complexField&, const labelList&, "
58 ) <<
"number of elements in direction " << idim
59 <<
" is not a power of 2" <<
endl
60 <<
" Number of elements in each direction = " << nn
65 const label ndim = nn.
size();
67 label i1, i2, i3, i2rev, i3rev, ip1, ip2, ip3, ifp1, ifp2;
68 label ibit, k1, k2, n, nprev, nrem, idim;
70 scalar theta, wi, wpi, wpr, wr, wtemp;
71 scalar* data =
reinterpret_cast<scalar*
>(field.
begin()) - 1;
76 if (isign == REVERSE_TRANSFORM)
91 for (idim=ndim; idim>=1; idim--)
94 nrem = ntot/(n*nprev);
100 for (i2=1; i2<=ip2; i2+=ip1)
104 for (i1=i2; i1<=i2 + ip1 - 2; i1+=2)
106 for (i3=i1; i3<=ip3; i3+=ip2)
108 i3rev = i2rev + i3 - i2;
109 SWAP(data[i3], data[i3rev]);
110 SWAP(data[i3 + 1], data[i3rev + 1]);
116 while (ibit >= ip1 && i2rev > ibit)
130 theta = isign*
TWOPI/(ifp2/ip1);
131 wtemp =
sin(0.5*theta);
132 wpr = -2.0*wtemp*wtemp;
137 for (i3 = 1; i3 <= ifp1; i3 += ip1)
139 for (i1 = i3; i1 <= i3 + ip1 - 2; i1 += 2)
141 for (i2 = i1; i2 <= ip3; i2 += ifp2)
145 tempr = scalar(wr*data[k2]) - scalar(wi*data[k2 + 1]);
146 tempi = scalar(wr*data[k2 + 1]) + scalar(wi*data[k2]);
147 data[k2] = data[k1] - tempr;
148 data[k2 + 1] = data[k1 + 1] - tempi;
150 data[k1 + 1] += tempi;
154 wr = (wtemp = wr)*wpr - wi*wpi + wr;
155 wi = wi*wpr + wtemp*wpi + wi;
165 if (isign == FORWARD_TRANSFORM)
173 scalar recRootN = 1.0/
sqrt(scalar(ntot));
177 field[i] *= recRootN;
197 transform(tfftField(), nn, FORWARD_TRANSFORM);
213 transform(tifftField(), nn, REVERSE_TRANSFORM);
237 tfftVectorField().replace
240 forwardTransform(tfield().
component(cmpt), nn)
246 return tfftVectorField;
266 tifftVectorField().replace
269 reverseTransform(tfield().
component(cmpt), nn)
275 return tifftVectorField;