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