Double To complex interpolating function with polar interpolation

App/src/TestMiscellaneous.cpp

@@ -7,7 +7,8 @@
 #include "TGraph.h"
 #include "TH2F.h"
 #include "TCanvas.h"
-
+#include "TGraphPolar.h"
+#include "TSystem.h"
 
 TestMiscellaneous::TestMiscellaneous()
 {
@@ -31,7 +32,7 @@ void TestMiscellaneous::test_DoubleToComplexInterpolatingFunction()
     MultiLayer *sample = dynamic_cast<MultiLayer *>(SampleFactory::instance().createItem("MultilayerOffspecTestcase1a"));
 
     OutputData<double > *data_alpha = new OutputData<double >;
-    data_alpha->addAxis(std::string("alpha_f"), 0.0*Units::degree, 2.0*Units::degree, 100);
+    data_alpha->addAxis(std::string("alpha_f"), 0.0*Units::degree, 2.0*Units::degree, 200);
 
     OpticalFresnel fresnelCalculator;
 
@@ -63,14 +64,14 @@ void TestMiscellaneous::test_DoubleToComplexInterpolatingFunction()
             R_map[angle] = R;
         }
         DoubleToComplexInterpolatingFunction T_function(T_map);
-        DoubleToComplexInterpolatingFunction R_function(R_map);
+        DoubleToComplexInterpolatingFunction R_function(R_map, DoubleToComplexInterpolatingFunction::Nearest);
 
         m_TT[i_layer] = T_function.clone();
         m_RR[i_layer] = R_function.clone();
     }
 
     double alpha_min(0), alpha_max(2.0*Units::degree);
-    int npoints = 99*5;
+    int npoints = 200*2;
 
     TGraph *gr1_exact = new TGraph(npoints);
     TGraph *gr2_interp = new TGraph(npoints);
@@ -78,35 +79,20 @@ void TestMiscellaneous::test_DoubleToComplexInterpolatingFunction()
 
     for(int i_point=0; i_point < npoints; i_point++){
         int i_layer_sel = 0;
-        double angle = alpha_min + i_point*(alpha_max-alpha_min)/double(npoints);
+        double angle = alpha_min + i_point*(alpha_max-alpha_min)/double(npoints-1);
         kvector_t kvec;
         kvec.setLambdaAlphaPhi(1.4*Units::angstrom, angle, 0.0);
         OpticalFresnel::MultiLayerCoeff_t coeffs;
         fresnelCalculator.execute(*sample, kvec, coeffs);
-//        hh2[0][0]->Fill();
         complex_t R = m_RR[i_layer_sel]->evaluate(angle);
         std::cout << i_point << " " << angle << " true R:" << coeffs[i_layer_sel].R << " interp:" << R << " " << std::abs(R - coeffs[i_layer_sel].R) << std::endl;
+        complex_t r = coeffs[i_layer_sel].R;
+//        std::cout << "RRR " << r << " abs:" << std::abs(r) << " arg:" << std::arg(r) << " " << Units::rad2deg(std::arg(r)) << std::endl << std::endl;
         gr1_exact->SetPoint(i_point, angle, std::abs(coeffs[i_layer_sel].R) );
         gr2_interp->SetPoint(i_point, angle, std::abs(R) );
         gr3_diff->SetPoint(i_point, angle, std::abs(R) - std::abs(coeffs[i_layer_sel].R) );
     }
 
-//    data_alpha->resetIndex();
-//    while (data_alpha->hasNext())
-//    {
-//        double angle = data_alpha->getCurrentValueOfAxis<double>("alpha_f") + ;
-//        kvector_t kvec;
-//        kvec.setLambdaAlphaPhi(1.4*Units::angstrom, angle, 0.0);
-//        OpticalFresnel::MultiLayerCoeff_t coeffs;
-//        fresnelCalculator.execute(*sample, kvec, coeffs);
-//        complex_t R = m_RR[0]->evaluate(angle);
-//        std::cout  << " " << angle << " true R:" << coeffs[0].R << " interp:" << R << " " << std::abs(R - coeffs[0].R) << std::endl;
-//        hh2[0][0]->Fill(angle, std::abs(R - coeffs[0].R) );
-
-//        data_alpha->next();
-//    }
-
-
 
     TCanvas *c1 = new TCanvas("c1","c1",1024, 768);
     c1->Divide(2,2);

Core/Tools/inc/DoubleToComplexInterpolatingFunction.h

