/* Code pour la simulation de l'équation d'ondes à d=1 */

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


#define sqr(z) (z)*(z)
#define delta(z) fabs((float)(z))==0 ? 1 : 0
#define theta(z) (z)>=0 ? 1 : 0
#define pulse(a,b,z) 1-(theta((a-1)-(z))+theta((z)-(b+1)))
#define pi 3.1415926535898 
#define v 10.0
#define v2 (v)*(v) 
#define L 20.0
#define N 2000
#define a (L)/(N)
#define M 10000
#define T 20.0
#define h (T)/(M)
#define Inst (int)(2.0*(L)/( (v)*(h) ))
#define sigma 1.0
#define alpha2 0.0

double laplacian_1d(double *X, int k){
    double temp;
    temp=(X[k+1]-2.0*X[k]+X[k-1])/sqr(a);
    return(temp);
}


void initial(double *X, double *Y){
    int n,m;
    double z,u,s;
    for(n=0;n<(N+1);n++){
	z=n*a;
	Y[n]=0.0;
	X[n]=sin(2.0*acos(-1.0L)*z/L);
//	X[n]=sin(pi*z/L);
//	X[n]= exp(-sqr(z-0.5*L)/sigma);
//	X[n]=pulse(0.25*L, 0.75*L, z);
    }

}



void bc(int i, double *X, double *Y){
    double t=i*h;
    Y[0]=0.0; Y[N]=0.0;
    X[0]=0.0; X[N]=0.0;
}

double Vder(double *X, double g, double b, int k){
    return(g*X[k]+b*X[k]*sqr(X[k]));
}

double F(double *X, double *Y, int k){
    double temp;
    temp=Y[k];
    return(temp);
}

double G(double *X,double *Y, int k){
    double temp;
    temp=v2*laplacian_1d(X,k)+Vder(X,alpha2, 0.0, k);
    return(temp);
}

// Runge Kutta for the wave equation
// write: dX/dt=F(X,Y)=Y, dY/dt=G(X,Y)=v^2*laplacian_1d(X)
// Use the K variables with X, the L variables with Y

void rk4(double *X, double *Y){
    int k;
    double temp_X[N+1], temp_Y[N+1], temp_K[4][N+1], temp_L[4][N+1];
 
    for(k=1;k<N;k++){
	temp_K[0][k]=h*F(X,Y,k);
	temp_L[0][k]=h*G(X,Y,k);
    }

    for(k=1;k<N;k++){
	temp_X[k]=X[k]+0.5*temp_K[0][k];
	temp_Y[k]=Y[k]+0.5*temp_L[0][k];
    }

    for(k=1;k<N;k++){
	temp_K[1][k]=h*F(temp_X, temp_Y, k);
	temp_L[1][k]=h*G(temp_X, temp_Y, k);
    }

    for(k=1;k<N;k++){
	temp_X[k]=X[k]+0.5*temp_K[1][k];
	temp_Y[k]=Y[k]+0.5*temp_L[1][k];
    }

    for(k=1;k<N;k++){
	temp_K[2][k]=h*F(temp_X, temp_Y, k);
	temp_L[2][k]=h*G(temp_X, temp_Y, k);
    }

    for(k=1;k<N;k++){
	temp_X[k]=X[k]+temp_K[2][k];
	temp_Y[k]=Y[k]+temp_L[2][k];
    }

    for(k=1;k<N;k++){
	temp_K[3][k]=h*F(temp_X, temp_Y, k);
	temp_L[3][k]=h*G(temp_X, temp_Y, k);
    }

    for(k=1;k<N;k++){
	X[k]=X[k]+(1./6.)*(temp_K[0][k]+2.0*( temp_K[1][k]+temp_K[2][k])+temp_K[3][k]);
	Y[k]=Y[k]+(1./6.)*(temp_L[0][k]+2.0*( temp_L[1][k]+temp_L[2][k])+temp_L[3][k]);
    }

}



void rk4_2(double *X, double *Y){
    int k;
    double temp_X[N+1], temp_Y[N+1], temp_K[4][N+1], temp_L[4][N+1];
 
    for(k=1;k<N;k++){
	temp_K[0][k]=h*Y[k];
	temp_L[0][k]=h*v2*laplacian_1d(X,k);
}
   
    for(k=1;k<N;k++){
	temp_X[k]=X[k]+0.5*temp_K[0][k];
	temp_Y[k]=Y[k]+0.5*temp_L[0][k];
    }


    for(k=1;k<N;k++){
	temp_K[1][k]=h*temp_Y[k];
	temp_L[1][k]=h*v2*laplacian_1d(temp_X,k);
    }


    for(k=1;k<N;k++){
	temp_X[k]=X[k]+0.5*temp_K[1][k];
	temp_Y[k]=Y[k]+0.5*temp_L[1][k];
    }


    for(k=1;k<N;k++){
	temp_K[2][k]=h*temp_Y[k];
	temp_L[2][k]=h*v2*laplacian_1d(temp_X,k);
    }

    for(k=1;k<N;k++){
	temp_X[k]=X[k]+temp_K[2][k];
	temp_Y[k]=Y[k]+temp_L[2][k];
    }

    for(k=1;k<N;k++){
	temp_K[3][k]=h*temp_Y[k];
	temp_L[3][k]=h*v2*laplacian_1d(temp_X,k);
    }

    for(k=1;k<N;k++){
	X[k]=X[k]+(1./6.)*(temp_K[0][k]+2.0*(temp_K[1][k]+temp_K[2][k])+temp_K[3][k]);
	Y[k]=Y[k]+(1./6.)*(temp_L[0][k]+2.0*(temp_L[1][k]+temp_L[2][k])+temp_L[3][k]);

    }


}





double simpson( double (*f)(double), double A, double B, int Ns){
    double h1=(B-A)/( (double)Ns +1.0);
    double temp=h1*(f(A)+f(B))/3.0;
    int j;
    double x;
    for(j=1;j<Ns;j=j+2){
	x=A+j*h1;
	temp=temp + (4.0*h1*f(x)/3.0);}

    for(j=2;j<Ns;j=j+2){
	x=A+j*h1;
	temp=temp + (2.0*h1*f(x)/3.0);} 
	return(temp);
    }

void output(int i, double *X){
    int n;
    for(n=0;n<(N+1);n++){
	printf("%g %g %g\n",i*h,n*a,X[n]);}
	
}


main(){
  double Phi[N+1], Psi[N+1];
    int n;
    int i; 




  
  /* condition initiale */

    initial(Phi,Psi);

/* conditions aux bords */

    bc(0,Phi,Psi);

/*    i=0;
      output(i,Phi);*/

/* évolution temporelle */

 for(i=1;i<M;i++){

	rk4(Phi,Psi);
	if(i%Inst==0){
	    output(i,Phi);
	    
	}
 }
  
   
}