D3TEST.C

/****************************************************************************
    (c) 1984 - 2003 by Scientific Endeavors Corporation.
    All rights reserved.

    This program does two plots. The first plot is a 3-D surface generated
    from a function using surfun().
****************************************************************************/

#include <graphic.h>                           /* Include all needed files */

#if TCQ                                 /* Set stack for Borland (Turbo) C */
extern unsigned _stklen = 0x3000;
#endif

double *carg, *sarg;                                   /* Define variables */
float *c, *s, *u, *one, *zero;
unsigned int memsize = 601;

int trapezoid(int npts, double h, double *integrand, float *integral);
float W_CDECL zfun(float x, float y);
int surface(void);
void fresnel(void);

/****************************************************************************
    Main program
****************************************************************************/
void GPC_MAIN(void)
{
                                            /* Allocate arrays */
    carg = (double *)gpcalloc(memsize, sizeof(double));
    sarg = (double *)gpcalloc(memsize, sizeof(double));
    c = (float *)gpcalloc(memsize, sizeof(double));
    s = (float *)gpcalloc(memsize, sizeof(double));
    u = (float *)gpcalloc(memsize, sizeof(double));
    one = (float *)gpcalloc(memsize, sizeof(double));
    zero = (float *)gpcalloc(memsize, sizeof(double));
    if (carg == (void *)NULL || sarg == (void *)NULL ||
        c == (void *)NULL || s == (void *)NULL ||
        u == (void *)NULL || one == (void *)NULL ||
        zero == (void *)NULL) {
        GPC_PUTS("D3TEST: Can't allocate data arrays");
        goto EndOfApp;
    }

    bgnplot(1, 'g', "d3.tkf");                   /* GraphiC initialization */
    metricunits(0);                                 /* Plots are in inches */

    if (surface())                                     /* Plot 3-D surface */
        fresnel();                               /* Plot Fresnel integrals */

    stopplot();                                    /* Close files and exit */

EndOfApp:
    gpcfree((void **)&carg);                                /* Free arrays */
    gpcfree((void **)&sarg);
    gpcfree((void **)&c);
    gpcfree((void **)&s);
    gpcfree((void **)&u);
    gpcfree((void **)&one);
    gpcfree((void **)&zero);
}

/****************************************************************************
    Plots a 3-D surface
****************************************************************************/
int surface(void)
{
    float x3axis, y3axis, z3axis, xvu, yvu, zvu;       /* Define variables */
    float xorig, xstp, xmax, yorig, ystp, ymax, zorig, zstep, zmax;
    float xdel, ydel;
    float pi = (float)PI;
    int ixpts, iypts;

    ixpts = 2;                                 /* Two segments across cell */
    iypts = 2;
    x3axis = 2.f;                  /* X axis is twice as long as other two */
    y3axis = 1.f;
    z3axis = 1.f;
    xvu = 8.8f;                  /* Distance to viewpoint in workbox units */
    yvu = 3.66f;
    zvu = 3.f;
    xorig = -pi;                                         /* X-axis scaling */
    xstp = 1.f;
    xdel = pi / 15.f;
    xmax = pi;
    yorig = 0.f;                                         /* Y-axis scaling */
    ystp =  1.0f;
    ydel = pi / 15.f;
    ymax = pi;
    zorig = -1.f;                                        /* Z-axis scaling */
    zstep = .5f;
    zmax = 1.f;

    startplot(WHITE);                           /* Initialize surface plot */
    font(2, "simplex.fnt", '\310', "triplex.fnt", '\311'); /* Select fonts */

    page(9.0f, 6.884f);                      /* Size of page and plot area */
    area2d(7.6f, 5.5f);

    color(BLACK);                                 /* Black box around plot */
    box();

    color(BROWN);                                      /* Brown axis names */
    x3name("X-axis");
    y3name("Y-axis");
    z3name("Z-axis");

    volm3d(x3axis, y3axis, z3axis);                  /* Define 3-D workbox */
    vuabs(xvu, yvu, zvu);                                 /* Set viewpoint */

    color(BLUE);                                   /* Blue axes and labels */
    tickout(1);                                      /* Ticks face outward */
    graf3d(xorig, xstp, xmax, yorig, ystp, ymax, zorig, zstep, zmax);

    tcolor(BROWN);        /* Top is brown, bottom is magenta, base is blue */
    bcolor(MAGENTA);
    base(1, 1);
    base(1, 0);
    survis("both");                /* Draw surface with both sides visible */
    surfun(zfun, ixpts, xdel, iypts, ydel);

    color(MAGENTA);                                  /* Magenta plot title */
    d3head("\311sin(y) * cos(x-y)", 't');

    return(!endplot());                          /* Terminate surface plot */
}

/****************************************************************************
    Define function to be plotted
****************************************************************************/
float W_CDECL zfun(float xx, float yy)
{
    return(float)(sin((double)yy)*cos((double)(xx - yy)));
}

