3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
9 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
10 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
11 * Copyright (c) 2001-2008, The GROMACS development team,
12 * check out http://www.gromacs.org for more information.
14 * This program is free software; you can redistribute it and/or
15 * modify it under the terms of the GNU General Public License
16 * as published by the Free Software Foundation; either version 2
17 * of the License, or (at your option) any later version.
19 * If you want to redistribute modifications, please consider that
20 * scientific software is very special. Version control is crucial -
21 * bugs must be traceable. We will be happy to consider code for
22 * inclusion in the official distribution, but derived work must not
23 * be called official GROMACS. Details are found in the README & COPYING
24 * files - if they are missing, get the official version at www.gromacs.org.
26 * To help us fund GROMACS development, we humbly ask that you cite
27 * the papers on the package - you can find them in the top README file.
29 * For more info, check our website at http://www.gromacs.org
32 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
46 #define SWAP(a, b) tempr = (a); (a) = (b); (b) = tempr
48 void four1(real data[], int nn, int isign)
50 int n, mmax, m, j, istep, i;
51 double wtemp, wr, wpr, wpi, wi, theta;
56 for (i = 1; i < n; i += 2)
60 SWAP(data[j], data[i]);
61 SWAP(data[j+1], data[i+1]);
64 while (m >= 2 && j > m)
75 theta = 6.28318530717959/(isign*mmax);
76 wtemp = sin(0.5*theta);
77 wpr = -2.0*wtemp*wtemp;
81 for (m = 1; m < mmax; m += 2)
83 for (i = m; i <= n; i += istep)
86 tempr = wr*data[j]-wi*data[j+1];
87 tempi = wr*data[j+1]+wi*data[j];
88 data[j] = data[i]-tempr;
89 data[j+1] = data[i+1]-tempi;
93 wr = (wtemp = wr)*wpr-wi*wpi+wr;
94 wi = wi*wpr+wtemp*wpi+wi;
102 static void realft(real data[], int n, int isign)
104 int i, i1, i2, i3, i4, n2p3;
105 real c1 = 0.5, c2, h1r, h1i, h2r, h2i;
106 double wr, wi, wpr, wpi, wtemp, theta;
108 theta = 3.141592653589793/(double) n;
119 wtemp = sin(0.5*theta);
120 wpr = -2.0*wtemp*wtemp;
125 for (i = 2; i <= n/2; i++)
127 i4 = 1+(i3 = n2p3-(i2 = 1+(i1 = i+i-1)));
128 h1r = c1*(data[i1]+data[i3]);
129 h1i = c1*(data[i2]-data[i4]);
130 h2r = -c2*(data[i2]+data[i4]);
131 h2i = c2*(data[i1]-data[i3]);
132 data[i1] = h1r+wr*h2r-wi*h2i;
133 data[i2] = h1i+wr*h2i+wi*h2r;
134 data[i3] = h1r-wr*h2r+wi*h2i;
135 data[i4] = -h1i+wr*h2i+wi*h2r;
136 wr = (wtemp = wr)*wpr-wi*wpi+wr;
137 wi = wi*wpr+wtemp*wpi+wi;
141 data[1] = (h1r = data[1])+data[2];
142 data[2] = h1r-data[2];
146 data[1] = c1*((h1r = data[1])+data[2]);
147 data[2] = c1*(h1r-data[2]);
152 static void twofft(real data1[], real data2[], real fft1[], real fft2[], int n)
155 real rep, rem, aip, aim;
157 nn3 = 1+(nn2 = 2+n+n);
158 for (j = 1, jj = 2; j <= n; j++, jj += 2)
160 fft1[jj-1] = data1[j];
165 fft1[2] = fft2[2] = 0.0;
166 for (j = 3; j <= n+1; j += 2)
168 rep = 0.5*(fft1[j]+fft1[nn2-j]);
169 rem = 0.5*(fft1[j]-fft1[nn2-j]);
170 aip = 0.5*(fft1[j+1]+fft1[nn3-j]);
171 aim = 0.5*(fft1[j+1]-fft1[nn3-j]);
183 void correl(real data1[], real data2[], int n, real ans[])
189 twofft(data1, data2, fft, ans, n);
191 for (i = 2; i <= n+2; i += 2)
194 ans[i-1] = (fft[i-1]*dum+fft[i]*ans[i])/no2;
195 ans[i] = (fft[i]*dum-fft[i-1]*ans[i])/no2;
198 realft(ans, no2, -1);