1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
|
/**
\file dbeta2d.hpp
\brief 2d double beta model
\author Junhua Gu
*/
#ifndef DBETA_MODEL2d_H_
#define DBETA_MODEL2d_H_
#define OPT_HEADER
#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);
}
const char* do_get_type_name()const
{
return "2d double beta model";
}
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
|
