git-svn-id: file:///home/svn/opt_utilities@1 ed2142bd-67ad-457f-ba7c-d818d4011675

28 files changed, 2144 insertions, 0 deletions

diff --git a/models/add_model.hpp b/models/add_model.hpp
new file mode 100644
index 0000000..ffac9dc
--- /dev/null
+++ b/models/add_model.hpp
@@ -0,0 +1,169 @@
+#ifndef ADD_MODEL_H_
+#define ADD_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+
+namespace opt_utilities
+{
+  template <typename Ty,typename Tx,typename Tp,typename Tstr>
+  class add_model
+    :public model<Ty,Tx,Tp,Tstr>
+  {
+  private:
+    model<Ty,Tx,Tp,Tstr>* do_clone()const
+    {
+      return new add_model<Ty,Tx,Tp,Tstr>(*this);
+    }
+  private:
+    add_model()
+    {
+    }
+
+  private:
+    model<Ty,Tx,Tp,Tstr>* pm1;
+    model<Ty,Tx,Tp,Tstr>* pm2;
+
+  public:
+    add_model(const model<Ty,Tx,Tp,Tstr>& m1,
+	      const model<Ty,Tx,Tp,Tstr>& m2)
+      :pm1(m1.clone()),pm2(m2.clone())
+    {
+      int np1=m1.get_num_params();
+      int np2=m2.get_num_params();
+      for(int i=0;i<np1;++i)
+	{
+	  param_info<Tp,Tstr> p(m1.get_param_info(i));
+	  param_info<Tp,Tstr> p1(p.get_name()+"1",p.get_default_value());
+	  this->push_param_info(p1);
+	}
+      for(int i=0;i<np2;++i)
+	{
+	  param_info<Tp,Tstr> p(m2.get_param_info(i));
+	  param_info<Tp,Tstr> p2(p.get_name()+"2",p.get_default_value());
+	  this->push_param_info(p2);
+	}
+    }
+
+    add_model(const add_model& rhs)
+      :pm1(NULL),pm2(NULL)
+    {
+      int np1(0),np2(0);
+      if(rhs.pm1)
+	{
+	  pm1=rhs.pm1->clone();
+	  np1=rhs.pm1->get_num_params();
+	  for(int i=0;i<np1;++i)
+	    {
+	      param_info<Tp,Tstr> p(rhs.pm1->get_param_info(i));
+	      param_info<Tp,Tstr> p1(p.get_name()+"1",p.get_default_value());
+	      this->push_param_info(p1);
+	    }
+	}
+      if(rhs.pm2)
+	{
+	  pm2=rhs.pm2->clone();
+	  np2=rhs.pm2->get_num_params();
+	  for(int i=0;i<np2;++i)
+	    {
+	      param_info<Tp,Tstr> p(rhs.pm2->get_param_info(i));
+	      param_info<Tp,Tstr> p2(p.get_name()+"2",p.get_default_value());
+	      this->push_param_info(p2);
+	    }
+	}
+    }
+
+    add_model& operator=(const add_model& rhs)
+    {
+      if(this==&rhs)
+	{
+	  return *this;
+	}
+      if(!pm1)
+	{
+	  //delete pm1;
+	  pm1->destroy();
+	}
+      if(!pm2)
+	{
+	  //delete pm2;
+	  pm2->destroy();
+	}
+      int np1(0),np2(0);
+      if(rhs.pm1)
+	{
+	  pm1=rhs.pm1->clone();
+	  np1=rhs.pm1->get_num_params();
+	  for(int i=0;i<np1;++i)
+	    {
+	      param_info<Tp,Tstr> p(rhs.pm1->get_param_info(i));
+	      param_info<Tp,Tstr> p1(p.get_name()+"1",p.get_default_value());
+	      this->push_param_info(p1);
+	    }
+	}
+      if(rhs.pm2)
+	{
+	  pm2=rhs.pm2->clone();
+	  np2=rhs.pm2->get_num_params();
+	  for(int i=0;i<np2;++i)
+	    {
+	      param_info<Tp,Tstr> p(rhs.pm2->get_param_info(i));
+	      param_info<Tp,Tstr> p2(p.get_name()+"2",p.get_default_value());
+	      this->push_param_info(p2);
+	    }
+	}
+      return *this;
+    }
+
+    ~add_model()
+    {
+      if(!pm1)
+	{
+	  //delete pm1;
+	  pm1->destroy();
+	}
+      if(!pm2)
+	{
+	  //delete pm2;
+	  pm2->destroy();
+	}
+    }
+
+  public:
+    Ty do_eval(const Tx& x,const Tp& param)
+    {
+      if(!pm1)
+	{
+	  throw opt_exception("incomplete model!");
+	}
+      if(!pm2)
+	{
+	  throw opt_exception("incomplete model!");
+	}
+      Tp p1(pm1->get_num_params());
+      Tp p2(pm2->get_num_params());
+      int i=0;
+      int j=0;
+      for(i=0;i<pm1->get_num_params();++i,++j)
+	{
+	  set_element(p1,i,get_element(param,j));
+	}
+      for(i=0;i<pm2->get_num_params();++i,++j)
+	{
+	  set_element(p2,i,get_element(param,j));
+	}
+      return pm1->eval(x,p1)+pm2->eval(x,p2);
+    }
+  };
+  
+  template <typename Ty,typename Tx,typename Tp,typename Tstr>
+  add_model<Ty,Tx,Tp,Tstr> operator+(const model<Ty,Tx,Tp,Tstr>& m1,
+				     const model<Ty,Tx,Tp,Tstr>& m2)
+  {
+    return add_model<Ty,Tx,Tp,Tstr>(m1,m2);
+  }
+};
+
+
+
+#endif
+//EOF
diff --git a/models/beta1d.hpp b/models/beta1d.hpp
new file mode 100644
index 0000000..e18aa34
--- /dev/null
+++ b/models/beta1d.hpp
@@ -0,0 +1,49 @@
+#ifndef BETA_MODEL_H_
+#define BETA_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+#include <iostream>
+
+namespace opt_utilities
+{
+  template <typename T>
+  class beta1d
+    :public model<T,T,std::vector<T>,std::string>
+  {
+  private:
+    model<T,T,std::vector<T> >* do_clone()const
+    {
+      return new beta1d<T>(*this);
+    }
+  public:
+    beta1d()
+    {
+      this->push_param_info(param_info<std::vector<T> >("S0",1));
+      this->push_param_info(param_info<std::vector<T> >("rc",10));
+      this->push_param_info(param_info<std::vector<T> >("beta",2./3.));
+      this->push_param_info(param_info<std::vector<T> >("bkg",0));
+    }
+
+
+    T do_eval(const T& x,const std::vector<T>& param)
+    {
+      T S0=get_element(param,0);
+      T r_c=get_element(param,1);
+      T beta=get_element(param,2);
+      T bkg=get_element(param,3);
+
+      return bkg+S0*pow(1+(x*x)/(r_c*r_c),-3*beta+static_cast<T>(.5));
+    }
+
+    std::string do_to_string()const
+    {
+      return "Beta model\n"
+	"S=S0*(1+(r/rc)^2)^(-3*beta+0.5)\n";
+    }
+  };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/beta2d.hpp b/models/beta2d.hpp
new file mode 100644
index 0000000..abb6da1
--- /dev/null
+++ b/models/beta2d.hpp
@@ -0,0 +1,60 @@
+#ifndef BETA_MODEL2d_H_
+#define BETA_MODEL2d_H_
+#include <core/fitter.hpp>
+#include <cmath>
+#include <cassert>
+#include "vecn.hpp"
+
+
+namespace opt_utilities
+{
+  
+  template <typename T>
+  class beta2d
+    :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 beta2d<T>(*this);
+    }
+  public:
+    beta2d()
+    {
+      this->push_param_info(param_info<std::vector<T> >("S0",1));
+      this->push_param_info(param_info<std::vector<T> >("rc1",100));
+      this->push_param_info(param_info<std::vector<T> >("rc2",100));
+      this->push_param_info(param_info<std::vector<T> >("rho",0));
+      this->push_param_info(param_info<std::vector<T> >("x0",100));
+      this->push_param_info(param_info<std::vector<T> >("y0",100));
+      this->push_param_info(param_info<std::vector<T> >("beta",2./3.));
+      this->push_param_info(param_info<std::vector<T> >("bkg",0));
+    }
+
+
+    T do_eval(const vecn<T,2>& xy,const std::vector<T>& param)
+    {
+      T S0=get_element(param,0);
+      T rc1=get_element(param,1);
+      T rc2=get_element(param,2);
+      T rho=get_element(param,3);
+      T x0=get_element(param,4);
+      T y0=get_element(param,5);
+      T beta=get_element(param,6);
+      T bkg=get_element(param,7);
+      T x=xy[0];
+      T y=xy[1];
+
+      T r=(x-x0)*(x-x0)/(rc1*rc1)+(y-y0)*(y-y0)/(rc2*rc2)
+	-2*rho*(x-x0)*(y-y0)/(rc1*rc2);
+      r=r<0?0:r;
+      
+      return bkg+S0*pow(1+r,-3*beta+static_cast<T>(.5));
+    }
+  };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/beta2d2.hpp b/models/beta2d2.hpp
new file mode 100644
index 0000000..bb473a3
--- /dev/null
+++ b/models/beta2d2.hpp
@@ -0,0 +1,73 @@
+#ifndef BETA_MODEL2d2_H_
+#define BETA_MODEL2d2_H_
+#include <core/fitter.hpp>
+#include <cmath>
+#include <cassert>
+#include "vecn.hpp"
+
+
+namespace opt_utilities
+{
+  
+  template <typename T>
+  class beta2d2
+    :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 beta2d2<T>(*this);
+    }
+  public:
+    beta2d2()
+    {
+      this->push_param_info(param_info<std::vector<T> >("r0",20));
+      this->push_param_info(param_info<std::vector<T> >("x0",100));
+      this->push_param_info(param_info<std::vector<T> >("y0",100));
+      this->push_param_info(param_info<std::vector<T> >("epsilon",0));
+      this->push_param_info(param_info<std::vector<T> >("theta",100));
+      this->push_param_info(param_info<std::vector<T> >("ampl",3));
+      this->push_param_info(param_info<std::vector<T> >("beta",2./3.));
+      this->push_param_info(param_info<std::vector<T> >("bkg",0));
+    }
+
+
+    T do_eval(const vecn<T,2>& xy,const std::vector<T>& param)
+    {
+      T x=xy[0];
+      T y=xy[1];
+
+      T r0=get_element(param,0);
+      T x0=get_element(param,1);
+      T y0=get_element(param,2);
+      T epsilon=get_element(param,3);
+      T theta=get_element(param,4);
+      T ampl=get_element(param,5);
+      T beta=get_element(param,6);
+      T bkg=get_element(param,7);
+
+      T x_new=(x-x0)*cos(theta)+(y-y0)*sin(theta);
+      T y_new=(y-y0)*cos(theta)-(x-x0)*sin(theta);
+      
+      //T _epsilon=sin(epsilon)-0.00001;
+      
+      T _epsilon=epsilon;
+
+      //      T r=sqrt(x_new*x_new*(1-_epsilon)*(1-_epsilon)
+      //	       + y_new*y_new);
+      
+      //r/=(1-_epsilon);
+      T r_r=x_new*x_new/exp(_epsilon)+y_new*y_new/exp(-_epsilon);
+      
+      
+      return bkg+ampl*pow(1+r_r/r0/r0,-3*beta+static_cast<T>(.5));
+    }
+
+
+  };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/bl1d.hpp b/models/bl1d.hpp
new file mode 100644
index 0000000..bd0ce7d
--- /dev/null
+++ b/models/bl1d.hpp
@@ -0,0 +1,56 @@
+#ifndef BROKEN_LINE_MODEL_H_
+#define BROKEN_LINE_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+
+namespace opt_utilities
+{
+  template <typename T>
+  class bl1d
+    :public model<T,T,std::vector<T>,std::string>
+  {
+  private:
+    model<T,T,std::vector<T> >* do_clone()const
+    {
+      return new bl1d<T>(*this);
+    }
+  public:
+    bl1d()
+    {
+      this->push_param_info(param_info<std::vector<T> >("break point y value",1));
+      this->push_param_info(param_info<std::vector<T> >("break point x value",1));
+      this->push_param_info(param_info<std::vector<T> >("slop 1",1));
+      this->push_param_info(param_info<std::vector<T> >("slop 2",1));
+    }
+
+  public:
+    T do_eval(const T& x,const std::vector<T>& param)
+    {
+      T x_b=get_element(param,0);
+      T f_b=get_element(param,1);
+      T k1=get_element(param,2);
+      T k2=get_element(param,3);
+      if(x<x_b)
+	{
+	  return k1*(x-x_b)+f_b;
+	}
+      else
+	{
+	  return k2*(x-x_b)+f_b;
+	}
+    }
+
+  private:
+    std::string do_to_string()const
+    {
+      return "broken linear model\n"
+	"y=k1*(x-x_b)+y_b for x<x_b\n"
+	"y=k2*(x-x_b)+y_b otherwise\n";
+    }
+  };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/bpl1d.hpp b/models/bpl1d.hpp
new file mode 100644
index 0000000..faf1336
--- /dev/null
+++ b/models/bpl1d.hpp
@@ -0,0 +1,56 @@
+#ifndef BROKEN_POWER_LAW_MODEL_H_
+#define BROKEN_POWER_LAW_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+
+namespace opt_utilities
+{
+  template <typename T>
+  class bpl1d
+    :public model<T,T,std::vector<T>,std::string>
+  {
+  private:
+    model<T,T,std::vector<T> >* do_clone()const
+    {
+      return new bpl1d<T>(*this);
+    }
+  public:
+    bpl1d()
+    {
+      this->push_param_info(param_info<std::vector<T> >("bpx",1));
+      this->push_param_info(param_info<std::vector<T> >("bpy",1));
+      this->push_param_info(param_info<std::vector<T> >("gamma1",1));
+      this->push_param_info(param_info<std::vector<T> >("gamma2",1));
+    }
+
+    T do_eval(const T& x,const std::vector<T>& param)
+    {
+      T x_b=get_element(param,0);
+      T f_b=get_element(param,1);
+      T gamma1=get_element(param,2);
+      T gamma2=get_element(param,3);
+      if(x<x_b)
+	{
+	  return f_b*pow(x,gamma1)/pow(x_b,gamma1);
+	}
+      else
+	{
+	  return f_b*pow(x,gamma2)/pow(x_b,gamma2);
+	}
+    }
+
+
+  private:
+    std::string do_to_string()const
+    {
+      return "broken power law\n"
+	"y=y_b*(x/x_b)^gamma1 for x<x_b\n"
+	"y=y_b*(x/x_b)^gamma2 otherwise\n";
+    }
+  };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/bremss.hpp b/models/bremss.hpp
new file mode 100644
index 0000000..9956a3e
--- /dev/null
+++ b/models/bremss.hpp
@@ -0,0 +1,44 @@
+#ifndef BREMSS_MODEL_H_
+#define BREMSS_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+
+namespace opt_utilities
+{
+  template <typename T>
+  class bremss
+    :public model<T,T,std::vector<T>,std::string>
+  {
+  private:
+    model<T,T,std::vector<T> >* do_clone()const
+    {
+      return new bremss<T>(*this);
+    }
+  public:
+    bremss()
+    {
+      this->push_param_info(param_info<std::vector<T> >("norm",1));
+      this->push_param_info(param_info<std::vector<T> >("kT",1));
+    }
+
+    T do_eval(const T& x,const std::vector<T>& param)
+    {
+      T norm=get_element(param,0);
+      T kT=get_element(param,1);
+      
+      return norm*sqrt(kT)*exp(-x/kT);
+    }
+
+  private:
+    std::string do_to_string()const
+    {
+      return "Simplified bremss model\n"
+	"flux=norm*kT^0.5*e^{-E/kT}\n";
+    }
+  };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/constant.hpp b/models/constant.hpp
new file mode 100644
index 0000000..8d06279
--- /dev/null
+++ b/models/constant.hpp
@@ -0,0 +1,42 @@
+#ifndef CONSTANT_MODEL_H_
+#define CONSTANT_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+
+namespace opt_utilities
+{
+  template <typename T>
+  class constant
+    :public model<T,T,std::vector<T>,std::string>
+  {
+  private:
+    model<T,T,std::vector<T> >* do_clone()const
+    {
+      return new constant<T>(*this);
+    }
+  public:
+    constant()
+    {
+      this->push_param_info(param_info<std::vector<T> >("c",1));
+    }
+
+  public:
+    T do_eval(const T& x,const std::vector<T>& param)
+    {
+      //return x*param[0]+param[1];
+      return get_element(param,0);
+    }
+
+  private:
+    std::string do_to_string()const
+    {
+      return "Constant\n"
+	"y=C\n";
+    }
+  };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/dbeta1d.hpp b/models/dbeta1d.hpp
new file mode 100644
index 0000000..71a1f52
--- /dev/null
+++ b/models/dbeta1d.hpp
@@ -0,0 +1,62 @@
+#ifndef DBETA_MODEL_H_
+#define DBETA_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+#include <iostream>
+namespace opt_utilities
+{
+  template <typename T>
+  class dbeta1d
+    :public model<T,T,std::vector<T>,std::string>
+  {
+  private:
+    model<T,T,std::vector<T> >* do_clone()const
+    {
+      return new dbeta1d<T>(*this);
+    }
+  public:
+    dbeta1d()
+    {
+      this->push_param_info(param_info<std::vector<T> >("S01",1));
+      this->push_param_info(param_info<std::vector<T> >("rc1",30));
+      this->push_param_info(param_info<std::vector<T> >("beta1",.6));
+
+      this->push_param_info(param_info<std::vector<T> >("S02",.5));
+      this->push_param_info(param_info<std::vector<T> >("rc2",20));
+      this->push_param_info(param_info<std::vector<T> >("beta2",.4));
+
+
+      this->push_param_info(param_info<std::vector<T> >("bkg",0));
+    }
+
+
+    T do_eval(const T& x,const std::vector<T>& param)
+    {
+      T S01=get_element(param,0);
+      T r_c1=get_element(param,1);
+      T beta1=get_element(param,2);
+
+      T S02=get_element(param,3);
+      T r_c2=get_element(param,4);
+      T beta2=get_element(param,5);
+
+
+      T bkg=get_element(param,6);
+
+      return bkg+S01*pow(1+(x*x)/(r_c1*r_c1),-3*beta1+static_cast<T>(.5))
+	+S02*pow(1+(x*x)/(r_c2*r_c2),-3*beta2+static_cast<T>(.5));
+    }
+
+  private:
+    std::string do_to_string()const
+    {
+      return "double 1d beta model\n"
+	"S=beta(S01,beta1,rc1)+beta(S02,beta2,rc2)\n";
+    }
+  };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/dbeta2d.hpp b/models/dbeta2d.hpp
new file mode 100644
index 0000000..6811270
--- /dev/null
+++ b/models/dbeta2d.hpp
@@ -0,0 +1,92 @@
+#ifndef DBETA_MODEL2d_H_
+#define DBETA_MODEL2d_H_
+#include <core/fitter.hpp>
+#include <cmath>
+#include <cassert>
+#include "vecn.hpp"
+
+namespace opt_utilities
+{
+  
+  template <typename T>
+  class dbeta2d
+    :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 dbeta2d<T>(*this);
+    }
+    
+    
+    
+  public:
+    dbeta2d()
+    {
+      push_param_info(param_info<std::vector<T> >("S01",1));
+      push_param_info(param_info<std::vector<T> >("rc11",50));
+      push_param_info(param_info<std::vector<T> >("rc21",50));
+      push_param_info(param_info<std::vector<T> >("rho1",0));
+      push_param_info(param_info<std::vector<T> >("x01",200));
+      push_param_info(param_info<std::vector<T> >("y01",200));
+      push_param_info(param_info<std::vector<T> >("beta1",2./3.));
+      
+      push_param_info(param_info<std::vector<T> >("S02",1));
+      push_param_info(param_info<std::vector<T> >("rc12",60));
+      push_param_info(param_info<std::vector<T> >("rc22",60));
+      push_param_info(param_info<std::vector<T> >("rho2",0));
+      push_param_info(param_info<std::vector<T> >("x02",200));
+      push_param_info(param_info<std::vector<T> >("y02",200));
+      push_param_info(param_info<std::vector<T> >("beta2",2./2.5));
+      push_param_info(param_info<std::vector<T> >("bkg",0));
+    }
+    
+    T do_eval(const vecn<T,2>& xy,const std::vector<T>& param)
+    {
+      T S01=get_element(param,0);
+      T rc11=get_element(param,1);
+      T rc21=get_element(param,2);
+      T rho1=get_element(param,3);
+      T x01=get_element(param,4);
+      T y01=get_element(param,5);
+      T beta1=get_element(param,6);
+      T S02=get_element(param,7);
+      T rc12=get_element(param,8);
+      T rc22=get_element(param,9);
+      T rho2=get_element(param,10);
+      T x02=get_element(param,11);
+      T y02=get_element(param,12);
+      T beta2=get_element(param,13);
+      T bkg=get_element(param,14);
+      
+      
+      T x=xy[0];
+      T y=xy[1];
+
+      rho1=rho1>1?1:rho1;
+      rho1=rho1<-1?-1:rho1;
+      rho2=rho2>1?1:rho2;
+      rho2=rho2<-1?-1:rho2;
+      
+      T r1=(x-x01)*(x-x01)/(rc11*rc11)+(y-y01)*(y-y01)/(rc21*rc21)
+	-2*rho1*(x-x01)*(y-y01)/(rc11*rc21);
+
+      T r2=(x-x02)*(x-x02)/(rc12*rc12)+(y-y02)*(y-y02)/(rc22*rc22)
+	-2*rho2*(x-x02)*(y-y02)/(rc12*rc22);
+      //      r1=r1<0?0:r1;
+      //r2=r2<0?0:r2;
+      assert(r1>=0);
+      assert(r2>=0);
+      
+
+      return bkg+S01*pow(1+r1,-3*beta1+static_cast<T>(.5))+
+	S02*pow(1+r2,-3*beta2+static_cast<T>(.5));
+    }
+  };
+};
+
+
+
+#endif
+//EOF
+
diff --git a/models/dbeta2d2.hpp b/models/dbeta2d2.hpp
new file mode 100644
index 0000000..d968dbc
--- /dev/null
+++ b/models/dbeta2d2.hpp
@@ -0,0 +1,100 @@
+#ifndef DDBETA_OOOMODEL2d2_H_
+#define DDBETA_OOOMODEL2d2_H_
+#include <core/fitter.hpp>
+#include <cmath>
+#include <cassert>
+#include "vecn.hpp"
+
+
+namespace opt_utilities
+{
+  
+  template <typename T>
+  class dbeta2d2
+    :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 dbeta2d2<T>(*this);
+    }
+  public:
+    dbeta2d2()
+    {
+      this->push_param_info(param_info<std::vector<T> >("ampl1",1));//0
+      this->push_param_info(param_info<std::vector<T> >("r01",1));//1
+      this->push_param_info(param_info<std::vector<T> >("x01",100));//2
+      this->push_param_info(param_info<std::vector<T> >("y01",100));//3      
+      this->push_param_info(param_info<std::vector<T> >("theta1",0));//4      
+      this->push_param_info(param_info<std::vector<T> >("beta1",2./3.));//5
+      this->push_param_info(param_info<std::vector<T> >("epsilon1",0));//6
+
+
+      this->push_param_info(param_info<std::vector<T> >("ampl2",1));//7
+      this->push_param_info(param_info<std::vector<T> >("r02",1));//8
+      this->push_param_info(param_info<std::vector<T> >("x02",100));//9
+      this->push_param_info(param_info<std::vector<T> >("y02",100));//10
+      this->push_param_info(param_info<std::vector<T> >("theta2",0));//11
+      this->push_param_info(param_info<std::vector<T> >("beta2",2./4.));//12
+      this->push_param_info(param_info<std::vector<T> >("epsilon2",0));//13
+
+      this->push_param_info(param_info<std::vector<T> >("bkg",0));//14
+    }
+
+
+    T do_eval(const vecn<T,2>& xy,const std::vector<T>& param)
+    {
+      T x=xy[0];
+      T y=xy[1];
+
+      T ampl1=get_element(param,0);
+      T r01=get_element(param,1);
+      T x01=get_element(param,2);
+      T y01=get_element(param,3);
+
+      T theta1=get_element(param,4);
+
+      T beta1=get_element(param,5);
+      T epsilon1=get_element(param,6);
+
+      T ampl2=get_element(param,7);      
+      T r02=get_element(param,8);
+      T x02=get_element(param,9);
+      T y02=get_element(param,10);
+   
+      T theta2=get_element(param,11);
+
+      T beta2=get_element(param,12);
+
+      T epsilon2=get_element(param,13);
+      T bkg=get_element(param,14);
+
+      T x_new1=(x-x01)*cos(theta1)+(y-y01)*sin(theta1);
+      T y_new1=(y-y01)*cos(theta1)-(x-x01)*sin(theta1);
+      
+      //T r1=sqrt(x_new1*x_new1*(1-epsilon1)*(1-epsilon1)
+      //       + y_new1*y_new1);
+      //r1/=(1-epsilon1);
+
+      T r1_r1=x_new1*x_new1/exp(epsilon1/30)+y_new1*y_new1/exp(-epsilon1/30);
+
+      T x_new2=(x-x02)*cos(theta2)+(y-y02)*sin(theta2);
+      T y_new2=(y-y02)*cos(theta2)-(x-x02)*sin(theta2);
+      
+      T r2_r2=x_new2*x_new2/exp(epsilon2/30)+y_new2*y_new2/exp(-epsilon2/30);
+      //   T r2=sqrt(x_new2*x_new2*(1-epsilon2)*(1-epsilon2)
+      //	       + y_new2*y_new2);
+      //r2/=(1-epsilon2);
+      
+      
+      
+      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));
+    }
+  };
+};
+
+
+
+#endif
+//EOF
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
diff --git a/models/dl_model.hpp b/models/dl_model.hpp
new file mode 100644
index 0000000..2fcf227
--- /dev/null
+++ b/models/dl_model.hpp
@@ -0,0 +1,178 @@
+#ifdef __linux__
+
+#ifndef DL_MODEL_H_
+#define DL_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+#include <iostream>
+#include <string>
+#include <sstream>
+#include <dlfcn.h>
+
+namespace opt_utilities
+{
+  template <typename T>
+  class dl_model
+    :public model<T,T,std::vector<T>,std::string>
+  {
+  private:
+    T (*calc_model)(T x,const T* p);
+    int nparams;
+    mutable void* handle;
+  private:
+    model<T,T,std::vector<T> >* do_clone()const
+    {
+      dl_model<T>* result=new dl_model<T>(*this);
+      this->handle=NULL;
+      return result;
+    }
+
+    //  public:
+  public:
+    dl_model()
+      :handle(NULL)
+    {}
+
+
+  public:
+    dl_model(const char* file_name)
+      :handle(NULL)
+    {
+      
+      handle=dlopen(file_name,RTLD_LAZY);
+      
+      if(!handle)
+	{
+	  throw opt_exception("faild loading object");
+	}
+      
+      calc_model=(T (*)(T,const T*))dlsym(handle,"calc_model");
+      
+      if(!calc_model)
+	{
+	  throw opt_exception("symble undefined");
+	}
+
+      const char* (*get_param_name)(int)
+	=(const char* (*)(int))dlsym(handle,"get_param_name");
+      
+      if(!get_param_name)
+	{
+	  throw opt_exception("symble undefined");
+	}
+
+      int (*get_num_params)()
+	=(int (*)())dlsym(handle,"get_num_params");
+      
+      if(!get_num_params)
+	{
+	  throw opt_exception("symble undefined");
+      if(!get_num_params)
+	{
+	  throw opt_exception("symble undefined");
+	}	}
+
+      T (*get_default_value)(int)
+	=(T (*)(int))dlsym(handle,"get_default_value");
+      
+      if(!get_default_value)
+	{
+	  throw opt_exception("symble undefined");
+	}
+
+      nparams=get_num_params();
+      
+      for(int i=0;i!=nparams;++i)
+	{
+	  this->push_param_info(param_info<std::vector<T> >(get_param_name(i),
+				get_default_value(i)));
+	}
+    }
+
+    ~dl_model()
+    {
+      if(handle)
+	{
+	  dlclose(handle);
+	}
+    }
+
+    void bind(const char* file_name)
+    {
+      if(handle)
+	{
+	  dlclose(handle);
+	}
+      this->clear_param_info();
+      handle=dlopen(file_name,RTLD_LAZY);
+
+      if(!handle)
+	{
+	  throw opt_exception("faild loading object");
+	}
+
+      calc_model=(T (*)(T,const T*))dlsym(handle,"calc_model");
+      
+      if(!calc_model)
+	{
+	  throw opt_exception("symble undefined");
+	}
+
+      const char* (*get_param_name)(int)
+	=(const char* (*)(int))dlsym(handle,"get_param_name");
+      
+      if(!get_param_name)
+	{
+	  throw opt_exception("symble undefined");
+	}
+
+
+      int (*get_num_params)()
+	=(int (*)())dlsym(handle,"get_num_params");
+      
+      if(!get_num_params)
+	{
+	  throw opt_exception("symble undefined");
+	}      
+
+
+      T (*get_default_value)(int)
+	=(T (*)(int))dlsym(handle,"get_default_value");
+     
+
+      if(!get_default_value)
+	{
+	  throw opt_exception("symble undefined");
+	}
+
+      
+      nparams=get_num_params();
+      for(int i=0;i!=nparams;++i)
+	{
+	  this->push_param_info(param_info<std::vector<T> >(get_param_name(i),
+				get_default_value(i)));
+	}
+    }
+
+    T do_eval(const T& x,const std::vector<T>& param)
+    {
+      if(handle==NULL)
+	{
+	  throw opt_exception("dl object unloaded");
+	}
+      return calc_model(x,&get_element(param,0));
+    }
+
+    std::string do_to_string()const
+    {
+      return "Dynamical load model\n"
+	"Should be loaded from an shared object file\n";
+    }
+  };
+};
+
+
+
+#endif
+#endif
+//EOF
diff --git a/models/dlmodel_template.c b/models/dlmodel_template.c
new file mode 100644
index 0000000..fd7aadb
--- /dev/null
+++ b/models/dlmodel_template.c
@@ -0,0 +1,36 @@
+#include <math.h>
+
+char p1name[2]="k";
+char p2name[2]="b";
+
+int main(int argc,char* argv[])
+{}
+
+
+int get_num_params()
+{
+  return 2;
+}
+
+const char* get_param_name(int n)
+{
+  if(n==0)
+    {
+      return p1name;
+    }
+  return p2name;
+}
+
+double get_default_value(int n)
+{
+  if(n==0)
+    return 1;
+  return 0;
+}
+
+double calc_model(double x,double p[])
+{
+  return p[0]*x+p[1];
+}
+
+
diff --git a/models/func_model.hpp b/models/func_model.hpp
new file mode 100644
index 0000000..4eda8d5
--- /dev/null
+++ b/models/func_model.hpp
@@ -0,0 +1,57 @@
+#ifndef FUNC_MODEL_H_
+#define FUNC_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+#include <iostream>
+#include <string>
+#include <sstream>
+namespace opt_utilities
+{
+  template <typename T>
+  class func_model
+    :public model<T,T,std::vector<T>,std::string>
+  {
+  private:
+    T (*func)(T x,const T* const& p);
+    int nparams;
+  private:
+    model<T,T,std::vector<T> >* do_clone()const
+    {
+      return new func_model<T>(*this);
+    }
+
+    //  public:
+  private:
+    func_model()
+    {}
+
+
+  public:
+    func_model(T (*_func)(T x,const T* const& p),int n)
+      :func(_func),nparams(n)
+    {
+      for(int i=0;i!=n;++i)
+	{
+	  std::ostringstream oss;
+	  oss<<i;
+	  this->push_param_info(param_info<std::vector<T> >(oss.str(),0));
+	}
+    }
+
+
+    T do_eval(const T& x,const std::vector<T>& param)
+    {
+      return func(x,&get_element(param,0));
+    }
+  private:
+    std::string do_to_string()const
+    {
+      return "Wrapper for necked C function\n";
+    }
+  };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/gauss1d.hpp b/models/gauss1d.hpp
new file mode 100644
index 0000000..958cec9
--- /dev/null
+++ b/models/gauss1d.hpp
@@ -0,0 +1,47 @@
+#ifndef GAUSS_MODEL_H_
+#define GAUSS_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+
+namespace opt_utilities
+{
+  template <typename T>
+  class gauss1d
+    :public model<T,T,std::vector<T>,std::string>
+  {
+  private:
+    model<T,T,std::vector<T> >* do_clone()const
+    {
+      return new gauss1d<T>(*this);
+    }
+  public:
+    gauss1d()
+    {
+      this->push_param_info(param_info<std::vector<T> >("N",1));
+      this->push_param_info(param_info<std::vector<T> >("x0",0));
+      this->push_param_info(param_info<std::vector<T> >("sigma",1));
+    }
+
+  public:
+    T do_eval(const T& x,const std::vector<T>& param)
+    {
+      T N=get_element(param,0);
+      T x0=get_element(param,1);
+      T sigma=get_element(param,2);
+      T y=(x-x0)/sigma;
+      return N*exp(-y*y);
+    }
+
+  private:
+    std::string do_to_string()const
+    {
+      return "Gaussian model\n"
+	"y=N*exp(-(x-x0)^2/sigma^2\n";
+    }
+  };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/lin1d.hpp b/models/lin1d.hpp
new file mode 100644
index 0000000..b0ddd31
--- /dev/null
+++ b/models/lin1d.hpp
@@ -0,0 +1,42 @@
+#ifndef LINEAR_MODEL_H_
+#define LINEAR_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+
+namespace opt_utilities
+{
+  template <typename T>
+  class lin1d
+    :public model<T,T,std::vector<T>,std::string>
+  {
+  private:
+    model<T,T,std::vector<T> >* do_clone()const
+    {
+      return new lin1d<T>(*this);
+    }
+  public:
+    lin1d()
+    {
+      this->push_param_info(param_info<std::vector<T> >("k",1));
+      this->push_param_info(param_info<std::vector<T> >("b",0));
+    }
+
+  public:
+    T do_eval(const T& x,const std::vector<T>& param)
+    {
+      return x*get_element(param,0)+get_element(param,1);
+    }
+
+  private:
+    std::string do_to_string()const
+    {
+      return "linear model\n"
+	"y=k*x+b\n";
+    }
+  };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/models.cc b/models/models.cc
new file mode 100644
index 0000000..7481b88
--- /dev/null
+++ b/models/models.cc
@@ -0,0 +1,176 @@
+#include "models.hpp"
+#include <core/opt_exception.hpp>
+#include "gauss1d.hpp"
+#include "bl1d.hpp"
+#include "nfw1d.hpp"
+#include "bpl1d.hpp"
+#include "beta1d.hpp"
+#include "nbeta1d.hpp"
+#include "dbeta1d.hpp"
+#include "lin1d.hpp"
+#include "pl1d.hpp"
+#include "poly1d.hpp"
+#include "bremss.hpp"
+#include "beta2d2.hpp"
+#include "beta2d.hpp"
+#include "dbeta2d2.hpp"
+#include "dbeta2d3.hpp"
+#include "dbeta2d.hpp"
+#include "dl_model.hpp"
+#include <iostream>
+using namespace std;
+
+
+namespace opt_utilities
+{
+  strmodel1d strm1d;
+  std::map<std::string,model<double,double,std::vector<double>,std::string>* > model_map;
+  std::map<std::string,model<double,vecn<double,2>,std::vector<double>,std::string>* > model2d_map;
+  std::list<std::string> get_model_name_list()
+  {
+    std::list<std::string> result;
+    for(std::map<std::string,model<double,double,std::vector<double>,std::string>* >::iterator i=model_map.begin();
+	i!=model_map.end();++i)
+      {
+	result.push_back(i->first);
+      }
+	return result;
+  }
+
+  std::list<std::string> get_model2d_name_list()
+  {
+    std::list<std::string> result;
+    for(map<std::string,model<double,vecn<double,2>,std::vector<double>,std::string>* > ::iterator i=model2d_map.begin();
+	i!=model2d_map.end();++i)
+      {
+	result.push_back(i->first);
+      }
+    return result;
+  }
+
+model<double,double,std::vector<double>,std::string>& get_1dmodel_by_name(const char* name)
+  {
+    std::map<std::string,model<double,double,std::vector<double>,std::string >* >::iterator iter;
+    iter=model_map.find(name);
+    if(iter==model_map.end()||iter->second==0)
+      {
+	throw opt_exception("model does not exist");
+      }
+    return *(iter->second);
+  }
+
+  strmodel1d& get_strm1d()
+  {
+    return strm1d;
+  }
+
+model<double,vecn<double,2>,std::vector<double>,std::string>& get_2dmodel_by_name(const char* name)
+  {
+    std::map<std::string,model<double,vecn<double,2>,std::vector<double>,std::string>* >::iterator iter;
+    iter=model2d_map.find(name);
+    if(iter==model2d_map.end()||iter->second==0)
+      {
+	throw opt_exception("model does not exist");
+      }
+    return *(iter->second);
+  }
+
+  int get_n_1dmodels()
+  {
+    return model_map.size();
+  }
+
+  int get_n_2dmodels()
+  {
+    return model2d_map.size();
+  }
+
+  class model_map_keeper_class
+  {
+  private:
+    void init_model_map()
+    {
+
+      //#define DECL_POLY(n) model_map["poly1d##n"]=new poly1d<double,n>;
+
+
+      model_map["lin1d"]=new lin1d<double>;
+      model_map["pl1d"]=new pl1d<double>;
+      model_map["bl1d"]=new bl1d<double>;
+      model_map["bpl1d"]=new bpl1d<double>;
+      model_map["beta1d"]=new beta1d<double>;  
+      model_map["bremss"]=new bremss<double>;
+      model_map["nbeta1d"]=new nbeta1d<double>;  
+      model_map["2beta1d"]=new dbeta1d<double>;
+      model_map["nfw"]=new nfw1d<double>;
+      model_map["gauss1d"]=new gauss1d<double>;  
+      model_map["poly1d2"]=new poly1d<double,2>;
+      model_map["poly1d3"]=new poly1d<double,3>;
+      model_map["poly1d4"]=new poly1d<double,4>;
+      model_map["poly1d5"]=new poly1d<double,5>;
+#ifdef __linux__
+      model_map["dlmodel"]=new dl_model<double>;
+#endif
+      //DECL_POLY(7)
+    }
+    
+    void release_model_map()
+    {
+      for(std::map<std::string,model<double,double,std::vector<double>,std::string>* >::iterator i=model_map.begin();
+	  i!=model_map.end();++i)
+	{
+	  delete i->second;
+	}
+    }
+  public:
+    model_map_keeper_class()
+    {
+      init_model_map();
+	  std::cerr<<"1d models Initialized"<<std::endl;
+    }
+    
+    ~model_map_keeper_class()
+    {
+      release_model_map();
+	  std::cerr<<"1d models Released"<<std::endl;
+    }
+
+  }model_map_keeper;
+
+  class model2d_map_keeper_class
+  {
+  private:
+    void init_model_map()
+    {
+      model2d_map["beta2d"]=new beta2d<double>;
+      model2d_map["beta2d2"]=new beta2d2<double>;
+      model2d_map["dbeta2d"]=new dbeta2d<double>;
+      model2d_map["dbeta2d2"]=new dbeta2d2<double>;
+      model2d_map["dbeta2d3"]=new dbeta2d3<double>;  
+    }
+    
+    void release_model_map()
+    {
+      for(std::map<std::string,model<double,vecn<double,2>,std::vector<double>,std::string>* >::iterator i=model2d_map.begin();
+	  i!=model2d_map.end();++i)
+	{
+	  delete i->second;
+	}
+    }
+  public:
+    model2d_map_keeper_class()
+    {
+      init_model_map();
+      std::cerr<<"2d models Initialized"<<std::endl;
+    }
+    
+    ~model2d_map_keeper_class()
+    {
+      release_model_map();
+      std::cerr<<"2d models Released"<<std::endl;
+    }
+
+  }model2d_map_keeper;
+
+
+};
diff --git a/models/models.hpp b/models/models.hpp
new file mode 100644
index 0000000..79201e8
--- /dev/null
+++ b/models/models.hpp
@@ -0,0 +1,36 @@
+#ifndef MODELS_HPP
+#define MODELS_HPP
+
+#include <core/fitter.hpp>
+#include <map>
+#include <string>
+#include <list>
+#include "vecn.hpp"
+#include "strmodel1d.hpp"
+
+
+
+namespace opt_utilities
+{
+  extern std::map<std::string,model<double,double,std::vector<double>,std::string>* > model_map;
+  extern std::map<std::string,model<double,vecn<double,2>,std::vector<double>,std::string >* > model2d_map;
+
+  extern strmodel1d strm1d;
+  extern std::list<std::string> get_model_name_list();
+  extern int get_n_1dmodels();
+  //  extern void init_model_map();
+  //  extern void release_model_map();
+  
+  extern std::list<std::string> get_model2d_name_list();
+  //  extern void init_model2d_map();
+  //  extern void release_model2d_map();
+  extern int get_n_2dmodels();
+  
+  extern model<double,double,std::vector<double>,std::string >& get_1dmodel_by_name(const char*);
+  extern  model<double,vecn<double,2>,std::vector<double>,std::string >& get_2dmodel_by_name(const char*);
+
+  extern strmodel1d& get_strm1d();
+}
+
+
+#endif
diff --git a/models/mul_model.hpp b/models/mul_model.hpp
new file mode 100644
index 0000000..0e64c1b
--- /dev/null
+++ b/models/mul_model.hpp
@@ -0,0 +1,171 @@
+#ifndef MUL_MODEL_H_
+#define MUL_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+
+namespace opt_utilities
+{
+  template <typename Ty,typename Tx,typename Tp,typename Tstr>
+  class mul_model
+    :public model<Ty,Tx,Tp,Tstr>
+  {
+  private:
+    model<Ty,Tx,Tp,Tstr>* do_clone()const
+    {
+      return new mul_model<Ty,Tx,Tp,Tstr>(*this);
+    }
+  private:
+    mul_model()
+    {
+    }
+
+  private:
+    model<Ty,Tx,Tp,Tstr>* pm1;
+    model<Ty,Tx,Tp,Tstr>* pm2;
+
+  public:
+    mul_model(const model<Ty,Tx,Tp,Tstr>& m1,
+		 const model<Ty,Tx,Tp,Tstr>& m2)
+      :pm1(m1.clone()),pm2(m2.clone())
+    {
+      int np1=m1.get_num_params();
+      int np2=m2.get_num_params();
+      for(int i=0;i<np1;++i)
+	{
+	  param_info<Tp,Tstr> p(m1.get_param_info(i));
+	  param_info<Tp,Tstr> p1(p.get_name()+"1",p.get_default_value());
+	  this->push_param_info(p1);
+	}
+      for(int i=0;i<np2;++i)
+	{
+	  param_info<Tp,Tstr> p(m2.get_param_info(i));
+	  param_info<Tp,Tstr> p2(p.get_name()+"2",p.get_default_value());
+	  this->push_param_info(p2);
+	}
+    }
+
+    mul_model(const mul_model& rhs)
+      :pm1(NULL),pm2(NULL)
+    {
+      int np1(0),np2(0);
+      if(rhs.pm1)
+	{
+	  pm1=rhs.pm1->clone();
+	  np1=rhs.pm1->get_num_params();
+	  for(int i=0;i<np1;++i)
+	    {
+	      param_info<Tp,Tstr> p(rhs.pm1->get_param_info(i));
+	      param_info<Tp,Tstr> p1(p.get_name()+"1",p.get_default_value());
+	      this->push_param_info(p1);
+	    }
+	}
+      if(rhs.pm2)
+	{
+	  pm2=rhs.pm2->clone();
+	  np2=rhs.pm2->get_num_params();
+	  for(int i=0;i<np2;++i)
+	    {
+	      param_info<Tp,Tstr> p(rhs.pm2->get_param_info(i));
+	      param_info<Tp,Tstr> p2(p.get_name()+"2",p.get_default_value());
+	      this->push_param_info(p2);
+	    }
+	}
+    }
+
+    mul_model& operator=(const mul_model& rhs)
+    {
+      if(this==&rhs)
+	{
+	  return *this;
+	}
+      if(!pm1)
+	{
+	  //delete pm1;
+	  pm1->destroy();
+	}
+      if(!pm2)
+	{
+	  //delete pm2;
+	  pm2->destroy();
+	}
+      int np1(0),np2(0);
+      if(rhs.pm1)
+	{
+	  pm1=rhs.pm1->clone();
+	  np1=rhs.pm1->get_num_params();
+	  for(int i=0;i<np1;++i)
+	    {
+	      param_info<Tp,Tstr> p(rhs.pm1->get_param_info(i));
+	      param_info<Tp,Tstr> p1(p.get_name()+"1",p.get_default_value());
+	      this->push_param_info(p1);
+	    }
+	}
+      if(rhs.pm2)
+	{
+	  pm2=rhs.pm2->clone();
+	  np2=rhs.pm2->get_num_params();
+	  for(int i=0;i<np2;++i)
+	    {
+	      param_info<Tp,Tstr> p(rhs.pm2->get_param_info(i));
+	      param_info<Tp,Tstr> p2(p.get_name()+"2",p.get_default_value());
+	      this->push_param_info(p2);
+	    }
+	}
+      return *this;
+    }
+
+
+    ~mul_model()
+    {
+      if(!pm1)
+	{
+	  //delete pm1;
+	  pm1->destroy();
+	}
+      if(!pm2)
+	{
+	  //delete pm2;
+	  pm2->destroy();
+	}
+    }
+
+  public:
+    Ty do_eval(const Tx& x,const Tp& param)
+    {
+      if(!pm1)
+	{
+	  throw opt_exception("incomplete model!");
+	}
+      if(!pm2)
+	{
+	  throw opt_exception("incomplete model!");
+	}
+      Tp p1(pm1->get_num_params());
+      Tp p2(pm2->get_num_params());
+      int i=0;
+      int j=0;
+      for(i=0;i<pm1->get_num_params();++i,++j)
+	{
+	  set_element(p1,i,get_element(param,j));
+	}
+      for(i=0;i<pm2->get_num_params();++i,++j)
+	{
+	  set_element(p2,i,get_element(param,j));
+	}
+      return pm1->eval(x,p1)*pm2->eval(x,p2);
+    }
+  };
+
+  template <typename Ty,typename Tx,typename Tp>
+  mul_model<Ty,Tx,Tp,Tstr> operator*(const model<Ty,Tx,Tp,Tstr>& m1,
+				const model<Ty,Tx,Tp,Tstr>& m2)
+  {
+    return mul_model<Ty,Tx,Tp,Tstr>(m1,m2);
+  }
+
+};
+
+
+
+#endif
+//EOF
diff --git a/models/nbeta1d.hpp b/models/nbeta1d.hpp
new file mode 100644
index 0000000..3806fe9
--- /dev/null
+++ b/models/nbeta1d.hpp
@@ -0,0 +1,49 @@
+#ifndef NBETA_MODEL_H_
+#define NBETA_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+#include <iostream>
+namespace opt_utilities
+{
+  template <typename T>
+  class nbeta1d
+    :public model<T,T,std::vector<T> ,std::string>
+  {
+  private:
+    model<T,T,std::vector<T> >* do_clone()const
+    {
+      return new nbeta1d<T>(*this);
+    }
+  public:
+    nbeta1d()
+    {
+      this->push_param_info(param_info<std::vector<T> >("S0",1));
+      this->push_param_info(param_info<std::vector<T> >("rc",10));
+      this->push_param_info(param_info<std::vector<T> >("beta",2./3.));
+      this->push_param_info(param_info<std::vector<T> >("bkg",0));
+    }
+
+
+    T do_eval(const T& x,const std::vector<T>& param)
+    {
+      T S0=get_element(param,0);
+      T r_c=get_element(param,1);
+      T beta=get_element(param,2);
+      T bkg=get_element(param,3);
+
+      return bkg+S0*pow(1+(x*x)/(r_c*r_c),-3./2.*beta);
+    }
+    
+  private:
+    std::string do_to_string()const
+    {
+      return "density beta model\n"
+	"n=n0*(1+(r/rc)^2)^(-1.5*beta)\n";
+    }
+  };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/nfw1d.hpp b/models/nfw1d.hpp
new file mode 100644
index 0000000..9c2dfd4
--- /dev/null
+++ b/models/nfw1d.hpp
@@ -0,0 +1,45 @@
+#ifndef NFW_MODEL_H_
+#define NFW_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+#include <iostream>
+namespace opt_utilities
+{
+  template <typename T>
+  class nfw1d
+    :public model<T,T,std::vector<T>,std::string>
+  {
+  private:
+    model<T,T,std::vector<T> >* do_clone()const
+    {
+      return new nfw1d<T>(*this);
+    }
+  public:
+    nfw1d()
+    {
+      this->push_param_info(param_info<std::vector<T> >("rho0",1.));
+      this->push_param_info(param_info<std::vector<T> >("rs",1.));
+    }
+
+
+    T do_eval(const T& x,const std::vector<T>& param)
+    {
+      T rho0=get_element(param,0);
+      T rs=get_element(param,1);
+      
+
+      return rho0/(x/rs*(1+x/rs)*(1+x/rs));
+    }
+
+    std::string do_to_string()const
+    {
+      return "NFW model\n"
+	"y=rho0/(r/rs*(1+r/rs)^2\n";
+    }
+  };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/pl1d.hpp b/models/pl1d.hpp
new file mode 100644
index 0000000..50b4c50
--- /dev/null
+++ b/models/pl1d.hpp
@@ -0,0 +1,43 @@
+#ifndef POWER_LAW_MODEL_H_
+#define POWER_LAW_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+
+namespace opt_utilities
+{
+  template <typename T>
+  class pl1d
+    :public model<T,T,std::vector<T>,std::string>
+  {
+  private:
+    model<T,T,std::vector<T> >* do_clone()const
+    {
+      return new pl1d<T>(*this);
+    }
+  public:
+    pl1d()
+    {
+      this->push_param_info(param_info<std::vector<T> >("Ampl",1));
+      this->push_param_info(param_info<std::vector<T> >("gamma",1));
+    }
+
+    T do_eval(const T& x,const std::vector<T>& param)
+    {
+      T A=get_element(param,0);
+      T gamma=get_element(param,1);
+      return A*pow(x,gamma);
+    }
+
+  private:
+    std::string do_to_string()const
+    {
+      return "Simple power law model\n"
+	"y=A*x^gamma\n";
+    }
+  };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/poly1d.hpp b/models/poly1d.hpp
new file mode 100644
index 0000000..a8c994b
--- /dev/null
+++ b/models/poly1d.hpp
@@ -0,0 +1,67 @@
+#ifndef POLY_MODEL_H_
+#define POLY_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+#include <sstream>
+#include <cassert>
+
+namespace opt_utilities
+{
+  template <typename T,int n>
+  class poly1d
+    :public model<T,T,std::vector<T>,std::string>
+  {
+  private:
+    model<T,T,std::vector<T> >* do_clone()const
+    {
+      return new poly1d<T,n>(*this);
+    }
+  public:
+    poly1d()
+    {
+      assert(n>=0);
+      for(int i=0;i<=n;++i)
+	{
+	  std::ostringstream ostr;
+	  ostr<<"a"<<i;
+	  this->push_param_info(param_info<std::vector<T> >(ostr.str().c_str(),1));  
+	}
+      
+    }
+
+  public:
+    T do_eval(const T& x,const std::vector<T>& param)
+    {
+      //      return x*get_element(param,0)+get_element(param,1);
+      T result(0);
+      for(int i=0;i<=n;++i)
+	{
+	  T xn(1);
+	  for(int j=0;j<i;++j)
+	    {
+	      xn*=x;
+	    }
+	  result+=get_element(param,i)*xn;
+	}
+      return result;
+    }
+
+  private:
+    std::string do_to_string()const
+    {
+      std::ostringstream ostr;
+      ostr<<n<<"-order polynorminal model\n";
+      for(int i=0;i<n;++i)
+	{
+	  ostr<<"a"<<i<<"*x^"<<i<<"+";
+	}
+      ostr<<"a"<<n<<"*x^"<<n;
+      return ostr.str();
+    }
+  };
+};
+
+
+
+#endif
+//EOF
diff --git a/models/pow_model.hpp b/models/pow_model.hpp
new file mode 100644
index 0000000..79baa4b
--- /dev/null
+++ b/models/pow_model.hpp
@@ -0,0 +1,120 @@
+#ifndef POW_MODEL_H_
+#define POW_MODEL_H_
+#include <core/fitter.hpp>
+#include <cmath>
+
+namespace opt_utilities
+{
+  template <typename Ty,typename Tx,typename Tp,typename Tstr>
+  class pow_model
+    :public model<Ty,Tx,Tp,Tstr>
+  {
+  private:
+    model<Ty,Tx,Tp,Tstr>* do_clone()const
+    {
+      return new pow_model<Ty,Tx,Tp,Tstr>(*this);
+    }
+  private:
+    pow_model()
+    {
+    }
+
+  private:
+    model<Ty,Tx,Tp,Tstr>* pm1;
+    typename value_type_trait<Tp>::value_type idx;
+
+  public:
+    pow_model(const model<Ty,Tx,Tp,Tstr>& m1,
+	      const typename value_type_trait<Tp>::value_type& index)
+      :pm1(m1.clone()),idx(index)
+    {
+      int np1=m1.get_num_params();
+      
+      for(int i=0;i<np1;++i)
+	{
+	  param_info<Tp,Tstr> p(m1.get_param_info(i));
+	  //param_info<Tp,Tstr> p1(p.get_name(),p.get_default_value());
+	  this->push_param_info(p);
+	}
+    }
+
+    pow_model(const pow_model& rhs)
+      :pm1(NULL),idx(0)
+    {
+      int np1(0);
+      if(rhs.pm1)
+	{
+	  pm1=rhs.pm1->clone();
+	  np1=rhs.pm1->get_num_params();
+	  for(int i=0;i<np1;++i)
+	    {
+	      param_info<Tp,Tstr> p(rhs.pm1->get_param_info(i));
+	      param_info<Tp,Tstr> p1(p.get_name()+"1",p.get_default_value());
+	      this->push_param_info(p1);
+	    }
+	}
+      idx=rhs.idx;
+    }
+
+    pow_model& operator=(const pow_model& rhs)
+    {
+      if(this==&rhs)
+	{
+	  return *this;
+	}
+      if(!pm1)
+	{
+	  //delete pm1;
+	  pm1->destroy();
+	}
+
+      int np1(0);
+      if(rhs.pm1)
+	{
+	  pm1=rhs.pm1->clone();
+	  np1=rhs.pm1->get_num_params();
+	  for(int i=0;i<np1;++i)
+	    {
+	      param_info<Tp,Tstr> p(rhs.pm1->get_param_info(i));
+	      // param_info<Tp,Tstr> p1(p.get_name()+"1",p.get_default_value());
+	      this->push_param_info(p);
+	    }
+	}
+      idx=rhs.idx;
+
+      return *this;
+    }
+
+    ~pow_model()
+    {
+      if(!pm1)
+	{
+	  //delete pm1;
+	  pm1->destroy();
+	}
+    }
+
+  public:
+    Ty do_eval(const Tx& x,const Tp& param)
+    {
+      if(!pm1)
+	{
+	  throw opt_exception("incomplete model!");
+	}
+
+      return std::pow(pm1->eval(x,param),idx);
+    }
+  };
+  
+  template <typename Ty,typename Tx,typename Tp,typename Tstr>
+  pow_model<Ty,Tx,Tp,Tstr> pow(const model<Ty,Tx,Tp,Tstr>& m1,
+			  const typename value_type_trait<Tp>::value_type& idx)
+  {
+    return pow_model<Ty,Tx,Tp,Tstr>(m1,idx);
+  }
+};
+
+
+
+#endif
+//EOF
diff --git a/models/strmodel1d.cc b/models/strmodel1d.cc
new file mode 100644
index 0000000..df3cd40
--- /dev/null
+++ b/models/strmodel1d.cc
@@ -0,0 +1,82 @@
+#include "strmodel1d.hpp"
+
+using namespace mu;
+using namespace std;
+using namespace opt_utilities;
+
+strmodel1d* strmodel1d::do_clone()const
+{
+  return new strmodel1d(*this);
+}
+
+strmodel1d::strmodel1d()
+{
+  set_buildin_fun();  
+}
+
+strmodel1d::strmodel1d(const strmodel1d& rhs)
+  :expr(rhs.expr),
+   par_names(rhs.par_names),
+   var_name(rhs.var_name),
+   par_vec(rhs.par_vec)
+{
+  set_buildin_fun();
+  set_expr(expr,par_names,var_name);
+}
+
+strmodel1d& strmodel1d::operator=(const strmodel1d& rhs)
+{
+  set_buildin_fun();
+  expr=rhs.expr;
+  par_names=rhs.par_names;
+  var_name=rhs.var_name;
+  par_vec=rhs.par_vec;
+  set_expr(expr,par_names,var_name);
+  return *this;
+}
+
+void strmodel1d::set_buildin_fun()
+{
+  mp.DefineFun("sin",sin);
+  mp.DefineFun("cos",cos);
+  mp.DefineFun("exp",exp);
+  mp.DefineFun("tan",tan);
+  mp.DefineFun("log",log);
+}
+
+
+void strmodel1d::set_expr(const string& _expr,
+			  const std::vector<std::string>& _par_names,
+			  const std::string& _var_name)
+{
+  expr=_expr;
+  par_names=_par_names;
+  var_name=_var_name;
+  this->clear_param_info();
+  par_vec.resize(par_names.size());
+  mp.ClearVar();
+  //  mp.ClearFun();
+  for(int i=0;i<par_vec.size();++i)
+    {
+      mp.DefineVar(par_names[i].c_str(),&par_vec[i]);
+      this->push_param_info(param_info<std::vector<double> >(par_names[i],0));
+    }
+  mp.DefineVar(var_name.c_str(),&x);
+  mp.SetExpr(expr.c_str());
+}
+
+
+double strmodel1d::do_eval(const double& _x,const vector<double>& p)
+{
+  for(int i=0;i<par_vec.size();++i)
+    {
+      //get_element(par_vec,i)=get_element(p,i);
+      set_element(par_vec,i,get_element(p,i));
+      //      cout<<par_vec[i]<<" ";
+    }
+  x=_x;
+  //  cout<<x<<" ";
+  double result(mp.Eval());
+  //  cout<<result<<endl;
+  return result;
+}
diff --git a/models/strmodel1d.hpp b/models/strmodel1d.hpp
new file mode 100644
index 0000000..8d992a0
--- /dev/null
+++ b/models/strmodel1d.hpp
@@ -0,0 +1,39 @@
+#ifndef STRMODEL1D_HPP
+#define STRMODEL1D_HPP
+#include <core/fitter.hpp>
+#include <cmath>
+#include <sstream>
+#include <cassert>
+#include <muparser/muParser.h>
+#include <vector>
+#include <string>
+
+class strmodel1d
+  :public opt_utilities::model<double,double,std::vector<double>,std::string>
+{
+private:
+  mu::Parser mp;
+  strmodel1d* do_clone()const;
+  std::vector<double> par_vec;
+  std::vector<std::string> par_names;
+  std::string expr;
+  std::string var_name;
+  double x;
+  void set_buildin_fun();
+public:
+  double do_eval(const double& x,const std::vector<double>& p);
+  strmodel1d();
+  strmodel1d(const strmodel1d& rhs);
+  strmodel1d& operator=(const strmodel1d& rhs);
+  
+
+  void set_expr(const std::string& _expr,
+		const std::vector<std::string>& _par_names,
+		const std::string& _var_name);
+
+};
+
+
+
+#endif
+//EOF
diff --git a/models/vecn.hpp b/models/vecn.hpp
new file mode 100644
index 0000000..cac1436
--- /dev/null
+++ b/models/vecn.hpp
@@ -0,0 +1,62 @@
+#ifndef VECN_HPP
+#define VECN_HPP
+#include <iostream>
+namespace opt_utilities
+{
+  template <typename T,int n>
+  class vecn
+  {
+  public:
+    T data[n];
+  public:
+    T& operator[](int i)
+    {
+      return data[i];
+    }
+    
+    const T& operator[](int i)const
+    {
+      return data[i];
+    }
+    
+    vecn()
+    {
+      for(int i=0;i<n;++i)
+	{
+	  data[i]=0;
+	}
+    }
+    
+  };
+
+  template <typename T,int n>
+  std::istream& operator>>(std::istream& is,vecn<T,n>& p)
+  {
+    for(int i=0;i<n;++i)
+      {
+	is>>p[i];
+	//	std::cout<<i<<std::endl;
+      }
+    return is;
+  }
+
+
+  template <typename T,int n>
+  std::ostream& operator<<(std::ostream& os,const vecn<T,n>& p)
+  {
+    os<<'[';
+    for(int i=0;i<n;++i)
+      {
+	os<<p[i]<<",";
+	//	std::cout<<i<<std::endl;
+      }
+    os<<']';
+    return os;
+  }
+
+
+  
+};
+
+#endif
+//EOF