471 lines
13 KiB
C
471 lines
13 KiB
C
|
|
/* Copyright (c) Mix Software 1988 */
|
|
/* ------------------------------------------------------------ */
|
|
/* bessel functions: j0, j1, jn */
|
|
/* y0, y1, yn */
|
|
double yn(n,x)
|
|
int n;
|
|
double x;
|
|
{
|
|
double y0(), y1();
|
|
double yz, yx, y, fx;
|
|
double _domerr();
|
|
extern double HUGE;
|
|
int negative = 0;
|
|
int i;
|
|
if (n == 0) return y0(x);
|
|
if (n == 1) return y1(x);
|
|
if (n < 0) {
|
|
n = -n;
|
|
if (n & 0x0001) negative = 1;
|
|
}
|
|
if (x <= 0.0) {_domerr("yn",x,x); return HUGE;}
|
|
yz = y0(x);
|
|
yx = y1(x);
|
|
fx = 2.0/x;
|
|
i = 1;
|
|
while (i < n) {
|
|
y = fx * i * yx - yz;
|
|
yz = yx;
|
|
yx = y;
|
|
i++;
|
|
}
|
|
if (negative) return -y;
|
|
return y;
|
|
} /* yn */
|
|
|
|
double jn(n,x)
|
|
int n;
|
|
double x;
|
|
{
|
|
double j0(), j1();
|
|
double jz, jx, j, fx;
|
|
int i;
|
|
if (n == 0) return j0(x);
|
|
if (n == 1) return j1(x);
|
|
if (n < 0) return jn(-n,-x);
|
|
jz = j0(x);
|
|
jx = j1(x);
|
|
fx = 2.0/x;
|
|
i = 1;
|
|
while (i < n) {
|
|
j = fx * i * jx - jz;
|
|
jz = jx;
|
|
jx = j;
|
|
i++;
|
|
}
|
|
return j;
|
|
} /* jn */
|
|
|
|
double y0(x)
|
|
double x;
|
|
{
|
|
double xsquare, xx, x8, y0;
|
|
double p, q;
|
|
double sin(), cos(), sqrt(), log(), j0();
|
|
double _domerr();
|
|
extern double HUGE;
|
|
void _pandq0();
|
|
static double pi = 3.1415926535897932;
|
|
static double piby4 = 0.78539816339744831;
|
|
static double twobypi = 0.63661977236758134;
|
|
static struct {
|
|
int terms;
|
|
double tbl[8];
|
|
} y0P = {8,
|
|
-0.40421653126104926e19,
|
|
0.96736537146262430e19,
|
|
-0.75173702582406051e18,
|
|
0.18702004643462448e17,
|
|
-0.20430408868635631e15,
|
|
0.10869647137874908e13,
|
|
-0.27725001364628094e10,
|
|
0.27295402301235549e7
|
|
};
|
|
static struct {
|
|
int terms;
|
|
double tbl[9];
|
|
} y0Q = {9,
|
|
0.54768700204477752e20,
|
|
0.72620404690767952e18,
|
|
0.49441312329373299e16,
|
|
0.22900727639481849e14,
|
|
0.80363630456265291e11,
|
|
0.22400416108570433e9,
|
|
0.50069757439413062e6,
|
|
0.86456632023518903e3,
|
|
0.10000000000000000e1
|
|
};
|
|
|
|
extern double (far *$0DPOWER)();
|
|
if (x <= 0.0) {_domerr("y0",x,x); return HUGE;}
|
|
if (x < 8.0) {
|
|
xsquare = x*x;
|
|
y0 = $0DPOWER(xsquare,&y0P)/$0DPOWER(xsquare,&y0Q);
|
|
return y0 + twobypi*j0(x)*log(x);
|
|
}
|
|
xx = x - piby4;
|
|
x8 = 8.0/x;
|
|
_pandq0(x8,&p,&q);
|
|
return sqrt(2.0/(pi*x)) * (p*sin(xx) + q*cos(xx));
|
|
} /* y0 */
|
|
|
|
double y1(x)
|
|
double x;
|
|
{
|
|
double xsquare, xx, x8, y1;
|
|
double p, q;
|
|
double sin(), cos(), sqrt(), log(), j1();
|
|
double _domerr();
|
|
extern double HUGE;
|
|
void _pandq1();
|
|
static double pi = 3.1415926535897932;
|
|
static double piby34 = 2.35619449019234493;
|
|
static double twobypi = 0.63661977236758134;
|
|
static struct {
|
|
int terms;
|
|
double tbl[8];
|
|
} y1P = {8,
|
|
-0.29238219615329625e20,
|
|
0.77485206821868396e19,
|
|
-0.34410480630841144e18,
|
|
0.59151607604900706e16,
|
|
-0.48633169425671751e14,
|
|
0.20496966737456622e12,
|
|
-0.42894719688552488e9,
|
|
0.35569240098305261e6
|
|
};
|
|
static struct {
|
|
int terms;
|
|
double tbl[9];
|
|
} y1Q = {9,
|
|
0.14913115113029204e21,
|
|
0.18186628417061350e19,
|
|
0.11316393826988845e17,
|
|
0.47551735888881377e14,
|
|
0.15002216991567090e12,
|
|
0.37166607986219303e9,
|
|
0.72691473071988846e6,
|
|
0.10726961437789255e4,
|
|
0.10000000000000000e1
|
|
};
|
|
|
|
extern double (far *$0DPOWER)();
|
|
if (x <= 0.0) {_domerr("y1",x,x); return HUGE;}
|
|
if (x < 8.0) {
|
|
xsquare = x*x;
|
|
y1 = x*$0DPOWER(xsquare,&y1P)/$0DPOWER(xsquare,&y1Q);
|
|
return y1 + twobypi*(j1(x)*log(x)-1.0/x);
|
|
}
|
|
xx = x - piby34;
|
|
x8 = 8.0/x;
|
|
_pandq1(x8,&p,&q);
|
|
return sqrt(2.0/(pi*x)) * (p*sin(xx) + q*cos(xx));
|
|
} /* y1 */
|
|
|
|
double j0(x)
|
|
double x;
|
|
{
|
|
double xsquare, xx, x8;
|
|
double p, q;
|
|
double sin(), cos(), sqrt(), fabs();
|
|
void _pandq0();
|
|
static double pi = 3.1415926535897932;
|
|
static double piby4 = 0.78539816339744831;
|
|
static struct {
|
|
int terms;
|
|
double tbl[12];
|
|
} j0P = {12,
|
|
0.18108389921416431e9,
|
|
-0.44464759561502412e8,
|
|
0.26292976097233912e7,
|
|
-0.66351202963350577e5,
|
|
0.90000108663859700e3,
|
|
-0.74117799951193489e1,
|
|
0.39771954892980604e-1,
|
|
-0.14463518634545936e-3,
|
|
0.36127859495623837e-6,
|
|
-0.60893587508087894e-9,
|
|
0.64247367849791009e-12,
|
|
-0.33133602728361709e-15
|
|
};
|
|
static struct {
|
|
int terms;
|
|
double tbl[4];
|
|
} j0Q = {4,
|
|
0.18108389921416430e9,
|
|
0.80621524203869792e6,
|
|
0.14154950117074173e4,
|
|
0.10000000000000000e1
|
|
};
|
|
static struct {
|
|
int terms;
|
|
double tbl[5];
|
|
} j0pP = {5,
|
|
0.21807736478830516e7,
|
|
0.30343163608475270e7,
|
|
0.11035444210405852e7,
|
|
0.11252515255664381e6,
|
|
0.22399036697750965e4
|
|
};
|
|
static struct {
|
|
int terms;
|
|
double tbl[6];
|
|
} j0pQ = {6,
|
|
0.21807736478830516e7,
|
|
0.30367122303337213e7,
|
|
0.11068209412295708e7,
|
|
0.11366275712613906e6,
|
|
0.23403140106394541e4,
|
|
0.10000000000000000e1
|
|
};
|
|
static struct {
|
|
int terms;
|
|
double tbl[5];
|
|
} j0qP = {5,
|
|
-0.17665456233082465e4,
|
|
-0.28720316121456664e4,
|
|
-0.12598828601325539e4,
|
|
-0.16263421062270593e3,
|
|
-0.44237584856933353e1
|
|
};
|
|
static struct {
|
|
int terms;
|
|
double tbl[6];
|
|
} j0qQ = {6,
|
|
0.11305891989633582e6,
|
|
0.18484510850351025e6,
|
|
0.82274660980144657e5,
|
|
0.11085805836751487e5,
|
|
0.35655140058576331e3,
|
|
0.10000000000000000e1
|
|
};
|
|
|
|
extern double (far *$0DPOWER)();
|
|
x = fabs(x);
|
|
if (x < 8.0) {
|
|
xsquare = x*x;
|
|
return $0DPOWER(xsquare,&j0P)/$0DPOWER(xsquare,&j0Q);
|
|
}
|
|
xx = x - piby4;
|
|
x8 = 8.0/x;
|
|
_pandq0(x8,&p,&q);
|
|
/*
|
|
xsquare = x8*x8;
|
|
p = $0DPOWER(xsquare,&j0pP)/$0DPOWER(xsquare,&j0pQ);
|
|
q = x8*$0DPOWER(xsquare,&j0qP)/$0DPOWER(xsquare,&j0qQ);
|
|
*/
|
|
return sqrt(2.0/(pi*x)) * (p*cos(xx) - q*sin(xx));
|
|
} /* j0 */
|
|
|
|
double j1(x)
|
|
double x;
|
|
{
|
|
double xsquare, xx, x8, value;
|
|
double p, q;
|
|
int negative = 0;
|
|
double sin(), cos(), sqrt();
|
|
void _pandq1();
|
|
static double pi = 3.1415926535897932;
|
|
static double piby34 = 2.35619449019234493;
|
|
extern double (far *$0DPOWER)();
|
|
static struct {
|
|
int terms;
|
|
double tbl[8];
|
|
} j1P = {8,
|
|
0.22148878804219631e18,
|
|
-0.25123742147032128e17,
|
|
0.84824207447812727e15,
|
|
-0.12498203672620249e14,
|
|
0.93165652967246732e11,
|
|
-0.36866689870229816e9,
|
|
0.74370238171199964e6,
|
|
-0.60795301796074136e3
|
|
};
|
|
static struct {
|
|
int terms;
|
|
double tbl[8];
|
|
} j1Q = {8,
|
|
0.44297757608439262e18,
|
|
0.51247127164848721e16,
|
|
0.29898363077254872e14,
|
|
0.11581921274668893e12,
|
|
0.32819403445341964e9,
|
|
0.69885861844850757e6,
|
|
0.10777412894333044e4,
|
|
0.10000000000000000e1
|
|
};
|
|
/* static struct {
|
|
int terms;
|
|
double tbl[5];
|
|
} j1pP = {5,
|
|
-0.17469576694409286e7,
|
|
-0.23669439140521428e7,
|
|
-0.82850363738723775e6,
|
|
-0.79574568713505959e5,
|
|
-0.14279066288270981e4
|
|
};
|
|
static struct {
|
|
int terms;
|
|
double tbl[6];
|
|
} j1pQ = {6,
|
|
-0.17469576694409286e7,
|
|
-0.23637451390226540e7,
|
|
-0.82423699065628189e6,
|
|
-0.78144050089391111e5,
|
|
-0.13084529833390797e4,
|
|
0.10000000000000000e1
|
|
};
|
|
static struct {
|
|
int terms;
|
|
double tbl[5];
|
|
} j1qP = {5,
|
|
0.14465282874995209e3,
|
|
0.17442916890924259e3,
|
|
0.51736532818365916e2,
|
|
0.37994453796980673e1,
|
|
0.36363466476034711e-1
|
|
};
|
|
static struct {
|
|
int terms;
|
|
double tbl[5];
|
|
} j1qQ = {5,
|
|
0.30859270133323172e4,
|
|
0.37343401060163018e4,
|
|
0.11191098527047487e4,
|
|
0.85223920643413404e2,
|
|
0.10000000000000000e1
|
|
};
|
|
*/
|
|
|
|
if (x < 0.0) {
|
|
negative = 1;
|
|
x = -x;
|
|
}
|
|
if (x < 8.0) {
|
|
xsquare = x*x;
|
|
value = x*$0DPOWER(xsquare,&j1P)/$0DPOWER(xsquare,&j1Q);
|
|
if (negative) return -value; else return value;
|
|
}
|
|
xx = x - piby34;
|
|
x8 = 8.0/x;
|
|
/*
|
|
xsquare = x8*x8;
|
|
p = $0DPOWER(xsquare,&j1pP)/$0DPOWER(xsquare,&j1pQ);
|
|
q = x8*$0DPOWER(xsquare,&j1qP)/$0DPOWER(xsquare,&j1qQ);
|
|
*/
|
|
_pandq1(x8,&p,&q);
|
|
value = sqrt(2.0/(pi*x)) * (p*cos(xx) - q*sin(xx));
|
|
if (negative) return -value; else return value;
|
|
} /* j1 */
|
|
|
|
void _pandq0(x, pptr, qptr)
|
|
double x;
|
|
double *pptr;
|
|
double *qptr;
|
|
{
|
|
double xsquare;
|
|
static struct {
|
|
int terms;
|
|
double tbl[5];
|
|
} P0p = {5,
|
|
0.21807736478830516e7,
|
|
0.30343163608475270e7,
|
|
0.11035444210405852e7,
|
|
0.11252515255664381e6,
|
|
0.22399036697750965e4
|
|
};
|
|
static struct {
|
|
int terms;
|
|
double tbl[6];
|
|
} P0q = {6,
|
|
0.21807736478830516e7,
|
|
0.30367122303337213e7,
|
|
0.11068209412295708e7,
|
|
0.11366275712613906e6,
|
|
0.23403140106394541e4,
|
|
0.10000000000000000e1
|
|
};
|
|
static struct {
|
|
int terms;
|
|
double tbl[5];
|
|
} Q0p = {5,
|
|
-0.17665456233082465e4,
|
|
-0.28720316121456664e4,
|
|
-0.12598828601325539e4,
|
|
-0.16263421062270593e3,
|
|
-0.44237584856933353e1
|
|
};
|
|
static struct {
|
|
int terms;
|
|
double tbl[6];
|
|
} Q0q = {6,
|
|
0.11305891989633582e6,
|
|
0.18484510850351025e6,
|
|
0.82274660980144657e5,
|
|
0.11085805836751487e5,
|
|
0.35655140058576331e3,
|
|
0.10000000000000000e1
|
|
};
|
|
extern double (far *$0DPOWER)();
|
|
|
|
xsquare = x*x;
|
|
*pptr = $0DPOWER(xsquare,&P0p)/$0DPOWER(xsquare,&P0q);
|
|
*qptr = x*$0DPOWER(xsquare,&Q0p)/$0DPOWER(xsquare,&Q0q);
|
|
}
|
|
|
|
void _pandq1(x, pptr, qptr)
|
|
double x;
|
|
double *pptr;
|
|
double *qptr;
|
|
{
|
|
double xsquare;
|
|
static struct {
|
|
int terms;
|
|
double tbl[5];
|
|
} P1p = {5,
|
|
-0.17469576694409286e7,
|
|
-0.23669439140521428e7,
|
|
-0.82850363738723775e6,
|
|
-0.79574568713505959e5,
|
|
-0.14279066288270981e4
|
|
};
|
|
static struct {
|
|
int terms;
|
|
double tbl[6];
|
|
} P1q = {6,
|
|
-0.17469576694409286e7,
|
|
-0.23637451390226540e7,
|
|
-0.82423699065628189e6,
|
|
-0.78144050089391111e5,
|
|
-0.13084529833390797e4,
|
|
0.10000000000000000e1
|
|
};
|
|
static struct {
|
|
int terms;
|
|
double tbl[5];
|
|
} Q1p = {5,
|
|
0.14465282874995209e3,
|
|
0.17442916890924259e3,
|
|
0.51736532818365916e2,
|
|
0.37994453796980673e1,
|
|
0.36363466476034711e-1
|
|
};
|
|
static struct {
|
|
int terms;
|
|
double tbl[5];
|
|
} Q1q = {5,
|
|
0.30859270133323172e4,
|
|
0.37343401060163018e4,
|
|
0.11191098527047487e4,
|
|
0.85223920643413404e2,
|
|
0.10000000000000000e1
|
|
};
|
|
|
|
extern double (far *$0DPOWER)();
|
|
|
|
xsquare = x*x;
|
|
*pptr = $0DPOWER(xsquare,&P1p)/$0DPOWER(xsquare,&P1q);
|
|
*qptr = x*$0DPOWER(xsquare,&Q1p)/$0DPOWER(xsquare,&Q1q);
|
|
}
|