#ifndef lint

static char *RCSid() { return RCSid("$Id: plot3d.c,v 1.278 2017-09-18 22:24:03 sfeam Exp $"); }

#include "plot3d.h"

#include "gp_types.h"

#include "alloc.h"

#include "axis.h"

#include "command.h"

#include "contour.h"

#include "datafile.h"

#include "eval.h"

#include "graph3d.h"

#include "hidden3d.h"

#include "misc.h"

#include "parse.h"

#include "pm3d.h"

#include "setshow.h"

#include "term_api.h"

#include "tabulate.h"

#include "util.h"

#include "variable.h" /* For locale handling */

#include "plot2d.h" /* Only for store_label() */

#include "matrix.h" /* Used by thin-plate-splines in dgrid3d */

#include "help.h"

t_data_mapping mapping3d = MAP3D_CARTESIAN;

int dgrid3d_row_fineness = 10;

int dgrid3d_col_fineness = 10;

int dgrid3d_norm_value = 1;

int dgrid3d_mode = DGRID3D_QNORM;

double dgrid3d_x_scale = 1.0;

double dgrid3d_y_scale = 1.0;

TBOOLEAN dgrid3d = FALSE;

TBOOLEAN dgrid3d_kdensity = FALSE;

static void calculate_set_of_isolines __PROTO((AXIS_INDEX value_axis, TBOOLEAN cross, struct iso_curve **this_iso,

AXIS_INDEX iso_axis, double iso_min, double iso_step, int num_iso_to_use,

AXIS_INDEX sam_axis, double sam_min, double sam_step, int num_sam_to_use

));

static int get_3ddata __PROTO((struct surface_points * this_plot));

static void eval_3dplots __PROTO((void));

static void grid_nongrid_data __PROTO((struct surface_points * this_plot));

static void parametric_3dfixup __PROTO((struct surface_points * start_plot, int *plot_num));

static struct surface_points * sp_alloc __PROTO((int num_samp_1, int num_iso_1, int num_samp_2, int num_iso_2));

static void sp_replace __PROTO((struct surface_points *sp, int num_samp_1, int num_iso_1, int num_samp_2, int num_iso_2));

static double splines_kernel __PROTO((double h));

static void thin_plate_splines_setup __PROTO(( struct iso_curve *old_iso_crvs, double **p_xx, int *p_numpoints ));

static double qnorm __PROTO(( double dist_x, double dist_y, int q ));

static double pythag __PROTO(( double dx, double dy ));

static void count_3dpoints __PROTO((struct surface_points *plot, int *nt, int *ni, int *nu));

static void load_contour_label_options __PROTO((struct text_label *contour_label));

struct surface_points *first_3dplot = NULL;

static struct udft_entry plot_func;

int plot3d_num=0;

static struct surface_points *

sp_alloc(int num_samp_1, int num_iso_1, int num_samp_2, int num_iso_2)

{

struct lp_style_type default_lp_properties = DEFAULT_LP_STYLE_TYPE;

struct surface_points *sp = gp_alloc(sizeof(*sp), "surface");

memset(sp,0,sizeof(struct surface_points));

/* Initialize various fields */

sp->lp_properties = default_lp_properties;

sp->fill_properties = default_fillstyle;

if (sp->fill_properties.fillstyle == FS_EMPTY)

sp->fill_properties.fillstyle = FS_SOLID;

default_arrow_style(&(sp->arrow_properties));

if (num_iso_2 > 0 && num_samp_1 > 0) {

int i;

struct iso_curve *icrv;

for (i = 0; i < num_iso_1; i++) {

icrv = iso_alloc(num_samp_2);

icrv->next = sp->iso_crvs;

sp->iso_crvs = icrv;

}

for (i = 0; i < num_iso_2; i++) {

icrv = iso_alloc(num_samp_1);

icrv->next = sp->iso_crvs;

sp->iso_crvs = icrv;

}

}

return (sp);

}

static void

sp_replace(

struct surface_points *sp,

int num_samp_1, int num_iso_1, int num_samp_2, int num_iso_2)

