/*
 * plotting.c
 *
 * plotting routines.
 *
 */
#include <stdio.h>
#include <errno.h>
extern int errno;
#include <assert.h>
#include <string.h>
#include <ctype.h>
#include <signal.h>
#include "datatypes.h"

void	ErrorExit();

#ifndef	min
#	define	min(A,B)	A<B? A : B
#endif
#ifndef	max
#	define	max(A,B)	A>B? A : B
#endif

extern bool	GraphicsInProgress;

/* GoldenRatio = (1+Sqrt(5))/2	*/
#define	GOLDENRATIO	1.6180339887499
#define ASPECTRATIO	1.0/GOLDENRATIO
extern VALUE	V_AngularUnit;
extern VALUE	V_AnnotationScale;
extern VALUE	V_AspectRatio;
extern VALUE	V_Boxed;
extern VALUE	V_Domain;
extern VALUE	V_Gridding;
extern VALUE	V_LabelScale;
extern VALUE	V_LineColor;
extern VALUE	V_LineStyle;
extern VALUE	V_Orientation;
extern VALUE	V_Origin;
extern VALUE	V_PlotColor;
extern VALUE	V_PlotDevice;
extern VALUE	V_PlotJoined;
extern VALUE	V_PlotPoints;
extern VALUE	V_PlotTitle;
extern VALUE	V_PlotType;
extern VALUE	V_PointScale;
extern VALUE	V_PointSymbol;
extern VALUE	V_PolarVariable;
extern VALUE	V_Range;
extern VALUE	V_SubPages;
extern VALUE	V_SupplyAbscissa;
extern VALUE	V_TitleScale;
extern VALUE	V_UseInputFile;
extern VALUE	V_Verbose;
extern VALUE	V_ViewPort;
extern VALUE	V_XLabel;
extern VALUE	V_XTick;
extern VALUE	V_YLabel;
extern VALUE	V_YTick;

char	*DeviceTypes[] = {
	"amiga",
	"printer",
	"iff",
	"hp",
	"aegis",
	"postscript",
	(char *)NULL
};

typedef struct line_style {
		int	ls_type;
		int	ls_nelms;
		int	*ls_space;
		int	*ls_mark;
	} LINESTYLE;
LIST	LineStyleList;
LINESTYLE	*LineStyles; int NLineStyles;
int		*PointCodes, NPointCodes;
int	space0 = 0, mark0 = 0, space1 = 1500, mark1 = 1500;

RECT	UserRect;
RECT	ViewPort = { 0.1,0.1, 0.9,0.9 };
RECT	WorldRect = {0.0,0.0, 1.0,1.0};


void
ListPlot(Data, NTuples, TupleSize)
FLOAT	**Data;
int	NTuples, TupleSize;
{
	void InitializeDevice();
	void InitializePlottingEnv();
#ifdef	ANSI_C
	void Plot2D(FLOAT **, int , int );
	void PolarPlot(FLOAT **, int , int );
#else
	void Plot2D();
	void PolarPlot();
#endif

	InitializeDevice();
	GraphicsInProgress = TRUE;

	InitializePlottingEnv(Data, NTuples, TupleSize);

	if (Vstrcmp(V_PlotType,"polar")) {
		PolarPlot(Data, NTuples, TupleSize);
	} else {
		Plot2D(Data, NTuples, TupleSize);
	}
	plend();
}


void
InitializeDevice()
{
	register int	i;
	register char	*DeviceStr;


	plorient(Vstrcmp(V_Orientation, "portrait")? 1 : 0);
	plselect(stdout);	/* write to stdout */
		

	DeviceStr = V_PlotDevice.val_u.udt_string;
	for (i=0; DeviceTypes[i]; i++) 
		if (strcmp(DeviceStr, DeviceTypes[i]) == 0) {
			plbeg(
				i+1,	/* device code	*/
				(int)V_SubPages.val_u.udt_interval.int_lo, /* nx */
				(int)V_SubPages.val_u.udt_interval.int_hi /* ny */
			);
			return;
		}
	fprintf(stderr,"(InitializeDevice) Unknown device type: %s\n", DeviceStr);
	ErrorExit();
}

