Refactored interference function creation and added isgisaxs6 sample

App/src/StandardSamples.cpp

@@ -13,6 +13,7 @@
 #include "Crystal.h"
 #include "MesoCrystal.h"
 #include "InterferenceFunction1DParaCrystal.h"
+#include "InterferenceFunction2DLattice.h"
 #include "InterferenceFunction2DParaCrystal.h"
 #include "FormFactorWeighted.h"
 #include "StochasticGaussian.h"
@@ -574,9 +575,27 @@ ISample *StandardSamples::IsGISAXS6_lattice()
     air_layer.setMaterial(p_air_material);
     Layer substrate_layer;
     substrate_layer.setMaterial(p_substrate_material);
-    IInterferenceFunction *p_interference_function = new InterferenceFunction1DParaCrystal(20.0*Units::nanometer,7*Units::nanometer, 1e3*Units::nanometer);
-    ParticleDecoration particle_decoration( new Particle(n_particle, new FormFactorCylinder(5*Units::nanometer, 5*Units::nanometer)));
+    Lattice2DIFParameters lattice_params = {
+            10.0*Units::nanometer,       // L1
+            10.0*Units::nanometer,       // L2
+            90.0*Units::degree,          // lattice angle
+            0.0*Units::degree,           // lattice orientation
+            20000.0*Units::nanometer,    // domain size 1
+            20000.0*Units::nanometer,    // domain size 2
+            300.0*Units::nanometer/2.0/M_PI, // correlation length 1
+            100.0*Units::nanometer/2.0/M_PI  // correlation length 2
+    };
+    InterferenceFunction2DLattice *p_interference_function = new InterferenceFunction2DLattice(lattice_params);
+    FTDistribution2DCauchy pdf(300.0*Units::nanometer/2.0/M_PI, 100.0*Units::nanometer/2.0/M_PI);
+    p_interference_function->setProbabilityDistribution(pdf);
+
+    FormFactorCylinder ff_cyl(5.0*Units::nanometer, 5.0*Units::nanometer);
+    ParticleDecoration particle_decoration( new Particle(n_particle, ff_cyl.clone()));
     particle_decoration.addInterferenceFunction(p_interference_function);
+    kvector_t position_2(5.0*Units::nanometer, 5.0*Units::nanometer, 0.0);
+    Particle *p_particle_2 = new Particle(n_particle, new FormFactorDecoratorPositionFactor(ff_cyl, position_2));
+    particle_decoration.addParticle( p_particle_2, 0.0, 1.0 );
+    particle_decoration.addInterferenceFunction(p_interference_function->clone());
     LayerDecorator air_layer_decorator(air_layer, particle_decoration);
 
     p_multi_layer->addLayer(air_layer_decorator);

Core/Algorithms/inc/IInterferenceFunctionStrategy.h

@@ -19,29 +19,28 @@
 #include "IInterferenceFunction.h"
 #include "Bin.h"
 #include "SafePointerVector.h"
+#include "StrategyBuilder.h"
 
 #include <vector>
 
 class IInterferenceFunctionStrategy
 {
 public:
-    virtual ~IInterferenceFunctionStrategy();
+    virtual ~IInterferenceFunctionStrategy() {};
     virtual void init(const std::vector<IFormFactor *> &form_factors,
             const std::vector<double> &fractions,
             const std::vector<IInterferenceFunction *> &interference_functions);
+    virtual void init(const SafePointerVector<FormFactorInfo> &form_factor_infos,
+            const SafePointerVector<IInterferenceFunction> &ifs);
     virtual double evaluate(const cvector_t &k_i, const Bin1DCVector &k_f_bin,
             double alpha_i, double alpha_f) const=0;
 protected:
-    void deleteVectors();
-    SafePointerVector<IFormFactor> m_form_factors; //!< Includes Scattering Length Density
-    std::vector<double> m_fractions;
-    SafePointerVector<IInterferenceFunction> m_interference_functions;
+    SafePointerVector<FormFactorInfo> m_ff_infos; //!< Form factor info
+    SafePointerVector<IFormFactor> m_form_factors; //TODO: remove
+    std::vector<double> m_fractions; //TODO: remove
+    SafePointerVector<IInterferenceFunction> m_ifs; //!< Interference functions
 };
 
