Merge branch 'master' into develop

App/App.pro

@@ -33,6 +33,7 @@ SOURCES += \
     src/TestIsGISAXS1.cpp \
     src/TestIsGISAXS2.cpp \
     src/TestIsGISAXS3.cpp \
+    src/TestIsGISAXS4.cpp \
     src/TestIsGISAXS9.cpp \
     src/TestIsGISAXS10.cpp \
     src/TestIsGISAXS11.cpp \
@@ -67,6 +68,7 @@ HEADERS += \
     inc/TestIsGISAXS1.h \
     inc/TestIsGISAXS2.h \
     inc/TestIsGISAXS3.h \
+    inc/TestIsGISAXS4.h \
     inc/TestIsGISAXS9.h \
     inc/TestIsGISAXS10.h \
     inc/TestIsGISAXS11.h \

App/inc/StandardSamples.h

@@ -19,6 +19,8 @@ ISample *IsGISAXS2_CylindersMixture();
 ISample *IsGISAXS3_CylinderDWBA();
 ISample *IsGISAXS3_CylinderBA();
 ISample *IsGISAXS3_CylinderBASize();
+ISample *IsGISAXS4_1DDL();
+ISample *IsGISAXS4_2DDL();
 ISample *IsGISAXS9_Pyramid();
 ISample *IsGISAXS9_RotatedPyramid();
 ISample *IsGISAXS10_CylindersParacrystal1D();

App/inc/TestIsGISAXS3.h

@@ -33,7 +33,7 @@ public:
     virtual void finalise();
 
 private:
-    // structure to holed info over several compare cases
+    // structure to hold info over several compare cases
     struct CompareStruct
     {
         CompareStruct(std::string _isginame, std::string _thisname, std::string _descr) : isginame(_isginame), thisname(_thisname), descr(_descr){}

App/inc/TestIsGISAXS4.h

+#ifndef TESTISGISAXS4_H_
+#define TESTISGISAXS4_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   TestIsGISAXS4.h
+//! @brief  Definition of TestIsGISAXS4 class for IsGISAXS validation
+//! @author Scientific Computing Group at FRM II
+//! @date   Oct 17, 2012
+
+#include "IFunctionalTest.h"
+
+//- -------------------------------------------------------------------
+//! @class TestIsGISAXS4
+//! @brief Comparison with IsGISAXS ex-4: cylinder with interference 1DDL or 2DDL
+//- -------------------------------------------------------------------
+class TestIsGISAXS4 : public IFunctionalTest
+{
+public:
+    TestIsGISAXS4();
+    virtual ~TestIsGISAXS4(){}
+
+    virtual void execute();
+    virtual void finalise();
+private:
+    // structure to hold info over several compare cases
+    struct CompareStruct
+    {
+        CompareStruct(std::string _isginame, std::string _thisname, std::string _descr) : isginame(_isginame), thisname(_thisname), descr(_descr){}
+        std::string isginame;
+        std::string thisname;
+        std::string descr;
+    };
+
+    std::string m_data_path;
+};
+
+
+#endif /* TESTISGISAXS4_H_ */

App/src/FunctionalTestFactory.cpp

@@ -6,6 +6,7 @@
 #include "TestIsGISAXS1.h"
 #include "TestIsGISAXS2.h"
 #include "TestIsGISAXS3.h"
+#include "TestIsGISAXS4.h"
 #include "TestIsGISAXS9.h"
 #include "TestIsGISAXS10.h"
 #include "TestIsGISAXS11.h"
@@ -42,6 +43,8 @@ FunctionalTestFactory::FunctionalTestFactory() : m_benchmark(0)
                  "functional test: isgisaxs ex-2 (mean form factors for particles with shape size distribution)");
     registerItem("isgisaxs3",   IFactoryCreateFunction<TestIsGISAXS3, IFunctionalTest>,
                  "functional test: isgisaxs ex-3 (cylinder FF)");