{

int i;

struct iso_curve *icrv, *icrvs = sp->iso_crvs;

while (icrvs) {

icrv = icrvs;

icrvs = icrvs->next;

iso_free(icrv);

}

sp->iso_crvs = NULL;

if (num_iso_2 > 0 && num_samp_1 > 0) {

for (i = 0; i < num_iso_1; i++) {

icrv = iso_alloc(num_samp_2);

icrv->next = sp->iso_crvs;

sp->iso_crvs = icrv;

}

for (i = 0; i < num_iso_2; i++) {

icrv = iso_alloc(num_samp_1);

icrv->next = sp->iso_crvs;

sp->iso_crvs = icrv;

}

} else

sp->iso_crvs = NULL;

}

sp_free(struct surface_points *sp)

{

while (sp) {

struct surface_points *next = sp->next_sp;

free(sp->title);

free(sp->title_position);

sp->title_position = NULL;

while (sp->contours) {

struct gnuplot_contours *next_cntrs = sp->contours->next;

free(sp->contours->coords);

free(sp->contours);

sp->contours = next_cntrs;

}

while (sp->iso_crvs) {

struct iso_curve *next_icrvs = sp->iso_crvs->next;

iso_free(sp->iso_crvs);

sp->iso_crvs = next_icrvs;

}

if (sp->labels) {

free_labels(sp->labels);

sp->labels = (struct text_label *)NULL;

}

free(sp);

sp = next;

}

}

plot3drequest()

{

int dummy_token0 = -1, dummy_token1 = -1;

AXIS_INDEX axis, u_axis, v_axis;

is_3d_plot = TRUE;

if (parametric && strcmp(set_dummy_var[0], "t") == 0) {

strcpy(set_dummy_var[0], "u");

strcpy(set_dummy_var[1], "v");

}

/* put stuff into arrays to simplify access */

AXIS_INIT3D(FIRST_X_AXIS, 0, 0);

AXIS_INIT3D(FIRST_Y_AXIS, 0, 0);

AXIS_INIT3D(FIRST_Z_AXIS, 0, 1);

AXIS_INIT3D(U_AXIS, 1, 0);

AXIS_INIT3D(V_AXIS, 1, 0);

AXIS_INIT3D(COLOR_AXIS, 0, 1);

/* Nonlinear mapping of x or y via linkage to a hidden primary axis. */

/* The user set autoscale for the visible axis; apply it also to the hidden axis. */

for (axis = 0; axis < NUMBER_OF_MAIN_VISIBLE_AXES; axis++) {

AXIS *secondary = &axis_array[axis];

if (axis == SAMPLE_AXIS)

continue;

if (secondary->linked_to_primary

&& secondary->linked_to_primary->index == -secondary->index) {

AXIS *primary = secondary->linked_to_primary;

primary->set_autoscale = secondary->set_autoscale;

axis_init(primary, 1);

}

}

if (!term) /* unknown */

int_error(c_token, "use 'set term' to set terminal type first");

/* Range limits for the entire plot are optional but must be given */

/* in a fixed order. The keyword 'sample' terminates range parsing. */

u_axis = (parametric ? U_AXIS : FIRST_X_AXIS);

v_axis = (parametric ? V_AXIS : FIRST_Y_AXIS);

dummy_token0 = parse_range(u_axis);

dummy_token1 = parse_range(v_axis);

if (parametric) {

parse_range(FIRST_X_AXIS);

parse_range(FIRST_Y_AXIS);

}

parse_range(FIRST_Z_AXIS);

check_axis_reversed(FIRST_X_AXIS);

check_axis_reversed(FIRST_Y_AXIS);

check_axis_reversed(FIRST_Z_AXIS);

if (equals(c_token,"sample") && equals(c_token+1,"["))

c_token++;

/* Clear out any tick labels read from data files in previous plot */

for (axis=0; axis<AXIS_ARRAY_SIZE; axis++) {

struct ticdef *ticdef = &axis_array[axis].ticdef;

if (ticdef->def.user)

ticdef->def.user = prune_dataticks(ticdef->def.user);

if (!ticdef->def.user && ticdef->type == TIC_USER)

ticdef->type = TIC_COMPUTED;

}

/* use the default dummy variable unless changed */

if (dummy_token0 > 0)

copy_str(c_dummy_var[0], dummy_token0, MAX_ID_LEN);

else

strcpy(c_dummy_var[0], set_dummy_var[0]);

if (dummy_token1 > 0)

copy_str(c_dummy_var[1], dummy_token1, MAX_ID_LEN);

else

strcpy(c_dummy_var[1], set_dummy_var[1]);

/* In "set view map" mode the x2 and y2 axes are legal */

/* but must be linked to the respective primary axis. */

if (splot_map) {

if ((axis_array[SECOND_X_AXIS].ticmode && !axis_array[SECOND_X_AXIS].linked_to_primary)

|| (axis_array[SECOND_Y_AXIS].ticmode && !axis_array[SECOND_Y_AXIS].linked_to_primary))

int_error(NO_CARET,

"Secondary axis must be linked to primary axis in order to draw tics");

}

eval_3dplots();

}

