NFFT Logo 3.0 API Reference
Main Page | Modules | Data Structures | Directories | File List | Data Fields | Globals

construct_data_2d1d.c

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;                /* some variables */
00017   double tmp;             /* a placeholder */
00018   nfft_plan my_plan;      /* plan for the two dimensional nfft  */
00019   FILE* fp;
00020     
00021   /* initialise my_plan */
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); /* execute the fft */
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 /* \} */

Generated on 1 Nov 2006 by Doxygen 1.4.4