+    registerItem("isgisaxs4",   IFactoryCreateFunction<TestIsGISAXS4, IFunctionalTest>,
+                 "functional test: isgisaxs ex-4 (paracrystal structure factors)");
     registerItem("isgisaxs9",   IFactoryCreateFunction<TestIsGISAXS9, IFunctionalTest>,
                  "functional test: isgisaxs ex-9 (rotated pyramid FF)");
     registerItem("isgisaxs10",  IFactoryCreateFunction<TestIsGISAXS10, IFunctionalTest>,

App/src/IsGISAXSTools.cpp

@@ -201,7 +201,6 @@ void IsGISAXSTools::drawOutputDataDifference1D(const OutputData<double> &left, c
 
     *left_clone -= *right_clone;
     *left_clone /= *right_clone;
-    left_clone->scaleAll(100.0);
 
     std::string histo_name = histogram_title;
     if (histo_name.empty()) {
@@ -209,7 +208,7 @@ void IsGISAXSTools::drawOutputDataDifference1D(const OutputData<double> &left, c
     }
 
     TH1D h1_spect("difference", histo_name.c_str(), 40, -20.0, 20.0);
-    h1_spect.GetXaxis()->SetTitle("log10( 100*(we-isgi)/isgi )");
+    h1_spect.GetXaxis()->SetTitle("log10( (we-isgi)/isgi )");
 
     left_clone->resetIndex();
     while (left_clone->hasNext())

App/src/SampleFactory.cpp

@@ -37,6 +37,10 @@ SampleFactory::SampleFactory()
     registerItem("IsGISAXS3_CylinderBA", StandardSamples::IsGISAXS3_CylinderBA);
     registerItem("IsGISAXS3_CylinderBASize", StandardSamples::IsGISAXS3_CylinderBASize);
 
+    // IsGISAXS4 example: cylinders on top of substrate with paracrystal structure factors
+    registerItem("IsGISAXS4_1DDL", StandardSamples::IsGISAXS4_1DDL);
+    registerItem("IsGISAXS4_2DDL", StandardSamples::IsGISAXS4_2DDL);
+
     // IsGISAXS9 example: pyramid and rotated pyramid
     registerItem("IsGISAXS9_Pyramid", StandardSamples::IsGISAXS9_Pyramid);
     registerItem("IsGISAXS9_RotatedPyramid", StandardSamples::IsGISAXS9_RotatedPyramid);

App/src/StandardSamples.cpp

@@ -13,6 +13,7 @@
 #include "Crystal.h"
 #include "MesoCrystal.h"
 #include "InterferenceFunction1DParaCrystal.h"
+#include "InterferenceFunction2DParaCrystal.h"
 #include "FormFactorWeighted.h"
 #include "StochasticGaussian.h"
 #include "Numeric.h"
@@ -389,6 +390,54 @@ ISample *StandardSamples::IsGISAXS3_CylinderBASize()
     return p_multi_layer;
 }
 
+/* ************************************************************************* */
+// IsGISAXS4 functional test: cylinders with 1DDL structure factor
+/* ************************************************************************* */
+ISample *StandardSamples::IsGISAXS4_1DDL()
+{
+    MultiLayer *p_multi_layer = new MultiLayer();
+    complex_t n_air(1.0, 0.0);
+    complex_t n_substrate(1.0-6e-6, 2e-8);
+    complex_t n_particle(1.0-6e-4, 2e-8);
+    const IMaterial *p_air_material = MaterialManager::instance().addHomogeneousMaterial("Air", n_air);
+    const IMaterial *p_substrate_material = MaterialManager::instance().addHomogeneousMaterial("Substrate", n_substrate);
+    Layer air_layer;
+    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)));
+    particle_decoration.addInterferenceFunction(p_interference_function);
+    LayerDecorator air_layer_decorator(air_layer, particle_decoration);
+
+    p_multi_layer->addLayer(air_layer_decorator);
+    p_multi_layer->addLayer(substrate_layer);
+    return p_multi_layer;
+}
+
+// IsGISAXS4 functional test: cylinders with 2DDL structure factor
+ISample *StandardSamples::IsGISAXS4_2DDL()
+{
+    MultiLayer *p_multi_layer = new MultiLayer();
+    complex_t n_air(1.0, 0.0);
+    complex_t n_substrate(1.0-6e-6, 2e-8);
+    complex_t n_particle(1.0-6e-4, 2e-8);
+    const IMaterial *p_air_material = MaterialManager::instance().addHomogeneousMaterial("Air", n_air);
+    const IMaterial *p_substrate_material = MaterialManager::instance().addHomogeneousMaterial("Substrate", n_substrate);
+    Layer air_layer;
+    air_layer.setMaterial(p_air_material);
+    Layer substrate_layer;
+    substrate_layer.setMaterial(p_substrate_material);
+    IInterferenceFunction *p_interference_function = InterferenceFunction2DParaCrystal::createHexagonal(20.0*Units::nanometer,1.0*Units::nanometer, 0.0,
+            20.0*Units::micrometer, 20.0*Units::micrometer);
+    ParticleDecoration particle_decoration( new Particle(n_particle, new FormFactorCylinder(5*Units::nanometer, 5*Units::nanometer)));
+    particle_decoration.addInterferenceFunction(p_interference_function);
+    LayerDecorator air_layer_decorator(air_layer, particle_decoration);
+
+    p_multi_layer->addLayer(air_layer_decorator);
+    p_multi_layer->addLayer(substrate_layer);
+    return p_multi_layer;
+}
 
 /* ************************************************************************* */
 // IsGISAXS9 functional test: pyramid
@@ -466,9 +515,9 @@ ISample *StandardSamples::IsGISAXS10_CylindersParacrystal1D()
     air_layer.setMaterial(p_air_material);
     Layer substrate_layer;
     substrate_layer.setMaterial(p_substrate_material);
-    IInterferenceFunction *p_interference_funtion = new InterferenceFunction1DParaCrystal(20.0*Units::nanometer,7*Units::nanometer, 1e7*Units::nanometer);
+    IInterferenceFunction *p_interference_function = new InterferenceFunction1DParaCrystal(20.0*Units::nanometer,7*Units::nanometer, 1e7*Units::nanometer);
     ParticleDecoration particle_decoration(new Particle(n_particle, new FormFactorCylinder(5*Units::nanometer, 5*Units::nanometer)));
-    particle_decoration.addInterferenceFunction(p_interference_funtion);
+    particle_decoration.addInterferenceFunction(p_interference_function);
     //    particle_decoration.setTotalParticleSurfaceDensity(1.0/(20.0*Units::nanometer*20.0*Units::nanometer));
     LayerDecorator air_layer_decorator(air_layer, particle_decoration);
 
@@ -541,13 +590,13 @@ ISample *StandardSamples::MesoCrystal1()
     air_layer.setMaterial(p_air_material);
     Layer substrate_layer;
     substrate_layer.setMaterial(p_substrate_material);
