| 1 | /* |
| 2 | floating point Bessel's function |
| 3 | of the first and second kinds |
| 4 | of order one |
| 5 | |
| 6 | j1(x) returns the value of J1(x) |
| 7 | for all real values of x. |
| 8 | |
| 9 | There are no error returns. |
| 10 | Calls sin, cos, sqrt. |
| 11 | |
| 12 | There is a niggling bug in J1 which |
| 13 | causes errors up to 2e-16 for x in the |
| 14 | interval [-8,8]. |
| 15 | The bug is caused by an inappropriate order |
| 16 | of summation of the series. rhm will fix it |
| 17 | someday. |
| 18 | |
| 19 | Coefficients are from Hart & Cheney. |
| 20 | #6050 (20.98D) |
| 21 | #6750 (19.19D) |
| 22 | #7150 (19.35D) |
| 23 | |
| 24 | y1(x) returns the value of Y1(x) |
| 25 | for positive real values of x. |
| 26 | For x<=0, error number EDOM is set and a |
| 27 | large negative value is returned. |
| 28 | |
| 29 | Calls sin, cos, sqrt, log, j1. |
| 30 | |
| 31 | The values of Y1 have not been checked |
| 32 | to more than ten places. |
| 33 | |
| 34 | Coefficients are from Hart & Cheney. |
| 35 | #6447 (22.18D) |
| 36 | #6750 (19.19D) |
| 37 | #7150 (19.35D) |
| 38 | */ |
| 39 | |
| 40 | #include <math.h> |
| 41 | #include <errno.h> |
| 42 | |
| 43 | int errno; |
| 44 | static double pzero, qzero; |
| 45 | static double tpi = .6366197723675813430755350535e0; |
| 46 | static double pio4 = .7853981633974483096156608458e0; |
| 47 | static double p1[] = { |
| 48 | 0.581199354001606143928050809e21, |
| 49 | -.6672106568924916298020941484e20, |
| 50 | 0.2316433580634002297931815435e19, |
| 51 | -.3588817569910106050743641413e17, |
| 52 | 0.2908795263834775409737601689e15, |
| 53 | -.1322983480332126453125473247e13, |
| 54 | 0.3413234182301700539091292655e10, |
| 55 | -.4695753530642995859767162166e7, |
| 56 | 0.2701122710892323414856790990e4, |
| 57 | }; |
| 58 | static double q1[] = { |
| 59 | 0.1162398708003212287858529400e22, |
| 60 | 0.1185770712190320999837113348e20, |
| 61 | 0.6092061398917521746105196863e17, |
| 62 | 0.2081661221307607351240184229e15, |
| 63 | 0.5243710262167649715406728642e12, |
| 64 | 0.1013863514358673989967045588e10, |
| 65 | 0.1501793594998585505921097578e7, |
| 66 | 0.1606931573481487801970916749e4, |
| 67 | 1.0, |
| 68 | }; |
| 69 | static double p2[] = { |
| 70 | -.4435757816794127857114720794e7, |
| 71 | -.9942246505077641195658377899e7, |
| 72 | -.6603373248364939109255245434e7, |
| 73 | -.1523529351181137383255105722e7, |
| 74 | -.1098240554345934672737413139e6, |
| 75 | -.1611616644324610116477412898e4, |
| 76 | 0.0, |
| 77 | }; |
| 78 | static double q2[] = { |
| 79 | -.4435757816794127856828016962e7, |
| 80 | -.9934124389934585658967556309e7, |
| 81 | -.6585339479723087072826915069e7, |
| 82 | -.1511809506634160881644546358e7, |
| 83 | -.1072638599110382011903063867e6, |
| 84 | -.1455009440190496182453565068e4, |
| 85 | 1.0, |
| 86 | }; |
| 87 | static double p3[] = { |
| 88 | 0.3322091340985722351859704442e5, |
| 89 | 0.8514516067533570196555001171e5, |
| 90 | 0.6617883658127083517939992166e5, |
| 91 | 0.1849426287322386679652009819e5, |
| 92 | 0.1706375429020768002061283546e4, |
| 93 | 0.3526513384663603218592175580e2, |
| 94 | 0.0, |
| 95 | }; |
| 96 | static double q3[] = { |
| 97 | 0.7087128194102874357377502472e6, |
| 98 | 0.1819458042243997298924553839e7, |
| 99 | 0.1419460669603720892855755253e7, |
| 100 | 0.4002944358226697511708610813e6, |
| 101 | 0.3789022974577220264142952256e5, |
| 102 | 0.8638367769604990967475517183e3, |
| 103 | 1.0, |
| 104 | }; |
| 105 | static double p4[] = { |
| 106 | -.9963753424306922225996744354e23, |
| 107 | 0.2655473831434854326894248968e23, |
| 108 | -.1212297555414509577913561535e22, |
| 109 | 0.2193107339917797592111427556e20, |
| 110 | -.1965887462722140658820322248e18, |
| 111 | 0.9569930239921683481121552788e15, |
| 112 | -.2580681702194450950541426399e13, |
| 113 | 0.3639488548124002058278999428e10, |
| 114 | -.2108847540133123652824139923e7, |
| 115 | 0.0, |
| 116 | }; |
| 117 | static double q4[] = { |
| 118 | 0.5082067366941243245314424152e24, |
| 119 | 0.5435310377188854170800653097e22, |
| 120 | 0.2954987935897148674290758119e20, |
| 121 | 0.1082258259408819552553850180e18, |
| 122 | 0.2976632125647276729292742282e15, |
| 123 | 0.6465340881265275571961681500e12, |
| 124 | 0.1128686837169442121732366891e10, |
| 125 | 0.1563282754899580604737366452e7, |
| 126 | 0.1612361029677000859332072312e4, |
| 127 | 1.0, |
| 128 | }; |
| 129 | |
| 130 | double |
| 131 | j1(arg) double arg;{ |
| 132 | double xsq, n, d, x; |
| 133 | double sin(), cos(), sqrt(); |
| 134 | int i; |
| 135 | |
| 136 | x = arg; |
| 137 | if(x < 0.) x = -x; |
| 138 | if(x > 8.){ |
| 139 | asympt(x); |
| 140 | n = x - 3.*pio4; |
| 141 | n = sqrt(tpi/x)*(pzero*cos(n) - qzero*sin(n)); |
| 142 | if(arg <0.) n = -n; |
| 143 | return(n); |
| 144 | } |
| 145 | xsq = x*x; |
| 146 | for(n=0,d=0,i=8;i>=0;i--){ |
| 147 | n = n*xsq + p1[i]; |
| 148 | d = d*xsq + q1[i]; |
| 149 | } |
| 150 | return(arg*n/d); |
| 151 | } |
| 152 | |
| 153 | double |
| 154 | y1(arg) double arg;{ |
| 155 | double xsq, n, d, x; |
| 156 | double sin(), cos(), sqrt(), log(), j1(); |
| 157 | int i; |
| 158 | |
| 159 | errno = 0; |
| 160 | x = arg; |
| 161 | if(x <= 0.){ |
| 162 | errno = EDOM; |
| 163 | return(-HUGE); |
| 164 | } |
| 165 | if(x > 8.){ |
| 166 | asympt(x); |
| 167 | n = x - 3*pio4; |
| 168 | return(sqrt(tpi/x)*(pzero*sin(n) + qzero*cos(n))); |
| 169 | } |
| 170 | xsq = x*x; |
| 171 | for(n=0,d=0,i=9;i>=0;i--){ |
| 172 | n = n*xsq + p4[i]; |
| 173 | d = d*xsq + q4[i]; |
| 174 | } |
| 175 | return(x*n/d + tpi*(j1(x)*log(x)-1./x)); |
| 176 | } |
| 177 | |
| 178 | static |
| 179 | asympt(arg) double arg;{ |
| 180 | double zsq, n, d; |
| 181 | int i; |
| 182 | zsq = 64./(arg*arg); |
| 183 | for(n=0,d=0,i=6;i>=0;i--){ |
| 184 | n = n*zsq + p2[i]; |
| 185 | d = d*zsq + q2[i]; |
| 186 | } |
| 187 | pzero = n/d; |
| 188 | for(n=0,d=0,i=6;i>=0;i--){ |
| 189 | n = n*zsq + p3[i]; |
| 190 | d = d*zsq + q3[i]; |
| 191 | } |
| 192 | qzero = (8./arg)*(n/d); |
| 193 | } |