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
00034 void construct(char * file, int N, int M, int Z)
00035 {
00036 int j,k,l;
00037 double real;
00038 nfft_plan my_plan;
00039 FILE* fp,*fk;
00040 int my_N[3],my_n[3];
00041
00042
00043
00044
00045 my_N[0]=Z; my_n[0]=ceil(Z*1.2);
00046 my_N[1]=N; my_n[1]=ceil(N*1.2);
00047 my_N[2]=N; my_n[2]=ceil(N*1.2);
00048 nfft_init_guru(&my_plan, 3, my_N, M, my_n, 6,
00049 PRE_PHI_HUT| PRE_PSI |MALLOC_X| MALLOC_F_HAT|
00050 MALLOC_F| FFTW_INIT| FFT_OUT_OF_PLACE,
00051 FFTW_MEASURE| FFTW_DESTROY_INPUT);
00052
00053 fp=fopen("knots.dat","r");
00054
00055 for(j=0;j<M;j++)
00056 fscanf(fp,"%le %le %le",&my_plan.x[3*(j)+1],
00057 &my_plan.x[3*(j)+2],&my_plan.x[3*(j)+0]);
00058
00059 fclose(fp);
00060
00061 fp=fopen("input_f.dat","r");
00062 fk=fopen(file,"w");
00063
00064 for(l=0;l<Z;l++) {
00065 for(j=0;j<N;j++)
00066 {
00067 for(k=0;k<N;k++)
00068 {
00069
00070 fscanf(fp,"%le ",&real);
00071 my_plan.f_hat[(N*N*l+N*j+k)] = real;
00072 }
00073 }
00074 }
00075
00076 if(my_plan.nfft_flags & PRE_PSI)
00077 nfft_precompute_psi(&my_plan);
00078
00079 nfft_trafo(&my_plan);
00080
00081
00082 for(j=0;j<my_plan.M_total;j++)
00083 fprintf(fk,"%le %le %le %le %le\n",my_plan.x[3*j+1],
00084 my_plan.x[3*j+2],my_plan.x[3*j+0],creal(my_plan.f[j]),cimag(my_plan.f[j]));
00085
00086
00087
00088 fclose(fk);
00089 fclose(fp);
00090
00091 nfft_finalize(&my_plan);
00092 }
00093
00094 int main(int argc, char **argv)
00095 {
00096 if (argc <= 4) {
00097 printf("usage: ./construct_data FILENAME N M Z\n");
00098 return 1;
00099 }
00100
00101 construct(argv[1], atoi(argv[2]),atoi(argv[3]),atoi(argv[4]));
00102
00103 return 1;
00104 }
00105