JCONT.C

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

    This program plots a family of proton tokamak orbits. If `_'= sub, and
    `^' = sup, the orbit equation for a family of orbits can be determined
    from the conservation of the second adiabatic invariant, J. For a simple
    model field, this is given by:
        J = 1-(r^2-1)^2 + eps*r*cos(th)
    eps is the inverse aspect ratio, a/R_0. All coordinates are normalized.
****************************************************************************/

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

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

float eps = .4f;
int lnstyle[ ] = {
    211, 221, 231, 241, 251, 211, 221, 231, 241, 251, 211
};
int lnstyle2[ ] = {
    211, 221, 231, 241
};
float x[121], y[121], level[1] = {1.15f};


float W_CDECL zfun1(float x, float y);
void jorbit(void);

/****************************************************************************
    Main program
****************************************************************************/
void GPC_MAIN(void)
{
    bgnplot(1, 'g', "jcont.tkf");       /* Initializes and names disk file */
    jorbit();
    stopplot();                                    /* Close files and exit */
}

/****************************************************************************
    Plot tokamak orbits
****************************************************************************/
void jorbit(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 = 2;                                 /* Two segments across cell */
    iypts = 2;
    x3axis = 5.5f;                       /* x and y axes are the same size */
    y3axis = 5.5f;
    z3axis = 5.0f;
    xvu = -3.f;                  /* Distance to viewpoint in workbox units */
    yvu = 24.f;
    zvu = 15.f;
    xorig = -1.5f;                               /* Set left end of x axis */
    xstp = .5f;
    xdel = .1f;                                      /* Set size of x step */
    xmax = 1.5f;                                 /* x ranges from -2 to +2 */
    yorig = -1.5f;                                    /* Set y-axis limits */
    ystp =  0.5f;
    ydel = .1f;
    ymax = 1.5f;
    zorig = -1.0f;                                    /* Set z-axis limits */
    zstep = .5f;
    zmax = 3.0f;

    startplot(BLACK);                                  /* Black background */
    font(3, "swissbld.fnt", '\310', "news.fnt", '\311',    /* Select fonts */
        "newsgrm.fnt", '\312');
    fillfont(1);
    metricunits(0);                                 /* Plots are in inches */

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

    volm3d(x3axis, y3axis, z3axis);                  /* Input workbox size */
    vuabs(xvu, yvu, zvu);                               /* Input viewpoint */

    color(RED);                                       /* Axis names in red */
    x3name("\310R-R]0[");
    y3name("\310y");
    z3name("\310J");

    color(WHITE);                                        /* White 3-D axes */
    tickout(1);                                      /* Ticks face outward */
    numht(.15f);
    graf3d(xorig, xstp, xmax, yorig, ystp, ymax, zorig, zstep, zmax);

    survis("both");                   /* Both top and bottom surfaces seen */
    tcolor(GREEN);               /* Surface with top green and bottom cyan */
    bcolor(CYAN);
    surfun(zfun1, ixpts, xdel, iypts, ydel);

    plane3d(1, zmax, 0);                      /* Plot contour in x-y plane */
    for(i = 0; i < 121; i++) {
        x[i] = y[i] = -1.5f + .025f * i;
    }
    contour(zfun1, x, 121, y, 121, 0.0f, 0.2f, 2.0f, lnstyle, 0);

    stack(1);                                   /* Plot contour on surface */
    tcurve(10);
    contour(zfun1, x, 121, y, 121, 1.0f, .2f, 1.6f, lnstyle2, 0);
    tcurve(0);
    plane3d(0, 0.f, 0);                         /* Turn off 2-D projection */

    fontfill(2, -RED, BRT_WHITE);        /* Plot title and function in red */
    titlht(.2f);
    d3head("\310Orbits correspond to J=constant contours", 't');
    prtfnt(1.0f, 0.5f, "\311J=1-(r[2]-1)[2]+\312e\311cos(\312q\311)", .2f,0);

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

/****************************************************************************
    Function to be plotted
****************************************************************************/
float W_CDECL zfun1(float xx, float yy)
{
    double rsq;
    float fn;

    rsq = xx * xx + yy * yy;
    fn = (float)(1. - (rsq - 1.) * (rsq - 1.) - eps * xx);
    if (fn < -1.0f)
        fn = -1.0f;
    return(fn);
}