2.9BSD - /usr/src/lib/m/j0.c

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

Defined functions

Defined variables