3 files changed, 139 insertions, 3 deletions

diff --git a/dynamical_fit/dynamical_fit.cpp b/dynamical_fit/dynamical_fit.cpp
index 2ab33ce..196df0d 100644
--- a/dynamical_fit/dynamical_fit.cpp
+++ b/dynamical_fit/dynamical_fit.cpp
@@ -15,9 +15,9 @@ using namespace opt_utilities;
 
 int main(int argc,char* argv[])
 {
-  if(argc!=2)
+  if(argc!=3)
     {
-      cerr<<"Usage: "<<argv[0]<<" <config file>"<<endl;
+      cerr<<"Usage: "<<argv[0]<<" <config file> <model dump file>"<<endl;
       exit(-1);
     }
   ifstream cfg_file(argv[1]);
@@ -57,7 +57,30 @@ int main(int argc,char* argv[])
 
   fit.load_data(dl.get_data_set());
 
+  fit.set_precision(1e-5);
+  fit.fit();
+  fit.fit();
+
   vector<double> p=fit.fit();
-  cout<<p[0]<<"\t"<<p[1]<<endl;
+  for(int i=0;i<fit.get_num_params();++i)
+    {
+      cerr<<fit.get_param_info(i).get_name()<<"\t"<<fit.get_param_info(i).get_value()<<endl;
+    }
   
+  ofstream ofs(argv[2]);
+  ofs<<"no no no\n";
+  double xmin,xmax;
+  xmin=fit.get_data_set().get_data(0).get_x();
+  xmax=fit.get_data_set().get_data(0).get_x();
+  for(int i=0;i<fit.get_data_set().size();++i)
+    {
+      xmin=std::min(xmin,fit.get_data_set().get_data(i).get_x());
+      xmax=std::max(xmax,fit.get_data_set().get_data(i).get_x());
+    }
+  
+  for(double x=xmin;x<xmax;x+=(xmax-xmin)/1000.)
+    {
+      ofs<<x<<"\t0\t0\t"<<fit.eval_model_raw(x,fit.get_all_params())<<"\t0\t0"<<endl;
+    }
+
 }
diff --git a/dynamical_models/lin1d.cpp b/dynamical_models/lin1d.cpp
new file mode 100644
index 0000000..de43a42
--- /dev/null
+++ b/dynamical_models/lin1d.cpp
@@ -0,0 +1,49 @@
+#include <core/fitter.hpp>
+#include <cmath>
+#include <vector>
+using namespace std;
+using 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);
+  }
+  
+  const char* do_get_type_name()const
+  {
+    return "1d linear model";
+  }
+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_get_information()const
+  {
+    return "<math><mrow> <mtext>f(x;k,b)=k x+b</mtext>		\
+    </mrow>							\
+</math>";
+  }
+};
+
+static lin1d<double> _model;
+
+
+extern "C" model<double,double,vector<double>,std::string>*
+create_model_object()
+{
+  return &_model;
+}
diff --git a/dynamical_models/vt_temperature.cpp b/dynamical_models/vt_temperature.cpp
new file mode 100644
index 0000000..39b474c
--- /dev/null
+++ b/dynamical_models/vt_temperature.cpp
@@ -0,0 +1,64 @@
+#include <core/fitter.hpp>
+#include <vector>
+#include <string>
+#include <cmath>
+
+using namespace std;
+using namespace opt_utilities;
+class vk_temperature
+  :public opt_utilities::model<double,double,std::vector<double>,std::string>
+{
+private:
+  opt_utilities::model<double,double,std::vector<double>,std::string>* do_clone()const;
+public:
+  vk_temperature();
+  double do_eval(const double& r,const std::vector<double>& param);
+  
+};
+
+
+model<double,double,vector<double>,std::string>* vk_temperature::do_clone()const
+{
+  return new vk_temperature(*this);
+}
+
+vk_temperature::vk_temperature()
+{
+  this->push_param_info(param_info<std::vector<double> >("T0",2,0.1,10));
+  this->push_param_info(param_info<std::vector<double> >("rcool",1e22,0.001,10));
+  this->push_param_info(param_info<std::vector<double> >("acool",1));
+  this->push_param_info(param_info<std::vector<double> >("Tmin",1,0.1,10));
+  this->push_param_info(param_info<std::vector<double> >("rt",1e22,0.001,10));
+  this->push_param_info(param_info<std::vector<double> >("a",1,-10,10));
+  this->push_param_info(param_info<std::vector<double> >("c",1,-10,10));
+  this->push_param_info(param_info<std::vector<double> >("b",1,-10,10));
+}
+
+
+
+
+
+double vk_temperature::do_eval(const double& r,const std::vector<double>& param)
+{
+  double T0=param[0];
+  double rcool=param[1];
+  double acool=param[2];
+  double Tmin=param[3];
+  double rt=param[4];
+  double a=param[5];
+  double c=param[6];
+  double b=param[7];
+  double x=pow(r/rcool,acool);
+  double tcool=(x+Tmin/T0)/(x+1);
+  double t=pow(r/rt,-a)/pow(1+pow(r/rt,b),c/b);
+  return T0*t*tcool;
+}
+
+
+static vk_temperature _model;
+
+extern "C" model<double,double,vector<double>,std::string>*
+create_model_object()
+{
+  return &_model;
+}