#include "HeatComp.h"

using namespace arma;

HeatComp::HeatComp(double dr0, double dz0, double pdphi, double pR, double pL, double pTaumax, double pT_in, char *const_file, double pvz, double pT_inf,double pdtau,double htc)
{
	SetConstM(const_file);

	R = pR;
	L = pL;
	Taumax = pTaumax;
	T_in = pT_in;
	vz = pvz;
	T_inf = pT_inf;
	p=0;

	SetDisc(dr0,dz0,pdphi);
	SetDtau(pdtau);
	SetCoolM(htc);
	
	CrystallizerLength=1;
	htc_k=2000;
	T_k=80;

	TertiaryLength=5;


	T.set_size(ndr,ndphi,ndz);

	cube T1;
	char fname[40]="HeatCompOut/T_0.cub";
	if (std::ifstream(fname).good()) //if there is former computed matrix
	{
		T1.load(fname,arma_ascii);
		if ((T.n_rows==T1.n_rows) && (T.n_cols==T1.n_cols) && (T.n_slices==T1.n_slices))
		{
			T=T1;
			T1.reset();
		}
		else T.fill(T_in);
	}
	else T.fill(T_in);

	for (int k=0;k<T.n_slices;k++)
	middle(k)=getConstT(T(0,0,k),'h');

}


HeatComp::HeatComp(double dr0, double dz0, double dphi, arma::cube T_old, char *const_file, double pvz, double pT_inf, double pdtau,double htc,double pTertiaryLength, double htcl,int n_z)
{
	//SetConstM(const_file);
	ConstM.load("HeatCompOut/ConstM.mat"); //this is workaround avoiding memory adressing issues, downside is, we need to generate this matrix by using first constructor

	Taumax = 0;
	T_in = T_old(0,0,0);
	vz = pvz;
	T_inf = pT_inf;
	p=0;

	CrystallizerLength=1; //NEEDS TO BE CHANGED
	htc_k=2000;
	T_k=80;

	TertiaryLength=pTertiaryLength;
	if (pTertiaryLength<0)
		TertiaryLength=0;

	SetDisc(dr0,dphi,dz0,T_old);
	SetDtau(pdtau);
	SetCoolM(htc,htcl,n_z);

	T=T_old;
	//T.fill(T_in); //ONLY FOR TESTING

	R = dr.n_rows*dr0; //ONLY constant!!!!
	L = dz.n_rows*dz0;
}


void HeatComp::SetDisc(double dr0, double dz0, double phi)
{
	//SetSqR(dr0,R); //dont need it ATM
	ndr=ceil(R/dr0);
	dr.set_size(ndr);
	dr.fill(dr0);

	dphi=phi;
	ndphi=ceil(2*datum::pi/dphi);

	ndz=ceil(L/dz0);
	dz.set_size(ndz);
	dz.fill(dz0);

	middle.set_size(ndz);

}

void HeatComp::SetDisc(double pdr0, double pdphi, double pdz0, cube &T_old)
{
	//SetSqR(dr0,R); //dont need it ATM
	ndr=T_old.n_rows;
	dr.set_size(ndr);
	dr.fill(pdr0);

	ndphi=T_old.n_cols;
	dphi=2*datum::pi/ndphi;


	ndz=T_old.n_slices;
	dz.set_size(ndz);
	dz.fill(pdz0);

	middle.set_size(ndz);

}

double HeatComp::getR(int ind)
{
	double out=dr(0);
	for(int j = 0; j < ind; j++)
		out+=dr(j);
	return out;
}


void HeatComp::SetSqR(double dr0 ,double pR)
{
	double x;
	x=(1.0-(dr0/pR))*(1.0-(dr0/pR));

	ndr=ceil(1.0/(1.0-x));
	dr.set_size(ndr);
	for(int i = 0; i < ndr; i++)
		dr(i)=(sqrt((double)(i+1)/ndr)-sqrt((double)(i)/ndr))*R; //prepsano o jedna
}

