From 5b8619b50d177af3710ce15c354d1c9dd1eb6e9a Mon Sep 17 00:00:00 2001
From: Walter Van Herck <w.van.herck@fz-juelich.de>
Date: Fri, 31 Aug 2012 17:59:38 +0200
Subject: [PATCH] Diffuse terms for nano particle crystal (1)

---
 App/src/TestMesoCrystal.cpp                   | 14 +++++++----
 .../inc/FormFactorDecoratorFactor.h           | 12 ++++++----
 .../src/DecouplingApproximationStrategy.cpp   |  2 +-
 ...LocalMonodisperseApproximationStrategy.cpp |  3 ++-
 Core/Samples/inc/IFormFactor.h                |  1 -
 .../inc/NanoParticleCrystalFormFactor.h       |  3 +++
 .../src/NanoParticleCrystalFormFactor.cpp     | 24 +++++++++++++++++++
 Core/Tools/inc/MathFunctions.h                | 13 ++++++++++
 8 files changed, 60 insertions(+), 12 deletions(-)

diff --git a/App/src/TestMesoCrystal.cpp b/App/src/TestMesoCrystal.cpp
index 7213bb049f7..e4d98743b4e 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_width = 600*Units::nanometer;
+    double meso_radius = 300*Units::nanometer;
     double surface_filling_ratio = 0.25;
-    double surface_density = surface_filling_ratio/meso_width/meso_width;
+    double surface_density = surface_filling_ratio/M_PI/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);
-    FormFactorGauss ff_meso(0.2*Units::micrometer, meso_width);
+    FormFactorCylinder 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
@@ -101,6 +101,7 @@ void TestMesoCrystal::initializeSample()
     size_t n_alpha_rotation_steps = 1;
     size_t n_np_size_steps = 1;
     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;
     double alpha_start = - (n_alpha_rotation_steps/2)*alpha_step;
     double np_size_step = 0.3*Units::nanometer/n_np_size_steps;
