1 #ifndef lint 2 static char sccsid[] = 3 "@(#)j1.c 4.2 (Berkeley) 8/21/85; 1.2 (ucb.elefunt) 09/11/85"; 4 #endif not lint 5 6 /* 7 floating point Bessel's function 8 of the first and second kinds 9 of order one 10 11 j1(x) returns the value of J1(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 J1 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 #6050 (20.98D) 26 #6750 (19.19D) 27 #7150 (19.35D) 28 29 y1(x) returns the value of Y1(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, j1. 36 37 The values of Y1 have not been checked 38 to more than ten places. 39 40 Coefficients are from Hart & Cheney. 41 #6447 (22.18D) 42 #6750 (19.19D) 43 #7150 (19.35D) 44 */ 45 46 #include <math.h> 47 #ifdef VAX 48 #include <errno.h> 49 #else /* IEEE double */ 50 static double zero = 0.e0; 51 #endif 52 static double pzero, qzero; 53 static double tpi = .6366197723675813430755350535e0; 54 static double pio4 = .7853981633974483096156608458e0; 55 static double p1[] = { 56 0.581199354001606143928050809e21, 57 -.6672106568924916298020941484e20, 58 0.2316433580634002297931815435e19, 59 -.3588817569910106050743641413e17, 60 0.2908795263834775409737601689e15, 61 -.1322983480332126453125473247e13, 62 0.3413234182301700539091292655e10, 63 -.4695753530642995859767162166e7, 64 0.2701122710892323414856790990e4, 65 }; 66 static double q1[] = { 67 0.1162398708003212287858529400e22, 68 0.1185770712190320999837113348e20, 69 0.6092061398917521746105196863e17, 70 0.2081661221307607351240184229e15, 71 0.5243710262167649715406728642e12, 72 0.1013863514358673989967045588e10, 73 0.1501793594998585505921097578e7, 74 0.1606931573481487801970916749e4, 75 1.0, 76 }; 77 static double p2[] = { 78 -.4435757816794127857114720794e7, 79 -.9942246505077641195658377899e7, 80 -.6603373248364939109255245434e7, 81 -.1523529351181137383255105722e7, 82 -.1098240554345934672737413139e6, 83 -.1611616644324610116477412898e4, 84 0.0, 85 }; 86 static double q2[] = { 87 -.4435757816794127856828016962e7, 88 -.9934124389934585658967556309e7, 89 -.6585339479723087072826915069e7, 90 -.1511809506634160881644546358e7, 91 -.1072638599110382011903063867e6, 92 -.1455009440190496182453565068e4, 93 1.0, 94 }; 95 static double p3[] = { 96 0.3322091340985722351859704442e5, 97 0.8514516067533570196555001171e5, 98 0.6617883658127083517939992166e5, 99 0.1849426287322386679652009819e5, 100 0.1706375429020768002061283546e4, 101 0.3526513384663603218592175580e2, 102 0.0, 103 }; 104 static double q3[] = { 105 0.7087128194102874357377502472e6, 106 0.1819458042243997298924553839e7, 107 0.1419460669603720892855755253e7, 108 0.4002944358226697511708610813e6, 109 0.3789022974577220264142952256e5, 110 0.8638367769604990967475517183e3, 111 1.0, 112 }; 113 static double p4[] = { 114 -.9963753424306922225996744354e23, 115 0.2655473831434854326894248968e23, 116 -.1212297555414509577913561535e22, 117 0.2193107339917797592111427556e20, 118 -.1965887462722140658820322248e18, 119 0.9569930239921683481121552788e15, 120 -.2580681702194450950541426399e13, 121 0.3639488548124002058278999428e10, 122 -.2108847540133123652824139923e7, 123 0.0, 124 }; 125 static double q4[] = { 126 0.5082067366941243245314424152e24, 127 0.5435310377188854170800653097e22, 128 0.2954987935897148674290758119e20, 129 0.1082258259408819552553850180e18, 130 0.2976632125647276729292742282e15, 131 0.6465340881265275571961681500e12, 132 0.1128686837169442121732366891e10, 133 0.1563282754899580604737366452e7, 134 0.1612361029677000859332072312e4, 135 1.0, 136 }; 137 138 double 139 j1(arg) double arg;{ 140 double xsq, n, d, x; 141 double sin(), cos(), sqrt(); 142 int i; 143 144 x = arg; 145 if(x < 0.) x = -x; 146 if(x > 8.){ 147 asympt(x); 148 n = x - 3.*pio4; 149 n = sqrt(tpi/x)*(pzero*cos(n) - qzero*sin(n)); 150 if(arg <0.) n = -n; 151 return(n); 152 } 153 xsq = x*x; 154 for(n=0,d=0,i=8;i>=0;i--){ 155 n = n*xsq + p1[i]; 156 d = d*xsq + q1[i]; 157 } 158 return(arg*n/d); 159 } 160 161 double 162 y1(arg) double arg;{ 163 double xsq, n, d, x; 164 double sin(), cos(), sqrt(), log(), j1(); 165 int i; 166 167 x = arg; 168 if(x <= 0.){ 169 #ifdef VAX 170 extern double infnan(); 171 return(infnan(EDOM)); /* NaN */ 172 #else /* IEEE double */ 173 return(zero/zero); /* IEEE machines: invalid operation */ 174 #endif 175 } 176 if(x > 8.){ 177 asympt(x); 178 n = x - 3*pio4; 179 return(sqrt(tpi/x)*(pzero*sin(n) + qzero*cos(n))); 180 } 181 xsq = x*x; 182 for(n=0,d=0,i=9;i>=0;i--){ 183 n = n*xsq + p4[i]; 184 d = d*xsq + q4[i]; 185 } 186 return(x*n/d + tpi*(j1(x)*log(x)-1./x)); 187 } 188 189 static 190 asympt(arg) double arg;{ 191 double zsq, n, d; 192 int i; 193 zsq = 64./(arg*arg); 194 for(n=0,d=0,i=6;i>=0;i--){ 195 n = n*zsq + p2[i]; 196 d = d*zsq + q2[i]; 197 } 198 pzero = n/d; 199 for(n=0,d=0,i=6;i>=0;i--){ 200 n = n*zsq + p3[i]; 201 d = d*zsq + q3[i]; 202 } 203 qzero = (8./arg)*(n/d); 204 } 205