-//    IInterferenceFunction *p_interference_funtion = new InterferenceFunctionNone();
-    IInterferenceFunction *p_interference_funtion = new InterferenceFunction1DParaCrystal(800.0*Units::nanometer,
+//    IInterferenceFunction *p_interference_function = new InterferenceFunctionNone();
+    IInterferenceFunction *p_interference_function = new InterferenceFunction1DParaCrystal(800.0*Units::nanometer,
         50*Units::nanometer, 1e7*Units::nanometer);
     ParticleDecoration particle_decoration;
     particle_decoration.addParticle(meso.clone(), 0.0, 0.5);
     particle_decoration.addParticle(meso2.clone(), 0.0, 0.5);
-    particle_decoration.addInterferenceFunction(p_interference_funtion);
+    particle_decoration.addInterferenceFunction(p_interference_function);
     LayerDecorator air_layer_decorator(air_layer, particle_decoration);
 
     p_multi_layer->addLayer(air_layer_decorator);
@@ -587,7 +636,7 @@ ISample *StandardSamples::MesoCrystal2()
     avg_layer.setThickness(0.2*Units::micrometer);
     Layer substrate_layer;
     substrate_layer.setMaterial(p_substrate_material);
-    IInterferenceFunction *p_interference_funtion = new InterferenceFunctionNone();
+    IInterferenceFunction *p_interference_function = new InterferenceFunctionNone();
     ParticleDecoration particle_decoration;
 
     //
@@ -612,7 +661,7 @@ ISample *StandardSamples::MesoCrystal2()
     particle_decoration.addParticle(meso_crystal, 0.2*Units::micrometer);
 
     particle_decoration.setTotalParticleSurfaceDensity(surface_density);
-    particle_decoration.addInterferenceFunction(p_interference_funtion);
+    particle_decoration.addInterferenceFunction(p_interference_function);
     LayerDecorator avg_layer_decorator(avg_layer, particle_decoration);
 
     p_multi_layer->addLayer(air_layer);

App/src/TestIsGISAXS4.cpp

+#include "TestIsGISAXS4.h"
+#include "IsGISAXSTools.h"
+#include "Units.h"
+#include "Utils.h"
+#include "MultiLayer.h"
+#include "GISASExperiment.h"
+#include "SampleFactory.h"
+#include "DrawHelper.h"
+
+#include "TCanvas.h"
+#include <gsl/gsl_errno.h>
+
+TestIsGISAXS4::TestIsGISAXS4() : IFunctionalTest("TestIsGISAXS4")
+{
+    m_data_path = std::string(Utils::FileSystem::GetHomePath()+"./Examples/IsGISAXS_examples/ex-4/");
+}
+
+void TestIsGISAXS4::execute()
+{
+    gsl_set_error_handler_off();
+
+    GISASExperiment experiment(mp_options);
+    experiment.setDetectorParameters(100, 0.0*Units::degree, 2.0*Units::degree, 100, 0.0*Units::degree, 2.0*Units::degree, true);
+    experiment.setBeamParameters(1.0*Units::angstrom, -0.2*Units::degree, 0.0*Units::degree);
+
+    MultiLayer *p_sample(0);
+
+//    // 1DDL
+    p_sample = dynamic_cast<MultiLayer *>(SampleFactory::instance().createItem("IsGISAXS4_1DDL"));
+    experiment.setSample(*p_sample);
+    experiment.runSimulation();
+    IsGISAXSTools::writeOutputDataToFile(*experiment.getOutputData(), m_data_path+"this_1DDL.ima");
+    delete p_sample;
+
+    // 2DDL
+    p_sample = dynamic_cast<MultiLayer *>(SampleFactory::instance().createItem("IsGISAXS4_2DDL"));
+    experiment.setSample(*p_sample);
+    experiment.runSimulation();
+    IsGISAXSTools::writeOutputDataToFile(*experiment.getOutputData(), m_data_path+"this_2DDLh.ima");
+    delete p_sample;
+}
+
+
+void TestIsGISAXS4::finalise()
+{
+    std::vector< CompareStruct > tocompare;
+    tocompare.push_back( CompareStruct("isgi_1DDL.ima",      "this_1DDL.ima",      "Cylinder 1DDL") );
+    tocompare.push_back( CompareStruct("isgi_2DDLh.ima",      "this_2DDLh.ima",      "Cylinder 2DDL") );
+
+    for(size_t i=0; i<tocompare.size(); ++i) {
+        OutputData<double> *isgi_data = IsGISAXSTools::readOutputDataFromFile( m_data_path+tocompare[i].isginame );
+        OutputData<double> *our_data = IsGISAXSTools::readOutputDataFromFile( m_data_path+tocompare[i].thisname );
+
+        std::ostringstream os;
+        os<<i;
+        std::string cname = getName()+"_c"+os.str();
+        TCanvas *c1 = DrawHelper::instance().createAndRegisterCanvas(cname.c_str(), tocompare[i].descr);
+        c1->Divide(2,2);
+
+        IsGISAXSTools::setMinimum(1.);
+        // our calculations
+        c1->cd(1); gPad->SetLogz();
+        IsGISAXSTools::drawOutputDataInPad(*our_data, "CONT4 Z", "Our cylinder FF");
+
+        // isgisaxs data
+        c1->cd(2); gPad->SetLogz();
+        IsGISAXSTools::drawOutputDataInPad(*isgi_data, "CONT4 Z", "IsGisaxs mean FF");
+
+        // difference
+        c1->cd(3);
+        IsGISAXSTools::setMinimum(-0.0001);
+        IsGISAXSTools::setMaximum(0.0001);
+        IsGISAXSTools::drawOutputDataDifference2D(*our_data, *isgi_data, "CONT4 Z", "2D Difference map");
+
+        // difference
+        c1->cd(4);
+        IsGISAXSTools::resetMinimumAndMaximum();
+        //IsGISAXSTools::setMinimum(1);
+        IsGISAXSTools::drawOutputDataDifference1D(*our_data, *isgi_data, "", "Difference spectra");
+
+        IsGISAXSTools::resetMinimum(); IsGISAXSTools::resetMaximum();
+        delete isgi_data;
+        delete our_data;
+    }
+}

App/src/TestMesoCrystal1.cpp

@@ -45,7 +45,7 @@ void TestMesoCrystal1::execute()
            , 256, -0.4*Units::degree, 0.066);
 //    experiment.setDetectorParameters(2, 0.96*Units::degree, 0.962*Units::degree
 //           , 2 , 0.376*Units::degree, 0.378*Units::degree);
-    experiment.setDetectorResolutionFunction(new ResolutionFunction2DSimple(0.0002, 0.0002));
+    experiment.setDetectorResolutionFunction(new ResolutionFunction2DSimple(0.0004, 0.0004));
     experiment.setBeamParameters(1.77*Units::angstrom, -0.4*Units::degree, 0.0*Units::degree);
     experiment.setBeamIntensity(8e12);
 
@@ -79,16 +79,16 @@ MesoCrystalBuilder::MesoCrystalBuilder()
 //, m_sigma_nanoparticle_radius(0.14*Units::nanometer)
 //, m_sigma_lattice_length_a(1.5*Units::nanometer)
 //, m_roughness(1.0*Units::nanometer)