-inline IInterferenceFunctionStrategy::~IInterferenceFunctionStrategy()
-{
-}
-
 inline void IInterferenceFunctionStrategy::init(
         const std::vector<IFormFactor*>& form_factors,
         const std::vector<double>& fractions,
@@ -52,10 +51,18 @@ inline void IInterferenceFunctionStrategy::init(
     for (size_t i=0; i<form_factors.size(); ++i) {
         m_form_factors.push_back(form_factors[i]->clone());
     }
-    m_interference_functions.clear();
+    m_ifs.clear();
     for (size_t i=0; i<interference_functions.size(); ++i) {
-        m_interference_functions.push_back(interference_functions[i]->clone());
+        m_ifs.push_back(interference_functions[i]->clone());
     }
 }
 
+inline void IInterferenceFunctionStrategy::init(
+        const SafePointerVector<FormFactorInfo> &form_factor_infos,
+        const SafePointerVector<IInterferenceFunction> &ifs)
+{
+    m_ff_infos = form_factor_infos;
+    m_ifs = ifs;
+}
+
 #endif /* IINTERFERENCEFUNCTIONSTRATEGY_H_ */

Core/Algorithms/inc/SimulationParameters.h

+#ifndef SIMULATIONPARAMETERS_H_
+#define SIMULATIONPARAMETERS_H_
+// ********************************************************************
+// * The BornAgain project                                            *
+// * Simulation of neutron and x-ray scattering at grazing incidence  *
+// *                                                                  *
+// * LICENSE AND DISCLAIMER                                           *
+// * Lorem ipsum dolor sit amet, consectetur adipiscing elit.  Mauris *
+// * eget quam orci. Quisque  porta  varius  dui,  quis  posuere nibh *
+// * mollis quis. Mauris commodo rhoncus porttitor.                   *
+// ********************************************************************
+//! @file   SimulationParameters.h
+//! @brief  Definition of SimulationParameters class
+//! @author Scientific Computing Group at FRM II
+//! @date   Jan 24, 2013
+
+//- -------------------------------------------------------------------
+//! @class SimulationParameters
+//! @brief Definition of parameter class collecting the different
+//! options for simulation
+//- -------------------------------------------------------------------
+class SimulationParameters
+{
+public:
+    SimulationParameters();
+    enum EFramework { DWBA, BA } me_framework;
+    enum EInterferenceApproximation { DA, LMA, SSCA } me_if_approx;
+    enum ELatticeType { NONE, LATTICE, PARA1D, PARA1DFINITE } me_lattice_type;
+};
+
+inline SimulationParameters::SimulationParameters()
+: me_framework(DWBA)
+, me_if_approx(DA)
+, me_lattice_type(NONE)
+{
+}
+
+#endif /* SIMULATIONPARAMETERS_H_ */

Core/Algorithms/inc/StrategyBuilder.h

+#ifndef STRATEGYBUILDER_H_
+#define STRATEGYBUILDER_H_
+// ********************************************************************
+// * The BornAgain project                                            *
+// * Simulation of neutron and x-ray scattering at grazing incidence  *
+// *                                                                  *
+// * LICENSE AND DISCLAIMER                                           *
+// * Lorem ipsum dolor sit amet, consectetur adipiscing elit.  Mauris *
+// * eget quam orci. Quisque  porta  varius  dui,  quis  posuere nibh *
+// * mollis quis. Mauris commodo rhoncus porttitor.                   *
+// ********************************************************************
+//! @file   LayerDecoratorStrategyBuilder.h
+//! @brief  Definition of LayerDecoratorStrategyBuilder class
+//! @author Scientific Computing Group at FRM II
+//! @date   Jan 24, 2013
+
+#include "SimulationParameters.h"
+#include "SafePointerVector.h"
+
+struct FormFactorInfo;
+class IInterferenceFunctionStrategy;
+class LayerDecorator;
+class Experiment;
+class IDoubleToPairOfComplexMap;
+class ParticleInfo;
+
+//- -------------------------------------------------------------------
+//! @class LayerDecoratorStrategyBuilder
+//! @brief Definition of methods to generate a simulation strategy from
+//! a LayerDecorator and SimulationParameters
+//- -------------------------------------------------------------------
+class LayerDecoratorStrategyBuilder
+{
+public:
+    LayerDecoratorStrategyBuilder(const LayerDecorator &decorated_layer,
+            const Experiment &experiment, const SimulationParameters &sim_params);
+    virtual ~LayerDecoratorStrategyBuilder();
+
+    //! set R and T coefficient map for DWBA simulation
+    void setReflectionTransmissionFunction(const IDoubleToPairOfComplexMap &rt_map);
+
+    //! create a strategy object which is able to calculate the scattering for fixed k_f
+    virtual IInterferenceFunctionStrategy *createStrategy();
+protected:
+    LayerDecorator *mp_layer_decorator;         //!< decorated layer
+    Experiment *mp_experiment;                  //!< experiment
+    SimulationParameters m_sim_params;          //!< simulation parameters
+    IDoubleToPairOfComplexMap *mp_RT_function;  //!< R and T coefficients for DWBA
+private:
+    //! collect the formfactor info of all particles in the decoration and decorate
+    //! these for DWBA when needed
+    void collectFormFactorInfos();
+    //! collect the interference functions
+    void collectInterferenceFunctions();
+    //! retrieve wavelength from experiment
+    double getWavelength();
+    //! create formfactor info for single particle
+    FormFactorInfo *createFormFactorInfo(const ParticleInfo *p_particle_info,
+            complex_t n_ambient_refractive_index, complex_t factor) const;
+
+    SafePointerVector<FormFactorInfo> m_ff_infos;
+    SafePointerVector<IInterferenceFunction> m_ifs;
+};
+
+struct FormFactorInfo : public ICloneable
+{
+    FormFactorInfo();
+    ~FormFactorInfo();
+    virtual FormFactorInfo *clone() const;
+    IFormFactor *mp_ff;
+    double m_pos_x, m_pos_y;
+    double m_abundance;
+};
+
+#endif /* STRATEGYBUILDER_H_ */

Core/Algorithms/src/DecouplingApproximationStrategy.cpp

@@ -28,7 +28,7 @@ double DecouplingApproximationStrategy::evaluate(const cvector_t& k_i,
         intensity += fraction*(std::norm(ff));
     }
     double amplitude_norm = std::norm(amplitude);
-    double itf_function = m_interference_functions[0]->evaluate(k_i-k_f_bin.getMidPoint());
+    double itf_function = m_ifs[0]->evaluate(k_i-k_f_bin.getMidPoint());
     return intensity + amplitude_norm*(itf_function-1.0);
 }
 
