// Résolution de l'équation de Schrödinger par une méthode spectrale
// On développe sur la base de fonctions propres de -d²/dx² = H_0
// Conventions: hbar = 1 et m = 0.5

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

#define sqr(x) ((x)*(x))
#define L 1.0
#define T 256.0 
#define PI acos(-1.0L)

#define LPOINTS 20
#define TPOINTS 40

#define Nsimp 127

#define ax  ((L)/(LPOINTS))
#define at  ((T)/(TPOINTS))

#define sigma (0.1*(L)) // largeur de la condition initiale
#define pos_moy (0.5*(L))
#define pos_moy1 (0.25*(L)) // centre de la condition initiale
#define pos_moy2 (0.75*(L)) // centre de la condition initiale

double simpson( double (*f)(double), double A, double B, int N){
//Méthode de Simpson pour l'intégration numérique 
    double h=(B-A)/( (double)N +1.0);
    double temp=h*(f(A)+f(B))/3.0;
    int j;
    double x;
    for(j=1;j<N;j=j+2){
	x=A+j*h;
	temp=temp + (4.0*h*f(x)/3.0);}

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



double eigenvalue(int n){// Valeurs propres indexées par un entier n
    return(sqr(n*PI/L));}

double eigenfunction(int n, double x){// fonctions propres indexées par n
    return(sqrt(2.0/L)*sin((double)(n*PI*x/L)));}

double initialcondition(double x){// condition initiale
  return(exp(-(sqr(x-pos_moy))/(2.0*sqr(sigma)))/sqrt(2.0*PI*sqr(sigma)));

  /*  double temp1 = exp(-(sqr(x-pos_moy1))/(2.0*sqr(sigma)))/sqrt(2.0*PI*sqr(sigma));

  double temp2 = exp(-(sqr(x-pos_moy2))/(2.0*sqr(sigma)))/sqrt(2.0*PI*sqr(sigma));

  return(0.5*(temp1+temp2));*/
}

double overlap(int n){/* calcule les coefficients du développement 
				   en fonctions propres de la solution */
  double f(double x){
    return(initialcondition(x)*eigenfunction(n,x));}
  
  return( simpson(f,0.0, L, Nsimp) ); // varier le nombre de points intermédiaires (ici Nsimp)
}


double complex psi(double x, double t){
  double complex somme = 0.0;
  int n, N = 16; 
  double cn; 

  for(n=1;n<N;n++){
    cn = overlap(n);
    somme = somme + cn*eigenfunction(n, x)*cexp(-I*eigenvalue(n)*t);
  }
  return(somme);
}

double density(double x, double t){
  return( sqr(cabs(psi(x,t))) );
}

double norm(double t){

  double f(double x){
    return(density(x,t));}

  return(simpson(f,0.0,L, Nsimp));
}


double cdm(double t){//centre de masse de la densité
  double f(double x){
    return(x*density(x,t));}
  return(simpson(f,0.0,L,Nsimp)/norm(t));
}


double largeur(double t){
  double f(double x){
    return(sqr(x-cdm(t))*density(x,t));
  }
  return( simpson(f,0.0,L,Nsimp)/norm(t) );
}



main(){
    int i, j;
    double x,t;

    
    

    for(i=0;i<=TPOINTS;i=i+(TPOINTS/8)){

	t = i*at;
	//		for(j=0;j<=LPOINTS;j++){
	//	    x = j*ax;
	    //	    	    printf("%g\t%g\t%g\n",t,x, density(x,t));
	//	}
	

	//	printf("%g\t%g\n", t, cdm(t));
	//	printf("%g\t%g\n", t, norm(t));

	printf("%g\t%g\n",t, largeur(t)); 
    }


}