/***************************************************************************
 *   Copyright (C) 2014 by Jordi Farjas					   *
 ***************************************************************************/

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

// iteration to solve the implicit scheme
// n_elem_tot:		total number of points of the space grid
// 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(double *alfa_n,double *alfa_a,double *Temp_n,double *Temp_a)
{
 int ii;
 double beta,tmp,delta,rel_cond,d_iM,d_im;

 if(Boundary_cond & RGT_PERIODIC)
	{	// coefficients A1 and A2 for the right boundary condition dT/dz=0
	 A_1[n_elem_tot-1]=1.0;
	 A_2[n_elem_tot-1]=0.0;	
	}
 else if(Boundary_cond & RGT_FLOW_CT)
	{	// coefficients A1 and A2 for the right boundary condition dT/dz=f_const
	 A_1[n_elem_tot-1]=1.0;
	 A_2[n_elem_tot-1]=f_const;	
	}
 else 
	{	// coefficients A1 and A2 for the right boundary condition T=T_furnace
	 A_1[n_elem_tot-1]=0.0;
	 A_2[n_elem_tot-1]=T_furnace;	
	}
	
// substrate coefficients A1 and A2
 for(ii=n_elem_tot-2;ii>=n_elem_sample;ii--) 
	{
	 d_iM=1.0+0.5/ii;
	 d_im=1.0-0.5/ii;
	 tmp=1.0+2.0*gamma_subs-gamma_subs*d_iM*A_1[ii+1];
	 A_2[ii]=(gamma_subs*d_iM*A_2[ii+1]+Temp_a[ii])/tmp;
	 A_1[ii]=gamma_subs*d_im/tmp;
	}

// coefficients A1 and A2 at the interface
 if(n_elem_subs>0)
	{ 
	 d_iM=1.0+0.5/(n_elem_sample-1.0);
	 d_im=1.0-0.5/(n_elem_sample-1.0);
	 rel_cond=Cond_subs/Cond_sample;
	 delta=2.0/(1.0+rel_cond/Dif_subs_adim);
	 if(Reaction_exists) beta=(alfa_n[n_elem_sample-1]-alfa_a[n_elem_sample-1])/H_cap_samp_adim/2.0*delta;
	 else beta=0.0;
	 beta+=Temp_a[n_elem_sample-1];
	 beta+=rel_cond*delta*gamma*d_iM*A_2[n_elem_sample];
	 tmp=1.0-delta*gamma*rel_cond*d_iM*A_1[n_elem_sample]+delta*gamma*(d_im+rel_cond*d_iM);
	 A_2[n_elem_sample-1]=beta/tmp;
	 A_1[n_elem_sample-1]=delta*d_im*gamma/tmp;
	}
// substrate coefficients A1 and A2
 if(Reaction_exists)
	{
	 for(ii=n_elem_sample-2;ii>0;ii--) 
		{
		 d_iM=1.0+0.5/ii;
		 d_im=1.0-0.5/ii;
		 beta=(alfa_n[ii]-alfa_a[ii])/H_cap_samp_adim;
		 tmp=1.0+2.0*gamma-gamma*d_iM*A_1[ii+1];
		 A_2[ii]=(gamma*d_iM*A_2[ii+1]+Temp_a[ii]+beta)/tmp;
		 A_1[ii]=gamma*d_im/tmp;
		}
	}
 else
	{
	 for(ii=n_elem_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-gamma*d_iM*A_1[ii+1];
		 A_2[ii]=(gamma*d_iM*A_2[ii+1]+Temp_a[ii])/tmp;
		 A_1[ii]=gamma*d_im/tmp;
		}
	}
 
 // calculation of new temperature at the bottom boundary

// Due to symmetry at the z axe dT/dr=0
 Temp_n[0]=A_2[1]/(1.0-A_1[1]);	
  
 // calculation of new temperatures  

 for(ii=1;ii<n_elem_tot;ii++) Temp_n[ii]=A_1[ii]*Temp_n[ii-1]+A_2[ii];
}
// calculation of the transformation degree aassuming constant temperature
void Sample::CalculaAlfaIsoterm(double *alfa_n,double *alfa_a,double *Temp)
{
 int ii;
 double tmp,tmp2;

 if(Reaction_exists==false) return;
 for(ii=0;ii<n_elem_sample;ii++)
	{
	 if(alfa_a[ii]>=1.0) alfa_n[ii]=1.0,alfa_a[ii]=1.0;
	 else
		{
		 tmp=H_cap_samp_adim*(exp(-1.0/Temp[ii]))*Inc_t;
#ifdef __F1__
		 alfa_n[ii]=1.0-(1.0-alfa_a[ii])*exp(-tmp);
#endif
#ifdef __F2__
		 alfa_n[ii]=1.0-1.0/(tmp+1.0/(1.0-alfa_a[ii]));
#endif
#ifdef __FN__
		 tmp2=pow((1.0-alfa_a[ii]),_cof_n);
		 tmp2=(1.0+tmp*(_cof_n-1.0)*tmp2-alfa_a[ii])/tmp2;
		 tmp=1.0/(1.0-_cof_n);
		 alfa_n[ii]=1.0-pow(tmp2,tmp);
#endif
#ifdef __KJMA__
		 tmp2=pow(-log(1.0-alfa_a[ii]),1.0/_cof_n);
		 tmp2=pow(tmp+tmp2,_cof_n);
		 alfa_n[ii]=1.0-exp(-tmp2);
#endif
#ifdef __PN__
		 tmp2=pow(alfa_a[ii],1.0/_cof_n);
		 alfa_n[ii]=pow(tmp+tmp2,_cof_n);
#endif
#ifdef __F0__
		 alfa_n[ii]=alfa_a[ii]+tmp;
#endif
#ifdef __D1__
		 alfa_n[ii]=sqrt(tmp+alfa_a[ii]*alfa_a[ii]);
#endif
#ifdef __D3_J__
		 tmp2=pow(1-alfa_a[ii],1.0/3.0);
		 tmp2=1.0-tmp2;
		 tmp2=tmp+tmp2*tmp2;
		 tmp2=1.0-sqrt(tmp2);
		 alfa_n[ii]=1.0-tmp2*tmp2*tmp2;
#endif
		 if(alfa_n[ii]>=1.0) alfa_n[ii]=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;
 double tmp,tmp2;
 
 if(Reaction_exists==false) return;
 for(ii=0;ii<n_elem_sample;ii++)	
	{
	 if(alfa_a[ii]>=1.0) alfa_n[ii]=1.0,alfa_a[ii]=1.0;
	 else
		{
		 tmp=Temp_n[ii]-Temp_a[ii];
		 if(fabs(tmp)<1.0e-14)	// T=constant, isothermal
			tmp=H_cap_samp_adim*(exp(-1.0/Temp_a[ii]))*Inc_t;
		 else			 // Temperature linear interpolation
			tmp=H_cap_samp_adim*Inc_t/tmp*(Int_Temp_Apr(1.0/Temp_n[ii])-Int_Temp_Apr(1.0/Temp_a[ii]));
#ifdef __F1__
		 alfa_n[ii]=1.0-(1.0-alfa_a[ii])*exp(-tmp);
#endif
#ifdef __F2__
		 alfa_n[ii]=1.0-1.0/(tmp+1.0/(1.0-alfa_a[ii]));
#endif
#ifdef __FN__
		 tmp2=pow((1.0-alfa_a[ii]),_cof_n);
		 tmp2=(1.0+tmp*(_cof_n-1.0)*tmp2-alfa_a[ii])/tmp2;
		 tmp=1.0/(1.0-_cof_n);
		 alfa_n[ii]=1.0-pow(tmp2,tmp);
#endif
#ifdef __KJMA__
		 tmp2=pow(-log(1.0-alfa_a[ii]),1.0/_cof_n);
		 tmp2=pow(tmp+tmp2,_cof_n);
		 alfa_n[ii]=1.0-exp(-tmp2);
#endif
#ifdef __PN__
		 tmp2=pow(alfa_a[ii],1.0/_cof_n);
		 alfa_n[ii]=pow(tmp+tmp2,_cof_n);
#endif
#ifdef __F0__
		 alfa_n[ii]=alfa_a[ii]+tmp;
#endif
#ifdef __D1__
		 alfa_n[ii]=sqrt(tmp+alfa_a[ii]*alfa_a[ii]);
#endif
#ifdef __D3_J__
		 tmp2=pow(1-alfa_a[ii],1.0/3.0);
		 tmp2=1.0-tmp2;
		 tmp2=tmp+tmp2*tmp2;
		 tmp2=1.0-sqrt(tmp2);
		 alfa_n[ii]=1.0-tmp2*tmp2*tmp2;
#endif
		 if(alfa_n[ii]>=1.0) alfa_n[ii]=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;
}