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