-: m_meso_radius(570.103*Units::nanometer)
-, m_surface_filling_ratio(0.113898)
-, m_meso_height(423.33*Units::nanometer)
-, m_sigma_meso_height(10.0608*Units::nanometer)
-, m_sigma_meso_radius(10.0882*Units::nanometer)
-, m_lattice_length_a(6.21014*Units::nanometer)
-, m_nanoparticle_radius(5.99176*Units::nanometer)
-, m_sigma_nanoparticle_radius(0.0681535*Units::nanometer)
-, m_sigma_lattice_length_a(2.61389*Units::nanometer)
-, m_roughness(28.0626*Units::nanometer)
+: m_meso_radius(1370*Units::nanometer)
+, m_surface_filling_ratio(0.114748)
+, m_meso_height(265.276*Units::nanometer)
+, m_sigma_meso_height(10.8148*Units::nanometer)
+, m_sigma_meso_radius(14.0738*Units::nanometer)
+, m_lattice_length_a(6.22525*Units::nanometer)
+, m_nanoparticle_radius(5.05257*Units::nanometer)
+, m_sigma_nanoparticle_radius(0.0877905*Units::nanometer)
+, m_sigma_lattice_length_a(1.95504*Units::nanometer)
+, m_roughness(0.13464*Units::nanometer)
 {
     init_parameters();
 }

App/src/TestMesoCrystal2.cpp

@@ -68,46 +68,46 @@ void TestMesoCrystal2::execute()
     OutputDataReader *reader = OutputDataIOFactory::instance().getReader(file_name);
     OutputData<double > *real_data = reader->getOutputData();
     OutputData<double > *real_data_half = doubleBinSize(*real_data);
-//    OutputData<double > *real_data_quarter = doubleBinSize(*real_data_half);
+    OutputData<double > *real_data_quarter = doubleBinSize(*real_data_half);
 //    OutputData<double > *real_data_eighth = doubleBinSize(*real_data_quarter);
     delete reader;
     c1->cd(1); gPad->SetLogz();
-    IsGISAXSTools::drawOutputDataInPad(*real_data_half, "CONT4 Z", "experiment");
+    IsGISAXSTools::drawOutputDataInPad(*real_data_quarter, "CONT4 Z", "experiment");
     c1->Update();
 
     // initializing experiment using real data to have detector axises like in real_data
-    initializeExperiment(real_data_half);
+    initializeExperiment(real_data_quarter);
     mp_experiment->printParameters();
 
     // setting fitSuite
     FitSuite *fitSuite = new FitSuite();
     fitSuite->setExperiment(mp_experiment);
-    fitSuite->setRealData(*real_data_half);
+    fitSuite->setRealData(*real_data_quarter);
     fitSuite->setMinimizer( new ROOTMinimizer("Minuit2", "Migrad") );
-    fitSuite->addFitParameter("*/lattice_length_a", 6.2*Units::nanometer, 1.0*Units::nanometer,
-            TRange<double>(4.0*Units::nanometer, 8.0*Units::nanometer) );
+//    fitSuite->addFitParameter("*/lattice_length_a", 6.2*Units::nanometer, 1.0*Units::nanometer,
+//            TRange<double>(4.0*Units::nanometer, 8.0*Units::nanometer) );
     fitSuite->addFitParameter("*/nanoparticle_radius", 5.7*Units::nanometer, 1.0*Units::nanometer,
             TRange<double>(2.0*Units::nanometer, 8.0*Units::nanometer) );
     fitSuite->addFitParameter("*/sigma_nanoparticle_radius", 0.1*Units::nanometer, 0.05*Units::nanometer,
             TRange<double>(0.01*Units::nanometer, 2.0*Units::nanometer) );
-    fitSuite->addFitParameter("*/meso_height", 500.0*Units::nanometer, 100.0*Units::nanometer,
-            TRange<double>(100.0*Units::nanometer, 2000.0*Units::nanometer) );
-    fitSuite->addFitParameter("*/meso_radius", 1000.0*Units::nanometer, 100.0*Units::nanometer,
-            TRange<double>(100.0*Units::nanometer, 5000.0*Units::nanometer) );
-    fitSuite->addFitParameter("*/sigma_meso_height", 5.0*Units::nanometer, 1.0*Units::nanometer,
-            TRange<double>(10.0*Units::nanometer, 200.0*Units::nanometer) );
-    fitSuite->addFitParameter("*/sigma_meso_radius", 50.0*Units::nanometer, 10.0*Units::nanometer,
-            TRange<double>(10.0*Units::nanometer, 500.0*Units::nanometer) );
-    fitSuite->addFitParameter("*/sigma_lattice_length_a", 1.5*Units::nanometer, 0.5*Units::nanometer,
-            TRange<double>(0.01*Units::nanometer, 4.0*Units::nanometer) );
-    fitSuite->addFitParameter("*/surface_filling_ratio", 0.25, 0.1,
-            TRange<double>(0.1, 0.4) );
-    fitSuite->addFitParameter("*/roughness", 1.0*Units::nanometer, 0.1*Units::nanometer,
-            TRange<double>(0.01*Units::nanometer, 50.0*Units::nanometer) );
-    fitSuite->addFitParameter("*/ResolutionFunction2D/sigma_x", 0.0002, 0.00001,
-            TRange<double>(0.0, 0.002) );
-    fitSuite->addFitParameter("*/ResolutionFunction2D/sigma_y", 0.0002, 0.00001,
-            TRange<double>(0.0, 0.002) );
+//    fitSuite->addFitParameter("*/meso_height", 500.0*Units::nanometer, 100.0*Units::nanometer,
+//            TRange<double>(100.0*Units::nanometer, 2000.0*Units::nanometer) );
+//    fitSuite->addFitParameter("*/meso_radius", 1000.0*Units::nanometer, 100.0*Units::nanometer,
+//            TRange<double>(100.0*Units::nanometer, 5000.0*Units::nanometer) );
+//    fitSuite->addFitParameter("*/sigma_meso_height", 5.0*Units::nanometer, 1.0*Units::nanometer,
+//            TRange<double>(10.0*Units::nanometer, 200.0*Units::nanometer) );
+//    fitSuite->addFitParameter("*/sigma_meso_radius", 50.0*Units::nanometer, 10.0*Units::nanometer,
+//            TRange<double>(10.0*Units::nanometer, 500.0*Units::nanometer) );
+//    fitSuite->addFitParameter("*/sigma_lattice_length_a", 1.5*Units::nanometer, 0.5*Units::nanometer,
+//            TRange<double>(0.01*Units::nanometer, 4.0*Units::nanometer) );
+//    fitSuite->addFitParameter("*/surface_filling_ratio", 0.25, 0.1,
+//            TRange<double>(0.1, 0.4) );
+//    fitSuite->addFitParameter("*/roughness", 1.0*Units::nanometer, 0.1*Units::nanometer,
+//            TRange<double>(0.01*Units::nanometer, 50.0*Units::nanometer) );
+//    fitSuite->addFitParameter("*/ResolutionFunction2D/sigma_x", 0.0002, 0.00001,
+//            TRange<double>(0.0, 0.002) );
+//    fitSuite->addFitParameter("*/ResolutionFunction2D/sigma_y", 0.0002, 0.00001,
+//            TRange<double>(0.0, 0.002) );
 
     IsGISAXSTools::setMinimum(1e2);
     std::string tree_file_name = Utils::FileSystem::GetHomePath()+"Examples/MesoCrystals/ex02_fitspheres/mesofit.tree";

