Homec4science
Diffusion tamaas (ci-fix-static-checks) Browse

tamaas/src/surface_generator_filter_bessel.cpp
bd979ec5f10fci-fix-static-checks

Last Change
Raw File

surface_generator_filter_bessel.cpp
View Options

/**
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de
* Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
* Solides)
*
* Tamaas is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
#include "surface_generator_filter_bessel.hh"
#include <cmath>
#include <fstream>
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
using namespace std;
SurfaceGeneratorFilterBessel::SurfaceGeneratorFilterBessel():
SurfaceGeneratorFilter(){
}
/* -------------------------------------------------------------------------- */
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);
}
}
}
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
c4science · Help