← Back to branch summary

~ubuntu-branches/ubuntu/saucy/deal.ii/saucy

« back to all changes in this revision

Viewing changes to contrib/boost/include/boost/math/special_functions/detail/bessel_yn.hpp

  • Committer: Bazaar Package Importer
  • Author(s): Adam C. Powell, IV
  • Date: 2009-05-08 23:13:50 UTC
  • Revision ID: james.westby@ubuntu.com-20090508231350-rrh1ltgi0tifabwc
Tags: upstream-6.2.0
Import upstream version 6.2.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
contrib/boost/include/boost/math/special_functions/detail/bessel_yn.hpp
 
1
//  Copyright (c) 2006 Xiaogang Zhang
 
2
//  Use, modification and distribution are subject to the
 
3
//  Boost Software License, Version 1.0. (See accompanying file
 
4
//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
 
5
 
 
6
#ifndef BOOST_MATH_BESSEL_YN_HPP
 
7
#define BOOST_MATH_BESSEL_YN_HPP
 
8
 
 
9
#ifdef _MSC_VER
 
10
#pragma once
 
11
#endif
 
12
 
 
13
#include <boost/math/special_functions/detail/bessel_y0.hpp>
 
14
#include <boost/math/special_functions/detail/bessel_y1.hpp>
 
15
#include <boost/math/policies/error_handling.hpp>
 
16
 
 
17
// Bessel function of the second kind of integer order
 
18
// Y_n(z) is the dominant solution, forward recurrence always OK (though unstable)
 
19
 
 
20
namespace boost { namespace math { namespace detail{
 
21
 
 
22
template <typename T, typename Policy>
 
23
T bessel_yn(int n, T x, const Policy& pol)
 
24
{
 
25
    T value, factor, current, prev;
 
26
 
 
27
    using namespace boost::math::tools;
 
28
 
 
29
    static const char* function = "boost::math::bessel_yn<%1%>(%1%,%1%)";
 
30
 
 
31
    if ((x == 0) && (n == 0))
 
32
    {
 
33
       return -policies::raise_overflow_error<T>(function, 0, pol);
 
34
    }
 
35
    if (x <= 0)
 
36
    {
 
37
       return policies::raise_domain_error<T>(function,
 
38
            "Got x = %1%, but x must be > 0, complex result not supported.", x, pol);
 
39
    }
 
40
 
 
41
    //
 
42
    // Reflection comes first:
 
43
    //
 
44
    if (n < 0)
 
45
    {
 
46
        factor = (n & 0x1) ? -1 : 1;  // Y_{-n}(z) = (-1)^n Y_n(z)
 
47
        n = -n;
 
48
    }
 
49
    else
 
50
    {
 
51
        factor = 1;
 
52
    }
 
53
 
 
54
    if (n == 0)
 
55
    {
 
56
        value = bessel_y0(x, pol);
 
57
    }
 
58
    else if (n == 1)
 
59
    {
 
60
        value = factor * bessel_y1(x, pol);
 
61
    }
 
62
    else
 
63
    {
 
64
       prev = bessel_y0(x, pol);
 
65
       current = bessel_y1(x, pol);
 
66
       int k = 1;
 
67
       BOOST_ASSERT(k < n);
 
68
       do
 
69
       {
 
70
           value = 2 * k * current / x - prev;
 
71
           prev = current;
 
72
           current = value;
 
73
           ++k;
 
74
       }
 
75
       while(k < n);
 
76
       value *= factor;
 
77
    }
 
78
    return value;
 
79
}
 
80
 
 
81
}}} // namespaces
 
82
 
 
83
#endif // BOOST_MATH_BESSEL_YN_HPP
 
84