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