// ................................................................. 3dpolar.cc
//   16.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 i,ii,ix,ix0,ixold,ixstart,iy0,iy,iyold,iystart,j,nComm,nPar,nPhi,nR,redraw;
	int *ParFont;
	unsigned long color;
	double alpha,alphadeg,beta,betadeg,DeltaPhi,DeltaR,dist,gamma,gammadeg,
	length,phi,r,radius,x,y;
	const double deg2rad=M_PI/180.0;
	double rnew[3],rold[3];
	double *par;
	double **rotmat,**zz;
	char **comandi,**ParName,**ParVal;
	clock_t musec1,musec2;
	timespec ts1,ts2;
	FILE *hnd;
	
	// ..........................................................................
	color=xvd.blue[0];
	// ............................................................. modo grafico
	StartXWindow(800,600,0,0,"3D Polar",2);
	// .......................................................... numero di punti
	nPhi=50;
	nR=50;
	length=800.0;
	radius=350.0;
	dist=200000.0;
	DeltaR=radius/double(nR);
	DeltaPhi=2.0*M_PI/double(nPhi);
	// ...................................................... allocazione memoria
	rotmat=new double * [3];
	for (i=0;i<3;i++) rotmat[i]=new double[3];
	zz=new double * [nR];
	for (i=0;i<nPhi;i++) zz[i]=new double [nPhi];
	// ............................................. allocazione memoria per menu'
	nPar=7;
	ParName=new char * [nPar];
	ParVal=new char * [nPar];
	for (i=0;i<nPar;i++) ParVal[i]=new char [24];
	ParName[0]="a";
	ParName[1]="b";
	ParName[2]="g";
	ParName[3]="a";
	ParName[4]="b";
	ParName[5]="c";
	ParName[6]="Distanza";
	nComm=1;
	comandi=new char * [nComm];
	comandi[0]="Fine";
	ParFont=new int [nPar];
	for(i=0;i<3;i++) ParFont[i]=5;
	for (i=3;i<7;i++) ParFont[i]=1;
	// .......................................................... allocate arrays
	par=new double [3];
	// ................................................................. costanti
	par[0]=1.0; //      a
	par[1]=0.0; //      b
	par[2]=0.0; //      c
	// ................................................................... angoli
	alphadeg=0.0;
	betadeg=0.0;
	gammadeg=0.0;
	// ....................................................... parametri per menu
	sprintf(ParVal[0],"%.4lf",alphadeg);
	sprintf(ParVal[1],"%.4lf",betadeg);
	sprintf(ParVal[2],"%.4lf",gammadeg);
	sprintf(ParVal[3],"%.4lf",par[0]);
	sprintf(ParVal[4],"%.4lf",par[1]);
	sprintf(ParVal[5],"%.4lf",par[2]);
	sprintf(ParVal[6],"%.0lf",dist);
	// ................................................. origine delle coordinate
	ix0=400;
	iy0=500;
	// .......................................................... disegna il menu
	XWinPermMenu(1,"3D Polar",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 Polar",ParFont,ParName,ParVal,nPar,16,comandi,nComm,1);
		if (ii==0) break;
		gammadeg+=0.1;
		// ........................................................... leggi i dati
		alphadeg=atof(ParVal[0]);
		betadeg=atof(ParVal[1]);
		//gammadeg=atof(ParVal[2]);
		par[0]=atof(ParVal[3]);
		par[1]=atof(ParVal[4]);
		par[2]=atof(ParVal[5]);
		dist=atof(ParVal[6]);
		// .................................................... gradi -> radianti
		alpha=alphadeg*deg2rad;
		beta=betadeg*deg2rad;
		gamma=gammadeg*deg2rad;
		// ................................................. matrice di rotazione
		rotmat[0][0]=cos(gamma)*cos(beta);
		rotmat[0][1]=-sin(gamma)*cos(beta);
		rotmat[0][2]=sin(beta);
		rotmat[1][0]=cos(gamma)*sin(beta)*sin(alpha)+sin(gamma)*cos(alpha);
		rotmat[1][1]=cos(gamma)*cos(alpha)-sin(gamma)*sin(beta)*sin(alpha);
		rotmat[1][2]=-cos(beta)*sin(alpha);
		rotmat[2][0]=sin(gamma)*sin(alpha)-cos(gamma)*sin(beta)*cos(alpha);
		rotmat[2][1]=sin(gamma)*sin(beta)*cos(alpha)+sin(alpha)*cos(gamma);
		rotmat[2][2]=cos(beta)*cos(alpha);
		// ............................................... coordinate da plottare
		for (i=0;i<nR;i++)
		{
			r=double(i)*DeltaR;
			for (j=0;j<nPhi;j++)
			{
				phi=double(j)*DeltaPhi;
				x=r*cos(phi);
				y=r*sin(phi);
				zz[i][j]=func(x,y,par);
			}
		}
		// .................................................. pulisci la finestra
		XSetForeground(xvd.display,xvd.gc,xvd.white);
		XFillRectangle(xvd.display,xvd.backpix[0],xvd.gc,0,0,800,600);
		// ..................................................... disegna gli assi
		XSetForeground(xvd.display,xvd.gc,xvd.blue[2]);
		// ............................................................... asse x
		rold[0]=500.0;rold[1]=0.0;rold[2]=0.0;
		MatVectMult(rotmat,rold,rnew,3);
		ix=ix0+int(rnew[0]+0.5);
		iy=iy0-int(rnew[2]+0.5);
		XDrawLine(xvd.display,xvd.backpix[0],xvd.gc,ix,iy,ix0,iy0);
		// ............................................................... asse y
		rold[0]=0.0;rold[1]=500.0;rold[2]=0.0;
		MatVectMult(rotmat,rold,rnew,3);
		ix=ix0+int(rnew[0]+0.5);
		iy=iy0-int(rnew[2]+0.5);
		XDrawLine(xvd.display,xvd.backpix[0],xvd.gc,ix,iy,ix0,iy0);
		// ............................................................... asse z
		rold[0]=0.0;rold[1]=0.0;rold[2]=500.0;
		MatVectMult(rotmat,rold,rnew,3);
		ix=ix0+int(rnew[0]+0.5);
		iy=iy0-int(rnew[2]+0.5);
		XDrawLine(xvd.display,xvd.backpix[0],xvd.gc,ix,iy,ix0,iy0);
		// ....................................................... linee r=const.
		for (i=1;i<nR;i++)
		{
			r=double(i)*DeltaR;
			phi=0;
			rold[0]=r*cos(phi);
			rold[1]=r*sin(phi);
			rold[2]=zz[i][0];
			MatVectMult(rotmat,rold,rnew,3);
			ixold=ixstart=ix0+int(rnew[0]*dist/(dist+rnew[1])+0.5);
			iyold=iystart=iy0-int(rnew[2]*dist/(dist+rnew[1])+0.5);
			for (j=1;j<=nPhi;j++)
			{
				phi=double(j%nPhi)*DeltaPhi;
				rold[0]=r*cos(phi);
				rold[1]=r*sin(phi);
				rold[2]=zz[i][j%nPhi];
				MatVectMult(rotmat,rold,rnew,3);
				ix=ix0+int(rnew[0]*dist/(dist+rnew[1])+0.5);
				iy=iy0-int(rnew[2]*dist/(dist+rnew[1])+0.5);
				XSetForeground(xvd.display,xvd.gc,color);
				XDrawLine(xvd.display,xvd.backpix[0],xvd.gc,ixold,iyold,ix,iy);
				ixold=ix;
				iyold=iy;
			}
		}
		// ................................................... righe phi = const.
		for (i=0;i<nPhi;i++)
		{
			phi=double(i)*DeltaPhi;
			rold[0]=0;
			rold[1]=0;
			rold[2]=zz[0][0];
			MatVectMult(rotmat,rold,rnew,3);
			ixold=ix0+int(rnew[0]*dist/(dist+rnew[1])+0.5);
			iyold=iy0-int(rnew[2]*dist/(dist+rnew[1])+0.5);
			for (j=1;j<nR;j++)
			{
				r=double(j)*DeltaR;
				rold[0]=r*cos(phi);
				rold[1]=r*sin(phi);
				rold[2]=zz[j][i];
				MatVectMult(rotmat,rold,rnew,3);
				ix=ix0+int(rnew[0]*dist/(dist+rnew[1])+0.5);
				iy=iy0-int(rnew[2]*dist/(dist+rnew[1])+0.5);
				XSetForeground(xvd.display,xvd.gc,color);
				XDrawLine(xvd.display,xvd.backpix[0],xvd.gc,ixold,iyold,ix,iy);
				ixold=ix;
				iyold=iy;
			}
		}
		// ...................................................... 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 Polar",ParFont,ParName,ParVal,nPar,16,comandi,nComm,2);
	// ............................................................. clear memory
	delete [] par;
	delete [] comandi;
	for (i=nPar-1;i>=0;i--) delete ParVal[i];
	delete [] ParVal;
	delete [] ParName;
	for (i=nR-1;i>=0;i--) delete zz[i];
	delete [] zz;
	for (i=2;i>=0;i--) delete [] rotmat[i];
	delete [] rotmat;
	// ....................................................... 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 0.005*(a*x*x+b*y*y);
}