#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

/***********************declaratino des variables globales***********************************************/
/********************************************************************************************************/

FILE *pf1;
FILE *pf2;

//Définition pointeurs

double *C,*C1,*C2,*DTC;	
double *pC,*pDTC,*pC1,*pC2;
	
//parametres du problème

double mu,gama;

//parametres d'espace

int ne,n;					//longueurs discrétisées
double dl;
int longueur;

//parametres de temps

int nt;						//temps_max discrétisé
int temps_max,temps_init;
double dt;

//parametres de phi

int phi_1,phi_2;
double alpha,beta;
double borne_phi;

//parametres de phi discrétisés

int dphi_1,dphi_2;
int borne_integ;

//partie dialogue

char choix1[20],choix2[20];
int choix0;

/***********************routine d'initialisation des paramètres******************************************/
/********************************************************************************************************/

void initialisation(void)

{

mu=1.;
gama=3.;

temps_max=5;
longueur=15;

dl=0.1;
dt=0.001;

phi_1=1;
phi_2=2;
alpha=2.;
beta=-1.;

borne_phi=2.;				//borne d'integration de phi non discretisée
borne_integ=ceil(borne_phi/dl);		//borne d'integration de phi discretisée

ne=ceil((double)longueur/dl);
n=ne+(2*borne_integ);
nt=ceil((double)temps_max/dt);

printf("%d  %d %d\n",ne,n,nt);
printf("%d  %d\n",dphi_1,dphi_2);
printf("%d\n",borne_integ);

}

/***********************routine d'allocation des vecteurs************************************************/
/********************************************************************************************************/

void allocations(void)
{

C= (double*) malloc((n)*sizeof(double)); 
C=C+borne_integ;
C1=(double*) malloc((n)*sizeof(double));
C1=C1+borne_integ;
C2=(double*) malloc((n)*sizeof(double));
C2=C2+borne_integ;
DTC=(double*) malloc((n)*sizeof(double));
DTC=DTC+borne_integ;

}

/***********************fonction poids phi***************************************************************/
/********************************************************************************************************/

double phi(double s)
{

double res;

if(s>=-phi_2 && s<=-phi_1)
	{
	res=beta;	
	}
else if(s>-phi_1 && s<phi_1)
	{
	res=alpha;	
	}
else if(s>=phi_1 && s<=phi_2)
	{
	res=beta;
	}
else
	{
	res=0.;
	}

return res;

}

/*double phi(double s)
{

double res;

if(s==0)
	{
	res=-2.;	
	}
else if(s==1)
	{
	res=1.;	
	}
else if(s==-1)
	{
	res=1.;
	}
else
	{
	res=0.;
	}

return res;

}*/

/*double phi(double s)
{

return (-11.719)*(s*s)+(1.87507);

}*/

/***********************fonction d'integration de c(x+y)*phi(y)dy****************************************/
/********************************************************************************************************/

double integ(double infinteg,double supinteg)

{
	double res;
	int a,b;
	int j,jmax;
	
	res=0.;
	jmax=(supinteg-infinteg);

	
	for(j=0;j<jmax;j++)
	{
	a=-borne_integ+j;
	b=a+1;
	//res=res+(((*(pC+a))*phi(a)+(*(pC+b))*phi(b))/2.);
	res=res+(*(pC1+a))*phi(a*dl);
	}
	
	res=res*dl;
	
	return res;
}

/***********************calcul des conditions aux limites periodiques de c*******************************/
/********************************************************************************************************/

void conditionslimites(void)
{

int l;
pC1=C1;
for(l=0;l<borne_integ;l++)     *(pC1+l-borne_integ)=*(pC1+ne+l-borne_integ);
for(l=(n-borne_integ);l<n;l++) *(pC1+l-borne_integ)=*(pC1+l+borne_integ-n);

}


/***********************calcul de la derivee de c par rapport a t****************************************/
/********************************************************************************************************/

double dt_c(void)

{

int i;

pC1=C1;
pDTC=DTC;

for(i=0;i<ne;i++)
{

(*pDTC)=-mu*(*pC1)+gama*integ(i-borne_integ,i+borne_integ)-(*pC1)*(*pC1)*(*pC1);

pC1=pC1+1;
pDTC=pDTC+1;	

}

pC1=C1;
pDTC=DTC;

}


/***********************runge kutta*********************************************************************/
/********************************************************************************************************/

