00001 #include <stdlib.h>
00002 #include <math.h>
00003 #include "util.h"
00004 #include "nfft3.h"
00005
00012 void construct(char * file, int N, int M, int Z)
00013 {
00014 int j,k,l;
00015 double real;
00016 nfft_plan my_plan;
00017 FILE* fp,*fk;
00018 int my_N[3],my_n[3];
00019
00020
00021
00022
00023 my_N[0]=Z; my_n[0]=ceil(Z*1.2);
00024 my_N[1]=N; my_n[1]=ceil(N*1.2);
00025 my_N[2]=N; my_n[2]=ceil(N*1.2);
00026 nfft_init_guru(&my_plan, 3, my_N, M, my_n, 6,
00027 PRE_PHI_HUT| PRE_PSI |MALLOC_X| MALLOC_F_HAT|
00028 MALLOC_F| FFTW_INIT| FFT_OUT_OF_PLACE,
00029 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00030
00031 fp=fopen("knots.dat","r");
00032
00033 for(j=0;j<M;j++)
00034 fscanf(fp,"%le %le %le",&my_plan.x[3*(j)+1],
00035 &my_plan.x[3*(j)+2],&my_plan.x[3*(j)+0]);
00036
00037 fclose(fp);
00038
00039 fp=fopen("input_f.dat","r");
00040 fk=fopen(file,"w");
00041
00042 for(l=0;l<Z;l++) {
00043 for(j=0;j<N;j++)
00044 {
00045 for(k=0;k<N;k++)
00046 {
00047
00048 fscanf(fp,"%le ",&real);
00049 my_plan.f_hat[(N*N*l+N*j+k)] = real;
00050 }
00051 }
00052 }
00053
00054 if(my_plan.nfft_flags & PRE_PSI)
00055 nfft_precompute_psi(&my_plan);
00056
00057 nfft_trafo(&my_plan);
00058
00059
00060 for(j=0;j<my_plan.M_total;j++)
00061 fprintf(fk,"%le %le %le %le %le\n",my_plan.x[3*j+1],
00062 my_plan.x[3*j+2],my_plan.x[3*j+0],creal(my_plan.f[j]),cimag(my_plan.f[j]));
00063
00064
00065
00066 fclose(fk);
00067 fclose(fp);
00068
00069 nfft_finalize(&my_plan);
00070 }
00071
00072 int main(int argc, char **argv)
00073 {
00074 if (argc <= 4) {
00075 printf("usage: ./construct_data FILENAME N M Z\n");
00076 return 1;
00077 }
00078
00079 construct(argv[1], atoi(argv[2]),atoi(argv[3]),atoi(argv[4]));
00080
00081 return 1;
00082 }
00083