void
InitializePlottingEnv(Data, NTuples, TupleSize)
FLOAT	**Data;
int	NTuples;
int	TupleSize;
{
	int	i, j;
	FLOAT	*Ptr;
	void	InitLineStyles();
	void	InitPointCodes();

	/* viewport	*/
	ViewPort = VtoRect(V_ViewPort);

	/* user window		*/
	UserRect.XMin = MAX_FLOAT;
	UserRect.YMin = MAX_FLOAT;
	UserRect.XMax = -MAX_FLOAT;
	UserRect.YMax = -MAX_FLOAT;
	
	if (VtoBoolean(V_SupplyAbscissa)) {
		UserRect.XMin = 0.0;
		UserRect.XMax = (FLOAT)(NTuples-1);
	} else {
		for (i=0; i<NTuples; i++) {
			if (Data[i][0] < UserRect.XMin)
				UserRect.XMin = Data[i][0];
			if (Data[i][0] > UserRect.XMax)
				UserRect.XMax = Data[i][0];
		}
	}
	for (i=0; i<NTuples; i++) 
		for (j=TupleSize-1,Ptr= &Data[i][1]; j>0; --j,Ptr++) {
			if (*Ptr < UserRect.YMin)
				UserRect.YMin = *Ptr;
			if (*Ptr > UserRect.YMax)
				UserRect.YMax = *Ptr;
		}
		
	
	/* font scaling		*/

	/* plot type		*/

	/* grid			*/

	/* line styles		*/
	InitLineStyles();
	

	/* line colors		*/

	/* point symbols	*/
	InitPointCodes();

	/* window type		*/

	/* X label		*/

	/* Y label		*/

	/* Plot title		*/

	/* aspect ratio		*/
}

#define	isMark(C)	(((C) == 'M' || (C) == 'm')? TRUE : FALSE)
#define	isSpace(C)	(((C) == 'S' || (C) == 's')? TRUE : FALSE)
#define	MarkLength(C)	(C) == 'M'? 1000 : 250
#define	SpaceLength(C)	(C) == 'S'? 1000 : 250
#define	MAX_MARKS	32

void
DecodeLineStyle(StyleStr, Mark, Space, N)
char	*StyleStr;
int	**Mark;
int	**Space;
int	*N;
{
	register char	*S;
	int	macc, sacc;
	int	i, j, k;
	int	mark[MAX_MARKS];
	int	space[MAX_MARKS];
	bool	ProcessingMark;

	/*
	 * 'S' = 1000   micron space
	 * 's' =  250   micron space
	 * 'M' = 1000 	micron line
	 * 'm' =  250	micron line
	 */
	ProcessingMark = TRUE; /* to satisfy some compilers */
	if (isMark(*StyleStr))
		ProcessingMark = TRUE;
	else if (isSpace(*StyleStr))
		ProcessingMark = FALSE;
	else {
		fprintf(stderr,"(DecodeLineStyle) Invalid Line style: \"%s\"\n", StyleStr);
		ErrorExit();
	}
	memset((char *)mark, 0, MAX_MARKS * sizeof(int));
	memset((char *)space, 0, MAX_MARKS * sizeof(int));
	for (macc=0,sacc=0,i=0,j=0,S=StyleStr ; *S && i<MAX_MARKS && j<MAX_MARKS; ++S) {
		if (isMark(*S) && ProcessingMark) {
			macc += MarkLength(*S);
		} else if (isMark(*S) && !ProcessingMark) {
			space[j++] = sacc;
			macc = MarkLength(*S);
			ProcessingMark = TRUE;
		} else if (isSpace(*S) && ProcessingMark) {
			mark[i++] = macc;
			sacc = SpaceLength(*S);
			ProcessingMark = FALSE;
		} else if (isSpace(*S) && !ProcessingMark) {
			sacc += SpaceLength(*S);
		} else if (isspace(*S)) {
			continue;
		} else {
			fprintf(stderr,"(DecodeLineStyle) Invalid line style code '%c' in \"%s\"\n", *S, StyleStr);
			ErrorExit();
		}
	}
	if (ProcessingMark)
		mark[i++] = macc;
	else
		space[j++] = sacc;

	k = max(i,j);

	if (((*Mark) = (int *)calloc(k, sizeof(int))) == (int *)NULL) {
		perror("(DecodeLineStyle) ");
		ErrorExit();
	}
	if (((*Space) = (int *)calloc(k, sizeof(int))) == (int *)NULL) {
		perror("(DecodeLineStyle) ");
		ErrorExit();
	}
	for (i=0; i<k; i++) {
		(*Mark)[i] = mark[i];
		(*Space)[i] = space[i];
	}
	*N = k;
}

void
InitLineStyles()
{
	int		i;
	LATOM		*A;
	LINESTYLE	*LS;
	int		*Mark, *Space;
	int		N;
#ifdef	ANSI_C
	LINESTYLE	*LSAlloc(int , int *, int *);
#else
	LINESTYLE	*LSAlloc();
#endif

	LineStyleList.list_n = 0;
	LineStyleList.list_head = (LATOM *)NULL;
	LineStyleList.list_tail = (LATOM *)NULL;

	LS = LSAlloc(0, (int *)NULL, (int *)NULL);
	Append(&LineStyleList, (generic *)LS);

	for (i=1,A=V_LineStyle.val_u.udt_set.list_head; A; A=A->la_next,i++) {
		DecodeLineStyle(
			(char *)A->la_p,
			&Mark,
			&Space,
			&N
		);
		LS = LSAlloc(N, Mark, Space);
		Append(&LineStyleList, (generic *)LS);
	}

	NLineStyles = i;
	if ((LineStyles = (LINESTYLE *)calloc(NLineStyles, sizeof(LINESTYLE))) == (LINESTYLE *)NULL) {
		perror("(InitLineStyles) ");
		ErrorExit();
	}
	for (i=0,A=LineStyleList.list_head; A; A=A->la_next,i++) {
		LineStyles[i].ls_mark = ((LINESTYLE *)(A->la_p))->ls_mark;
		LineStyles[i].ls_space = ((LINESTYLE *)(A->la_p))->ls_space;
		LineStyles[i].ls_nelms = ((LINESTYLE *)(A->la_p))->ls_nelms;
	}
}