Core/Core.pro

@@ -63,6 +63,7 @@ SOURCES += \
     Samples/src/ICompositeSample.cpp \
     Samples/src/IMaterial.cpp \
     Samples/src/InterferenceFunction1DParaCrystal.cpp \
+    Samples/src/InterferenceFunction2DParaCrystal.cpp \
     Samples/src/IRoughness.cpp \
     Samples/src/ISample.cpp \
     Samples/src/Lattice.cpp \
@@ -180,6 +181,7 @@ HEADERS += \
     Samples/inc/IInterferenceFunction.h \
     Samples/inc/IMaterial.h \
     Samples/inc/InterferenceFunction1DParaCrystal.h \
+    Samples/inc/InterferenceFunction2DParaCrystal.h \
     Samples/inc/InterferenceFunctionNone.h \
     Samples/inc/IRoughness.h \
     Samples/inc/ISample.h \

Core/Samples/inc/InterferenceFunction1DParaCrystal.h

@@ -11,7 +11,7 @@
 // ********************************************************************
 //! @file   InterferenceFunction1DParaCrystal.h
 //! @brief  Definition of InterferenceFunction1DParaCrystal class
-//! @author herck
+//! @author Scientific Computing Group at FRM II
 //! @date   19.06.2012
 
 #include "IInterferenceFunction.h"
@@ -39,7 +39,7 @@ private:
     //! initialize pool parameters, i.e. register some of class members for later access via parameter pool
     virtual void init_parameters();
 
-	// replace these with strategy pattern for different algorithms
+	//TODO: replace these with strategy pattern for different algorithms
 	complex_t FTGaussianCorrLength(double qpar) const;
 };
 

Core/Samples/inc/InterferenceFunction2DParaCrystal.h

+#ifndef INTERFERENCEFUNCTION2DPARACRYSTAL_H_
+#define INTERFERENCEFUNCTION2DPARACRYSTAL_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   InterferenceFunction2DParaCrystal.h
+//! @brief  Definition of InterferenceFunction2DParaCrystal class
+//! @author Scientific Computing Group at FRM II
+//! @date   Oct 17, 2012
+
+#include "IInterferenceFunction.h"
+
+class InterferenceFunction2DParaCrystal : public IInterferenceFunction
+{
+public:
+    InterferenceFunction2DParaCrystal(double peak_distance, double alpha_lattice, double width, double corr_length=0.0);
+    virtual ~InterferenceFunction2DParaCrystal() {}
+    virtual InterferenceFunction2DParaCrystal *clone() const {
+        InterferenceFunction2DParaCrystal *p_new = new InterferenceFunction2DParaCrystal(m_peak_distance, m_alpha_lattice, m_width, m_corr_length);
+        p_new->setDomainSizes(m_domain_size_1, m_domain_size_2);
+        return p_new;
+    }
+    static InterferenceFunction2DParaCrystal *createSquare(double peak_distance, double width, double corr_length=0.0,
+            double domain_size_1=0.0, double domain_size_2=0.0) {
+        InterferenceFunction2DParaCrystal *p_new = new InterferenceFunction2DParaCrystal(peak_distance, M_PI/2.0, width, corr_length);
+        p_new->setDomainSizes(domain_size_1, domain_size_2);
+        return p_new;        return new InterferenceFunction2DParaCrystal(peak_distance, M_PI/2.0, width, corr_length);
+    }
+    static InterferenceFunction2DParaCrystal *createHexagonal(double peak_distance, double width, double corr_length=0.0,
+            double domain_size_1=0.0, double domain_size_2=0.0) {
+        InterferenceFunction2DParaCrystal *p_new = new InterferenceFunction2DParaCrystal(peak_distance, 2.0*M_PI/3.0, width, corr_length);
+        p_new->setDomainSizes(domain_size_1, domain_size_2);
+        return p_new;
+    }
+
+    void setDomainSizes(double size_1, double size_2) {
+        m_domain_size_1 = size_1;
+        m_domain_size_2 = size_2;
+    }
+
+    virtual double evaluate(const cvector_t &q) const;
+
+protected:
+    double m_peak_distance;
+    double m_alpha_lattice; //!< Angle between lattice basis vectors
+    double m_width;
+    double m_corr_length;
+    bool m_use_corr_length;
+    double m_domain_size_1;
+    double m_domain_size_2;
+private:
+    //! copy constructor and assignment operator are hidden since there is a clone method
+    InterferenceFunction2DParaCrystal(const InterferenceFunction2DParaCrystal &);
+    InterferenceFunction2DParaCrystal &operator=(const InterferenceFunction2DParaCrystal &);
+
+    //! initialize pool parameters, i.e. register some of class members for later access via parameter pool
+    virtual void init_parameters();
+
+    //! Calculate interference function for fixed rotation xi
+    double interferenceForXi(double xi) const;
+
+    //! calculate interference function for fixed xi in 1d
+    double interference1D(double qpar, double xi, double domain_size) const;
+
+    //TODO: replace these with strategy pattern for different algorithms
+    complex_t FTCauchyDistribution(double qpar, double xi) const;
+
+    mutable double m_qpar;
+};
+
+double wrapFunctionForGSL(double xi, void* p_unary_function);
+
+
+
+#endif /* INTERFERENCEFUNCTION2DPARACRYSTAL_H_ */

