00001 #include <stdlib.h>
00002 #include "util.h"
00003 #include "nfft3.h"
00004
00014 void construct(char * file, int N, int M, int Z, fftw_complex *mem)
00015 {
00016 int j,z;
00017 double tmp;
00018 nfft_plan my_plan;
00019 FILE* fp;
00020
00021
00022 nfft_init_2d(&my_plan,N,N,M/Z);
00023
00024 fp=fopen("knots.dat","r");
00025
00026 for(j=0;j<my_plan.M_total;j++)
00027 {
00028 fscanf(fp,"%le %le %le",&my_plan.x[2*j+0],&my_plan.x[2*j+1],&tmp);
00029 }
00030 fclose(fp);
00031
00032 fp=fopen(file,"w");
00033
00034 for(z=0;z<Z;z++) {
00035 tmp = (double) z;
00036
00037 for(j=0;j<N*N;j++)
00038 my_plan.f_hat[j] = mem[(z*N*N+N*N*Z/2+j)%(N*N*Z)];
00039
00040 if(my_plan.nfft_flags & PRE_PSI)
00041 nfft_precompute_psi(&my_plan);
00042
00043 nfft_trafo(&my_plan);
00044
00045 for(j=0;j<my_plan.M_total;j++)
00046 {
00047 fprintf(fp,"%le %le %le %le %le\n",my_plan.x[2*j+0],my_plan.x[2*j+1],tmp/Z-0.5,
00048 creal(my_plan.f[j]),cimag(my_plan.f[j]));
00049 }
00050 }
00051 fclose(fp);
00052
00053 nfft_finalize(&my_plan);
00054 }
00055
00060 void fft(int N,int M,int Z, fftw_complex *mem)
00061 {
00062 fftw_plan plan;
00063 plan = fftw_plan_many_dft(1, &Z, N*N,
00064 mem, NULL,
00065 N*N, 1,
00066 mem, NULL,
00067 N*N,1 ,
00068 FFTW_FORWARD, FFTW_ESTIMATE);
00069
00070 fftw_execute(plan);
00071 fftw_destroy_plan(plan);
00072 }
00073
00078 void read_data(int N,int M,int Z, fftw_complex *mem)
00079 {
00080 int i,z;
00081 double real;
00082 FILE* fin;
00083 fin=fopen("input_f.dat","r");
00084
00085 for(z=0;z<Z;z++)
00086 {
00087 for(i=0;i<N*N;i++)
00088 {
00089 fscanf(fin,"%le ",&real );
00090 mem[(z*N*N+N*N*Z/2+i)%(N*N*Z)]=real;
00091 }
00092 }
00093 fclose(fin);
00094 }
00095
00096 int main(int argc, char **argv)
00097 {
00098 fftw_complex *mem;
00099
00100 if (argc <= 4) {
00101 printf("usage: ./construct_data FILENAME N M Z\n");
00102 return 1;
00103 }
00104
00105 mem = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * atoi(argv[2]) * atoi(argv[2]) * atoi(argv[4]));
00106
00107 read_data(atoi(argv[2]),atoi(argv[3]),atoi(argv[4]), mem);
00108
00109 fft(atoi(argv[2]),atoi(argv[3]),atoi(argv[4]), mem);
00110
00111 construct(argv[1],atoi(argv[2]),atoi(argv[3]),atoi(argv[4]), mem);
00112
00113 fftw_free(mem);
00114
00115 return 1;
00116 }
00117