LINESTYLE	*
LSAlloc(n, Mark, Space)
int	n;
int	*Mark, *Space;
{
	register LINESTYLE	*LS;

	if ((LS = (LINESTYLE *)calloc(1, sizeof(LINESTYLE))) == (LINESTYLE *)NULL) {
		perror("(LSAlloc) ");
		assert(LS != (LINESTYLE *)NULL);
		exit(errno);
	}
	LS->ls_mark = Mark;
	LS->ls_space = Space;
	LS->ls_nelms = n;
	return(LS);
}

void
InitPointCodes()
{
	register LATOM	*A;
	register int	i;
	extern int	*PointCodes, NPointCodes;
	extern VALUE	V_PointSymbol;

	if (!isSet(V_PointSymbol)) {
		NPointCodes = 0;
		return;
	}

	/* count # codes	*/
	for (i=0,A=V_PointSymbol.val_u.udt_set.list_head; A; A=A->la_next, i++)
		;

	NPointCodes = i;
	if (NPointCodes == 0)
		return;

	if ((PointCodes = (int *)calloc(NPointCodes, sizeof(int))) == (int *)NULL) {
		perror("(InitPointCodes) ");
		assert(PointCodes != (int *)NULL);
		ErrorExit();
	}
	for (i=0,A=V_PointSymbol.val_u.udt_set.list_head; A; A=A->la_next, i++) {
		PointCodes[i] = atoi((char *)A->la_p);
	}
}