@@ -36,7 +36,7 @@ bool DecouplingApproximationStrategy::checkVectorSizes() const
 {
     size_t n_ffs = m_form_factors.size();
     size_t n_frs = m_fractions.size();
-    size_t n_ifs = m_interference_functions.size();
+    size_t n_ifs = m_ifs.size();
     return (n_ffs==n_frs && n_ifs==1);
 }
 

Core/Algorithms/src/IsGISAXSMorphologyFileStrategy.cpp

@@ -88,7 +88,7 @@ bool IsGISAXSMorphologyFileStrategy::checkVectorSizes()
 {
     size_t n_ffs = m_form_factors.size();
     size_t n_frs = m_fractions.size();
-    size_t n_ifs = m_interference_functions.size();
+    size_t n_ifs = m_ifs.size();
     return (n_ffs==n_frs && n_ifs==0);
 }
 

Core/Algorithms/src/LocalMonodisperseApproximationStrategy.cpp

@@ -22,7 +22,7 @@ double LocalMonodisperseApproximationStrategy::evaluate(const cvector_t& k_i,
     double intensity = 0.0;
     for (size_t i=0; i<m_form_factors.size(); ++i) {
         complex_t ff = m_form_factors[i]->evaluate(k_i, k_f_bin, alpha_i, alpha_f);
-        double itf_function = m_interference_functions[i]->evaluate(k_i-k_f_bin.getMidPoint());
+        double itf_function = m_ifs[i]->evaluate(k_i-k_f_bin.getMidPoint());
         double fraction = m_fractions[i];
         intensity += fraction*(itf_function*std::norm(ff));
     }
@@ -33,6 +33,6 @@ bool LocalMonodisperseApproximationStrategy::checkVectorSizes()
 {
     size_t n_ffs = m_form_factors.size();
     size_t n_frs = m_fractions.size();
-    size_t n_ifs = m_interference_functions.size();
+    size_t n_ifs = m_ifs.size();
     return (n_ffs==n_frs && n_ifs==n_ffs);
 }

Core/Algorithms/src/SizeSpacingCorrelationApproximationStrategy.cpp

@@ -42,7 +42,7 @@ bool SizeSpacingCorrelationApproximationStrategy::checkVectorSizes() const
 {
     size_t n_ffs = m_form_factors.size();
     size_t n_frs = m_fractions.size();
-    size_t n_ifs = m_interference_functions.size();
+    size_t n_ifs = m_ifs.size();
     return (n_ffs==n_frs && n_ifs==1);
 }
 
@@ -76,7 +76,7 @@ complex_t SizeSpacingCorrelationApproximationStrategy::getCharacteristicDistribu
         double qp) const
 {
     const InterferenceFunction1DParaCrystal *p_iff = dynamic_cast<const InterferenceFunction1DParaCrystal *>(
-            m_interference_functions[0]);
+            m_ifs[0]);
     if (p_iff==0) {
         throw ClassInitializationException("Wrong interference function for SSCA");
     }

Core/Algorithms/src/StrategyBuilder.cpp

