diff --git a/App/src/TestMesoCrystal.cpp b/App/src/TestMesoCrystal.cpp
index e4d98743b4ec2cb9c4274d2e594235589fad1052..4c81803a276b314da17663b5ee5e8f37a196ca9e 100644
--- a/App/src/TestMesoCrystal.cpp
+++ b/App/src/TestMesoCrystal.cpp
@@ -68,14 +68,14 @@ void TestMesoCrystal::initializeSample()
{
delete mp_sample;
// create mesocrystal
- double meso_radius = 300*Units::nanometer;
+ double meso_radius = 600*Units::nanometer;
double surface_filling_ratio = 0.25;
- double surface_density = surface_filling_ratio/M_PI/meso_radius/meso_radius;
+ double surface_density = surface_filling_ratio/meso_radius/meso_radius;
complex_t n_particle(1.0-1.55e-5, 1.37e-6);
complex_t avg_n_squared_meso = 0.7886*n_particle*n_particle + 0.2114;
complex_t n_avg = std::sqrt(surface_filling_ratio*avg_n_squared_meso + 1.0 - surface_filling_ratio);
complex_t n_particle_adapted = std::sqrt(n_avg*n_avg + n_particle*n_particle - 1.0);
- FormFactorCylinder ff_meso(0.2*Units::micrometer, meso_radius);
+ FormFactorLorentz ff_meso(0.2*Units::micrometer, meso_radius);
// MesoCrystal meso2(npc.clone(), new FormFactorPyramid(0.2*Units::micrometer, meso_radius, 84*Units::degree));
// Create multilayer
@@ -97,9 +97,9 @@ void TestMesoCrystal::initializeSample()
// IInterferenceFunction *p_interference_funtion = new InterferenceFunction1DParaCrystal(800.0*Units::nanometer,
// 50*Units::nanometer, 1e7*Units::nanometer);
NanoParticleDecoration particle_decoration;
- size_t n_phi_rotation_steps = 1;
- size_t n_alpha_rotation_steps = 1;
- size_t n_np_size_steps = 1;
+ size_t n_phi_rotation_steps = 13;
+ size_t n_alpha_rotation_steps = 5;
+ size_t n_np_size_steps = 3;
double phi_step = 2*M_PI/3.0/n_phi_rotation_steps;
double phi_start = 0.0;
double alpha_step = 4*Units::degree/n_alpha_rotation_steps;
@@ -149,7 +149,7 @@ MesoCrystal* createMesoCrystal(double nanoparticle_radius, complex_t n_particle,
pos_vector.push_back(position_2);
LatticeBasis basis(particle, pos_vector);
NanoParticleCrystal npc(basis, lat);
- double relative_sigma_np_radius = 0.0;
+ double relative_sigma_np_radius = 0.15;
double dw_factor = relative_sigma_np_radius*relative_sigma_np_radius*nanoparticle_radius*nanoparticle_radius/6.0;
npc.setDWFactor(dw_factor);
return new MesoCrystal(npc.clone(), p_meso_form_factor->clone());
diff --git a/Core/Samples/inc/FormFactorCylinder.h b/Core/Samples/inc/FormFactorCylinder.h
index fb98d84b9f37e6e94fd5f278d69ccfcd98b63f6b..d60d2b40b9bdf5d47be72b12db3899a2c9c39ecd 100644
--- a/Core/Samples/inc/FormFactorCylinder.h
+++ b/Core/Samples/inc/FormFactorCylinder.h
@@ -30,7 +30,6 @@ public:
protected:
virtual complex_t evaluate_for_q(cvector_t q) const;
- virtual complex_t evaluate_for_complex_qz(kvector_t q, complex_t qz) const;
private:
//! copy constructor and assignment operator are hidden since there is a clone method
diff --git a/Core/Samples/inc/FormFactorGauss.h b/Core/Samples/inc/FormFactorGauss.h
index 92a6efb075a09e7ebb22f6ce5aa40c6a2aeb8231..c389b5ea6e7b5897e9de5e48d546d7b279cf7a03 100644
--- a/Core/Samples/inc/FormFactorGauss.h
+++ b/Core/Samples/inc/FormFactorGauss.h
@@ -30,7 +30,6 @@ public:
protected:
virtual complex_t evaluate_for_q(cvector_t q) const;
- virtual complex_t evaluate_for_complex_qz(kvector_t q, complex_t qz) const;
private:
//! copy constructor and assignment operator are hidden since there is a clone method
diff --git a/Core/Samples/inc/FormFactorLorentz.h b/Core/Samples/inc/FormFactorLorentz.h
index 945f7f89e448121bfa1245a6df458dc7b3a8fd24..76b42d8be01ccdd9b8f80b458fdc15c645ce8a62 100644
--- a/Core/Samples/inc/FormFactorLorentz.h
+++ b/Core/Samples/inc/FormFactorLorentz.h
@@ -30,7 +30,6 @@ public:
protected:
virtual complex_t evaluate_for_q(cvector_t q) const;
- virtual complex_t evaluate_for_complex_qz(kvector_t q, complex_t qz) const;
private:
//! copy constructor and assignment operator are hidden since there is a clone method
diff --git a/Core/Samples/src/FormFactorCylinder.cpp b/Core/Samples/src/FormFactorCylinder.cpp
index eb17e8cf97295677103e4f7edfa659b6f4f8c08b..86824b4714b27d854d41a5a50c3d5dabb60ab808 100644
--- a/Core/Samples/src/FormFactorCylinder.cpp
+++ b/Core/Samples/src/FormFactorCylinder.cpp
@@ -44,23 +44,6 @@ FormFactorCylinder* FormFactorCylinder::clone() const
// return new FormFactorCylinder(mp_height->clone(), mp_radius->clone());
}
-complex_t FormFactorCylinder::evaluate_for_complex_qz(kvector_t q, complex_t qz) const
-{
-// double R = mp_radius->getCurrent();
-// double H = mp_height->getCurrent();
- double R = m_radius;
- double H = m_height;
-
- complex_t qzH_half = qz*H/2.0;
- complex_t z_part = H*MathFunctions::Sinc(qzH_half)*std::exp(complex_t(0.0, 1.0)*qzH_half);
-
- double qrR = q.magxy()*R;
- double J1_qrR_div_qrR = qrR > Numeric::double_epsilon ? MathFunctions::Bessel_J1(qrR)/qrR : 0.5;
- double radial_part = 2*M_PI*R*R*J1_qrR_div_qrR;
-
- return radial_part*z_part;
-}
-
complex_t FormFactorCylinder::evaluate_for_q(cvector_t q) const
{
// double R = mp_radius->getCurrent();
@@ -72,7 +55,7 @@ complex_t FormFactorCylinder::evaluate_for_q(cvector_t q) const
complex_t z_part = H*MathFunctions::Sinc(qzH_half)*std::exp(complex_t(0.0, 1.0)*qzH_half);
complex_t qrR = q.magxy()*R;
- complex_t J1_qrR_div_qrR = std::abs(qrR) > Numeric::double_epsilon ? MathFunctions::Bessel_J1(qrR)/qrR : 0.5;
+ complex_t J1_qrR_div_qrR = std::abs(qrR) > Numeric::double_epsilon ? MathFunctions::Bessel_J1(std::abs(qrR))/qrR : 0.5;
complex_t radial_part = 2*M_PI*R*R*J1_qrR_div_qrR;
return radial_part*z_part;
diff --git a/Core/Samples/src/FormFactorGauss.cpp b/Core/Samples/src/FormFactorGauss.cpp
index 3cef83a4c0d47ef57300d84f7aaa5132d8a69ea6..306b74a1033df7356d9e1170b00121a500730fde 100644
--- a/Core/Samples/src/FormFactorGauss.cpp
+++ b/Core/Samples/src/FormFactorGauss.cpp
@@ -54,32 +54,16 @@ FormFactorGauss* FormFactorGauss::clone() const
// return new FormFactorGauss(mp_height->clone(), mp_radius->clone());
}
-complex_t FormFactorGauss::evaluate_for_complex_qz(kvector_t q, complex_t qz) const
-{
-// double R = mp_radius->getCurrent();
-// double H = mp_height->getCurrent();
- double R = m_width;
- double H = m_height;
-
- complex_t z_part = H*std::exp(-qz*qz*H*H/4.0/M_PI);
-
- double qrR = q.magxy()*R;
- double radial_part = R*R*std::exp(-qrR*qrR/4.0/M_PI);
-
- return radial_part*z_part;
-}
-
complex_t FormFactorGauss::evaluate_for_q(cvector_t q) const
{
-// double R = mp_radius->getCurrent();
+// double R = mp_width->getCurrent();
// double H = mp_height->getCurrent();
double R = m_width;
double H = m_height;
- complex_t z_part = H*std::exp(-q.z()*q.z()*H*H/4.0/M_PI);
+ complex_t z_part = H*std::exp(-q.z()*q.z()*H*H/8.0/M_PI);
- complex_t qrR = q.magxy()*R;
- complex_t radial_part = R*R*std::exp(-qrR*qrR/4.0/M_PI);
+ complex_t radial_part = R*R*std::exp(-(q.x()*q.x()+q.y()*q.y())*R*R/8.0/M_PI);
return radial_part*z_part;
}
diff --git a/Core/Samples/src/FormFactorLorentz.cpp b/Core/Samples/src/FormFactorLorentz.cpp
index 506690a8dc2b99f0dfd6dbcf405727acf1d596ec..702c09099eb04154b5c6f0df54aa7edd06532661 100644
--- a/Core/Samples/src/FormFactorLorentz.cpp
+++ b/Core/Samples/src/FormFactorLorentz.cpp
@@ -54,20 +54,6 @@ FormFactorLorentz* FormFactorLorentz::clone() const
// return new FormFactorLorentz(mp_height->clone(), mp_radius->clone());
}
-complex_t FormFactorLorentz::evaluate_for_complex_qz(kvector_t q, complex_t qz) const
-{
-// double R = mp_radius->getCurrent();
-// double H = mp_height->getCurrent();
- double R = m_width;
- double H = m_height;
-
- complex_t z_part = H/(1.0 + (H*qz/2.0)*(H*qz/2.0));
-
- double radial_part = R*R/(1.0 + (R*q.x()/2.0)*(R*q.x()/2.0))/(1.0 + (R*q.y()/2.0)*(R*q.y()/2.0));
-
- return radial_part*z_part;
-}
-
complex_t FormFactorLorentz::evaluate_for_q(cvector_t q) const
{
// double R = mp_radius->getCurrent();
@@ -75,9 +61,9 @@ complex_t FormFactorLorentz::evaluate_for_q(cvector_t q) const
double R = m_width;
double H = m_height;
- complex_t z_part = H/(1.0 + (H*q.z()/2.0)*(H*q.z()/2.0));
+ complex_t z_part = H/std::sqrt(1.0 + (H*q.z()/2.0)*(H*q.z()/2.0));
- complex_t radial_part = R*R/(1.0 + (R*q.x()/2.0)*(R*q.x()/2.0))/(1.0 + (R*q.y()/2.0)*(R*q.y()/2.0));
+ complex_t radial_part = R*R/std::sqrt(1.0 + (R*q.x()/2.0)*(R*q.x()/2.0))/std::sqrt(1.0 + (R*q.y()/2.0)*(R*q.y()/2.0));
return radial_part*z_part;
}
diff --git a/Core/Tools/inc/MathFunctions.h b/Core/Tools/inc/MathFunctions.h
index 7670c9d574f3d349cec9285305a2b8163280fcff..777315e9574adc399fb5c040874720f5b34a63c9 100644
--- a/Core/Tools/inc/MathFunctions.h
+++ b/Core/Tools/inc/MathFunctions.h
@@ -38,7 +38,7 @@ double GenerateUniformRandom();
double Bessel_J1(double value);
-complex_t Bessel_J1(complex_t value);
+//complex_t Bessel_J1(complex_t value);
double Sinc(double value);
@@ -81,8 +81,7 @@ inline complex_t MathFunctions::Sinc(complex_t value) // Sin(x)/x
if(std::abs(value)<Numeric::double_epsilon) {
return complex_t(1.0, 0.0);
}
- return (std::exp(complex_t(0.0, 1.0)*value) - std::exp(complex_t(0.0, -1.0)*value))
- /(complex_t(0.0, 2.0)*value);
+ return std::sin(value)/value;
}
inline complex_t MathFunctions::Laue(complex_t value, size_t N) // Exp(iNx/2)*Sin((N+1)x)/Sin(x)
diff --git a/Core/Tools/src/MathFunctions.cpp b/Core/Tools/src/MathFunctions.cpp
index d69d26bbd55762cc12a4e51e398f19c005d908b8..dffe15451b0966adb75882bc7d86e9f21d006850 100644
--- a/Core/Tools/src/MathFunctions.cpp
+++ b/Core/Tools/src/MathFunctions.cpp
@@ -96,81 +96,81 @@ std::vector<complex_t > MathFunctions::FastFourierTransform(const std::vector<do
return MathFunctions::FastFourierTransform(cdata, ftCase);
}
-complex_t MathFunctions::Bessel_J1(complex_t value)
-{
- complex_t z1,z2,cr,cp,cs,cp0,cq0,cp1,cq1,ct1,ct2,cu;
- complex_t result;
- double a0;
- int k,kz;
-
- static complex_t czero(0.0, 0.0);
- static complex_t cone(1.0, 0.0);
- static double a1[] = {
- 0.1171875,
- -0.1441955566406250,
- 0.6765925884246826,
- -6.883914268109947,
- 1.215978918765359e2,
- -3.302272294480852e3,
- 1.276412726461746e5,
- -6.656367718817688e6,
- 4.502786003050393e8,
- -3.833857520742790e10,
- 4.011838599133198e12,
- -5.060568503314727e14,
- 7.572616461117958e16,
- -1.326257285320556e19};
- static double b1[] = {
- -0.1025390625,
- 0.2775764465332031,
- -1.993531733751297,
- 2.724882731126854e1,
- -6.038440767050702e2,
- 1.971837591223663e4,
- -8.902978767070678e5,
- 5.310411010968522e7,
- -4.043620325107754e9,
- 3.827011346598605e11,
- -4.406481417852278e13,
- 6.065091351222699e15,
- -9.833883876590679e17,
- 1.855045211579828e20};
-
- a0 = abs(value);
- z2 = value*value;
- z1 = value;
- if (a0 == 0.0) {
- result = czero;
- return result;
- }
- if (real(value) < 0.0) z1 = -value;
- if (a0 <= 12.0) {
- result = cone;
- cr = cone;
- for (k=1;k<=40;k++) {
- cr *= -0.25*z2/(k*(k+1.0));
- result += cr;
- if (abs(cr) < abs(result)*Numeric::double_epsilon) break;
- }
- result *= 0.5*z1;
- return result;
- }
- else {
- if (a0 >= 50.0) kz = 8; // can be changed to 10
- else if (a0 >= 35.0) kz = 10; // " " " 12
- else kz = 12; // " " " 14
- cp1 = cone;
- for (k=0;k<kz;k++) {
- cp1 += a1[k]*pow(z1,-2.0*k-2.0);
- }
- cq1 = 0.375/z1;
- for (k=0;k<kz;k++) {
- cq1 += b1[k]*pow(z1,-2.0*k-3.0);
- }
- result = cu*(cp1*cos(ct2)-cq1*sin(ct2));
- return result;
- }
-}
+//complex_t MathFunctions::Bessel_J1(complex_t value)
+//{
+// complex_t z1,z2,cr,cp,cs,cp0,cq0,cp1,cq1,ct1,ct2,cu;
+// complex_t result;
+// double a0;
+// int k,kz;
+//
+// static complex_t czero(0.0, 0.0);
+// static complex_t cone(1.0, 0.0);
+// static double a1[] = {
+// 0.1171875,
+// -0.1441955566406250,
+// 0.6765925884246826,
+// -6.883914268109947,
+// 1.215978918765359e2,
+// -3.302272294480852e3,
+// 1.276412726461746e5,
+// -6.656367718817688e6,
+// 4.502786003050393e8,
+// -3.833857520742790e10,
+// 4.011838599133198e12,
+// -5.060568503314727e14,
+// 7.572616461117958e16,
+// -1.326257285320556e19};
+// static double b1[] = {
+// -0.1025390625,
+// 0.2775764465332031,
+// -1.993531733751297,
+// 2.724882731126854e1,
+// -6.038440767050702e2,
+// 1.971837591223663e4,
+// -8.902978767070678e5,
+// 5.310411010968522e7,
+// -4.043620325107754e9,
+// 3.827011346598605e11,
+// -4.406481417852278e13,
+// 6.065091351222699e15,
+// -9.833883876590679e17,
+// 1.855045211579828e20};
+//
+// a0 = abs(value);
+// z2 = value*value;
+// z1 = value;
+// if (a0 == 0.0) {
+// result = czero;
+// return result;
+// }
+// if (real(value) < 0.0) z1 = -value;
+// if (a0 <= 12.0) {
+// result = cone;
+// cr = cone;
+// for (k=1;k<=40;k++) {
+// cr *= -0.25*z2/(k*(k+1.0));
+// result += cr;
+// if (abs(cr) < abs(result)*Numeric::double_epsilon) break;
+// }
+// result *= 0.5*z1;
+// return result;
+// }
+// else {
+// if (a0 >= 50.0) kz = 8; // can be changed to 10
+// else if (a0 >= 35.0) kz = 10; // " " " 12
+// else kz = 12; // " " " 14
+// cp1 = cone;
+// for (k=0;k<kz;k++) {
+// cp1 += a1[k]*pow(z1,-2.0*k-2.0);
+// }
+// cq1 = 0.375/z1;
+// for (k=0;k<kz;k++) {
+// cq1 += b1[k]*pow(z1,-2.0*k-3.0);
+// }
+// result = cu*(cp1*cos(ct2)-cq1*sin(ct2));
+// return result;
+// }
+//}
/* ************************************************************************* */
// convolution of two real vectors of equal size