void
Plot2D(Data, NTuples, TupleSize)
FLOAT	**Data;
int	NTuples;
int	TupleSize;
{
	int	i, j;
	FLOAT	*X, *Y;
	FLOAT	x, y;
	FLOAT	x0, y0;
	FLOAT	xtick, ytick;
	FLOAT	Xmin,Xmax, Ymin,Ymax;
	FLOAT	MedianX, StdDevX;
	FLOAT	MedianY, StdDevY;
	int	nxsub, nysub;
	bool	AutoRange, AutoDomain;
	bool	LogX, LogY;
	char	Xopt[16], Yopt[16];
#ifdef	ANSI_C
	void	StatisticsOfX(
			FLOAT **,
			int,
			int ,
			FLOAT *,
			FLOAT *
		);
	void	StatisticsOfY(
			FLOAT **,
			int,
			int ,
			FLOAT *,
			FLOAT *
		);
	double	log10(double);
#else
	void	StatisticsOfX();
	void	StatisticsOfY();
	double	log10();
#endif

	pladv(0);
	ViewPort = VtoRect(V_ViewPort);

	if (isString(V_AspectRatio)) {
		/* automatic --> STRETCH */
		plvpor(ViewPort.XMin, ViewPort.XMax, ViewPort.YMin, ViewPort.YMax);
/*debug
		plvsta();
debug*/
	} else {
		assert(isDbl(V_AspectRatio));
		plvasp(VtoDbl(V_AspectRatio));
	}
	
	LogX = (Vstrcmp(V_PlotType, "loglin") || Vstrcmp(V_PlotType,"loglog"))?
			TRUE : FALSE;
	LogY = (Vstrcmp(V_PlotType, "linlog") || Vstrcmp(V_PlotType,"loglog"))?
			TRUE : FALSE;

	/* window bounds */
	AutoDomain = FALSE;
	if (isString(V_Domain) && Vstrcmp(V_Domain, "All")) {
		/* all points	*/
		Xmin = UserRect.XMin;
		Xmax = UserRect.XMax;
	} else if (isString(V_Domain) && Vstrcmp(V_Domain,"Automatic")) {
		/* automatic 	*/
		AutoDomain = TRUE;
		StatisticsOfX(Data, NTuples, TupleSize, &MedianX, &StdDevX);
		x = MedianX - 3.0*StdDevX;
		Xmin = UserRect.XMin < x? x : UserRect.XMin;
		x = MedianX + 3.0*StdDevX;
		Xmax = UserRect.XMax > x? x : UserRect.XMax;
	} else {
		assert(isInterval(V_Domain));
		Xmin = V_Domain.val_u.udt_interval.int_lo;
		Xmax = V_Domain.val_u.udt_interval.int_hi;
	}
	if (LogX) {
		if (Xmin == 0.0 || Xmax == 0.0) {
			/* error */
			fprintf(stderr,"ListPlot: Trying to find log10(0) in data! Stop.\n");
			ErrorExit();
		}
		Xmin = log10((double)Xmin);
		Xmax = log10((double)Xmax);
	}

	AutoRange = FALSE;
	if (isString(V_Range) && Vstrcmp(V_Range, "All")) {
		/* all */
		Ymin = UserRect.YMin;
		Ymax = UserRect.YMax;
	} else if (isString(V_Range) && Vstrcmp(V_Range,"Automatic")) {
		/* automatic 	*/
		AutoRange = TRUE;
		StatisticsOfY(Data, NTuples, TupleSize, &MedianY, &StdDevY);
		y = MedianY - 3.0*StdDevY;
		Ymin = UserRect.YMin < y? y : UserRect.YMin;
		y = MedianY + 3.0*StdDevY;
		Ymax = UserRect.YMax > y? y : UserRect.YMax;
	} else {
		assert(isInterval(V_Range));
		Ymin = V_Range.val_u.udt_interval.int_lo;
		Ymax = V_Range.val_u.udt_interval.int_hi;
	}
	if (LogY) {
		if (Ymin == 0.0 || Ymax == 0.0) {
			/* error */
			fprintf(stderr,"ListPlot: Trying to find log10(0) in data! Stop.\n");
			ErrorExit();
		}
		Ymin = log10((double)Ymin);
		Ymax = log10((double)Ymax);
	}
	plwind(
		Xmin,
		Xmax,
		Ymin,
		Ymax
	);
	plschr(0., VtoDbl(V_AnnotationScale));


	/* X and Y ticks */
	if (isString(V_XTick)) {
		/* automatic */
		xtick = 0.;
		nxsub = 0;
	} else {
		xtick = V_XTick.val_u.udt_interval.int_lo;
		if (LogX) xtick = log10((double)xtick);
		nxsub = V_XTick.val_u.udt_interval.int_hi;
	}
	if (isString(V_YTick)) {
		/* automatic */
		ytick = 0.;
		nysub = 0;
	} else {
		ytick = V_YTick.val_u.udt_interval.int_lo;
		if (LogY) ytick = log10((double)ytick);
		nysub = V_YTick.val_u.udt_interval.int_hi;
	}

	/* Axis options	*/
	Xopt[0] = (char)NULL;
	Yopt[0] = (char)NULL;
	if (isBoolean(V_Boxed) && VtoBoolean(V_Boxed)) {
		strcat(Xopt, "bc");
		strcat(Yopt, "bc");
	} else if (isString(V_Boxed)) {
		if (Vstrchr(V_Boxed, 'b'))
			strcat(Xopt, "b");
		if (Vstrchr(V_Boxed, 't'))
			strcat(Xopt, "c");
		if (Vstrchr(V_Boxed, 'r'))
			strcat(Yopt, "c");
		if (Vstrchr(V_Boxed, 'l'))
			strcat(Yopt, "b");
	} else {
		;
	}
	if (LogX)
		strcat(Xopt, "l");
	if (LogY)
		strcat(Yopt, "l");
	strcat(Xopt, "nst");
	strcat(Yopt, "nstv");

	if (isString(V_Origin) && Vstrcmp(V_Origin, "Median")) {
		if (AutoDomain == FALSE) {
			StatisticsOfX(Data,NTuples,TupleSize,&MedianX,&StdDevX);
		}
		if (AutoRange == FALSE) {
			StatisticsOfY(Data,NTuples,TupleSize,&MedianY,&StdDevY);
		}
		plaxes(
			MedianX,MedianY,
			Xopt, xtick, nxsub,
			Yopt, ytick,nysub
		);
	} else if (isInterval(V_Origin)) {
		x0 = V_Origin.val_u.udt_interval.int_lo;
		y0 = V_Origin.val_u.udt_interval.int_hi;
		plaxes(
			x0,y0,
			Xopt, xtick, nxsub,
			Yopt, ytick,nysub
		);
	} else if (isString(V_Origin) && Vstrcmp(V_Origin, "Automatic")) {
		plbox(Xopt, xtick, nxsub, Yopt, ytick,nysub);
	} else {
		plbox(Xopt, xtick, nxsub, Yopt, ytick,nysub);
	}

	/* gridding	*/
	if (VtoBoolean(V_Gridding)) {
		/* gridding on */
		plstyl(1, &mark1, &space1);
		plbox("g", xtick, nxsub, "g", ytick,nysub);
		plstyl(0, &mark0, &space0);
	}

	/* plot curves	*/
	if ((X = (FLOAT *)calloc(NTuples, sizeof(FLOAT))) == (FLOAT *)NULL) {
		perror("(Plot2D) ");
		ErrorExit();
	}
	if ((Y = (FLOAT *)calloc(NTuples, sizeof(FLOAT))) == (FLOAT *)NULL) {
		perror("(Plot2D) ");
		ErrorExit();
	}
	if (VtoBoolean(V_SupplyAbscissa)) {
		for (i=0; i<NTuples; i++)
			X[i] = LogX? log10((double)(i+1)) : (FLOAT)i;
	} else {
		for (i=0; i<NTuples; i++)
			X[i] = LogX? log10((double)Data[i][0]) : Data[i][0];
	}
	for (i=1; i<TupleSize; i++) {
		if (NLineStyles > 0) {
			/* use user-supplied patterns */
			plstyl(
				LineStyles[(i-1)%NLineStyles].ls_nelms,
				LineStyles[(i-1)%NLineStyles].ls_mark,
				LineStyles[(i-1)%NLineStyles].ls_space
			);
		} else {
			pllsty((i-1)%8 + 1);
		}
		plcol(i);
		for (j=0; j<NTuples; j++)
			Y[j] = LogY? log10((double)Data[j][i]) : Data[j][i];

		if (VtoBoolean(V_PlotJoined))
			plline(NTuples, X, Y);

		if (VtoBoolean(V_PlotPoints)) {
			plschr(0., VtoDbl(V_PointScale));
			pllsty(1);
			if (NPointCodes > 0) {
				plpoin(NTuples, X, Y, PointCodes[(i-1)%NPointCodes]);
			} else {
				plpoin(NTuples, X, Y, i);
			}
		}
	}


	/* restore line styles etc... */
	pllsty(1);
	plcol(1);
	plschr(0., VtoDbl(V_LabelScale));
	pllab(VtoString(V_XLabel), VtoString(V_YLabel), "");
	plschr(0., VtoDbl(V_TitleScale));
	pllab("","", VtoString(V_PlotTitle));
}

