special_functions/detail/bessel_jy_derivatives_asym.hpp

0001 //  Copyright (c) 2013 Anton Bikineev
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 //
0007 // This is a partial header, do not include on it's own!!!
0008 //
0009 // Contains asymptotic expansions for derivatives of Bessel J(v,x) and Y(v,x)
0010 // functions, as x -> INF.
0011 #ifndef BOOST_MATH_SF_DETAIL_BESSEL_JY_DERIVATIVES_ASYM_HPP
0012 #define BOOST_MATH_SF_DETAIL_BESSEL_JY_DERIVATIVES_ASYM_HPP
0013
0014 #ifdef _MSC_VER
0015 #pragma once
0016 #endif
0017
0018 namespace boost{ namespace math{ namespace detail{
0019
0020 template <class T>
0021 inline T asymptotic_bessel_derivative_amplitude(T v, T x)
0022 {
0023    // Calculate the amplitude for J'(v,x) and I'(v,x)
0024    // for large x: see A&S 9.2.30.
0025    BOOST_MATH_STD_USING
0026    T s = 1;
0027    const T mu = 4 * v * v;
0028    T txq = 2 * x;
0029    txq *= txq;
0030
0031    s -= (mu - 3) / (2 * txq);
0032    s -= ((mu - 1) * (mu - 45)) / (txq * txq * 8);
0033
0034    return sqrt(s * 2 / (boost::math::constants::pi<T>() * x));
0035 }
0036
0037 template <class T>
0038 inline T asymptotic_bessel_derivative_phase_mx(T v, T x)
0039 {
0040    // Calculate the phase of J'(v, x) and Y'(v, x) for large x.
0041    // See A&S 9.2.31.
0042    // Note that the result returned is the phase less (x - PI(v/2 - 1/4))
0043    // which we'll factor in later when we calculate the sines/cosines of the result:
0044    const T mu = 4 * v * v;
0045    const T mu2 = mu * mu;
0046    const T mu3 = mu2 * mu;
0047    T denom = 4 * x;
0048    T denom_mult = denom * denom;
0049
0050    T s = 0;
0051    s += (mu + 3) / (2 * denom);
0052    denom *= denom_mult;
0053    s += (mu2 + (46 * mu) - 63) / (6 * denom);
0054    denom *= denom_mult;
0055    s += (mu3 + (185 * mu2) - (2053 * mu) + 1899) / (5 * denom);
0056    return s;
0057 }
0058
0059 template <class T, class Policy>
0060 inline T asymptotic_bessel_y_derivative_large_x_2(T v, T x, const Policy& pol)
0061 {
0062    // See A&S 9.2.20.
0063    BOOST_MATH_STD_USING
0064    // Get the phase and amplitude:
0065    const T ampl = asymptotic_bessel_derivative_amplitude(v, x);
0066    const T phase = asymptotic_bessel_derivative_phase_mx(v, x);
0067    BOOST_MATH_INSTRUMENT_VARIABLE(ampl);
0068    BOOST_MATH_INSTRUMENT_VARIABLE(phase);
0069    //
0070    // Calculate the sine of the phase, using
0071    // sine/cosine addition rules to factor in
0072    // the x - PI(v/2 - 1/4) term not added to the
0073    // phase when we calculated it.
0074    //
0075    const T cx = cos(x);
0076    const T sx = sin(x);
0077    const T vd2shifted = (v / 2) - 0.25f;
0078    const T ci = cos_pi(vd2shifted, pol);
0079    const T si = sin_pi(vd2shifted, pol);
0080    const T sin_phase = sin(phase) * (cx * ci + sx * si) + cos(phase) * (sx * ci - cx * si);
0081    BOOST_MATH_INSTRUMENT_CODE(sin(phase));
0082    BOOST_MATH_INSTRUMENT_CODE(cos(x));
0083    BOOST_MATH_INSTRUMENT_CODE(cos(phase));
0084    BOOST_MATH_INSTRUMENT_CODE(sin(x));
0085    return sin_phase * ampl;
0086 }
0087
0088 template <class T, class Policy>
0089 inline T asymptotic_bessel_j_derivative_large_x_2(T v, T x, const Policy& pol)
0090 {
0091    // See A&S 9.2.20.
0092    BOOST_MATH_STD_USING
0093    // Get the phase and amplitude:
0094    const T ampl = asymptotic_bessel_derivative_amplitude(v, x);
0095    const T phase = asymptotic_bessel_derivative_phase_mx(v, x);
0096    BOOST_MATH_INSTRUMENT_VARIABLE(ampl);
0097    BOOST_MATH_INSTRUMENT_VARIABLE(phase);
0098    //
0099    // Calculate the sine of the phase, using
0100    // sine/cosine addition rules to factor in
0101    // the x - PI(v/2 - 1/4) term not added to the
0102    // phase when we calculated it.
0103    //
0104    BOOST_MATH_INSTRUMENT_CODE(cos(phase));
0105    BOOST_MATH_INSTRUMENT_CODE(cos(x));
0106    BOOST_MATH_INSTRUMENT_CODE(sin(phase));
0107    BOOST_MATH_INSTRUMENT_CODE(sin(x));
0108    const T cx = cos(x);
0109    const T sx = sin(x);
0110    const T vd2shifted = (v / 2) - 0.25f;
0111    const T ci = cos_pi(vd2shifted, pol);
0112    const T si = sin_pi(vd2shifted, pol);
0113    const T sin_phase = cos(phase) * (cx * ci + sx * si) - sin(phase) * (sx * ci - cx * si);
0114    BOOST_MATH_INSTRUMENT_VARIABLE(sin_phase);
0115    return sin_phase * ampl;
0116 }
0117
0118 template <class T>
0119 inline bool asymptotic_bessel_derivative_large_x_limit(const T& v, const T& x)
0120 {
0121    BOOST_MATH_STD_USING
0122    //
0123    // This function is the copy of math::asymptotic_bessel_large_x_limit
0124    // It means that we use the same rules for determining how x is large
0125    // compared to v.
0126    //
0127    // Determines if x is large enough compared to v to take the asymptotic
0128    // forms above.  From A&S 9.2.28 we require:
0129    //    v < x * eps^1/8
0130    // and from A&S 9.2.29 we require:
0131    //    v^12/10 < 1.5 * x * eps^1/10
0132    // using the former seems to work OK in practice with broadly similar
0133    // error rates either side of the divide for v < 10000.
0134    // At double precision eps^1/8 ~= 0.01.
0135    //
0136    return (std::max)(T(fabs(v)), T(1)) < x * sqrt(boost::math::tools::forth_root_epsilon<T>());
0137 }
0138
0139 }}} // namespaces
0140
0141 #endif // BOOST_MATH_SF_DETAIL_BESSEL_JY_DERIVATIVES_ASYM_HPP