/***************************************************************************
 *   Copyright (C) 2014 by Jordi Farjas					   *
 ***************************************************************************/
// Alternating direction implicit method.

#include <math.h>
#include "Combustion.h"

// iteration to solve the alternating implicit scheme: radial implicit
// Temp_a,alfa_a:	temperatures, K, and degree of transformation of previous iteration
// Temp_n,alfa_n:	new evaluated temperatures, K, and degree of transformation
// A1, A2:			coeficients used in the resolution of the implicit equation

void Sample::CalculaImp_r(double **alfa_n,double **alfa_a,double **Temp_n,double **Temp_a)
{
 int ii,jj;
 double beta,tmp,a_temp,rel_cond,rel_condm1d2,delta,deltap,ipsil,ipsilp,ipsilc,d_iM,d_im;

 rel_cond=Cond_subs/Cond_sample;
 rel_condm1d2=(1.0+rel_cond)/2.0;
 delta=2.0/(1.0+rel_cond/Dif_subs_adim);
 deltap=4.0/(1.0+3.0*rel_cond/Dif_subs_adim);
 ipsil=rel_condm1d2*gamma_2*delta;
 ipsilp=gamma_2*deltap*rel_condm1d2;
 ipsilc=deltap*gamma_2*rel_cond;

  //sample + substrat
 for(jj=1;jj<n_z_sample-1;jj++)
	{
	 ii=n_r_tot-1;
	 if(Boundary_cond & RGT_PERIODIC)
		{	// coefficients A1 and A2 for the right boundary condition dT/dr=0
		 A_1[ii]=1.0;
		 A_2[ii]=0.0;	
		}
	 else if(Boundary_cond & RGT_FLOW_CT)
		{	// coefficients A1 and A2 for the right boundary condition dT/dr=f_const
		 A_1[ii]=1.0;
		 A_2[ii]=f_const;	
		}
	 else 
		{	// coefficients A1 and A2 for the right boundary condition dT/dr=f_const
		 A_1[ii]=0.0;
		 A_2[ii]=T_furnace;	
		}
	 // right vessel wall
	 for(ii=n_r_tot-2;ii>=n_r_sample;ii--) 
		{
		 d_iM=1.0+0.5/ii;
		 d_im=1.0-0.5/ii;
		 tmp=1.0+2.0*gamma_subs_2-gamma_subs_2*d_iM*A_1[ii+1];
		 a_temp=d_iM*A_2[ii+1]+Temp_a[ii][jj+1]+Temp_a[ii][jj-1]-2.0*Temp_a[ii][jj];
		 A_2[ii]=(gamma_subs_2*a_temp+Temp_a[ii][jj])/tmp;
		 A_1[ii]=gamma_subs_2*d_im/tmp;
		}
     if(n_r_subs>0)	// vertical interface
		{
		 ii=n_r_sample-1;
		 d_iM=1.0+0.5/ii;
		 d_im=1.0-0.5/ii;
		 if(Reaction_exists) beta=(alfa_n[ii][jj]-alfa_a[ii][jj])/H_cap_samp_adim/2.0*delta;
		 else beta=0.0;
		 beta+=Temp_a[ii][jj];
		 beta+=rel_cond*delta*gamma_2*d_iM*A_2[ii+1];
	 	 tmp=1.0-delta*gamma_2*rel_cond*d_iM*A_1[ii+1]+delta*gamma_2*(d_im+rel_cond*d_iM);
		 a_temp=gamma_2*delta*rel_condm1d2*(Temp_a[ii][jj+1]+Temp_a[ii][jj-1]-2.0*Temp_a[ii][jj]);
		 A_2[ii]=(a_temp+beta)/tmp;
		 A_1[ii]=delta*d_im*gamma_2/tmp;
		}
	 // sample calculation
	 if(Reaction_exists)
		{
		for(ii=n_r_sample-2;ii>0;ii--) 
			{
			 d_iM=1.0+0.5/ii;
			 d_im=1.0-0.5/ii;
			 beta=(alfa_n[ii][jj]-alfa_a[ii][jj])/H_cap_samp_adim;
			 tmp=1.0+2.0*gamma_2-gamma_2*d_iM*A_1[ii+1];
			 a_temp=d_iM*A_2[ii+1]+Temp_a[ii][jj+1]+Temp_a[ii][jj-1]-2.0*Temp_a[ii][jj];
			 A_2[ii]=(gamma_2*a_temp+beta+Temp_a[ii][jj])/tmp;
			 A_1[ii]=gamma_2*d_im/tmp;
			}
		}
	 else
		{
		 for(ii=n_r_sample-2;ii>0;ii--) 
			{
			 d_iM=1.0+0.5/ii;
			 d_im=1.0-0.5/ii;
			 tmp=1.0+2.0*gamma_2-gamma_2*d_iM*A_1[ii+1];
 			 a_temp=d_iM*A_2[ii+1]+Temp_a[ii][jj+1]+Temp_a[ii][jj-1]-2.0*Temp_a[ii][jj];
			 A_2[ii]=(gamma_2*a_temp+Temp_a[ii][jj])/tmp;
			 A_1[ii]=gamma_2*d_im/tmp;
			}
		}
	// Due to symmetry at the z axe dT/dr=0
	 Temp_n[0][jj]=A_2[1]/(1.0-A_1[1]);	
	 
	 for(ii=1;ii<n_r_tot;ii++) Temp_n[ii][jj]=A_1[ii]*Temp_n[ii-1][jj]+A_2[ii];
	}	// end sample+vessel

 // bottom vessel
 for(jj=n_z_sample;jj<n_z_tot-1;jj++)
	{
	 ii=n_r_tot-1;
	 if(Boundary_cond & RGT_PERIODIC)
		{	// coefficients A1 and A2 for the right boundary condition dT/dr=0
		 A_1[ii]=1.0;
		 A_2[ii]=0.0;	
		}
	 else if(Boundary_cond & RGT_FLOW_CT)
		{	// coefficients A1 and A2 for the right boundary condition dT/dr=f_const
		 A_1[ii]=1.0;
		 A_2[ii]=f_const;	
		}
	 else 
		{	// coefficients A1 and A2 for the right boundary condition T=T_furnace
		 A_1[ii]=0.0;
		 A_2[ii]=T_furnace;	
		}
	 for(ii=n_r_tot-2;ii>0;ii--) 
		{
		 d_iM=1.0+0.5/ii;
		 d_im=1.0-0.5/ii;
		 tmp=1.0+2.0*gamma_subs_2-gamma_subs_2*d_iM*A_1[ii+1];
		 a_temp=d_iM*A_2[ii+1]+Temp_a[ii][jj+1]+Temp_a[ii][jj-1]-2.0*Temp_a[ii][jj];
		 A_2[ii]=(gamma_subs_2*a_temp+Temp_a[ii][jj])/tmp;
		 A_1[ii]=gamma_subs_2*d_im/tmp;
		}
	 	// Due to symmetry at the z axe dT/dr=0
	 Temp_n[0][jj]=A_2[1]/(1.0-A_1[1]);	
	 
	 for(ii=1;ii<n_r_tot;ii++) Temp_n[ii][jj]=A_1[ii]*Temp_n[ii-1][jj]+A_2[ii];
	}	// end bottom vessel

 if(n_z_subs>0)// horizontal interface
	{
	 jj=n_z_sample-1;
	 ii=n_r_tot-1;
	 if(Boundary_cond & RGT_PERIODIC)
		{	// coefficients A1 and A2 for the right boundary condition dT/dr=0
		 A_1[ii]=1.0;
		 A_2[ii]=0.0;	
		}
	 else if(Boundary_cond & RGT_FLOW_CT)
		{	// coefficients A1 and A2 for the right boundary condition dT/dr=f_const
		 A_1[ii]=1.0;
		 A_2[ii]=f_const;	
		}
	 else 
		{	// coefficients A1 and A2 for the right boundary condition T=T_furnace
		 A_1[ii]=0.0;
		 A_2[ii]=T_furnace;	
		}
     // right vessel wall
	 for(ii=n_r_tot-2;ii>=n_r_sample;ii--) 
		{
		 d_iM=1.0+0.5/ii;
		 d_im=1.0-0.5/ii;
		 tmp=1.0+2.0*gamma_subs_2-gamma_subs_2*d_iM*A_1[ii+1];
		 a_temp=d_iM*A_2[ii+1]+Temp_a[ii][jj+1]+Temp_a[ii][jj-1]-2.0*Temp_a[ii][jj];
		 A_2[ii]=(gamma_subs_2*a_temp+Temp_a[ii][jj])/tmp;
		 A_1[ii]=gamma_subs_2*d_im/tmp;
		}

	 if(n_r_subs>0)	// interface at the sample's corner
		{
		 jj=n_z_sample-1;
		 ii=n_r_sample-1;
		 d_iM=1.0+0.5/ii;
		 d_im=1.0-0.5/ii;
		 if(Reaction_exists) beta=(alfa_n[ii][jj]-alfa_a[ii][jj])/H_cap_samp_adim/4.0*deltap;
		 else beta=0.0;
		 tmp=1.0+ipsilp*d_im+ipsilc*d_iM*(1.0-A_1[ii+1]);
		 a_temp=ipsilp*(Temp_a[ii][jj-1]-Temp_a[ii][jj]);
		 a_temp+=ipsilc*(Temp_a[ii][jj+1]-Temp_a[ii][jj]);
		 A_2[ii]=(a_temp+beta+ipsilc*d_iM*A_2[ii+1]+Temp_a[ii][jj])/tmp;
		 A_1[ii]=d_im*ipsilp/tmp;
		}	
	 if(Reaction_exists) // horizontal interface
		{
		 for(ii=n_r_sample-2;ii>0;ii--) 
			{
			 d_iM=1.0+0.5/ii;
			 d_im=1.0-0.5/ii;
			 beta=(alfa_n[ii][jj]-alfa_a[ii][jj])/H_cap_samp_adim/2.0*delta;
			 tmp=1.0-ipsil*(d_iM*A_1[ii+1]-2.0);
			 a_temp=gamma_2*delta*(Temp_a[ii][jj-1]-Temp_a[ii][jj]);
			 a_temp+=gamma_2*delta*rel_cond*(Temp_a[ii][jj+1]-Temp_a[ii][jj]);
			 A_2[ii]=(a_temp+beta+ipsil*d_iM*A_2[ii+1]+Temp_a[ii][jj])/tmp;
			 A_1[ii]=ipsil*d_im/tmp;
			}
		}
	 else
		{
		 for(ii=n_r_sample-2;ii>0;ii--) 
			{
	 	 	 d_iM=1.0+0.5/ii;
			 d_im=1.0-0.5/ii;
			 tmp=1.0-ipsil*(d_iM*A_1[ii+1]-2.0);
 			 a_temp=gamma_2*delta*(Temp_a[ii][jj-1]-Temp_a[ii][jj]);
			 a_temp+=gamma_2*delta*rel_cond*(Temp_a[ii][jj+1]-Temp_a[ii][jj]);
			 A_2[ii]=(a_temp+ipsil*d_iM*A_2[ii+1]+Temp_a[ii][jj])/tmp;
			 A_1[ii]=ipsil*d_im/tmp;
			}
		}
		// Due to symmetry at the z axe dT/dr=0
	 Temp_n[0][jj]=A_2[1]/(1.0-A_1[1]);	

	 for(ii=1;ii<n_r_tot;ii++) Temp_n[ii][jj]=A_1[ii]*Temp_n[ii-1][jj]+A_2[ii];
	}	// end interface

// top boundary conditions
  if(Boundary_cond & TOP_PERIODIC) // top boundary condition dT/dz=0
	 for(ii=0;ii<n_r_tot;ii++)	Temp_n[ii][0]=Temp_n[ii][1];		
  else if(Boundary_cond & TOP_FLOW_CT)	//  top boundary condition dT/dz=f_const
	 for(ii=0;ii<n_r_tot;ii++)	Temp_n[ii][0]=Temp_n[ii][1]-f_const;		
  else //  top boundary condition T = T_furnace
	 for(ii=0;ii<n_r_tot;ii++) Temp_n[ii][0]=T_furnace;

  // bottom boundary conditions
  if(Boundary_cond & BTM_PERIODIC) // bottom boundary condition dT/dz=0
	 for(ii=0;ii<n_r_tot;ii++)	Temp_n[ii][n_z_tot-1]=Temp_n[ii][n_z_tot-2];		
  else if(Boundary_cond & BTM_FLOW_CT)	// bottom boundary condition dT/dz=f_const
	 for(ii=0;ii<n_r_tot;ii++)	Temp_n[ii][n_z_tot-1]=Temp_n[ii][n_z_tot-2]+f_const;		
  else //  bottom boundary condition T= T_furnace
	 for(ii=0;ii<n_r_tot;ii++) Temp_n[ii][n_z_tot-1]=T_furnace;
}