void HeatComp::SetConstM(char * filename)
{

	ConstM0.load(filename);
	ConstM.set_size(ConstM0(ConstM0.n_rows-1),3);

	double c=0;

	for (int i = 0; i < ConstM0.n_rows; i++) //set exact values
	{
		c=floor(ConstM0(i,0)+0.5); //missing function round();
		ConstM(c-1,0)=ConstM0(i,0);
		ConstM(c-1,1)=1000*ConstM0(i,1); //recaculate to basic units
		ConstM(c-1,2)=ConstM0(i,2);
	}

	int j=0,br=0,i_u=0,i_l=0;

	for (int i = 25; i < ConstM.n_rows - 1; i++) //interpolate
	{
		if (ConstM(i,0)<0.1)
		{
			ConstM(i,0)=i+1;
			i_u=i-1;
			while (br<1)
			{
				j++;

				if (ConstM(i+j,0)>0)
				{
					i_l=i+j;
					br=10;
				}
			}
			br=j=0;

			for (int l = 1; l < 3; l++)
			{
				ConstM(i,l) = ConstM(i_u,l) + ((ConstM(i_l,l) - ConstM(i_u,l))*(i - i_u)) / (i_l - i_u);
			}
		}
	}
	i_u=i_l=0;

	ConstM.save("HeatCompOut/ConstM.mat",raw_ascii);
	ConstM0.save("HeatCompOut/ConstMorig.mat",raw_ascii);
}

double HeatComp::getConstT(double &T, char swtch)
{
	vec cnst(3);
	cnst(1)=ConstM(floor(T+0.5)-1,1);
	cnst(2)=ConstM(floor(T+0.5)-1,2);
	switch (swtch) {
	case 'h' :
		return cnst(1);

	case 'k' :
		return cnst(2);
	default :
		exit (EXIT_FAILURE);
	}
}

double HeatComp::getConstH(double &H, char swtch) //accuracy +- 1 degree °C
{
	vec cnst(3);
	uvec c=find(ConstM.col(1) > H,1);
//	uvec c=find(ConstM0.col(1) <= (H/1000),1,"last"); //UPRAvit?
//	int a=c(0);
//	for (int col=0; col<3;col++)
//		{
//		      cnst(col)=( (ConstM0(a+1,col)-ConstM0(c(0),col))*( (H/1000)-ConstM0(c(0),1) ) ) / ( ConstM0(c(0)+1,1)-ConstM0(c(0),1) ) + ConstM0(c(0),col);
//		}

	cnst(0)=ConstM(c(0)-1,0);
	cnst(1)=ConstM(c(0)-1,1);
	cnst(2)=ConstM(c(0)-1,2);

	switch (swtch) {
	case 't' :
		return cnst(0);
	case 'h' :
		return cnst(1);
	case 'k' :
		return cnst(2);
	default :
		exit (EXIT_FAILURE);
	}
}

cube HeatComp::TempToEn(cube &T)
{
	cube H(T.n_rows,T.n_cols,T.n_slices);

	for (int k = 0; k < T.n_slices; k++)
	{
		for (int j = 0; j < T.n_cols; j++)
		{
			for (int i = 0; i < T.n_rows; i++)
				H(i,j,k)=getConstT(T(i,j,k),'h');
		}
	}

	return H;
}

cube HeatComp::EnToTemp(cube &H)
{
	cube pT;
	pT.set_size(H.n_rows, H.n_cols, H.n_slices);

	for (int k = 0; k < H.n_slices; k++)
	{
		for (int j = 0; j < H.n_cols; j++)
		{
			for (int i = 0; i < H.n_rows; i++)
				pT(i,j,k)=getConstH(H(i,j,k),'t');
		}
	}

	return pT;

}

void HeatComp::SetCoolM(double htc) //TBD
{
	CoolM.set_size(ndphi,ndz);
	CoolM.fill(htc);
//	for (int i=0;i<50;i++)
//	CoolM(span::all,span(i,i+1)).fill(500);

}

