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