TZTEST.C

/*********************************************************************
This program will generate plots of air traffic over the United
States superimposed on a map which is expressed in terms of a series
of line segments.  This version will only plot that map in the input
coordinates (assumed to be lat-long) with no projection.  The basic
map data is the file "state.chn"
**********************************************************************/

#include <GraphiC.h>

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

#define LATN 50.
#define LATS 24.
#define LATP 53
#define LONP 121
#define LONGW 126.
#define LONGE 66.
#define DELTA_DEGREE 0.5

FILE *in_ptr;                                 /* pointer to the input file */
char density_file_name[] = "tztest.dat";
char output_file_name[] = "tztest.tkf";
int clr[] = {
    1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15
};

/* Color plot suggested values */
/* For contour plots */
int intensityc[] = {
    174, 173, 171, 170, 145, 147, 120, 117, 122,
    95, 96, 97, 98, 99, 124, 118, 119, 74,
    73, 72, 71, 70, 91, 90, 93, 92, 68, 87,
    94, 69, 58, 64, -1, 63, 59, 84, 83, 113, 103,
    109, 77, 104, 108, 133, 134, 101, 139, 138, 169,
    168, 164, 159, 158, 154, 153, 152, 151, 150
};
int nintensityc = 58;

/* For 3-d plots */
int intensity[] = {
    174, 173, 170, 145, 120, -2, 123, 98, 74, 69,
    -1, 144, 77, 108, -5, 159, 164, 150, 101, 107
};
int nintensity = 20;

/* The following Gray-scale settings work well for black-and-white plots */
/*
int intensity[] = {
    231, 230, 229, 228, 227, 226, 225, 224,
    223, 222, 221, 220, 219, 218, 217, 216
};
int nintensity = 16;
*/

/* Patch plot arguments */
int ppsmooth[2] = {0, 2};
int pptype[2] = {'c', 'a'};

char hctype = 0;
float density_max;
int error_number = 0;

/**************** Function Declarations  ******************/
void map_x_y(float *parea_x, float *parea_y);
void block_r(float**px, float**py, int *pi_flag, int*pn_pt);
void plot_block(float*x, float*y, int n_pt);
void plot_density(int nplot);