@@ -109,8 +110,8 @@ void TestMesoCrystal::initializeSample()
         for (size_t j=0; j<n_alpha_rotation_steps; ++j) {
             for (size_t k=0; k<n_np_size_steps; ++k) {
                 double R = np_size_start + k*np_size_step;
-                Geometry::RotateZ3D transform1(i*phi_step);
-                Geometry::RotateY3D transform2(alpha_start+j*alpha_step);
+                Geometry::RotateZ3D transform1(phi_start + i*phi_step);
+                Geometry::RotateY3D transform2(alpha_start + j*alpha_step);
                 Geometry::Transform3D *p_total_transform = new Geometry::Transform3D(transform1*transform2);
                 particle_decoration.addNanoParticle(createMesoCrystal(R, n_particle_adapted, &ff_meso), p_total_transform, 0.2*Units::micrometer);
     //            particle_decoration.addNanoParticle(meso2, transform1, 0.2*Units::micrometer, 0.5);
@@ -148,5 +149,8 @@ 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 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/Algorithms/inc/FormFactorDecoratorFactor.h b/Core/Algorithms/inc/FormFactorDecoratorFactor.h
index 7bb5d8e82f4..1cf3e21eac4 100644
--- a/Core/Algorithms/inc/FormFactorDecoratorFactor.h
+++ b/Core/Algorithms/inc/FormFactorDecoratorFactor.h
@@ -25,6 +25,8 @@ public:
     virtual ~FormFactorDecoratorFactor();
 
     virtual complex_t evaluate(cvector_t k_i, cvector_t k_f, double alpha_i, double alpha_f) const;
+    virtual double evaluateDiffuse(cvector_t k_i, cvector_t k_f, double alpha_i, double alpha_f) const;
+
     virtual int getNumberOfStochasticParameters() const;
 
 protected:
@@ -54,13 +56,15 @@ inline complex_t FormFactorDecoratorFactor::evaluate(cvector_t k_i,
     return m_factor*mp_form_factor->evaluate(k_i, k_f, alpha_i, alpha_f);
 }
 
+inline double FormFactorDecoratorFactor::evaluateDiffuse(cvector_t k_i,
+        cvector_t k_f, double alpha_i, double alpha_f) const
+{
+    return std::norm(m_factor)*mp_form_factor->evaluateDiffuse(k_i, k_f, alpha_i, alpha_f);
+}
+
 inline int FormFactorDecoratorFactor::getNumberOfStochasticParameters() const
 {
     return mp_form_factor->getNumberOfStochasticParameters();
 }
 
-
-
-
-
 #endif /* FORMFACTORDECORATORFACTOR_H_ */
diff --git a/Core/Algorithms/src/DecouplingApproximationStrategy.cpp b/Core/Algorithms/src/DecouplingApproximationStrategy.cpp
index 06f1c9679d7..265fa86f242 100644
--- a/Core/Algorithms/src/DecouplingApproximationStrategy.cpp
+++ b/Core/Algorithms/src/DecouplingApproximationStrategy.cpp
@@ -25,7 +25,7 @@ double DecouplingApproximationStrategy::evaluate(cvector_t k_i,
         complex_t ff = m_form_factors[i]->evaluate(k_i, k_f, alpha_i, alpha_f);
         double fraction = m_fractions[i];
         amplitude += fraction*ff;
-        intensity += fraction*(std::norm(ff) + m_form_factors[i]->evaluateDiffuse(k_i, k_f, alpha_i, alpha_f));
+        intensity += fraction*(std::norm(ff)); // + m_form_factors[i]->evaluateDiffuse(k_i, k_f, alpha_i, alpha_f));
     }
     double amplitude_norm = std::norm(amplitude);
     double itf_function = m_interference_functions[0]->evaluate(k_i-k_f);
diff --git a/Core/Algorithms/src/LocalMonodisperseApproximationStrategy.cpp b/Core/Algorithms/src/LocalMonodisperseApproximationStrategy.cpp
index 6126de9a32f..4b44010fc31 100644
--- a/Core/Algorithms/src/LocalMonodisperseApproximationStrategy.cpp
+++ b/Core/Algorithms/src/LocalMonodisperseApproximationStrategy.cpp
@@ -24,7 +24,8 @@ double LocalMonodisperseApproximationStrategy::evaluate(
         complex_t ff = m_form_factors[i]->evaluate(k_i, k_f, alpha_i, alpha_f);
         double itf_function = m_interference_functions[i]->evaluate(k_i-k_f);
         double fraction = m_fractions[i];
-        intensity += fraction*itf_function*std::norm(ff);
+        intensity += fraction*(itf_function*std::norm(ff)
+                                + m_form_factors[i]->evaluateDiffuse(k_i, k_f, alpha_i, alpha_f));
     }
     return intensity;
 }
diff --git a/Core/Samples/inc/IFormFactor.h b/Core/Samples/inc/IFormFactor.h
index dc0eeccbeef..67859245f24 100644
--- a/Core/Samples/inc/IFormFactor.h
+++ b/Core/Samples/inc/IFormFactor.h
@@ -55,7 +55,6 @@ inline double IFormFactor::evaluateDiffuse(cvector_t k_i, cvector_t k_f, double
     return 0.0;
 }
 
-
 class IFormFactorDecorator : public IFormFactor
 {
 public:
diff --git a/Core/Samples/inc/NanoParticleCrystalFormFactor.h b/Core/Samples/inc/NanoParticleCrystalFormFactor.h
index 45962f300d7..408f636e536 100644
--- a/Core/Samples/inc/NanoParticleCrystalFormFactor.h
+++ b/Core/Samples/inc/NanoParticleCrystalFormFactor.h
@@ -33,16 +33,19 @@ public:
 
     virtual void setAmbientRefractiveIndex(complex_t refractive_index);
 
+//    virtual double evaluateDiffuse(cvector_t k_i, cvector_t k_f, double alpha_i, double alpha_f) const;
 protected:
     virtual complex_t evaluate_for_q(cvector_t q) const;
 private:
     void calculateLargestReciprocalDistance();
+//    void initializeDiffuseNanoparticleFormfactors();
     Lattice m_lattice;
     NanoParticle *mp_nano_particle;
     IFormFactor *mp_basis_form_factor;
     IFormFactor *mp_meso_form_factor;
     complex_t m_ambient_refractive_index;
     double m_max_rec_length;
+//    std::vector<IFormFactor *> m_diffuse_nanoparticle_ffs;
 };
 
 
diff --git a/Core/Samples/src/NanoParticleCrystalFormFactor.cpp b/Core/Samples/src/NanoParticleCrystalFormFactor.cpp
index 24fc901077c..a1d533babce 100644
--- a/Core/Samples/src/NanoParticleCrystalFormFactor.cpp
+++ b/Core/Samples/src/NanoParticleCrystalFormFactor.cpp
@@ -62,6 +62,16 @@ complex_t NanoParticleCrystalFormFactor::evaluate_for_q(cvector_t q) const
     return 8.0*pi3*result/volume;
 }
 
