#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define N 4		

const double g 	= 9.81; 	//!< constante de gravité
double       l 	= 1.0;  	//!< longueur du pendule
double	 alpha	= 3.14/4.0;	//!angle d'inclinaison du systeme(ou du vecteur gravite)

double  *a1, *a2, *a3, *a4,*ec,*ep,*em,*s,*S1,*S2,*S3,*S4;


/* conditions initiales de notre probleme */
void init(double *S) 
{
	
	S[0]	=	1.0;//------->	theta		
  	S[1]	=	0.0;//------->	vitesse en theta	
	S[2]  	= 	0.0;//------->	phi		
  	S[3]	=	0.0;//------->  vitesse en phi 		

}

void derive (double *S ,double* dS)
{

		dS[0]=	S[1] ;
		dS[1]=	S[3]*S[3]*(sin(2*S[0])/2)-(g/l)*sin(S[0])*cos(alpha) ; 	//------->acceleration en theta
		dS[2]=	S[3] ;
		dS[3]=	-2*S[1]*S[3]/tan(S[0]) ;		   		//------->acceleration en phi
}


void RK2(double *S, double h)
{
	int i;
	derive(S, a1);
	for(i=0;i<N;i++){

		S1[i]=S[i]+a1[i]*h/2;
		derive(S1, a2);

		S2[i]=S1[i]+a2[i]*h/2;
		derive(S2, a3);
		
		S3[i]=S2[i]+a3[i]*h;
		derive(S3,a4);

		S[i]=S[i]+(h/6)*(a1[i]+2*a2[i]+2*a3[i]+a4[i]);

			}		
}


int main ()
{

	double  h, T,emoy,cpt,m=1.0;
	int n;
	int  i,j,k;
	double ecart;
	// On ouvre un fichier dans lequel on va écrire

	//FILE *fp;
	FILE *PHI;
	//FILE *nrj;
	
	// On détermine le plafond n

	T	=	250;					//Temps de simulation
	h 	= 	1e-5;					//pas de temps
	n 	= 	 T/h;					//nbre d'iterations temporelles
	
	s 	=(double*)malloc((N)*sizeof(double));
	S1	=(double*)malloc((N)*sizeof(double));
	S2	=(double*)malloc((N)*sizeof(double));
	S3	=(double*)malloc((N)*sizeof(double));
	
	a1	=(double*)malloc((N)*sizeof(double));
	a2	=(double*)malloc((N)*sizeof(double));
	a3	=(double*)malloc((N)*sizeof(double));
	a4	=(double*)malloc((N)*sizeof(double));
	
	
	ec	=(double*)malloc((n+1)*sizeof(double));   	// energie cinetique du systeme
	ep	=(double*)malloc((n+1)*sizeof(double));   	// energie potentielle du systeme
	em	=(double*)malloc((n+1)*sizeof(double));   	// energie mecanique du systeme

//conditions initiales
for(k=6;k<7;k++){

		char  filename[80];
		
		//----------------------------------------------------------------------
		//-- creation des noms des fichiers multiples 
		//----------------------------------------------------------------------
		snprintf( filename, sizeof(filename), "projetp%d.dat", k );
		fp  = fopen( filename, "w" );

		snprintf( filename, sizeof(filename), "phip%d.dat", k );
		PHI = fopen( filename, "w" );

		//Valeur des diferentes energies

		snprintf( filename, sizeof(filename), "nrj%d.dat", k );
		nrj = fopen( filename, "w" );

		double pi=2*asin(1);

		em[0] = 9.0;
				/*		
			theta	<-> s[0]
  			vtheta	<-> s[1]
			phi	<-> s[2]
  			vphi 	<-> s[3]    */
		s[0]  = pi/k;
		s[1]  = sqrt((2*em[0]/(m*l*l))-2*(g*(1-cos(s[0])*cos(alpha))/l)-s[3]*s[3]*s[0]*s[0]);
		s[2]  = 0.0;
		s[3]  = sqrt(g/l);
	fprintf(stderr,"numero de k=%d\n",k);
		
	// CALCUL Ep ,Ec et Em

	ec[0]	=	0.5*m*l*l*(s[3]*s[3]*sin(s[0])*sin(s[0])+s[1]*s[1]);
	ep[0]	=	em[0]-ec[0];

	fprintf(stderr,"%0.1lf  %0.1lf  %0.1lf\n",ec[0],ep[0],em[0]);
		
			
	for (i=0; i<n; i++) {	
		
		double phiv	=	sin(s[2]);
 	
		RK2(s,h);
		if(phiv*sin(s[2])<0)	
		{
			if(phiv<0)
				fprintf (PHI, "%0.16lf  %0.16lf  %0.16lf  %0.16lf \n", s[0] , s[1] , s[2]/acos(-1.0),s[3]);
		}

		ec[i]	=	0.5*m*l*l*(s[3]*s[3]*sin(s[0])*sin(s[0])+s[1]*s[1]);       	// calcul de l'energie cinetique a chaque pas de temps
		ep[i]	=	m*g*l*(1-cos(s[0]))*cos(alpha);                  		// calcul de l'energie potentiel a chaque pas de temp
		em[i]	=	ec[i]+ep[i];  							// calul de l'energie mecanique du systeme
	
			if(i%100==0)
	
				fprintf (fp, "%0.16lf  %0.16lf  %0.16lf  %0.16lf \n", s[0] ,s[1], s[2], s[3]);
		

	}
	
	
		
//Calcul de l'energie moyenne de la valeur d'energie mecanique et de son ecart type.
	cpt	= 0.0;
	for (j=0 ; j<n ;j++)
	{
		cpt	=  	cpt+em[j];
		
	}
	emoy=cpt/n;
	
	ecart=0;
	for (j=0 ; j<n ;j++)
	{
		ecart += (em[j]-emoy)*(em[j]-emoy);
	}
	
	ecart=sqrt(ecart/n);

	
	/*Impression dans les fichiers de l'energie mecanique moyenne et l'ecart type de celle ci*/
	
	printf ("energie moyenne %0.16lf\n",emoy);
	printf ("ecart-type      %0.16lf\n", ecart);




	}

	free 	(s );
	free 	(a1);
	free 	(a2);
	free 	(a3);
	free 	(a4);
	free 	(S1);
	free 	(S2);
	free 	(S3);
	free 	(S4);
	free 	(ec);
	free 	(ep);
	free 	(em);
	fclose	(fp);
	fclose 	(PHI);
	fclose 	(nrj);


	return 0;//fin du programme.
}