LMU/musrfit
LMU/musrfit
5
0
Fork 0
Code Issues 6 Pull Requests Actions Packages Projects Releases 3 Activity

Merged muonspin/musrfit:root6 into master

Browse Source
This commit is contained in:
Zaher Salman 2019-10-28 09:21:05 +01:00
commit 4ef4316b97
1 changed files with 68 additions and 37 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

105
src/tests/fftw3_bare/fftw3_test.cpp
View File

@ -3,41 +3,72 @@
#include <fftw3.h> #include <fftw3.h>
int main(int argc, char *argv[]) { void fftw3_test_syntax() {
const unsigned int N=8192; std::cout << std::endl;
std::cout << "fftw3_test [N [offset]]" << std::endl;
fftw_complex *in, *out; std::cout << " N: 2^N Fourier Power" << std::endl;
fftw_plan my_plan; std::cout << " default value if not given: N=13 (2^13 = 8192)" << std::endl;
std::cout << " offset: allows to have a length != 2^N, namely 2^N + offset" << std::endl;
in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*N); std::cout << " default value of offset if not given is 0." << std::endl;
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*N); std::cout << std::endl;
}
// feed input
double x=0.0; int main(int argc, char *argv[]) {
for (unsigned int i=0; i<N; i++) { unsigned int N=8192;
x = (double)i / (double)N; unsigned int offset=0;
in[i][0] = sin(100.0*x)*exp(-x/8.0);
in[i][1] = 0.0; if (argc > 3) {
//as35 if ((i%100)==0) std::cout << i << ": x=" << x << ", in=" << in[i][0] << std::endl; fftw3_test_syntax();
} return 1;
}
my_plan = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
int a=-1;
fftw_execute(my_plan); for (int i=1; i<argc; i++) {
int status = sscanf(argv[i], "%d", &a);
double sumRe = 0.0, sumIm = 0.0; if (status != 1) {
for (unsigned int i=0; i<N/2; i++) { fftw3_test_syntax();
//as35 if ((i % 10) == 0) std::cout << i << ": re=" << out[i][0] << ", im=" << out[i][1] << std::endl; return 2;
sumRe += out[i][0]; }
sumIm += out[i][1]; if (i==1) // N
} N = (unsigned int)pow(2.0, a);
if (i==2) // offset
std::cout << "sumRe=" << sumRe << std::endl; offset = (unsigned int)a;
std::cout << "sumIm=" << sumIm << std::endl; }
std::cout << "debug> N=" << N << ", offset=" << offset << std::endl;
fftw_destroy_plan(my_plan);
fftw_free(in); fftw_complex *in, *out;
fftw_free(out); fftw_plan my_plan;
return 0; N += offset;
in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*N);
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*N);
// feed input
double x=0.0;
for (unsigned int i=0; i<N; i++) {
x = (double)i / (double)N;
in[i][0] = sin(100.0*x)*exp(-x/8.0);
in[i][1] = 0.0;
//as35 if ((i%100)==0) std::cout << i << ": x=" << x << ", in=" << in[i][0] << std::endl;
}
my_plan = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(my_plan);
double sumRe = 0.0, sumIm = 0.0;
for (unsigned int i=0; i<N/2; i++) {
//as35 if ((i % 10) == 0) std::cout << i << ": re=" << out[i][0] << ", im=" << out[i][1] << std::endl;
sumRe += out[i][0];
sumIm += out[i][1];
}
std::cout << "sumRe=" << sumRe << std::endl;
std::cout << "sumIm=" << sumIm << std::endl;
fftw_destroy_plan(my_plan);
fftw_free(in);
fftw_free(out);
return 0;
} }