00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #include "config.h"
00021
00022 #include <stdlib.h>
00023 #include <math.h>
00024 #ifdef HAVE_COMPLEX_H
00025 #include <complex.h>
00026 #endif
00027
00028 #include "nfft3util.h"
00029 #include "nfft3.h"
00030
00037 static void construct(char * file, int N, int M, int Z)
00038 {
00039 int j,k,l;
00040 double real;
00041 nfft_plan my_plan;
00042 FILE* fp,*fk;
00043 int my_N[3],my_n[3];
00044
00045
00046
00047
00048 my_N[0]=Z; my_n[0]=ceil(Z*1.2);
00049 my_N[1]=N; my_n[1]=ceil(N*1.2);
00050 my_N[2]=N; my_n[2]=ceil(N*1.2);
00051 nfft_init_guru(&my_plan, 3, my_N, M, my_n, 6,
00052 PRE_PHI_HUT| PRE_PSI |MALLOC_X| MALLOC_F_HAT|
00053 MALLOC_F| FFTW_INIT| FFT_OUT_OF_PLACE,
00054 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00055
00056 fp=fopen("knots.dat","r");
00057
00058 for(j=0;j<M;j++)
00059 fscanf(fp,"%le %le %le",&my_plan.x[3*(j)+1],
00060 &my_plan.x[3*(j)+2],&my_plan.x[3*(j)+0]);
00061
00062 fclose(fp);
00063
00064 fp=fopen("input_f.dat","r");
00065 fk=fopen(file,"w");
00066
00067 for(l=0;l<Z;l++) {
00068 for(j=0;j<N;j++)
00069 {
00070 for(k=0;k<N;k++)
00071 {
00072
00073 fscanf(fp,"%le ",&real);
00074 my_plan.f_hat[(N*N*l+N*j+k)] = real;
00075 }
00076 }
00077 }
00078
00079 if(my_plan.nfft_flags & PRE_PSI)
00080 nfft_precompute_psi(&my_plan);
00081
00082 nfft_trafo(&my_plan);
00083
00084
00085 for(j=0;j<my_plan.M_total;j++)
00086 fprintf(fk,"%le %le %le %le %le\n",my_plan.x[3*j+1],
00087 my_plan.x[3*j+2],my_plan.x[3*j+0],creal(my_plan.f[j]),cimag(my_plan.f[j]));
00088
00089
00090
00091 fclose(fk);
00092 fclose(fp);
00093
00094 nfft_finalize(&my_plan);
00095 }
00096
00097 int main(int argc, char **argv)
00098 {
00099 if (argc <= 4) {
00100 printf("usage: ./construct_data FILENAME N M Z\n");
00101 return 1;
00102 }
00103
00104 construct(argv[1], atoi(argv[2]),atoi(argv[3]),atoi(argv[4]));
00105
00106 return 1;
00107 }
00108