FERMI.C

/* This program gives an example of doing simultaneous plots with GraphiC.
Custom y-axis labels are used.

The example is a simple model of Fermi acceleration. Fermi claimed that
cosmic rays could be accelerated by a similar method to very high velocities.
However, as you can see from the results, the method does not work well
because the stochastic acceleration regions are divided by nonstochastic
regions at M, 2M, etc. The brown points start on this separatrix. They can
sneak through the region at M, but it is unlikely that the trajectory will
also sneak through the separatrix at 2M.

A good reference for this topic is Joseph Ford's "A Picture Book of
Stochasticity" on page 140 in the book Topics in Nonlinear Dynamics: A
Tribute to Sir Edward Bullard, edited by Siebe Jorna, American Institute of
Physics, 1978.

A light ball bounces in one dimension between two massive walls. One wall
oscillates with a velocity given by
        V(t) = (V/4)[1 - 2{t}]
where {t} is the fractional part of t (denoted by ft in the code).

      V(t)
        |
        |
     V/4|\    |\    |\    |\    |\    |\
        | \   | \   | \   | \   | \   | \   |
      --|--\--|--\--|--\--|--\--|--\--|--\--|- t
        |0  \ |1  \ |2  \ |3  \ |4  \ |5  \ |6
        |    \|    \|    \|    \|    \|    \|
        |

The velocity of the ball just before the nth collision with the oscillating
wall is

    Un+1 = |u + {tn} - 1/2|

and if the space between the walls (d) is much greater than the maximun
amplitude of oscillation (V/16 = 2a),

    {tn+1} = {{tn} + M/Un+1}

where u = the velocity of the particle/v, and M = d/16a.
This program simultaneously plots the velocity of the particle and its
phase space trajectory for several starting points.
*/

#include <graphic.h>

#define NCURVE 7

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

float M = 10.0f;
int nmax = 4000;
double ustart[NCURVE];
double tstart[NCURVE];
int stepmax[NCURVE];

/*******************************************************************/
void gpc_init(void);
void GPC_MAIN(void)
{
    int i, ncolor = 14;
    int ch;
    int n = 0;                                         /* Iteration number */
    float u[2];                      /* Relative velocity before collision */
    float ft[2];                                  /* Fractional parts of t */
    float nstep[2];                     /* Steps (need float for the plot) */

    bgnplot(1, 'g', "notek");      /* Suppress the very large TKF file */
    plotabsolute(1);                 /* Prevent GraphiC from resizing axes */
    metricunits(0);                /* Ensure scaling is done in inch units */
    startplot(BLACK);
    font(3,"simplex.fnt",'\310',"swiss.fnt",'\311',"simgrma.fnt",'\312');
    fillfont(1);
    gpc_init();
    ustart[0] = ustart[1] = M;         /* Pick interesting starting points */
    stepmax[0] = stepmax[1] = 500;
    tstart[0] = .1f;
    tstart[1] = .4f;
    ustart[2] = .4f*sqrt((double)M);
    tstart[2] = .5f;
    stepmax[2] = 4000;
    ustart[3] = 1.45f*M;
    tstart[3] = .5f;
    stepmax[3] = 2000;
    ustart[4] = 1.5f*M;
    tstart[4] = 0.0f;
    stepmax[4] = 4000;
    ustart[5] = 2.0f*M;
    tstart[5] = 0.5f;
    stepmax[5] = 4000;
    ustart[6] = M;
    tstart[6] = 0.0f;
    stepmax[6] = 4000;
    for (i = 0; i < NCURVE; i++) {
        u[0] = (float)ustart[i];
        ft[0] = (float)tstart[i];
        nstep[0] = 0.0f;
        context(0);
        symbol(ft[0], u[0], 16);
        for (n = 1; n < stepmax[i]; n++) {
            nstep[1] = (float)n;
            u[1] = u[0] + ft[0] - .5f;
            if (u[1] < 0.0f)
                u[1] = -u[1];
            ft[1] = ft[0] + (M/u[1]);
            while (ft[1] > 1.0f)
                ft[1] -= 1.0f;
            context(0);                                   /* First context */
            symbol(ft[1], u[1], 16);
            context(1);                                  /* Second context */
            curve(nstep, u, 2, 0);
            u[0] = u[1];
            ft[0] = ft[1];
            nstep[0] = nstep[1];
            if ((ch = cstsq()) == 'q')
                goto end;
        }
        ncolor++;
        if (ncolor==15)
            ncolor = 1;
        if (ncolor==7)
            ncolor = 9;
        context(0);            /* Change colors for the different contexts */
        color(ncolor);
        context(1);
        color(ncolor);
    }
end:
    simuplot(0);   /* Must end the simultaneous plots to make the TKF file */
    endplot();
    stopplot();
}

/***************************************************************************/
/* Set up the two contexts for the simultaneous plots */
void gpc_init(void)
{
    simuplot(2);                                           /* Two contexts */
    context(0);/* Establish page and plot parameters for the first context */
    color(BRT_WHITE);
    frame(1, 1);
    tcurve(2);
    page(4.0f, 6.884f);
    pgshift(0.5f, 0.0f);
    symht(.05f);                         /* This call must come after page */
    area2d(3.5f, 5.7f);               /* Set the area of the Poincar‚ plot */
    xname("\311{t} = Phase\310");
    yname("\311u = Relative Velocity\310");
    graf("%-1.2f", 0.0f, 1.0f, 1.0f, " ", 0.0f, 0.5f*M, 2.0f*M);
    pltfnt(-.02f, 0.0f, "\310""0", .1f, 90);   /* Custom labels for y axis */
    pltfnt(-.02f, M, "\310M", .1f, 90);
    pltfnt(-.02f, 2.f*M, "\310""2M", .1f, 90);
    pltfnt(-.02f, (float)(sqrt((double)M))*.5f, "\310(\312""1\310M)/2", .1f, 90);

/* LABEL THE PLOT */
    color(BLUE);                /* Outline the title in blue (for printer) */
    prtfnt(5.0f, 1.5f, "\311|14|Fermi Acceleration", 0.2f, 0);
    color(BRT_WHITE);
    tfont(3);
/* The `#' toggles the thick fonts */
    prtfnt(5.0f, 1.0f, "\310#u]n+1[ = @|u]n[ + {t]n[} - \312A\310 @|", 0.18f, 0);
    prtfnt(5.0f, 0.5f, "\310{t]n+1[} = \312{\310{t]n[} + M/u]n+1[\312}#\310", 0.18f, 0);
    tfont(0);
    tcurve(0);
    color(YELLOW);                                    /* Set context color */

/* START VELOCITY PLOT */
    context(1);/* Establish page and plot parameters for the second context */
    color(BRT_WHITE);
    symht(.01f);                         /* This call must come after page */
    page(4.5f, 4.0f);
    pgshift(4.5f,2.3f);
    axnamht(.3f);
    xname("\311Iteration");
    yname("\311u = Relative Velocity");
    area2d(3.5f, 4.0f);
                   /* The `-' is to left align the axis labels */
    graf("%-1.0f", 0.0f, 0.1f*nmax, (float)nmax,
        " ", 0.0f, .5f*M, 2.f*M);
                               /* Custom labels for y-axis */
    pltfnt(-.02f*nmax, -0.0f, "\310""0", .08f, 90);
    pltfnt(-.02f*nmax, M, "\310M", .08f, 90);
    pltfnt(-.02f*nmax, 2.f*M, "\310""2M", .08f, 90);
    pltfnt(-.02f*nmax, (float)(sqrt((double)M))*.5f, "\310(\312""1\310M)/2", .08f, 90);
    color(YELLOW);
}