20 std::complex<double> *endpts,std::complex<double> *endcor,
double *corfac);
21 void dft_cplx(
double w,
int n,
double a,
double b,std::complex<double> *f,
22 std::complex<double> &res, std::complex<double> &err);
25 #define PI 3.14159265358979323846 26 #define ADFT_MINPTS 16 44 typedef std::complex<double>
cplx;
75 assert(size>0 && nmax>0);
76 data_=
new cplx [size];
96 std::cerr <<
"= not implemented for adft" << std::endl;
101 std::cerr <<
" copy constructor not implemented for for adft" << std::endl;
122 void resize(
int size,
int nmax){
123 if(size<0 || nmax<0) {std::cerr <<
"adft: size<0 or nmax<0" << std::endl;abort();}
129 data_=
new cplx [size];
134 a_=
new double [nmax];
135 b_=
new double [nmax];
155 cplx *interval_data(
int i){
156 if(i>=n_) {std::cerr <<
"interval i="<<i<<
" >n="<< n_<<
" in adft" << std::endl;abort();}
181 template <
class function>
182 void init(
double a,
double b,
int n,
function &fz){
187 if(a>=b){std::cerr <<
"adft_func: a>=b" << std::endl;abort();}
188 if(n<1){std::cerr <<
"adft_func: n<1" << std::endl;abort();}
189 if(nmax_<1 || size_<n) resize(n+1,1);
196 dx=(b-a)/((
double) n);
221 template <
class function>
222 void split(
int interv,
function &fz){
224 int j,j0,j1,jm,len0,n0=n_,i1=interv+1,tot0;
225 double a0,b0,mid,x,dx;
229 if(interv>=n0 || interv<0) {std::cerr <<
"adft_cplx_split: interv out of range" << std::endl; abort();}
237 tot_ += len_[interv] + 1 ;
238 if(n_ > nmax_ || tot_ > size_ ) {std::cerr <<
"adft_cplx_split: something too small" << std::endl; abort();}
241 for(j=n0-1;j>interv;j--){
255 f0=interval_data(interv);
258 dx=(b0-a0)/(2.0*len0);
263 for(j0=len0-1;j0>=jm;j0--){
271 for(j0=jm;j0>0;j0--){
281 for(j0=len0;j0>=jm;j0--){
289 for(j0=jm-1;j0>=0;j0--){
320 cplx res_i,err_i,*f_i;
321 if(n_<=0) {std::cerr <<
"adft_cplx: no interval in f" << std::endl; abort();}
326 f_i=interval_data(i);
365 template <
class function>
366 void sample(
double w,
double a,
double b,
function &fz,
int n,
int limit){
370 double abserr,abserr_i;
371 cplx z,res_i,err_i,res_l,err_l,res_r,err_r;
372 std::vector<cplx> results(limit),errors(limit);
374 if(n<ADFT_MINPTS || n%2!=0){ std::cerr <<
"adft_cplx_sample: n<16 || n%2!=0" << std::endl; abort();}
375 if(limit<1){std::cerr <<
"adft_cplx_sample: limit<1" << std::endl; abort();}
376 if(a>=b){ std::cerr <<
"adft_cplx_sample: a>=b" << std::endl;abort();}
378 if(limit != nmax_ || (n+1)*limit != size_ ) resize((n+1)*limit,limit);
392 for(i=0;i<n_cur;i++){
394 abserr_i=sqrt(z.real()*z.real()+z.imag()*z.imag());
407 fourier::dft_cplx(w,n,a_[isplit+1],b_[isplit+1],interval_data(isplit+1),res_r,err_r);
409 for(i=n_cur-1;i>isplit;i--){
410 results[i+1]=results[i];
411 errors[i+1]=errors[i];
413 results[isplit+1]=res_r;
414 errors[isplit+1]=err_r;
415 results[isplit]=res_l;
416 errors[isplit]=err_l;
std::complex< double > cplx
Class adft_func contains functions for computing accurate Fourier transforms by adaptive splitting o...
std::complex< double > cplx
void get_dftcorr_linear(double th, double *corfac, cplx *endcor)
Returns correction factor and boundary correction term needed for linearly corrected DFT...
void sample(double w, double a, double b, function &fz, int n, int limit)
Fourier samples the given function to determine the individual intervals.
void dft_cplx(double w, int n, double a, double b, cplx *f, cplx &res, cplx &err)
Computes the Fourier integral using cubically corrected DFT.
void complex_dftcor_cubic(double w, double delta, double a, double b, cplx *endpts, cplx *endcor, double *corfac)
Returns correction factor and evaluate boundary correction terms needed for cubically corrected DFT...
void dft(double w, cplx &result, cplx &err)
Computes the cubically corrected DFT.
void get_dftcorr_cubic(double th, double *corfac, cplx *endcor)
Returns correction factor and boundary correction terms needed for cubically corrected DFT...