+//double NanoParticleCrystalFormFactor::evaluateDiffuse(cvector_t k_i,
+//        cvector_t k_f, double alpha_i, double alpha_f) const
+//{
+//    // TODO: remove fixed params
+//    cvector_t k_zero;
+//    double nbr_nanoparticles = mp_meso_form_factor->evaluate(k_zero, k_zero, 0.0, 0.0)/m_lattice.getVolume();
+//
+//
+//}
+
 void NanoParticleCrystalFormFactor::calculateLargestReciprocalDistance()
 {
     kvector_t a1 = m_lattice.getBasisVectorA();
@@ -71,3 +81,17 @@ void NanoParticleCrystalFormFactor::calculateLargestReciprocalDistance()
     m_max_rec_length = std::max(M_PI/a1.mag(), M_PI/a2.mag());
     m_max_rec_length = std::max(m_max_rec_length, M_PI/a3.mag());
 }
+
+//void NanoParticleCrystalFormFactor::initializeDiffuseNanoparticleFormfactors()
+//{
+//    // TODO: remove fixed params
+//    NanoParticle *p_nano = mp_nano_particle->clone();
+//    double heigth = 0.2*Units::micrometer;
+//    size_t nbr_heigths = 10;
+//    size_t nbr_shapes = 5;
+//    for (size_t i=0; i<nbr_shapes; ++i) {
+//        for (size_t j=0; j< nbr_heigths; ++j) {
+//
+//        }
+//    }
+//}
diff --git a/Core/Tools/inc/MathFunctions.h b/Core/Tools/inc/MathFunctions.h
index 376cc478238..7670c9d574f 100644
--- a/Core/Tools/inc/MathFunctions.h
+++ b/Core/Tools/inc/MathFunctions.h
@@ -44,6 +44,8 @@ double Sinc(double value);
 
 complex_t Sinc(complex_t value);
 
+complex_t Laue(complex_t value, size_t N);
+
 enum TransformCase { ForwardFFT, BackwardFFT };
 std::vector<complex_t > FastFourierTransform(const std::vector<complex_t > &data, TransformCase tcase);
 
@@ -83,4 +85,15 @@ inline complex_t MathFunctions::Sinc(complex_t value)  // Sin(x)/x
     		/(complex_t(0.0, 2.0)*value);
 }
 
+inline complex_t MathFunctions::Laue(complex_t value, size_t N) // Exp(iNx/2)*Sin((N+1)x)/Sin(x)
+{
+    if (N==0) {
+        return complex_t(1.0, 0.0);
+    }
+    if(std::abs(value)<Numeric::double_epsilon) {
+        return complex_t(N+1.0, 0.0);
+    }
+    return std::exp(complex_t(0.0, 1.0)*value*(double)N/2.0)*std::sin(value*(N+1.0)/2.0)/std::sin(value/2.0);
+}
+
 #endif // MATHFUNCTIONS_H
-- 
GitLab