42 Real start,
Real end,
Size size,
const std::pair<Real, Real>& cPoints,
43 const bool requireCPoint)
46 QL_REQUIRE(end > start,
"end must be larger than start");
48 const Real cPoint = cPoints.first;
53 "cPoint must be between start and end");
55 "density > 0 required");
57 "density must be given if cPoint is given");
59 "cPoint is required in grid but not given");
64 std::vector<Real> u, z;
65 ext::shared_ptr<Interpolation> transform;
66 const Real c1 = std::asinh((start - cPoint) / density);
67 const Real c2 = std::asinh((end - cPoint) / density);
71 if (!
close(cPoint, start) && !
close(cPoint, end)) {
72 const Real z0 = -c1 / (c2 - c1);
74 std::max(std::min(std::lround(z0 * (
size - 1)),
75 static_cast<long>(
size) - 2),
83 transform = ext::shared_ptr<Interpolation>(
87 for (
Size i = 1; i <
size - 1; ++i) {
88 const Real li = requireCPoint ? (*transform)(i*dx) : i*dx;
90 + density*std::sinh(c1*(1.0 - li) + c2*li);
94 for (
Size i = 1; i <
size - 1; ++i) {
102 for (
Size i = 0; i <
size - 1; ++i) {
109 class OdeIntegrationFct {
111 OdeIntegrationFct(
const std::vector<Real>& points,
112 const std::vector<Real>& betas,
118 return rk_(odeFct, y0, x0, x1);
127 return a/std::sqrt(
s);
134 bool equal_on_first(
const std::pair<Real, Real>& p1,
135 const std::pair<Real, Real>& p2) {
143 const std::vector<ext::tuple<Real, Real, bool> >& cPoints,
147 QL_REQUIRE(end > start,
"end must be larger than start");
149 std::vector<Real> points, betas;
150 for (
const auto& cPoint : cPoints) {
151 points.push_back(ext::get<0>(cPoint));
152 betas.push_back(
squared(ext::get<1>(cPoint) * (end - start)));
157 for (
Size i=0; i < points.size(); ++i) {
158 const Real c1 = std::asinh((start-points[i])/betas[i]);
159 const Real c2 = std::asinh((end-points[i])/betas[i]);
160 aInit+=(c2-c1)/points.size();
163 OdeIntegrationFct fct(points, betas, tol);
165 [&](
Real x) ->
Real {
return fct.solve(x, start, 0.0, 1.0) - end; },
166 tol, aInit, 0.1*aInit);
170 x[0] = 0.0;
y[0] = start;
174 y[i] = fct.solve(a,
y[i-1], x[i-1], x[i]);
178 const Real dy =
y.back() - end;
186 std::vector<std::pair<Real, Real> > w(1, std::make_pair(0.0, 0.0));
188 for (
Size i=0; i < points.size(); ++i) {
189 if (ext::get<2>(cPoints[i]) && points[i] > start && points[i] < end) {
191 const Size j = std::distance(
y.begin(),
192 std::lower_bound(
y.begin(),
y.end(), points[i]));
195 [&](
Real x) ->
Real {
return odeSolution(x,
true) - points[i]; },
198 w.emplace_back(std::min(x[
size - 2], x[j]), e);
201 w.emplace_back(1.0, 1.0);
202 std::sort(w.begin(), w.end());
203 w.erase(std::unique(w.begin(), w.end(), equal_on_first), w.end());
205 std::vector<Real> u(w.size()), z(w.size());
206 for (
Size i=0; i < w.size(); ++i) {
Runge-Kutta ODE integration.
1-D array used in linear algebra.
ext::function< T(const Real, const T)> OdeFct1d
1-D array used in linear algebra.
Concentrating1dMesher(Real start, Real end, Size size, const std::pair< Real, Real > &cPoints=(std::pair< Real, Real >(Null< Real >(), Null< Real >())), bool requireCPoint=false)
std::vector< Real > locations_
std::vector< Real > dplus_
std::vector< Real > dminus_
Linear interpolation between discrete points
template class providing a null value for a given type.
Real solve(const F &f, Real accuracy, Real guess, Real step) const
floating-point comparisons
const std::vector< Real > & betas_
const std::vector< Real > & points_
One-dimensional grid mesher concentrating around critical points.
Classes and functions for error handling.
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
std::size_t Size
size of a container
linear interpolation between discrete points
functionals and combinators not included in the STL
bool close(const Quantity &m1, const Quantity &m2, Size n)
bool close_enough(const Quantity &m1, const Quantity &m2, Size n)