Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 09:39:47 

0001 /*
0002  * Copyright Nick Thompson, 2020
0003  * Use, modification and distribution are subject to the
0004  * Boost Software License, Version 1.0. (See accompanying file
0005  * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0006  */
0007 #ifndef BOOST_MATH_INTERPOLATORS_DETAIL_SEPTIC_HERMITE_DETAIL_HPP
0008 #define BOOST_MATH_INTERPOLATORS_DETAIL_SEPTIC_HERMITE_DETAIL_HPP
0009 #include <algorithm>
0010 #include <stdexcept>
0011 #include <sstream>
0012 #include <limits>
0013 #include <cmath>
0014 
0015 namespace boost {
0016 namespace math {
0017 namespace interpolators {
0018 namespace detail {
0019 
0020 template<class RandomAccessContainer>
0021 class septic_hermite_detail {
0022 public:
0023     using Real = typename RandomAccessContainer::value_type;
0024     septic_hermite_detail(RandomAccessContainer && x, RandomAccessContainer && y, RandomAccessContainer && dydx, RandomAccessContainer && d2ydx2, RandomAccessContainer && d3ydx3) 
0025     : x_{std::move(x)}, y_{std::move(y)}, dydx_{std::move(dydx)}, d2ydx2_{std::move(d2ydx2)}, d3ydx3_{std::move(d3ydx3)}
0026     {
0027         if (x_.size() != y_.size())
0028         {
0029             throw std::domain_error("Number of abscissas must = number of ordinates.");
0030         }
0031         if (x_.size() != dydx_.size())
0032         {
0033             throw std::domain_error("Numbers of derivatives must = number of abscissas.");
0034         }
0035         if (x_.size() != d2ydx2_.size())
0036         {
0037             throw std::domain_error("Number of second derivatives must equal number of abscissas.");
0038         }
0039         if (x_.size() != d3ydx3_.size())
0040         {
0041             throw std::domain_error("Number of third derivatives must equal number of abscissas.");
0042         }
0043 
0044         if (x_.size() < 2)
0045         {
0046             throw std::domain_error("At least 2 abscissas are required.");
0047         }
0048         Real x0 = x_[0];
0049         for (decltype(x_.size()) i = 1; i < x_.size(); ++i)
0050         {
0051             Real x1 = x_[i];
0052             if (x1 <= x0)
0053             {
0054                 throw std::domain_error("Abscissas must be sorted in strictly increasing order x0 < x1 < ... < x_{n-1}");
0055             }
0056             x0 = x1;
0057         }
0058     }
0059 
0060     void push_back(Real x, Real y, Real dydx, Real d2ydx2, Real d3ydx3)
0061     {
0062         using std::abs;
0063         using std::isnan;
0064         if (x <= x_.back()) {
0065              throw std::domain_error("Calling push_back must preserve the monotonicity of the x's");
0066         }
0067         x_.push_back(x);
0068         y_.push_back(y);
0069         dydx_.push_back(dydx);
0070         d2ydx2_.push_back(d2ydx2);
0071         d3ydx3_.push_back(d3ydx3);
0072     }
0073 
0074     Real operator()(Real x) const
0075     {
0076         if  (x < x_[0] || x > x_.back())
0077         {
0078             std::ostringstream oss;
0079             oss.precision(std::numeric_limits<Real>::digits10+3);
0080             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
0081                 << x_[0] << ", " << x_.back() << "]";
0082             throw std::domain_error(oss.str());
0083         }
0084         // t \in [0, 1)
0085         if (x == x_.back())
0086         {
0087             return y_.back();
0088         }
0089 
0090         auto it = std::upper_bound(x_.begin(), x_.end(), x);
0091         auto i = std::distance(x_.begin(), it) -1;
0092         Real x0 = *(it-1);
0093         Real x1 = *it;
0094         Real dx = (x1-x0);
0095         Real t = (x-x0)/dx;
0096 
0097         // See: 
0098         // http://seisweb.usask.ca/classes/GEOL481/2017/Labs/interpolation_utilities_matlab/shermite.m
0099         Real t2 = t*t;
0100         Real t3 = t2*t;
0101         Real t4 = t3*t;
0102         Real dx2 = dx*dx/2;
0103         Real dx3 = dx2*dx/3;
0104 
0105         Real s = t4*(-35 + t*(84 + t*(-70 + 20*t)));
0106         Real z4 = -s;
0107         Real z0 = s + 1;
0108         Real z1 = t*(1 + t3*(-20 + t*(45 + t*(-36 + 10*t))));
0109         Real z2 = t2*(1 + t2*(-10 + t*(20 + t*(-15 + 4*t))));
0110         Real z3 = t3*(1 + t*(-4 + t*(6 + t*(-4 + t))));
0111         Real z5 = t4*(-15 + t*(39 + t*(-34 + 10*t)));
0112         Real z6 = t4*(5 + t*(-14 + t*(13 - 4*t)));
0113         Real z7 = t4*(-1 + t*(3 + t*(-3+t)));
0114 
0115         Real y0 = y_[i];
0116         Real y1 = y_[i+1];
0117         // Velocity:
0118         Real v0 = dydx_[i];
0119         Real v1 = dydx_[i+1];
0120         // Acceleration:
0121         Real a0 = d2ydx2_[i];
0122         Real a1 = d2ydx2_[i+1];
0123         // Jerk:
0124         Real j0 = d3ydx3_[i];
0125         Real j1 = d3ydx3_[i+1];
0126 
0127         return z0*y0 + z4*y1 + (z1*v0 + z5*v1)*dx + (z2*a0 + z6*a1)*dx2 + (z3*j0 + z7*j1)*dx3;
0128     }
0129 
0130     Real prime(Real x) const
0131     {
0132         if  (x < x_[0] || x > x_.back())
0133         {
0134             std::ostringstream oss;
0135             oss.precision(std::numeric_limits<Real>::digits10+3);
0136             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
0137                 << x_[0] << ", " << x_.back() << "]";
0138             throw std::domain_error(oss.str());
0139         }
0140         if (x == x_.back())
0141         {
0142             return dydx_.back();
0143         }
0144 
0145         auto it = std::upper_bound(x_.begin(), x_.end(), x);
0146         auto i = std::distance(x_.begin(), it) -1;
0147         Real x0 = *(it-1);
0148         Real x1 = *it;
0149         Real y0 = y_[i];
0150         Real y1 = y_[i+1];
0151         Real v0 = dydx_[i];
0152         Real v1 = dydx_[i+1];
0153         Real a0 = d2ydx2_[i];
0154         Real a1 = d2ydx2_[i+1];
0155         Real j0 = d3ydx3_[i];
0156         Real j1 = d3ydx3_[i+1];
0157         Real dx = x1 - x0;
0158         Real t = (x-x0)/dx;
0159         Real t2 = t*t;
0160         Real t3 = t2*t;
0161         Real z0 = 140*t3*(1 + t*(-3 + t*(3 - t)));
0162         Real z1 = 1 + t3*(-80 + t*(225 + t*(-216 + 70*t)));
0163         Real z2 = t3*(-60 + t*(195 + t*(-204 + 70*t)));
0164         Real z3 = 1 + t2*(-20 + t*(50 + t*(-45 + 14*t)));
0165         Real z4 = t2*(10 + t*(-35 + t*(39 - 14*t)));
0166         Real z5 = 3 + t*(-16 + t*(30 + t*(-24 + 7*t)));
0167         Real z6 = t*(-4 + t*(15 + t*(-18 + 7*t)));
0168 
0169         Real dydx = z0*(y1-y0)/dx;
0170         dydx += z1*v0 + z2*v1;
0171         dydx += (x-x0)*(z3*a0 + z4*a1);
0172         dydx += (x-x0)*(x-x0)*(z5*j0 + z6*j1)/6;
0173         return dydx;
0174     }
0175 
0176     inline Real double_prime(Real) const
0177     {
0178         return std::numeric_limits<Real>::quiet_NaN();
0179     }
0180 
0181     friend std::ostream& operator<<(std::ostream & os, const septic_hermite_detail & m)
0182     {
0183         os << "(x,y,y') = {";
0184         for (size_t i = 0; i < m.x_.size() - 1; ++i) {
0185             os << "(" << m.x_[i] << ", " << m.y_[i] << ", " << m.dydx_[i] << ", " << m.d2ydx2_[i] <<  ", " << m.d3ydx3_[i] << "),  ";
0186         }
0187         auto n = m.x_.size()-1;
0188         os << "(" << m.x_[n] << ", " << m.y_[n] << ", " << m.dydx_[n] << ", " << m.d2ydx2_[n] << m.d3ydx3_[n] << ")}";
0189         return os;
0190     }
0191 
0192     int64_t bytes()
0193     {
0194         return 5*x_.size()*sizeof(Real) + 5*sizeof(x_);
0195     }
0196 
0197     std::pair<Real, Real> domain() const
0198     {
0199         return {x_.front(), x_.back()};
0200     }
0201 
0202 private:
0203     RandomAccessContainer x_;
0204     RandomAccessContainer y_;
0205     RandomAccessContainer dydx_;
0206     RandomAccessContainer d2ydx2_;
0207     RandomAccessContainer d3ydx3_;
0208 };
0209 
0210 template<class RandomAccessContainer>
0211 class cardinal_septic_hermite_detail {
0212 public:
0213     using Real = typename RandomAccessContainer::value_type;
0214     cardinal_septic_hermite_detail(RandomAccessContainer && y, RandomAccessContainer && dydx, RandomAccessContainer && d2ydx2, RandomAccessContainer && d3ydx3, Real x0, Real dx)
0215     : y_{std::move(y)}, dy_{std::move(dydx)}, d2y_{std::move(d2ydx2)}, d3y_{std::move(d3ydx3)}, x0_{x0}, inv_dx_{1/dx}
0216     {
0217         if (y_.size() != dy_.size())
0218         {
0219             throw std::domain_error("Numbers of derivatives must = number of ordinates.");
0220         }
0221         if (y_.size() != d2y_.size())
0222         {
0223             throw std::domain_error("Number of second derivatives must equal number of ordinates.");
0224         }
0225         if (y_.size() != d3y_.size())
0226         {
0227             throw std::domain_error("Number of third derivatives must equal number of ordinates.");
0228         }
0229         if (y_.size() < 2)
0230         {
0231             throw std::domain_error("At least 2 abscissas are required.");
0232         }
0233 
0234         if (dx <= 0)
0235         {
0236             throw std::domain_error("dx > 0 is required.");
0237         }
0238 
0239         for (auto & dy : dy_)
0240         {
0241             dy *= dx;
0242         }
0243         for (auto & d2y : d2y_)
0244         {
0245             d2y *= (dx*dx/2);
0246         }
0247         for (auto & d3y : d3y_)
0248         {
0249             d3y *= (dx*dx*dx/6);
0250         }
0251 
0252     }
0253 
0254     inline Real operator()(Real x) const
0255     {
0256         Real xf = x0_ + (y_.size()-1)/inv_dx_;
0257         if  (x < x0_ || x > xf)
0258         {
0259             std::ostringstream oss;
0260             oss.precision(std::numeric_limits<Real>::digits10+3);
0261             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
0262                 << x0_ << ", " << xf << "]";
0263             throw std::domain_error(oss.str());
0264         }
0265         if (x == xf)
0266         {
0267             return y_.back();
0268         }
0269         return this->unchecked_evaluation(x);
0270     }
0271 
0272     inline Real unchecked_evaluation(Real x) const {
0273         using std::floor;
0274         Real s3 = (x-x0_)*inv_dx_;
0275         Real ii = floor(s3);
0276         auto i = static_cast<decltype(y_.size())>(ii);
0277         Real t = s3 - ii;
0278         if (t == 0) {
0279             return y_[i];
0280         }
0281         // See:
0282         // http://seisweb.usask.ca/classes/GEOL481/2017/Labs/interpolation_utilities_matlab/shermite.m
0283         Real t2 = t*t;
0284         Real t3 = t2*t;
0285         Real t4 = t3*t;
0286 
0287         Real s = t4*(-35 + t*(84 + t*(-70 + 20*t)));
0288         Real z4 = -s;
0289         Real z0 = s + 1;
0290         Real z1 = t*(1 + t3*(-20 + t*(45 + t*(-36+10*t))));
0291         Real z2 = t2*(1 + t2*(-10 + t*(20 + t*(-15+4*t))));
0292         Real z3 = t3*(1 + t*(-4+t*(6+t*(-4+t))));
0293         Real z5 = t4*(-15 + t*(39 + t*(-34 + 10*t)));
0294         Real z6 = t4*(5 + t*(-14 + t*(13-4*t)));
0295         Real z7 = t4*(-1 + t*(3+t*(-3+t)));
0296 
0297         Real y0 = y_[i];
0298         Real y1 = y_[i+1];
0299         Real dy0 = dy_[i];
0300         Real dy1 = dy_[i+1];
0301         Real a0 = d2y_[i];
0302         Real a1 = d2y_[i+1];
0303         Real j0 = d3y_[i];
0304         Real j1 = d3y_[i+1];
0305 
0306         return z0*y0 + z1*dy0 + z2*a0 + z3*j0 + z4*y1 + z5*dy1 + z6*a1 + z7*j1;
0307     }
0308 
0309     inline Real prime(Real x) const
0310     {
0311         Real xf = x0_ + (y_.size()-1)/inv_dx_;
0312         if  (x < x0_ || x > xf)
0313         {
0314             std::ostringstream oss;
0315             oss.precision(std::numeric_limits<Real>::digits10+3);
0316             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
0317                 << x0_ << ", " << xf << "]";
0318             throw std::domain_error(oss.str());
0319         }
0320         if (x == xf)
0321         {
0322             return dy_.back()/inv_dx_;
0323         }
0324 
0325         return this->unchecked_prime(x);
0326     }
0327 
0328     inline Real unchecked_prime(Real x) const
0329     {
0330         using std::floor;
0331         Real s3 = (x-x0_)*inv_dx_;
0332         Real ii = floor(s3);
0333         auto i = static_cast<decltype(y_.size())>(ii);
0334         Real t = s3 - ii;
0335         if (t==0)
0336         {
0337             return dy_[i]/inv_dx_;
0338         }
0339  
0340         Real y0 = y_[i];
0341         Real y1 = y_[i+1];
0342         Real dy0 = dy_[i];
0343         Real dy1 = dy_[i+1];
0344         Real a0 = d2y_[i];
0345         Real a1 = d2y_[i+1];
0346         Real j0 = d3y_[i];
0347         Real j1 = d3y_[i+1];
0348         Real t2 = t*t;
0349         Real t3 = t2*t;
0350         Real z0 = 140*t3*(1 + t*(-3 + t*(3 - t)));
0351         Real z1 = 1 + t3*(-80 + t*(225 + t*(-216 + 70*t)));
0352         Real z2 = t3*(-60 + t*(195 + t*(-204 + 70*t)));
0353         Real z3 = 1 + t2*(-20 + t*(50 + t*(-45 + 14*t)));
0354         Real z4 = t2*(10 + t*(-35 + t*(39 - 14*t)));
0355         Real z5 = 3 + t*(-16 + t*(30 + t*(-24 + 7*t)));
0356         Real z6 = t*(-4 + t*(15 + t*(-18 + 7*t)));
0357 
0358         Real dydx = z0*(y1-y0)*inv_dx_;
0359         dydx += (z1*dy0 + z2*dy1)*inv_dx_;
0360         dydx += 2*t*(z3*a0 + z4*a1)*inv_dx_;
0361         dydx += t*t*(z5*j0 + z6*j1);
0362         return dydx;
0363     }
0364 
0365     inline Real double_prime(Real x) const
0366     {
0367         Real xf = x0_ + (y_.size()-1)/inv_dx_;
0368         if  (x < x0_ || x > xf)
0369         {
0370             std::ostringstream oss;
0371             oss.precision(std::numeric_limits<Real>::digits10+3);
0372             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
0373                 << x0_ << ", " << xf << "]";
0374             throw std::domain_error(oss.str());
0375         }
0376         if (x == xf)
0377         {
0378             return d2y_.back()*2*inv_dx_*inv_dx_;
0379         }
0380 
0381         return this->unchecked_double_prime(x);
0382     }
0383 
0384     inline Real unchecked_double_prime(Real x) const
0385     {
0386         using std::floor;
0387         Real s3 = (x-x0_)*inv_dx_;
0388         Real ii = floor(s3);
0389         auto i = static_cast<decltype(y_.size())>(ii);
0390         Real t = s3 - ii;
0391         if (t==0)
0392         {
0393             return d2y_[i]*2*inv_dx_*inv_dx_;
0394         }
0395 
0396         Real y0 = y_[i];
0397         Real y1 = y_[i+1];
0398         Real dy0 = dy_[i];
0399         Real dy1 = dy_[i+1];
0400         Real a0 = d2y_[i];
0401         Real a1 = d2y_[i+1];
0402         Real j0 = d3y_[i];
0403         Real j1 = d3y_[i+1];
0404         Real t2 = t*t;
0405 
0406         Real z0 = 420*t2*(1 + t*(-4 + t*(5 - 2*t)));
0407         Real z1 = 60*t2*(-4 + t*(15 + t*(-18 + 7*t)));
0408         Real z2 = 60*t2*(-3 + t*(13 + t*(-17 + 7*t)));
0409         Real z3 = (1 + t2*(-60 + t*(200 + t*(-225 + 84*t))));
0410         Real z4 = t2*(30 + t*(-140 + t*(195 - 84*t)));
0411         Real z5 = t*(1 + t*(-8 + t*(20 + t*(-20 + 7*t))));
0412         Real z6 = t2*(-2 + t*(10 + t*(-15 + 7*t)));
0413 
0414         Real d2ydx2 = z0*(y1-y0)*inv_dx_*inv_dx_;
0415         d2ydx2 += (z1*dy0 + z2*dy1)*inv_dx_*inv_dx_;
0416         d2ydx2 += (z3*a0 + z4*a1)*2*inv_dx_*inv_dx_;
0417         d2ydx2 += 6*(z5*j0 + z6*j1)/(inv_dx_*inv_dx_);
0418 
0419         return d2ydx2;
0420     }
0421 
0422     int64_t bytes() const
0423     {
0424         return 4*y_.size()*sizeof(Real) + 2*sizeof(Real) + 4*sizeof(y_);
0425     }
0426 
0427     std::pair<Real, Real> domain() const
0428     {
0429         return {x0_, x0_ + (y_.size()-1)/inv_dx_};
0430     }
0431 
0432 private:
0433     RandomAccessContainer y_;
0434     RandomAccessContainer dy_;
0435     RandomAccessContainer d2y_;
0436     RandomAccessContainer d3y_;
0437     Real x0_;
0438     Real inv_dx_;
0439 };
0440 
0441 
0442 template<class RandomAccessContainer>
0443 class cardinal_septic_hermite_detail_aos {
0444 public:
0445     using Point = typename RandomAccessContainer::value_type;
0446     using Real = typename Point::value_type;
0447     cardinal_septic_hermite_detail_aos(RandomAccessContainer && dat, Real x0, Real dx)
0448     : data_{std::move(dat)}, x0_{x0}, inv_dx_{1/dx}
0449     {
0450         if (data_.size() < 2) {
0451             throw std::domain_error("At least 2 abscissas are required.");
0452         }
0453         if (data_[0].size() != 4) {
0454             throw std::domain_error("There must be 4 data items per struct.");
0455         }
0456 
0457         for (auto & datum : data_)
0458         {
0459             datum[1] *= dx;
0460             datum[2] *= (dx*dx/2);
0461             datum[3] *= (dx*dx*dx/6);
0462         }
0463     }
0464 
0465     inline Real operator()(Real x) const
0466     {
0467         Real xf = x0_ + (data_.size()-1)/inv_dx_;
0468         if  (x < x0_ || x > xf)
0469         {
0470             std::ostringstream oss;
0471             oss.precision(std::numeric_limits<Real>::digits10+3);
0472             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
0473                 << x0_ << ", " << xf << "]";
0474             throw std::domain_error(oss.str());
0475         }
0476         if (x == xf)
0477         {
0478             return data_.back()[0];
0479         }
0480         return this->unchecked_evaluation(x);
0481     }
0482 
0483     inline Real unchecked_evaluation(Real x) const
0484     {
0485         using std::floor;
0486         Real s3 = (x-x0_)*inv_dx_;
0487         Real ii = floor(s3);
0488         auto i = static_cast<decltype(data_.size())>(ii);
0489         Real t = s3 - ii;
0490         if (t==0)
0491         {
0492             return data_[i][0];
0493         }
0494         Real t2 = t*t;
0495         Real t3 = t2*t;
0496         Real t4 = t3*t;
0497 
0498         Real s = t4*(-35 + t*(84 + t*(-70 + 20*t)));
0499         Real z4 = -s;
0500         Real z0 = s + 1;
0501         Real z1 = t*(1 + t3*(-20 + t*(45 + t*(-36+10*t))));
0502         Real z2 = t2*(1 + t2*(-10 + t*(20 + t*(-15+4*t))));
0503         Real z3 = t3*(1 + t*(-4+t*(6+t*(-4+t))));
0504         Real z5 = t4*(-15 + t*(39 + t*(-34 + 10*t)));
0505         Real z6 = t4*(5 + t*(-14 + t*(13-4*t)));
0506         Real z7 = t4*(-1 + t*(3+t*(-3+t)));
0507 
0508         Real y0 = data_[i][0];
0509         Real dy0 = data_[i][1];
0510         Real a0 = data_[i][2];
0511         Real j0 = data_[i][3];
0512         Real y1 = data_[i+1][0];
0513         Real dy1 = data_[i+1][1];
0514         Real a1 = data_[i+1][2];
0515         Real j1 = data_[i+1][3];
0516 
0517         return z0*y0 + z1*dy0 + z2*a0 + z3*j0 + z4*y1 + z5*dy1 + z6*a1 + z7*j1;
0518     }
0519 
0520     inline Real prime(Real x) const
0521     {
0522         Real xf = x0_ + (data_.size()-1)/inv_dx_;
0523         if  (x < x0_ || x > xf)
0524         {
0525             std::ostringstream oss;
0526             oss.precision(std::numeric_limits<Real>::digits10+3);
0527             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
0528                 << x0_ << ", " << xf << "]";
0529             throw std::domain_error(oss.str());
0530         }
0531         if (x == xf)
0532         {
0533             return data_.back()[1]*inv_dx_;
0534         }
0535 
0536         return this->unchecked_prime(x);
0537     }
0538 
0539     inline Real unchecked_prime(Real x) const
0540     {
0541         using std::floor;
0542         Real s3 = (x-x0_)*inv_dx_;
0543         Real ii = floor(s3);
0544         auto i = static_cast<decltype(data_.size())>(ii);
0545         Real t = s3 - ii;
0546         if (t == 0)
0547         {
0548             return data_[i][1]*inv_dx_;
0549         }
0550 
0551         Real y0 = data_[i][0];
0552         Real y1 = data_[i+1][0];
0553         Real dy0 = data_[i][1];
0554         Real dy1 = data_[i+1][1];
0555         Real a0 = data_[i][2];
0556         Real a1 = data_[i+1][2];
0557         Real j0 = data_[i][3];
0558         Real j1 = data_[i+1][3];
0559         Real t2 = t*t;
0560         Real t3 = t2*t;
0561         Real z0 = 140*t3*(1 + t*(-3 + t*(3 - t)));
0562         Real z1 = 1 + t3*(-80 + t*(225 + t*(-216 + 70*t)));
0563         Real z2 = t3*(-60 + t*(195 + t*(-204 + 70*t)));
0564         Real z3 = 1 + t2*(-20 + t*(50 + t*(-45 + 14*t)));
0565         Real z4 = t2*(10 + t*(-35 + t*(39 - 14*t)));
0566         Real z5 = 3 + t*(-16 + t*(30 + t*(-24 + 7*t)));
0567         Real z6 = t*(-4 + t*(15 + t*(-18 + 7*t)));
0568 
0569         Real dydx = z0*(y1-y0)*inv_dx_;
0570         dydx += (z1*dy0 + z2*dy1)*inv_dx_;
0571         dydx += 2*t*(z3*a0 + z4*a1)*inv_dx_;
0572         dydx += t*t*(z5*j0 + z6*j1);
0573         return dydx;
0574     }
0575 
0576     inline Real double_prime(Real x) const
0577     {
0578         Real xf = x0_ + (data_.size()-1)/inv_dx_;
0579         if  (x < x0_ || x > xf)
0580         {
0581             std::ostringstream oss;
0582             oss.precision(std::numeric_limits<Real>::digits10+3);
0583             oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
0584                 << x0_ << ", " << xf << "]";
0585             throw std::domain_error(oss.str());
0586         }
0587         if (x == xf)
0588         {
0589             return data_.back()[2]*2*inv_dx_*inv_dx_;
0590         }
0591 
0592         return this->unchecked_double_prime(x);
0593     }
0594 
0595     inline Real unchecked_double_prime(Real x) const
0596     {
0597         using std::floor;
0598         Real s3 = (x-x0_)*inv_dx_;
0599         Real ii = floor(s3);
0600         auto i = static_cast<decltype(data_.size())>(ii);
0601         Real t = s3 - ii;
0602         if (t == 0)
0603         {
0604             return data_[i][2]*2*inv_dx_*inv_dx_;
0605         }
0606         Real y0 = data_[i][0];
0607         Real y1 = data_[i+1][0];
0608         Real dy0 = data_[i][1];
0609         Real dy1 = data_[i+1][1];
0610         Real a0 = data_[i][2];
0611         Real a1 = data_[i+1][2];
0612         Real j0 = data_[i][3];
0613         Real j1 = data_[i+1][3];
0614         Real t2 = t*t;
0615 
0616         Real z0 = 420*t2*(1 + t*(-4 + t*(5 - 2*t)));
0617         Real z1 = 60*t2*(-4 + t*(15 + t*(-18 + 7*t)));
0618         Real z2 = 60*t2*(-3 + t*(13 + t*(-17 + 7*t)));
0619         Real z3 = (1 + t2*(-60 + t*(200 + t*(-225 + 84*t))));
0620         Real z4 = t2*(30 + t*(-140 + t*(195 - 84*t)));
0621         Real z5 = t*(1 + t*(-8 + t*(20 + t*(-20 + 7*t))));
0622         Real z6 = t2*(-2 + t*(10 + t*(-15 + 7*t)));
0623 
0624         Real d2ydx2 = z0*(y1-y0)*inv_dx_*inv_dx_;
0625         d2ydx2 += (z1*dy0 + z2*dy1)*inv_dx_*inv_dx_;
0626         d2ydx2 += (z3*a0 + z4*a1)*2*inv_dx_*inv_dx_;
0627         d2ydx2 += 6*(z5*j0 + z6*j1)/(inv_dx_*inv_dx_);
0628 
0629         return d2ydx2;
0630     }
0631 
0632     int64_t bytes() const
0633     {
0634         return data_.size()*data_[0].size()*sizeof(Real) + 2*sizeof(Real) + sizeof(data_);
0635     }
0636 
0637     std::pair<Real, Real> domain() const
0638     {
0639         return {x0_, x0_ + (data_.size() -1)/inv_dx_};
0640     }
0641 
0642 private:
0643     RandomAccessContainer data_;
0644     Real x0_;
0645     Real inv_dx_;
0646 };
0647 
0648 }
0649 }
0650 }
0651 }
0652 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!