source: sasview/src/sas/models/c_extension/cephes/yn.c @ c81f9b6

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since c81f9b6 was 79492222, checked in by krzywon, 10 years ago

Changed the file and folder names to remove all SANS references.

  • Property mode set to 100644
File size: 1.8 KB
Line 
1/*                                                      yn.c
2 *
3 *      Bessel function of second kind of integer order
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * double x, y, yn();
10 * int n;
11 *
12 * y = yn( n, x );
13 *
14 *
15 *
16 * DESCRIPTION:
17 *
18 * Returns Bessel function of order n, where n is a
19 * (possibly negative) integer.
20 *
21 * The function is evaluated by forward recurrence on
22 * n, starting with values computed by the routines
23 * y0() and y1().
24 *
25 * If n = 0 or 1 the routine for y0 or y1 is called
26 * directly.
27 *
28 *
29 *
30 * ACCURACY:
31 *
32 *
33 *                      Absolute error, except relative
34 *                      when y > 1:
35 * arithmetic   domain     # trials      peak         rms
36 *    DEC       0, 30        2200       2.9e-16     5.3e-17
37 *    IEEE      0, 30       30000       3.4e-15     4.3e-16
38 *
39 *
40 * ERROR MESSAGES:
41 *
42 *   message         condition      value returned
43 * yn singularity   x = 0              MAXNUM
44 * yn overflow                         MAXNUM
45 *
46 * Spot checked against tables for x, n between 0 and 100.
47 *
48 */
49
50/*
51Cephes Math Library Release 2.8:  June, 2000
52Copyright 1984, 1987, 2000 by Stephen L. Moshier
53*/
54
55#include "mconf.h"
56#ifdef ANSIPROT
57extern double y0 ( double );
58extern double y1 ( double );
59extern double log ( double );
60#else
61double y0(), y1(), log();
62#endif
63extern double MAXNUM, MAXLOG;
64
65double yn( n, x )
66int n;
67double x;
68{
69double an, anm1, anm2, r;
70int k, sign;
71
72if( n < 0 )
73        {
74        n = -n;
75        if( (n & 1) == 0 )      /* -1**n */
76                sign = 1;
77        else
78                sign = -1;
79        }
80else
81        sign = 1;
82
83
84if( n == 0 )
85        return( sign * y0(x) );
86if( n == 1 )
87        return( sign * y1(x) );
88
89/* test for overflow */
90if( x <= 0.0 )
91        {
92        mtherr( "yn", SING );
93        return( -MAXNUM );
94        }
95
96/* forward recurrence on n */
97
98anm2 = y0(x);
99anm1 = y1(x);
100k = 1;
101r = 2 * k;
102do
103        {
104        an = r * anm1 / x  -  anm2;
105        anm2 = anm1;
106        anm1 = an;
107        r += 2.0;
108        ++k;
109        }
110while( k < n );
111
112
113return( sign * an );
114}
Note: See TracBrowser for help on using the repository browser.