Commit 838609de03e740c2d2ee1c191d1fc42aaae3a938

  • avatar
  • Heikki Mäntysaari <heikki.mantysaari @j…u.fi> (Committer)
  • Fri Aug 19 10:44:36 EEST 2011
  • avatar
  • Heikki Mäntysaari <heikki.mantysaari @j…u.fi> (Author)
  • Fri Aug 19 10:44:36 EEST 2011
fourier.c BesselJ transformation handles situation where
integral over the invervals starts to give zero as a result
after a few nonzero terms.
src/fourier/fourier.c
(42 / 13)
  
4646
4747static int Nzeros;
4848//static gsl_sum_levin_u_workspace *workspace;
49static double epsilon1=1.0e-12;
50static double epsilon=1.0e-12;
49static double epsilon1=1.0e-7; //1.0e-12;
50static double epsilon=1.0e-7; //1.0e-12;
5151
5252
5353/* Table with the first 1000 zeros of J0(x) */
5454
55//double J0zero[2001]={0.0,\ ...
55/* double J0zero[2001]={0.0,\ ... */
5656#include "besselj0zero.h"
5757
5858
5959
60//double J1zero[1001]={0.0,\ ...
60/* double J1zero[1001]={0.0,\ ... */
6161#include "besselj1zero.h"
6262
6363
64// Generic integration routine in an interval
64/* Generic integration routine in an interval */
6565
6666double integrate(double a,double b,gsl_function F,\
6767 gsl_integration_workspace *gsl_wksp){
6868 double error,term_n;
6969
70 gsl_integration_qag(&F,a,b,0.0,epsilon1,4000,6,gsl_wksp,&term_n,&error);
70 // Error handling added by H.M.
71 int status
72 = gsl_integration_qag(&F,a,b,0.0,epsilon1,4000,6,gsl_wksp,&term_n,&error);
7173
74 if (status && fabs(term_n)>epsilon1)
75 {
76 fprintf(stderr, "Integration failed at fourier.c:integrate interval "
77 " [%f, %f], result %E, relerr %E\n", a,b,term_n, fabs(error/term_n));
78 }
79
7280 return term_n;
7381}
7482
8787 fourier_int_data *dt=(fourier_int_data*)param;
8888
8989#ifndef FORTRAN
90 return j0(k*(dt->x))*(dt->function)(k,dt->params);
90 /* changed by H.M. */
91 double result = j0(k*(dt->x))*(dt->function)(k,dt->params);
92 if (isinf(result) || isnan(result))
93 fprintf(stderr, "Inf/nan result %f at k=%f\n", result, k);
94 return result;
9195#else
9296 return j0(k*(dt->x))*(dt->function)(&k,dt->params);
9397#endif
126126 int N;
127127 int dN=20;
128128 int precision_OK=0;
129 int j;
129 int i,j;
130130 fourier_int_data params;
131131 gsl_sum_levin_u_workspace *workspace=gsl_sum_levin_u_alloc (1+Nzeros);
132132 gsl_integration_workspace *gsl_wksp=gsl_integration_workspace_alloc(4000);
156156 do {
157157 for(;j<N+Ni;j++){
158158 tmp[j-Ni]=integrate(zeros[j]/x,zeros[j+1]/x,F,gsl_wksp);
159
160 /* All terms given to gls_sum_levin_u_accel must be nonzero.
161 * If we just got 0, most likely all the other terms in the series
162 * are zero (within the numerical accuracy), thus compute the sum of
163 * the nonzero terms and return it.
164 * Added by H.M. */
165 if (tmp[j-Ni]==0)
166 {
167 sum=0;
168 for (i=0; i<=j-Ni; i++)
169 sum += tmp[i];
170 gsl_integration_workspace_free(gsl_wksp);
171 free(tmp);
172 gsl_sum_levin_u_free(workspace);
173 return sum;
174 }
159175 }
160176 gsl_sum_levin_u_accel(tmp,N,workspace,&sum,&error);
161177 if (workspace->terms_used<N) {
162178 precision_OK=1;
163 } else {
179 }
180 else {
164181 if (fabs(error/sum)<=epsilon) {
165 precision_OK=1;
166 } else {
167 N+=dN;
182 precision_OK=1;
168183 }
184 else {
185 N+=dN;
186 }
169187 }
170188 } while((precision_OK==0)&&(N<=Nzeros));
171189
296296
297297 // We have to sum a finite number of terms ==> we must not use sum
298298 // acceleration algorithms !!
299
300299 for(j=Ni;j<Nf;j++){
301300 sum+=integrate(zeros[j]/x,zeros[j+1]/x,F,gsl_wksp);
302301 }