/****************************************************************************
    Plot Fresnel integrals
****************************************************************************/
void fresnel()
{
    double t;                                          /* Define variables */
    double piby2 = (float)PI_OVER_TWO, tsq;
    unsigned int i;

    for(i = 0; i < memsize; i++) {                        /* Generate data */
        t = 1.e-2 * (double)i;
        tsq = t * t;
        u[i] = (float)t;
        one[i] = 1.0f;
        zero[i] = 0.0f;
        carg[i] = cos(piby2 * tsq);
        sarg[i] = sin(piby2 * tsq);
    }
    trapezoid(memsize, .01, carg, c);
    trapezoid(memsize, .01, sarg, s);

    startplot(WHITE);                           /* Initialize Fresnel plot */
    font(4, "microb.fnt", '\310', "duplex.fnt", '\311',    /* Select fonts */
        "complex.fnt", '\312', "compgrma.fnt", '\313');

    page(9.0f, 6.884f);                      /* Size of page and plot area */
    physor(1.5f, -.5f);
    area2d(7.6f, 5.5f);

    color(BLACK);                                 /* Black box around plot */
    box();

    volm3d(7.0f, 4.0f, 4.0f);                  /* Define 3-D workbox units */
    vuabs(30.f, -20.f, 30.f);                             /* Set viewpoint */

    color(GREEN);                                      /* Green axis names */
    fntchg((Uchar)'\311');
    titlht(.25f);
    numht(.15f);
    axnamht(.20f);
    tfont(5);
    x3name("#z#");
    y3name("#C(z)#");
    z3name("#S(z)#");

    graf3d(0.0f, 1.0f, 6.0f, 0.0f, 0.5f, 1.0f, 0.0f, 0.5f, 1.0f);

    dashf(1);                                           /* Draw grid lines */
    d3line(0.0f, 0.0f, 0.0f, 0.0f, 1.0f, 0.0f);
    d3line(0.0f, 1.0f, 0.0f, 0.0f, 1.0f, 1.0f);
    d3line(0.0f, 1.0f, 1.0f, 0.0f, 0.0f, 1.0f);
    d3line(0.0f, 1.0f, 0.0f, 6.0f, 1.0f, 0.0f);
    d3line(6.0f, 1.0f, 0.0f, 6.0f, 1.0f, 1.0f);
    d3line(6.0f, 1.0f, 1.0f, 0.0f, 1.0f, 1.0f);
    dashf(2);
    for(i = 1; i < 6; i++){
        d3line((float)i, 1.0f, 0.0f, (float)i, 1.0f, 1.0f);
        d3line((float)i, 0.0f, 0.0f, (float)i, 1.0f, 0.0f);
        if (i < 5)
            d3line(0.0f, (float)i * .2f, 0.0f, 0.0f, (float)i * .2f, 1.0f);
    }
    dashf(1);

    color(GREEN);                                      /* Green plot title */
    tfont(3);
    skew(.3f);
    prtfnt(0.2f, 6.1f, "`\310#Fresnel Integrals#`", .3f, 0);
    tfont(0);

    color(MAGENTA);                            /* Magenta Fresnel integral */
    tcurve(5);
    curv3d(u, c, s, memsize);
    tcurve(3);

    color(BLUE);                                /* Blue x-z component C(z) */
    charspc(2);
    prtfnt(0.2f, 5.6f,
        "\312#C(z)=\313,\312[[ z]]]]0[[cos(\313Ap\312t[2])dt", .15f, 0);
    charspc(0);
    curv3d(u, one, s, memsize);

    color(RED);                                  /* Red x-y component S(z) */
    charspc(2);
    prtfnt(0.2f, 5.1f,
        "\312S(z)=\313,\312[[ z]]]]0[[sin(\313Ap\312t[2])dt\310", .15f,0);
    charspc(0);
    curv3d(u, c, zero, memsize);

    color(LGT_CYAN);                         /* Bright cyan Cornu's spiral */
    tfont(2);
    prtfnt(0.5f, .3f, "\311`#Cornu\'s Spiral#`", .2f, 0);
    tfont(0);
    tcurve(5);
    curv3d(zero, c, s, memsize);

    for(i = 0; i < memsize; i++) {
        c[i] = -c[i];
        s[i] = -s[i];
    }
    curv3d(zero, c, s, memsize);

    for(i = 0; i < memsize; i++) {
        c[i] = -c[i];
        s[i] = -s[i];
    }

    dateit((Uchar)'\311');

    endplot();                                   /* Terminate Fresnel plot */
}

/****************************************************************************
    Integrate equally spaced sample points with the
    extended trapezoidal rule:
    h[f(0)/2+f(1)+.....+f(m-1)+f(m)/2]
****************************************************************************/
int trapezoid(int npts, double h, double *integrand, float *integral)
{
    double hby2, sum, newterm;
    int i;

    if(npts < 2)
        return(0);
    hby2 = h * .5;
    integral[0] = 0;
    sum = hby2 * integrand[0];
    for(i = 1; i < npts; i++) {
        newterm = hby2 * integrand[i];
        sum += newterm;
        integral[i] = (float)sum;
        sum += newterm;
    }
    return(1);
}