refresh_3dbounds(struct surface_points *first_plot, int nplots)

{

struct surface_points *this_plot = first_plot;

int iplot; /* plot index */

for (iplot = 0; iplot < nplots; iplot++, this_plot = this_plot->next_sp) {

int i; /* point index */

struct axis *x_axis = &axis_array[FIRST_X_AXIS];

struct axis *y_axis = &axis_array[FIRST_Y_AXIS];

struct axis *z_axis = &axis_array[FIRST_Z_AXIS];

struct iso_curve *this_curve;

/* IMAGE clipping is done elsewhere, so we don't need INRANGE/OUTRANGE checks */

if (this_plot->plot_style == IMAGE

|| this_plot->plot_style == RGBIMAGE

|| this_plot->plot_style == RGBA_IMAGE) {

if (x_axis->set_autoscale)

process_image(this_plot, IMG_UPDATE_AXES);

continue;

}

for ( this_curve = this_plot->iso_crvs; this_curve; this_curve = this_curve->next) {

for (i=0; i<this_curve->p_count; i++) {

struct coordinate GPHUGE *point = &this_curve->points[i];

if (point->type == UNDEFINED)

continue;

else

point->type = INRANGE;

/* If the state has been set to autoscale since the last plot,

/* This autoscaling logic is parallel to that in

if (!this_plot->noautoscale) {

autoscale_one_point(x_axis, point->x);

autoscale_one_point(y_axis, point->y);

}

if (!inrange(point->x, x_axis->min, x_axis->max)

|| !inrange(point->y, y_axis->min, y_axis->max)) {

point->type = OUTRANGE;

continue;

}

if (!this_plot->noautoscale) {

autoscale_one_point(z_axis, point->z);

}

if (!inrange(point->z, z_axis->min, z_axis->max)) {

point->type = OUTRANGE;

continue;

}

} /* End of this curve */

}

} /* End of this plot */

/* handle 'reverse' ranges */

axis_check_range(FIRST_X_AXIS);

axis_check_range(FIRST_Y_AXIS);

axis_check_range(FIRST_Z_AXIS);

/* Make sure the bounds are reasonable, and tweak them if they aren't */

axis_checked_extend_empty_range(FIRST_X_AXIS, NULL);

axis_checked_extend_empty_range(FIRST_Y_AXIS, NULL);

axis_checked_extend_empty_range(FIRST_Z_AXIS, NULL);

}

static double

splines_kernel(double h)

{

if (h > 0.0) { return h * h * log(h); }

return 0.0;

}

static void

thin_plate_splines_setup( struct iso_curve *old_iso_crvs,

double **p_xx, int *p_numpoints )

{

int i, j, k;

double *xx, *yy, *zz, *b, **K, d;

int numpoints, *indx;

struct iso_curve *oicrv;

numpoints = 0;

for (oicrv = old_iso_crvs; oicrv != NULL; oicrv = oicrv->next) {

numpoints += oicrv->p_count;

}

xx = gp_alloc(sizeof(xx[0]) * (numpoints + 3) * (numpoints + 8),

"thin plate splines in dgrid3d");

/* the memory needed is not really (n+3)*(n+8) for now,

K = gp_alloc(sizeof(K[0]) * (numpoints + 3),

"matrix : thin plate splines 2d");

yy = xx + numpoints;

zz = yy + numpoints;

b = zz + numpoints;

/* HBB 20010424: Count actual input points without the UNDEFINED

numpoints = 0;

for (oicrv = old_iso_crvs; oicrv != NULL; oicrv = oicrv->next) {

struct coordinate GPHUGE *opoints = oicrv->points;

for (k = 0; k < oicrv->p_count; k++, opoints++) {

/* HBB 20010424: avoid crashing for undefined input */

if (opoints->type == UNDEFINED)

continue;

xx[numpoints] = opoints->x;

yy[numpoints] = opoints->y;

zz[numpoints] = opoints->z;

numpoints++;

}

}