@@ -21,8 +21,10 @@
 class DoubleToComplexInterpolatingFunction : public IDoubleToComplexFunction
 {
 public:
+    enum InterpolatingMode { Nearest, Linear, Polar };
+
 	virtual ~DoubleToComplexInterpolatingFunction();
-	DoubleToComplexInterpolatingFunction(const std::map<double, complex_t> &value_map);
+    DoubleToComplexInterpolatingFunction(const std::map<double, complex_t> &value_map, InterpolatingMode imode=Nearest);
 	virtual DoubleToComplexInterpolatingFunction *clone() const;
 
     virtual complex_t evaluate(double value);
@@ -34,6 +36,8 @@ protected:
 	double m_low_step;
 	double m_high_step;
 
+    InterpolatingMode m_interpolating_mode;
+
 private:
     //! copy constructor and assignment operator are hidden since there is a clone method
     DoubleToComplexInterpolatingFunction(const DoubleToComplexInterpolatingFunction &);

Core/Tools/src/DoubleToComplexInterpolatingFunction.cpp

@@ -9,12 +9,14 @@
 #include "Exceptions.h"
 #include <sstream>
 
+
+
 DoubleToComplexInterpolatingFunction::~DoubleToComplexInterpolatingFunction()
 {
 }
 
-DoubleToComplexInterpolatingFunction::DoubleToComplexInterpolatingFunction(const std::map<double, complex_t> &value_map)
-    : m_value_map(value_map)
+DoubleToComplexInterpolatingFunction::DoubleToComplexInterpolatingFunction(const std::map<double, complex_t> &value_map, InterpolatingMode imode)
+    : m_value_map(value_map), m_interpolating_mode(imode)
 {
 	m_lower_limit = (*m_value_map.begin()).first;
 	m_upper_limit = (*m_value_map.rbegin()).first;
@@ -24,7 +26,7 @@ DoubleToComplexInterpolatingFunction::DoubleToComplexInterpolatingFunction(const
 
 DoubleToComplexInterpolatingFunction* DoubleToComplexInterpolatingFunction::clone() const
 {
-    DoubleToComplexInterpolatingFunction *p_new = new DoubleToComplexInterpolatingFunction(m_value_map);
+    DoubleToComplexInterpolatingFunction *p_new = new DoubleToComplexInterpolatingFunction(m_value_map, m_interpolating_mode);
     return p_new;
 }
 
@@ -47,8 +49,24 @@ complex_t DoubleToComplexInterpolatingFunction::evaluate(double value)
     std::map<double, complex_t>::const_iterator lower_it = m_value_map.lower_bound(value);
     --lower_it;
     std::map<double, complex_t>::const_iterator upper_it = m_value_map.upper_bound(value);
-    // Linear interpolation:
-    complex_t result_difference = (*upper_it).second - (*lower_it).second;
+
     double interpolating_factor = (value - (*lower_it).first)/((*upper_it).first-(*lower_it).first);
-    return (*lower_it).second + result_difference * interpolating_factor;
+    if(m_interpolating_mode == Nearest) {
+        if( interpolating_factor < 0.5 ) {
+            return (*lower_it).second;
+        } else {
+            return (*upper_it).second;
+        }
+    } else if ( m_interpolating_mode == Linear ) {
+        complex_t result_difference = (*upper_it).second - (*lower_it).second;
+        return (*lower_it).second + result_difference * interpolating_factor;
+    } else if (m_interpolating_mode == Polar) {
+        double cabs_difference = std::abs((*upper_it).second) - std::abs((*lower_it).second);
+        double cabs = std::abs((*lower_it).second) + cabs_difference * interpolating_factor;
+        double carg_difference = std::arg((*upper_it).second) - std::arg((*lower_it).second);
+        double carg = std::arg((*lower_it).second) + carg_difference * interpolating_factor;
+        return std::polar(cabs, carg);
+    } else {
+        throw DomainErrorException("DoubleToComplexInterpolatingFunction::evaluate() -> Error. Unknown interpolation mode");
+    }
 }

Examples/Performance/perf_history.txt

@@ -33,3 +33,4 @@
 # debug is on again
 2012-08-03 09:21:48 | jcnsopc73  | macosx64, 2800 MHz      | 77519.4 | 4.04858 | 3.89864 | 0.17574 | 
 2012-08-08 11:41:47 | jcnsopc73  | macosx64, 2800 MHz      | 81967.2 | 4.04858 | 3.89105 | 0.26246 | 
+2012-08-08 17:27:53 | jcnsopc73  | macosx64, 2800 MHz      | 83333.3 | 4.07332 | 3.89864 | 0.26143 |