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

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


double L = 1.0, lambda, kappa, ell, thetan; //lambda = 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 eigenfm1(double x){
  return(sqrt( (2.0/lambda)/(1.0-exp(-(2.0/lambda))) )*exp(-x/lambda) );
}



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

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

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


}

double initialcondition(double x){// condition initiale
  double norm = sqrt(2.0*PI*sqr(sigma));
  double expo = sqr(x-0.5*L)/(2.0*sqr(sigma));

  return(eigenfm1(x));
  //  return( exp(-expo)/norm ); //la gaussienne, condition initiale sur rho

  //    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( eigenfunction(1,x) +  eigenfunction(2,x) ); //testing an  eigenfunction as initial condition
}







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

    double x,t;
    
    //scanf("%lf", &L);

    double beta; 
    scanf("%lf", &lambda);
    scanf("%lf", &beta);
    scanf("%d", &N);
    int Npas; 
    scanf("%d",&Npas);
    double Dt;
    scanf("%lf",&Dt);

    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
    
    //Evolution dans le temps pour l'équation d'ondes

    //Calcul des coefficients à partir de la condition initiale
    //===========================================================
    double C0, K[N+2]; //cas spécial, pour dphi(u,0)/dtau=0
    for(i=0;i<(N+2);i++){
      x = i*(L/(N+1));
      f[i] = initialcondition(x)*eigenfm1(x);
    }
    //    C0=integrate(f,s,0.0,L,N)/2.0;//expression valable pour beta=0
    C0=integrate(f,s,0.0,L,N);//expression valable pour beta= non-nul

    //    C0 = 1.0; // On choisit comme condiiton initiale eigenfm1(x)

    for(n=1;n<(N+2);n++){
      for(i=0;i<(N+2);i++){
	x = i*(L/(N+1));
	f[i] = initialcondition(x)*eigenfunction(n,x);
      }
      K[n]=integrate(f,s,0.0,L,N);

      //      K[n]=0.0; // La condition initiale a récouvrement zéro avec les fn
    }

    //    K[1] = 1.0;

    //Calcul des coefficients de la condition initiale pour la partie non-linéaire

    double b0, b[N+2];

    for(i=0;i<(N+2);i++){
      x = i*(L/(N+1));
      f[i] = sin(initialcondition(x))*eigenfm1(x);
    }
    //    C0=integrate(f,s,0.0,L,N)/2.0;//expression valable pour beta=0
    b0=integrate(f,s,0.0,L,N);//expression valable pour beta= non-nul

    //    C0 = 1.0; // On choisit comme condiiton initiale eigenfm1(x)

    for(n=1;n<(N+2);n++){
      for(i=0;i<(N+2);i++){
	x = i*(L/(N+1));
	f[i] = sin(initialcondition(x))*eigenfunction(n,x);
      }
      b[n]=integrate(f,s,0.0,L,N);

      //      K[n]=0.0; // La condition initiale a récouvrement zéro avec les fn
    }

    //    K[1] = 1.0;

    //Reconstruction de sin(phi(u,0)) à partir de la série 
    for(i=0;i<(N+2);i++){
      x = i*(L/(N+1));
      somme = b0*eigenfm1(x);
      for(n=1;n<(N+2);n++){
	somme = somme + b[n]*eigenfunction(n,x);
      }
      //      printf("%g\t%g\t%g\n",x,sin(initialcondition(x)),somme);
    }


    //==========================================================
    //Evolution dans le temps
    double Dx = L/(N+1);

    double Omega, wn;
    double somme0, somme1, somme_dot0, somme_dot1;
    double Jbord; //le courant au bord
    double C0_temp, C0_hat, K_temp[N+2], K_hat[N+2];

    //    Omega=sqrt(-beta -sqr(1.0/lambda));

   

    //    for(j=0;j<(21*N);j=j+80){
    for(j=0;j<Npas;j++){
      t = j*Dt;
 
      // Mise à jour des C0 et C_hat0;
      C0_temp = C0 + Dt*C0_hat; 
      C0_hat = C0 + Dt*( (1.0/sqr(lambda))*C0-beta*b0 );
       
      for(n=1;n<(N+2);n++){
	K_temp[n] = K[n] + Dt*K_hat[n];
	K_hat[n] = K_hat[n] +  
	  Dt*(-eigenvalue(n)*K[n]-beta*b[n] );
      }
   
      //Mise à jour des coefficients
      C0 = C0_temp; 
      for(n=1;n<(N+2);n++){
	K[n] = K_temp[n];
      }

      //Calcul de phi(u) au nouvel instant
      for(i=0;i<(N+1);i++){
	x = i*Dx;
	somme = C0*eigenfm1(x);
	for(n=1;n<(N+1);n++){
	  somme = somme + K[n]*eigenfunction(n,x);
	}
      }

      //Calcul des nouvelles valeurs pour b0 et b[n]

      for(i=0;i<(N+2);i++){
	x = i*(L/(N+1));
	f[i] = sin(initialcondition(x))*eigenfm1(x);
      }
      b0=integrate(f,s,0.0,L,N);//expression valable pour beta= non-nul
      
      for(n=1;n<(N+2);n++){
	for(i=0;i<(N+2);i++){
	  x = i*(L/(N+1));
	  f[i] = sin(initialcondition(x))*eigenfunction(n,x);
	}
	b[n]=integrate(f,s,0.0,L,N);
      }

      //Calcul du courant de bord
      x = 0.0;
      somme0 = C0*eigenfm1(x);
      somme_dot0 = C0_hat*eigenfm1(x);
      for(n=1;n<(N+1);n++){
	somme0 = somme0 + K[n]*eigenfunction(n,x);
	somme_dot0 = somme_dot0 + K_hat[n]*eigenfunction(n,x);
      }

      x = 1.0;
      somme1 = C0*eigenfm1(x);
      somme_dot1 = C0_hat*eigenfm1(x);
      for(n=1;n<(N+1);n++){
	somme1 = somme1 + K[n]*eigenfunction(n,x);
	somme_dot1 = somme_dot1 + K_hat[n]*eigenfunction(n,x);
      }

      Jbord = (somme_dot1*somme1-somme_dot0*somme0)/lambda;
      printf("%g\t%g\n",t,Jbord);



      /*

   //      for(i=0;i<=(N+1);i++){
      //	x = i*Dx;

      //le courant à x = 0
      x = 0.0;
      //	somme = C0*cos(t*Omega)*eigenfm1(x);
      somme0 = C0*cos(t*Omega)*eigenfm1(x);
      somme_dot0 = -C0*Omega*sin(Omega*t)*eigenfm1(x);

	for(n=1;n<(N+1);n++){
	  wn = sqrt(eigenvalue(n)-beta);
	  //	  somme = somme + K[n]*cos(wn*t)*eigenfunction(n,x);
	  somme0 = somme0 + K[n]*cos(wn*t)*eigenfunction(n,x);
	  somme_dot0 = somme_dot0 -K[n]*wn*sin(wn*t)*eigenfunction(n,x);
	}

      x = 1.0;
      //	somme = C0*cos(t*Omega)*eigenfm1(x);
      somme1 = C0*cos(t*Omega)*eigenfm1(x);
      somme_dot1 = -C0*Omega*sin(Omega*t)*eigenfm1(x);

	for(n=1;n<(N+1);n++){
	  wn = sqrt(eigenvalue(n)-beta);
	  //	  somme = somme + K[n]*cos(wn*t)*eigenfunction(n,x);
	  somme1 = somme1 + K[n]*cos(wn*t)*eigenfunction(n,x);
	  somme_dot1 = somme_dot1 -K[n]*wn*sin(wn*t)*eigenfunction(n,x);
	}

	Jbord = (somme_dot1*somme1-somme_dot0*somme0)/lambda;
	
	//	printf("%g\t%g\t%g\n",t,x,sqr(somme));

	//      } //fin de boucle sur l'espace, étiquetté par i

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

    } //fin de boucle sur le temps, étiquetté par j
	    
}