1 files changed, 91 insertions, 0 deletions
diff --git a/models/dbeta2d3.hpp b/models/dbeta2d3.hpp
new file mode 100644
index 0000000..e2edb78
--- /dev/null
+++ b/models/dbeta2d3.hpp
@@ -0,0 +1,91 @@
+#ifndef DDBETA3_MODEL2d2_H_
+#define DDBETA3_MODEL2d2_H_
+#include <core/fitter.hpp>
+#include <cmath>
+#include <cassert>
+#include "vecn.hpp"
+
+
+namespace opt_utilities
+{
+  
+  template <typename T>
+  class dbeta2d3
+    :public model<T,vecn<T,2>,std::vector<T>,std::string>
+  {
+  private:
+    model<T,vecn<T,2>,std::vector<T> >* do_clone()const
+    {
+      return new dbeta2d3<T>(*this);
+    }
+  public:
+    dbeta2d3()
+    {
+
+      this->push_param_info(param_info<std::vector<T> >("x0",256));//1
+      this->push_param_info(param_info<std::vector<T> >("y0",256));//2      
+      this->push_param_info(param_info<std::vector<T> >("epsilon",0));//3
+      this->push_param_info(param_info<std::vector<T> >("theta",0));//4      
+      
+      this->push_param_info(param_info<std::vector<T> >("ampl1",1));//5
+      this->push_param_info(param_info<std::vector<T> >("beta1",0.6));//6
+      this->push_param_info(param_info<std::vector<T> >("r01",30));//7
+      
+      this->push_param_info(param_info<std::vector<T> >("ampl2",.5));//8
+      this->push_param_info(param_info<std::vector<T> >("beta2",.4));//9
+      this->push_param_info(param_info<std::vector<T> >("r02",20));//10
+      
+      this->push_param_info(param_info<std::vector<T> >("bkg",0));//11
+
+
+    }
+
+
+    T do_eval(const vecn<T,2>& xy,const std::vector<T>& param)
+    {
+      T x=xy[0];
+      T y=xy[1];
+      
+      T x0=get_element(param,0);
+      T y0=get_element(param,1);
+      T epsilon=get_element(param,2);
+      T theta=get_element(param,3);
+ 
+      T ampl1=(get_element(param,4));
+      T beta1=(get_element(param,5));
+      T r01=(get_element(param,6));
+
+      T ampl2=(get_element(param,7));
+      T beta2=(get_element(param,8));
+      T r02=(get_element(param,9));
+
+      T bkg=get_element(param,10);
+
+      T x_new1=(x-x0)*cos(theta)+(y-y0)*sin(theta);
+      T y_new1=(y-y0)*cos(theta)-(x-x0)*sin(theta);
+      
+ 
+      T r1_r1=x_new1*x_new1/exp(epsilon/30)+y_new1*y_new1/exp(-epsilon/30);
+      //T r1=sqrt(x_new1*x_new1*(1-epsilon)*(1-epsilon)+y_new1*y_new1)/(1-epsilon);
+
+      T x_new2=(x-x0)*cos(theta)+(y-y0)*sin(theta);
+      T y_new2=(y-y0)*cos(theta)-(x-x0)*sin(theta);
+      
+      T r2_r2=x_new2*x_new2/exp(epsilon/30)+y_new2*y_new2/exp(-epsilon/30);
+      //T r2=sqrt(x_new2*x_new2*(1-epsilon)*(1-epsilon)+y_new2*y_new2)/(1-epsilon);
+
+      
+      
+      return bkg+ampl1*pow(1+(r1_r1/r01/r01),-3*beta1+static_cast<T>(.5))
+      	+ampl2*pow(1+(r2_r2/r02/r02),-3*beta2+static_cast<T>(.5));
+
+      //return bkg+pow(1+r1*r1/r01/r01,-3*beta1+static_cast<T>(.5))+
+      //pow(1+r2*r2/r02/r02,-3*beta2+static_cast<T>(.5));
+    }
+  };
+};
+
+
+
+#endif
+//EOF