// .................................................................... buca.cc
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <corso.h>

extern XVideoData xvd;

void func(double *x,double *y,double *dydx);

main()
{
	int DeltaSave,i,j,nEigen,nPoints,iy1,iy2,yEigv;
	unsigned long color [10];
	double DeltaX,DeltaX2,EigvStart,EigvStep,tolerance,xx,xMax,yy;
	double par[4];
	double *eigv,*x;
	double **y;
	
	// ............................................................. modo grafico
	StartXWindow(800,600,0,0,"Buca",1);
	XSelectInput(xvd.display,xvd.win[0],ButtonPressMask|ButtonReleaseMask|
			PointerMotionMask|KeyPressMask|KeyReleaseMask|ExposureMask|
					StructureNotifyMask);
	// ..........................................................................
	color[0]=xvd.black;
	color[1]=xvd.red[0];
	color[2]=xvd.green[0];
	color[3]=xvd.blue[0];
	color[4]=xvd.cyan[1];
	// ..........................................................................
	DeltaSave=4;
	// ..........................................................................
	nEigen=5;
	eigv=new double [nEigen];
	// ..........................................................................
	nPoints=500;
	x=new double [nPoints+8];
	y=new double * [2];
	for (i=0;i<2;i++) y[i]=new double [nPoints+8];
	// ..........................................................................
	EigvStep=0.005;
	tolerance=1.0e-12;
	par[2]=1.0;
	par[3]=2.0;
	xMax=10.0;
	// .................................................... disegna il potenziale
	DeltaX=double(xMax)/double(nPoints);
	DeltaX2=DeltaX*DeltaX;
	XSetForeground(xvd.display,xvd.gc,xvd.black);
	iy1=0;
	for (i=0;i<800;i++)
	{ 
		yy=DeltaX*double((i-399)*(i-399));
		iy2=int(599.0-yy);
		XDrawLine(xvd.display,xvd.backpix[0],xvd.gc,i-1,iy1,i,iy2);
		iy1=iy2;
	}
	XDrawLine(xvd.display,xvd.backpix[0],xvd.gc,399,0,399,599);
	XRedraw();
	// ..........................................................................
	EigvStart=0.0;
	for (i=0;i<nEigen;i++)
	{
		eigv[i]=SymmWell(par,x,y,nPoints,DeltaSave,xMax,i,EigvStart,EigvStep,
										 tolerance,func);
		printf("%4d%15.5le\n",i,eigv[i]);
		EigvStart=eigv[i]+EigvStep;
		yEigv=599-int(2.0*eigv[i]/DeltaX);
		PlotPsi(y[0],nPoints,(i+1)%2,400.0,yEigv,color[i]);
		iy1=599-int(2.0*eigv[i]/DeltaX);
		XDrawLine(xvd.display,xvd.backpix[0],xvd.gc,0,iy1,799,iy1);
		XRedraw();
	}
	// ..........................................................................
	XWait(0);
	CloseXWindow();
	// ...........................................................................
	for (i=1;i>=0;i--) delete [] y[i];
	delete [] y;
	delete [] x;
	delete [] eigv;
}

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

void func(double *x,double *y,double *dydx)
{
	double A,B,E;
	
	E=x[1];
	A=x[2];
	B=x[3];
	// ................................................................ velocita'
	dydx[0]=y[1];
	// ............................................................ accelerazione
	dydx[1]=(A*x[0]*x[0]-B*E)*y[0];
}