/*********************************************************************/
/*
First, we make 2-D plots of the data. In this case, the data
represents aircraft locations reported at 5-minute intervals.
Therefore, it does not (and should not) represent a smooth function.
Data may vary dramatically from one grid point to the next. For data
of this sort, a "patch plot" is the best way of representing this
discontinuous nature. In the second plot, we want to make a contour
plot of the data, but it must be smoothed first in order that the
contour-plotting algorithm can work properly. The call to smear()
eliminates the jagged peaks and cliffs but also changes the
appearance of the plot. You, the scientist user, must evaluate the
tradeoffs between these two approaches.

Finally, we make a 3-d surface plot with contour levels. This sort of
plot must be even smmother than the 2-D contour plot, so the data is
smeared some more.
*/
void GPC_MAIN(void)
{
    int i_flag;/* line segments will be read until this flag is set to EOF */
    int n_pt;                          /* number of segments in a map line */
    float area_x, area_y;              /* plot area to preserve dimensions */
    char *str = "\302Aircraft density over the U.S. at 20:00-21:00 GMT";
    float *x, *y;                       /* Arrays for generating map lines */
    int nplot;                                              /* Plot number */

    bgnplot(1,  'g',  output_file_name);
    font(3, "simplex.fnt", '\301', "swiss.fnt", '\302',    /* Select fonts */
        "newsgrm.fnt", '\303');
    fillfont(1);
    physor(0.0f, 1.0f);
    for (nplot = 0; nplot <= 1; nplot++) {
        startplot(WHITE);                     /* Get GraphiC to setup plot */
        color(0);
        rotate(0);
        page(9.00f,  6.844f);
        tmargin(0.0f);
        ctline(str, .2f);
        linesp(2.0f);

        ctline("\302""11/5/91", .15f);
        ctline("\302(planes/hour/10,000 mi[2])", .15f);
        map_x_y(&area_x, &area_y);
        plotabsolute(1);
        area2d(area_x, area_y);
        area_x = (float)((9.0 - area_x) / 2.0);
        area_y = (float)((6.884 - area_y) / 2.0);
        physor(area_x, area_y);
        axesoff(AXESOFF);
        graf("%-4.1f", -(float)LONGW, 5.f, -(float)LONGE, "%-4.1f",
            (float)LATS, 4.5f, (float)LATN);

        plot_density(nplot);
        if (error_number)
            goto Error;

        /* NOW READ IN MAP DATA AND CYCLE THROUGH THE PLOT ROUTINE */
        color(BLACK);
        i_flag = 0;
        n_pt = 0;
        tcurve(3);                           /* Thick state outlines lines */
        while (i_flag != EOF) {                         /* Plot each state */
            block_r(&x, &y, &i_flag, &n_pt);
            if (error_number)
                goto Error;
            if (i_flag != EOF) {
                plot_block(x, y, n_pt);
                if (error_number)
                        goto Error;
            }
            else {
                block_r(&x, &y, &i_flag, &n_pt); /* Free memory */
                if (error_number)
                        goto Error;
                break;
            }
        }
        tcurve(1);                             /* Reset the line thickness */
        if (hctype) {              /* Allow unattended plotting if desired */
            hardcopy(hctype);
        }
        if (endplot())
            goto EndOfApp;
    }
    /* Make a 3-d surface plot with contour levels */
    /* The next call to physor resets it to its default */
    physor(-1000.0f, -1000.0f);
    axesoff(AXESON);                                  /* Turn axes back on */
    startplot(WHITE);                         /* Get GraphiC to setup plot */
    color(BLUE);
    rotate(0);
    font(3, "microb.fnt", '\301', "swiss.fnt", '\302',     /* Select fonts */
        "newsgrm.fnt", '\303');
    page(9.00f,  6.844f);
    area2d(8.5f, 6.5f);
    tmargin(0.5f);
    linesp(1.5f);
    ctline(str, .2f);
    ctline("\302""11/5/91", .2f);
    color(BLACK);
    ctline("\301(planes/hour/10,000 mi[2])", .15f);
    volm3d(80.5f, 50.f, 50.f);
    vuabs(110.f, -500.f, 1200.f);
    color(BLACK);
    x3name("\302Longitude");
    y3name("\302North Latitude");
    z3name("\302Aircraft Density");
    fntchg((Uchar)'\301');
    graf3d( -(float)LONGW, 10.f, -(float)LONGE,
            (float)LATS, 9.0f, (float)LATN, 0.f, 20.f, 40.f );
    d3color(intensity, nintensity, FILL_ON_CONTOURS);
    plot_density(2);
    if (error_number)
        goto Error;

    color(BLACK);
    plane3d(1, 0.f, 0);

    /* NOW READ IN MAP DATA AND CYCLE THROUGH THE PLOT ROUTINE */
    color(WHITE); /* Erase state outlines first to get good printer output */
    i_flag = 0;
    tcurve(3);
    n_pt = 0;
    while (i_flag != EOF) {                             /* Plot each state */
        block_r(&x, &y, &i_flag, &n_pt);
        if (error_number)
            goto Error;
        if (i_flag != EOF) {
            plot_block(x, y, n_pt);
            if (error_number)
                goto Error;
        }
        else {
            block_r(&x, &y, &i_flag, &n_pt);                /* Free memory */
            if (error_number)
                goto Error;
            break;
        }
    }
    color(BLACK);
    i_flag = 0;
    n_pt = 0;
    while (i_flag != EOF) {                             /* Plot each state */
        block_r(&x, &y, &i_flag, &n_pt);
        if (error_number)
            goto Error;
        if (i_flag != EOF) {
            plot_block(x, y, n_pt);
            if (error_number)
                goto Error;
        }
        else {
            block_r(&x, &y, &i_flag, &n_pt);                /* Free memory */
            if (error_number)
                goto Error;
            break;
        }
    }
    tcurve(1);
Error:
    endplot();
EndOfApp:
    stopplot();

    if (error_number) {
        switch (error_number) {
            case 1:
                GPC_PUTS("TZTEST: Could not open map data file");
                break;

            case 2:
                GPC_PUTS("TZTEST: Could not allocate map data arrays");
                break;

            case 3:
                GPC_PUTS("TZTEST: Not enough points in a segment");
                break;

            case 4:
                GPC_PUTS("TZTEST: Could not open density data file");
                break;

            case 5:
                GPC_PUTS("TZTEST: Could not allocate density data arrays");
                break;
        }
    }
}

/*********************************************************************/
/* Calculate map parameters */
void map_x_y(float *parea_x, float *parea_y)
{
    float a_ratio;                                  /* aspect ratio of map */
    float area_x, area_y;              /* Plot area to preserve dimensions */
    float pagexy;                                     /* Page aspect ratio */

    pagexy = (float)(9.00 / 6.884);
    a_ratio = (float)((LONGW - LONGE) * cos((LATS + (LATN - LATS) / 2.)
        * 2. * PI / 360.) / (LATN-LATS));
    if (a_ratio >= pagexy) {
        area_x = 8.5f;
        area_y = area_x/a_ratio;
    }
    else {
        area_y = 6.884f;
        area_x = area_y * a_ratio;
    }
    *parea_x = area_x;
    *parea_y = area_y;
}

