19using namespace CurveFitting;
46 vector<double> coeff(
m_nlayer * 3 + 7, 0.0);
48 for (
int i = 0; i < 7; ++i) {
61 std::vector<double> dn(
m_nlayer + 2);
62 std::vector<double> rnbn(
m_nlayer + 2);
63 std::vector<double> zbb(
m_nlayer + 2);
65 std::vector<double> cy(nData);
67 std::complex<double> c0(0.0, 0.0);
68 std::complex<double> ci(0.0, 1.0);
69 std::complex<double> cr(1.0, 0.0);
70 std::vector<std::complex<double>> rnfn(
m_nlayer + 2);
71 std::vector<std::complex<double>> pfn(
m_nlayer + 2);
72 std::vector<std::complex<double>> betan(
m_nlayer + 2);
73 std::complex<double> rnf(0.0, 0.0);
74 std::complex<double> rnf1(0.0, 0.0);
75 std::complex<double> ac1(0.0, 0.0);
76 std::complex<double> ac2(0.0, 0.0);
77 std::complex<double> ac3(0.0, 0.0);
78 std::complex<double> ac4(0.0, 0.0);
79 std::complex<double> a12t(0.0, 0.0);
80 std::complex<double> a22t(0.0, 0.0);
81 std::complex<double> btm(0.0, 0.0);
82 std::complex<double> btm1(0.0, 0.0);
83 std::complex<double> cbtm(0.0, 0.0);
84 std::complex<double> cbtm1(0.0, 0.0);
86 std::complex<double> a111(0.0, 0.0);
87 std::complex<double> a112(0.0, 0.0);
88 std::complex<double> a121(0.0, 0.0);
89 std::complex<double> a122(0.0, 0.0);
90 std::complex<double> a211(0.0, 0.0);
91 std::complex<double> a212(0.0, 0.0);
92 std::complex<double> a221(0.0, 0.0);
93 std::complex<double> a222(0.0, 0.0);
94 std::complex<double> a311(0.0, 0.0);
95 std::complex<double> a312(0.0, 0.0);
96 std::complex<double> a321(0.0, 0.0);
97 std::complex<double> a322(0.0, 0.0);
99 double theta0 = coeff[0] * M_PI / 180.0;
100 double scalefac = coeff[1];
103 zbb[0] = coeff[4] * coeff[4];
105 double bgd = coeff[5];
106 double pthet = coeff[6];
112 for (
int i = 0; i <
m_nlayer; ++i) {
113 rnbn[i + 1] = coeff[7 + i * 3];
114 dn[i + 1] = coeff[8 + i * 3];
115 zbb[i + 1] = coeff[9 + i * 3];
119 double tmax, tmin,
x, st0, ct0, ans, gauss, f;
121 std::vector<double> xnit(nit1);
124 double dthet = theta0 * pthet / 100.0;
125 dthet = dthet / 2.35;
126 double dthetr = dthet * 2.51;
127 tmax = theta0 + dthet * 3;
128 tmin = theta0 - dthet * 3;
129 double dt = (tmax - tmin) / nit;
130 for (
int i = 0; i < nit1; i++) {
131 xnit[i] = tmin + dt * i;
136 for (
size_t j = 0; j < nData; ++j) {
137 double lambda = 4 * M_PI * sin(theta0) / xValues[j];
140 double tlc = 8.0 * M_PI * M_PI / tl;
141 double con = tl / (2.0 * M_PI);
143 for (k = 0; k <
m_nlayer + 2; k++) {
144 rnfn[k] = (1.0 - con * rnbn[k]);
147 for (ii = 0; ii < nit1; ii++) {
151 pfn[0] = rnfn[0] * st0;
153 for (
int i = 1; i <
m_nlayer + 1; i++) {
154 rnf = (rnfn[i] * rnfn[i]) * cr;
155 rnf1 = (rnfn[0] * rnfn[0]) * cr;
156 pfn[i] = sqrt(rnf - (rnf1 * ct0 * ct0));
160 rnf1 = (rnfn[0] * rnfn[0]) * cr;
161 pfn[
m_nlayer + 1] = sqrt(rnf - (rnf1 * ct0 * ct0));
163 for (
int i = 1; i <
m_nlayer + 1; i++) {
164 betan[i] = 2.0 * M_PI * dn[i] * pfn[i] /
lambda;
168 if ((pfn[0] + pfn[1]) != c0) {
169 a12t = (pfn[0] - pfn[1]) / (pfn[0] + pfn[1]);
173 a112 = a12t * exp(-tlc * zbb[0] * pfn[0] * pfn[1]);
177 for (
int i = 1; i <
m_nlayer + 1; i++) {
179 btm1 = -1.0 * betan[i] * ci;
180 if (imag(btm1) == 0.0)
181 btm = betan[i] * 2.0 * ci;
183 if (real(btm1) == 0.0)
185 if (imag(btm1) == 0.0)
188 if ((pfn[i] + pfn[i + 1]) != c0)
189 a22t = (pfn[i] - pfn[i + 1]) / (pfn[i] + pfn[i + 1]);
190 if ((pfn[i] + pfn[i + 1]) == c0)
192 a22t = a22t * exp(-tlc * zbb[i] * pfn[i] * pfn[i + 1]);
197 a311 = a111 * a211 + a112 * a221;
198 a312 = a111 * a212 + a112 * a222;
199 a321 = a121 * a211 + a122 * a221;
200 a322 = a121 * a212 + a122 * a222;
215 ans = real((ac1 * ac2) / (ac3 * ac4));
218 gauss = (1.0 / dthetr) * exp(-0.5 * ((theta0 -
x) / dthet) * (theta0 -
x) / dthet);
220 cy[j] = cy[j] + f * dt;
222 out[j] = (cy[j] + bgd) * scalefac;
238 if (attName ==
"nlayer") {
241 throw std::invalid_argument(
"ReflectivityMulf: reflectivity number of "
242 "layers cannot be negative.");
246 vector<double> coeff(
m_nlayer * 3 + 7, 0.0);
247 for (
int i = 0; i < 7; ++i)
258 coeff[7 + i * 3] = 0.0;
259 coeff[8 + i * 3] = 0.0;
260 coeff[9 + i * 3] = 0.0;
263 for (
int i = 0; i <
m_nlayer; ++i) {
285 for (
int i = 0; i <
m_nlayer; ++i) {
const std::vector< double > * lambda
#define DECLARE_FUNCTION(classname)
Macro for declaring a new type of function to be used with the FunctionFactory.
Attribute is a non-fitting parameter.
int asInt() const
Returns int value if attribute is a int, throws exception otherwise.
void declareAttribute(const std::string &name, const API::IFunction::Attribute &defaultValue)
Declare a single attribute.
void storeAttributeValue(const std::string &name, const API::IFunction::Attribute &value)
Store an attribute's value.
void clearAllParameters()
Nonvirtual member which removes all declared parameters.
void declareParameter(const std::string &name, double initValue=0, const std::string &description="") override
Declare a new parameter.
double getParameter(size_t i) const override
Get i-th parameter.
ReflectivityMulf : Calculate the ReflectivityMulf from a simple layer model.
void function1D(double *out, const double *xValues, const size_t nData) const override
Function to calcualte reflectivity.
void setAttribute(const std::string &attName, const Attribute &) override
Set a value to attribute attName.
int m_nlayer
ReflectivityMulf layers.
void init() override
Function initialization. Declare function parameters in this method.
std::string to_string(const wide_integer< Bits, Signed > &n)