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 zero 7 8 j0(x) returns the value of J0(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 J0 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 #5849 (19.22D) 23 #6549 (19.25D) 24 #6949 (19.41D) 25 26 y0(x) returns the value of Y0(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, j0. 32 33 The values of Y0 have not been checked 34 to more than ten places. 35 36 Coefficients are from Hart & Cheney. 37 #6245 (18.78D) 38 #6549 (19.25D) 39 #6949 (19.41D) 40 */ 41 42 static double pzero, qzero; 43 static double tpi = .6366197723675813430755350535e0; 44 static double pio4 = .7853981633974483096156608458e0; 45 static double p1[] = { 46 0.4933787251794133561816813446e21, 47 -.1179157629107610536038440800e21, 48 0.6382059341072356562289432465e19, 49 -.1367620353088171386865416609e18, 50 0.1434354939140344111664316553e16, 51 -.8085222034853793871199468171e13, 52 0.2507158285536881945555156435e11, 53 -.4050412371833132706360663322e8, 54 0.2685786856980014981415848441e5, 55 }; 56 static double q1[] = { 57 0.4933787251794133562113278438e21, 58 0.5428918384092285160200195092e19, 59 0.3024635616709462698627330784e17, 60 0.1127756739679798507056031594e15, 61 0.3123043114941213172572469442e12, 62 0.6699987672982239671814028660e9, 63 0.1114636098462985378182402543e7, 64 0.1363063652328970604442810507e4, 65 1.0 66 }; 67 static double p2[] = { 68 0.5393485083869438325262122897e7, 69 0.1233238476817638145232406055e8, 70 0.8413041456550439208464315611e7, 71 0.2016135283049983642487182349e7, 72 0.1539826532623911470917825993e6, 73 0.2485271928957404011288128951e4, 74 0.0, 75 }; 76 static double q2[] = { 77 0.5393485083869438325560444960e7, 78 0.1233831022786324960844856182e8, 79 0.8426449050629797331554404810e7, 80 0.2025066801570134013891035236e7, 81 0.1560017276940030940592769933e6, 82 0.2615700736920839685159081813e4, 83 1.0, 84 }; 85 static double p3[] = { 86 -.3984617357595222463506790588e4, 87 -.1038141698748464093880530341e5, 88 -.8239066313485606568803548860e4, 89 -.2365956170779108192723612816e4, 90 -.2262630641933704113967255053e3, 91 -.4887199395841261531199129300e1, 92 0.0, 93 }; 94 static double q3[] = { 95 0.2550155108860942382983170882e6, 96 0.6667454239319826986004038103e6, 97 0.5332913634216897168722255057e6, 98 0.1560213206679291652539287109e6, 99 0.1570489191515395519392882766e5, 100 0.4087714673983499223402830260e3, 101 1.0, 102 }; 103 static double p4[] = { 104 -.2750286678629109583701933175e20, 105 0.6587473275719554925999402049e20, 106 -.5247065581112764941297350814e19, 107 0.1375624316399344078571335453e18, 108 -.1648605817185729473122082537e16, 109 0.1025520859686394284509167421e14, 110 -.3436371222979040378171030138e11, 111 0.5915213465686889654273830069e8, 112 -.4137035497933148554125235152e5, 113 }; 114 static double q4[] = { 115 0.3726458838986165881989980e21, 116 0.4192417043410839973904769661e19, 117 0.2392883043499781857439356652e17, 118 0.9162038034075185262489147968e14, 119 0.2613065755041081249568482092e12, 120 0.5795122640700729537480087915e9, 121 0.1001702641288906265666651753e7, 122 0.1282452772478993804176329391e4, 123 1.0, 124 }; 125 126 static void 127 asympt(double arg) 128 { 129 double zsq, n, d; 130 int i; 131 132 zsq = 64 / (arg*arg); 133 for(n=0,d=0,i=6;i>=0;i--) { 134 n = n*zsq + p2[i]; 135 d = d*zsq + q2[i]; 136 } 137 pzero = n/d; 138 for(n=0,d=0,i=6;i>=0;i--) { 139 n = n*zsq + p3[i]; 140 d = d*zsq + q3[i]; 141 } 142 qzero = (8/arg)*(n/d); 143 } 144 145 double 146 j0(double arg) 147 { 148 double argsq, n, d; 149 int i; 150 151 if(arg < 0) 152 arg = -arg; 153 if(arg > 8) { 154 asympt(arg); 155 n = arg - pio4; 156 return sqrt(tpi/arg)*(pzero*cos(n) - qzero*sin(n)); 157 } 158 argsq = arg*arg; 159 for(n=0,d=0,i=8;i>=0;i--) { 160 n = n*argsq + p1[i]; 161 d = d*argsq + q1[i]; 162 } 163 return n/d; 164 } 165 166 double 167 y0(double arg) 168 { 169 double argsq, n, d; 170 int i; 171 172 errno = 0; 173 if(arg <= 0) { 174 errno = EDOM; 175 return(-HUGE_VAL); 176 } 177 if(arg > 8) { 178 asympt(arg); 179 n = arg - pio4; 180 return sqrt(tpi/arg)*(pzero*sin(n) + qzero*cos(n)); 181 } 182 argsq = arg*arg; 183 for(n=0,d=0,i=8;i>=0;i--) { 184 n = n*argsq + p4[i]; 185 d = d*argsq + q4[i]; 186 } 187 return n/d + tpi*j0(arg)*log(arg); 188 } 189