source: sasview/src/sas/models/c_extension/cephes/iv.c @ 79492222

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 79492222 was 79492222, checked in by krzywon, 9 years ago

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

  • Property mode set to 100644
File size: 2.2 KB
Line 
1/*                                                      iv.c
2 *
3 *      Modified Bessel function of noninteger order
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * double v, x, y, iv();
10 *
11 * y = iv( v, x );
12 *
13 *
14 *
15 * DESCRIPTION:
16 *
17 * Returns modified Bessel function of order v of the
18 * argument.  If x is negative, v must be integer valued.
19 *
20 * The function is defined as Iv(x) = Jv( ix ).  It is
21 * here computed in terms of the confluent hypergeometric
22 * function, according to the formula
23 *
24 *              v  -x
25 * Iv(x) = (x/2)  e   hyperg( v+0.5, 2v+1, 2x ) / gamma(v+1)
26 *
27 * If v is a negative integer, then v is replaced by -v.
28 *
29 *
30 * ACCURACY:
31 *
32 * Tested at random points (v, x), with v between 0 and
33 * 30, x between 0 and 28.
34 *                      Relative error:
35 * arithmetic   domain     # trials      peak         rms
36 *    DEC       0,30          2000      3.1e-15     5.4e-16
37 *    IEEE      0,30         10000      1.7e-14     2.7e-15
38 *
39 * Accuracy is diminished if v is near a negative integer.
40 *
41 * See also hyperg.c.
42 *
43 */
44/*                                                     iv.c    */
45/*      Modified Bessel function of noninteger order            */
46/* If x < 0, then v must be an integer. */
47
48
49/*
50Cephes Math Library Release 2.8:  June, 2000
51Copyright 1984, 1987, 1988, 2000 by Stephen L. Moshier
52*/
53
54
55#include "mconf.h"
56#ifdef ANSIPROT
57extern double hyperg ( double, double, double );
58extern double exp ( double );
59extern double gamma ( double );
60extern double log ( double );
61extern double fabs ( double );
62extern double floor ( double );
63#else
64double hyperg(), exp(), gamma(), log(), fabs(), floor();
65#endif
66extern double MACHEP, MAXNUM;
67
68double iv( v, x )
69double v, x;
70{
71int sign;
72double t, ax;
73
74/* If v is a negative integer, invoke symmetry */
75t = floor(v);
76if( v < 0.0 )
77        {
78        if( t == v )
79                {
80                v = -v; /* symmetry */
81                t = -t;
82                }
83        }
84/* If x is negative, require v to be an integer */
85sign = 1;
86if( x < 0.0 )
87        {
88        if( t != v )
89                {
90                mtherr( "iv", DOMAIN );
91                return( 0.0 );
92                }
93        if( v != 2.0 * floor(v/2.0) )
94                sign = -1;
95        }
96
97/* Avoid logarithm singularity */
98if( x == 0.0 )
99        {
100        if( v == 0.0 )
101                return( 1.0 );
102        if( v < 0.0 )
103                {
104                mtherr( "iv", OVERFLOW );
105                return( MAXNUM );
106                }
107        else
108                return( 0.0 );
109        }
110
111ax = fabs(x);
112t = v * log( 0.5 * ax )  -  x;
113t = sign * exp(t) / gamma( v + 1.0 );
114ax = v + 0.5;
115return( t * hyperg( ax,  2.0 * ax,  2.0 * x ) );
116}
Note: See TracBrowser for help on using the repository browser.