1 #ifndef lint 2 static char sccsid[] = 3 "@(#)j0.c 4.2 (Berkeley) 8/21/85; 5.1 (ucb.elefunt) 11/30/87"; 4 #endif /* not lint */ 5 6 /* 7 floating point Bessel's function 8 of the first and second kinds 9 of order zero 10 11 j0(x) returns the value of J0(x) 12 for all real values of x. 13 14 There are no error returns. 15 Calls sin, cos, sqrt. 16 17 There is a niggling bug in J0 which 18 causes errors up to 2e-16 for x in the 19 interval [-8,8]. 20 The bug is caused by an inappropriate order 21 of summation of the series. rhm will fix it 22 someday. 23 24 Coefficients are from Hart & Cheney. 25 #5849 (19.22D) 26 #6549 (19.25D) 27 #6949 (19.41D) 28 29 y0(x) returns the value of Y0(x) 30 for positive real values of x. 31 For x<=0, if on the VAX, error number EDOM is set and 32 the reserved operand fault is generated; 33 otherwise (an IEEE machine) an invalid operation is performed. 34 35 Calls sin, cos, sqrt, log, j0. 36 37 The values of Y0 have not been checked 38 to more than ten places. 39 40 Coefficients are from Hart & Cheney. 41 #6245 (18.78D) 42 #6549 (19.25D) 43 #6949 (19.41D) 44 */ 45 46 #include <math.h> 47 #if defined(vax)||defined(tahoe) 48 #include <errno.h> 49 #else /* defined(vax)||defined(tahoe) */ 50 static double zero = 0.e0; 51 #endif /* defined(vax)||defined(tahoe) */ 52 static double pzero, qzero; 53 static double tpi = .6366197723675813430755350535e0; 54 static double pio4 = .7853981633974483096156608458e0; 55 static double p1[] = { 56 0.4933787251794133561816813446e21, 57 -.1179157629107610536038440800e21, 58 0.6382059341072356562289432465e19, 59 -.1367620353088171386865416609e18, 60 0.1434354939140344111664316553e16, 61 -.8085222034853793871199468171e13, 62 0.2507158285536881945555156435e11, 63 -.4050412371833132706360663322e8, 64 0.2685786856980014981415848441e5, 65 }; 66 static double q1[] = { 67 0.4933787251794133562113278438e21, 68 0.5428918384092285160200195092e19, 69 0.3024635616709462698627330784e17, 70 0.1127756739679798507056031594e15, 71 0.3123043114941213172572469442e12, 72 0.6699987672982239671814028660e9, 73 0.1114636098462985378182402543e7, 74 0.1363063652328970604442810507e4, 75 1.0 76 }; 77 static double p2[] = { 78 0.5393485083869438325262122897e7, 79 0.1233238476817638145232406055e8, 80 0.8413041456550439208464315611e7, 81 0.2016135283049983642487182349e7, 82 0.1539826532623911470917825993e6, 83 0.2485271928957404011288128951e4, 84 0.0, 85 }; 86 static double q2[] = { 87 0.5393485083869438325560444960e7, 88 0.1233831022786324960844856182e8, 89 0.8426449050629797331554404810e7, 90 0.2025066801570134013891035236e7, 91 0.1560017276940030940592769933e6, 92 0.2615700736920839685159081813e4, 93 1.0, 94 }; 95 static double p3[] = { 96 -.3984617357595222463506790588e4, 97 -.1038141698748464093880530341e5, 98 -.8239066313485606568803548860e4, 99 -.2365956170779108192723612816e4, 100 -.2262630641933704113967255053e3, 101 -.4887199395841261531199129300e1, 102 0.0, 103 }; 104 static double q3[] = { 105 0.2550155108860942382983170882e6, 106 0.6667454239319826986004038103e6, 107 0.5332913634216897168722255057e6, 108 0.1560213206679291652539287109e6, 109 0.1570489191515395519392882766e5, 110 0.4087714673983499223402830260e3, 111 1.0, 112 }; 113 static double p4[] = { 114 -.2750286678629109583701933175e20, 115 0.6587473275719554925999402049e20, 116 -.5247065581112764941297350814e19, 117 0.1375624316399344078571335453e18, 118 -.1648605817185729473122082537e16, 119 0.1025520859686394284509167421e14, 120 -.3436371222979040378171030138e11, 121 0.5915213465686889654273830069e8, 122 -.4137035497933148554125235152e5, 123 }; 124 static double q4[] = { 125 0.3726458838986165881989980e21, 126 0.4192417043410839973904769661e19, 127 0.2392883043499781857439356652e17, 128 0.9162038034075185262489147968e14, 129 0.2613065755041081249568482092e12, 130 0.5795122640700729537480087915e9, 131 0.1001702641288906265666651753e7, 132 0.1282452772478993804176329391e4, 133 1.0, 134 }; 135 136 double 137 j0(arg) double arg;{ 138 double argsq, n, d; 139 double sin(), cos(), sqrt(); 140 int i; 141 142 if(arg < 0.) arg = -arg; 143 if(arg > 8.){ 144 asympt(arg); 145 n = arg - pio4; 146 return(sqrt(tpi/arg)*(pzero*cos(n) - qzero*sin(n))); 147 } 148 argsq = arg*arg; 149 for(n=0,d=0,i=8;i>=0;i--){ 150 n = n*argsq + p1[i]; 151 d = d*argsq + q1[i]; 152 } 153 return(n/d); 154 } 155 156 double 157 y0(arg) double arg;{ 158 double argsq, n, d; 159 double sin(), cos(), sqrt(), log(), j0(); 160 int i; 161 162 if(arg <= 0.){ 163 #if defined(vax)||defined(tahoe) 164 extern double infnan(); 165 return(infnan(EDOM)); /* NaN */ 166 #else /* defined(vax)||defined(tahoe) */ 167 return(zero/zero); /* IEEE machines: invalid operation */ 168 #endif /* defined(vax)||defined(tahoe) */ 169 } 170 if(arg > 8.){ 171 asympt(arg); 172 n = arg - pio4; 173 return(sqrt(tpi/arg)*(pzero*sin(n) + qzero*cos(n))); 174 } 175 argsq = arg*arg; 176 for(n=0,d=0,i=8;i>=0;i--){ 177 n = n*argsq + p4[i]; 178 d = d*argsq + q4[i]; 179 } 180 return(n/d + tpi*j0(arg)*log(arg)); 181 } 182 183 static 184 asympt(arg) double arg;{ 185 double zsq, n, d; 186 int i; 187 zsq = 64./(arg*arg); 188 for(n=0,d=0,i=6;i>=0;i--){ 189 n = n*zsq + p2[i]; 190 d = d*zsq + q2[i]; 191 } 192 pzero = n/d; 193 for(n=0,d=0,i=6;i>=0;i--){ 194 n = n*zsq + p3[i]; 195 d = d*zsq + q3[i]; 196 } 197 qzero = (8./arg)*(n/d); 198 } 199