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