#include "backTime_solver.h"





CBackTimeSolve :: CBackTimeSolve()

{

	pTable = NULL;

	analytic = NULL;

	right_border = NULL;

}



CBackTimeSolve :: ~CBackTimeSolve()

{

	clearArray();

	analytic = NULL;

	right_border = NULL;

}





double ** CBackTimeSolve::getTable()

{

	return pTable;

}





double * CBackTimeSolve::getGradient(double * prev_grad)

{

	for (int k =0; k < K; k++)

	{

		prev_grad[k] = -pTable[k][0];

	}

	return prev_grad;



}





/*

void ReverseSolve::CheckSolve(string f_name)

{

	ofstream out(f_name.c_str());

	out << fixed;

	float temp;

	for (int k = 0; k < K; k++)

	{

		for (int i = 0; i < N; i++)

		{

			temp = r(stX + i*h, k*tau);

			out << setprecision (5) << "t = " << k*tau << "   x = " << stX + i*h << "   u(x,y) = " << u_array[k][i] << "   u*(x,y) = " << temp << "  Difference = " << fabs(u_array[k][i] - temp) << endl;

		}

		out << endl << endl;

	}

}*/











void CBackTimeSolve:: clearArray()

{

	if (pTable != NULL)

	{

		for (int k = 0; k < K; k++)

		{

			delete [] pTable[k];

		}

		delete [] pTable;

	}

	pTable = NULL;

}





/*

void ReverseSolve::CheckResult(double **correct, string f_name)

{

	ofstream out(f_name.c_str());

	out << fixed;

	double temp;

	for (int k = 0; k < K; k++)

	{

		for (int i = 0; i < N; i++)

		{

			temp = fabs(correct[k][i] - u_array[k][i]);

			out << setprecision (5) << "t = " << k*tau << "   x = " << stX + i*h << "   u(x,y) = " << u_array[k][i] << "   u*(x,y) = " << correct[k][i] << "  Difference = " << temp << endl;

		}

		out << endl << endl;

	}

}*/













void CBackTimeSolve::checkResult( string f_name)

{

	ofstream out(f_name.c_str());

	out << fixed;

	for (int k = K-1; k >= 0; k--)

	{

		for (int i = 0; i < N; i++)

		{

			out << setprecision (5) << "t = " << T - k*tau << "   x = " << stX + i*h << "   u(x,y) = " << pTable[k][i] << endl;

		}

		out << endl << endl;

	}

}



void CBackTimeSolve::setTask(double _a, double _b, double _c, 

						   double _alpha, double _beta, double _gamma, double _delta, 

						    double _stX, double _endX, double _Tr, int _N, int _K, 

							double *_right_border, double *_analitic)

{

	clearArray();

	a = _a;

	b = _b;

	c = _c;

	alpha = _alpha;

	beta = _beta;

	gamma = _gamma;

	delta = _delta;

	l = _endX;

	T = _Tr;

	stX = _stX;

	N = _N + 1;

	K = _K + 1;

	h = _endX/_N;

	tau = _Tr/_K;

	pTable = createTable();

	right_border = _right_border;

	analytic = _analitic;

}



double CBackTimeSolve::r(double x, double t)

{

	    return exp(-1*a*t)*sin(x);

	//return sin(t) * cos(x);

	//return exp(- a * a * t) * cos(x + b * t);

	//return exp((c - a * a) * t) * sin(x + b * t);

}



double ** CBackTimeSolve::createTable()

{

	double ** new_array = new double*[K];

	for (int i =0; i < K ; i++)

	{

		new_array[i] = new double[N];

	}

	return new_array;

}



void CBackTimeSolve::UseSweepMethod(double *d, double *A, double *B, double *C, int N, double *res)

{

	double* P = new double[N];

	double* Q = new double[N];

	double a,b,c;	

	b = B[0];

	c = C[0];

	P[0] = -c / b;

	Q[0] = d[0] / b;



	for (int i=1; i < N-1; i++)

	{

		a = A[i];

		b = B[i];

		c = C[i];

		P[i] = -c / (b + a*P[i-1]);

		Q[i] = (d[i] - a*Q[i-1]) / (b + a*P[i-1]);

	}



	a = A[N-1];

	b = B[N-1];

	P[N-1] = 0;

	Q[N-1] = (d[N-1] - a * Q[N-2]) / (b + a * P[N-2]);



	res[N-1] = Q[N-1];					//SweapMethod

	for (int i=N-2; i >= 0; i--)

	{

		res[i] = P[i]*res[i+1] + Q[i];

	}

}







double CBackTimeSolve:: psiT(double x, double t)

{

	return 0.0;

}



// ????????????

//double CBackTimeSolve:: psiT(double x, double t)

//{

//	return sin (x);

//}







double CBackTimeSolve:: f(double x, double t)

{

	return 0.0;

}



double CBackTimeSolve:: phi0(double x, double t)

{

	return 0.0;

}





double CBackTimeSolve:: phiL(double x, int k)

{ 

	return (2.0 * ( right_border[k] - analytic[k])) ; 

	// ?? ??????? ?????? ?????? ??? ????????? ?? ?????? ??????? ??????? 

	// ????????????? ??????? ????????(??????????????) ??????  

}



// ????????????

//double CBackTimeSolve:: phiL(double x, double t)

//{

//	//return (2.0 * ( border[k] - tStar[k])); 

//	return -1*exp(-1*a*t);

//	//return -sin(t);

//	//return (  exp(-a * a * t) * (cos(b * t) + sin (b*t))  );

//	//return (  -exp((c - a*a)*t) * (cos (b * t) + sin (b * t))  );

//}







/*

void CBackTimeSolve::setBorder(double * _right_border)

{

	right_border = _right_border;

}

*/



void CBackTimeSolve::getSolve()

{





	double * A = new double[N];

	double * B = new double[N];

	double * C = new double[N];

	double * d = new double[N];



	double temp_a, temp_b, temp_c;

	temp_a = (a/h/h - b/h/2.0)*tau;

	temp_b = -(2*a*tau/h/h + 1 - c*tau);

	temp_c = (a/h/h + b/2/h)*tau;



	A[0] = 0;

	B[0]= ((-3 + temp_a/temp_c )*alpha + 2* beta*h);

	C[0]= (4 + temp_b/temp_c)*alpha;



	A[N-1] = - (4 + temp_b/temp_a)*gamma;

	B[N-1] = (3 - temp_c/temp_a)*gamma + 2*h*delta;

	C[N-1] = 0;



	for (int i = 1; i < N-1; i++)

	{

		A[i] = temp_a;

		B[i] = temp_b;

		C[i] = temp_c;

	}



	for (int i = 0; i < N; i++)

	{

		pTable[K-1][i] = psiT(stX + i*h, 0 ); //psiT(i ,0);

	}



	for (int k = K-2; k >= 0; k-- )

	{

		for (int i = 1; i < N; i++ )

		{



			d[i] = -pTable[k+1][i] - f(stX + i* h, ( k*tau) )*tau;



		}



		d[0] = phi0(0, ( k*tau) )*2*h + ( d[1] + f(0, (k*tau)) )*alpha/temp_c ;

		d[N-1] = phiL(l, k )*2*h - (d[N-2] + f(l, k*tau) )*gamma/temp_a;

		// ??? ? ?????? ???? ??????? k- ?? ???????? ?? ???????.

		

		//d[N-1] = phiL(l,( k* tau) )*2*h - (d[N-2] + f(l,( k*tau)) )*gamma/temp_a;      ????????????



		UseSweepMethod(d, A, B, C, N, pTable[k]); 

	}

}