// .............................................................. 3daxonplot.cc
//   26.JUL.2007

#include <curses.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <corso.h>

extern XVideoData xvd;

double func(double x,double y,double *par);

main()
{
	int di,dj,i,ii,im1,ip1,ix,ix0,ixold,ixstart,iy0,iy,iyold,iystart,j,
	jm1,jp1,l,k,kk,nComm,nPar,nX,nY,nTot,skip;
	int *ind,*ParFont;
	unsigned long color;
	double alpha,alphadeg,beta,betadeg,DeltaX,DeltaY,depth,
	length,x,xLen,xScale,y,yLen,yScale,z,zLen,zScale;
	double rnew[3],rold[3];
	const double deg2rad=M_PI/180.0;
	double *par,*zz;
	double **Point3d;
	char **comandi,**ParName,**ParVal;
	clock_t musec1,musec2;
	timespec ts1,ts2;
	XPoint ptswap;
	XPoint pt2d[5];
	
	// ............................................................. modo grafico
	StartXWindow(800,600,0,0,"3D Axon Plot",2);
	color=xvd.blue[0];
	// .......................................................... numero di punti
	nX=50;
	nY=50;
	nTot=nX*nY;
	xLen=400.0;
	yLen=400.0;
	zLen=300.0;
	DeltaX=2.0*xLen/double(nX);
	DeltaY=2.0*yLen/double(nY);
	// ...................................................... allocazione memoria
	ind=new int [nTot];
	zz=new double [nTot];
	Point3d=new double * [nTot];
	for (i=0;i<nTot;i++) Point3d[i]=new double [3];
	// ............................................. allocazione memoria per menu'
	nPar=10;
	ParName=new char * [nPar];
	ParVal=new char * [nPar];
	for (i=0;i<nPar;i++) ParVal[i]=new char [24];
	ParName[0]="Origine X";
	ParName[1]="Origine Y";
	ParName[2]="a";
	ParName[3]="b";
	ParName[4]="a";
	ParName[5]="b";
	ParName[6]="c";
	ParName[7]="Scala X";
	ParName[8]="Scala Y";
	ParName[9]="Scala Z";
	nComm=1;
	comandi=new char * [nComm];
	comandi[0]="Fine";
	ParFont=new int [nPar];
	for (i=0;i<nPar;i++) ParFont[i]=1;
	ParFont[2]=ParFont[3]=5;
	// .......................................................... allocate arrays
	par=new double [3];
	// ................................................................. costanti
	par[0]=1.0; //      a
	par[1]=1.0; //      b
	par[2]=1.0; //      c
	// ................................................. origine delle coordinate
	ix0=100;
	iy0=500;
	xScale=1.0;
	yScale=1.0;
	zScale=1.0;
	// ................................................................... angoli
	alphadeg=110.0;
	betadeg=70.0;
	// ....................................................... parametri per menu
	sprintf(ParVal[0],"%.d",ix0);
	sprintf(ParVal[1],"%.d",xvd.height[0]-iy0);
	sprintf(ParVal[2],"%.4lf",alphadeg);
	sprintf(ParVal[3],"%.4lf",betadeg);
	sprintf(ParVal[4],"%.4lf",par[0]);
	sprintf(ParVal[5],"%.4lf",par[1]);
	sprintf(ParVal[6],"%.4lf",par[2]);
	sprintf(ParVal[7],"%.4lf",xScale);
	sprintf(ParVal[8],"%.4lf",yScale);
	sprintf(ParVal[9],"%.4lf",zScale);
	// .......................................................... disegna il menu
	XWinPermMenu(1,"3D Axon Plot",ParFont,ParName,ParVal,nPar,16,comandi,nComm,0);
	// .............................................................. sincronizza
	XSync(xvd.display,TRUE);
	// .................................................................... ciclo
	ts1.tv_sec=0;
	for (;;)
	{
		XSync(xvd.display,FALSE);
		musec1=clock();
		// .......................................................... leggi il menu
		ii=XWinPermMenu(1,"3D Axon Plot",ParFont,ParName,ParVal,nPar,16,comandi,nComm,1);
		if (ii==0) break;
		//if (ii==nComm)
		if (TRUE)
		{ // ......................................................... leggi i dati
			ix0=atoi(ParVal[0]);
			iy0=xvd.height[0]-atoi(ParVal[1]);
			alphadeg=atof(ParVal[2]);
			betadeg=atof(ParVal[3]);
			par[0]=atof(ParVal[4]);
			par[1]=atof(ParVal[5]);
			par[2]=atof(ParVal[6]);
			xScale=atof(ParVal[7]);
			yScale=atof(ParVal[8]);
			zScale=atof(ParVal[9]);
			// .................................................... gradi -> radianti
			alpha=alphadeg*deg2rad;
			beta=betadeg*deg2rad;
			// ........................................... coordinate tridimensionali
			kk=0;
			for (j=0;j<nY;j++) for (i=0;i<nX;i++)
			{
				x=double(i)*DeltaX*xScale;
				y=double(j)*DeltaY*yScale;
				zz[kk]=func(x,y,par)*zScale;
				kk++;
			}
			// .................................................. pulisci la finestra
			XSetForeground(xvd.display,xvd.gc,xvd.white);
			XFillRectangle(xvd.display,xvd.backpix[0],xvd.gc,0,0,800,600);
			XSetForeground(xvd.display,xvd.gc,xvd.black);
			// .............................................................. disegno
			for (k=nTot-1;k>=0;k--)
			{
				j=k/nX;
				if (j==0) continue;
				i=k%nX;
				if (i==0) continue;
				kk=0;
				for (di=0;di>=-1;di--) for (dj=0;dj>=-1;dj--)
				{
					x=double(i+di)*DeltaX;
					y=double(j+dj)*DeltaY;
					ii=(j+dj)*nX+i+di;
					pt2d[kk].x=ix0+int(x*sin(alpha)+y*sin(beta)+0.5);
					pt2d[kk].y=iy0-int(y*cos(beta)+x*cos(alpha)+zz[ii]+0.5);
					kk++;
				}
				ptswap=pt2d[2];pt2d[2]=pt2d[3];pt2d[3]=ptswap;
				pt2d[4]=pt2d[0];
				// ................................................ pulisci il poligono
				XSetForeground(xvd.display,xvd.gc,xvd.white);
				XFillPolygon(xvd.display,xvd.backpix[0],xvd.gc,pt2d,4,Convex,
										 CoordModeOrigin);
				// ................................................ disegna il reticolo
				XSetForeground(xvd.display,xvd.gc,xvd.black);
				XDrawLines(xvd.display,xvd.backpix[0],xvd.gc,pt2d,5,
									 CoordModeOrigin);
			}
			// ..................................................... disegna gli assi
			XSetForeground(xvd.display,xvd.gc,xvd.blue[2]);
			// ............................................................... asse x
			ix=ix0+int(xLen*sin(alpha)+0.5);
			iy=iy0-int(xLen*cos(alpha)+0.5);
			XDrawLine(xvd.display,xvd.backpix[0],xvd.gc,ix,iy,ix0,iy0);
			// ............................................................... asse y
			ix=ix0+int(yLen*sin(beta)+0.5);
			iy=iy0-int(yLen*cos(beta)+0.5);
			XDrawLine(xvd.display,xvd.backpix[0],xvd.gc,ix,iy,ix0,iy0);
			// ............................................................... asse z
			XDrawLine(xvd.display,xvd.backpix[0],xvd.gc,ix0,iy0-int(zLen+0.5),ix0,iy0);
		}
		// ...................................................... copia sul monitor
		XCopyArea(xvd.display,xvd.backpix[0],xvd.win[0],xvd.gc,0,0,800,600,0,0);
		// ................................................................ aspetta
		musec2=clock();
		ts1.tv_nsec=10000000-(musec2-musec1)*100000; 
		if (ts1.tv_nsec>0) nanosleep(&ts1,&ts2);
	}
	// ............................................................... erase menu
	XWinPermMenu(1,"3D Axon Plot",ParFont,ParName,ParVal,nPar,16,comandi,nComm,2);
	// ............................................................. clear memory
	delete [] ParFont;
	delete [] par;
	delete [] comandi;
	for (i=nPar-1;i>=0;i--) delete ParVal[i];
	delete [] ParVal;
	delete [] ParName;
	for (i=nTot-1;i>=0;i--) delete [] Point3d[i];
	delete [] Point3d;
	delete [] zz;
	delete [] ind;
	// ....................................................... chiudi le finestre
	CloseXWindow();
}

// ----------------------------------------------------------------------------

double func(double x,double y,double *par)
{
	double a,b,c;
	
	a=par[0];
	b=par[1];
	c=par[2];
	
	//return 0.005*(a*(x*x+y*y)+b*x*y*y-c*x*x*x);
	return 100.0*a*b*b*c*c/((x*x+b*b)*(y*y+c*c))+
			100.0*a*b*b*c*c/(((x-100.0)*(x-100.)+b*b)*((y-150.0)*(y-150.0)+c*c));
}