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

#define sqr(x) (x)*(x)
#define sigma 0.1
#define PI acos(-1.0)




void simpson(double s[], int N){//set the coefficients of the numerical integration procedure for Simpson's rule: N=odd
  int i;
  s[0] = 1./3.;
  for(i=1;i<=N;i=i+2){
    s[i]=4./3.;
  }
  for(i=2;i<N;i=i+2){
    s[i]=2./3.;
  }
  s[N+1]=1./3.;
}

void trapezes(double s[], int N){// set the coefficients of the numerical integration procedure for the trapezoidal rule; N can be even or odd
  int i;
  s[0]=0.5;
  for(i=1;i<(N+1);i++){
    s[i]=1.0;
  }
  s[N+1]=0.5;
}

double integrate(double f[], double s[], double A, double B, int N){//numerical integration: 
  int i; double somme=s[0]*f[0];

  for(i=1;i<=N;i++){
    somme = somme + s[i]*f[i];
  }
  somme = somme + s[N+1]*f[N+1];
  return( ((B-A)/(N+1))*somme);
}


double valpropre(int n){// Valeurs propres de la seconde dérivée spatiale
    return(sqr(n*PI));}


double fonpropre(int n, double x){// fonctions propres de la seconde dérivée spatiale avec conditions Robin
  return(sqrt(2.0)*sin(n*PI*x)); // eigenfunctions for Dirichlet boundary conditions

}

double initialcondition(double x){// condition initiale
  //  return(exp(-sqr(x-0.5*L)/(2.0*sqr(sigma)))/sqrt(2.0*PI*sqr(sigma))); //une gaussienne
  //  return(eigenfunction(1,x)); //testing an  eigenfunction as initial condition
}







main(){ 
    double complex somme;
    int n, l, j;
    int N=128;
    
    double c[N+1];

    //initialisation des coefficients
    for(n=0;n<=N;n++){
      c[n]=0.0;
    }

    //on se donne la condition initiale, en donnant des valeurs non-nulles à certains coefficients

     c[4]=1.0;
 
    double u,tau;
    double Du = 1.0/N;
    double Dtau = 1.0/N;
    
    int Tmax=sqr(N);

    for(j=0;j<=Tmax; j = j+10){ //boucle sur le temps

      tau = j*Dtau;

      for(l=0;l<=N;l++){ //boucle sur l'espace
	u = l*Du;

	somme = 0.0;

	for(n=1;n<=N;n++){// la boucle sur la somme sur les modes
	  somme = somme + c[n]*cexp(-I*valpropre(n)*tau)*fonpropre(n,u);
	
	}// on ferme la boucle sur la somme sur les modes

	printf("%g\t%g\t%g\n",tau,u,sqr(cabs(somme)));

      } // on ferme la boucle sur l'espace

    } // on ferme la boucle sur le temps
    

    




      
    
}