BESS.C
/****************************************************************************
(c) 1984 - 2003 by Scientific Endeavors Corporation.
All rights reserved.
This program plots circular membrane vibration modes.
Illustrates a 2-D projection with a grid.
****************************************************************************/
#include <graphic.h> /* Include all needed files */
#if 1-MSQ
#ifndef fsign
#define fsign(x) (((x) < (float)0.) ? (float)(-1.) : (float)1.)
#endif
#endif
#if TCQ
extern unsigned stklen=6144; /* Set stack for Borland (Turbo) C */
#endif
float eps = .4f;
char *larray[ ] = {""};
int lnstyle[ ] = {
211, 221, 231, 241, 251, 211,
221, 231, 241, 251, 211
};
int lnstyle2[ ] = {
211, 221, 231, 241, 201, 211, 221, 231, 241, 201,
211, 221, 231, 241, 201, 211, 221, 231, 241, 201,
211, 221, 231, 241, 201, 211, 221, 231, 241, 201,
211, 221, 231, 241, 201, 211, 221, 231, 241, 201,
211, 221, 231, 241, 201, 211, 221, 231, 241, 201
};
float x[121], y[121], level[1] = {1.15};
/*int clrs[] = {447,417,414,412,418,421,427,442,445,447};*/
int clrs[] = {406,407,408,409,410,411,375,339,303,267,261,255};
float W_CDECL zfun1(float x, float y);
float W_CDECL zfun2(float r, float theta);
#if 1 - (MSQ | IBMQ)
double jn(int n, double x);
double yn(int n, double x);
double bessj0q(double x);
double bessj1q(double x);
double bessy0q(double x);
double bessy1q(double x);
#endif
/****************************************************************************
Main program
****************************************************************************/
void GPC_MAIN(void)
{
float x3axis, y3axis, z3axis, xvu, yvu, zvu, xorig, xstp, xmax;
float yorig, ystp, ymax, zorig, zstep, zmax;
float xdel, ydel;
int i, ixpts, iypts;
ixpts = 1; /* Two segments across cell */
iypts = 1;
x3axis = 5.5f; /* Define workbox area */
y3axis = 5.5f;
z3axis = 5.0f;
xvu = -48.f; /* Position of view point in workbox units */
yvu = 48.f;
zvu = 15.f;
xorig = -1.0f; /* Set x-axis limits */
xstp = .5f;
xdel = .1f;
xmax = 1.0f;
yorig = -1.0f; /* Set y-axis limits */
ystp = 0.5f;
ydel = .1f;
ymax = 1.0f;
zorig = -1.f; /* Set z-axis limits */
zstep = .5f;
zmax = .5f;
for(i = 0; i < 81; i++) {
x[i] = .0125f * i;
y[i] = (float)PI * i / 40.f;
}
bgnplot(1, 'g', "bessel.tkf"); /* GraphiC initialization */
startplot(BRT_WHITE);
metricunits(0); /* Plots are in inches */
font(3, "swiss.fnt", '\310', "swissbld.fnt", '\311',
"newsgrm.fnt", '\312' ); /* Select fonts */
fillfont(1);
page(9.0f, 6.884f); /* Size of page and plot area */
area2d(7.4f, 6.0f);
plotabsolute(1);
volm3d(x3axis,y3axis,z3axis); /* Set workbox size */
vuabs(xvu,yvu,zvu); /* Specify viewpoint */
axnamht(.25f);
color(BLACK); /* Invisible 3-D axes */
numht(.18f);
tickout(1); /* Ticks face outward */
polar(1);
axesoff(GPC_AXESOFF);
graf3d(xorig, xstp, xmax, yorig, ystp, ymax, zorig, zstep, zmax);
d3color(clrs, 12, GPC_DRAW_MESH);
survis("both"); /* Both top and bottom surfaces seen */
surfun(zfun1, ixpts, xdel, iypts, ydel); /* Plot surface */
tline(8);
tcurve(8);
plane3d(1, zorig, 0); /* Plot 2-D projection */
polgrid(1.f, 0.f, 1,larray, 0, 0, 0);
contour(zfun2, (float *)x, 81, (float *)y,
81,-0.4f, 0.2f, 1.0f, lnstyle2, 0);
tline(3);
color(BLACK); /* Draw grids in back of plot */
axesoff(GPC_AXESON);
grid(1);
plane3d(3,xmax,1); /* y-z plane */
y3name("\310amplitude"); /* Vertical axis name */
x3name("\310y"); /* Horizontal axis name */
d3grid(1, 1);
plane3d(2,yorig,1); /* x-z plane */
x3name("\310x"); /* Horizontal axis name */
d3grid(1, 1);
plane3d(0,0.0f,0); /* Restore plane to 2-D */
tline(1);
titlht(.2f); /* Plot title */
color(RED);
skew(.1f);
d3head("\311Bessel Function Plot: \310 `J]21[cos(2\312q\310)",'b');
endplot(); /* Plot closing routines */
stopplot();
}
/****************************************************************************
Define function for 3-D surface plot
****************************************************************************/
float W_CDECL zfun1(float xx, float yy)
{
double r;
double theta;
float fn;
r = sqrt(xx * xx + yy * yy);
if (r >= 1.0)
fn = 0.0;
else {
if (fabs(xx) < 1.0e-20)
theta = PI / 2.0;
else
theta = atan(yy / xx);
fn = (float)jn(2, 5.1356 * r) * cos(2. * theta);
}
return (fn);
}
/****************************************************************************
Define function for contour plot
****************************************************************************/
float W_CDECL zfun2(float r, float theta)
{
float fn;
fn = jn(2, 5.1356 * r) * cos(2. * (double)theta);
return (fn);
}
#if 1 - (MSQ | IBMQ)
/********************************************************************
This following code contains the routines for implementing
the bessel functions.
********************************************************************/
/********************************************************************
FUNCTION NAME:
ARGUMENTS:
PURPOSE:
********************************************************************/
double jn(int n, double x)
{
int j, m, jsum;
int iacc = 40;
double bigno = 1.e10;
double bigni = 1.e-10;
double bjm, bjp, bj, tox, sum, bessj;
if(n == 0)
return(bessj0q(x));
if(n == 1)
return(bessj1q(x));
if(x == 0.)
return(0.);
tox = 2. / x;
if(x > (double)n){
bjm = bessj0q(x);
bj = bessj1q(x);
for(j = 1; j < n; j++){
bjp = j * tox * bj - bjm;
bjm = bj;
bj = bjp;
}
return(bj);
}
else{
m = 2 * ((n + (int)(sqrt((double)(iacc * n)))) / 2);
bessj = 0.;
jsum = 0;
sum = 0.;
bjp = 0.;
bj = 1.;
for(j = m; j > 0; j--){
bjm = j * tox * bj - bjp;
bjp = bj;
bj = bjm;
if(fabs(bj) > bigno){
bj *= bigni;
bjp *= bigni;
bessj *= bigni;
sum *= bigni;
}
if(jsum != 0)
sum += bj;
jsum = 1 - jsum;
if(j == n)
bessj = bjp;
}
sum = 2. * sum - bj;
return(bessj / sum);
}
}
/********************************************************************
FUNCTION NAME:
ARGUMENTS:
PURPOSE:
********************************************************************/
double yn(int n,double x)
{
int j;
double by, bym, byp, tox;
if(n == 0)
return(bessy0q(x));
if (n == 1)
return(bessy1q(x));
if(x == 0.)
return(0.);
tox = 2. / x;
by = bessy1q(x);
bym = bessy0q(x);
for(j = 1; j < n; j++){
byp = j * tox * by - bym;
bym = by;
by = byp;
}
return(by);
}
/********************************************************************
FUNCTION NAME:
ARGUMENTS:
PURPOSE:
********************************************************************/
double bessj0q(double x)
{
double y, ax, xx, z;
double p1 = 1.0;
double p2 = -.1098628627e-2;
double p3 = .2734510407e-4;
double p4 = -.2073370639e-5;
double p5 = .2093887211e-6;
double q1 = -.1562499995e-1;
double q2 = .1430488765e-3;
double q3 = -.6911147651e-5;
double q4 = .7621095161e-6;
double q5 = -.934945152e-7;
double r1 = 57568490574.;
double r2 = -13362590354.;
double r3 = 651619640.7;
double r4 = -11214424.18;
double r5 = 77392.33017;
double r6 = -184.9052456;
double s1 = 57568490411.;
double s2 = 1029532985.;
double s3 = 9494680.718;
double s4 = 59272.64853;
double s5 = 267.8532712;
double s6 = 1.;
if(fabs(x) < 8.){
y = x * x;
return((r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * r6))))) /
(s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6))))));
}
else{
ax = fabs(x);
z = 8. / ax;
y = z * z;
xx = ax - .785398164;
return(sqrt(.636619772 / ax) *
(cos(xx) * (p1 + y * (p2 + y * (p3 + y * (p4 + y * p5)))) -
z * sin(xx) * (q1 + y * (q2 + y * (q3 + y * (q4 * y + q5))))));
}
}
/********************************************************************
FUNCTION NAME:
ARGUMENTS:
PURPOSE:
********************************************************************/
double bessj1q(double x)
{
double y, xx, ax, z;
double r1 = 72362614232.;
double r2 = -7895059235.;
double r3 = 242396853.1;
double r4 = -2972611.439;
double r5 = 15704.48260;
double r6 = -30.16036606;
double s1 = 144725228442.;
double s2 = 2300535178.;
double s3 = 18583304.74;
double s4 = 99447.43394;
double s5 = 376.9991397;
double s6 = 1.;
double p1 = 1.;
double p2 = .183105e-2;
double p3 = -.3516396496e-4;
double p4 = .2457520174e-5;
double p5 = -.240337019e-6;
double q1 = .04687499995;
double q2 = -.2002690873e-3;
double q3 = .8449199096e-5;
double q4 = -.88228987e-6;
double q5 = .105787412e-6;
if(fabs(x) < 8.){
y = x * x;
return(x * (r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * r6)))))/
(s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6))))));
}
else {
ax = fabs(x);
z = 8. / ax;
y = z * z;
xx = ax - 2.356194491;
return(sqrt(.636619772 / ax) *
(cos(xx) * (p1 + y * (p2 + y * (p3 + y * (p4 + y * p5)))) -
z * sin(xx) * (q1 + y * (q2 + y * (q3 + y * (q4 + y * q5))))) *
fsign(x));
}
}
/********************************************************************
FUNCTION NAME:
ARGUMENTS:
PURPOSE:
********************************************************************/
double bessy0q(double x)
{
double y, z, xx;
double p1 = 1.;
double p2 = -.1098628627e-2;
double p3 = .2734510407e-4;
double p4 = -.2073370639e-5;
double p5 = .2093887211e-6;
double q1 = -.1562499995e-1;
double q2 = .1430488765e-3;
double q3 = -.6911147651e-5;
double q4 = .7621095161e-6;
double q5 = -.934945152e-7;
double r1 = -2957821389.;
double r2 = 7062834065.;
double r3 = -512359803.6;
double r4 = 10879881.29;
double r5 = -86327.92757;
double r6 = 228.4622733;
double s1 = 40076544269.;
double s2 = 745249964.8;
double s3 = 7189466.438;
double s4 = 47447.26470;
double s5 = 226.1030244;
double s6 = 1.;
if(x < 8.){
y = x * x;
return((r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * r6))))) /
(s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6))))) +
.63661972 * bessj0q(x) * log(x));
}
else {
z = 8. / x;
y = z * z;
xx = x - .785398164;
return(sqrt(.636619772 / x) * (sin(xx) *
(p1 + y * (p2 + y * (p3 + y * (p4 + y * p5)))) + z * cos(xx) *
(q1 + y * (q2 + y * (q3 + y * (q4 + y * q5))))));
}
}
/********************************************************************
FUNCTION NAME:
ARGUMENTS:
PURPOSE:
********************************************************************/
double bessy1q(double x)
{
double y, z, xx;
double p1 = 1.;
double p2 = .183105e-2;
double p3 = -.3516396496e-4;
double p4 = .2457520174e-5;
double p5 = -.241337019e-6;
double q1 = .04687499995;
double q2 = -.2002690873e-3;
double q3 = .8449199096e-5;
double q4 = -.88228987e-6;
double q5 = .105787312e-6;
double r1 = -.4900604943e13;
double r2 = .1275274390e13;
double r3 = -.5153438139e11;
double r4 = .7349264551e9;
double r5 = -.4237922726e7;
double r6 = .8511937935e4;
double s1 = .2499580570e14;
double s2 = .4244419664e12;
double s3 = .3733650367e10;
double s4 = .2245904002e8;
double s5 = .1020426050e6;
double s6 = .3549632885e3;
double s7 = 1.;
if(x < 8.){
y = x * x;
return(x * (r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y *
r6))))) / (s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y *
(s6 + y * s7)))))) + .636619772 * (bessj1q(x) * log(x) - 1. / x));
}
else {
z = 8. / x;
y = z * z;
xx = x - 2.356194491;
return(sqrt(.636619772 / x) * (sin(xx) * (p1 + y * (p2 + y * (p3 +
y * (p4 + y * p5)))) + z * cos(xx) * (q1 + y * (q2 + y * (q3 +
y * (q4 + y * q5))))));
}
}
#endif