4.3BSD - /usr/src/usr.lib/libom/j0.c

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

Defined functions

Defined variables