/*
 * %G%	%W%
 */
#include <stdio.h>
#include "xplot.h"
#include <sys/errno.h>


extern int errno;
/*extern  int cmadr[256];*/
/*extern int docolor;*/
/*extern int nulwindow;*/

int 	dontplot = 0;
int 	dodraw = 1.0;
int 	doaxis = 1.0;

/*
 * Draw plot in a window
 *
 */

XTextItem	c[3];			/* Labels for the axis */

xplot(xp)
	struct xpinfo *xp;
{
	register float *pvec;
	register XPoint *xptr;
	register int i;
	register int j;
	register int z;
	register float ival;
	register float jval;
	register float half_window_x;
	register float half_window_y;
	register float half_window_z;
	float	scale_mat[4][4];		/* Scale matrix */
	float	trans_mat[4][4];		/* trans_mat */
	float	back_trans[4][4];		/* restoring trans_mat */
	float 	taxis[3][2][4];
	float	int_mat[4][4];			/* intermediate matrix */
	float	axis[3][2][4];			/* matrix of axis */
	float	*Plot[MAXGRAPHS];		/* actual matrix */
	float	*tplot[MAXGRAPHS];		/* temp matrix */
	XPoint	*thead;
	int 	iinc;
	int		jinc;
	int		origin_x;
	int		origin_y;
	float	xnorm;
	float	ynorm;
	float	znorm;
	int		pix1,pix2;


	if(dontplot) return;

	/*if(docolor < 0 ) {*/
		pix1 =  xp->forePixel;
		pix2 =  xp->forePixel;
#ifdef notdef
	}else {
		pix1 = cmadr[96];
		pix2 = cmadr[48];
	}