+#include "StrategyBuilder.h"
+
+#include <cmath>
+
+LayerDecoratorStrategyBuilder::LayerDecoratorStrategyBuilder(
+        const LayerDecorator &decorated_layer, const Experiment &experiment,
+        const SimulationParameters &sim_params)
+: mp_layer_decorator(decorated_layer.clone())
+, mp_experiment(experiment.clone())
+, m_sim_params(sim_params)
+, mp_RT_function(0)
+{
+}
+
+LayerDecoratorStrategyBuilder::~LayerDecoratorStrategyBuilder()
+{
+    delete mp_layer_decorator;
+    delete mp_experiment;
+    delete mp_RT_function;
+}
+
+void LayerDecoratorStrategyBuilder::setReflectionTransmissionFunction(
+        const IDoubleToPairOfComplexMap& rt_map)
+{
+    if (mp_RT_function != &rt_map) {
+        delete mp_RT_function;
+        mp_RT_function = rt_map.clone();
+    }
+}
+
+IInterferenceFunctionStrategy* LayerDecoratorStrategyBuilder::createStrategy()
+{
+    collectFormFactorInfos();
+    collectInterferenceFunctions();
+    size_t n_particles = m_ff_infos.size();
+    size_t n_ifs = m_ifs.size();
+    IInterferenceFunctionStrategy *p_result(0);
+    switch (m_sim_params.me_if_approx)
+    {
+    case SimulationParameters::DA:
+        p_result = new DecouplingApproximationStrategy();
+        break;
+    case SimulationParameters::LMA:
+        p_result = new LocalMonodisperseApproximationStrategy();
+        break;
+    case SimulationParameters::SSCA:
+        if (n_ifs<2) {
+            throw Exceptions::ClassInitializationException(
+                    "SSCA requires an interference function");
+        }
+        double kappa = m_ifs[0]->getKappa();
+        if (kappa<=0.0) {
+            throw Exceptions::ClassInitializationException(
+                    "SSCA requires a strictly positive coupling value");
+        }
+        p_result = new SizeSpacingCorrelationApproximationStrategy(kappa);
+        break;
+    default:
+        throw Exceptions::ClassInitializationException(
+                "Unknown interference function approximation");
+    }
+    if (!p_result) {
+        throw Exceptions::ClassInitializationException(
+                "Could not create appropriate strategy");
+    }
+    p_result->init(m_ff_infos, m_ifs);
+    return p_result;
+}
+
+void LayerDecoratorStrategyBuilder::collectFormFactorInfos()
+{
+    m_ff_infos.clear();
+    const IDecoration *p_decoration = mp_layer_decorator->getDecoration();
+    complex_t n_layer = mp_layer_decorator->getRefractiveIndex();
+    double wavelength = getWavelength();
+    complex_t wavevector_scattering_factor = M_PI/wavelength/wavelength;
+    size_t number_of_particles = p_decoration->getNumberOfParticles();
+    for (size_t particle_index=0; particle_index<number_of_particles; ++particle_index) {
+        const ParticleInfo *p_particle_info = p_decoration->getParticleInfo(particle_index);
+        FormFactorInfo *p_ff_info = createFormFactorInfo(p_particle_info, n_layer,
+                wavevector_scattering_factor);
+        m_ff_infos.push_back(p_ff_info);
+    }
+    return;
+}
+
+void LayerDecoratorStrategyBuilder::collectInterferenceFunctions()
+{
+    m_ifs.clear();
+    if (mp_layer_decorator->getDecoration()->getNumberOfInterferenceFunctions()) {
+        m_ifs = mp_layer_decorator->getDecoration()->getInterferenceFunctions();
+    }
+    else m_ifs.push_back(new InterferenceFunctionNone);
+}
+
+double LayerDecoratorStrategyBuilder::getWavelength()
+{
+    cvector_t ki = mp_experiment->getBeam().getCentralK();
+    kvector_t ki_real(ki.x().real(), ki.y().real(), ki.z().real());
+    return 2.0*M_PI/ki_real.mag();
+}
+
+FormFactorInfo *LayerDecoratorStrategyBuilder::createFormFactorInfo(
+        const ParticleInfo *p_particle_info, complex_t n_ambient_refractive_index,
+        complex_t factor) const
+{
+    FormFactorInfo *p_result = new FormFactorInfo;
+    Particle *p_particle_clone = p_particle_info->getParticle()->clone();
+    const Geometry::Transform3D *p_transform = p_particle_info->getTransform3D();
+
+    // formfactor
+    p_particle_clone->setAmbientRefractiveIndex(n_ambient_refractive_index);
+    IFormFactor *ff_particle = p_particle_clone->createFormFactor();
+    delete p_particle_clone;
+    IFormFactor *ff_transformed(0);
+    if(p_transform) {
+        ff_transformed = new FormFactorDecoratorTransformation(ff_particle, new Geometry::Transform3D(*p_transform));
+    } else{
+        ff_transformed = ff_particle;
+    }
+    IFormFactor *p_ff_framework(0);
+    switch (m_sim_params.me_framework)
+    {
+    case SimulationParameters::BA:    // Born Approximation
+        p_ff_framework = ff_transformed;
+        break;
+    case SimulationParameters::DWBA:  // Distorted Wave Born Approximation
+        if (mp_RT_function==0) {
+            throw Exceptions::ClassInitializationException(
+                    "R and T coefficients are necessary for DWBA");
+        }
+        double depth = p_particle_info->getDepth();
+        FormFactorDWBAConstZ *p_dwba_ff = new FormFactorDWBAConstZ(ff_transformed, depth);
+        p_dwba_ff->setReflectionTransmissionFunction(*mp_RT_function);
+        p_ff_framework = p_dwba_ff;
+        break;
+    default:
+        throw Exceptions::RuntimeErrorException("Framework must be BA or DWBA");
+    }
+    FormFactorDecoratorFactor *p_ff = new FormFactorDecoratorFactor(p_ff_framework, factor);
+    p_result->mp_ff = p_ff;
+    // Other info (position and abundance
+    const PositionParticleInfo *p_pos_particle_info = dynamic_cast<const PositionParticleInfo *>(p_particle_info);
+    if (p_pos_particle_info) {
+        kvector_t position = p_pos_particle_info->getPosition();
+        p_result->m_pos_x = position.x();
+        p_result->m_pos_y = position.y();
+    }
+    p_result->m_abundance = p_particle_info->getAbundance();
+    return p_result;
+}
+
+FormFactorInfo::FormFactorInfo()
+: mp_ff(0)
+, m_pos_x(0.0)
+, m_pos_y(0.0)
+, m_abundance(0.0)
+{
+}
+
+FormFactorInfo::~FormFactorInfo()
+{
+    delete mp_ff;
+}
+
+FormFactorInfo* FormFactorInfo::clone() const
+{
+    FormFactorInfo *p_result = new FormFactorInfo();
+    p_result->m_abundance = m_abundance;
+    p_result->m_pos_x = m_pos_x;
+    p_result->m_pos_y = m_pos_y;
+    p_result->mp_ff = mp_ff->clone();
+    return p_result;
+}

