//                                22.MAY.2006                   filoelastico.cc

#include <curses.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.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,yScale,yScale2;
	double *tt,*y;
	char answ[122],buff[40];
	char **comandi,**ParName,**ParVal;
	FILE *hnd;
	// ............................................................. modo grafico
	StartXWindow(800,600,0,0,"Filo elastico");
	// ...................................................... allocazione memoria
	nPar=3;
	ParName=new (char *) [nPar];
	ParVal=new (char *) [nPar];
	for (i=0;i<nPar;i++) ParVal[i]=new char [24];
	ParName[0]="Lunghezza";
	ParName[1]="Costante di Hooke";
	ParName[2]="Smorzamento";
	// .................................................................. comandi
	nComm=4;
	comandi=new (char *) [nComm];
	comandi[0]="caduta";
	comandi[1]="calcione";
	comandi[2]="ferma";
	comandi[3]="fine";
	// .......................................................... allocate arrays
	tt=new double [4];
	y=new double [4];
	// ................................................................. costanti
	tt[1]=300; //           lunghezza
	tt[2]=20; //    costante di Hooke
	tt[3]=0.005; //       smorzamento
	sprintf(ParVal[0],"%.4lf",tt[1]);
	sprintf(ParVal[1],"%.4lf",tt[2]);
	sprintf(ParVal[2],"%.4lf",tt[3]);
	// ...................................................... 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
	XPermMenu(1,"Filo elastico",ParName,ParVal,nPar,16,comandi,nComm,20,400,&mx,&my,0);
	// ............................................... primo passo di Runge-Kutta
	tt[0]=0;
	RungeStep(tt,y,4,step,0,func);
	XSync(xvd.display,TRUE);
	// ..................................................... ciclo di Runge Kutta
	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,xvd.gc,0,0,800,400);
		// ..................................................... disegna il pendolo
		XSetForeground(xvd.display,xvd.gc,xvd.black);
		XFillArc(xvd.display,xvd.backpix,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,xvd.gc,ix0,iy0,ix,iy);
		XSetForeground(xvd.display,xvd.gc,xvd.red[1]);
		ix-=10;
		iy-=10;
		XFillArc(xvd.display,xvd.backpix,xvd.gc,ix,iy,20,20,0,23040);
		XCopyArea(xvd.display,xvd.backpix,xvd.win,xvd.gc,0,0,800,400,0,0);
		XSync(xvd.display,FALSE);
		// .......................................................... leggi il menu
		ii=XPermMenu(1,"Filo elastico",ParName,ParVal,nPar,16,comandi,nComm,20,
								 400,&mx,&my,1);
		if (ii==0)
		{ // ............................................................... caduta
			tt[0]=0;
			tt[1]=atof(ParVal[0]); //         lunghezza
			tt[2]=atof(ParVal[1]); // costante di Hooke
			tt[3]=atof(ParVal[2]); //       smorzamento
			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;
			tt[1]=atof(ParVal[0]); //           lunghezza
			tt[2]=atof(ParVal[1]); //   costante di Hooke
			tt[3]=atof(ParVal[2]); //         smorzamento
			y[0]=0; //                               x0=0
			y[1]=sqrt(9.8*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
			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]);
			tt[2]=atof(ParVal[1]);
		}
		XSync(xvd.display,TRUE);
	}
	// .................................................... last Runge-Kutta step
	RungeStep(tt,y,4,step,2,func);
	// ............................................................... erase menu
	ii=XPermMenu(1,"Filo elastico",ParName,ParVal,nPar,16,comandi,nComm,20,
							 400,&mx,&my,2);
	// ..........................................................................
	delete [] y;
	delete [] tt;
	CloseXWindow();
}

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

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