Core/Samples/src/InterferenceFunction1DParaCrystal.cpp

@@ -25,7 +25,9 @@ void InterferenceFunction1DParaCrystal::init_parameters()
 
 double InterferenceFunction1DParaCrystal::evaluate(const cvector_t &q) const
 {
-	double qpar = q.magxy().real();
+    double qxr = q.x().real();
+    double qyr = q.y().real();
+    double qpar = std::sqrt(qxr*qxr + qyr*qyr);
 	complex_t p_transformed = FTGaussianCorrLength(qpar);
 	double interference_function = 1.0 + 2*(p_transformed/(complex_t(1.0, 0.0)-p_transformed)).real();
 	return interference_function;

Core/Samples/src/InterferenceFunction2DParaCrystal.cpp

+#include "InterferenceFunction2DParaCrystal.h"
+#include "MathFunctions.h"
+
+#include <functional>
+
+double wrapFunctionForGSL(double xi, void* p_unary_function)
+{
+    std::binder1st<std::const_mem_fun1_t<double, InterferenceFunction2DParaCrystal, double> > *p_f =
+            (std::binder1st<std::const_mem_fun1_t<double, InterferenceFunction2DParaCrystal, double> > *)p_unary_function;
+    return (*p_f)(xi);
+}
+
+InterferenceFunction2DParaCrystal::InterferenceFunction2DParaCrystal(double peak_distance, double alpha_lattice, double width, double corr_length)
+: m_peak_distance(peak_distance)
+, m_alpha_lattice(alpha_lattice)
+, m_width(width)
+, m_corr_length(corr_length)
+, m_use_corr_length(true)
+, m_domain_size_1(0.0)
+, m_domain_size_2(0.0)
+{
+    setName("InterferenceFunction2DParaCrystal");
+    if (m_corr_length==0.0) {
+        m_use_corr_length = false;
+    }
+    init_parameters();
+}
+
+double InterferenceFunction2DParaCrystal::evaluate(const cvector_t &q) const
+{
+    double qxr = q.x().real();
+    double qyr = q.y().real();
+    m_qpar = std::sqrt(qxr*qxr + qyr*qyr);
+    std::binder1st<std::const_mem_fun1_t<double, InterferenceFunction2DParaCrystal, double> > f_base = std::bind1st(std::mem_fun(&InterferenceFunction2DParaCrystal::interferenceForXi), this);
+    gsl_function f;
+    f.function = &wrapFunctionForGSL;
+    f.params = &f_base;
+    double result = MathFunctions::Integrate1D(&f, 0.0, M_PI)/M_PI;
+    return result;
+}
+
+void InterferenceFunction2DParaCrystal::init_parameters()
+{
+    getParameterPool()->clear();
+    getParameterPool()->registerParameter("peak_distance", &m_peak_distance);
+    getParameterPool()->registerParameter("width", &m_width);
+    getParameterPool()->registerParameter("corr_length", &m_corr_length);
+}
+
+
+double InterferenceFunction2DParaCrystal::interferenceForXi(double xi) const
+{
+    double result = interference1D(m_qpar, xi, m_domain_size_1)*interference1D(m_qpar, xi + m_alpha_lattice, m_domain_size_2);
+    return result;
+}
+
+double InterferenceFunction2DParaCrystal::interference1D(double qpar,
+        double xi, double domain_size) const
+{
+    double result;
+    int n = (int)std::abs(domain_size/m_peak_distance);
+    double nd = (double)n;
+    complex_t fp = FTCauchyDistribution(qpar, xi);
+    if (n<1) {
+        result = ((1.0 + fp)/(1.0 - fp)).real();
+    }
+    else {
+        if (std::abs(1.0-fp) < Numeric::double_epsilon) {
+            result = nd;
+        }
+        else {
+            complex_t tmp;
+            if (std::log(std::abs(fp)+Numeric::double_min)*nd < std::log(Numeric::double_min)) {
+                tmp = complex_t(0.0, 0.0);
+            }
+            else {
+                tmp = std::pow(fp,n-1);
+            }
+            // TODO: remove factor 2.0 in second term
+            double intermediate = ((1.0-1.0/nd)*fp/(1.0-fp) - fp*fp*(1.0-tmp)/nd/(1.0-fp)/(1.0-fp)).real();
+            result = 1.0 + 2.0*intermediate;
+        }
+    }
+    return result;
+}
+
+complex_t InterferenceFunction2DParaCrystal::FTCauchyDistribution(double qpar, double xi) const
+{
+    complex_t phase = std::exp(complex_t(0.0, 1.0)*qpar*m_peak_distance*std::cos(xi));
+    double amplitude = std::pow(1.0+qpar*qpar*m_width*m_width, -1.5);
+    complex_t result = phase*amplitude;
+    if (m_use_corr_length) {
+        result *= std::exp(-m_peak_distance/m_corr_length);
+    }
+    return result;
+}
+

Core/Tools/inc/MathFunctions.h

@@ -14,14 +14,16 @@
 //! @author Scientific Computing Group at FRM II
 //! @date   20.04.2012
 
-#include <cstdlib>
-#include <vector>
-#include "gsl/gsl_sf_bessel.h"
-#include "gsl/gsl_sf_trig.h"
 #include "Types.h"
 #include "Units.h"
 #include "Numeric.h"
 
+#include <cstdlib>
+#include <vector>
+
+#include "gsl/gsl_sf_bessel.h"
+#include "gsl/gsl_sf_trig.h"
+#include "gsl/gsl_integration.h"
 
 namespace MathFunctions
 {
@@ -55,53 +57,22 @@ std::vector<complex_t > FastFourierTransform(const std::vector<double > &data, T
 std::vector<complex_t> ConvolveFFT(const std::vector<double> &signal, const std::vector<double> &resfunc);
 
 //! fast sine calculations (not actually fast)
-inline double FastSin(const double& x) {
-    const double P = 0.225f;
-    const double A = 16 * std::sqrt(P);
-    const double B = (1 - P) / std::sqrt(P);
-    double y = x / (2 * M_PI);
-    y = y - std::floor(y + 0.5);  // y in range -0.5..0.5
-    y = A * y * (0.5 - std::abs(y));
-    return y * (B + std::abs(y));
-}
+double FastSin(const double& x);
 
 //! fast cosine calculation  (not actually fast)
-inline double FastCos(const double& x) {
-    const double P = 0.225f;
-    const double A = 16 * std::sqrt(P);
-    const double B = (1 - P) / std::sqrt(P);
-    double y = x / (2 * M_PI)+0.25; // pi/2. shift
-    y = y - std::floor(y + 0.5);  // y in range -0.5..0.5
-    y = A * y * (0.5 - std::abs(y));
-    return y * (B + std::abs(y));
-}
+double FastCos(const double& x);
 
 //! fast complex sine calculation
-inline complex_t FastSin(const complex_t &x) {
-    // sin(a+bi) = sin(a)cosh(b) + i*cos(a)*sinh(b);
-    //return complex_t( FastSin(x.real())*std::cosh(x.imag()), FastCos(x.real())*std::sinh(x.imag()));
-    return complex_t( std::sin(x.real())*std::cosh(x.imag()), std::cos(x.real())*std::sinh(x.imag()));
-}
+complex_t FastSin(const complex_t &x);
 
 //! fast complex cosine calculation
-inline complex_t FastCos(const complex_t &x) {
-    // cos(a+bi) = cos(a)cosh(b) - i*sin(a)*sinh(b);
-    //return complex_t( FastCos(x.real())*std::cosh(x.imag()), -1*FastSin(x.real())*std::sinh(x.imag()));
-    return complex_t( std::cos(x.real())*std::cosh(x.imag()), -1*std::sin(x.real())*std::sinh(x.imag()));
-}
+complex_t FastCos(const complex_t &x);
 
 //! simultaneous complex sine and cosine calculations
-inline void FastSinCos(const complex_t &x, complex_t &xsin, complex_t &xcos)
-{
-    double sina = FastSin(x.real());
-    double cosa = std::sqrt(1-sina*sina);
-    double sinhb = std::sinh(x.imag());
-    double coshb = std::sqrt(1-sinhb*sinhb);
-    xsin.real() =sina*coshb;
-    xsin.imag() =cosa*sinhb;
-    xcos.real() =cosa*coshb;
-    xcos.imag() =-sina*sinhb;
-}
+void FastSinCos(const complex_t &x, complex_t &xsin, complex_t &xcos);
+
+//! numerical integration
+double Integrate1D(gsl_function *p_function, double start, double end);
 
 } // Namespace MathFunctions
 
@@ -145,7 +116,53 @@ inline complex_t MathFunctions::Laue(complex_t value, size_t N) // Exp(iNx/2)*Si
     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);
 }
 
