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));
}