Core/Core.pro

@@ -76,6 +76,7 @@ SOURCES += \
     Samples/src/ICompositeSample.cpp \
     Samples/src/IMaterial.cpp \
     Samples/src/InterferenceFunction1DParaCrystal.cpp \
+    Samples/src/InterferenceFunction2DLattice.cpp \
     Samples/src/InterferenceFunction2DParaCrystal.cpp \
     Samples/src/ISample.cpp \
     Samples/src/IsGISAXSMorphologyFileDecoration.cpp \
@@ -163,7 +164,9 @@ HEADERS += \
     Algorithms/inc/MultiLayerRoughnessDWBASimulation.h \
     Algorithms/inc/OpticalFresnel.h \
     Algorithms/inc/ResolutionFunction2DSimple.h \
+    Algorithms/inc/SimulationParameters.h \
     Algorithms/inc/SizeSpacingCorrelationApproximationStrategy.h \
+    Algorithms/inc/StrategyBuilder.h \
     Algorithms/inc/ThreadInfo.h \
     \
     FormFactors/inc/FormFactorBox.h \
@@ -213,6 +216,7 @@ HEADERS += \
     Samples/inc/IInterferenceFunction.h \
     Samples/inc/IMaterial.h \
     Samples/inc/InterferenceFunction1DParaCrystal.h \
+    Samples/inc/InterferenceFunction2DLattice.h \
     Samples/inc/InterferenceFunction2DParaCrystal.h \
     Samples/inc/InterferenceFunctionNone.h \
     Samples/inc/IRoughness.h \
@@ -220,6 +224,7 @@ HEADERS += \
     Samples/inc/ISelectionRule.h \
     Samples/inc/IsGISAXSMorphologyFileDecoration.h \
     Samples/inc/Lattice.h \
+    Samples/inc/Lattice2DIFParameters.h \
     Samples/inc/LatticeBasis.h \
     Samples/inc/Layer.h \
     Samples/inc/LayerDecorator.h \

Core/Samples/inc/IDecoration.h

@@ -15,9 +15,12 @@
 //! @date   Jun 22, 2012
 
 #include "ICompositeSample.h"
-#include "IInterferenceFunctionStrategy.h"
 #include "IFormFactor.h"
 #include "ParticleInfo.h"
