From aebdd8834e038cfa5c3c78eeeeaa92682a0c7f60 Mon Sep 17 00:00:00 2001
From: Celine Durniak <c.durniak@fz-juelich.de>
Date: Wed, 8 May 2013 09:30:07 +0200
Subject: [PATCH] Modifications of mathematical operations with complex vectors
(BasicVector3D.cpp)
---
App/inc/TestFormFactor.h | 2 +-
App/inc/TestMiscellaneous.h | 3 +
App/src/SampleFactory.cpp | 2 +-
App/src/TestFormFactor.cpp | 77 ++++++++-
App/src/TestFormFactors.cpp | 25 ++-
App/src/TestMiscellaneous.cpp | 157 +++++++++++++++++-
Core/Core.pro | 2 +
Core/FormFactors/inc/FormFactorTethraedron.h | 53 ++++++
.../FormFactors/src/FormFactorTethraedron.cpp | 103 ++++++++++++
Core/Geometry/src/BasicVector3D.cpp | 45 +++--
10 files changed, 443 insertions(+), 26 deletions(-)
create mode 100644 Core/FormFactors/inc/FormFactorTethraedron.h
create mode 100644 Core/FormFactors/src/FormFactorTethraedron.cpp
diff --git a/App/inc/TestFormFactor.h b/App/inc/TestFormFactor.h
index 6d5eeee9cfa..b231e15e4b7 100644
--- a/App/inc/TestFormFactor.h
+++ b/App/inc/TestFormFactor.h
@@ -26,7 +26,7 @@ class TestFormFactor : public IFunctionalTest
TestFormFactor();
~TestFormFactor();
virtual void execute();
- void draw();
+ void draw(), draw4();
private:
OutputData<double> *mp_intensity_output;
diff --git a/App/inc/TestMiscellaneous.h b/App/inc/TestMiscellaneous.h
index 189c0547cd3..4c5752611e8 100644
--- a/App/inc/TestMiscellaneous.h
+++ b/App/inc/TestMiscellaneous.h
@@ -34,6 +34,9 @@ class TestMiscellaneous : public IFunctionalTest
//! form factor as a function of qx,qy,qz
void test_FormFactor();
+ //! Re, Im or Amp, phase of form factors as functions of qx,qy,qz
+ void test_FormFactor1();
+
//! opengl mesocrystal drawing
void test_DrawMesocrystal();
diff --git a/App/src/SampleFactory.cpp b/App/src/SampleFactory.cpp
index acf4d59fc24..d28e60a9ab4 100644
--- a/App/src/SampleFactory.cpp
+++ b/App/src/SampleFactory.cpp
@@ -93,7 +93,7 @@ SampleFactory::SampleFactory()
registerItem("FormFactor_Ellipsoid", StandardSamples::FormFactor_Ellipsoid);
registerItem("FormFactor_FullSpheroid", StandardSamples::FormFactor_FullSpheroid);
registerItem("FormFactor_HemiSpheroid", StandardSamples::FormFactor_HemiSpheroid);
- registerItem("FormFactor_Parallelpiped", StandardSamples::FormFactor_Parallelpiped);
+ registerItem("FormFactor_Parallelepiped", StandardSamples::FormFactor_Parallelpiped);
registerItem("FormFactor_Cylinder", StandardSamples::FormFactor_Cylinder);
registerItem("FormFactor_Pyramid", StandardSamples::FormFactor_Pyramid);
registerItem("FormFactor_FullSphere", StandardSamples::FormFactor_FullSphere);
diff --git a/App/src/TestFormFactor.cpp b/App/src/TestFormFactor.cpp
index c44069b51ed..c6193f1d46f 100644
--- a/App/src/TestFormFactor.cpp
+++ b/App/src/TestFormFactor.cpp
@@ -15,6 +15,7 @@
#include "TestFormFactor.h"
#include "Types.h"
+#include "DrawHelper.h"
#include "TCanvas.h"
#include "TH2.h"
@@ -60,9 +61,83 @@ void TestFormFactor::execute()
*it = std::pow(std::abs(m_ff.evaluate(k_i, k_f_zero_bin, alpha_i, alpha_f)),2);
++it;
}
- draw();
+ draw4();
}
+
+void TestFormFactor::draw4()
+{
+ // creation of 2D histogram from calculated intensities
+ // TCanvas *c1 = new TCanvas("c1_test_formfactor", "Cylinder Formfactor", 0, 0, 1024, 768);
+ // (void)c1;
+
+ TCanvas *c1_xy = new TCanvas("c1_test_formfactor","Cylinder Formfactor",1024,768);
+ gPad->SetRightMargin(0.11);
+ gStyle->SetPalette(1);
+ gStyle->SetOptStat(0);
+ DrawHelper::SetMagnifier(c1_xy);
+ c1_xy->Divide(2,2);
+
+ const IAxis *p_y_axis = mp_intensity_output->getAxis("detector y-axis");
+ const IAxis *p_z_axis = mp_intensity_output->getAxis("detector z-axis");
+ size_t y_size = p_y_axis->getSize();
+ size_t z_size = p_z_axis->getSize();
+ double y_start = (*p_y_axis)[0];
+ double y_end = (*p_y_axis)[y_size-1];
+ double z_start = (*p_z_axis)[0];
+ double z_end = (*p_z_axis)[z_size-1];
+
+ std::vector<TH2D *> p_hist2Da;
+ p_hist2Da.resize(4, 0);
+ /* p_hist2Da[0] = new TH2D("p_hist2D", "Cylinder Formfactor Ampl", (int)y_size,
+ y_start, y_end, (int)z_size, z_start, z_end);
+ p_hist2Da[1] = new TH2D("p_hist2D", "Cylinder Formfactor Re", (int)y_size,
+ y_start, y_end, (int)z_size, z_start, z_end);
+ p_hist2Da[2] = new TH2D("p_hist2D", "Cylinder Formfactor Im", (int)y_size,
+ y_start, y_end, (int)z_size, z_start, z_end);*/
+
+ for (int i=0; i<4;i++) {
+ //p_hist2Da[i]->UseCurrentStyle();
+ p_hist2Da[i] = new TH2D("p_hist2D", "Cylinder Formfactor", (int)y_size,
+ y_start, y_end, (int)z_size, z_start, z_end);
+ p_hist2Da[i]->GetXaxis()->SetTitle("phi_{f} (rad)");
+ p_hist2Da[i]->GetYaxis()->SetTitle("alpha_{f} (rad)");
+ }
+
+ TH2D *p_hist2D = new TH2D("p_hist2D", "Cylinder Formfactor", (int)y_size, y_start, y_end, (int)z_size, z_start, z_end);
+ //p_hist2D->UseCurrentStyle();
+ p_hist2D->GetXaxis()->SetTitle("phi_{f} (rad)");
+ p_hist2D->GetYaxis()->SetTitle("alpha_{f} (rad)");
+
+ OutputData<double>::const_iterator it = mp_intensity_output->begin();
+ //OutputData<double>::const_iterator it = mp_intensity_output->begin();
+ while (it != mp_intensity_output->end())
+ {
+ size_t index_y = mp_intensity_output->getIndexOfAxis("detector y-axis", it.getIndex());
+ size_t index_z = mp_intensity_output->getIndexOfAxis("detector z-axis", it.getIndex());
+ double x_value = (*p_y_axis)[index_y];
+ double y_value = (*p_z_axis)[index_z];
+ double z_value = std::log(*it);
+ p_hist2D->Fill(x_value, y_value, z_value);
+ p_hist2Da[0]->Fill(x_value, y_value, z_value);
+ ++it;
+ }
+
+
+ c1_xy->cd(1);
+ p_hist2D->SetContour(50);
+ p_hist2D->Draw("CONT4 z");
+
+ c1_xy->cd(2);
+
+ p_hist2Da[0]->SetContour(99);
+ p_hist2Da[0]->Draw("cont4 z");
+
+
+
+}
+
+
void TestFormFactor::draw()
{
// creation of 2D histogram from calculated intensities
diff --git a/App/src/TestFormFactors.cpp b/App/src/TestFormFactors.cpp
index 904577aede9..fc086adaa90 100644
--- a/App/src/TestFormFactors.cpp
+++ b/App/src/TestFormFactors.cpp
@@ -84,6 +84,13 @@ void TestFormFactors::execute()
simulation.runSimulation();
OutputDataIOFactory::writeOutputData(*simulation.getOutputData(), Utils::FileSystem::GetHomePath()+"./Examples/FormFactors/this_FullSphere.ima");
+ //Parallelepiped
+ sample = dynamic_cast<MultiLayer *>(SampleFactory::createSample("FormFactor_Parallelepiped"));
+ simulation.setSample(*sample);
+ simulation.runSimulation();
+ OutputDataIOFactory::writeOutputData(*simulation.getOutputData(), Utils::FileSystem::GetHomePath()+"./Examples/FormFactors/this_Parallelepiped.ima");
+
+
//Prism3
sample = dynamic_cast<MultiLayer *>(SampleFactory::createSample("FormFactor_Prism3"));
simulation.setSample(*sample);
@@ -95,21 +102,21 @@ void TestFormFactors::finalise()
{
std::vector<std::string > this_files;
// this_files.push_back(Utils::FileSystem::GetHomePath()+"./Examples/FormFactors/this_Box.ima");
-// this_files.push_back(Utils::FileSystem::GetHomePath()+"./Examples/FormFactors/this_Cone.ima");
- this_files.push_back(Utils::FileSystem::GetHomePath()+"./Examples/FormFactors/this_Sphere.ima");
-// this_files.push_back(Utils::FileSystem::GetHomePath()+"./Examples/FormFactors/this_Ellipsoid.ima");
-// this_files.push_back(Utils::FileSystem::GetHomePath()+"./Examples/FormFactors/this_FullSpheroid.ima");
-// this_files.push_back(Utils::FileSystem::GetHomePath()+"./Examples/FormFactors/this_HemiSpheroid.ima");
-// this_files.push_back(Utils::FileSystem::GetHomePath()+"./Examples/FormFactors/this_Parallelpiped.ima");
+ this_files.push_back(Utils::FileSystem::GetHomePath()+"./Examples/FormFactors/this_Cone.ima");
+// this_files.push_back(Utils::FileSystem::GetHomePath()+"./Examples/FormFactors/this_Sphere.ima");
+ this_files.push_back(Utils::FileSystem::GetHomePath()+"./Examples/FormFactors/this_Ellipsoid.ima");
+ this_files.push_back(Utils::FileSystem::GetHomePath()+"./Examples/FormFactors/this_FullSpheroid.ima");
+ this_files.push_back(Utils::FileSystem::GetHomePath()+"./Examples/FormFactors/this_HemiSpheroid.ima");
+// this_files.push_back(Utils::FileSystem::GetHomePath()+"./Examples/FormFactors/this_Parallelepiped.ima");
// this_files.push_back(Utils::FileSystem::GetHomePath()+"./Examples/FormFactors/this_Cylinder.ima");
// this_files.push_back(Utils::FileSystem::GetHomePath()+"./Examples/FormFactors/this_Pyramid.ima");
- this_files.push_back(Utils::FileSystem::GetHomePath()+"./Examples/FormFactors/this_FullSphere.ima");
+// this_files.push_back(Utils::FileSystem::GetHomePath()+"./Examples/FormFactors/this_FullSphere.ima");
// this_files.push_back(Utils::FileSystem::GetHomePath()+"./Examples/FormFactors/this_Prism3.ima");
int ncomparison = (int)this_files.size();
TCanvas *c1 = DrawHelper::createAndRegisterCanvas("Form Factors", "TestFormFactors");
- c1->Divide(4,5);
-
+ // c1->Divide(4,5);
+ c1->Divide(2,2);
for(int i=0; i<ncomparison; i++) {
OutputData<double> *our_data = OutputDataIOFactory::getOutputData(this_files[i]);
diff --git a/App/src/TestMiscellaneous.cpp b/App/src/TestMiscellaneous.cpp
index 98b5d6daf8b..3e97fa52741 100644
--- a/App/src/TestMiscellaneous.cpp
+++ b/App/src/TestMiscellaneous.cpp
@@ -41,6 +41,7 @@
#include "TGraphPolar.h"
#include "TRandom.h"
#include "TBenchmark.h"
+#include "TStyle.h"
TestMiscellaneous::TestMiscellaneous()
{
@@ -48,12 +49,13 @@ TestMiscellaneous::TestMiscellaneous()
void TestMiscellaneous::execute()
{
- test_LogSystem();
+ //test_LogSystem();
//test_OutputDataTo2DArray();
//test_KVectorContainer();
//test_OutputDataIOFactory();
//test_FastSin();
//test_DoubleToComplexInterpolatingFunction();
+ test_FormFactor1();
//test_FormFactor();
//test_DrawMesocrystal();
}
@@ -323,6 +325,159 @@ void TestMiscellaneous::test_FormFactor()
}
}
+/* ************************************************************************* */
+// plots of form factor :
+// contourplot (amp & phase or Re & Im) and log(|F|**2) vs. log(q)
+/* ************************************************************************* */
+void TestMiscellaneous::test_FormFactor1()
+{
+ FormFactorFullSphere ff_fullsphere(5.*Units::nanometer);
+
+ FormFactorCylinder ff_cylinder(10.*Units::nanometer,
+ 5.*Units::nanometer);
+ // IFormFactor& ff = ff_cylinder;
+
+ FormFactorParallelepiped ff_para(7.*Units::nanometer,
+ 6.*Units::nanometer);
+ // IFormFactor& ff = ff_para;
+
+ FormFactorPyramid ff_pyramid(10.*Units::nanometer,
+ 5.*Units::nanometer,
+ Units::deg2rad(54.73 ));
+ // IFormFactor& ff = ff_pyramid;
+
+ FormFactorPrism3 ff_prism3(5.*Units::nanometer,
+ 5.*Units::nanometer);
+ // IFormFactor& ff = ff_prism3;
+
+ FormFactorSphere ff_sphere(5.*Units::nanometer,
+ 5.*Units::nanometer);
+ // IFormFactor& ff = ff_sphere;
+
+ FormFactorBox ff_box(5*Units::nanometer,
+ 5*Units::nanometer,
+ 5*Units::nanometer);
+ // IFormFactor& ff = ff_box;
+
+ IFormFactor& ff = ff_fullsphere;
+
+ double qmin(-4.0), qmax(4.0);
+ double lambda = 1.0;
+ double alpha_i = 0.2*M_PI/180.0;
+ cvector_t k_i;
+ k_i.setLambdaAlphaPhi(lambda, -alpha_i, 0.0);
+ //double phi_fmin(-4.0), phi_fmax(4.0), alpha_fmin(-4.), alpha_fmax(4.);
+ int nbins(101);
+ double dq =(qmax-qmin)/(nbins-1);
+ //double dphi_f = (phi_fmax-phi_fmin)/(nbins-1);
+ //double dalpha_f = (alpha_fmax-alpha_fmin)/(nbins-1);
+
+ TH2D *vh2_xy = new TH2D("vh2_xy","vh2_xy;q_{x};q_{y};qz",nbins, qmin-dq/2., qmax+dq/2., nbins, qmin-dq/2., qmax+dq/2.);
+
+ OutputData<double> *p_data = new OutputData<double>();
+ p_data->addAxis(std::string("qx"), nbins, qmin, qmax);
+ p_data->addAxis(std::string("qy"), nbins, qmin, qmax);
+ p_data->addAxis(std::string("qz"), 1, qmin, qmax);
+ OutputData<double>::const_iterator it = p_data->begin();
+ double z = p_data->getValueOfAxis("qz", it.getIndex());
+
+ while (it != p_data->end()) {
+ double x = p_data->getValueOfAxis("qx", it.getIndex());
+ double y = p_data->getValueOfAxis("qy", it.getIndex());
+
+ cvector_t q(x,y,z);
+ cvector_t q0(0,0,0);
+ Bin1DCVector q0_bin(q0, q0);
+ double value = std::abs(ff.evaluate(q,q0_bin, 0.0, 0.0));
+ //double valuep = std::abs(ff.evaluate(q,q0_bin, 0.0, 0.0));
+ //double valuer = std::abs(ff.evaluate(q,q0_bin, 0.0, 0.0));
+ //double valuei = std::abs(ff.evaluate(q,q0_bin, 0.0, 0.0));
+
+ vh2_xy->Fill(x,y,value);
+
+ ++it;
+ }
+
+ TCanvas *c1_xy = new TCanvas("c1_xy","c1_xy",1024,768);
+ DrawHelper::SetMagnifier(c1_xy);
+ //c1_xy->Divide(2,2);
+ c1_xy->cd(1);
+ gPad->SetRightMargin(0.11);
+ gPad->SetLogz();
+ vh2_xy->GetXaxis()->SetNdivisions(510);
+ vh2_xy->GetYaxis()->SetNdivisions(510);
+ vh2_xy->SetContour(99);
+ gStyle->SetPalette(1);
+ gStyle->SetOptStat(0);
+ vh2_xy->Draw("cont4 z");
+ c1_xy->Print("test.eps");
+
+ /*
+ OutputData<double>::iterator it = mp_intensity_output->begin();
+ const IAxis *p_y_axis = mp_intensity_output->getAxis("detector y-axis");
+ const IAxis *p_z_axis = mp_intensity_output->getAxis("detector z-axis");
+ cvector_t k_i;
+ k_i.setLambdaAlphaPhi(lambda, -alpha_i, 0.0);
+ while (it != mp_intensity_output->end())
+ {
+ size_t index_y = mp_intensity_output->getIndexOfAxis("detector y-axis", it.getIndex());
+ size_t index_z = mp_intensity_output->getIndexOfAxis("detector z-axis", it.getIndex());
+ double phi_f = M_PI*(*p_y_axis)[index_y]/180.0;
+ double alpha_f = M_PI*(*p_z_axis)[index_z]/180.0;
+ cvector_t k_f;
+ k_f.setLambdaAlphaPhi(lambda, alpha_f, phi_f);
+ Bin1DCVector k_f_zero_bin(k_f, k_f);
+ *it = std::pow(std::abs(m_ff.evaluate(k_i, k_f_zero_bin, alpha_i, alpha_f)),2);
+ ++it;
+ }
+
+ const IAxis *p_y_axis = mp_intensity_output->getAxis("detector y-axis");
+ const IAxis *p_z_axis = mp_intensity_output->getAxis("detector z-axis");
+ size_t y_size = p_y_axis->getSize();
+ size_t z_size = p_z_axis->getSize();
+ double y_start = (*p_y_axis)[0];
+ double y_end = (*p_y_axis)[y_size-1];
+ double z_start = (*p_z_axis)[0];
+ double z_end = (*p_z_axis)[z_size-1];
+ p_hist2D->UseCurrentStyle();
+ p_hist2D->GetXaxis()->SetTitle("phi_f");
+ p_hist2D->GetYaxis()->SetTitle("alpha_f");
+*/
+
+ /*
+ c1->Divide(2,2);
+
+ c1->cd(1); gPad->SetLogz();
+ IsGISAXSTools::setMinimum(hmin);
+ if(hmax>0) IsGISAXSTools::setMaximum(hmax);
+
+ TH1 *hist(0);
+ hist = IsGISAXSTools::getOutputDataTH2D(output, "p_hist1D");
+ if( hasMinimum() ) hist->SetMinimum(m_hist_min);
+ if( hasMaximum() ) hist->SetMaximum(m_hist_max);
+ hist->SetTitle(histogram_title.c_str());
+ hist->DrawCopy(draw_options.c_str());
+
+ // dealing with masks
+ if(output.getMask()) {
+ TPolyMarker *poly = new TPolyMarker();
+ const IAxis *p_axis0 = output.getAxis(0);
+ const IAxis *p_axis1 = output.getAxis(1);
+ int i_point(0);
+ for(OutputData<double>::const_iterator it = output.begin();
+ it!= output.end(); ++it) {
+ size_t axis0_index = output.toCoordinate(it.getIndex(), 0);
+ size_t axis1_index = output.toCoordinate(it.getIndex(), 1);
+ double axis0_value = (*p_axis0)[axis0_index];
+ double axis1_value = (*p_axis1)[axis1_index];
+ poly->SetPoint(i_point++, axis0_value, axis1_value);
+ }
+ poly->Draw("same");
+ }
+*/
+
+}
+
/* ************************************************************************* */
// test double to complex interpolating function
/* ************************************************************************* */
diff --git a/Core/Core.pro b/Core/Core.pro
index 4e5d541cffe..2b9e9efccf3 100644
--- a/Core/Core.pro
+++ b/Core/Core.pro
@@ -92,6 +92,7 @@ SOURCES += \
FormFactors/src/FormFactorPrism6.cpp \
FormFactors/src/FormFactorPyramid.cpp \
FormFactors/src/FormFactorSphere.cpp \
+ FormFactors/src/FormFactorTethraedron.cpp \
FormFactors/src/FormFactorWeighted.cpp \
FormFactors/src/IFormFactorBorn.cpp \
\
@@ -246,6 +247,7 @@ HEADERS += \
FormFactors/inc/FormFactorPyramid.h \
FormFactors/inc/FormFactorSphere.h \
FormFactors/inc/FormFactorSphereGaussianRadius.h \
+ FormFactors/inc/FormFactorTethraedron.h \
FormFactors/inc/FormFactorWeighted.h \
FormFactors/inc/FormFactors.h \
FormFactors/inc/IFormFactor.h \
diff --git a/Core/FormFactors/inc/FormFactorTethraedron.h b/Core/FormFactors/inc/FormFactorTethraedron.h
new file mode 100644
index 00000000000..23f647fcbbb
--- /dev/null
+++ b/Core/FormFactors/inc/FormFactorTethraedron.h
@@ -0,0 +1,53 @@
+// ************************************************************************** //
+//
+// BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file FormFactors/inc/FormFactorTethraedron.h
+//! @brief Defines class FormFactorTethraedron
+//!
+//! @homepage http://apps.jcns.fz-juelich.de/BornAgain
+//! @license GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2013
+//! @authors Scientific Computing Group at MLZ Garching
+//! @authors C. Durniak, G. Pospelov, W. Van Herck, J. Wuttke
+//
+// ************************************************************************** //
+
+#ifndef FORMFACTORTETHRAEDRON_H
+#define FORMFACTORTETHRAEDRON_H
+
+
+#include "IFormFactorBorn.h"
+#include "IStochasticParameter.h"
+
+//! Form factor of tethraedron.
+
+class FormFactorTethraedron : public IFormFactorBorn
+{
+ public:
+ //! @brief tethraedron constructor
+ //! @param height of tethraedron
+ //! @param half_side half of tethraedron's base
+ //! @param angle in radians between base and facet
+ FormFactorTethraedron(double height, double half_side, double alpha);
+
+ ~FormFactorTethraedron() {}
+ virtual FormFactorTethraedron *clone() const;
+
+ virtual int getNumberOfStochasticParameters() const { return 3; }
+
+ virtual double getHeight() const { return m_height; }
+
+ virtual complex_t evaluate_for_q(const cvector_t& q) const;
+
+ protected:
+ virtual void init_parameters();
+
+ private:
+ double m_height;
+ double m_half_side;
+ double m_alpha;
+ double m_root3; // Cached value of square root of 3
+};
+
+#endif // FORMFACTORTETHRAEDRON_H
diff --git a/Core/FormFactors/src/FormFactorTethraedron.cpp b/Core/FormFactors/src/FormFactorTethraedron.cpp
new file mode 100644
index 00000000000..f24cbc0b594
--- /dev/null
+++ b/Core/FormFactors/src/FormFactorTethraedron.cpp
@@ -0,0 +1,103 @@
+// ************************************************************************** //
+//
+// BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file FormFactors/src/FormFactorTethraedron.cpp
+//! @brief Implements class FormFactorTethraedron.
+//!
+//! @homepage http://apps.jcns.fz-juelich.de/BornAgain
+//! @license GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2013
+//! @authors Scientific Computing Group at MLZ Garching
+//! @authors C. Durniak, G. Pospelov, W. Van Herck, J. Wuttke
+//
+// ************************************************************************** //
+
+#include "FormFactorTethraedron.h"
+#include "StochasticDiracDelta.h"
+#include "MathFunctions.h"
+
+FormFactorTethraedron::FormFactorTethraedron(
+ double height, double half_side, double alpha)
+{
+ setName("FormFactorTethraerdon");
+ m_height = height;
+ m_half_side = half_side;
+ m_alpha = alpha;
+ m_root3 = std::sqrt(3.0);
+ init_parameters();
+}
+
+void FormFactorTethraedron::init_parameters()
+{
+ getParameterPool()->clear();
+ getParameterPool()->registerParameter("height", &m_height);
+ getParameterPool()->registerParameter("half_side", &m_half_side);
+ getParameterPool()->registerParameter("alpha", &m_alpha);
+}
+
+FormFactorTethraedron* FormFactorTethraedron::clone() const
+{
+ FormFactorTethraedron *p_clone =
+ new FormFactorTethraedron(m_height, m_half_side, m_alpha);
+ return p_clone;
+}
+
+complex_t FormFactorTethraedron::evaluate_for_q(const cvector_t& q) const
+{
+
+ double H = m_height;
+ double R = m_half_side;
+ double tga = std::tan(m_alpha);
+
+ complex_t qx = q.x();
+ complex_t qy = q.y();
+ complex_t qz = q.z();
+
+ complex_t F;
+ const complex_t im(0,1);
+
+ if (std::abs(qx)==0.0 && std::abs(qy)==0.0) {
+ complex_t qzH_half = qz*H/2.0;
+ F = m_root3*R*R*H*std::exp(im*qzH_half)*MathFunctions::Sinc(qzH_half);
+ }
+ else {
+ if (std::abs(qx*qx-3.0*qy*qy)==0.0) {
+ complex_t qa = 2.*qy/tga + qz/2.;
+ complex_t qb = - qy/tga + qz/2.;
+ F = H/2.*m_root3*std::exp(im*2.*qy*R/m_root3)*(
+ - std::exp(im*qa*H-im*2.*m_root3*qy*R)*MathFunctions::Sinc(qa*H)
+ + std::exp(im*qb*H)*MathFunctions::Sinc(qb*H)*
+ ( 1.- 3.*qy/(tga*qb) - im*2.*m_root3*qy*R )
+ + 3.*qy/(qb*tga)*std::exp(im*2.*qb*H)
+ )/( qx*qx ) ;
+ } else {
+
+ complex_t q1, q2, q3;
+ double L;
+ L = 2.*tga*R/m_root3-H;
+ q3 = (qy/tga - qz/2.);
+
+ if (std::abs(qx)==0.0) {
+
+ F = - 2.*H*MathFunctions::Sinc(q3*H)*
+ std::exp(im*(q3*L+qz*R*tga/m_root3))/((m_root3*qy*qy));
+
+ } else {
+
+ q1=(1./2.)*((m_root3*qx-qy)/tga - qz);
+ q2=(1./2.)*((m_root3*qx+qy)/tga + qz);
+
+ F = -(1. + m_root3*qy/qx)*MathFunctions::Sinc(q1*H)*std::exp(im*q1*L)
+ -(1. - m_root3*qy/qx)*MathFunctions::Sinc(q2*H)*std::exp(-im*q2*L)
+ + 2.*MathFunctions::Sinc(q3*H)*std::exp(im*q3*L);
+
+ F = F*H*std::exp(im*qz*R*tga/m_root3)*m_root3/((qx*qx-3.*qy*qy));}
+ }
+ }
+
+ return F;
+}
+
+
+
diff --git a/Core/Geometry/src/BasicVector3D.cpp b/Core/Geometry/src/BasicVector3D.cpp
index 2597f2838e6..0476fda2d16 100644
--- a/Core/Geometry/src/BasicVector3D.cpp
+++ b/Core/Geometry/src/BasicVector3D.cpp
@@ -33,11 +33,11 @@ double BasicVector3D<double>::mag2() const
return x()*x()+y()*y()+z()*z();
}
-//! @TODO eliminate this, it is plain wrong
+//! Returns squared magnitude of the vector.
template<>
complex_t BasicVector3D<complex_t>::mag2() const
{
- return x()*x()+y()*y()+z()*z();
+ return x()*std::conj(x())+y()*std::conj(y())+z()*std::conj(z());
}
//! Returns magnitude of the vector.
@@ -65,11 +65,11 @@ double BasicVector3D<double>::magxy2() const
return x()*x()+y()*y();
}
-//! @TODO eliminate this, it is plain wrong
+//! Returns squared distance from z axis.
template<>
complex_t BasicVector3D<complex_t>::magxy2() const
{
- return x()*x()+y()*y();
+ return x()*std::conj(x())+y()*std::conj(y());
}
//! Returns distance from z axis.
@@ -79,7 +79,7 @@ double BasicVector3D<double>::magxy() const
return std::sqrt(magxy2());
}
-//! @TODO eliminate this, it is plain wrong
+//! Returns distance from z axis.
template<>
complex_t BasicVector3D<complex_t>::magxy() const
{
@@ -109,28 +109,47 @@ double BasicVector3D<double>::cosTheta() const
return std::abs(ma) == 0 ? 1 : z()/ma;
}
+
//! @TODO eliminate this, it is plain wrong
-template<>
+ template<>
complex_t BasicVector3D<complex_t>::cosTheta() const
-{
+ {
complex_t ma = mag();
- return std::abs(ma) == 0 ? 1 : z()/ma;
+ return std::abs(ma) == 0 ? 1 : abs(z())/ma;
}
+// //! @TODO eliminate this, it is plain wrong
+// template<>
+ //complex_t BasicVector3D<complex_t>::cosTheta() const
+// {
+// complex_t ma = mag();
+// return std::abs(ma) == 0 ? 1 : z()/ma;
+//}
+
// -----------------------------------------------------------------------------
// Combine two vectors
// -----------------------------------------------------------------------------
-//! @TODO check usage: unlikely to be correct
-template<>
-BasicVector3D<complex_t> BasicVector3D<complex_t>::cross(
- const BasicVector3D<complex_t>& v) const
+ //! Returns cross product of vectors
+ template<>
+ BasicVector3D<double> BasicVector3D<double>::cross(
+ const BasicVector3D<double>& v) const
{
- return BasicVector3D<complex_t>(y()*v.z()-v.y()*z(),
+ return BasicVector3D<double>(y()*v.z()-v.y()*z(),
z()*v.x()-v.z()*x(),
x()*v.y()-v.x()*y());
}
+// //! @TODO check usage: unlikely to be correct
+// template<>
+// BasicVector3D<complex_t> BasicVector3D<complex_t>::cross(
+// const BasicVector3D<complex_t>& v) const
+//{
+// return BasicVector3D<complex_t>(y()*v.z()-v.y()*z(),
+// z()*v.x()-v.z()*x(),
+// x()*v.y()-v.x()*y());
+//}
+
//! Returns square of transverse component with respect to given axis.
template<>
double BasicVector3D<double>::perp2(const BasicVector3D<double>& v) const
--
GitLab