for (i = 0; i < numpoints + 3; i++) {

K[i] = b + (numpoints + 3) * (i + 1);

}

for (i = 0; i < numpoints; i++) {

for (j = i + 1; j < numpoints; j++) {

double dx = xx[i] - xx[j], dy = yy[i] - yy[j];

K[i][j] = K[j][i] = -splines_kernel(sqrt(dx * dx + dy * dy));

}

K[i][i] = 0.0; /* here will come the weights for errors */

b[i] = zz[i];

}

for (i = 0; i < numpoints; i++) {

K[i][numpoints] = K[numpoints][i] = 1.0;

K[i][numpoints + 1] = K[numpoints + 1][i] = xx[i];

K[i][numpoints + 2] = K[numpoints + 2][i] = yy[i];

}

b[numpoints] = 0.0;

b[numpoints + 1] = 0.0;

b[numpoints + 2] = 0.0;

K[numpoints][numpoints] = 0.0;

K[numpoints][numpoints + 1] = 0.0;

K[numpoints][numpoints + 2] = 0.0;

K[numpoints + 1][numpoints] = 0.0;

K[numpoints + 1][numpoints + 1] = 0.0;

K[numpoints + 1][numpoints + 2] = 0.0;

K[numpoints + 2][numpoints] = 0.0;

K[numpoints + 2][numpoints + 1] = 0.0;

K[numpoints + 2][numpoints + 2] = 0.0;

indx = gp_alloc(sizeof(indx[0]) * (numpoints + 3), "indexes lu");

/* actually, K is *not* positive definite, but

lu_decomp(K, numpoints + 3, indx, &d);

lu_backsubst(K, numpoints + 3, indx, b);

free( K );

free( indx );

*p_xx = xx;

*p_numpoints = numpoints;

}

static double

qnorm( double dist_x, double dist_y, int q )

{

double dist = 0.0;

switch (q) {

case 1:

dist = dist_x + dist_y;

break;

case 2:

dist = dist_x * dist_x + dist_y * dist_y;

break;

case 4:

dist = dist_x * dist_x + dist_y * dist_y;

dist *= dist;

break;

case 8:

dist = dist_x * dist_x + dist_y * dist_y;

dist *= dist;

dist *= dist;

break;

case 16:

dist = dist_x * dist_x + dist_y * dist_y;

dist *= dist;

dist *= dist;

dist *= dist;

break;

default:

dist = pow(dist_x, (double)q ) + pow(dist_y, (double)q );

break;

}

return dist;

}

static double

pythag( double dx, double dy )

{

double x, y;

x = fabs(dx);

y = fabs(dy);

if (x > y) { return x*sqrt(1.0 + (y*y)/(x*x)); }

if (y==0.0) { return 0.0; }

return y*sqrt(1.0 + (x*x)/(y*y));

}

static void

grid_nongrid_data(struct surface_points *this_plot)

{

int i, j, k;

double x, y, z, w, dx, dy, xmin, xmax, ymin, ymax;

struct iso_curve *old_iso_crvs = this_plot->iso_crvs;

struct iso_curve *icrv, *oicrv, *oicrvs;

/* these are only needed for thin_plate_splines */

double *xx = NULL, *yy = NULL, *zz = NULL, *b = NULL;

int numpoints = 0;

xx = NULL; /* save to call free() on NULL if xx has never been used */

/* Compute XY bounding box on the original data. */

/* FIXME HBB 20010424: Does this make any sense? Shouldn't we just

xmin = xmax = old_iso_crvs->points[0].x;

ymin = ymax = old_iso_crvs->points[0].y;

for (icrv = old_iso_crvs; icrv != NULL; icrv = icrv->next) {

struct coordinate GPHUGE *points = icrv->points;

for (i = 0; i < icrv->p_count; i++, points++) {

/* HBB 20010424: avoid crashing for undefined input */

if (points->type == UNDEFINED)

continue;

if (xmin > points->x)

xmin = points->x;