void HeatComp::SetCoolM(double htc0, double htcl, int n_z)
{
	CoolM.set_size(ndphi,ndz);
	CoolM.fill(htc0);
	CoolM(span::all,span(0,n_z)).fill(htcl);

}

void HeatComp::SetDtau(double pdtau) //TBD
{
	(*this).dtau=pdtau;
}

HeatComp::~HeatComp(void) //vz, T_inf, dtau
{
	erase();
}

void HeatComp::erase(void)
{
	T.reset();
		ConstM.reset(); CoolM.reset();
		dphi=0;
		R=0;
		L=0;
		vz=0;
		Taumax=0;
		T_inf=0;
		T_in=0;
		TertiaryLength=0;
		CrystallizerLength=0;
		htc_k=0;
		T_k=0;
		dr.reset();dz.reset();
		ndr=ndphi=ndz=0;
		middle.reset();
}

cube HeatComp::Compute(int pmax, char swtch) //predelat i,j,k na k,i,j!
{
	double epsilon,htc_n,htc_r,QdrOUT,Qdr,Qdphi,Qdz,Tout;
	int size_m=T.n_rows;
	int size_n=T.n_cols;
	int size_o=T.n_slices;
	int Qdr_mid=0;
	double c=0;
//	int pmax=7000;//(int)(Taumax/dtau);
//	int delta=1000;
//	double c,c1,c2,c3,c4,c5;
	cube H;
	H.set_size(size_m,size_n,size_o);
	
	cube T_old=T;
//	cube dT;
	mat TempProfile(size_m,size_o);
	vec middle_old=middle;
	cube H_old = TempToEn(T);

	bool use_middle=false;
	char fnameT[]="HeatCompOut/T_0_final.cub";
	char fnameH[]="HeatCompOut/H_0.cub";
	char fnameOctave[]="TempProfile_final.mat";
	p=0;

	while (p<pmax)
	{
		for(int k = 0; k < size_o; k++) //z coordinate
		{
			for(int j = 0; j < size_n; j++) // angle phi
			{
				for(int i = 0; i < size_m; i++) // r coordinate
				{

					if ((i==ndr-1) && (j>0) && (j<ndphi-1) && (k>0) && (k<ndz-1)) //BC condition on the surface
					{
						//epsilon=0.85;

						c=42.68-0.0268*(T_old(ndr-1,j,k)+273.15);
						epsilon=0.85/pow(exp(1+c),0.0115);

						htc_n = 0.84*pow(T_old(ndr-1,j,k) - T_inf,1.0/3) + CoolM(j,k);
						if (k>ndz-TertiaryLength/((dz(k)+dz(k+1))/2))
							htc_n = 0.84*pow(T_old(ndr-1,j,k) - T_inf,1.0/3);
						//cooling with jets modelled as aditional htc
						htc_r = epsilon*sigma*(pow(T_old(ndr-1,j,k) + 273.15,2)+ pow(T_inf + 273.15,2))*((T_old(ndr-1,j,k) + 273.15) + (T_inf + 273.15));
						QdrOUT = (htc_r + htc_n)*(T_old(ndr-1,j,k) - T_inf);
						//QdrOUT=0;
						if (k<CrystallizerLength/((dz(k)+dz(k+1))/2))
							QdrOUT=htc_k*(T_old(ndr-1,j,k)-T_k);
					//	QdrOUT=8000;

						// Tout: computing of surface temperature
						//double acha=;
//						cout <<dr(ndr-1)<< endl;
//						cout <<getR(ndr-1)<< endl;
//						cout <<QdrOUT<< endl;
//
//						cout << getConstT(T_old(ndr-1,j,k),'k') << endl;
//						cout <<T(ndr-2,j,k)<< endl;

						Tout = -(2 * dr(ndr-1)*getR(ndr-1) * QdrOUT)/getConstT(T_old(ndr-1,j,k), 'k') + T(ndr-2,j,k);
//						if (p % (int)Tout==500)
//							cout <<Tout<<'\n';
						Qdr= 2*((Tout-T(ndr-1,j,k))/dr(ndr-1) - (T(ndr-1,j,k)-T(ndr-2,j,k))/dr(ndr-2))/((dr(ndr-1)+dr(ndr-2)));
						Qdphi=(T(ndr-1,j+1,k)-2*T(ndr-1,j,k)+T(ndr-1,j-1,k))/((getR(ndr-1)*dphi)*(getR(ndr-1)*dphi));
						Qdz= 2*((T(ndr-1,j,k+1)-T(ndr-1,j,k))/dz(k) - (T(ndr-1,j,k)-T(ndr-1,j,k-1))/dz(k-1))/(dz(k)+dz(k-1));

						H(ndr-1,j,k)=H_old(ndr-1,j,k)+dtau*(getConstT(T_old(ndr-1,j,k),'k'))*(Qdr+Qdphi+Qdz)+vz*dtau*((H_old(ndr-1,j,k-1)-H_old(ndr-1,j,k))/dz(k-1));
					}

					if ((i>0) && (i<ndr-1) && (j>0) && (j<ndphi-1) && (k>0) && (k<ndz-1)) //Inner computation
					{
						
						Qdr= 2*((T(i+1,j,k)-T(i,j,k))/dr(i+1) - (T(i,j,k)-T(i-1,j,k))/dr(i))/((dr(i+1)+dr(i)));
						Qdphi=(T(i,j-1,k)-2*T(i,j,k)+T(i,j+1,k))/((getR(i)*dphi)*(getR(i)*dphi));
						Qdz= 2*((T(i,j,k+1)-T(i,j,k))/dz(k) - (T(i,j,k)-T(i,j,k-1))/dz(k-1))/(dz(k)+dz(k-1));

						H(i,j,k)=H_old(i,j,k)+dtau*getConstT(T_old(i,j,k),'k')*(Qdr+Qdphi+Qdz)+vz*dtau*((H_old(i,j,k-1)-H_old(i,j,k))/dz(k-1));
					}

					//BC conditions

					switch (swtch)
					{
					case 'y':
						H(i,j,0) = getConstT(T_in, 'h'); //temp in crystallizer
						break;
					case 'n':
						H(i,j,0) = H(i,j,1); //zero flux
						break;
					}
					middle(0)=getConstT(T_in, 'h');
					H(i,0,k)=H(i,1,k);
					H(i,ndphi-1,k)=H(i,ndphi-2,k);
					H(i,j,ndz-1)=H(i,j,ndz-2);
					middle(ndz-1)=middle(ndz-2);

					if ((i==0) && (use_middle) && (j>0) && (j<ndphi-1) && (k>0) && (k<ndz-1))//to compute elements next to the axis
					{
						Qdr= 2*((T(i+1,j,k)-T(i,j,k))/dr(i) - (T(i,j,k)-getConstH(middle(k),'t'))/dr(0))/((dr(0)+dr(0)));
						Qdphi=(T(i,j-1,k)-2*T(i,j,k)+T(i,j+1,k))/((getR(i)*dphi)*(getR(i)*dphi));
						Qdz= 2*((T(i,j,k+1)-T(i,j,k))/dz(k) - (T(i,j,k)-T(i,j,k-1))/dz(k-1))/(dz(k)+dz(k-1));

						H(i,j,k)=H_old(i,j,k)+dtau*getConstT(T_old(i,j,k),'k')*(Qdr+Qdphi+Qdz)+vz*dtau*((H_old(i,j,k-1)-H_old(i,j,k))/dz(k-1));

					}

					if (!use_middle)
						H(0,j,k)=H(1,j,k);

				}
				if ((use_middle) && (k>0) && (k<ndz-1))
				{
					Qdr_mid+= 2*((T(0,j,k)-getConstH(middle(k),'t'))/dr(0) - (getConstH(middle(k),'t')-T(1,j,k))/dr(0))/((dr(0)+dr(0)));
					Qdphi=0;
					Qdz= 2*((getConstH(middle(k+1),'t')-getConstH(middle(k),'t'))/dz(k) - (getConstH(middle(k),'t')-getConstH(middle(k-1),'t'))/dz(k-1))/(dz(k)+dz(k-1));
				}
			}

			if ((use_middle) && (k>0) && (k<ndz-1))
			{
				middle(k)=middle_old(k)+dtau*getConstH(middle_old(k),'k')*(Qdr_mid+Qdphi+Qdz)+vz*dtau*((middle_old(k-1)-middle_old(k))/dz(k-1));
//				c=middle(k);
				Qdr_mid=0;
			}
		}
	//	H.save("H_0.3.cub",arma_ascii);
	//	T.save("T_1.cub",arma_ascii);
	//	H.save("H_err.cub",arma_ascii);
		T=EnToTemp(H); //getConstT a EnToEnt NAPROSTO kriticke v rychlosti vypcotu
		middle_old=middle;

		if (p % 500 == 0)
		{
			cout << p << "          out of " << pmax << "\n";
//			H.save(fnameH,arma_ascii);
//			T.save(fnameT,arma_ascii);
//			cout <<a/(ndphi-2) <<'\n';
		}
		p++;
//		dT=abs(T-T_old);
//		delta=dT.slice(ndz-1).max();

		H_old=H;
		T_old=T;
	}
	
	H.save(fnameH,arma_ascii);
	T.save(fnameT,arma_ascii);
//	T.save("HeatCompOut/Thdf5.cub",hdf5_binary);
	for ( int k=0; k<T.n_slices; k++)
	{
			TempProfile(span::all,k)=T.slice(k)(span::all,1);
	}
	char str0[50];
	sprintf(str0,"HeatCompOut/%s",fnameOctave);

	TempProfile.save(str0,raw_ascii);
	char str[100];
	sprintf(str,"cd \HeatCompOut; octave -q --eval \"visualiseF('%s')\";exit",fnameOctave);
	system(str);

	return T;
}

