00001
00007
00008 #include <math.h>
00009 #include <stdlib.h>
00010 #include <string.h>
00011 #include <string.h>
00012
00013
00014 #include "nfft3.h"
00015
00016
00017 #include "legendre.h"
00018
00019
00020 #include "api.h"
00021
00022
00023 #include "util.h"
00024
00034 #define NFSFT_DEFAULT_NFFT_CUTOFF 6
00035
00041 #define NFSFT_DEFAULT_THRESHOLD 1000
00042
00048 #define NFSFT_BREAK_EVEN 5
00049
00056 static struct nfsft_wisdom wisdom = {false,0U};
00057
00080 inline void c2e(nfsft_plan *plan)
00081 {
00082 int k;
00083 int n;
00084 double complex last;
00085 double complex act;
00086 double complex *xp;
00087 double complex *xm;
00088 int low;
00089 int up;
00090 int lowe;
00091 int upe;
00093
00094 memset(plan->f_hat_intern,0U,(2*plan->N+2)*sizeof(double complex));
00095
00096
00097 lowe = -plan->N + (plan->N%2);
00098 upe = -lowe;
00099
00100
00101 for (n = lowe; n <= upe; n += 2)
00102 {
00103
00104
00105 xm = &(plan->f_hat_intern[NFSFT_INDEX(-1,n,plan)]);
00106 xp = &(plan->f_hat_intern[NFSFT_INDEX(+1,n,plan)]);
00107 for(k = 1; k <= plan->N; k++)
00108 {
00109 *xp *= 0.5;
00110 *xm-- = *xp++;
00111 }
00112
00113
00114 *xm = 0.0;
00115 }
00116
00117
00118 low = -plan->N + (1-plan->N%2);
00119 up = -low;
00120
00121
00122
00123 for (n = low; n <= up; n += 2)
00124 {
00125
00126
00127
00128 plan->f_hat_intern[NFSFT_INDEX(0,n,plan)] *= 2.0;
00129 xp = &(plan->f_hat_intern[NFSFT_INDEX(-plan->N-1,n,plan)]);
00130
00131
00132 *xp++ = 0.0;
00133 xm = &(plan->f_hat_intern[NFSFT_INDEX(plan->N,n,plan)]);
00134 last = *xm;
00135 *xm = 0.5 * I * (0.5*xm[-1]);
00136 *xp++ = -(*xm--);
00137 for (k = plan->N-1; k > 0; k--)
00138 {
00139 act = *xm;
00140 *xm = 0.5 * I * (0.5*(xm[-1] - last));
00141 *xp++ = -(*xm--);
00142 last = act;
00143 }
00144 *xm = 0.0;
00145 }
00146 }
00147
00158 inline void c2e_transposed(nfsft_plan *plan)
00159 {
00160 int k;
00161 int n;
00162 double complex last;
00163 double complex act;
00164 double complex *xp;
00165 double complex *xm;
00166 int low;
00167 int up;
00168 int lowe;
00169 int upe;
00171
00172 lowe = -plan->N + (plan->N%2);
00173 upe = -lowe;
00174
00175
00176 for (n = lowe; n <= upe; n += 2)
00177 {
00178
00179
00180 xm = &(plan->f_hat[NFSFT_INDEX(-1,n,plan)]);
00181 xp = &(plan->f_hat[NFSFT_INDEX(+1,n,plan)]);
00182 for(k = 1; k <= plan->N; k++)
00183 {
00184 *xp += *xm--;
00185 *xp++ *= 0.5;;
00186 }
00187 }
00188
00189
00190 low = -plan->N + (1-plan->N%2);
00191 up = -low;
00192
00193
00194 for (n = low; n <= up; n += 2)
00195 {
00196
00197
00198 xm = &(plan->f_hat[NFSFT_INDEX(-1,n,plan)]);
00199 xp = &(plan->f_hat[NFSFT_INDEX(+1,n,plan)]);
00200 for(k = 1; k <= plan->N; k++)
00201 {
00202 *xp++ -= *xm--;
00203 }
00204
00205 plan->f_hat[NFSFT_INDEX(0,n,plan)] =
00206 -0.25*I*plan->f_hat[NFSFT_INDEX(1,n,plan)];
00207 last = plan->f_hat[NFSFT_INDEX(1,n,plan)];
00208 plan->f_hat[NFSFT_INDEX(1,n,plan)] =
00209 -0.25*I*plan->f_hat[NFSFT_INDEX(2,n,plan)];
00210
00211 xp = &(plan->f_hat[NFSFT_INDEX(2,n,plan)]);
00212 for (k = 2; k < plan->N; k++)
00213 {
00214 act = *xp;
00215 *xp = -0.25 * I * (xp[1] - last);
00216 xp++;
00217 last = act;
00218 }
00219 *xp = 0.25 * I * last;
00220
00221 plan->f_hat[NFSFT_INDEX(0,n,plan)] *= 2.0;
00222 }
00223 }
00224
00225 void nfsft_init(nfsft_plan *plan, int N, int M)
00226 {
00227
00228 nfsft_init_advanced(plan, N, M, NFSFT_MALLOC_X | NFSFT_MALLOC_F |
00229 NFSFT_MALLOC_F_HAT);
00230 }
00231
00232 void nfsft_init_advanced(nfsft_plan* plan, int N, int M,
00233 unsigned int flags)
00234 {
00235
00236 nfsft_init_guru(plan, N, M, flags, PRE_PHI_HUT | PRE_PSI | FFTW_INIT |
00237 FFT_OUT_OF_PLACE, NFSFT_DEFAULT_NFFT_CUTOFF);
00238 }
00239
00240 void nfsft_init_guru(nfsft_plan *plan, int N, int M, unsigned int flags,
00241 int nfft_flags, int nfft_cutoff)
00242 {
00243 int *nfft_size;
00244 int *fftw_size;
00245
00246
00247 plan->flags = flags;
00248
00249
00250 plan->N = N;
00251 plan->M_total = M;
00252
00253
00254
00255
00256
00257
00258
00259 plan->N_total = (2*plan->N+2)*(2*plan->N+2);
00260
00261
00262
00263 if (plan->flags & NFSFT_PRESERVE_F_HAT)
00264 {
00265 plan->f_hat_intern = (double complex*) calloc(plan->N_total,
00266 sizeof(double complex));
00267 }
00268
00269
00270 if (plan->flags & NFSFT_MALLOC_F_HAT)
00271 {
00272 plan->f_hat = (double complex*) calloc(plan->N_total,
00273 sizeof(double complex));
00274 }
00275
00276
00277 if (plan->flags & NFSFT_MALLOC_F)
00278 {
00279 plan->f = (double complex*) calloc(plan->M_total,sizeof(double complex));
00280 }
00281
00282
00283 if (plan->flags & NFSFT_MALLOC_X)
00284 {
00285 plan->x = (double*) calloc(plan->M_total,2*sizeof(double));
00286 }
00287
00288
00289 if (plan->flags & NFSFT_NO_FAST_ALGORITHM)
00290 {
00291 }
00292 else
00293 {
00294 nfft_size = (int*)malloc(2*sizeof(int));
00295 fftw_size = (int*)malloc(2*sizeof(int));
00296
00298 nfft_size[0] = 2*plan->N+2;
00299 nfft_size[1] = 2*plan->N+2;
00300 fftw_size[0] = 4*plan->N;
00301 fftw_size[1] = 4*plan->N;
00302
00304 nfft_init_guru(&plan->plan_nfft, 2, nfft_size, plan->M_total, fftw_size,
00305 nfft_cutoff, nfft_flags,
00306 FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
00307
00308
00309 plan->plan_nfft.x = plan->x;
00310
00311 plan->plan_nfft.f = plan->f;
00312
00313 plan->plan_nfft.f_hat = plan->f_hat;
00314
00317
00318
00319
00320
00321 free(nfft_size);
00322 free(fftw_size);
00323 }
00324 }
00325
00326 void nfsft_precompute(int N, double kappa, unsigned int nfsft_flags,
00327 unsigned int fpt_flags)
00328 {
00329 int n;
00330
00331
00332 if (wisdom.initialized == true)
00333 {
00334 return;
00335 }
00336
00337
00338 wisdom.flags = nfsft_flags;
00339
00340
00341
00342 nfft_next_power_of_2_exp(N,&wisdom.N_MAX,&wisdom.T_MAX);
00343
00344
00345 if (wisdom.flags & NFSFT_NO_DIRECT_ALGORITHM)
00346 {
00347 wisdom.alpha = NULL;
00348 wisdom.beta = NULL;
00349 wisdom.gamma = NULL;
00350 }
00351 else
00352 {
00353
00354 wisdom.alpha = (double*) malloc((wisdom.N_MAX+1)*(wisdom.N_MAX+2)*
00355 sizeof(double));
00356 wisdom.beta = (double*) malloc((wisdom.N_MAX+1)*(wisdom.N_MAX+2)*
00357 sizeof(double));
00358 wisdom.gamma = (double*) malloc((wisdom.N_MAX+1)*(wisdom.N_MAX+2)*
00359 sizeof(double));
00361
00362
00363 alpha_al_all(wisdom.alpha,wisdom.N_MAX);
00364 beta_al_all(wisdom.beta,wisdom.N_MAX);
00365 gamma_al_all(wisdom.gamma,wisdom.N_MAX);
00366 }
00367
00368
00369 if (wisdom.flags & NFSFT_NO_FAST_ALGORITHM)
00370 {
00371 }
00372 else if (wisdom.N_MAX >= NFSFT_BREAK_EVEN)
00373 {
00374
00375
00376
00377 if (wisdom.alpha != NULL)
00378 {
00379
00380
00381 wisdom.set = fpt_init(wisdom.N_MAX, wisdom.T_MAX,
00382 fpt_flags | FPT_AL_SYMMETRY | FPT_PERSISTENT_DATA);
00383 for (n = 0; n <= wisdom.N_MAX; n++)
00384 {
00385
00386
00387
00388 fpt_precompute(wisdom.set,n,&wisdom.alpha[ROW(n)],&wisdom.beta[ROW(n)],
00389 &wisdom.gamma[ROW(n)],n,kappa);
00390 }
00391 }
00392 else
00393 {
00394
00395 wisdom.alpha = (double*) malloc((wisdom.N_MAX+2)*sizeof(double));
00396 wisdom.beta = (double*) malloc((wisdom.N_MAX+2)*sizeof(double));
00397 wisdom.gamma = (double*) malloc((wisdom.N_MAX+2)*sizeof(double));
00398 wisdom.set = fpt_init(wisdom.N_MAX, wisdom.T_MAX,
00399 fpt_flags | FPT_AL_SYMMETRY);
00400 for (n = 0; n <= wisdom.N_MAX; n++)
00401 {
00402
00403
00404
00405
00406 alpha_al_row(wisdom.alpha,wisdom.N_MAX,n);
00407 beta_al_row(wisdom.beta,wisdom.N_MAX,n);
00408 gamma_al_row(wisdom.gamma,wisdom.N_MAX,n);
00409
00410
00411 fpt_precompute(wisdom.set,n,wisdom.alpha,wisdom.beta,wisdom.gamma,n,
00412 kappa);
00413 }
00414
00415 free(wisdom.alpha);
00416 free(wisdom.beta);
00417 free(wisdom.gamma);
00418 wisdom.alpha = NULL;
00419 wisdom.beta = NULL;
00420 wisdom.gamma = NULL;
00421 }
00422 }
00423
00424
00425 wisdom.initialized = true;
00426 }
00427
00428 void nfsft_forget()
00429 {
00430
00431 if (wisdom.initialized == false)
00432 {
00433
00434 return;
00435 }
00436
00437
00438 if (wisdom.flags & NFSFT_NO_DIRECT_ALGORITHM)
00439 {
00440 }
00441 else
00442 {
00443
00444 free(wisdom.alpha);
00445 free(wisdom.beta);
00446 free(wisdom.gamma);
00447 wisdom.alpha = NULL;
00448 wisdom.beta = NULL;
00449 wisdom.gamma = NULL;
00450 }
00451
00452
00453 if (wisdom.flags & NFSFT_NO_FAST_ALGORITHM)
00454 {
00455 }
00456 else if (wisdom.N_MAX >= NFSFT_BREAK_EVEN)
00457 {
00458
00459 fpt_finalize(wisdom.set);
00460 }
00461
00462
00463 wisdom.initialized = false;
00464 }
00465
00466
00467 void nfsft_finalize(nfsft_plan *plan)
00468 {
00469
00470 nfft_finalize(&plan->plan_nfft);
00471
00472
00473
00474 if (plan->flags & NFSFT_PRESERVE_F_HAT)
00475 {
00476 free(plan->f_hat_intern);
00477 }
00478
00479
00480 if (plan->flags & NFSFT_MALLOC_F_HAT)
00481 {
00482
00483 free(plan->f_hat);
00484 }
00485
00486
00487 if (plan->flags & NFSFT_MALLOC_F)
00488 {
00489
00490 free(plan->f);
00491 }
00492
00493
00494 if (plan->flags & NFSFT_MALLOC_X)
00495 {
00496
00497 free(plan->x);
00498 }
00499 }
00500
00501 void ndsft_trafo(nfsft_plan *plan)
00502 {
00503 int m;
00504 int k;
00505 int n;
00506 int n_abs;
00507 double *alpha;
00508
00509
00510 double *gamma;
00511
00512
00513 double complex *a;
00514 double complex it1;
00515 double complex it2;
00516 double complex temp;
00517 double complex f_m;
00518
00519 double stheta;
00520 double sphi;
00521
00522 if (wisdom.flags & NFSFT_NO_DIRECT_ALGORITHM)
00523 {
00524 return;
00525 }
00526
00527
00528 if (plan->flags & NFSFT_PRESERVE_F_HAT)
00529 {
00530 memcpy(plan->f_hat_intern,plan->f_hat,plan->N_total*
00531 sizeof(double complex));
00532 }
00533 else
00534 {
00535 plan->f_hat_intern = plan->f_hat;
00536 }
00537
00538
00539
00540
00541 if (plan->flags & NFSFT_NORMALIZED)
00542 {
00543
00544 for (k = 0; k <= plan->N; k++)
00545 {
00546 for (n = -k; n <= k; n++)
00547 {
00548
00549 plan->f_hat_intern[NFSFT_INDEX(k,n,plan)] *=
00550 sqrt((2*k+1)/(4.0*PI));
00551 }
00552 }
00553 }
00554
00555
00556 if (plan->N == 0)
00557 {
00558
00559
00560
00561 for (m = 0; m < plan->M_total; m++)
00562 {
00563 plan->f[m] = plan->f_hat_intern[NFSFT_INDEX(0,0,plan)];
00564 }
00565 }
00566 else
00567 {
00568
00569
00570
00571
00572
00573
00574
00575 for (m = 0; m < plan->M_total; m++)
00576 {
00577
00578 stheta = cos(2.0*PI*plan->x[2*m+1]);
00579
00580 sphi = 2.0*PI*plan->x[2*m];
00581
00582
00583 f_m = 0.0;
00584
00585
00586
00587
00588
00589 for (n = -plan->N; n <= plan->N; n++)
00590 {
00591
00592 a = &(plan->f_hat_intern[NFSFT_INDEX(0,n,plan)]);
00593
00594
00595 n_abs = abs(n);
00596
00597
00598 alpha = &(wisdom.alpha[ROW(n_abs)]);
00599 gamma = &(wisdom.gamma[ROW(n_abs)]);
00600
00601
00602 it2 = a[plan->N];
00603 it1 = a[plan->N-1];
00604 for (k = plan->N; k > n_abs + 1; k--)
00605 {
00606 temp = a[k-2] + it2 * gamma[k];
00607 it2 = it1 + it2 * alpha[k] * stheta;
00608 it1 = temp;
00609 }
00610
00611
00612 if (n_abs < plan->N)
00613 {
00614 it2 = it1 + it2 * wisdom.alpha[ROWK(n_abs)+1] * stheta;
00615 }
00616
00617
00618
00619
00620
00621
00622
00623 f_m += it2 * wisdom.gamma[ROW(n_abs)] *
00624 pow(1- stheta * stheta, 0.5*n_abs) * cexp(I*n*sphi);
00625 }
00626
00627
00628 plan->f[m] = f_m;
00629 }
00630 }
00631 }
00632
00633 void ndsft_adjoint(nfsft_plan *plan)
00634 {
00635 int m;
00636 int k;
00637 int n;
00638 int n_abs;
00639 double *alpha;
00640
00641
00642 double *gamma;
00643
00644
00645 double complex it1;
00646 double complex it2;
00647 double complex temp;
00648 double stheta;
00649 double sphi;
00650
00651 if (wisdom.flags & NFSFT_NO_DIRECT_ALGORITHM)
00652 {
00653 return;
00654 }
00655
00656
00657 memset(plan->f_hat,0U,plan->N_total*sizeof(double complex));
00658
00659
00660 if (plan->N == 0)
00661 {
00662
00663
00664
00665 for (m = 0; m < plan->M_total; m++)
00666 {
00667 plan->f_hat[NFSFT_INDEX(0,0,plan)] += plan->f[m];
00668 }
00669 }
00670 else
00671 {
00672
00673
00674
00675 for (m = 0; m < plan->M_total; m++)
00676 {
00677
00678 stheta = cos(2.0*PI*plan->x[2*m+1]);
00679
00680 sphi = 2.0*PI*plan->x[2*m];
00681
00682
00683 for (n = -plan->N; n <= plan->N; n++)
00684 {
00685
00686 n_abs = abs(n);
00687
00688
00689 alpha = &(wisdom.alpha[ROW(n_abs)]);
00690 gamma = &(wisdom.gamma[ROW(n_abs)]);
00691
00692
00693
00694
00695 it1 = plan->f[m] * wisdom.gamma[ROW(n_abs)] *
00696 pow(1 - stheta * stheta, 0.5*n_abs) * cexp(-I*n*sphi);
00697 plan->f_hat[NFSFT_INDEX(n_abs,n,plan)] += it1;
00698
00699 if (n_abs < plan->N)
00700 {
00701 it2 = it1 * wisdom.alpha[ROWK(n_abs)+1] * stheta;
00702 plan->f_hat[NFSFT_INDEX(n_abs+1,n,plan)] += it2;
00703 }
00704
00705
00706 for (k = n_abs+2; k <= plan->N; k++)
00707 {
00708 temp = it2;
00709 it2 = alpha[k] * stheta * it2 + gamma[k] * it1;
00710 it1 = temp;
00711 plan->f_hat[NFSFT_INDEX(k,n,plan)] += it2;
00712 }
00713 }
00714 }
00715 }
00716
00717
00718
00719
00720 if (plan->flags & NFSFT_NORMALIZED)
00721 {
00722
00723 for (k = 0; k <= plan->N; k++)
00724 {
00725 for (n = -k; n <= k; n++)
00726 {
00727
00728 plan->f_hat[NFSFT_INDEX(k,n,plan)] *=
00729 sqrt((2*k+1)/(4.0*PI));
00730 }
00731 }
00732 }
00733
00734
00735 if (plan->flags & NFSFT_ZERO_F_HAT)
00736 {
00737 for (n = -plan->N; n <= plan->N+1; n++)
00738 {
00739 memset(&plan->f_hat[NFSFT_INDEX(-plan->N-1,n,plan)],0U,
00740 (plan->N+1+abs(n))*sizeof(double complex));
00741 }
00742 }
00743 }
00744
00745 void nfsft_trafo(nfsft_plan *plan)
00746 {
00747 int k;
00748 int n;
00749 #ifdef DEBUG
00750 double t, t_pre, t_nfft, t_fpt, t_c2e, t_norm;
00751 t_pre = 0.0;
00752 t_norm = 0.0;
00753 t_fpt = 0.0;
00754 t_c2e = 0.0;
00755 t_nfft = 0.0;
00756 #endif
00757
00758 if (wisdom.flags & NFSFT_NO_FAST_ALGORITHM)
00759 {
00760 return;
00761 }
00762
00763
00764
00765
00766 if (wisdom.initialized == 0 || plan->N > wisdom.N_MAX)
00767 {
00768 return;
00769 }
00770
00771
00772 if (plan->N < NFSFT_BREAK_EVEN)
00773 {
00774
00775 ndsft_trafo(plan);
00776 }
00777
00778
00779 else if (plan->N <= wisdom.N_MAX)
00780 {
00781
00782 if (plan->flags & NFSFT_PRESERVE_F_HAT)
00783 {
00784 memcpy(plan->f_hat_intern,plan->f_hat,plan->N_total*
00785 sizeof(double complex));
00786 }
00787 else
00788 {
00789 plan->f_hat_intern = plan->f_hat;
00790 }
00791
00792
00793
00794
00795 plan->plan_nfft.x = plan->x;
00796 plan->plan_nfft.f = plan->f;
00797 plan->plan_nfft.f_hat = plan->f_hat_intern;
00798
00799
00800
00801
00802 if (plan->flags & NFSFT_NORMALIZED)
00803 {
00804
00805 for (k = 0; k <= plan->N; k++)
00806 {
00807 for (n = -k; n <= k; n++)
00808 {
00809
00810 plan->f_hat_intern[NFSFT_INDEX(k,n,plan)] *=
00811 sqrt((2*k+1)/(4.0*PI));
00812 }
00813 }
00814 }
00815
00816
00817 if (plan->flags & NFSFT_USE_DPT)
00818 {
00819
00820 for (n = -plan->N; n <= plan->N; n++)
00821 {
00822
00823 fflush(stderr);
00824 dpt_trafo(wisdom.set,abs(n),
00825 &plan->f_hat_intern[NFSFT_INDEX(abs(n),n,plan)],
00826 &plan->f_hat_intern[NFSFT_INDEX(0,n,plan)],
00827 plan->N,0U);
00828 }
00829 }
00830 else
00831 {
00832
00833 for (n = -plan->N; n <= plan->N; n++)
00834 {
00835
00836 fflush(stderr);
00837 fpt_trafo(wisdom.set,abs(n),
00838 &plan->f_hat_intern[NFSFT_INDEX(abs(n),n,plan)],
00839 &plan->f_hat_intern[NFSFT_INDEX(0,n,plan)],
00840 plan->N,0U);
00841 }
00842 }
00843
00844
00845 c2e(plan);
00846
00847
00848
00849
00850 if (plan->flags & NFSFT_USE_NDFT)
00851 {
00852
00853 ndft_trafo(&plan->plan_nfft);
00854 }
00855 else
00856 {
00857
00858
00859 nfft_trafo(&plan->plan_nfft);
00860 }
00861 }
00862 }
00863
00864 void nfsft_adjoint(nfsft_plan *plan)
00865 {
00866 int k;
00867 int n;
00868
00869 if (wisdom.flags & NFSFT_NO_FAST_ALGORITHM)
00870 {
00871 return;
00872 }
00873
00874
00875
00876
00877 if (wisdom.initialized == 0 || plan->N > wisdom.N_MAX)
00878 {
00879 return;
00880 }
00881
00882
00883 if (plan->N < NFSFT_BREAK_EVEN)
00884 {
00885
00886 ndsft_adjoint(plan);
00887 }
00888
00889 else if (plan->N <= wisdom.N_MAX)
00890 {
00891
00892
00893
00894
00895
00896 plan->plan_nfft.x = plan->x;
00897 plan->plan_nfft.f = plan->f;
00898 plan->plan_nfft.f_hat = plan->f_hat;
00899
00900
00901
00902
00903 if (plan->flags & NFSFT_USE_NDFT)
00904 {
00905
00906
00907
00908 ndft_adjoint(&plan->plan_nfft);
00909 }
00910 else
00911 {
00912
00913
00914
00915
00916 nfft_adjoint(&plan->plan_nfft);
00917 }
00918
00919
00920
00921
00922 c2e_transposed(plan);
00923
00924
00925 if (plan->flags & NFSFT_USE_DPT)
00926 {
00927
00928 for (n = -plan->N; n <= plan->N; n++)
00929 {
00930
00931
00932 dpt_transposed(wisdom.set,abs(n),
00933 &plan->f_hat[NFSFT_INDEX(abs(n),n,plan)],
00934 &plan->f_hat[NFSFT_INDEX(0,n,plan)],
00935 plan->N,0U);
00936 }
00937 }
00938 else
00939 {
00940
00941
00942 for (n = -plan->N; n <= plan->N; n++)
00943 {
00944
00945
00946 fpt_transposed(wisdom.set,abs(n),
00947 &plan->f_hat[NFSFT_INDEX(abs(n),n,plan)],
00948 &plan->f_hat[NFSFT_INDEX(0,n,plan)],
00949 plan->N,0U);
00950 }
00951 }
00952
00953
00954
00955
00956 if (plan->flags & NFSFT_NORMALIZED)
00957 {
00958
00959
00960
00961 for (k = 0; k <= plan->N; k++)
00962 {
00963 for (n = -k; n <= k; n++)
00964 {
00965
00966 plan->f_hat[NFSFT_INDEX(k,n,plan)] *=
00967 sqrt((2*k+1)/(4.0*PI));
00968 }
00969 }
00970 }
00971
00972
00973 if (plan->flags & NFSFT_ZERO_F_HAT)
00974 {
00975
00976
00977 for (n = -plan->N; n <= plan->N+1; n++)
00978 {
00979 memset(&plan->f_hat[NFSFT_INDEX(-plan->N-1,n,plan)],0U,
00980 (plan->N+1+abs(n))*sizeof(double complex));
00981 }
00982 }
00983
00984
00985 }
00986 }
00987
00988 void nfsft_precompute_x(nfsft_plan *plan)
00989 {
00990
00991 plan->plan_nfft.x = plan->x;
00992
00993
00994 if (plan->plan_nfft.nfft_flags & PRE_ONE_PSI)
00995 nfft_precompute_one_psi(&plan->plan_nfft);
00996 }
00997