if (xmax < points->x)

xmax = points->x;

if (ymin > points->y)

ymin = points->y;

if (ymax < points->y)

ymax = points->y;

}

}

dx = (xmax - xmin) / (dgrid3d_col_fineness - 1);

dy = (ymax - ymin) / (dgrid3d_row_fineness - 1);

/* Create the new grid structure, and compute the low pass filtering from

this_plot->iso_crvs = NULL;

this_plot->num_iso_read = dgrid3d_col_fineness;

this_plot->has_grid_topology = TRUE;

if (dgrid3d_mode == DGRID3D_SPLINES) {

thin_plate_splines_setup( old_iso_crvs, &xx, &numpoints );

yy = xx + numpoints;

zz = yy + numpoints;

b = zz + numpoints;

}

for (i = 0, x = xmin; i < dgrid3d_col_fineness; i++, x += dx) {

struct coordinate GPHUGE *points;

icrv = iso_alloc(dgrid3d_row_fineness + 1);

icrv->p_count = dgrid3d_row_fineness;

icrv->next = this_plot->iso_crvs;

this_plot->iso_crvs = icrv;

points = icrv->points;

for (j=0, y=ymin; j<dgrid3d_row_fineness; j++, y+=dy, points++) {

z = w = 0.0;

/* as soon as ->type is changed to UNDEFINED, break out of

points->type = INRANGE;

if (dgrid3d_mode == DGRID3D_SPLINES) {

z = b[numpoints];

for (k = 0; k < numpoints; k++) {

double dx = xx[k] - x, dy = yy[k] - y;

z = z - b[k] * splines_kernel(sqrt(dx * dx + dy * dy));

}

z = z + b[numpoints + 1] * x + b[numpoints + 2] * y;

} else { /* everything, except splines */

for (oicrv = old_iso_crvs; oicrv != NULL; oicrv = oicrv->next) {

struct coordinate GPHUGE *opoints = oicrv->points;

for (k = 0; k < oicrv->p_count; k++, opoints++) {

if (dgrid3d_mode == DGRID3D_QNORM) {

double dist = qnorm( fabs(opoints->x - x),

fabs(opoints->y - y),

dgrid3d_norm_value );

if (dist == 0.0) {

/* HBB 981209: revised flagging as undefined */

/* Supporting all those infinities on various

points->type = UNDEFINED;

z = opoints->z;

w = 1.0;

break; /* out of inner loop */

} else {

z += opoints->z / dist;

w += 1.0/dist;

}

} else { /* ALL else: not spline, not qnorm! */

double weight = 0.0;

double dist=pythag((opoints->x-x)/dgrid3d_x_scale,

(opoints->y-y)/dgrid3d_y_scale);

if (dgrid3d_mode == DGRID3D_GAUSS) {

weight = exp( -dist*dist );

} else if (dgrid3d_mode == DGRID3D_CAUCHY) {

weight = 1.0/(1.0 + dist*dist );

} else if (dgrid3d_mode == DGRID3D_EXP) {

weight = exp( -dist );

} else if (dgrid3d_mode == DGRID3D_BOX) {

weight = (dist<1.0) ? 1.0 : 0.0;

} else if (dgrid3d_mode == DGRID3D_HANN) {

if (dist < 1.0) {

weight = 0.5*(1-cos(2.0*M_PI*dist));

}

}

z += opoints->z * weight;

w += weight;

}

}

/* PKJ: I think this is only relevant for qnorm */

if (points->type != INRANGE)

break; /* out of the second-inner loop as well ... */

}

} /* endif( dgrid3d_mode == DGRID3D_SPLINES ) */

/* Now that we've escaped the loops safely, we know that we

points->type = INRANGE;

points->x = x;

points->y = y;

/* Honor requested x and y limits */

/* Historical note: This code was not in 4.0 or 4.2. It imperfectly */

/* restores the clipping behaviour of version 3.7 and earlier. */

if ((x < axis_array[x_axis].min && !(axis_array[x_axis].autoscale & AUTOSCALE_MIN))

|| (x > axis_array[x_axis].max && !(axis_array[x_axis].autoscale & AUTOSCALE_MAX))

|| (y < axis_array[y_axis].min && !(axis_array[y_axis].autoscale & AUTOSCALE_MIN))

|| (y > axis_array[y_axis].max && !(axis_array[y_axis].autoscale & AUTOSCALE_MAX)))

points->type = OUTRANGE;

if (dgrid3d_mode != DGRID3D_SPLINES && !dgrid3d_kdensity)

z = z / w;

STORE_AND_UPDATE_RANGE(points->z, z, points->type, z_axis,

this_plot->noautoscale, continue);

if (this_plot->pm3d_color_from_column)

int_error(NO_CARET,

"Gridding of the color column is not implemented");

else {

coord_type dummy_type = points->type;

STORE_AND_UPDATE_RANGE( points->CRD_COLOR, z, dummy_type,

COLOR_AXIS,

this_plot->noautoscale, continue);

}

}

}

free(xx); /* save to call free on NULL pointer if splines not used */

/* Delete the old non grid data. */

for (oicrvs = old_iso_crvs; oicrvs != NULL;) {

oicrv = oicrvs;

oicrvs = oicrvs->next;

iso_free(oicrv);

}

}

static int

get_3ddata(struct surface_points *this_plot)

{

int xdatum = 0;

int ydatum = 0;

int j;

int pt_in_iso_crv = 0;

struct iso_curve *this_iso;

int retval = 0;

double v[MAXDATACOLS];

/* Initialize the space that will hold input data values */

memset(v, 0, sizeof(v));

if (mapping3d == MAP3D_CARTESIAN) {

/* do this check only, if we have PM3D / PM3D-COLUMN not compiled in */

if (df_no_use_specs == 2)

int_error(this_plot->token, "Need 1 or 3 columns for cartesian data");

/* HBB NEW 20060427: if there's only one, explicit using

if (df_no_use_specs == 1)

df_axis[0] = FIRST_Z_AXIS;

} else {

if (df_no_use_specs == 1)

int_error(this_plot->token, "Need 2 or 3 columns for polar data");

}

this_plot->num_iso_read = 0;

this_plot->has_grid_topology = TRUE;

this_plot->pm3d_color_from_column = FALSE;

/* we ought to keep old memory - most likely case

if (this_plot->iso_crvs != NULL) {

struct iso_curve *icrv, *icrvs = this_plot->iso_crvs;

while (icrvs) {

icrv = icrvs;

icrvs = icrvs->next;

iso_free(icrv);

}

this_plot->iso_crvs = NULL;

}

/* data file is already open */

if (df_matrix)

this_plot->has_grid_topology = TRUE;

{

/*{{{ read surface from text file */

struct iso_curve *local_this_iso = iso_alloc(samples_1);

struct coordinate GPHUGE *cp;

struct coordinate GPHUGE *cphead = NULL; /* Only for VECTOR plots */

double x, y, z;

double xtail, ytail, ztail;

double zlow = VERYLARGE, zhigh = -VERYLARGE;

double color = VERYLARGE;

int pm3d_color_from_column = FALSE;

#define color_from_column(x) pm3d_color_from_column = x

if (this_plot->plot_style == LABELPOINTS)

expect_string( 4 );

if (this_plot->plot_style == VECTOR) {

local_this_iso->next = iso_alloc(samples_1);

local_this_iso->next->p_count = 0;

}

/* If the user has set an explicit locale for numeric input, apply it */

/* here so that it affects data fields read from the input file. */

set_numeric_locale();

/* Initial state */

df_warn_on_missing_columnheader = TRUE;

while ((retval = df_readline(v,MAXDATACOLS)) != DF_EOF) {

j = retval;

if (j == 0) /* not blank line, but df_readline couldn't parse it */

int_warn(NO_CARET, "Bad data on line %d of file %s",

df_line_number, df_filename ? df_filename : "");

if (j == DF_SECOND_BLANK)

break; /* two blank lines */

if (j == DF_FIRST_BLANK) {

/* Images are in a sense similar to isocurves.

if ((this_plot->plot_style == IMAGE)

|| (this_plot->plot_style == RGBIMAGE)

|| (this_plot->plot_style == RGBA_IMAGE))

continue;

if (this_plot->plot_style == VECTOR)

continue;

/* one blank line */

if (pt_in_iso_crv == 0) {

if (xdatum == 0)

continue;

pt_in_iso_crv = xdatum;

}

if (xdatum > 0) {

local_this_iso->p_count = xdatum;

local_this_iso->next = this_plot->iso_crvs;

this_plot->iso_crvs = local_this_iso;

this_plot->num_iso_read++;

if (xdatum != pt_in_iso_crv)

this_plot->has_grid_topology = FALSE;

local_this_iso = iso_alloc(pt_in_iso_crv);

xdatum = 0;

ydatum++;

}

continue;

}

else if (j == DF_FOUND_KEY_TITLE){

/* only the shared part of the 2D and 3D headers is used */

df_set_key_title((struct curve_points *)this_plot);

continue;

}

else if (j == DF_KEY_TITLE_MISSING){

fprintf(stderr,

"get_data: key title not found in requested column\n"

);

continue;

}

else if (j == DF_COLUMN_HEADERS) {

continue;

}

/* its a data point or undefined */

if (xdatum >= local_this_iso->p_max) {

/* overflow about to occur. Extend size of points[]

iso_extend(local_this_iso, xdatum + xdatum + 1000);

if (this_plot->plot_style == VECTOR) {

iso_extend(local_this_iso->next, xdatum + xdatum + 1000);

local_this_iso->next->p_count = 0;

}

}

cp = local_this_iso->points + xdatum;

if (this_plot->plot_style == VECTOR) {

if (j < 6) {

cp->type = UNDEFINED;

continue;

}

cphead = local_this_iso->next->points + xdatum;

}

if (j == DF_UNDEFINED || j == DF_MISSING) {

cp->type = UNDEFINED;

goto come_here_if_undefined;

}

cp->type = INRANGE; /* unless we find out different */

/* EAM Oct 2004 - Substantially rework this section */

/* now that there are many more plot types. */

x = y = z = 0.0;

xtail = ytail = ztail = 0.0;

/* The x, y, z coordinates depend on the mapping type */

switch (mapping3d) {

case MAP3D_CARTESIAN:

if (j == 1) {

x = xdatum;

y = ydatum;

z = v[0];

j = 3;

break;

}

if (j == 2) {

if (PM3DSURFACE != this_plot->plot_style)

int_error(this_plot->token,

"2 columns only possible with explicit pm3d style (line %d)",

df_line_number);

x = xdatum;

y = ydatum;

z = v[0];

color_from_column(TRUE);

color = v[1];

j = 3;

break;

}

/* Assume everybody agrees that x,y,z are the first three specs */

if (j >= 3) {

x = v[0];

y = v[1];

z = v[2];

break;

}

break;

case MAP3D_SPHERICAL:

if (j < 2)

int_error(this_plot->token, "Need 2 or 3 columns");

if (j < 3) {

v[2] = 1; /* default radius */

j = 3;

}

/* Convert to radians. */

v[0] *= ang2rad;

v[1] *= ang2rad;

x = v[2] * cos(v[0]) * cos(v[1]);

y = v[2] * sin(v[0]) * cos(v[1]);

z = v[2] * sin(v[1]);

break;

case MAP3D_CYLINDRICAL:

if (j < 2)

int_error(this_plot->token, "Need 2 or 3 columns");

if (j < 3) {

v[2] = 1; /* default radius */

j = 3;

}

/* Convert to radians. */

v[0] *= ang2rad;

x = v[2] * cos(v[0]);

y = v[2] * sin(v[0]);

z = v[1];

break;

default:

int_error(NO_CARET, "Internal error: Unknown mapping type");

return retval;

}

if (j < df_no_use_specs)

int_error(this_plot->token,

"Wrong number of columns in input data - line %d",

df_line_number);

/* Work-around for hidden3d, which otherwise would use the */

/* color of the vector midpoint rather than the endpoint. */

if (this_plot->plot_style == IMPULSES) {

if (this_plot->lp_properties.pm3d_color.type == TC_Z) {

color = z;

color_from_column(TRUE);

}

}

/* After the first three columns it gets messy because */

/* different plot styles assume different contents in the columns */

if (( this_plot->plot_style == POINTSTYLE || this_plot->plot_style == LINESPOINTS)) {

int varcol = 3;

if (this_plot->lp_properties.p_size == PTSZ_VARIABLE)

cp->CRD_PTSIZE = v[varcol++];

if (this_plot->lp_properties.p_type == PT_VARIABLE)