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