BAR3D.C

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

        This program plots 3-D bars.
****************************************************************************/

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

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

#define DIM 31                                   /* Define array dimension */

float a = .30f;                                         /* Plot parameters */
float R0 = 1.0f;
float Te0 = 850.f;                                                /* In eV */
float Eb = 25000.f;                                               /* In eV */
float ne0 = 10.e19f;                                             /* In m-3 */
float Ec;
float taus;
float vDrift;
double lambda0;
float s[DIM],                              /* Distance along the beam path */
    Nb[DIM],                    /* Number of beam particles along the path */
    Ni[DIM];                    /* Number of fast ions born along the path */

double xsect(double EeV);
float W_CDECL shine(float E, float n0);
float W_CDECL r0q(float E, float n0);

/****************************************************************************
    Main program
****************************************************************************/
void GPC_MAIN(void)
{
    bgnplot(1, 'g', "bar3d.tkf");                /* GraphiC initialization */
    startplot(WHITE);
    metricunits(0);                                 /* Plots are in inches */
    font(1, "swiss.fnt", '\310');                          /* Select fonts */
    sclfix(1);
    fillfont(1);

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

    color(BLACK);                        /* Draw 3-D axes with axis titles */
    volm3d(2.0f, 3.0f, 1.0f);
    vuabs(200.f, -100.f, 100.f);
    x3name("\310beam energy (keV)");
    y3name("\310density x 10[19] m[-3]");
    z3name("\310shine-through");
    graf3d(5.f, 5.f, 40.f, 1.f, 1.f, 10.f, 0.0f, .2f, 1.0f);

    color(BLUE);                                /* Plot base grid (z == 0) */
    surfun(r0q, 0, 7.f, 0, 1.f);

    axesoff(1);                           /* Redraw axes to reset horizons */
    graf3d(5.f, 5.f, 40.f, 1.f, 1.f, 10.f, 0.0f, .2f, 1.0f);

    color(RED);                                           /* Plot 3-D bars */
    d3bars(.2f);
    surfun(shine, 0, 7.f, 0, 1.f);
    d3bars(0.f);

    color(BLACK);                                            /* Plot title */
    d3head("\310Beam Shine Through in a Tokamak", 't');

    endplot();                                    /* Plot closing routines */
    stopplot();
}

/****************************************************************************
    Uses Riviere's Beam deposition cross sections  in m2
****************************************************************************/
double xsect(double EeV)
{
    double scx, si, lte;

    lte = log10(EeV);
    scx = (0.6937e-14 * pow((1. - 0.155 * lte), 2.)) /
        (1. + 0.1112e-14 * pow(EeV, 3.3));
    si = pow(10.0, -.8712 * lte * lte + 8.156 * lte - 34.833);
    return(1.e-4 * (si + scx));
}

/****************************************************************************
    Return the beam shine-through
****************************************************************************/
float W_CDECL shine(float E, float n0)
{
    double lam, n, Ed;
    float rv;

    n = 1.e19 * n0;
    Ed = 1000. * E;
    lam = 1. / (n * xsect(Ed));
    rv = (float)exp(-(double)(1.33 * a) / lam);
    return(rv);
}

/****************************************************************************
    Dummy routine (always returns zero)
****************************************************************************/
float W_CDECL r0q(float E, float n0)
{
    E = 0.f;
    n0 = 0.f;
    return((float)(E + n0));
}