#include <iostream>
#include <fstream>
#include <math.h>
#include "prototypes.h"
using namespace std;

//          Parameters in the PDE

#define D       .5 /* The diffusion coefficient */

//          Solution values to be computed

#define XMAX     1 /* The maximum value of x, for plotting */
#define XMIN    -1 /* The minimum value of x, for plotting */
#define T        1 /* The final time of the computation    */

//          Parameters for the PDE solver

#define DX      .2 /* The step size in x for PDE solving, might not get exactly this many   */
#define DXRATIO .5 /* A parameter for the PDE solver       */

int main() {

	ofstream OutCake("C:\\Documents and Settings\\CIMS\\Desktop\\outcake.csv");
	if ( !OutCake ) cout << "No OutCake" << endl;

	OutCake << " D " << " , " << " T " << endl
		    <<   D   << " , " <<   T   << endl;

	double* u1; // Arrays to store the computed PDE solution
	double* u2; // Arrays to store the computed PDE solution
	int     nx;      // The number of x values between XMIN and XMAX
	int     nu;      // The size of the u array, for the requested values and those outside.
	int     nt;      // The number of time steps.

	double xrange = (double) XMAX - (double) XMIN;

	nx = (int) ( ( xrange /DX ) + 1 ); // Compute the number of grid points inside
	                                                    // the region asked for.
	double dx = xrange /( nx - 1 );  //The actual dx possibly not DX.
	double dt = DXRATIO*dx*dx/(2*D) ; // The time step, computed by the formula
	nt = (int) ( ( (double) T ) / dt ) +1;
	dt  = ( (double) T )/ nt;   //  Possibly different from the original dt by a little.
	nu  = nx + 2*nt; // Number of u values is the number asked for plus a buffer of size nt on
		            // either side.

	OutCake << " nt " << " , " << " nx " << " , " << " dt " << " , " << " dx " << " , " << " xrange " << endl
            <<   nt   << " , " <<   nx   << " , " <<   dt   << " , " <<   dx  << " , " <<    xrange   << endl;

    u1 = new double[nu];
	u2 = new double[nu];

    double xMin = (double) XMIN;
	int ec;  // the error code = the return value of the major routines.
	   ec =
	FinalValues( u1, xMin-(dx*nt), dx, nu);

	for ( int iu = 0; iu < nu; iu++ ) OutCake << u1[iu] << " , " ;
	OutCake << endl;

// The time stepping loop, which is the point of all this

	int it = 0;
	while ( it < nt ) {       // The complicated logic alternates between (u2 <- updated u1 ) and
		   ec =               // ( u1 <- updated u2 ).
	    TimeStep( &u2[it], &u1[it], dx, dt, D, nu-it );
	    it++;
		if ( it == nt ) break;
		   ec =
	    TimeStep( &u1[it], &u2[it], dx, dt, D, nu-it );
		it++; }

	// put the results into the .csv file.  First the x values then the uvalues.

	double x = xMin;
	int ixMax = nu - nt;  // The largest x index that has a computed value at the final time.
	for ( int ix = nt; ix < ixMax; ix++ ){
		OutCake << x << " , ";
		x += dx; }
	OutCake << endl;
	for ( int ix = nt; ix < ixMax; ix++ )
		OutCake << u1[ix] << " , ";
	OutCake << endl;

	OutCake << "all done" << endl;

	return 0;
}