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
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
|
// Copyright (C) 2016 The Qt Company Ltd.
// SPDX-License-Identifier: LicenseRef-Qt-Commercial OR GPL-3.0-only
/*
* Galaxy creation code obtained from http://beltoforion.de/galaxy/galaxy_en.html
* Thanks to Ingo Berg, great work.
* Licensed under a Creative Commons Attribution 3.0 License
* http://creativecommons.org/licenses/by/3.0/
*/
#include "cumulativedistributor.h"
#include <QDebug>
#include <math.h>
CumulativeDistributor::CumulativeDistributor()
: m_pDistFun(NULL),
m_vM1(),
m_vY1(),
m_vX1(),
m_vM2(),
m_vY2(),
m_vX2()
{
}
CumulativeDistributor::~CumulativeDistributor()
{
}
void CumulativeDistributor::setupRealistic(qreal I0, qreal k, qreal a,
qreal RBulge, qreal min, qreal max,
int nSteps)
{
m_fMin = min;
m_fMax = max;
m_nSteps = nSteps;
m_I0 = I0;
m_k = k;
m_a = a;
m_RBulge = RBulge;
m_pDistFun = &CumulativeDistributor::Intensity;
// build the distribution function
buildCDF(m_nSteps);
}
void CumulativeDistributor::buildCDF(int nSteps)
{
qreal h = (m_fMax - m_fMin) / nSteps;
qreal x = 0, y = 0;
m_vX1.clear();
m_vY1.clear();
m_vX2.clear();
m_vM1.clear();
m_vM1.clear();
// Simpson rule for integration of the distribution function
m_vY1.push_back(0.0);
m_vX1.push_back(0.0);
for (int i = 0; i < nSteps; i += 2) {
x = (i + 2) * h;
y += h/3 * ((this->*m_pDistFun)(m_fMin + i*h) + 4*(this->*m_pDistFun)(m_fMin + (i+1)*h) + (this->*m_pDistFun)(m_fMin + (i+2)*h) );
m_vM1.push_back((y - m_vY1.back()) / (2*h));
m_vX1.push_back(x);
m_vY1.push_back(y);
}
// normieren
for (int i = 0; i < m_vM1.size(); i++) {
m_vY1[i] /= m_vY1.back();
m_vM1[i] /= m_vY1.back();
}
m_vY2.clear();
m_vM2.clear();
m_vY2.push_back(0.0);
qreal p = 0;
h = 1.0 / nSteps;
for (int i = 1, k = 0; i < nSteps; ++i) {
p = (qreal)i * h;
for (; m_vY1[k+1] <= p; ++k) {
}
y = m_vX1[k] + (p - m_vY1[k]) / m_vM1[k];
m_vM2.push_back( (y - m_vY2.back()) / h);
m_vX2.push_back(p);
m_vY2.push_back(y);
}
}
qreal CumulativeDistributor::valFromProp(qreal fVal)
{
if (fVal < 0.0 || fVal > 1.0)
throw std::runtime_error("out of range");
qreal h = 1.0 / m_vY2.size();
int i = int(fVal / h);
qreal remainder = fVal - i*h;
int min = qMin(m_vY2.size(), m_vM2.size());
if (i >= min)
i = min - 1;
return (m_vY2[i] + m_vM2[i] * remainder);
}
qreal CumulativeDistributor::IntensityBulge(qreal R, qreal I0, qreal k)
{
return I0 * exp(-k * pow(R, 0.25));
}
qreal CumulativeDistributor::IntensityDisc(qreal R, qreal I0, qreal a)
{
return I0 * exp(-R/a);
}
qreal CumulativeDistributor::Intensity(qreal x)
{
return (x < m_RBulge) ? IntensityBulge(x, m_I0, m_k) : IntensityDisc(x - m_RBulge, IntensityBulge(m_RBulge, m_I0, m_k), m_a);
}
|