+//! fast sine calculations (not actually fast)
+inline double MathFunctions::FastSin(const double& x) {
+    const double P = 0.225f;
+    const double A = 16 * std::sqrt(P);
+    const double B = (1 - P) / std::sqrt(P);
+    double y = x / (2 * M_PI);
+    y = y - std::floor(y + 0.5);  // y in range -0.5..0.5
+    y = A * y * (0.5 - std::abs(y));
+    return y * (B + std::abs(y));
+}
+
+//! fast cosine calculation  (not actually fast)
+inline double MathFunctions::FastCos(const double& x) {
+    const double P = 0.225f;
+    const double A = 16 * std::sqrt(P);
+    const double B = (1 - P) / std::sqrt(P);
+    double y = x / (2 * M_PI)+0.25; // pi/2. shift
+    y = y - std::floor(y + 0.5);  // y in range -0.5..0.5
+    y = A * y * (0.5 - std::abs(y));
+    return y * (B + std::abs(y));
+}
+
+//! fast complex sine calculation
+inline complex_t MathFunctions::FastSin(const complex_t &x) {
+    // sin(a+bi) = sin(a)cosh(b) + i*cos(a)*sinh(b);
+    //return complex_t( FastSin(x.real())*std::cosh(x.imag()), FastCos(x.real())*std::sinh(x.imag()));
+    return complex_t( std::sin(x.real())*std::cosh(x.imag()), std::cos(x.real())*std::sinh(x.imag()));
+}
 
+//! fast complex cosine calculation
+inline complex_t MathFunctions::FastCos(const complex_t &x) {
+    // cos(a+bi) = cos(a)cosh(b) - i*sin(a)*sinh(b);
+    //return complex_t( FastCos(x.real())*std::cosh(x.imag()), -1*FastSin(x.real())*std::sinh(x.imag()));
+    return complex_t( std::cos(x.real())*std::cosh(x.imag()), -1*std::sin(x.real())*std::sinh(x.imag()));
+}
 
