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

// implicit method
// 1D axial coordinate
// substrate located on top of the system

#include <stdio.h>
#include <math.h>
#include <conio.h>
#include <string.h>
#include <stdlib.h>
#include "Combustion.h"

void SaltaLinia(FILE *PtrF);

// void constructor
Sample::Sample()
{
 Reaction_exists=true;	
 A_1=NULL;				
 A_2=NULL;				
 n_elem_tot=701;		
 n_elem_subs=0;	
 n_elem_sample=701;
 gamma=100.0;		
 gamma_subs=100.0;		
 Tot_time=200;		
 Cond_sample=0.05;		
 Cond_subs=0.05;		
 H_cap_sample=1000;	     	
 H_cap_subs=1000;	     	
 dens_sample=1000;		
 dens_subs=1000;		
 Ent=2.6e5;        
 K_0=1.05e13;  	
 E_A=160000.0;	        
 Length_sample=0.004;		
 Length_subs=0.0;		
 Avg_alpha=0.0;		
 elap_time=0.0;			
 Vel_esc=20.0/60.0;		
 h_conve=0.0;		
 T_room=293.15;		
 Length_trans=-1.0; 		
 T_min=373.0;	 	
 T_Subs_Ini=373.0;	 	
 N_Front=50;		
 N_Point_Calc=100;	
 N_Point_z=2000;	
 Err_Rel_Max=1e-8;		
 h_div_k=h_conve/Cond_subs;	
 
 gamma_max=gamma;
 gamma=1.0;
 Dif_subs_adim=1.0;
 Boundary_cond= TOP_T_CONST | BTM_FLOW_CT;
 Com_Front_Created=false;
}

Sample::~Sample()
{
 if(A_1!=NULL) delete [] A_1;		
 if(A_2!=NULL) delete [] A_2;
}

// read data from a file
int Sample::ImputData(char *File_Name)
{
 int ii;
 FILE *PF;
 errno_t err;

 err=fopen_s(&PF,File_Name,"r");
 if(err)
  {
   printf("\n\n\n\t\t\tUnable to open the file: %s\n\t\t\tProgram ended",File_Name);
   Help();
   return 0;
  }
 SaltaLinia(PF);
 fscanf_s(PF,"%d",&n_elem_sample);
 SaltaLinia(PF);
 fscanf_s(PF,"%lf",&gamma);
 if(gamma<0.01) gamma=0.01;
 gamma_max=gamma;
 SaltaLinia(PF);
 fscanf_s(PF,"%lf",&Tot_time);
 SaltaLinia(PF);
 fscanf_s(PF,"%lf",&Cond_sample);
 SaltaLinia(PF);
 fscanf_s(PF,"%lf",&H_cap_sample);
 SaltaLinia(PF);
 fscanf_s(PF,"%lf",&dens_sample);
 SaltaLinia(PF);
 fscanf_s(PF,"%lf",&Ent);
 SaltaLinia(PF);
 fscanf_s(PF,"%lf",&K_0);
 SaltaLinia(PF);
 fscanf_s(PF,"%lf",&E_A);
 SaltaLinia(PF);
 fscanf_s(PF,"%lf",&Length_sample);
 SaltaLinia(PF);
 fscanf_s(PF,"%lf",&Cond_subs);
 SaltaLinia(PF);
 fscanf_s(PF,"%lf",&H_cap_subs);
 SaltaLinia(PF);
 fscanf_s(PF,"%lf",&dens_subs);
 SaltaLinia(PF);
 fscanf_s(PF,"%lf",&Length_subs);
 SaltaLinia(PF);
 fscanf_s(PF,"%lf",&Vel_esc);
 SaltaLinia(PF);
 Vel_esc/=60.0;	// heating rate in K/s
 fscanf_s(PF,"%lf",&h_conve);
 SaltaLinia(PF);
 h_div_k=h_conve/Cond_subs;	
 fscanf_s(PF,"%lf",&T_room);
 SaltaLinia(PF);
 fscanf_s(PF,"%lf",&Length_trans);
 SaltaLinia(PF);
 fscanf_s(PF,"%lf",&T_Subs_Ini);
 SaltaLinia(PF);
 fscanf_s(PF,"%lf",&T_min);
 SaltaLinia(PF);
 fscanf_s(PF,"%d",&ii);
 SaltaLinia(PF);
 if(ii==0) Int_Fronts=false;
 else Int_Fronts=true;
 fscanf_s(PF,"%d",&N_Front);
 if(N_Front<10) N_Front=10;
 if(N_Front>999) N_Front=999;
 SaltaLinia(PF);
 fscanf_s(PF,"%d",&N_Point_z);
 if(N_Point_z<100) N_Point_z=100;
 SaltaLinia(PF);
 fscanf_s(PF,"%d",&N_Point_Calc);
 if(N_Point_Calc<1) N_Point_Calc=1;
 SaltaLinia(PF);
 fscanf_s(PF,"%lf",&Err_Rel_Max);
 if(Err_Rel_Max<0.0) Err_Rel_Max*=-1.0;
 if(Err_Rel_Max>LIMIT) Err_Rel_Max=LIMIT;
 SaltaLinia(PF);
 fscanf_s(PF,"%d",&Boundary_cond);
 if(Boundary_cond<0 || Boundary_cond>2) Boundary_cond=0;
 Boundary_cond<<=4;
 SaltaLinia(PF);
 fscanf_s(PF,"%d",&ii);
 if(ii<0 || ii>2) ii=0;
 Boundary_cond|=ii;
 return 1;
}