/*********************************************************************/
/* Read a state outline from the map data file */
void block_r(float **px, float **py, int *pi_flag, int *pn_pt)
{
    int i, j;                            /* dummy ints to clear state id's */
    int i_flag, n_pt;
    static float *x = NULL, *y = NULL;

    if (*pi_flag == 0) {                             /* First time through */
        in_ptr = gfopen("state.chn", "r", 1);
        if (in_ptr == NULL) {
            error_number = 1;
            return;
        }
    }
    else {
        gpcfree((void **)&x);
        gpcfree((void **)&y);
    }
    if (*pi_flag == EOF)                               /* Free memory only */
        return;

    i_flag = fscanf(in_ptr, "%d %d %d", &i, &j, &n_pt);
    *pi_flag = i_flag;
    if (i_flag == EOF) {
        gfclose(in_ptr);
        return;
    }
    x = (float *)gpcalloc(n_pt, sizeof(float));
    y = (float *)gpcalloc(n_pt, sizeof(float));
    *px = x;
    *py = y;
    *pn_pt = n_pt;
    if (x == NULL || y == NULL) {
        error_number = 2;
        return;
    }
    for (i = 0;i < n_pt;i++) {
        fscanf(in_ptr, "%e %e", &x[i], &y[i]);
        x[i] = -x[i];                             /* Longitude is negative */
    }
}

/*********************************************************************/
/* Plot a state outline */
void plot_block(float *x, float *y, int n_pt)
{
    if (n_pt < 2) {
        error_number = 3;
        return;
    }
    curve(x, y, n_pt, 0);
}

/*********************************************************************/
void plot_density(int nplot)
{
    int lon_pts, lat_pts, i, j;
    FILE *in_data;
    char density_file[80];
    char str[351], *p, seps[]=" ,\n\r\t";
    float *px = NULL, *py = NULL, *pdensity = NULL;

/*
    First open input file and read in max and number of points.
    Then allocate memory.
*/
    in_data = gfopen(density_file_name, "r", 1);
    if (in_data == NULL) {
        error_number = 4;
        return;
    }
    lon_pts = LONP;
    lat_pts = LATP;
    pdensity = (float *)gpcalloc(LATP*LONP,  sizeof(float));
    px = (float *)gpcalloc(LONP, sizeof(float));
    py = (float *)gpcalloc(LATP, sizeof(float));
    if (pdensity == NULL || px == NULL || py == NULL) {
        error_number = 5;
        return;
    }

/*
    Now read in data -- BE CAREFUL because plotting routine requires x
    (longitude) to vary most rapidly, then fill x, y vectors.
*/
    density_max = 0.0f;
    for (i = 0; i < lon_pts; i++) {
        *(px+i) = (float) (-LONGW + i * DELTA_DEGREE);
       j = 0;
        fgets(str, 350, in_data);
        p = strtok(str, seps);
        while (p != NULL) {
            *(pdensity + i + j * lon_pts) = (float) atof(p);
            if (*(pdensity+i+j*lon_pts) > density_max)
                density_max = *(pdensity+i+j*lon_pts);
            *(py+j) = (float)(LATS + j * DELTA_DEGREE);
            j++;
            p = strtok(NULL, seps);
        }
    }
    if (nplot < 2) {
        patchmat(pdensity, px, lon_pts, py, lat_pts, 0.0f, 0.0f,
            pptype[nplot], intensityc, nintensityc, 0, 0.5f, 0.2f, 8.5f,
            0.8f, -5, "%-2.1f", '\301', 0.1f, ppsmooth[nplot]);
    }
    else {  /* 3-D plot */
        /* SMEAR THE DATA TWO TIMES TO MAKE A SMOOTHER FUNCTION */
        /*  NOTE: The above call to patchmat has already smoothed the
            pdensity array two times. The 3-D surface is smoothed
            two times more.
        */
        smear(pdensity, lon_pts, lat_pts, 2);
        color(BLACK);
        surmatc(pdensity, lon_pts, lat_pts, 5, 5);
    }
    color(RED);
    sprintf(density_file, "\302Maximum unsmoothed density = %3.1f",
        density_max);
    if (nplot < 2)
        pltfnt((float)(-LONGW + .5), (float)(LATS + .5), density_file, .15f, 0);
    else
        prtfnt(5.f, 5.f, density_file, .15f, 0);
    gpcfree((void **)&px);
    gpcfree((void **)&py);
    gpcfree((void **)&pdensity);
}