cube HeatComp::ComputeOpt(int ind,int pmax, cube T_L) //predelat i,j,k na k,i,j!
{
	double epsilon,htc_n,htc_r,QdrOUT,Qdr,Qdphi,Qdz,Tout;
	int size_m=T.n_rows;
	int size_n=T.n_cols;
	int size_o=T.n_slices;
	int Qdr_mid=0;

//	int pmax=7000;//(int)(Taumax/dtau);
//	int delta=1000;
//	double c,c1,c2,c3,c4,c5;
	cube H;
	H.set_size(size_m,size_n,size_o);

	cube T_old=T;
//	cube dT;
	mat TempProfile(size_m,size_o);
	vec middle_old=middle;
	cube H_old = TempToEn(T);

	bool use_middle=false;
	char fnameT[]="HeatCompOut/T_0.cub";
	char fnameH[]="HeatCompOut/H_0.cub";
	p=0;
	int c;



	while (p<pmax)
	{
		for(int k = 0; k < size_o; k++) //z coordinate
		{
			for(int j = 0; j < size_n; j++) // angle phi
			{
				for(int i = 0; i < size_m; i++) // r coordinate
				{

					if ((i==ndr-1) && (j>0) && (j<ndphi-1) && (k>0) && (k<ndz-1)) //BC condition on the surface
					{
						//epsilon=0.85;
						c=42.68-0.0268*(T_old(ndr-1,j,k)+273.15);
						epsilon=0.85/pow(exp(1+c),0.0115);

						htc_n = 0.84*pow(T_old(ndr-1,j,k) - T_inf,1.0/3) + CoolM(j,k); //cooling with jets modelled as aditional convection htc
						if (k>ndz-TertiaryLength/((dz(k)+dz(k+1))/2))//TERTIARY COOLING NEEDS WORK
							htc_n = 0.84*pow(T_old(ndr-1,j,k) - T_inf,1.0/3);

						htc_r = epsilon*sigma*(pow(T_old(ndr-1,j,k) + 273.15,2) + pow(T_inf + 273.15,2))*((T_old(ndr-1,j,k) + 273.15) + (T_inf + 273.15));

						QdrOUT = (htc_r + htc_n)*(T_old(ndr-1,j,k) - T_inf);
						if ((k<CrystallizerLength/((dz(k)+dz(k+1))/2)) && (ind==0) ) //part with crystallizer
							QdrOUT=htc_k*(T_old(ndr-1,j,k)-T_k);
					//	QdrOUT=8000;

						// Tout: computing of surface temperature
						//double acha=;
//						cout <<dr(ndr-1)<< endl;
//						cout <<getR(ndr-1)<< endl;
//						cout <<QdrOUT<< endl;
//
//						cout << getConstT(T_old(ndr-1,j,k),'k') << endl;
//						cout <<T(ndr-2,j,k)<< endl;

						Tout = -(2 * dr(ndr-1)*getR(ndr-1) * QdrOUT)/getConstT(T_old(ndr-1,j,k), 'k') + T(ndr-2,j,k); //as T(ndr,j,k)
//						if (p % (int)Tout==500)
//							cout <<Tout<<'\n';
						Qdr= 2*((Tout-T(ndr-1,j,k))/dr(ndr-1) - (T(ndr-1,j,k)-T(ndr-2,j,k))/dr(ndr-2))/((dr(ndr-1)+dr(ndr-2)));
						Qdphi=(T(ndr-1,j+1,k)-2*T(ndr-1,j,k)+T(ndr-1,j-1,k))/((getR(ndr-1)*dphi)*(getR(ndr-1)*dphi));
						Qdz= 2*((T(ndr-1,j,k+1)-T(ndr-1,j,k))/dz(k) - (T(ndr-1,j,k)-T(ndr-1,j,k-1))/dz(k-1))/(dz(k)+dz(k-1));

						H(ndr-1,j,k)=H_old(ndr-1,j,k)+dtau*(getConstT(T_old(ndr-1,j,k),'k'))*(Qdr+Qdphi+Qdz)+vz*dtau*((H_old(ndr-1,j,k-1)-H_old(ndr-1,j,k))/dz(k-1));
					}

					if ((i>0) && (i<ndr-1) && (j>0) && (j<ndphi-1) && (k>0) && (k<ndz-1)) //Inner computation
					{

						Qdr= 2*((T(i+1,j,k)-T(i,j,k))/dr(i+1) - (T(i,j,k)-T(i-1,j,k))/dr(i))/((dr(i+1)+dr(i)));
						Qdphi=(T(i,j-1,k)-2*T(i,j,k)+T(i,j+1,k))/((getR(i)*dphi)*(getR(i)*dphi));
						Qdz= 2*((T(i,j,k+1)-T(i,j,k))/dz(k) - (T(i,j,k)-T(i,j,k-1))/dz(k-1))/(dz(k)+dz(k-1));

						H(i,j,k)=H_old(i,j,k)+dtau*getConstT(T_old(i,j,k),'k')*(Qdr+Qdphi+Qdz)+vz*dtau*((H_old(i,j,k-1)-H_old(i,j,k))/dz(k-1));
					}

					//BC conditions

					switch (ind)
					{
					case 0:
						H(i,j,0) = getConstT(T_in, 'h'); //temp in crystallizer
						break;
					default:
						H(i,j,0) = getConstT(T_L(i,j,0),'h');
						//H(i,j,0) = H(i,j,1);
//						if ((i>0) && (i<ndr-1) && (j>0) && (j<ndphi-1) && (k==0)) //Inner computation for 0 cut using -1 element (T_L)
//						{
//							Qdr= 2*((T(i+1,j,0)-T(i,j,0))/dr(i+1) - (T(i,j,0)-T(i-1,j,0))/dr(i))/((dr(i+1)+dr(i)));
//							Qdphi=(T(i,j-1,0)-2*T(i,j,0)+T(i,j+1,0))/((getR(i)*dphi)*(getR(i)*dphi));
//							Qdz= 2*((T(i,j,0+1)-T(i,j,0))/dz(0) - (T(i,j,0)-T_L(i,j,0))/dz(0))/(dz(0)+dz(0));
//							H(i,j,0)=H_old(i,j,0)+dtau*getConstT(T_old(i,j,0),'k')*(Qdr+Qdphi+Qdz)+vz*dtau*((getConstT(T_L(i,j,0),'h')-H_old(i,j,k))/dz(k));
//						}

						break;
					}
					middle(0)=getConstT(T_in, 'h');
					H(i,0,k)=H(i,1,k);
					H(i,ndphi-1,k)=H(i,ndphi-2,k);
					H(i,j,ndz-1)=H(i,j,ndz-2);
					middle(ndz-1)=middle(ndz-2);

					if ((i==0) && (use_middle) && (j>0) && (j<ndphi-1) && (k>0) && (k<ndz-1))//to compute elements next to the axis
					{
						Qdr= 2*((T(i+1,j,k)-T(i,j,k))/dr(i) - (T(i,j,k)-getConstH(middle(k),'t'))/dr(0))/((dr(0)+dr(0)));
						Qdphi=(T(i,j-1,k)-2*T(i,j,k)+T(i,j+1,k))/((getR(i)*dphi)*(getR(i)*dphi));
						Qdz= 2*((T(i,j,k+1)-T(i,j,k))/dz(k) - (T(i,j,k)-T(i,j,k-1))/dz(k-1))/(dz(k)+dz(k-1));

						H(i,j,k)=H_old(i,j,k)+dtau*getConstT(T_old(i,j,k),'k')*(Qdr+Qdphi+Qdz)+vz*dtau*((H_old(i,j,k-1)-H_old(i,j,k))/dz(k-1));

					}

					if (!use_middle)
						H(0,j,k)=H(1,j,k);

				}
				if ((use_middle) && (k>0) && (k<ndz-1))
				{
					Qdr_mid+= 2*((T(0,j,k)-getConstH(middle(k),'t'))/dr(0) - (getConstH(middle(k),'t')-T(1,j,k))/dr(0))/((dr(0)+dr(0)));
					Qdphi=0;
					Qdz= 2*((getConstH(middle(k+1),'t')-getConstH(middle(k),'t'))/dz(k) - (getConstH(middle(k),'t')-getConstH(middle(k-1),'t'))/dz(k-1))/(dz(k)+dz(k-1));
				}
			}

			if ((use_middle) && (k>0) && (k<ndz-1))
			{
				middle(k)=middle_old(k)+dtau*getConstH(middle_old(k),'k')*(Qdr_mid+Qdphi+Qdz)+vz*dtau*((middle_old(k-1)-middle_old(k))/dz(k-1));
//				c=middle(k);
//				cout <<c <<'\n';
				Qdr_mid=0;
			}
		}

	//	H.save("H_0.3.cub",arma_ascii);
	//	T.save("T_1.cub",arma_ascii);
	//	H.save("H_err.cub",arma_ascii);
		T=EnToTemp(H); //getConstT a EnToEnt NAPROSTO kriticke v rychlosti vypcotu
		middle_old=middle;

		if (p % 500 == 0)
		{
			//cout << p << "          out of " << pmax << "\n";
//			H.save(fnameH,arma_ascii);
//
		}
		p++;
//		dT=abs(T-T_old);
//		delta=dT.slice(ndz-1).max();

		H_old=H;
		T_old=T;
		//T.save("T.cub",arma_ascii);
	}

//	H.save(fnameH,arma_ascii);
//	T.save(fnameT,arma_ascii);
//	T.save("HeatCompOut/Thdf5.cub",hdf5_binary);
	for ( int k=0; k<T.n_slices; k++)
	{
			TempProfile(span::all,k)=T.slice(k)(span::all,1);
	}

//	TempProfile.save("HeatCompOut/TempProfile_part.mat",raw_ascii);
	//cout << 'Total number of iterations: ' << p <<'\n';
	//cout << 'Real time: '<<(int)(p*dtau/60) <<'min.';
	//fprintf('\n Real time computed: %d min \n',p*(Taumax/dtau)/60);
	return T;
}


