00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include <stdlib.h>
00022 #include <math.h>
00023 #include <complex.h>
00024
00025 #include "nfft3util.h"
00026 #include "nfft3.h"
00027
00037 void construct(char * file, int N, int M, int Z, fftw_complex *mem)
00038 {
00039 int j,z;
00040 double tmp;
00041 nfft_plan my_plan;
00042 FILE* fp;
00043
00044
00045 nfft_init_2d(&my_plan,N,N,M/Z);
00046
00047 fp=fopen("knots.dat","r");
00048
00049 for(j=0;j<my_plan.M_total;j++)
00050 {
00051 fscanf(fp,"%le %le %le",&my_plan.x[2*j+0],&my_plan.x[2*j+1],&tmp);
00052 }
00053 fclose(fp);
00054
00055 fp=fopen(file,"w");
00056
00057 for(z=0;z<Z;z++) {
00058 tmp = (double) z;
00059
00060 for(j=0;j<N*N;j++)
00061 my_plan.f_hat[j] = mem[(z*N*N+N*N*Z/2+j)%(N*N*Z)];
00062
00063 if(my_plan.nfft_flags & PRE_PSI)
00064 nfft_precompute_psi(&my_plan);
00065
00066 nfft_trafo(&my_plan);
00067
00068 for(j=0;j<my_plan.M_total;j++)
00069 {
00070 fprintf(fp,"%le %le %le %le %le\n",my_plan.x[2*j+0],my_plan.x[2*j+1],tmp/Z-0.5,
00071 creal(my_plan.f[j]),cimag(my_plan.f[j]));
00072 }
00073 }
00074 fclose(fp);
00075
00076 nfft_finalize(&my_plan);
00077 }
00078
00083 void fft(int N,int M,int Z, fftw_complex *mem)
00084 {
00085 fftw_plan plan;
00086 plan = fftw_plan_many_dft(1, &Z, N*N,
00087 mem, NULL,
00088 N*N, 1,
00089 mem, NULL,
00090 N*N,1 ,
00091 FFTW_FORWARD, FFTW_ESTIMATE);
00092
00093 fftw_execute(plan);
00094 fftw_destroy_plan(plan);
00095 }
00096
00101 void read_data(int N,int M,int Z, fftw_complex *mem)
00102 {
00103 int i,z;
00104 double real;
00105 FILE* fin;
00106 fin=fopen("input_f.dat","r");
00107
00108 for(z=0;z<Z;z++)
00109 {
00110 for(i=0;i<N*N;i++)
00111 {
00112 fscanf(fin,"%le ",&real );
00113 mem[(z*N*N+N*N*Z/2+i)%(N*N*Z)]=real;
00114 }
00115 }
00116 fclose(fin);
00117 }
00118
00119 int main(int argc, char **argv)
00120 {
00121 fftw_complex *mem;
00122
00123 if (argc <= 4) {
00124 printf("usage: ./construct_data FILENAME N M Z\n");
00125 return 1;
00126 }
00127
00128 mem = (fftw_complex*) nfft_malloc(sizeof(fftw_complex) * atoi(argv[2]) * atoi(argv[2]) * atoi(argv[4]));
00129
00130 read_data(atoi(argv[2]),atoi(argv[3]),atoi(argv[4]), mem);
00131
00132 fft(atoi(argv[2]),atoi(argv[3]),atoi(argv[4]), mem);
00133
00134 construct(argv[1],atoi(argv[2]),atoi(argv[3]),atoi(argv[4]), mem);
00135
00136 nfft_free(mem);
00137
00138 return 1;
00139 }
00140