+//! simultaneous complex sine and cosine calculations
+inline void MathFunctions::FastSinCos(const complex_t &x, complex_t &xsin, complex_t &xcos)
+{
+    double sina = FastSin(x.real());
+    double cosa = std::sqrt(1-sina*sina);
+    double sinhb = std::sinh(x.imag());
+    double coshb = std::sqrt(1-sinhb*sinhb);
+    xsin.real() =sina*coshb;
+    xsin.imag() =cosa*sinhb;
+    xcos.real() =cosa*coshb;
+    xcos.imag() =-sina*sinhb;
+}
 
 #endif // MATHFUNCTIONS_H

Core/Tools/src/MathFunctions.cpp

@@ -192,3 +192,13 @@ std::vector<complex_t> MathFunctions::ConvolveFFT(const std::vector<double> &sig
     std::vector<complex_t > result = MathFunctions::FastFourierTransform(fft_prod, MathFunctions::BackwardFFT);
     return result;
 }
+
+double MathFunctions::Integrate1D(gsl_function *p_function, double start,
+        double end)
+{
+    gsl_integration_workspace *ws = gsl_integration_workspace_alloc(200);
+    double result, error;
+    gsl_integration_qag(p_function, start, end, 1e-10, 1e-8, 50, 1, ws, &result, &error);
+    gsl_integration_workspace_free(ws);
+    return result;
+}

Core/test001.h

+for testing

Examples/IsGISAXS_examples/ex-4/1DDL.inp

+##########################################
+#   GISAXS SIMULATIONS : INPUT PARAMETERS   
+###########################################
+
+# Base filename
+1DDL
+############################  Framework and beam  parameters ############################################
+# Framework   Diffuse,  Multilayer, Number of index slices, Polarization 
+   DWBA         LMA              0             25              ss 
+# Beam Wavelenght  :  Lambda(nm),  Wl_distribution,  Sigma_Wl/Wl,  Wl_min(nm), Wl_max(nm), nWl,  xWl
+                      0.1            none               0.3           0.08         0.12     20    3
+# Beam Alpha_i     :  Alpha_i(deg), Ai_distribution, Sigma_Ai(deg), Ai_min(deg), Ai_max(deg), nAi, xAi
+                       0.2            none           0.1         0.15              0.25    30   2
+# Beam 2Theta_i     :  2Theta_i(deg), Ti_distribution, Sigma_Ti(deg), Ti_min(deg), Ti_max(deg), nTi, XTi
+                         0              none             0.5            -0.5            0.5      10   2
+# Substrate :  n-delta_S,   n-beta_S,   Layer thickness(nm), n-delta_L,   n-beta_L,  RMS roughness(nm) 
+	         6.E-06        2.e-8          0.               1.E-05      5.E-07         0.
+# Particle : n-delta_I,     n-beta_I,     Depth(nm),  n-delta_SH,   n-beta_SH        
+            6.E-04          2.e-8         0         8.E-04          2.e-8
+################################# Grid parameters ######################################################
+# Ewald mode
+  T
+# Output angle (deg) :   Two theta min-max, Alphaf min-max,  n(1),   n(2)				
+                             0.     2         0.  2           100      100
+# Output q(nm-1) :  Qx min-max, Qy min-max, Qz min-max,  n(1), n(2), n(3)
+                     -1  1      -1   1      -2   0         200   200    1
+##################################   Particle parameters #################################################
+# Number of different particle types
+1
+# Particle type,     Probability
+cylinder                1
+# Geometrical parameters : Base angle (deg),  Height ratio, Flattening,  FS-radii/R
+	                       54.73           1.		          1.          0.8    0.8
+# Shell thicknesses (nm) : dR,   dH,  dW
+                           0    0     0
+# H_uncoupled, W_uncoupled
+       T         T
+# Size of particle	: Radius(nm), R_distribution, SigmaR/R,	Rmin(nm),	Rmax(nm), nR, xR
+                           5          none            0.01       0.1         11      100   4
+# Height aspect ratio	: Height/R,   H_distribution, SigmaH/H,    Hmin/R,	Hmax/R,   nH, xH,   rho_H
+                           1          none            0.1       0.1         11       25   2      0
+# Width aspect ratio	: Width/R,   W_distribution, SigmaW/W,    Wmin/R,	Wmax/R,   nW, xW,   rho_W
+                            2          none                0.4        1         300       15   2   0
+# Orientation of particle : Zeta(deg), Z_distribution, SigmaZ(deg), Zmin(deg), Zmax(deg),	nZ, xZ
+                             0          none       20.            0             120     30   2
+##################################### Lattice parameters #################################################
+# Lattice type
+          1DDL 
+#  Interference function :    Peak position D(nm),   w(nm), Statistics,  Eta_Voigt, Size-Distance coupling, Cut-off
+                               20                     7      gau	      0.5            0            1000
+# Pair correlation function :   Density(nm-2), D1(nm),  Hard core coverage, CxR
+                                   0.007         25          0.3             1.  
+# Lattice parameters : L(1)(nm), L(2)(nm), Angle(deg,  Xi_fixed 
+                        10       10          90.           F
+                       Xi(deg), Xi_distribution, SigmaXi(deg), Ximin(deg),  Ximax(deg), nXi, xXi
+                        0           gate          20           0.            240.       3      -2
+                       Domain sizes DL(nm), DL_distribution, SigmaDL/DL, DLmin(nm), DLmax(nm), nDL, XDL
+                         20000   20000        none         0.2  0.2   200  200   10000  10000   10 10 -2 -2
+# Imperfect lattice :  Rod description,  Rod shape,  
+                        rec_ellip         cau cau
+                       Correlation lenghts(nm),  Rod orientation(deg)
+                         3000        1000            0   90                
+# Paracrystal :    Probability description
+                        ellip
+                   Disorder factors w(nm), DL-statistical distribution and rod orientation (deg)
+                       0.5  0.5   0.5  0.5
+                       cau   cau   cau   cau
+                        0   90     0    90
+# Pattern  :     Regular pattern content,   Number of particles per pattern
+	                      F                          2	
+                 Positions xp/L, Debye-Waller factors B11/L1 B22/L1 B12/L1 
+                   0.    0.     0.     0.   0.
+                   0.5  0.5     0.     0.   0.
+