void Sample::ListData(FILE *PF2)
{
 fputs("\n\t\t\tFull implicit method",PF2);
 fputs("\n\t\t\tSimulation parameters",PF2);
 fprintf(PF2,"\n\tNumber of points of the sample grid %d",n_elem_sample);
 fprintf(PF2,"\n\tMaximum gamma parameter allowed, dimensionless %lg",gamma_max);
 fprintf(PF2,"\n\tmaximun time allowed for integration %lg s",Tot_time);
 fprintf(PF2,"\n\tSample thermal conductivity %lg W/(m K)",Cond_sample);
 fprintf(PF2,"\n\tSample Heat capacity %lg J/(kg K)",H_cap_sample);
 fprintf(PF2,"\n\tSample density %lg kg/m3",dens_sample);
 fprintf(PF2,"\n\tReaction enthalpy %lg J/m3",Ent);
 fprintf(PF2,"\n\tPreexponential term s %lg 1/s",K_0);
 fprintf(PF2,"\n\tActivation energy %lg J/mol",E_A);
 fprintf(PF2,"\n\tSample height %lg m",Length_sample);
#ifdef __F1__
 fputs("\n\tModel: First order reaction",PF2);
#endif
#ifdef __F2__
 fputs("\n\tModel: Second order reaction",PF2);
#endif
#ifdef __FN__
 fprintf(PF2,"\n\tModel: reaction of order %d",_cof_n);
#endif
#ifdef __KJMA__
 fprintf(PF2,"\n\tModel: KJMA of order %d",_cof_n);
#endif
#ifdef __PN__
 fprintf(PF2,"\n\tModel: Power law of order %d",_cof_n);
#endif
#ifdef __D1__
 fputs("\n\tModel: 1D diffusion",PF2);
#endif
#ifdef __D3_J__
 fputs("\n\tModel: Jander's 3D diffusion",PF2);
#endif

 fprintf(PF2,"\n\tSubstrate thermal conductivity %lg W/(m K)",Cond_subs);
 fprintf(PF2,"\n\tSubstrate Heat capacity %lg J/(kg K)",H_cap_subs);
 fprintf(PF2,"\n\tSubstrate density %lg kg/m3",dens_subs);
 fprintf(PF2,"\n\tSubstrate height %lg m",Length_subs);
 fprintf(PF2,"\n\tConstant temperature rise %lg K/min",Vel_esc*60.0);
 fprintf(PF2,"\n\tConvection coefficient %lg W/m^2/K",h_conve);
 fprintf(PF2,"\n\tSourrounding gas temperature %lg K",T_room*f_Temp);
 fprintf(PF2,"\n\tInitial sample temperature: %lg K",T_min);
 if(Int_Fronts) fputs("\n\tSaves profiles at constant degree of transformation intervals.",PF2);
 else fputs("\n\tSaves profiles at constant tempeterature intervals.",PF2);
 fprintf(PF2,"\n\tNumber of profiles recorded: %d",N_Front);
 fprintf(PF2,"\n\tNumber of points of the time evolution file: %d",N_Point_z);
 fprintf(PF2,"\n\tMaximum relative error allowed between iterations %lg\n",Err_Rel_Max);
 fputs("\t\tBoundary conditions\n",PF2);
 if(Boundary_cond & TOP_PERIODIC) fputs("\tTop side, zero heat flow (periodic)\n",PF2);
 else if(Boundary_cond & TOP_FLOW_CT) fputs("\tTop side, constant heat flow (convection)\n",PF2);
 else fputs("\tTop side, constant temperature= T furnace\n",PF2);
 if(Boundary_cond & BTM_PERIODIC) fputs("\tBottom side, zero heat flow (periodic)\n",PF2);
 else if(Boundary_cond & BTM_FLOW_CT) fputs("\tBottom side, constant heat flow (convection)\n",PF2);
 else fputs("\tBottom side, constant temperature= T furnace\n",PF2);

 if(Reaction_exists)
	{
	 fprintf(PF2,"\tInitial transformed fraction: %lg",alfa_in);
	 fprintf(PF2,"\n\tMaximun and initial sample temperature: %lg i %lg [K]",T_max,T_min);
	 fprintf(PF2,"\n\tDimensionless heat capacity: %lg",H_cap_samp_adim);
 	 fprintf(PF2,"\n\tDimensionless heating rate: %lg",Vel_esc_adim);
	}
 else
	 fputs("\tTher is no chemical reaction.",PF2);
 fprintf(PF2,"\n\tTemperature scale factor: %lg [K]",f_Temp);
 fprintf(PF2,"\n\tTime scale factor: %lg [s]",f_time);
 fprintf(PF2,"\n\tSpace scale factor: %lg [m]",f_space);
 fprintf(PF2,"\n\tSpace and time steps: %lg [m] i %lg [s]\n\n",Inc_z*f_space,Inc_t*f_time);
 if(n_elem_0>0)
	{
	 fputs("\n\tSHS, ignitinon induced\n\n",PF2);
	 fprintf(PF2,"\n\tIgnition region %lg (dimensionless)",Length_trans);
	}
}

