ORBITS.C

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

    This program plots Tokamak Banana Orbits.
****************************************************************************/

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

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

double  q_0 = 3.0,
        R_0 = 1.0,
        e = 1.602e-19,
        m = 1.67e-27,
        RM_T = 0.75,
        r_T = 0.3,
        PtoP = 0.8,
        v = 2.0e6;
float C;
char print=0;

float W_CDECL zfun1(float x, float y);
float W_CDECL zfun2(float x, float y);
void orbit(void);

/****************************************************************************
    Main program
****************************************************************************/
void GPC_MAIN(void)
{
    bgnplot(1, 'g', "orbits.tkf");        /* Initialize and name disk file */
    orbit();
    stopplot();                                    /* Close files and exit */
}

/****************************************************************************
    Plot banana orbits
****************************************************************************/
void orbit(void)
{
    float x3axis, y3axis, z3axis, xvu, yvu, zvu, xorig, xstp, xmax;
    float yorig, ystp, ymax, zorig, zstep, zmax;
    float xdel, ydel, x[21], yp[21], ym[21], z[21], t1;
    int i, np1, np2, ixpts, iypts;
    double temp;

    C = (float)(2. * q_0 * m * v * sqrt(PtoP) / e);
    ixpts = 2;                                 /* Two segments across cell */
    iypts = 2;
    x3axis = 3.0f;                 /* x, y, and z axes are the same length */
    y3axis = 3.0f;
    z3axis = 3.0f;
    xvu =-3.f;                /* x distance to view point in workbox units */
    yvu = 24.f;               /* y distance to view point in workbox units */
    zvu = 7.f;                /* z distance to view point in workbox units */
    xorig = (float)RM_T;                        /* Sets left end of x axis */
    xstp = .25f;
    xdel = .05f;                                    /* Sets size of x step */
    xmax = 1.5f;                              /* x ranges from 0.75 to 1.5 */
    yorig = 0.0f;                                     /* Set y-axis limits */
    ystp =  0.1f;
    ydel = .05f;
    ymax = .5f;
    zorig = -.5f;                                     /* Set z-axis limits */
    zstep = .25f;
    zmax = .25f;
    startplot(WHITE);                                  /* White background */
    font(3, "swissbld.fnt", '\310', "news.fnt", '\311',    /* Select fonts */
        "newsgrm.fnt",'\312');
    fillfont(1);                                    /* Enable font filling */
    metricunits(0);                                 /* Plots are in inches */
    sclfix(1);                               /* All plots use same scaling */

    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);                              /* Inputs viewpoint */

    color(RED);                                      /* Axis labels in red */
    x3name("\310R (meters)");
    y3name("\310z (meters)");
    z3name("\310P]\312f[ \310-P]\312f]\310TIP[[");

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

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

    axesoff(1);                     /* Second call resets the 3-D horizons */
    graf3d(xorig, xstp, xmax, yorig, ystp, ymax, zorig, zstep, zmax);

    tcolor(RED);                /* Surface with red top and magenta bottom */
    bcolor(MAGENTA);
    surfun(zfun1, ixpts, xdel, iypts, ydel);

    for(i = 0; i < 21; i++) {
        x[i] = xorig + (xmax - xorig) * .05f * i;       /* x scans R major */
        z[i] = 0.0f;
    }
    for(i = 0; i < 21; i++) {
        temp = R_0 / RM_T - R_0 / (double)x[i];
        if (temp > 0)
                temp = sqrt(temp);
        t1 = (C) * x[i] * (float)temp;
        yp[i] = (float)(r_T * r_T + t1 - (x[i] - R_0) * (x[i] - R_0));
        ym[i] = t1;
        if(yp[i] <= 0.0f)
                break;
        yp[i] = (float)sqrt((double)yp[i]);
    }
    np1 = i;
    tcurve(12);                           /* Curve will be 12 pixels thick */
    color(LGT_RED);                                     /* Light red curve */
    curv3d(x, yp, z, i);


    for(i = 0; i < 21; i++) {
        ym[i] = (float)(r_T * r_T - ym[i] - (x[i] - R_0) * (x[i] - R_0) );
        if(ym[i] <= 0.0f)
                break;
        ym[i] = (float)sqrt((double)ym[i]);
    }
    np2 = i;
    color(LGT_GREEN);                                 /* Light green curve */
    curv3d(x, ym, z, i);

    plane3d(1, zorig, 0);         /* Establish a plane for 2-D projections */
    tcurve(16);                           /* Curve will be 16 pixels thick */
    color(LGT_RED);                                     /* Light red curve */
    curve(x, yp, np1, 0);

    color(LGT_GREEN);                                 /* Light green curve */
    curve(x, ym, np2, 0);
    tcurve(0);                                     /* Reset line thickness */

/*
    It is best to put the curve labels in the X-Z plane to 1) prevent
    distortion of the characters, and 2) so that the labels move
    automatically as the view point is changed. z=.1 is selected for the
    plane location because this is close to the ends of the curves.
*/

    plane3d(2, .1f, 0);                      /* x-z plane for curve labels */
    color(RED);                                    /* Red vector and label */
    uvector(1.25f, -.45f, x[np1-1]-.01f, -0.5f, 2.f, .1f, "01");
    pltfnt(1.25f, -.45f, "\311v\312]7[\311>0", .2f, 0);

    color(GREEN);                                /* Green vector and label */
    uvector(.75f, -.40f, x[np2-1]-.05f, -0.5f, 2.f, .1f, "01");
    pltfnt(.75f, -.40f, "\311v\312]7[\311<0", .2f, 0);
    plane3d(0, 0.f, 0);                         /* Turn off 2-D projection */

    color(BLUE);                                        /* Blue plot title */
    titlht(.2f);
    d3head("\310Topology of Tokamak Banana Orbits", 't');

    endplot();                                           /* Terminate plot */
}

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

    RM = xx - R_0;
    rsq = RM * RM + yy * yy;
    if ((double)xx > RM_T)
        t1 = (C) * xx * sqrt((R_0 / RM_T) - (R_0 / xx));
    else
        t1 = (double)0.0;
    t2 = (rsq - r_T * r_T);
    return((float)(t1 - t2));
}

/****************************************************************************
    Function to be plotted
****************************************************************************/
float W_CDECL zfun2(float xx, float yy)
{
    double rsq, RM, t1, t2;

    RM = xx - R_0;
    rsq = RM * RM + yy * yy;
    if ((double)xx >= RM_T)
        t1 = (C) * xx * sqrt((R_0 / RM_T) - (R_0 / xx));
    else
        t1 = (double)0.0;
    t2 = (rsq - r_T * r_T);
    return((float)(-t1 - t2));
}