#define	DTR(D)	(double)((D)*3.1415927/180.0)

void
PolarPlot(Data, NTuples, TupleSize)
FLOAT	**Data;
int	NTuples;
int	TupleSize;
{
	int	i, j;
	FLOAT	*A, *R;
	FLOAT	x, y;
	FLOAT	xtick, ytick;
	FLOAT	atick, rtick;
	FLOAT	Xmin,Xmax, Ymin,Ymax;
	FLOAT	r, rmax, ang, amax;
	FLOAT	MedianX, StdDevX;
	FLOAT	MedianY, StdDevY;
	FLOAT	da;
	FLOAT	x0,y0, xn,yn;
	double	pi;
	char	s[32];
	int	nxsub, nysub, nrsub, nasub;
	bool	InRadians, XisAngular;
	bool	AutoDomain, AutoRange;
#ifdef	ANSI_C
	void	StatisticsOfX(
			FLOAT **,
			int ,
			int ,
			FLOAT *,
			FLOAT *
		);
	void	StatisticsOfY(
			FLOAT **,
			int ,
			int ,
			FLOAT *,
			FLOAT *
		);
	double	log10(double); double fabs();
	double	sqrt(), sin(), cos(), atan();
#else
	void	StatisticsOfX();
	void	StatisticsOfY();
	double	log10(); double fabs();
	double	sqrt(), sin(), cos(), atan();
#endif

	pi = 4.0 * atan((double)1.0);
	InRadians = Vstrcmp(V_AngularUnit, "degrees")? FALSE : TRUE;
	XisAngular = Vstrcmp(V_PolarVariable,"angle")? TRUE : FALSE;

	pladv(0);

	if (isString(V_AspectRatio)) {
		/* automatic --> STRETCH */
		plvasp(1.0);
	} else {
		assert(isDbl(V_AspectRatio));
		plvasp(VtoDbl(V_AspectRatio));
	}

	/* window bounds */
	AutoDomain = FALSE;
	if (isString(V_Domain) && Vstrcmp(V_Domain, "All")) {
		/* all points	*/
		Xmin = UserRect.XMin;
		Xmax = UserRect.XMax;
	} else if (isString(V_Domain) && Vstrcmp(V_Domain,"Automatic")) {
		/* automatic 	*/
		AutoDomain = TRUE;
		StatisticsOfX(Data, NTuples, TupleSize, &MedianX, &StdDevX);
		x = MedianX - 3.0*StdDevX;
		Xmin = UserRect.XMin < x? x : UserRect.XMin;
		x = MedianX + 3.0*StdDevX;
		Xmax = UserRect.XMax > x? x : UserRect.XMax;
	} else {
		assert(isInterval(V_Domain));
		Xmin = V_Domain.val_u.udt_interval.int_lo;
		Xmax = V_Domain.val_u.udt_interval.int_hi;
	}

	AutoRange = FALSE;
	if (isString(V_Range) && Vstrcmp(V_Range, "All")) {
		/* all */
		Ymin = UserRect.YMin;
		Ymax = UserRect.YMax;
	} else if (isString(V_Range) && Vstrcmp(V_Range,"Automatic")) {
		/* automatic 	*/
		AutoRange = TRUE;
		StatisticsOfY(Data, NTuples, TupleSize, &MedianY, &StdDevY);
		y = MedianY - 3.0*StdDevY;
		Ymin = UserRect.YMin < y? y : UserRect.YMin;
		y = MedianY + 3.0*StdDevY;
		Ymax = UserRect.YMax > y? y : UserRect.YMax;
	} else {
		assert(isInterval(V_Range));
		Ymin = V_Range.val_u.udt_interval.int_lo;
		Ymax = V_Range.val_u.udt_interval.int_hi;
	}
	rmax = XisAngular? Ymax : Xmax;
	plwind(
		-1.3*rmax,
		1.3*rmax,
		-1.3*rmax,
		1.3*rmax
	);


	/* X and Y ticks */
	if (isString(V_XTick)) {
		/* automatic */
		xtick = XisAngular? (InRadians? pi/6.0 : 30.0) : 0.25 * rmax;
		nxsub = 4;
	} else {
		xtick = V_XTick.val_u.udt_interval.int_lo;
		nxsub = V_XTick.val_u.udt_interval.int_hi;
	}
	if (isString(V_YTick)) {
		/* automatic */
		ytick = XisAngular? 0.25 * rmax : (InRadians? pi/6.0 : 30.0);
		nysub = 4;
	} else {
		ytick = V_YTick.val_u.udt_interval.int_lo;
		nysub = V_YTick.val_u.udt_interval.int_hi;
	}
	rtick = XisAngular? ytick : xtick;
	atick = XisAngular? xtick : ytick;
	nasub = XisAngular? nxsub : nysub;
	nrsub = XisAngular? nysub : nxsub;
	amax = InRadians? 2.0 * pi : 360.0;

	/* Axis options	*/
	if (VtoBoolean(V_Gridding)) {
		/* draw circles	*/
		plstyl(0, &mark1, space1);
		/*  circular ticks	*/
		for (i=0,r=(rtick/nrsub); r<rmax; r += (rtick/nrsub),i++) {
			da = i%nrsub? 0.05 * atick : 0.1 * atick;
			for (ang=0; ang<amax; ang += atick) {
				x0 = r * cos(InRadians? ang-da : DTR(ang-da)); 
				y0 = r * sin(InRadians? ang-da : DTR(ang-da)); 
				xn = r * cos(InRadians? ang+da : DTR(ang+da)); 
				yn = r * sin(InRadians? ang+da : DTR(ang+da)); 
				pljoin(x0,y0,xn,yn);
			}
		}
#ifdef	CIRCULAR_TICKS
		/* major circles	*/
		for (r=rtick; r<rmax; r += rtick) {
			x0 = r; y0 = 0.0;
			for (ang=1.0; ang<=360.0; ang++) {
				xn = r * cos((double)(DTR(ang))); 
				yn = r * sin((double)(DTR(ang))); 
				pljoin(x0,y0,xn,yn);
				x0 = xn; y0 = yn;
			}
		}
#endif
		/* draw spokes	*/
		for (i=0,ang=0.0; ang<amax; ang += (atick/nasub),i++) {
			xn = rmax * cos(InRadians? ang : DTR(ang));
			yn = rmax * sin(InRadians? ang : DTR(ang));
			if (i%nasub == 0) {
				pljoin((FLOAT)0.0,(FLOAT)0.0,xn,yn);
			} else {
				pljoin(xn,yn,1.05*xn,1.05*yn);
			}
		}
	}
	/* write angle labels	*/
	plschr(0., VtoDbl(V_AnnotationScale));
	j = amax / atick;
	for (ang=0.0,i=0; i<j; ang += atick,i++) {
		if (InRadians) {
			xn = rmax * cos((double)ang);
			yn = rmax * sin((double)ang);
			if (fabs(ang-(2.0*pi)) > 1.0e-2)
				sprintf(s, "%.2f #gp", (float)ang/pi);
		} else {
			xn = rmax * cos((double)DTR(ang));
			yn = rmax * sin((double)DTR(ang));
			sprintf(s, "%d", (int)(ang+0.5));
		}
		if (xn >= 0)
			plptex(xn,yn,xn,yn,-0.15, s);
		else
			plptex(xn,yn,-xn,-yn,1.15, s);
	}
	plstyl(0, &mark0, &space0);

	/* plot curves	*/
	if ((A = (FLOAT *)calloc(NTuples, sizeof(FLOAT))) == (FLOAT *)NULL) {
		perror("(Plot2D) ");
		ErrorExit();
	}
	if ((R = (FLOAT *)calloc(NTuples, sizeof(FLOAT))) == (FLOAT *)NULL) {
		perror("(Plot2D) ");
		ErrorExit();
	}
	if (VtoBoolean(V_SupplyAbscissa)) {
		for (i=0; i<NTuples; i++)
			if (XisAngular)
				A[i] = i*2.0*pi/(NTuples-1);
			else
				R[i] = i*rmax/(NTuples-1);
	} else {
		for (i=0; i<NTuples; i++)
			if (XisAngular)
				A[i] = InRadians? Data[i][0] :  DTR(Data[i][0]);
			else
				R[i] = Data[i][0];
	}
	for (i=1; i<TupleSize; i++) {
		if (NLineStyles > 0) {
			/* use user-supplied patterns */
			plstyl(
				LineStyles[(i-1)%NLineStyles].ls_nelms,
				LineStyles[(i-1)%NLineStyles].ls_mark,
				LineStyles[(i-1)%NLineStyles].ls_space
			);
		} else {
			pllsty((i-1)%8);
		}
		plcol(i);

		/* get and convert data if necessary	*/
		if (InRadians && !XisAngular)
			for (j=0; j<NTuples; j++)
				A[j] = Data[j][i];
		else if ((!InRadians) && !XisAngular)
			for (j=0; j<NTuples; j++)
				A[j] = DTR(Data[j][i]);
		else
			for (j=0; j<NTuples; j++)
				R[j] = Data[j][i];

		/* convert to cartesian	A<->X  R<->Y*/
		for (j=0; j<NTuples; j++) {
			x0 = R[j] * cos(A[j]);
			y0 = R[j] * sin(A[j]);
			A[j] = x0;
			R[j] = y0;
		}
		if (VtoBoolean(V_PlotJoined)) {
			x0 = A[0];
			y0 = R[0];
			for (j=1; j<NTuples; j++) {
				pljoin(A[j-1], R[j-1], A[j],R[j]);
			}
		}

		if (VtoBoolean(V_PlotPoints)) {
			plschr(0., VtoDbl(V_PointScale));
			pllsty(1);
			if (NPointCodes > 0) {
				plpoin(NTuples, A, R, PointCodes[(i-1)%NPointCodes]);
			} else {
				plpoin(NTuples, A, R, i);
			}
		}
	}

	/* restore line styles etc... */
	pllsty(1);
	plcol(1);
	plschr(0., VtoDbl(V_LabelScale));
	pllab(VtoString(V_XLabel), VtoString(V_YLabel), "");
	plschr(0., VtoDbl(V_TitleScale));
	pllab("","", VtoString(V_PlotTitle));
}


