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