//																02.APR.2009									 xpendolomolla.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 *t,double *y,double *dydt);

main()
{
	int DeltaSave,i,ii,ix,ix0,iy,iy0,iyf,iyf0,mx,my,nComm,nPar;
	double alpha,dt,step,tMax,tMin,tWait,yScale,yScale2;
	double *tt,*y;
	char answ[122],buff[40],title[40];
	char **comandi,**ParName,**ParVal;
	const char *ParFixName[6]=
	{
		"Lunghezza",//               0
		"Massa",//                   1
		"g",//                       2
		"Costante di Hooke",//       3
		"Smorzamento",//             4
		"Intervallo/ms"//            5
	};
	const char *ComFixName[4]=
	{
		"Caduta",//                  0
		"Calcione",//                1
		"Ferma",//                   2
		"Fine"//                     3
	};
	clock_t musec1,musec2;
	timespec ts1,ts2;
	FILE *hnd;
	// ............................................................. modo grafico
	strcpy(title,"Pendolo-molla");
	StartXWindow(800,600,0,0,title,2);
	XSelectInput(xvd.display,xvd.win[0],ButtonPressMask|ButtonReleaseMask|
			PointerMotionMask|KeyPressMask|KeyReleaseMask|ExposureMask|
					StructureNotifyMask);
	// ...................................................... allocazione memoria
	nPar=6;
	ParName=new char * [nPar];
	ParVal=new char * [nPar];
	for (i=0;i<nPar;i++) {ParName[i]=new char [24];ParVal[i]=new char [24];}
	for (i=0;i<nPar;i++) strcpy(ParName[i],ParFixName[i]);
	// .................................................................. comandi
	nComm=4;
	comandi=new char * [nComm];
	for (i=0;i<nComm;i++) {comandi[i]=new char [16];strcpy(comandi[i],ComFixName[i]);}
	// .......................................................... allocate arrays
	tt=new double [6];
	y=new double [4];
	// ............................................................... intervallo
	tWait=5.0;
	// ................................................................. costanti
	tt[1]=300.0; //           lunghezza
	tt[2]=1.0; //             massa
	tt[3]=9.81; //            g
	tt[4]=20.0; //           costante di Hooke
	tt[5]=0.005; //           smorzamento
	sprintf(ParVal[0],"%.4lf",tt[1]);
	sprintf(ParVal[1],"%.4lf",tt[2]);
	sprintf(ParVal[2],"%.4lf",tt[3]);
	sprintf(ParVal[3],"%.4lf",tt[4]);
	sprintf(ParVal[4],"%.4lf",tt[5]);
	sprintf(ParVal[5],"%.4lf",tWait);
	// ...................................................... condizioni iniziali
	y[0]=tt[1]; //               x0=lunghezza
	y[1]=0.0; //         velocita' x iniziale
	y[2]=0.0; //                         y0=0
	y[3]=0.0; //         velocita' y iniziale
	// ...................................................... intervallo di tempo
	step=0.02;
	// ............................................................ campionamento
	DeltaSave=1;
	// ........................................................... punti iniziali
	ix0=400;
	iy0=20;
	// .......................................................... disegna il menu
	XWinPermMenu(0,title,ParName,ParVal,nPar,16,comandi,nComm,0);
	// ............................................... primo passo di Runge-Kutta
	tt[0]=0;
	RungeStep(tt,y,4,step,0,func);
	XSync(xvd.display,TRUE);
	// ................................................ 0 secondi nell'intervallo
	ts1.tv_sec=0;
	// ..................................................... ciclo di Runge Kutta
	i=0;
	musec1=clock();
	for (;;)
	{
		RungeStep(tt,y,4,step,1,func);
		if (i%DeltaSave) continue;
		// ................................ seleziona e pulisci il contesto grafico
		XSetForeground(xvd.display,xvd.gc,xvd.white);
		XFillRectangle(xvd.display,xvd.backpix[0],xvd.gc,0,0,xvd.width[0],
									 xvd.height[0]);
		// ............................................... disegna la circonferenza
		ii=int(tt[1]+0.5);
		XSetForeground(xvd.display,xvd.gc,xvd.yellow[0]);
		XDrawArc(xvd.display,xvd.backpix[0],xvd.gc,ix0-ii,iy0-ii,2*ii,2*ii,0,23040);
		// ..................................................... disegna il pendolo
		XSetForeground(xvd.display,xvd.gc,xvd.black);
		XFillArc(xvd.display,xvd.backpix[0],xvd.gc,ix0-2,iy0-2,4,4,0,23040);
		ix=ix0+int(y[0]+0.5);
		iy=iy0+int(y[2]+0.5);
		XDrawLine(xvd.display,xvd.backpix[0],xvd.gc,ix0,iy0,ix,iy);
		XSetForeground(xvd.display,xvd.gc,xvd.red[1]);
		ix-=10;
		iy-=10;
		XFillArc(xvd.display,xvd.backpix[0],xvd.gc,ix,iy,20,20,0,23040);
		XRedraw();
		musec2=clock();
		ts1.tv_nsec=int(tWait*1.0e6)-(musec2-musec1)*1000;
		if (ts1.tv_nsec>0) nanosleep(&ts1,&ts2);
		musec1=clock();
		XSync(xvd.display,FALSE);
		// .......................................................... leggi il menu
		ii=XWinPermMenu(0,title,ParName,ParVal,nPar,16,
										comandi,nComm,1);
		if (ii==0)
		{ // ............................................................... caduta
			tt[0]=0;
			y[0]=tt[1]; //                 x0=lunghezza
			y[1]=0.0; //           velocita' x iniziale
			y[2]=0.0; //                           y0=0
			y[3]=0.0; //           velocita' y iniziale
		}
		else if (ii==1)
		{ // ............................................................. calcione
			tt[0]=0;
			y[0]=0; //                               x0=0
			y[1]=sqrt(tt[3]*tt[1]); // velocita' x iniziale
			y[2]=tt[1]; //                   y0=lunghezza
			y[3]=0.0; //             velocita' y iniziale
		}
		else if (ii==2)
		{ // ..................................................... ferma il pendolo
			tt[0]=0;
			tt[1]=atof(ParVal[0]); //           lunghezza
			tt[2]=atof(ParVal[1]); //   costante di Hooke
			tt[3]=atof(ParVal[2]); //         smorzamento
			tt[4]=atof(ParVal[4]); //                peso
			y[0]=0; //                               x0=0
			y[1]=0;; //              velocita' x iniziale
			y[2]=tt[1]; //                   y0=lunghezza
			y[3]=0.0; //             velocita' y iniziale
		}
		else if (ii==3) break;
		else if (ii==nComm)
		{ // ........................................ riparti con i nuovi parametri
			tt[1]=atof(ParVal[0]);//         lunghezza
			tt[2]=atof(ParVal[1]);//         massa
			tt[3]=atof(ParVal[2]);//         g
			tt[4]=atof(ParVal[3]);//         Hooke
			tt[5]=atof(ParVal[4]);//         smorzamento
			tWait=atof(ParVal[5]);//         intervallo
			sprintf(ParVal[0],"%.4lf",tt[1]);
			sprintf(ParVal[1],"%.4lf",tt[2]);
			sprintf(ParVal[2],"%.4lf",tt[3]);
			sprintf(ParVal[3],"%.4lf",tt[4]);
			sprintf(ParVal[4],"%.4lf",tt[5]);
			sprintf(ParVal[5],"%.4lf",tWait);
		}
		XSync(xvd.display,TRUE);
		i++;
	}
	// .................................................... last Runge-Kutta step
	RungeStep(tt,y,4,step,2,func);
	// ............................................................... erase menu
	XWinPermMenu(1,title,ParName,ParVal,nPar,16,comandi,nComm,2);
	// ..........................................................................
	delete [] y;
	delete [] tt;
	CloseXWindow();
}

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

void func(double *t,double *y,double *dydt)
{
	double delta,hooke,g,k,length,m,restlength;
	
	restlength=t[1];
	m=t[2];
	g=t[3];
	hooke=t[4];
	k=t[5];
	// .........................................................................
	length=sqrt(y[0]*y[0]+y[2]*y[2]);
	delta=length-restlength;
	// ................................................................ velocita'
	dydt[0]=y[1];
	dydt[2]=y[3];
	// .......................................................... accelerazione x
	dydt[1]=-(hooke*delta*y[0]/length+k*y[1])/m;
	// .......................................................... accelerazione y
	dydt[3]=g-(hooke*delta*y[2]/length+k*y[3])/m;
}