#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)
//#define lambdaplus 5.0 

double L = 1.0, D = 1.0, lambda, kappa, ell, thetan;


/*

double L, D = 0.01; double c1 = 1.0, c2 = 1.0;

void source(double b[], int N){
  
  int n,m1;
  for(n=1;n<N;n++){
    m1=1-2*(n%2);
    thetan = atan(L/(n*lambdaplus*PI));
    b[n]=-c1*(sqrt(2.*L*sqr(L))/(n*PI))*( sin(thetan)*m1+((1-m1)/(n*PI))*cos(thetan) );
  }
}

*/

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( 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);

    scanf("%lf", &lambda);
    scanf("%d", &N);

    double a[N+2],b[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
    
        //some tests

    //Testing the orthonormality of the eigenfunctions 

    /*    scanf("%d %d", &n, &m);

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

    printf("%d\t%d\t%g\n",n,m,integrate(f,s,0.0,L,N));
    */

    /*    for(i=0;i<(N+2);i++){
      x = i*(L/(N+1));
      printf("%g\t%g\n",x,eigenfunction(1,x));
    }
    */
    

    //afficher la condition initiale

    /*    for(i=0;i<(N+2);i++){
      x = i*(L/(N+1));
      printf("%g\t%g\n",x,initialcondition(x)); OK!
    }
    */


    //computing the coefficient of the initial condition, that belongs to finite spectrum

    double am1;

    for(i=0;i<(N+2);i++){
      x = i*(L/(N+1));
      f[i] = initialcondition(x)*eigenfm1(x);
    }
    am1 = integrate(f,s,0.0,L,N);
    
    //computing the coefficients of the initial condition, that belong to the infinite spectrum

            for(n=0;n<(N+2);n++){
	      for(i=0;i<(N+2);i++){
		x = i*(L/(N+1));
		f[i] = initialcondition(x)*eigenfunction(n,x);
	  }
	      a[n]=integrate(f,s,0.0,L,N);
      //      printf("%d\t%g\n",n,a[n]);
      
	}
    

   //check the reconstruction of the initial condition


	    /*	    
	    	   for(i=0;i<(N+2);i++){
	  x = i*(L/(N+1));
	  somme = am1*eigenfm1(x);
	  for(n=0;n<(N+2);n++){
	    somme = somme + a[n]*eigenfunction(n,x);
	  }
	  printf("%g\t%g\n",x,somme);
	}
	    */
	    
	    
    

    //Time evolution-toy model

	    double Dtau = 1e-3;
	    int jinit, jfinal, jincr;
	    scanf("%d %d %d", &jinit, &jfinal, &jincr);

	    double omegan, somme0, sommed0, sommeL, sommedL, JRobinbord; 
	    for(j=jinit;j<jfinal; j=j+jincr){
	      t = j*Dtau;
	      //	      for(i=0;i<(N+2);i++){
	      x = 0.0;

	      somme0 = am1*eigenfm1(x); sommed0 = 0.0;
		for(n=0;n<N;n++){
		  omegan=sqrt( eigenvalue(n)+sqr(1./lambda) );
		  somme0 = somme0 + a[n]*cos(omegan*t)*eigenfunction(n,x);
		  sommed0 = sommed0 - a[n]*omegan*sin(omegan*t)*eigenfunction(n,x);
		}

		x = L;

		sommeL = am1*eigenfm1(x); sommedL = 0.0;
		for(n=0;n<N;n++){
		  omegan = sqrt( eigenvalue(n)+sqr(1./lambda) );
		  sommeL = sommeL + a[n]*cos(omegan*t)*eigenfunction(n,x);
		  sommedL = sommedL - a[n]*omegan*sin(omegan*t)*eigenfunction(n,x);
		}
		JRobinbord=(sommeL*sommed0-somme0*sommedL)/lambda;
		printf("%g\t%g\n",t,JRobinbord);
		//		printf("%g\t%g\t%g\n",t,x,somme);
		//	      }

	    }

















    //    source(b,N);

	    /*
		   
	    	    double Dtau = 1e-3;

	   	    int jinit, jfinal, jincr;
	    scanf("%d %d %d", &jinit, &jfinal, &jincr);

	    for(j=jinit;j<jfinal;j=j+jincr){
	      t = j*Dtau;
	      for(i=0;i<(N+2);i++){
		x = i*(L/(N+1));
		somme = am1*eigenfm1(x);
		for(n=0;n<N;n++){
		  somme = somme + a[n]*exp(-eigenvalue(n)*t-(t/sqr(lambda)))*eigenfunction(n,x);
	}

		printf("%g\t%g\t%g\n",t,x,exp(-x/lambda)*somme);//afficher rho

      }

    }
	    */	   

	    /*	    scanf("%lf", &t);

	    for(i=0;i<(N+2);i++){
	      x = i*(L/(N+1));
	      somme = am1*eigenfm1(x);
	      for(n=0;n<N;n++){
		somme = somme + a[n]*exp(-eigenvalue(n)*t-(t/sqr(lambda)))*eigenfunction(n,x);
	      }

	      printf("%g\t%g\t%g\n",x,somme,am1*eigenfm1(x));

	    }
	    */

    
}