void
StatisticsOfX(Data, NTuples,TupleSize, Median, StdDev)
FLOAT	**Data;
int	NTuples;
int	TupleSize;
FLOAT	*Median;
FLOAT	*StdDev;
{
#ifdef	ANSI_C
	FLOAT	MedianOfX();
	FLOAT	StdDevOfX();
/*debug
	FLOAT	MedianOfX(FLOAT **Data,int NTuples,int TupleSize);
	FLOAT	StdDevOfX(FLOAT **Data,int NTuples,int TupleSize);
debug*/
#else
	FLOAT	MedianOfX();
	FLOAT	StdDevOfX();
#endif
	*StdDev = StdDevOfX(Data,NTuples,TupleSize);	
	*Median = MedianOfX(Data,NTuples,TupleSize);
}


void
StatisticsOfY(Data, NTuples,TupleSize, Median, StdDev)
FLOAT	**Data;
int	NTuples;
int	TupleSize;
FLOAT	*Median;
FLOAT	*StdDev;
{
#ifdef	ANSI_C
	FLOAT	MedianOfY();
	FLOAT	StdDevOfY();
/*debug
	FLOAT	MedianOfY(FLOAT **Data,int NTuples,int TupleSize);
	FLOAT	StdDevOfY(FLOAT **Data,int NTuples,int TupleSize);
debug*/
#else
	FLOAT	MedianOfY();
	FLOAT	StdDevOfY();
#endif
	*StdDev = StdDevOfY(Data,NTuples,TupleSize);	
	*Median = MedianOfY(Data,NTuples,TupleSize);
}