+#include "SafePointerVector.h"
+
+class IInterferenceFunctionStrategy;
+class IInterferenceFunction;
 
 class IDecoration : public ICompositeSample
 {
@@ -29,16 +32,22 @@ public:
     virtual IInterferenceFunctionStrategy *createStrategy(
             const std::vector<IFormFactor *> &form_factors) const=0;
 
-    /// Get number of particles
+    //! Get number of particles
     virtual size_t getNumberOfParticles() const=0;
 
-    /// get information about particle with index
+    //! get information about particle with index
     virtual const ParticleInfo *getParticleInfo(size_t index) const=0;
 
-    /// Get surface density of all particles
+    //! get number of interference functions
+    virtual size_t getNumberOfInterferenceFunctions() const { return 0; }
+
+    //! get interference functions
+    virtual SafePointerVector<IInterferenceFunction> getInterferenceFunctions() const=0;
+
+    //! get surface density of all particles
     double getTotalParticleSurfaceDensity() const { return m_total_particle_surface_density; }
 
-    /// Set surface density of all particles
+    //! set surface density of all particles
     void setTotalParticleSurfaceDensity(double surface_density) { m_total_particle_surface_density = surface_density; }
 
 private:

Core/Samples/inc/IInterferenceFunction.h

@@ -24,6 +24,7 @@ public:
 
     virtual double evaluate(const cvector_t &q) const=0;
 	virtual IInterferenceFunction *clone() const=0;
+    virtual double getKappa() const { return 0.0; }
 };
 
 

Core/Samples/inc/InterferenceFunction1DParaCrystal.h

@@ -28,7 +28,7 @@ public:
 	}
 
 	void setKappa(double kappa) { m_kappa = kappa; }
-	double getKappa() const { return m_kappa; }
+	virtual double getKappa() const { return m_kappa; }
     virtual double evaluate(const cvector_t &q) const;
     //TODO: replace these with strategy pattern for different algorithms
     complex_t FTGaussianCorrLength(double qpar) const;
@@ -43,6 +43,4 @@ private:
     virtual void init_parameters();
 };
 
-
-
 #endif /* INTERFERENCEFUNCTION1DPARACRYSTAL_H_ */

Core/Samples/inc/InterferenceFunction2DLattice.h

+#ifndef INTERFERENCEFUNCTION2DLATTICE_H_
+#define INTERFERENCEFUNCTION2DLATTICE_H_
+// ********************************************************************
+// * The BornAgain project                                            *
+// * Simulation of neutron and x-ray scattering at grazing incidence  *
+// *                                                                  *
+// * LICENSE AND DISCLAIMER                                           *
+// * Lorem ipsum dolor sit amet, consectetur adipiscing elit.  Mauris *
+// * eget quam orci. Quisque  porta  varius  dui,  quis  posuere nibh *
+// * mollis quis. Mauris commodo rhoncus porttitor.                   *
+// ********************************************************************
+//! @file   InterferenceFunction2DLattice.h
+//! @brief  Definition of InterferenceFunction2DLattice class
+//! @author Scientific Computing Group at FRM II
+//! @date   Jan 22, 2013
+
+#include "IInterferenceFunction.h"
+#include "Lattice2DIFParameters.h"
+#include "FTDistributions.h"
+
+class InterferenceFunction2DLattice : public IInterferenceFunction
+{
+public:
+    InterferenceFunction2DLattice(const Lattice2DIFParameters &lattice_params);
+    virtual ~InterferenceFunction2DLattice() {}
+
+    virtual InterferenceFunction2DLattice *clone() const {
+        InterferenceFunction2DLattice *p_clone = new InterferenceFunction2DLattice(m_lattice_params);
+        p_clone->setProbabilityDistribution(*mp_pdf);
+        return p_clone;
+    }
+
+    void setProbabilityDistribution(const IFTDistribution2D &pdf);
+
+    virtual double evaluate(const cvector_t &q) const;
+protected:
+    double interferenceAtOneRecLatticePoint(double qx, double qy) const;
+    void transformToPrincipalAxes(double qx, double qy, double gamma, double delta, double &q_pa_1, double &q_pa_2) const;
+    void calculateReciprocalVectorFraction(double qx, double qy, double &qx_frac, double &qy_frac) const;
+    Lattice2DIFParameters m_lattice_params;
+    IFTDistribution2D *mp_pdf;
+private:
+    //! initialize pool parameters, i.e. register some of class members for later access via parameter pool
+    virtual void init_parameters();
+    void initialize_rec_vectors();
+    double m_asx, m_asy;
+    double m_bsx, m_bsy;
+};
+
+#endif /* INTERFERENCEFUNCTION2DLATTICE_H_ */

Core/Samples/inc/IsGISAXSMorphologyFileDecoration.h

@@ -44,32 +44,42 @@ public:
     void addParticle(const Particle &particle, kvector_t position, double abundance=1.0);
     void addParticle(Particle *p_particle, kvector_t position, double abundance=1.0);
 
-    /// Add particle info
+    /// add particle info
     void addParticleInfo(const PositionParticleInfo &info);
 
