void timeStep()
{
	double t0 = t;
	do {
	//	int x = 31;
	//	int y = 61;
		cout << "t = " << t << endl;
		if(debug)
			cout << t0 << endl;
		for(int y = 0; y<sizeY; y++)
		{
			cout << t << endl;
			for(int x = 0; x<sizeX; x++)
			{
				if(debug)
					cout << "c" << endl;
				grid [x][y] = grid[x][y] * pie(x,y); //so P(t2) = P(t1) * somejunk
			}
		}
		t+=tau;
		cout << t << endl;
	} while (abs(c0 * (t - t)) > sizeX / T / framesPerSec); 
    glutPostRedisplay();
    glFlush();
}

using namespace std;
//using namespace itpp;

//////////////////////////////////VARIABLE DECLARATION////////////////////////////////////////////
double T = 5;
double k = 2; //determines the speed of motion, used in e^ikx
bool debug = false; //True for debugging, used mainly in glut but also gives nice console output every step of the way.
double t; //the time variable
double m = 1; //mass
double c0 = 1; //beginning speed
double xcenter = 32; //x coordinate for the center of the pressure wave
double ycenter = 64; //y coordinate for the center of the pressure wave
double sigmaXY = 5; //determines the width of the pressure wave in the x
double sigma0 = 6; //determines the sigma for the daf
double tau = 0; //time step as determined by the nyquist condition, meant to be less than nyquist
double temp; //Temp double for various calculations, resets as needed
double speed; //velocity profile
double delX = 1; //Grid spacing
double tempP[upperM]; //The values of the array for the Pk 
double bar = 0.1; //This is the difference between xi and xj, yi and yj, etc. Can be tuned for more accuracy.
double ** grid;
double tMax = 10000; //This is tMax, set to something needlessly big so it runs infinitely
double framesPerSec = 60;       // animation rate for screen redraws
bool running; //controls animation
double expDelT;
double * dafX;
double * dafY;

void init()
{	if(debug)
	cout << "a" << endl;
	t = 0; //resets time
	double upperBound = 20; //This is the upper Range for the daf calculations
	tau = (c0*c0)/(k*k*100);
	grid = new double*[sizeX];

	for (int asdf = 0; asdf < sizeX; ++asdf)
		grid[asdf] = new double[sizeY];
	
	//Generates your gaussian
	
	for(int y = 0; y<sizeY; y++)
	{
		for(int x = 0; x<sizeX; x++)
		{
			grid[x][y] = 0;
		}
	}
	
	expDelT = exp(-(tau*tau)); //this is the e^-deltaT^2

	//Generates your DAF Kernel and stores it as an array for the X and an array for the Y
	dafX = new double[sizeX];
	dafY = new double[sizeY];
	
	for(int iterDaf = 0; iterDaf<sizeX; iterDaf++)
	{
		dafX[iterDaf] = deltaX(upperBound, iterDaf);
	}
	
	for(int iterDaf = 0; iterDaf<sizeY; iterDaf++)
	{
		dafX[iterDaf] = deltaY(upperBound, iterDaf);
	}
	
	//Does all the clearing that is needed
	cout << "done finally" << endl;
    glClearColor(1.0, 1.0, 1.0, 0.0);//Define our background color
    gluOrtho2D(0,sizeY,0,sizeX);//(NEW)  Define our viewing area
	glDisable(GL_DEPTH_TEST);
	glEnable(GL_POINT_SMOOTH);
	glPointSize(5);
	//Generates your gaussian
	
	t = 0;
	
	for(int y = 0; y<sizeY; y++)
	{
		for(int x = 0; x<sizeX; x++)
		{
			grid[x][y] = exp(-(((x-xcenter)*(x-xcenter))/(2*sigmaXY*sigmaXY)+((y-ycenter)*(y-ycenter))/(2*sigmaXY*sigmaXY)));
		}
	}
	
}

//Definition:
T a=0;
//Declaration (will cause errors if in the same file as the definition):
extern T a;