FLOAT
StdDevOfX(Data,NTuples,TupleSize)
FLOAT	**Data;
int	NTuples;
int	TupleSize;
{
	int	i;
	FLOAT	StdDev, x;
	FLOAT	var, sum, mean;
	double sqrt(), fabs();

	if (VtoBoolean(V_SupplyAbscissa)) {
		return((FLOAT)(NTuples*NTuples)/12.0);
	}

	sum = 0.0;
	for (i=0; i<NTuples; i++)
		sum += Data[i][0];
	mean = sum / NTuples;

	var = 0.0;
	for (i=0; i<NTuples; i++) {
		x = fabs((double)(Data[i][0] -  mean));
		var += x * x;
	}
	var /= (NTuples-1);
	StdDev = sqrt((double)var);
	return(StdDev);	
}


FLOAT
StdDevOfY(Data,NTuples,TupleSize)
FLOAT	**Data;
int	NTuples;
int	TupleSize;
{
	int	i,j,D;
	FLOAT	StdDev, x;
	FLOAT	var, sum, mean;
	FLOAT	n;
	double sqrt(), fabs();

	D = VtoBoolean(V_SupplyAbscissa) ? 0 : 1;
	n = NTuples * (TupleSize-D);
	sum = 0.0;
	for (i=0; i<NTuples; i++)
		for (j=D; j<TupleSize; j++) {
			sum += Data[i][j];
		}
	mean = sum / n;

	var = 0.0;
	for (i=0; i<NTuples; i++)
		for (j=D; j<TupleSize; j++) {
			x = fabs((double)(Data[i][j] -  mean));
			var += x * x;
		}
	var /= (n-1);
	StdDev = sqrt((double)var);
	return(StdDev);	
}