int Sample::Data_Initialization()
{
	// calculation of the system dimensions
 Inc_z=Length_sample/(double)(n_elem_sample-1);
 Length_sample=Inc_z*(n_elem_sample-1.0);
 n_elem_subs=(int)(Length_subs/Inc_z)+1;
 if(n_elem_subs<2) 
	{
	 n_elem_subs=0;
	 Length_subs=0.0;
	 n_elem_tot=n_elem_sample;
	}
 else
	{
	 Length_subs=Inc_z*(n_elem_subs-1.0);
	 n_elem_tot=n_elem_subs+n_elem_sample-1;
	}
 	// memory allocation
 if(A_1!=NULL) delete [] A_1;		
 if(A_2!=NULL) delete [] A_2;
 A_1= new double  [n_elem_tot]; 
 A_2= new double  [n_elem_tot]; 
 
 if(A_1==NULL || A_2==NULL)
	{
	 puts("\n\n\n\t\t\tUnable to allocate the memory\n\t\t\tProgram ended");
 	 return 0;
	}

 // Initialization of dimensionless parameters
 Dif_sample=Cond_sample/dens_sample/H_cap_sample;
 Dif_subs=Cond_subs/dens_subs/H_cap_subs;
 h_div_k=h_conve/Cond_subs;	
 if(Ent>0.0)	// Normalization when a reaction exists
	{
 	 Reaction_exists=true;
	 f_time=H_cap_sample*E_A/R_G/Ent/K_0;	// time scale factor
	 f_Temp=E_A/R_G;	// temperature scale factor
	 f_space=sqrt(Cond_sample*E_A/R_G/Ent/K_0/dens_sample);	// space scale factor
	 H_cap_samp_adim=H_cap_sample*E_A/R_G/Ent;
	 T_min_adim=T_min/f_Temp;
	 T_room/=f_Temp;
	 if(Vel_esc>0.0) alfa_in=1-exp(-(H_cap_samp_adim/Vel_esc)*Int_Temp_Apr(1/T_min_adim));	
	 else alfa_in=0.0;
	 T_max_adim=T_min_adim+(1.0/H_cap_samp_adim)*(1.0-alfa_in);	// adiabatic temperature
	 T_max=T_max_adim*f_Temp;
	 T_furnace_in=T_Subs_Ini/f_Temp;
	 Vel_esc_adim=Vel_esc*f_time/f_Temp;
 	}
 else	// no chemical reaction
	{
 	 Reaction_exists=false;
	 Int_Fronts=false;
	 f_space=Length_sample;		// space scale factor
	 f_time=Length_sample*Length_sample/Cond_sample*dens_sample*H_cap_sample;	// time scale factor
	 f_Temp=1.0;	// temperature scale factor
	 T_max=T_Subs_Ini+Vel_esc*Tot_time;
	 alfa_in=0.0;
	 T_min_adim=T_min/f_Temp;
	 T_max_adim=T_max/f_Temp;
	 T_furnace_in=T_Subs_Ini/f_Temp;
	 Vel_esc_adim=Vel_esc*f_time/f_Temp;
	}
 
 // Dimensionless parameters
 Inc_z/=f_space;
 Length_sample_adim=Inc_z*(n_elem_sample-1.0);
 if(n_elem_subs>0) Length_subs_adim=Inc_z*(n_elem_subs-1.0);
 else Length_subs_adim=0.0;
 Length_trans/=f_space;
 h_div_k*=f_space;
 Tot_time_adim=Tot_time/f_time;
 Dif_subs_adim=Dif_subs/Dif_sample;
 
// simulation parameters initialization
 gamma=100.0;
 T_furnace=T_furnace_in;
 f_const=h_div_k*Inc_z*T_room;	// boundary condition flow term for convection
 Inc_t=gamma*Inc_z*Inc_z;
 gamma_subs=gamma*Dif_subs_adim;

 if(Length_trans>0.0)	// SHS, ignition induced
	{
	 n_elem_0=(int)(Length_trans/Inc_z)+1;
	 if(n_elem_0>n_elem_sample) n_elem_0=n_elem_sample;
	}
 else n_elem_0=0;

 if(Int_Fronts)
	{
	 Inc_Front=(1.0-alfa_in)/(double)N_Front;
	 Val_Front=alfa_in;
	 Inc_N_Point_z=(1.0-alfa_in)/(double)N_Point_z;
	 Val_N_Point_z=alfa_in;
	}
 else
	{
	 Inc_Front=(Tot_time)/(double)N_Front;
	 Val_Front=0.0;
	 Inc_N_Point_z=(Tot_time)/(double)N_Point_z;
	 Val_N_Point_z=0.0;
	}

 elap_time=0.0;
 elap_time_adim=0.0;
 Avg_alpha=0.0;
 k_front=0;
 Com_Front_Created=false;
 return 1;
}
/*---------------------------------------------------------------------------
    Some help on the screen
---------------------------------------------------------------------------*/