double rungekutta(void)
{

int i;

pC=C;
pC1=C1;
pC2=C2;
pDTC=DTC;


//c_sol -> C
//c_int -> C1
//c_old -> C2

for(i=0;i<ne;i++) //sauvegarde de c_old
{

(*pC2)=(*pC);
(*pC1)=(*pC);

pC=pC+1;
pC1=pC1+1;
pC2=pC2+1;
pDTC=pDTC+1;

}

pC1=C1;
pC2=C2;
pC=C;
pDTC=DTC;

conditionslimites();
dt_c();

	/**************etape1***************/

for(i=0;i<ne;i++)
{

(*pC)=(*pC2)+dt*((*pDTC)/6.);
(*pC1)=(*pC2)+dt*(*pDTC)*0.5;

pC=pC+1;
pC1=pC1+1;
pC2=pC2+1;
pDTC=pDTC+1;
}

pC1=C1;
pC2=C2;
pC=C;
pDTC=DTC;

conditionslimites();
dt_c();

	/**************etape2***************/

for(i=0;i<ne;i++)
{

(*pC)=(*pC)+dt*((*pDTC)/3.);
(*pC1)=(*pC2)+dt*(*pDTC)*0.5;

pC=pC+1;
pC1=pC1+1;
pC2=pC2+1;
pDTC=pDTC+1;
}

pC1=C1;
pC2=C2;
pC=C;
pDTC=DTC;

conditionslimites();
dt_c();

	/**************etape3***************/

for(i=0;i<ne;i++)
{

(*pC)=(*pC)+dt*((*pDTC)/3.);
(*pC1)=(*pC2)+dt*(*pDTC);

pC=pC+1;
pC1=pC1+1;
pC2=pC2+1;
pDTC=pDTC+1;
}

pC1=C1;
pC2=C2;
pC=C;
pDTC=DTC;

conditionslimites();
dt_c();

	/**************etape4***************/

for(i=0;i<ne;i++)
{

(*pC)=(*pC)+dt*((*pDTC)/6.);

pC=pC+1;
pC1=pC1+1;
pC2=pC2+1;
pDTC=pDTC+1;
}

pC1=C1;
pC2=C2;
pC=C;
pDTC=DTC;

}


/***********************boucle principale de calcul selon la methode d'Euler*****************************/
/********************************************************************************************************/

	/********** Creation  ********/

void boucle_temps_1(void)
{

int i;
double bruit;

srand48(time(NULL));

for(i=0;i<n;i++) 
{
bruit=drand48();
*(C+i-borne_integ)=0.01*bruit;

}

for(i=0;i<nt;i++)
{

rungekutta();
conditionslimites();

} // fin de la boucle sur le temps

}
	/********** RUNGEKUTTA II ********/

void boucle_temps_2(void)

{

int u,i,j,k,l;
double x;
char diese;

printf("nom du fichier initial:\n");
scanf("%s",&choix1);

pf1 = fopen(choix1,"r");

fscanf(pf1,"%c %d",&diese,&temps_init);
temps_init=temps_init+temps_max;

	for(u=0;u<ne;u++)
	{
		
	fscanf(pf1,"%lf %lf\n", &x ,(C+u));
	
	} // fin de la boucle d'écriture dans pf1

//conditionslimites();

for(i=0;i<nt;i++)
{

rungekutta();
conditionslimites();

}

} //fin de la boucle Euler 2

/***********************liberation de la memoire occupée par les vecteurs********************************/
/********************************************************************************************************/

void liberation(void)
{

C=C-borne_integ;
C1=C1-borne_integ;
C2=C2-borne_integ;
DTC=DTC-borne_integ;

free(C);
free(C1);
free(C2);
free(DTC);

}

/***********************fonction principale**************************************************************/
/********************************************************************************************************/

main()

{

double x;
int u;
int test;

// initialisations de paramètres

initialisation();

// allocations

allocations();

	/**** Demande de l'endroit de depart ****/

test=0;
while(test==0)

{

printf("\n");
printf("partir de t=0             ----> tapez 0\n");
printf("partir d'un fichier existant ----> tapez 1\n");
printf("\n");

scanf("%d",&choix0);

	if(choix0==0)
	{
	temps_init=0;
	
	boucle_temps_1();
	
	test=1;
	
	}
	
	else if(choix0==1)
	{
	
	boucle_temps_2();	
	
	test=1;
	
	}

	else
	{
	printf("erreur de frappe!!!\n");
	}

}

	/**** Affichage des resultats dans pf2 ****/

printf("nom du fichier final:\n");
scanf("%s",&choix2);

pf2 = fopen(choix2,"w");

temps_init=temps_init+temps_max;

fprintf(pf2,"# %d\n",temps_init);

for(u=0;u<ne;u++)
{
	x=u*dl;	
	
	fprintf(pf2,"%lf %lf\n", x , *(C+u));

}

if(choix0==1)
{

fclose(pf1);

}

fclose(pf2);

// liberation

liberation();

}// fin du main