#define	BIG	1.0e30
#define	AFAC	1.5
#define	AMP	1.5

FLOAT
MedianOfX(Data,NTuples,TupleSize)
FLOAT	**Data;
int	NTuples;
int	TupleSize;
/* an iterative computation of the median...
 * Code taken from Numerical Recipes..
 */
{
	int	np, nm, i;
	FLOAT	xx,xp,xm,sumx,sum,eps,stemp,dum,ap,am,aa,a;
	FLOAT	Median;
	double	fabs();

	if (VtoBoolean(V_SupplyAbscissa))
		return((FLOAT)(0.5*NTuples));

	/* first estimate	*/
	a = 0.5 * (Data[0][0] + Data[NTuples-1][0]);
	eps = fabs((double)(Data[NTuples-1][0] - Data[0][0]));

	am = -(ap=BIG);
	for (;;) {
		sum = sumx = 0.0;
		np = nm = 0;		/* # point above and below median */
		xm = -(xp = BIG);

		for (i=0; i<NTuples; i++) {
			xx = Data[i][0];
			if (xx != a) {
				if (xx > a) {
					++np;
					if (xx < xp) xp = xx;
				} else if (xx < a) {
					++nm;
					if (xx > xm) xm = xx;
				}
				sum += dum=1.0/(eps+fabs((double)(xx-a)));
				sumx += xx * dum;
			}
		}

		stemp = (sumx / sum) - a;
		if ((np-nm) >= 2) {
			/* guess is too low make another pass */
			am = a;
			aa = stemp < 0.0? xp : xp + stemp*AMP;
			if (aa > ap) aa = 0.5*(a+ap);
			eps = AFAC * fabs((double)(aa-a));
			a = aa;
		} else if ((nm-np) >= 2) {
			/* guess is too high make another pass */
			ap = a;
			aa = stemp > 0.0? xm : xm+stemp*AMP;
			if (aa < am) aa = 0.5*(a+am);
			eps = AFAC*fabs((double)(aa-a));
			a = aa;
		} else {	/* got it	*/
			if (NTuples%2 == 0) {
				Median = 0.5*(np==nm? xp+xm : np>nm? a+xp : xm+a);
			} else {
				Median = np == nm? a : np>nm? xp : xm;
			}
			return(Median);
		}
	}
}

FLOAT
MedianOfY(Data,N,M)
FLOAT	**Data;
int	N;
int	M;
/* an iterative computation of the median...
 * Code taken from Numerical Recipes..
 */
{
	int	n, np, nm, i,j,D;
	FLOAT	xx,xp,xm,sumx,sum,eps,stemp,dum,ap,am,aa,a;
	FLOAT	Median;
	double	fabs();

	D = VtoBoolean(V_SupplyAbscissa) ? 0 : 1;
	n = N * (M-D);

	/* first estimate	*/
	a = 0.5 * (Data[0][D] + Data[N-1][M-1]);
	eps = fabs((double)(Data[N-1][M-1] - Data[0][D]));

	am = -(ap=BIG);
	for (;;) {
		sum = sumx = 0.0;
		np = nm = 0;		/* # point above and below median */
		xm = -(xp = BIG);

		for (i=0; i<N; i++) {
			for (j=D; j<M; j++) {
				xx = Data[i][j];
				if (xx != a) {
					if (xx > a) {
						++np;
						if (xx < xp) xp = xx;
					} else if (xx < a) {
						++nm;
						if (xx > xm) xm = xx;
					}
					sum += dum=1.0/(eps+fabs((double)(xx-a)));
					sumx += xx * dum;
				}
			}
		}

		stemp = (sumx / sum) - a;
		if ((np-nm) >= 2) {
			/* guess is too low make another pass */
			am = a;
			aa = stemp < 0.0? xp : xp + stemp*AMP;
			if (aa > ap) aa = 0.5*(a+ap);
			eps = AFAC * fabs((double)(aa-a));
			a = aa;
		} else if ((nm-np) >= 2) {
			/* guess is too high make another pass */
			ap = a;
			aa = stemp > 0.0? xm : xm+stemp*AMP;
			if (aa < am) aa = 0.5*(a+am);
			eps = AFAC*fabs((double)(aa-a));
			a = aa;
		} else {	/* got it	*/
			if (n%2 == 0) {
				Median = 0.5*(np==nm? xp+xm : np>nm? a+xp : xm+a);
			} else {
				Median = np == nm? a : np>nm? xp : xm;
			}
			return(Median);
		}
	}
}