void Help()
{
 puts("\n\n\t1D simulation program, axial coordinate.");
 puts("\n\n\tInput format:");
 puts("\n\t\t1D_cart 'file''\n\t");
 puts("The file should contain the system parameters.\n");
  _getch();
}

void SaltaLinia(FILE *PtrF)
{
 int ch;

 for(ch=fgetc(PtrF);ch!='\n' && ch!='\r' && ch!=EOF;) ch=fgetc(PtrF);
}

void Sample::RegiFront(double *Temp,double *alfa,FILE *PF2,char *cami)
{
 char  buffer[_MAX_PATH],Nom[_MAX_PATH];
 int i;
 FILE *PF;
 errno_t err;
 double m_posadim;
 
 if(k_front>999) return;
 sprintf_s(Nom,"Front%03d.dat",k_front); 	
 strcpy_s(buffer,cami);
 strcat_s(buffer,Nom);
 err=fopen_s(&PF,buffer,"wt");
 if(err)
  {
   printf("\n\n\n\t\t\tUnable to create the file: %s\n\t\t\tProgram ended",buffer);
   return;
  }
 fprintf(PF2,"\tSaving the profile at %s time %lg [s] alpha=%lg\n",Nom,elap_time,Avg_alpha);
 printf("\tSaving the profile at %s\n",Nom);
 fprintf(PF,"Position_z[adim] Position_z[m] Temperature[adim] Temperature[K] alpha=%.3lf \n",Avg_alpha);
 
 // sample

 fprintf(PF,"%lg %lg %lg %lg %lg\n",0.0,0.0,Temp[0],f_Temp*Temp[0],alfa[0]);
 for(i=1;i<n_elem_sample-1;i++)
	{
 	 m_posadim=Inc_z*(double)i;
	 fprintf(PF,"%lg %lg %lg %lg %lg\n",m_posadim,f_space*m_posadim,Temp[i],f_Temp*Temp[i],alfa[i]);
	}
 
 // substrate
 for(i=n_elem_sample-1;i<n_elem_tot-1;i++)
	{
 	 m_posadim=Inc_z*(double)i;
	 fprintf(PF,"%lg %lg %lg %lg %lg\n",m_posadim,f_space*m_posadim,Temp[i],f_Temp*Temp[i],0.0);
	}
 m_posadim=Inc_z*(double)(n_elem_tot-1);
 fprintf(PF,"%lg %lg %lg %lg %lg\n",m_posadim,f_space*m_posadim,Temp[n_elem_tot-1],Temp[n_elem_tot-1]*f_Temp,0.0);

 fclose(PF);
 k_front++;	
 Val_Front+=Inc_Front;
}
// splits filename and extension
 int SeparaNomExt(char *Fitxer,char *Nom)
 {
  int i,j;

  j=strlen(Fitxer);
  if(j<0 || j>_MAX_PATH) return 0;
  i=j;
  do
	{
	 i--;
	}while(Fitxer[i]!='.' && i>-1);
 strcpy_s(Nom,_MAX_PATH,Fitxer);
 if(i==0) i=j;
 else Nom[i]='\x0';
 if(j-i!=4) return 1;
 if(Fitxer[i+1]=='L' || Fitxer[i+1]=='l')
	{
 	 if(Fitxer[i+1]=='O' || Fitxer[i+2]=='o')
		{
		if(Fitxer[i+3]=='G' || Fitxer[i+3]=='g')
			{
			 puts("\n\tThe file with parameter data cannot have extension .log\n");
			 return 0;
			}
		}
	}
  return 1;
 }
