Page MenuHomec4science
Files F69284232

surface_generator_filter_bessel.cpp
No OneTemporary
Actions

File Metadata

Created
Mon, Jul 1, 02:41

surface_generator_filter_bessel.cpp
View Options

/* -------------------------------------------------------------------------- */
#include "surface_generator_filter_bessel.hh"
#include <math.h>
#include <fstream>
/* -------------------------------------------------------------------------- */
SurfaceGeneratorFilterBessel::SurfaceGeneratorFilterBessel(const std::string & inputfile):
SurfaceGenerator(inputfile),SurfaceGeneratorFilter(inputfile){
}
/* -------------------------------------------------------------------------- */
Real SurfaceGeneratorFilterBessel::besselFunction1(Real x){
Real ax,z;
Real xx,y,ans,ans1,ans2;
if ((ax=fabs(x)) < 8.0) { //Direct rational approximation.
y=x*x;
ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1
+y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
ans2=144725228442.0+y*(2300535178.0+y*(18583304.74
+y*(99447.43394+y*(376.9991397+y*1.0))));
ans=ans1/ans2;
} else { //Fitting function (6.5.9).
z=8.0/ax;
y=z*z;
xx=ax-2.356194491;
ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
+y*(0.2457520174e-5+y*(-0.240337019e-6))));
ans2=0.04687499995+y*(-0.2002690873e-3
+y*(0.8449199096e-5+y*(-0.88228987e-6
+y*0.105787412e-6)));
ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
if (x < 0.0) ans = -ans;
}
return ans;
}
// Real SurfaceGeneratorFilterBessel::besselFunction1(Real x,int order){
// Real res = 0;
// Real sign = 1;
// Real x_2p = 1;
// Real p_factorial = 1;
// Real two_power2p = 1;
// for (int p = 0 ; p < order ; ++p){
// res += sign*x_2p/two_power2p/p_factorial/p_factorial/(p+1);
// /* fprintf(stderr,"(p = %d) adding %e = %e * %e / %e / %e ^2 / %d => %e (x=%e)\n",
// p,
// sign*x_2p/two_power2p/p_factorial/p_factorial/(p+1),
// sign,
// x_2p,
// two_power2p,
// p_factorial,
// p+1, res,x);
// */
// sign *= -1;
// x_2p *= x*x;
// two_power2p *= 4;
// p_factorial *= p+1;
// }
// if (isnan(res) || isinf(res)) {
// fprintf(stderr,"bessel function computation failed %.15e %d %.15e\n",x,order,res);
// exit(-1);
// }
// return res*(x/2);
// }
inline Real sinc(Real x){
return sin(M_PI*x)/M_PI/x;
}
/* -------------------------------------------------------------------------- */
Real SurfaceGeneratorFilterBessel::lanczosFunction(Real x,int N){
if (fabs(x) < N) return sinc(2*x/N);
return 0;
}
/* -------------------------------------------------------------------------- */
void SurfaceGeneratorFilterBessel::computeFilterCoefficients(){
int n = surface->size();
h_coeff = new Surface<Real>(n,1.);
for (int i = 0 ; i < n ; ++i)
for (int j = 0 ; j < n ; ++j){
Real A;
if (i == 0 && j == 0){
A = 1;
}
else
A = sqrt(i*i+j*j);
Real x = 2.*M_PI*alpha*A;
Real J1 = besselFunction1(x);
(*h_coeff)(i,j) = alpha*1.6*J1/A*lanczosFunction(A,n);//*alpha
// h[i+j*n] = lanczosFunction(A,n/2);
if (isnan((*h_coeff)(i,j)) || isinf((*h_coeff)(i,j)))
{
fprintf(stderr,"h computation failed %.15e %.15e %.15e\n",
alpha,J1,sqrt(i*i+j*j));
exit(-1);
}
}
}
/* -------------------------------------------------------------------------- */
void SurfaceGeneratorFilterBessel::parseInputFile(std::ifstream & fp) {
std::string line;
int gridsize;
while (! fp.eof() )
{
getline (fp,line);
size_t pos = line.find("#"); //remove the comment
line = line.substr(0, pos);
trim(line); // remove unnecessary spaces
/* Read all options one by one */
readInput(line, "seed" , random_seed);
readInput(line, "gridsize", gridsize );
readInput(line, "alpha" , alpha );
}
return;
}

Event Timeline

c4science · Help