sample()

template<class function >

void fourier::adft_func::sample ( double  w, double  a, double  b, function &  fz, int  n, int  limit  )

inline

Fourier samples the given function to determine the individual intervals.

Purpose

Computes the Fourier integral for fixed frequency \(\omega\). The interval \([a,b]\) is iteratively split until the error of the Fourier integral over each subinterval is smaller or equal than the smallest error.

Parameters

w	frequency \(\omega\) in the Fourier integral
a	lower interval bound
b	upper interval bound
fz	[template class function] scalar complex function
n	number of points in each interval
limit	max number of intervals

Definition at line 366 of file fourier.hpp.

References fourier::dft_cplx().

Referenced by cntr::green_equilibrium_mat(), and cntr::smooth_box::set_norm().

                                                                       {
       // n:      number of points in each interval
       // limit:  max number of intervals 
       int i,n_cur,isplit;
       double abserr,abserr_i;
       cplx z,res_i,err_i,res_l,err_l,res_r,err_r;
       std::vector<cplx> results(limit),errors(limit);
    
       if(n<ADFT_MINPTS || n%2!=0){ std::cerr << "adft_cplx_sample: n<16 || n%2!=0" << std::endl; abort();}
       if(limit<1){std::cerr << "adft_cplx_sample: limit<1" << std::endl; abort();}
       if(a>=b){ std::cerr << "adft_cplx_sample: a>=b" << std::endl;abort();}
       
       if(limit != nmax_ || (n+1)*limit != size_ ) resize((n+1)*limit,limit);  
       // Sample f in first interval
       init(a,b,n,fz);
       if(limit==1) return;
       // split first interval
       split(0,fz);
       // do the first fourier transform
       fourier::dft_cplx(w,n,a_[0],b_[0],interval_data(0),results[0],errors[0]);
       fourier::dft_cplx(w,n,a_[1],b_[1],interval_data(1),results[1],errors[1]);
       while(n_+1<limit){
          // look for largest error in res-table
          n_cur=n_;
          abserr=-1.0;
          isplit=-1;
          for(i=0;i<n_cur;i++){
             z=errors[i];
             abserr_i=sqrt(z.real()*z.real()+z.imag()*z.imag());
             if(abserr_i>abserr){
               isplit=i;
               abserr=abserr_i;
             }
          }
          // split interval isplit and integrate partial intervals

          //std::cout << "isplit= "<< isplit << "n= "<< n_ << std::endl; 
//          std::cout << "n " << n_ << "i " << isplit << "[ " << a_[isplit] << "," << b_[isplit] << "]" << std::endl; 
          split(isplit,fz);

          fourier::dft_cplx(w,n,a_[isplit],b_[isplit],interval_data(isplit),res_l,err_l);
          fourier::dft_cplx(w,n,a_[isplit+1],b_[isplit+1],interval_data(isplit+1),res_r,err_r);
          // update res-table
          for(i=n_cur-1;i>isplit;i--){
             results[i+1]=results[i];
             errors[i+1]=errors[i];
          }
          results[isplit+1]=res_r;
          errors[isplit+1]=err_r;
          results[isplit]=res_l;
          errors[isplit]=err_l;
       }
   }        

Here is the call graph for this function:

Here is the caller graph for this function: