/*this is file contour.c */
#include "ana_structures.h"
#include <stdio.h>
#include <string.h>
#include <math.h>
extern	struct	sym_desc sym[];
extern	int	iorder, symbol_context;
extern	float	current_pen;
extern	float	wxb, wxt, wyb, wyt;
extern	int	tkplot(float x, float y, int ip);
extern	float	float_arg();
/* contour parameters and flags */
int	autocon = 1, contour_mode, contour_box, contour_flag;
int	contour_nlev = 5, contour_sym, contour_ticks = 3;
float	contour_tick_pen = 1, contour_border = 11;
float	contour_dash_lev, contour_tick_fac = 0.5;
static	char	*level_name = {"$CONTOURS"};
static	int	nx, ny;
static	float	xa, xb, ya, yb;
static	float	xsc, ysc;
/*------------------------------------------------------------------------- */
int ana_contour(narg,ps)			/* contour routine */		
/* call is CONTOUR, image, nlev, xa, xb, ya, yb */
int	narg, ps[];
{
struct	ahead	*h;
int	iq, i, nc;
float	*pf, *cl, xmax, xmin, *pp, yq, xspace, xq;
iq = ps[0];
if ( sym[iq].class != 4 ) return execute_error(66);
iq = ana_float(1, &iq);
h = (struct ahead *) sym[iq].spec.array.ptr;
ny = 1;		nx = h->dims[0];	if (h->ndim != 1) ny = h->dims[1];
pf = (float *) ((char *) h + sizeof( struct ahead ));
if (narg > 1) contour_nlev = int_arg( ps[1]);
if (contour_nlev <= 0 || contour_nlev > 50)
  { printf("invalid # of contour levels = %d\n",contour_nlev); return -1; }
xa = wxb; xb = wxt; ya = wyb; yb = wyt;			/*default window */
if (narg > 2) xa = float_arg( ps[2]);
if (narg > 3) xb = float_arg( ps[3]);
if (narg > 4) ya = float_arg( ps[4]);
if (narg > 5) yb = float_arg( ps[5]);
/* check if $contours array exists yet */
if (contour_sym == 0)  { i = symbol_context;  symbol_context = 0;
	contour_sym = find_sym(level_name);	symbol_context = i; }
						/* check for auto mode */
if (autocon) {					/* compute levels */
	redef_array(contour_sym, 3, 1, &contour_nlev);
	h = (struct ahead *) sym[contour_sym].spec.array.ptr;
	cl = (float *) ((char *) h + sizeof( struct ahead ));
	pp = pf;	nc = nx * ny - 1;	xmin = xmax = *pp++;
	while (nc--) {	if (*pp > xmax)  xmax = *pp;
			if (*pp < xmin)  xmin = *pp;	pp++; }
	yq = xmax - xmin;	pp = cl;
	switch (contour_mode) {
	case 0:  					/* linear */
	xspace = yq / contour_nlev;	*pp++ = xmin + 0.5 * xspace;  break;
	case 1:  					/* sqrt */
	xspace = sqrt(yq) / contour_nlev;	*pp++ = 0.0;  break;
	case 2:  					/* log */
	xq = MAX( xmin, 1.0);	xspace = log(yq) / contour_nlev;
	*pp++ = log(xq) + 0.5 * xspace;  break;
	}
	nc = contour_nlev - 1;
	while (nc--) { *pp = *(pp - 1) + xspace;  pp++; }
	nc = contour_nlev;	pp = cl;
	switch (contour_mode) {
	case 1:  					/* sqrt */
	while (nc--) { *pp = pow(*pp, 2.0) + xmin;  pp++; } break;
	case 2:  					/* log */
	while (nc--) { *pp =exp(*pp) + xmin;  pp++; } break;
	}
} else {				/* not auto mode, check $contours */
	if ( sym[contour_sym].class != 4 ) return execute_error(109);
	iq = ana_float(1, &contour_sym);
	h = (struct ahead *) sym[iq].spec.array.ptr;
	if (h->ndim != 1) return execute_error(109);
	cl = (float *) ((char *) h + sizeof( struct ahead ));
	contour_nlev = MIN(contour_nlev, h->dims[0]);  /* min of these used */
}
xsc = (xb - xa)/nx;	ysc = (yb - ya)/ny;
					/* do we want border and ticks? */
if (contour_flag == 0)	{	/* set flag means none, set via tvcon */
if (contour_border > 0) box();
if (contour_ticks > 0) ticks();
}
#if NeXT
anacon( pf, &nx, &ny, cl, &contour_nlev, &iorder, &xa, &ya, &xb, &yb,
	&contour_dash_lev);
#else
anacon_( pf, &nx, &ny, cl, &contour_nlev, &iorder, &xa, &ya, &xb, &yb,
	&contour_dash_lev);
#endif
return 1;
}
/*------------------------------------------------------------------------- */
int box()
{
float	oldpen, pen;
oldpen = current_pen;
pen = contour_border;
set_pen(pen);
tkplot(xa, ya, 0);	tkplot(xb, ya, 1);	tkplot(xb, yb, 1);
tkplot(xa, yb, 1);	tkplot(xa, ya, 1);
set_pen(oldpen);
return 1;
}
/*------------------------------------------------------------------------- */
int ticks()
{
int	nxx, nyy;
float	oldpen, pen;
float	tx, ty, y, x, xs, ys;
if (contour_tick_pen <= 0) return 1;
oldpen = current_pen;
pen = contour_tick_pen;
set_pen(pen);
if (iorder%2) { xs = ysc; ys = xsc;} else { xs = xsc;  ys = ysc;}
tx = xsc * contour_tick_fac;		ty = ysc * contour_tick_fac;
tx = ty = MIN(tx, ty);				/* use smallest of the 2 */
							/* do the sides */
if ( iorder < 2 || iorder > 5) {
y = ya;
printf("xs, ys, iorder = %f %f %d\n", xs,ys,iorder);
while (y < yb)
 { tkplot(xa,y,0); tkplot(xa+tx,y,1); tkplot(xb,y,0); tkplot(xb-tx,y,1);
 y += ys; }
} else {
y = yb;
while (y > ya)
 { tkplot(xa,y,0); tkplot(xa+tx,y,1); tkplot(xb,y,0); tkplot(xb-tx,y,1);
 y -= ys; }
}
							/* top and bottom */
if (iorder < 4 ) {
x = xa;
while (x < xb)
 { tkplot(x,ya,0); tkplot(x,ya+ty,1); tkplot(x,yb,0); tkplot(x,yb-ty,1);
 x += xs; }
} else {
y = yb;
while (y > ya)
 { tkplot(x,ya,0); tkplot(x,ya+ty,1); tkplot(x,yb,0); tkplot(x,yb-ty,1);
 x -= xs; }
}
				/* now draw the origin box if wanted */
if ( contour_box) {
	switch	(iorder/2) {
	case 0: tkplot(xa,ya,-4); break;
	case 1: tkplot(xa,yb,-4); break;
	case 2: tkplot(xb,yb,-4); break;
	case 3: tkplot(xb,ya,-4); break;
}
}
set_pen(oldpen);
return 1;
}
/*------------------------------------------------------------------------- */
#if NeXT
int tkdash(aa,bb,ndsh,ntimes)
#else
int tkdash_(aa,bb,ndsh,ntimes)
#endif
/* actually no dashed line support yet 3/9/92 */
float	*aa,*bb;
int	*ndsh, *ntimes;
{
int	nc;
float	x, y;
nc = *ntimes;
while (nc-- > 0) {
x = *aa++;	y = *bb++;	tkplot(x,y,0);
x = *aa++;	y = *bb++;	tkplot(x,y,1);
}
}
/*------------------------------------------------------------------------- */