#endif
	z = 0;
	half_window_x = xp->Window_X * 0.5;
	half_window_y = xp->Window_Y * 0.5;
	half_window_z = xp->Window_Z * 0.5;

	xnorm = -xp->min_x;
	ynorm = -xp->min_y;
	znorm = -xp->min_z;
	if (!xp->wmap.first) {
 		identity(xp->wmap.world_mat);

		my_fill_mat(scale_mat,0.0,4,4);	
		scale_mat[0][0] =  .5;
		xp->wmap.Sx = (xp->Window_X/(xp->max_x - xp->min_x))*(1.0/xp->X_scale);	
 		scale_mat[1][1] =  .5;
		xp->wmap.Sy = (xp->Window_Y/(xp->max_y - xp->min_y))*(1.0/xp->Y_scale);	
		scale_mat[2][2] =  .5;
		xp->wmap.Sz = (xp->Window_Z/(xp->max_z - xp->min_z))*(1.0/xp->Z_scale);	
		scale_mat[3][3] = 1.0;

		/*
		 * Scale world coords 
		 */
		set_wcm(scale_mat,xp->wmap.world_mat);		/* update the world */

		xp->wmap.first = 1;

		for(i=0; i<3; i++ ) {
			c[i].nchars = 1;
			c[i].delta = 1;
			c[i].font = None;
		}
	}

	/*
	 * Make our axis matrix
	 */
	fill_axis(taxis, (xp->min_x+xnorm)*xp->wmap.Sx,
					 (xp->max_x+xnorm)*xp->wmap.Sx,
			         (xp->min_y+ynorm)*xp->wmap.Sy,
					 (xp->max_y+ynorm)*xp->wmap.Sy,
			  		 (xp->min_z+znorm)*xp->wmap.Sz,
					 (xp->max_z+znorm)*xp->wmap.Sz);
	/*
	 * Set wcm rotations.
	 */
	if (xp->Xtheta) 
		rotate_in_x(RAD(xp->Xtheta),xp->wmap.world_mat);	
	if (xp->Ytheta) 
		rotate_in_y(RAD(xp->Ytheta),xp->wmap.world_mat);		
	if (xp->Ztheta) 
		rotate_in_z(RAD(xp->Ztheta),xp->wmap.world_mat);		
	/*
	 * Translate so origin is in center of all plots, so we can 
	 * rotate about the center of the image.
	 */
	identity(trans_mat);
	trans_mat[3][0] = (-(half_window_x) * 1.0/xp->X_scale);
	trans_mat[3][1] = (-(half_window_y) * 1.0/xp->Y_scale);
	trans_mat[3][2] = (-(half_window_z) * 1.0/xp->Z_scale); 
	
	identity(back_trans);
	back_trans[3][0] = ((half_window_x) * 1.0/xp->X_scale); 
	back_trans[3][1] = ((half_window_y) * 1.0/xp->Y_scale); 
	back_trans[3][2] = ((half_window_z) * 1.0/xp->Z_scale); 
	
	/*
	 * To get our cartesian display, we need to mult our 
	 * picture by the world_mat, and then plot X and Y
	 */
	ival = 0;
	
	iinc = ((xp->Window_X * 1.0/(xp->X_scale))/(xp->numgraphs-1));
	jinc = ((xp->Window_Y * 1.0/(xp->Y_scale))/(xp->numz-1));
 
	my_mat_mult(trans_mat, xp->wmap.world_mat, int_mat, 4, 4);


	if(dodraw > 0) {
		for (i=0; i<xp->numgraphs; i++) {
			if(xp->doxyz) {
				Plot[i]  = (float *) malloc(3 * xp->numz*4*sizeof(float));
				tplot[i] = (float *) malloc(3 * xp->numz*4*sizeof(float));
			}else {
				Plot[i]  = (float *) malloc(xp->numz*4*sizeof(float));
				tplot[i] = (float *) malloc(xp->numz*4*sizeof(float));
			}
			pvec = tplot[i];
			jval = 0;
			for (j=0; j<xp->numz; j++) {
				if(xp->doxyz) { 		/* normalize */
					pvec[DX] = (xp->data[z++]+xnorm) * xp->wmap.Sx;
					pvec[DY] = (xp->data[z++]+ynorm) * xp->wmap.Sy;
				}else {
					pvec[DX] = ival;
					pvec[DY] = jval;
				}
				pvec[DZ] = (xp->data[z++]+znorm) * xp->wmap.Sz;
				pvec[DS] = 1.0;
				jval +=  jinc;
				pvec += 4;
			}
			my_mat_mult(tplot[i], int_mat, Plot[i],xp->numz,4);
			my_mat_mult(Plot[i], back_trans, tplot[i],xp->numz,4);
			ival += iinc;
		}
	}


	/*
	 * Put the axis at origin
	 */
	if(doaxis > 0) {
		for (i=0; i<3; i++) {
			my_mat_mult(taxis[i], int_mat, axis[i], 2, 4);
			my_mat_mult(axis[i], back_trans, taxis[i], 2, 4);
		}
	}

	if(dodraw || doaxis) {
		if(xp->doxyz) {
			if((thead=(XPoint *)malloc(3*xp->numgraphs*xp->numz * 
									xp->numz*sizeof(XPoint))) == NULL) {
				perror("doxyz :");
				exit(1);
			}
		}else {
			if((thead=(XPoint *)malloc(xp->numgraphs * xp->numz *
									xp->numz*sizeof(XPoint))) == NULL) {
				perror("nodoxyz");
				exit(1);
			}
		}
	}
	if(dodraw > 0) {
		xptr = thead;
		for(i = 0; i < xp->numgraphs; i++) {
			pvec = tplot[i];
			for (j = 0; j < xp->numz; j++) {
#ifdef shear
				xptr->x = pvec[DX] / pvec[DS] + 
							((xp->Window_X - (xp->Window_X * 
													1.0/xp->X_scale))/2);
				xptr->y = pvec[DY] / pvec[DS] + 
							((xp->Window_Y - (xp->Window_Y * 
													1.0/xp->Y_scale))/2);
#else
				xptr->x = pvec[DX] + ((xp->Window_X - (xp->Window_X * 
													1.0/xp->X_scale))/2);
				xptr->y = pvec[DY] + ((xp->Window_Y - (xp->Window_Y * 
													1.0/xp->Y_scale))/2);
#endif
				xptr++;
				pvec += 4;
			}
			XDrawLines(xp->dpy,xp->plot,xp->gc,thead,
				xp->numz, CoordModeOrigin);
			xptr = thead;
		}
/*		XDrawLines(xp->dpy,xp->plot,xp->gc,thead,xp->numgraphs*xp->numz,*/
/*						CoordModeOrigin);*/
		xptr = thead;
		for (i = 0; i < xp->numz; i++) {
			for (j = 0; j < xp->numgraphs; j++) {
				pvec = tplot[j] + i * 4;
#ifdef shear
				xptr->x = pvec[DX] / pvec[DS] + 
								   ((xp->Window_X - (xp->Window_X * 
													1.0/xp->X_scale))/2);
				xptr->y = pvec[DY] / pvec[DS] + 
								   ((xp->Window_Y - (xp->Window_Y * 
													1.0/xp->Y_scale))/2);
#else
				xptr->x = pvec[DX] + ((xp->Window_X - (xp->Window_X * 
													1.0/xp->X_scale))/2);
				xptr->y = pvec[DY] + ((xp->Window_Y - (xp->Window_Y * 
													1.0/xp->Y_scale))/2);
#endif
				xptr++;
			}
			XDrawLines(xp->dpy,xp->plot,xp->gc,thead,
				xp->numgraphs, CoordModeOrigin);
			xptr = thead;
		}
/*		XDrawLines(xp->dpy,xp->plot,xp->gc,thead,xp->numgraphs*xp->numz,*/
/*						CoordModeOrigin);*/
	}

	/*
	 * Lastly draw the axis
	 */
	if(doaxis > 0) {
		c[0].chars = "*";
		c[1].chars = "X";
		c[2].chars = "Y";
		c[3].chars = "Z";
		xptr = thead;
		pvec = taxis[0][0];
		origin_x =  pvec[DX] + ((xp->Window_X - (xp->Window_X * 
													1.0/xp->X_scale))/2); 
		origin_y =  pvec[DY] + ((xp->Window_Y - (xp->Window_Y * 
													1.0/xp->Y_scale))/2); 

		XDrawText(xp->dpy,xp->plot,xp->gc,origin_x,origin_y,&c[0],1);
		for(i=0; i<3; i++) {
			for (j=0; j<2; j++) {
#ifdef shear
				xptr->x = (int)(pvec[DX]/pvec[DS]) +
						 + ((xp->Window_X - (xp->Window_X * 
													1.0/xp->X_scale))/2);	
				xptr->y = (int)(pvec[DY]/pvec[DS]) +
						 + ((xp->Window_Y - (xp->Window_Y * 
													1.0/xp->Y_scale))/2);	
				
#else
				xptr->x = pvec[DX] + ((xp->Window_X - (xp->Window_X * 
														1.0/xp->X_scale))/2);
				xptr->y = pvec[DY] + ((xp->Window_Y - (xp->Window_Y * 
														1.0/xp->Y_scale))/2);
#endif
				if (j == 1)
					XDrawText(xp->dpy,xp->plot,xp->gc,xptr->x,
							xptr->y,&c[i+1],1);
				xptr++;
				pvec += 4;
			}
		}
		XDrawLines(xp->dpy,xp->plot,xp->gc,thead,3*2,CoordModeOrigin);
	}

	if(dodraw || doaxis) {
		free(thead);
	}
	for (i=0; i<xp->numgraphs; i++) {
		free(Plot[i]);
		free(tplot[i]);
	}
}