#include "parabolicVector_solver.h"



CParabolicVector_solver::CParabolicVector_solver()

{

	q_arr = NULL;

}



CParabolicVector_solver::~CParabolicVector_solver()

{

	if (pTable != NULL)

	{

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

		{

			delete [] pTable[i];

		}

		delete [] pTable;

		pTable = NULL;

	}

	delete[] q_arr;

	q_arr = NULL;

}



double **CParabolicVector_solver:: getTable()

{

	return pTable;

}



int CParabolicVector_solver::getK()

{

	return K;

}



int CParabolicVector_solver::getN()

{

	return N;

}



double * CParabolicVector_solver::getLeftBorder()

{

	double *temp_arr = new double[K];

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

	{

		temp_arr[k] = pTable[k][0];

	}

	return temp_arr;

}



double * CParabolicVector_solver::getRightBorder()

{

	double *temp_arr = new double[K];

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

	{

		temp_arr[k] = pTable[k][N-1];

	}

	return temp_arr;

}



void CParabolicVector_solver::setGrid(int _N, int _K, double _xb, double _xe, double _maxT)

{

	stepT = _maxT / _K;

	stepX = (_xe - _xb) / _N;



	if (pTable != NULL)

	{

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

		{

			delete [] pTable[i];

		}

		delete [] pTable;

		pTable = NULL;

	}

	if (q_arr != NULL)

	{

		delete [] q_arr;

		q_arr = NULL;

	}





	N = _N + 1;

	K = _K + 1;

	xb = _xb;

	xe = _xe;

	maxT = _maxT;



//



	pTable = new double* [K];

	q_arr = new double[K];

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

	{

		pTable[i] = new double[N];

		q_arr[i] = 0.0;

	}

	

}



void CParabolicVector_solver::setTask(	double _a, double _b, double _c, double(*_f)(double, double),

							   double _alpha1, double _betta1, double _alpha2, double _betta2,

							   double * _left_border, double (*_phi2)(double, double), 

							   double (*_psi)(double, double)) // Only after setGrid;

{

	a = _a;

	b = _b;

	c = _c;

	f = _f;

	alpha1 = _alpha1;

	alpha2 = _alpha2;

	betta1 = _betta1;

	betta2 = _betta2;

	phi2 = _phi2;

	psi = _psi;

	setBorder(_left_border);

}



void CParabolicVector_solver::solve_implicit()

{

	double sigma = (a * a * stepT) / (stepX * stepX);



	//????????? ?????? ?????? ?? ???????	psi(x, 0) = sin(x)

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

	{

		pTable[0][i] = psi(xb + i * stepX, 0.0);

	}



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

	double* A = new double[N];

	double* B = new double[N];

	double* C = new double[N];

	double* D = new double[N];

	double* P = new double[N];

	double* Q = new double[N];

	double* x = new double[N];



	double csigma = a * a * stepT / (stepX * stepX); // only constant, don't want to copy-paste

	double csigma1 = 2.0 * a * a * alpha1 / (2.0 * a * a * stepX - b * stepX * stepX);	// only constant, don't want to copy-paste

	double csigma2 = 2.0 * a * a * alpha2 / (2.0 * a * a * stepX + b * stepX * stepX);	// only constant, don't want to copy-paste



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

	{

		A[0] = 0.0;

		B[0] = csigma1 * ( ((c * stepX * stepX) / (2.0 * a * a)) - 1.0 - ((stepX * stepX) / (2.0 * a * a * stepT))) + betta1;

		C[0] = csigma1;

		D[0] = phi0(xb, i ) - pTable[i-1][0] * (csigma1 * ((stepX * stepX) / (2.0 * a * a * stepT))) - csigma1 * ((stepX * stepX) / (2.0 * a * a)) * f(xb, i * stepT);

	// ??? ? ?????? ???? - ??????????? i -?? ??????? ?? ???????.

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

		{

			A[j] = csigma - ((b * stepT) / (2.0 * stepX));

			B[j] = (-2.0 * csigma) - 1.0 + c * stepT;

			C[j] = csigma + ((b * stepT) / (2.0 * stepX));

			D[j] = -pTable[i-1][j] - f(xb + j * stepX, i * stepT) * stepT;



			if (fabs(B[i]) < fabs(A[i]) + fabs(C[i]))

			{

				printf("error in coeff\n");

			}

		}



		A[N - 1] = -csigma2;

		B[N - 1] = csigma2 * (1.0 + ((stepX * stepX) / (2.0 * a * a * stepT)) - ((c * stepX * stepX) / (2.0 * a * a))) + betta2;

		C[N - 1] = 0.0;

		D[N - 1] = phi2(xe, i * stepT) + pTable[i-1][N - 1] * (csigma2 * stepX * stepX / (2.0 * a * a * stepT)) + f(xe, i * stepT) * csigma2 * stepX * stepX / (2.0 * a * a);



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

		double aa, bb, cc;



		bb = B[0];

		cc = C[0];



		P[0] = -cc / bb;

		Q[0] = D[0] / bb;



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

		{

			aa = A[k];

			bb = B[k];

			cc = C[k];

			P[k] = -cc / (bb + aa*P[k-1]);

			Q[k] = (D[k] - aa*Q[k-1]) / (bb + aa*P[k-1]);

		}



		aa = A[N - 1];

		bb = B[N - 1];

		P[N - 1] = 0;

		Q[N - 1] = (D[N - 1] - aa * Q[N - 2]) / (bb + aa * P[N - 2]);



		x[N - 1] = Q[N - 1];

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

		{

			x[k] = P[k] * x[k+1] + Q[k];

		}

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



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

			pTable[i][j] = x[j];

	}



	delete [] A;

	delete [] B;

	delete [] C;

	delete [] D;

	delete [] P;

	delete [] Q;

	delete [] x;

}



double CParabolicVector_solver::phi0(double x, int i){

	return -q_arr[i];

}



void CParabolicVector_solver::setBorder(double *_arr){

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

		q_arr[i] = _arr[i];

	}

}