// iteration to solve the alternating implicit scheme: axial implicit
// Temp_a,alfa_a:	temperatures, K, and degree of transformation of previous iteration
// Temp_n,alfa_n:	new evaluated temperatures, K, and degree of transformation
// A1, A2:			coeficients used in the resolution of the implicit equation

void Sample::CalculaImp_z(double **alfa_n,double **alfa_a,double **Temp_n,double **Temp_a)
{
 int ii,jj;
 double beta,tmp,a_temp,rel_cond,rel_condm1d2,delta,deltap,ipsil,ipsilp,ipsilc,d_iM,d_im;

 rel_cond=Cond_subs/Cond_sample;
 rel_condm1d2=(1.0+rel_cond)/2.0;
 delta=2.0/(1.0+rel_cond/Dif_subs_adim);
 deltap=4.0/(1.0+3.0*rel_cond/Dif_subs_adim);
 ipsil=rel_condm1d2*gamma_2*delta;
 ipsilp=gamma_2*deltap*rel_condm1d2;
 ipsilc=deltap*gamma_2*rel_cond;

  // sample + substrate
 for(ii=1;ii<n_r_sample-1;ii++)
	{
	 d_iM=1.0+0.5/ii;
	 d_im=1.0-0.5/ii;
	 jj=n_z_tot-1;
	 if(Boundary_cond & BTM_PERIODIC)
		{	// coefficients A1 and A2 for the bottom boundary condition dT/dz=0
		 A_1[jj]=1.0;
		 A_2[jj]=0.0;	
		}
	 else if(Boundary_cond & BTM_FLOW_CT)
		{	// coefficients A1 and A2 for the bottom boundary condition dT/dz=f_const
		 A_1[jj]=1.0;
		 A_2[jj]=f_const;	
		}
	 else 
		{	// coefficients A1 and A2 for the bottom boundary condition T=T_furnace
		 A_1[jj]=0.0;
		 A_2[jj]=T_furnace;	
		}
	 // bottom vessel wall
	 for(jj=n_z_tot-2;jj>=n_z_sample;jj--)  
		{
		 tmp=1.0+2.0*gamma_subs_2-gamma_subs_2*A_1[jj+1];
 		 a_temp=A_2[jj+1]+d_iM*Temp_a[ii+1][jj]+d_im*Temp_a[ii-1][jj]-2.0*Temp_a[ii][jj];
		 A_2[jj]=(gamma_subs_2*a_temp+Temp_a[ii][jj])/tmp;
		 A_1[jj]=gamma_subs_2/tmp;
		}
	 if(n_z_subs>0)	// interface
		{
		 jj=n_z_sample-1;
		 if(Reaction_exists) beta=(alfa_n[ii][jj]-alfa_a[ii][jj])/H_cap_samp_adim/2.0*delta;
		 else beta=0.0;
		 beta+=Temp_a[ii][jj];
		 beta+=rel_cond*delta*gamma_2*A_2[jj+1];
 	 	 tmp=1.0-delta*gamma_2*rel_cond*A_1[jj+1]+delta*gamma_2*(1.0+rel_cond);
		 a_temp=gamma_2*delta*rel_condm1d2*(d_iM*Temp_a[ii+1][jj]+d_im*Temp_a[ii-1][jj]-2.0*Temp_a[ii][jj]);
		 A_2[jj]=(a_temp+beta)/tmp;
		 A_1[jj]=delta*gamma_2/tmp;
		}
	 //  sample
	 if(Reaction_exists)		
		{
		 for(jj=n_z_sample-2;jj>0;jj--)  
			{
		     beta=(alfa_n[ii][jj]-alfa_a[ii][jj])/H_cap_samp_adim;
			 tmp=1.0+2.0*gamma_2-gamma_2*A_1[jj+1];
			 a_temp=A_2[jj+1]+d_iM*Temp_a[ii+1][jj]+d_im*Temp_a[ii-1][jj]-2.0*Temp_a[ii][jj];
			 A_2[jj]=(gamma_2*a_temp+beta+Temp_a[ii][jj])/tmp;
			 A_1[jj]=gamma_2/tmp;
			}
		}
	 else
		{
		 for(jj=n_z_sample-2;jj>0;jj--)  
			{
			 tmp=1.0+2.0*gamma_2-gamma_2*A_1[jj+1];
 			 a_temp=A_2[jj+1]+d_iM*Temp_a[ii+1][jj]+d_im*Temp_a[ii-1][jj]-2.0*Temp_a[ii][jj];
			 A_2[jj]=(gamma_2*a_temp+Temp_a[ii][jj])/tmp;
			 A_1[jj]=gamma_2/tmp;
			}
		}
	 if(Boundary_cond & TOP_PERIODIC) // top boundary condition dT/dz=0
		 Temp_n[ii][0]=A_2[1]/(1.0-A_1[1]);	
	 else if(Boundary_cond & TOP_FLOW_CT)	// top boundary condition dT/dz=f_const
		 Temp_n[ii][0]=(A_2[1]-f_const)/(1.0-A_1[1]);
	 else Temp_n[ii][0]=T_furnace; //  top boundary condition T=T_furnace

	 for(jj=1;jj<n_z_tot;jj++) Temp_n[ii][jj]=A_1[jj]*Temp_n[ii][jj-1]+A_2[jj];
	}	// end sample + vessel
 // right vessel wall
 for(ii=n_r_sample;ii<n_r_tot-1;ii++)
	{
	 d_iM=1.0+0.5/ii;
	 d_im=1.0-0.5/ii;
	 jj=n_z_tot-1;
	 if(Boundary_cond & BTM_PERIODIC)
		{	// coefficients A1 and A2 for the bottom boundary condition dT/dz=0
		 A_1[jj]=1.0;
		 A_2[jj]=0.0;	
		}
	 else if(Boundary_cond & BTM_FLOW_CT)
		{	// coefficients A1 and A2 for the bottom boundary condition dT/dz=f_const
		 A_1[jj]=1.0;
		 A_2[jj]=f_const;	
		}
	 else 
		{	// coefficients A1 and A2 for the bottom boundary condition T=T_furnace
		 A_1[jj]=0.0;
		 A_2[jj]=T_furnace;	
		}
	 for(jj=n_z_tot-2;jj>0;jj--)  
		{
		 tmp=1.0+2.0*gamma_subs_2-gamma_subs_2*A_1[jj+1];
 		 a_temp=A_2[jj+1]+d_iM*Temp_a[ii+1][jj]+d_im*Temp_a[ii-1][jj]-2.0*Temp_a[ii][jj];
		 A_2[jj]=(gamma_subs_2*a_temp+Temp_a[ii][jj])/tmp;
		 A_1[jj]=gamma_subs_2/tmp;
		}
	 if(Boundary_cond & TOP_PERIODIC) // top boundary condition dT/dz=0
		 Temp_n[ii][0]=A_2[1]/(1.0-A_1[1]);	
	 else if(Boundary_cond & TOP_FLOW_CT)	// top boundary condition dT/dz=f_const
		 Temp_n[ii][0]=(A_2[1]-f_const)/(1.0-A_1[1]);
	 else Temp_n[ii][0]=T_furnace; //  top boundary condition T=T_furnace

	 for(jj=1;jj<n_z_tot;jj++) Temp_n[ii][jj]=A_1[jj]*Temp_n[ii][jj-1]+A_2[jj];
	}

 if(n_r_subs>0)// interface
	{
	 ii=n_r_sample-1;
	 d_iM=1.0+0.5/ii;
	 d_im=1.0-0.5/ii;
	 jj=n_z_tot-1;
	 if(Boundary_cond & BTM_PERIODIC)
		{	// coefficients A1 and A2 for the bottom boundary condition dT/dz=0
		 A_1[jj]=1.0;
		 A_2[jj]=0.0;	
		}
	 else if(Boundary_cond & BTM_FLOW_CT)
		{	// coefficients A1 and A2 for the bottom boundary condition dT/dz=f_const
		 A_1[jj]=1.0;
		 A_2[jj]=f_const;	
		}
	 else 
		{	/// coefficients A1 and A2 for the bottom boundary condition T=T_furnace
		 A_1[jj]=0.0;
		 A_2[jj]=T_furnace;	
		}
	 for(jj=n_z_tot-2;jj>=n_z_sample;jj--)  
		{
		 tmp=1.0+2.0*gamma_subs_2-gamma_subs_2*A_1[jj+1];
 		 a_temp=A_2[jj+1]+d_iM*Temp_a[ii+1][jj]+d_im*Temp_a[ii-1][jj]-2.0*Temp_a[ii][jj];
		 A_2[jj]=(gamma_subs_2*a_temp+Temp_a[ii][jj])/tmp;
		 A_1[jj]=gamma_subs_2/tmp;
		}
	if(n_z_subs>0)	// interficie
		{
		 jj=n_z_sample-1;
		 ii=n_r_sample-1;
		 d_iM=1.0+0.5/ii;
		 d_im=1.0-0.5/ii;
		 if(Reaction_exists) beta=(alfa_n[ii][jj]-alfa_a[ii][jj])/H_cap_samp_adim/4.0*deltap;
		 else beta=0.0;
		 tmp=1.0+ipsilp+ipsilc*(1.0-A_1[jj+1]);
		 a_temp=ipsilp*d_im*(Temp_a[ii-1][jj]-Temp_a[ii][jj]);
		 a_temp+=ipsilc*d_iM*(Temp_a[ii+1][jj]-Temp_a[ii][jj]);
		 A_2[jj]=(a_temp+beta+ipsilc*A_2[jj+1]+Temp_a[ii][jj])/tmp;
		 A_1[jj]=ipsilp/tmp;
		}	
	 if(Reaction_exists)
		{
		 for(jj=n_z_sample-2;jj>0;jj--) 
			{
			 beta=(alfa_n[ii][jj]-alfa_a[ii][jj])/H_cap_samp_adim/2.0*delta;
			 tmp=1.0-ipsil*(A_1[jj+1]-2.0);
			 a_temp=gamma_2*delta*d_im*(Temp_a[ii-1][jj]-Temp_a[ii][jj]);
			 a_temp+=gamma_2*delta*rel_cond*d_iM*(Temp_a[ii+1][jj]-Temp_a[ii][jj]);
			 A_2[jj]=(a_temp+beta+ipsil*A_2[jj+1]+Temp_a[ii][jj])/tmp;
			 A_1[jj]=ipsil/tmp;
			}
		}
	 else
		{
		 for(jj=n_z_sample-2;jj>0;jj--)  
			{
			 tmp=1.0-ipsil*(A_1[jj+1]-2.0);
 			 a_temp=gamma_2*delta*d_im*(Temp_a[ii-1][jj]-Temp_a[ii][jj]);
			 a_temp+=gamma_2*delta*rel_cond*d_iM*(Temp_a[ii+1][jj]-Temp_a[ii][jj]);
			 A_2[jj]=(a_temp+ipsil*A_2[jj+1]+Temp_a[ii][jj])/tmp;
			 A_1[jj]=ipsil/tmp;
			}
		}
	 if(Boundary_cond & TOP_PERIODIC) // top boundary condition dT/dz=0
		 Temp_n[ii][0]=A_2[1]/(1.0-A_1[1]);	
	 else if(Boundary_cond & TOP_FLOW_CT)	// top boundary condition dT/dz=f_const
		 Temp_n[ii][0]=(A_2[1]-f_const)/(1.0-A_1[1]);
	 else Temp_n[ii][0]=T_furnace; //  top boundary condition T=T_furnace

	 for(jj=1;jj<n_z_tot;jj++) Temp_n[ii][jj]=A_1[jj]*Temp_n[ii][jj-1]+A_2[jj];
	}	// end interface
  
  // boundary condtions 
// due to symetry  dT/dr=0
  for(jj=1;jj<n_z_tot-1;jj++) 	Temp_n[0][jj]=Temp_n[1][jj];		
 
  if(Boundary_cond & RGT_PERIODIC) // right boundary condition dT/dr=0
	 for(jj=1;jj<n_z_tot-1;jj++) 	Temp_n[n_r_tot-1][jj]=Temp_n[n_r_tot-2][jj];		
  else if(Boundary_cond & RGT_FLOW_CT)	//  right boundary condition dT/dr=f_const
	 for(jj=1;jj<n_z_tot-1;jj++) 	Temp_n[n_r_tot-1][jj]=Temp_n[n_r_tot-2][jj]+f_const;		
  else //  right boundary condition T= T_furnace
	 for(jj=1;jj<n_z_tot-1;jj++)  Temp_n[n_r_tot-1][jj]=T_furnace;


  // Boundary condition top corners
  if(Boundary_cond & TOP_PERIODIC) // top boundary condition dT/dz=0
		Temp_n[0][0]=Temp_n[0][1],Temp_n[n_r_tot-1][0]=Temp_n[n_r_tot-1][1];		
  else if(Boundary_cond & TOP_FLOW_CT)	//  top boundary condition dT/dz=f_const
	 Temp_n[0][0]=Temp_n[0][1]-f_const,Temp_n[n_r_tot-1][0]=Temp_n[n_r_tot-1][1]-f_const;		
  else //  top boundary condition T=T_furnace
	 Temp_n[0][0]=T_furnace,Temp_n[n_r_tot-1][0]=T_furnace;

  // Boundary condition bottom corners
  if(Boundary_cond & BTM_PERIODIC) // bottom boundary condition dT/dz=0
	 Temp_n[0][n_z_tot-1]=Temp_n[0][n_z_tot-2],Temp_n[n_r_tot-1][n_z_tot-1]=Temp_n[n_r_tot-1][n_z_tot-2];		
  else if(Boundary_cond & BTM_FLOW_CT)	//  bottom boundary condition dT/dz=f_const
	 Temp_n[0][n_z_tot-1]=Temp_n[0][n_z_tot-2]+f_const,Temp_n[n_r_tot-1][n_z_tot-1]=Temp_n[n_r_tot-1][n_z_tot-2]+f_const;		
  else //  bottom boundary condition T=T_furnace
	 Temp_n[0][n_z_tot-1]=T_furnace,Temp_n[n_r_tot-1][n_z_tot-1]=T_furnace;
}
// calculation of the transformation degree aassuming constant temperature
void Sample::CalculaAlfaIsoterm(double **alfa_n,double **alfa_a,double **Temp)
{
 int ii,jj;
 double tmp,tmp2;

 if(Reaction_exists==false) return;
 for(ii=0;ii<n_r_sample;ii++)	
	{
	 for(jj=0;jj<n_z_sample;jj++)
		{
		 if(alfa_a[ii][jj]>=1.0) alfa_n[ii][jj]=1.0,alfa_a[ii][jj]=1.0;
		 else
			{
			 tmp=H_cap_samp_adim*(exp(-1.0/Temp[ii][jj]))*Inc_t_2;
#ifdef __F1__
			 alfa_n[ii][jj]=1.0-(1.0-alfa_a[ii][jj])*exp(-tmp);
#endif
#ifdef __F2__
			 alfa_n[ii][jj]=1.0-1.0/(tmp+1.0/(1.0-alfa_a[ii][jj]));
#endif
#ifdef __FN__
			 tmp2=pow((1.0-alfa_a[ii][jj]),_cof_n);
			 tmp2=(1.0+tmp*(_cof_n-1.0)*tmp2-alfa_a[ii][jj])/tmp2;
			 tmp=1.0/(1.0-_cof_n);
			 alfa_n[ii][jj]=1.0-pow(tmp2,tmp);
#endif
#ifdef __KJMA__
			 tmp2=pow(-log(1.0-alfa_a[ii][jj]),1.0/_cof_n);
			 tmp2=pow(tmp+tmp2,_cof_n);
			 alfa_n[ii][jj]=1.0-exp(-tmp2);
#endif
#ifdef __PN__
			 tmp2=pow(alfa_a[ii][jj],1.0/_cof_n);
			 alfa_n[ii][jj]=pow(tmp+tmp2,_cof_n);
#endif
#ifdef __F0__
			 alfa_n[ii][jj]=alfa_a[ii][jj]+tmp;
#endif
#ifdef __D1__
			 alfa_n[ii][jj]=sqrt(tmp+alfa_a[ii][jj]*alfa_a[ii][jj]);
#endif
#ifdef __D3_J__
			 tmp2=pow(1-alfa_a[ii][jj],1.0/3.0);
			 tmp2=1.0-tmp2;
			 tmp2=tmp+tmp2*tmp2;
			 tmp2=1.0-sqrt(tmp2);
			 alfa_n[ii][jj]=1.0-tmp2*tmp2*tmp2;
#endif
			 if(alfa_n[ii][jj]>=1.0) alfa_n[ii][jj]=1.0;
			}
		}
	}
}
// calculation of the transformation degree aassuming a constant temperature rise
void Sample::CalculaAlfaNoIsoterm(double **alfa_n,double **alfa_a,double **Temp_n,double **Temp_a)
{
 int ii,jj;
 double tmp,tmp2;
 
 if(Reaction_exists==false) return;
 for(ii=0;ii<n_r_sample;ii++)	
	{
	 for(jj=0;jj<n_z_sample;jj++)
		{
		 if(alfa_a[ii][jj]>=1.0) alfa_n[ii][jj]=1.0,alfa_a[ii][jj]=1.0;
		 else
			{
			 tmp=Temp_n[ii][jj]-Temp_a[ii][jj];
			 if(fabs(tmp)<1.0e-14)	// T=constant, isothermal
				tmp=H_cap_samp_adim*(exp(-1.0/Temp_a[ii][jj]))*Inc_t_2;
			 else			 // Temperature linear interpolation
				tmp=H_cap_samp_adim*Inc_t_2/tmp*(Int_Temp_Apr(1.0/Temp_n[ii][jj])-Int_Temp_Apr(1.0/Temp_a[ii][jj]));
#ifdef __F1__
			 alfa_n[ii][jj]=1.0-(1.0-alfa_a[ii][jj])*exp(-tmp);
#endif
#ifdef __F2__
			 alfa_n[ii][jj]=1.0-1.0/(tmp+1.0/(1.0-alfa_a[ii][jj]));
#endif
#ifdef __FN__
			 tmp2=pow((1.0-alfa_a[ii][jj]),_cof_n);
			 tmp2=(1.0+tmp*(_cof_n-1.0)*tmp2-alfa_a[ii][jj])/tmp2;
			 tmp=1.0/(1.0-_cof_n);
			 alfa_n[ii][jj]=1.0-pow(tmp2,tmp);
#endif
#ifdef __KJMA__
			 tmp2=pow(-log(1.0-alfa_a[ii][jj]),1.0/_cof_n);
			 tmp2=pow(tmp+tmp2,_cof_n);
			 alfa_n[ii][jj]=1.0-exp(-tmp2);
#endif
#ifdef __PN__
			 tmp2=pow(alfa_a[ii][jj],1.0/_cof_n);
			 alfa_n[ii][jj]=pow(tmp+tmp2,_cof_n);
#endif
#ifdef __F0__
			 alfa_n[ii][jj]=alfa_a[ii][jj]+tmp;
#endif
#ifdef __D1__
			 alfa_n[ii][jj]=sqrt(tmp+alfa_a[ii][jj]*alfa_a[ii][jj]);
#endif
#ifdef __D3_J__
			 tmp2=pow(1-alfa_a[ii][jj],1.0/3.0);
			 tmp2=1.0-tmp2;
			 tmp2=tmp+tmp2*tmp2;
			 tmp2=1.0-sqrt(tmp2);
			 alfa_n[ii][jj]=1.0-tmp2*tmp2*tmp2;
#endif
			 if(alfa_n[ii][jj]>=1.0) alfa_n[ii][jj]=1.0;
			}
		}
	 }
}
// calculation of integral p(x)
double Int_Temp_Apr(double x)
{
 double num,denom;

 num=((((x+40.0)*x+552.0)*x+3168.0)*x+7092.0)*x+4320.0;
 denom=(((((x+42.0)*x+630.0)*x+4200.0)*x+12600.0)*x+15120.0)*x+5040.0;
 return exp(-x)/x*num/denom;
}