-    /// Get number of particles
+    /// get number of particles
     virtual size_t getNumberOfParticles() const { return m_particles.size(); }
 
     /// get information about particle with index
     virtual const PositionParticleInfo *getParticleInfo(size_t index) const;
 
-    /// Get abundance fraction of particle with index
+    /// get abundance fraction of particle with index
     double getAbundanceFractionOfParticle(size_t index) const;
 
-    /// Add interference function
+    /// add interference function
     void addInterferenceFunction(IInterferenceFunction* p_interference_function);
     void addInterferenceFunction(const IInterferenceFunction &interference_function);
 
-    /// Get interference function with index
+    //! get number of interference functions
+    virtual size_t getNumberOfInterferenceFunctions() const {
+        return m_interference_functions.size();
+    }
+
+    //! get interference functions
+    virtual SafePointerVector<IInterferenceFunction> getInterferenceFunctions() const {
+        return m_interference_functions;
+    }
+
+    /// get interference function with index
     const IInterferenceFunction* getInterferenceFunction(size_t index) const;
 
-    /// Create interference function strategy
+    /// create interference function strategy
     IInterferenceFunctionStrategy *createStrategy(const std::vector<IFormFactor *> &form_factors) const;
 
-    /// Get surface density of all particles
+    /// get surface density of all particles
     double getTotalParticleSurfaceDensity() const { return m_total_particle_surface_density; }
 
-    /// Set surface density of all particles
+    /// set surface density of all particles
     void setTotalParticleSurfaceDensity(double surface_density) { m_total_particle_surface_density = surface_density; }
 private:
     /// copy constructor and assignment operator are hidden since there is a clone method

Core/Samples/inc/Lattice2DIFParameters.h

+#ifndef LATTICE2DIFPARAMETERS_H_
+#define LATTICE2DIFPARAMETERS_H_
+// ********************************************************************
+// * The BornAgain project                                            *
+// * Simulation of neutron and x-ray scattering at grazing incidence  *
+// *                                                                  *
+// * LICENSE AND DISCLAIMER                                           *
+// * Lorem ipsum dolor sit amet, consectetur adipiscing elit.  Mauris *
+// * eget quam orci. Quisque  porta  varius  dui,  quis  posuere nibh *
+// * mollis quis. Mauris commodo rhoncus porttitor.                   *
+// ********************************************************************
+//! @file   Lattice2DIFParameters.h
+//! @brief  Definition of Lattice2DIFParameters
+//! @author Scientific Computing Group at FRM II
+//! @date   Jan 22, 2013
+
+#include "IParameterized.h"
+
+struct Lattice2DIFParameters //: public IParameterized
+{
+    double m_length_1, m_length_2;
+    double m_angle;
+    double m_xi;
+    double m_domain_size_1, m_domain_size_2;
+    double m_corr_length_1, m_corr_length_2;
+};
+
+
+#endif /* LATTICE2DIFPARAMETERS_H_ */

Core/Samples/inc/ParticleDecoration.h

@@ -55,6 +55,16 @@ public:
     /// Get abundance fraction of particle with index
     double getAbundanceFractionOfParticle(size_t index) const;
 
+    //! get number of interference functions
+    virtual size_t getNumberOfInterferenceFunctions() const {
+        return m_interference_functions.size();
+    }
+
+    //! get interference functions
+    virtual SafePointerVector<IInterferenceFunction> getInterferenceFunctions() const {
+        return m_interference_functions;
+    }
+
     /// Add interference function
     void addInterferenceFunction(IInterferenceFunction* p_interference_function);
     void addInterferenceFunction(const IInterferenceFunction &interference_function);
@@ -82,9 +92,10 @@ private:
         registerChild(child);
     }
 
+    //TODO: replace with SafePointerVector
     std::vector<ParticleInfo *> m_particles;
     ///< Vector of the types of particles
-    std::vector<IInterferenceFunction *> m_interference_functions;
+    SafePointerVector<IInterferenceFunction> m_interference_functions;
     ///< Currently only a scalar interference function (instead of matrix)
     double m_total_abundance;
 };

Core/Samples/src/InterferenceFunction2DLattice.cpp

