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

#define sqr(x) ( (x)*(x) )
//#define sigma 0.01
#define PI acos(-1.0)
//#define Gamma  5.0
#define Astar sqrt( 2.0/ ( Gamma*(1.0-exp( -2.0/Gamma ) ) ) ) 
#define L 1.0

double  thetan; 

double sigma, Gamma;

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 evaluestar(){
  return( -1.0/sqr(Gamma) );
}

double eigenfneg(double x){
  return(Astar*exp(-( (x/L)/Gamma )) );
}

double eigenfnegprime(double x){
  return(-(Astar/(L*Gamma))*exp(-x/(L*Gamma)) );
}

double eigenvalue(int n){// Valeurs propres de la seconde dérivée spatiale
  return(sqr(n*PI/L));} //choix d'unités: v=1 et L=1


double eigenfunction(int n, double x){// fonctions propres de la seconde dérivée spatiale avec conditions Robin
  
   thetan = 2.0*atan(Gamma*n*PI);
   return( sqrt(2.0/L)*sin( (n*PI*x/L)-0.5*thetan)  );
  

     //  return(sqrt(2.0/L)*sin(n*PI*x/L)); // eigenfunctions for Dirichlet boundary conditions

}

double eigenfunctionprime(int n, double x){

  thetan = 2.0*atan(Gamma*n*PI);
  return( sqrt(2.0/L)*(n*PI/L)*cos( (n*PI*x/L) - 0.5*thetan ) );
}

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(1.0); // une constante
  //  return(eigenfneg(x));//fonction propre de la valeur propre négative
  //   return(eigenfunction(1,x)); //testing an  eigenfunction as initial condition
}







main(){
    double somme;
    int n, i, j;
    int N;

    double u,tau;
    
    scanf("%d", &N);

    scanf("%lf %lf", &sigma, &Gamma);

    double cstar, c[N+2];

    double s[N+2], f[N+2];
    //     simpson(s,N); // set the coefficients of the numerical integration for Simpson's rule
        trapezes(s,N); // set the coefficients of the numerical integration for trapezes

    // Exercice 0: Retrouver que la fonction propre de valeur propre négative  est  correctement normalisée 

    for(i=0;i<=(N+1);i++){
      u = (float)i/(float)(N+1);
      //      f[i] = sqr(Astar)*exp(-2.0*u/Gamma);
      f[i] = sqr(eigenfneg(u));
    }

    //   printf("%17.14g\n",integrate(f,s,0.0,1.0,N));


// Exercice 1: Retrouver que les fonctions propres de valeur propre positive sont correctement normalisées 


    for(n=1;n<(N+1);n++){

    for(i=0;i<=(N+1);i++){
      u = (float)i/(float)(N+1);
      f[i] =  sqr(eigenfunction(n,u));
    }
    //    printf("%d\t%17.14g\n",n,integrate(f,s,0.0,1.0,N));

    }


    //Exercice 2: Orthogonalité des fonctions propres
    //(a) Orthoginalité entre la fstar et les fn

    for(n=1;n<(N+1);n++){
      for(i=0;i<=(N+1);i++){
	u = (float)i/(float)(N+1);
	f[i] = eigenfneg(u)*eigenfunction(n,u);
      }
      
      //      printf("%d\t%17.14g\n",n,integrate(f,s,0.0,1.0,N));

    }

    //Exercice 3: Réproduire la condition initiale

    //Calcul du coefficient cstar

    for(i=0;i<=(N+1);i++){
      u = (float)i/(float)(N+1);
      f[i] = initialcondition(u)*eigenfneg(u);
    }
    cstar = integrate(f,s,0.0,1.0,N);

    //Calculs des coefficients c[n] 

    for(n=1;n<(N+1);n++){
      for(i=0;i<=(N+1);i++){
	u = (float)i/(float)(N+1);
	f[i] = initialcondition(u)*eigenfunction(n,u);
      }
      c[n] = integrate(f,s,0.0,1.0,N);
    }

    //Test que la condition initiale est correctement réproduite

    for(i=0;i<=(N+1);i++){
      u = (float)i/(float)(N+1);
      somme = cstar*eigenfneg(u);
      for(n=1;n<(N+1);n++){
	somme = somme + c[n]*eigenfunction(n,u);
      }
      //            printf("%g\t%g\t%g\n",u,initialcondition(u),somme);
    }
    
    //Exercice 4: Evolution dans le temps, pour le cas spécial Omega = 0

    double w[N+2];
    for(n=1;n<(N+1);n++){
      w[n] = sqrt( sqr(n*PI) + sqr(1.0/Gamma) );
    }

    double T1 = 2.0*PI/w[1], TN = 2.0*PI/w[N]; 
    /* Le courant de bord de l'énergie pour le cas où la condition initiale est
       une fonction propre, appartenant au spectre positif

     */

    //Pour préciser un seul mode
    //    scanf("%d", &n); 

    double Jbord, dphidtau0,dphidtau1,dphidu0,dphidu1;

    for(tau=0.0; tau < 2*T1; tau = tau + 10*TN){
      //Jbord pour le cas de la condition initiale, qui serait un seul mode,  qu appartenant au spectre positif
      //      Jbord = (w[n]/Gamma)*sin(w[n]*tau)*( sqr(eigenfunction(n,0.0))-sqr(eigenfunction(n,1.0)) );

      //Le calcul de Jbord, dans le cas d'une condition initiale générale
      dphidtau0 = 0.0; dphidtau1 = 0.0;
      dphidu0 = cstar*eigenfneg(0.0); dphidu1 = cstar*eigenfneg(1.0);

      for(n=1;n<(N+1);n++){
	dphidtau0 = dphidtau0 + c[n]*w[n]*sin(w[n]*tau)*eigenfunction(n,0.0);
	dphidtau1 = dphidtau1 + c[n]*w[n]*sin(w[n]*tau)*eigenfunction(n,1.0);

	dphidu0 = dphidu0 + c[n]*cos(w[n]*tau)*eigenfunction(n,0.0);
	dphidu1 = dphidu1 + c[n]*cos(w[n]*tau)*eigenfunction(n,1.0);

      }

      dphidtau0 = -dphidtau0;
      dphidtau1 = -dphidtau1;

      dphidu0 = -dphidu0/Gamma;
      dphidu1 = -dphidu1/Gamma;


      Jbord = dphidu0*dphidtau0 - dphidu1*dphidtau1;

      //            printf("%g\t%g\n",tau, Jbord);
      
    }

    
    for(tau=0.0; tau <20*T1; tau = tau + 0.01*TN){

      /*
	Calculer la valeur de la solution à chaque bord, c.à.d. phi(0.,tau) et
	phi(1.,tau)
      */ 
      //=====================================================================
      u = 0.0;
      somme = cstar*eigenfneg(u);
      for(n=1;n<(N+1);n++){
	somme = somme + c[n]*cos(w[n]*tau)*eigenfunction(n,u);
      }
      printf("%g\t%g\t",tau, somme);

      u = 1.0;
      somme = cstar*eigenfneg(u);
      for(n=1;n<(N+1);n++){
	somme = somme + c[n]*cos(w[n]*tau)*eigenfunction(n,u);
      }
      printf("%g\n",somme);
      //=====================================================================
      /*      for(i=0;i<=(N+1);i++){
	u = (float)i/(float)(N+1);
	somme = cstar*eigenfneg(u);

	for(n=1;n<(N+1);n++){
	  somme = somme + c[n]*cos(w[n]*tau)*eigenfunction(n,u);
	}
	printf("%g\t%g\t%g\n",tau, u, somme);
      }
      */

    }
    

}