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