int TreuPath(char *Fitxer,char *Nom)
 {
  int i,j;

  j=strlen(Fitxer);
  if(j<0 || j>_MAX_PATH) return 0;
  i=j;
  do
	{
	 i--;
	}while(Fitxer[i]!='/' && Fitxer[i]!='\\' &&i>-1);
 strcpy_s(Nom,_MAX_PATH,Fitxer);
 if(i==0) 
	{
	 i=j;
	 Nom[0]='\x0';
	 return 0;
	}
 else Nom[i+1]='\x0';
 return 1;
 }
 
// Locates the position of the propagation fron if any
// Prop_Up=true searches a front propagating upwards, false downwards

void Sample::RegPosFront(FILE *PF,double *alfa,bool Prop_Up)
{
 int ii;
 double a_temp1,a_temp2,m_temp;
 bool Front_found,Begin_found;

 if(Reaction_exists==false) return;

 Front_found=false;
 Begin_found=false;
 if(Prop_Up)   // front propagating upwards
	{
	 ii=n_elem_sample-1;
	 do
		{
		 if(alfa[ii]>0.66) 
			{
			 if(ii==n_elem_sample-1) return;
			 m_temp=(alfa[ii]-alfa[ii+1]);
			 if(m_temp>6.0e-15)
				{
				 a_temp1=-(0.66-alfa[ii+1])/m_temp+(double)ii+1.0;
				 a_temp1*=Inc_z;
				 Begin_found=true;
				}
			}
		 ii--;
		}while(ii>=0 && Begin_found==false);
	 if(ii<1) return;
	 do
		{
		 if(alfa[ii]>0.99) 
			{
			 m_temp=(alfa[ii]-alfa[ii+1]);
			 if(m_temp>6.0e-15) a_temp2=ii;
			 else a_temp2=-(0.99-alfa[ii+1])/m_temp+(double)ii+1.0;
			 a_temp2*=Inc_z;
			 Front_found=true;
			}
		 ii--;
		}while(ii>=0 && Front_found==false);
	 }
 else 
	{		 // front propagating downwards
	 ii=0;
	 do
		{
		 if(alfa[ii]>0.66) 
			{
			 if(ii==0) return;
			 m_temp=(alfa[ii]-alfa[ii-1]);
			 if(m_temp>6.0e-15)
				{
			     a_temp2=(0.66-alfa[ii-1])/m_temp+(double)ii-1.0;
				 a_temp2*=Inc_z;
				 Begin_found=true;
				}
			}
		 ii++;
		}while(ii<n_elem_sample && Begin_found==false);
	 if(ii>=n_elem_sample) return;
	 do
		{
		 if(alfa[ii]>0.99) 
			{
			 m_temp=(alfa[ii]-alfa[ii-1]);
			 if(m_temp<6e-15) a_temp1=ii;
			 else a_temp1=(0.99-alfa[ii-1])/m_temp+(double)ii-1.0;
			 a_temp1*=Inc_z;
			 Front_found=true;
			}
		 ii++;
		}while(ii<n_elem_sample && Front_found==false);
	}
 if(Front_found)
	{
	 fprintf(PF,"%lf %lg %lg %lg %lg\n",elap_time,T_furnace*f_Temp,Avg_alpha,a_temp1*f_space,(a_temp1-a_temp2)*f_space);     
	 if(Com_Front_Created==false)
		{
		 Com_Front_Created=true;
		 puts("\n\n\tWarning!!!!!\n\n\tA combustion front has been created!!!\n\n\n");
		}
	}
}