//**  Hello,
//**  I put couple of things here.  You still need to make several other
//**  things.  I don't want you spend too much time in designing those
//**  functions.  Try and ask me for help if you can not finish it.
//**
//**
//**  I have corrected couple of things and now this file is complete.
//**  Of course, you need to make other things to complete the project.

//**************************************************************************
//****  Following is a sample program  named 'sample.C'.
//****  Compile the program  (Command is : g++ sample.C -lm )
//****  and run it  (Command is : a.out  )
//**************************************************************************


//**************************************************************************
//**** Some files to be included
//**************************************************************************
#include <iostream.h>
#include <fstream.h>
#include <stream.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>

//**************************************************************************
//**** Constants to be used in the scheme
//**************************************************************************
const double L1=-PI;
const double L2=PI;
const double alpha=0.2;                   // 
const double t1=alpha*alpha/4;

const double beta=2;                      // 
const double t2=beta*beta/4;

const double eta=0.25;

//**************************************************************************
//**** Exposure we want to have
//**************************************************************************
double E(double x){
  if(x<-0.5) return(-0.0);
  if(x<0.5)  return(1.0);
  return(0.0);
}

//**************************************************************************
//****  Now we have to compute the coefficients of Fourier seires 
//****  given by E(x).  They are given in terms of integrations.
//****  Numerically, integration is a summation of small rectangles.
//**************************************************************************

//************************************************************
//****  We introduce constants to be used in the integration.
//************************************************************
const int MAX=126;                        // MAX=(L2-L1)/dx
const double dx=(L2-L1)/(MAX);            // dx
const int N=11;

//************************************************************
//****  Coefficients of Fourier seires
//****  I have used notation a[n] in the class.
//****  In the book, p. 58, they use e[n].
//************************************************************
double a[N];  //these coefficients are computed in the main function.

double EN(double x){
  double result;
  int n;

  result=a[0]*0.5;
  for(n=1;n<N;++n)
    result=result+a[n]*cos(n*x);

  return(result);
}

//**************************************************************************
//****  Now you have to choose your own D(x).
//**************************************************************************
  //**************************************************************************
  //****  First coefficients.
  //**************************************************************************
double d[N]; //these coefficients are computed in the main function.

  //**************************************************************************
  //****  Now D(x).  Following is a simple case. 
  //**************************************************************************
double D(double x){
  double result;
  int n;

  result=d[0]*0.5;
  for(n=1;n<N;++n)
    result=result+d[n]*cos(n*x);

  return(result);
}

  //**************************************************************************
  //****  Now actual expose AE(x) : This is the function given by
  //****  the formulat (3.4)
  //**************************************************************************
double K(double x1,double x2,double index){
  return(exp(-(x1*x1+x2*x2)/index/index)/PI/index/index);
}

double AE(double x){    //  Remember that (3.4) in the book is double integration.
  double y1,y2,result;
  int j,l;

  result=0;
  for(l=-0*MAX;l<1*MAX;++l){
    y2=L1+(l+0.5)*dx;
    for(j=-0*MAX;j<1*MAX;++j){
      y1=L1+(j+0.5)*dx;
      result=result+((K(x-y1,0-y2,alpha)+eta*K(x-y1,0-y2,beta))*D(y1))*dx*dx/(1+eta);
    }
  }
  return(result);
}

//************************************************************
//**** In my main function I print out those functions.
//************************************************************
void main(){
  filebuf fi;
  int j,n;
  double x;

  ofstream output1("newE", ios::out);
  ofstream output2("newEN", ios::out);
  ofstream output3("newD", ios::out);
  ofstream output4("newAE", ios::out);

for(n=0;n<N;++n){
  a[n]=0;
  for(j=0;j<MAX;++j){
    x=L1+(j+0.5)*dx;
    a[n]=a[n]+E(x)*cos(n*x)*dx;
  }
    a[n]=a[n]/PI;
}

for(n=0;n<N;++n){
   d[n]=a[n]*(1+eta)/(exp(-(n*n*t1))+eta*exp(-(n*n*t2)));
}
  for(x=-2;x<2;x=x+dx){
     output1<<form("%f %f", x,E(x))<<endl;
     output2<<form("%f %f", x,EN(x))<<endl;
     output3<<form("%f %f", x,D(x))<<endl;
//     output4<<form("%f %f", x,AE(x))<<endl;     // Note that AE and EN should be same.
  }
cout<<form("t1=%f, t2=%f",t1,t2)<<endl;
}