+#include "InterferenceFunction2DLattice.h"
+
+InterferenceFunction2DLattice::InterferenceFunction2DLattice(
+        const Lattice2DIFParameters& lattice_params)
+: m_lattice_params(lattice_params)
+, mp_pdf(0)
+{
+    setName("InterferenceFunction2DLattice");
+    init_parameters();
+    initialize_rec_vectors();
+}
+
+void InterferenceFunction2DLattice::setProbabilityDistribution(
+        const IFTDistribution2D& pdf)
+{
+    if (mp_pdf) delete mp_pdf;
+    mp_pdf = pdf.clone();
+}
+
+double InterferenceFunction2DLattice::evaluate(const cvector_t& q) const
+{
+    double result = 0.0;
+    double qxr = q.x().real();
+    double qyr = q.y().real();
+    double qx_frac, qy_frac;
+    calculateReciprocalVectorFraction(qxr, qyr, qx_frac, qy_frac);
+    int na = 40;
+    int nb = 40;
+    for (int i=-na-1; i<na+2; ++i)
+    {
+        for (int j=-nb-1; j<nb+2; ++j)
+        {
+            double qx = qx_frac + i*m_asx + j*m_bsx;
+            double qy = qy_frac + i*m_asy + j*m_bsy;
+            result += interferenceAtOneRecLatticePoint(qx, qy);
+        }
+    }
+//    double prefactor = 2.0*M_PI*m_lattice_params.m_corr_length_1*m_lattice_params.m_corr_length_2;
+    return result;
+}
+
+double InterferenceFunction2DLattice::interferenceAtOneRecLatticePoint(
+        double qx, double qy) const
+{
+    double qp1, qp2;
+    double gamma = m_lattice_params.m_xi + mp_pdf->getGamma();
+    double delta = mp_pdf->getDelta();
+    transformToPrincipalAxes(qx, qy, gamma, delta, qp1, qp2);
+    return mp_pdf->evaluate(qp1, qp2);
+}
+
+void InterferenceFunction2DLattice::transformToPrincipalAxes(double qx,
+        double qy, double gamma, double delta, double& q_pa_1, double& q_pa_2) const
+{
+    q_pa_1 = qx*std::cos(gamma) + qy*std::sin(gamma);
+    q_pa_2 = qx*std::cos(gamma+delta) + qy*std::sin(gamma+delta);
+}
+
+void InterferenceFunction2DLattice::calculateReciprocalVectorFraction(double qx,
+        double qy, double& qx_frac, double& qy_frac) const
+{
+    double a = m_lattice_params.m_length_1;
+    double b = m_lattice_params.m_length_2;
+    double xi = m_lattice_params.m_xi;
+    double xialpha = xi + m_lattice_params.m_angle;
+    int qa_int = (int)( a*(qx*std::cos(xi)+qy*std::sin(xi))/(2.0*M_PI) );
+    int qb_int = (int)( b*(qx*std::cos(xialpha)+qy*std::sin(xialpha))/(2.0*M_PI) );
+    qx_frac = qx - qa_int*m_asx - qb_int*m_bsx;
+    qy_frac = qy - qa_int*m_asy - qb_int*m_bsy;
+}
+
+void InterferenceFunction2DLattice::init_parameters()
+{
+    getParameterPool()->clear();
+}
+
+void InterferenceFunction2DLattice::initialize_rec_vectors()
+{
+    double sinalpha = std::sin(m_lattice_params.m_angle);
+    double ainv = 2.0*M_PI/m_lattice_params.m_length_1/sinalpha;
+    double binv = 2.0*M_PI/m_lattice_params.m_length_2/sinalpha;
+    double xi = m_lattice_params.m_xi;
+    double xialpha = xi + m_lattice_params.m_angle;
+    m_asx = ainv*std::sin(xialpha);
+    m_asy = -ainv*std::cos(xialpha);
+    m_bsx = -binv*std::sin(xi);
+    m_bsy = binv*std::cos(xi);
+}

Core/Samples/src/ParticleDecoration.cpp

@@ -24,10 +24,6 @@ ParticleDecoration::~ParticleDecoration()
     for (size_t i=0; i<m_particles.size(); ++i) {
         delete m_particles[i];
     }
-
-    for (size_t i=0; i<m_interference_functions.size(); ++i) {
-        delete m_interference_functions[i];
-    }
 }
 
 
@@ -139,7 +135,7 @@ IInterferenceFunctionStrategy* ParticleDecoration::createStrategy(
     size_t n_particles = m_particles.size();
     size_t n_ifs = m_interference_functions.size();
     if (n_ifs==1) {
-        InterferenceFunction1DParaCrystal *p_iff = dynamic_cast<InterferenceFunction1DParaCrystal *>(
+        const InterferenceFunction1DParaCrystal *p_iff = dynamic_cast<const InterferenceFunction1DParaCrystal *>(
                 m_interference_functions[0]);
         if (p_iff == 0 || p_iff->getKappa() == 0.0) {
             p_strategy = new DecouplingApproximationStrategy();
@@ -157,7 +153,7 @@ IInterferenceFunctionStrategy* ParticleDecoration::createStrategy(
         ostr << "n_particles:" << n_particles << " n_interference_function:" << n_ifs << ".";
         throw ClassInitializationException(ostr.str());
     }
-    p_strategy->init(form_factors, fractions, m_interference_functions);
+    p_strategy->init(form_factors, fractions, m_interference_functions.getSTLVector());
     return p_strategy;
 }