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

// Alternating direction implicit method
// cylindrical coordinates
// Sourrounding vessel at the bottom and right sides

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

void SaltaLinia(FILE *PtrF);
bool CreateImageAlfa(int n_x,int n_y,double **alfa, BYTE *img_RGB,BYTE *img_BMP,long size_RGB);

// void constructor
Sample::Sample()
{
 Reaction_exists=true;	
 A_1=NULL;				
 A_2=NULL;				
 n_r_sample=701;		
 n_z_sample=351;		
 n_r_subs=0;		
 n_z_subs=0;		
 n_r_tot=701;		
 n_z_tot=351;		
 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;
 Radius_sample=0.004;
 Btm_wall_thk=0.0;
 Rgt_wall_thk=0.0;
 Avg_alpha=0.0;		
 elap_time=0.0;			
 Vel_esc=20.0/60.0;		
 h_conve=0.0;		
 T_room=293.15;		
 Radius_trans=-1.0;
 Height_trans=-1.0;
 T_min=373.0;	 	
 T_Subs_Ini=373.0;	 	
 col_z=0;			
 col_r=0;
 N_Front_bmp=300;
 N_Front_prof=50;
 N_Point_Calc=100;	
 N_Point_s=2000;	
 Err_Rel_Max=1e-14;		
 h_div_k=h_conve/Cond_subs;	
 img_BMP=NULL;
 img_RGB=NULL;
 
 gamma_max=gamma;
 gamma=1.0;
 Dif_subs_adim=1.0;
 Boundary_cond= TOP_PERIODIC | BTM_T_CONST | RGT_FLOW_CT;
 Com_Front_Created=false;
}

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

// 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_r_sample);
 SaltaLinia(PF);
 fscanf_s(PF,"%d",&n_z_sample);
 SaltaLinia(PF);
 fscanf_s(PF,"%lf",&gamma);
 if(gamma<0.01) gamma=0.01;
 gamma_max=gamma;
 gamma=100.0;
 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",&Radius_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",&Rgt_wall_thk);
 SaltaLinia(PF);
 fscanf_s(PF,"%lf",&Btm_wall_thk);
 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",&Radius_trans);
 SaltaLinia(PF);
  fscanf_s(PF,"%lf",&Height_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_bmp);
 if(N_Front_bmp<0) N_Front_bmp=0;
 if(N_Front_bmp>999) N_Front_bmp=999;
 SaltaLinia(PF);
 fscanf_s(PF,"%d",&N_Front_prof);
 if(N_Front_prof<0) N_Front_prof=0;
 if(N_Front_prof>99) N_Front_prof=99;
 SaltaLinia(PF);
 fscanf_s(PF,"%d",&N_Point_s);
 if(N_Point_s<100) N_Point_s=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;
 Boundary_cond<<=4;
 SaltaLinia(PF);
 fscanf_s(PF,"%d",&ii);
 if(ii<0 || ii>2) ii=0;
 Boundary_cond|=ii;
 Boundary_cond<<=4;
 SaltaLinia(PF);
 fscanf_s(PF,"%d,%d",&col_z,&col_r);
 if(col_r<0) col_r=0;
 if(col_z<0) col_z=0;
 if(col_r>n_r_sample-1) col_r=n_r_sample-1;
 if(col_z>n_z_sample-1) col_z=n_z_sample-1;
 fclose(PF);
 return 1;
}

void Sample::ListData(FILE *PF2)
{
 fputs("\n\tAlternating direction implicit method (cylindrical coordinates)",PF2);
 fputs("\n\t\t\tSimulation parameters",PF2);
 fprintf(PF2,"\n\tNumber of points of the radial sample grid  %d",n_r_sample);
 fprintf(PF2,"\n\tNumber of points of the axial sample grid  %d",n_z_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 radius and height %lg, %lg m",Radius_sample,Height_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\tVessel thermal conductivity %lg W/(m K)",Cond_subs);
 fprintf(PF2,"\n\tVessel Heat capacity %lg J/(kg K)",H_cap_subs);
 fprintf(PF2,"\n\tVessel density %lg kg/m3",dens_subs);
 fprintf(PF2,"\n\tVessel thickness (right and bottom) %lg, %lg m",Rgt_wall_thk,Btm_wall_thk);
 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_prof);
 fprintf(PF2,"\n\tNumber of images recorded: %d",N_Front_bmp);
 fprintf(PF2,"\n\tNumber of points of the time evolution file: %d",N_Point_s);
 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 substrate\n",PF2);
 if(Boundary_cond & RGT_PERIODIC) fputs("\tRight side, zero heat flow (periodic)\n",PF2);
 else if(Boundary_cond & RGT_FLOW_CT) fputs("\tRight side, constant heat flow (convection)\n",PF2);
 else fputs("\tRight side, constant temperature= T substrate\n",PF2);
 fprintf(PF2,"\tFile %d and column %d saved at 2D profiles\n",col_z,col_r);

 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_s*f_space,Inc_t*f_time);
 if(n_r_0>0)
	{
	 fputs("\n\tSHS, ignitinon induced\n\n",PF2);
	 fprintf(PF2,"\n\tRadius and height ignited %lg, %lg (m)",Radius_trans*f_space,Height_trans*f_space);
	}
}

int Sample::Data_Initialization()
{
	// calculation of the system dimensions
 Inc_s=Radius_sample/(double)(n_r_sample-1);
 Radius_sample=Inc_s*(n_r_sample-1.0);
 Height_sample=Inc_s*(n_z_sample-1.0);
 n_z_subs=(int)(Btm_wall_thk/Inc_s)+1;
 if(n_z_subs<2) 
	{
	 n_z_subs=0;
	 Btm_wall_thk=0.0;
	 n_z_tot=n_z_sample;
	}
 else
	{
	 Btm_wall_thk=Inc_s*(n_z_subs-1.0);
	 n_z_tot=n_z_subs+n_z_sample-1;
	}
 n_r_subs=(int)(Rgt_wall_thk/Inc_s)+1;
 if(n_r_subs<2) 
	{
	 n_r_subs=0;
	 Rgt_wall_thk=0.0;
	 n_r_tot=n_r_sample;
	}
 else
	{
	 Rgt_wall_thk=Inc_s*(n_r_subs-1.0);
	 n_r_tot=n_r_subs+n_r_sample-1;
	}
 	// memory allocation
 if(A_1!=NULL) delete [] A_1;		
 if(A_2!=NULL) delete [] A_2;
 if(n_r_tot>n_z_tot)
	{
	 A_1= new double  [n_r_tot]; 
 	 A_2= new double  [n_r_tot]; 
	}
 else
	{
	 A_1= new double  [n_z_tot]; 
 	 A_2= new double  [n_z_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;
	}
 if(img_BMP!=NULL) delete [] img_BMP;
 if(img_RGB!=NULL) delete [] img_RGB;
 size_RGB=n_r_tot*n_z_tot*3;
 img_RGB = new BYTE[size_RGB];
 SizeRGBToBMPB(n_r_tot, n_z_tot, size_BMP);
 img_BMP = new BYTE[size_BMP];
 if(img_RGB==NULL || img_BMP==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=Radius_sample;		// space scale factor
	 f_time=Radius_sample*Radius_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_s/=f_space;
 Radius_sample_adim=Inc_s*(n_r_sample-1.0);
 Height_sample_adim=Inc_s*(n_z_sample-1.0);
 if(n_r_subs>0) Rgt_wall_thk_adim=Inc_s*(n_r_subs-1.0);
 else Rgt_wall_thk_adim=0.0;
 if(n_z_subs>0) Btm_wall_thk_adim=Inc_s*(n_z_subs-1.0);
 else Btm_wall_thk_adim=0.0;
 Radius_trans/=f_space;
 Height_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_s*T_room;	// boundary condition flow term for convection
 gamma_2=gamma/2.0;
 Inc_t_2=gamma_2*Inc_s*Inc_s;
 Inc_t=2.0*Inc_t_2;
 gamma_subs=gamma*Dif_subs_adim;
 gamma_subs_2=gamma_subs/2.0;

 if(Radius_trans>0.0 && Height_trans>0.0)	// SHS, ignition induced
	{
	 n_r_0=(int)(Radius_trans/Inc_s)+1;
	 n_z_0=(int)(Height_trans/Inc_s)+1;
	 if(n_r_0>n_r_sample) n_r_0=n_r_sample;
	 if(n_z_0>n_z_sample) n_z_0=n_z_sample;
	}
 else n_r_0=0,n_z_0=0;

 if(Int_Fronts)
	{
	 Inc_Front_bmp=(1.0-alfa_in)/(double)N_Front_bmp;
	 Val_Front_bmp=alfa_in;
	 Inc_Front_prof=(1.0-alfa_in)/(double)N_Front_prof;
	 Val_Front_prof=alfa_in;
	 Inc_N_Point_s=(1.0-alfa_in)/(double)N_Point_s;
	 Val_N_Point_s=alfa_in;
	}
 else
	{
	 Inc_Front_bmp=(Tot_time)/(double)N_Front_bmp;
	 Val_Front_bmp=0.0;
	 Inc_Front_prof=(Tot_time)/(double)N_Front_prof;
	 Val_Front_prof=0.0;
	 Inc_N_Point_s=(Tot_time)/(double)N_Point_s;
	 Val_N_Point_s=0.0;
	}

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

void Help()
{
 puts("\n\n\t2D simulation program, cylindrical coordinates.");
 puts("\n\n\tInput format:");
 puts("\n\t\t2D_cyl '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_prof>999) return;
 sprintf_s(Nom,"Front_r%03d.dat",k_front_prof); 	
 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\n",Nom);
 printf("\tSaving the profile at %s\n",Nom);
 fprintf(PF,"Position_r[adim] Position_r[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][col_z],f_Temp*Temp[0][col_z],alfa[0][col_z]);
 for(i=1;i<n_r_sample-1;i++)
	{
 	 m_posadim=Inc_s*(double)i;
	 fprintf(PF,"%lg %lg %lg %lg %lg\n",m_posadim,f_space*m_posadim,Temp[i][col_z],f_Temp*Temp[i][col_z],alfa[i][col_z]);
	}
 
 // substrate
 for(i=n_r_sample-1;i<n_r_tot-1;i++)
	{
 	 m_posadim=Inc_s*(double)i;
 	 fprintf(PF,"%lg %lg %lg %lg %lg\n",m_posadim,f_space*m_posadim,Temp[i][col_z],f_Temp*Temp[i][col_z],0.0);
	}
 m_posadim=Inc_s*(double)(n_r_tot-1);
 fprintf(PF,"%lg %lg %lg %lg %lg\n",m_posadim,f_space*m_posadim,Temp[n_r_sample-1][col_z],Temp[n_r_sample-1][col_z]*f_Temp,0.0);
 fclose(PF);

 sprintf_s(Nom,"Front_z%03d.dat",k_front_prof); 	
 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\n",Nom);
 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[col_r][0],f_Temp*Temp[col_r][0],alfa[col_r][0]);
 for(i=1;i<n_z_sample-1;i++)
	{
 	 m_posadim=Inc_s*(double)i;
	 fprintf(PF,"%lg %lg %lg %lg %lg\n",m_posadim,f_space*m_posadim,Temp[col_r][i],f_Temp*Temp[col_r][i],alfa[col_r][i]);
	}
 
 // substrate
 for(i=n_z_sample-1;i<n_z_tot-1;i++)
	{
 	 m_posadim=Inc_s*(double)i;
 	 fprintf(PF,"%lg %lg %lg %lg %lg\n",m_posadim,f_space*m_posadim,Temp[col_r][i],f_Temp*Temp[col_r][i],0.0);
	}
 m_posadim=Inc_s*(double)(n_z_tot-1);
 fprintf(PF,"%lg %lg %lg %lg %lg\n",m_posadim,f_space*m_posadim,Temp[col_r][n_z_sample-1],Temp[col_r][n_z_sample-1]*f_Temp,0.0);
 fclose(PF);
 k_front_prof++;	
 Val_Front_prof+=Inc_Front_prof;
}
// 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;
 }
// This function creates a BMP image from the array "values" of size  n_r_tot n_z_tot
// colors are graded according to Max_Value and Min_Value
// Max_Value: maximum value of the array and Min_Value: minimum value of the array

bool Sample::CreateImatge(double **values,double Max_Value,double Min_Value)
{
 int ii,jj;
 BYTE vermell,vert,blau;
 double tmp;
	
 if (values==NULL)	return false;
 memset (img_RGB, 255, size_RGB );

 for (ii=0;ii<n_r_tot;ii++)
	 for (jj=0;jj<n_z_tot;jj++)
		{
		 tmp=255.0*(values[ii][jj]-Min_Value)/(Max_Value-Min_Value);
		 if(tmp>255.0) tmp=255.0;
		 else  if(tmp<0.0) tmp=0.0;
		 vermell=(BYTE)tmp;
		 blau=255-vermell;
		 if(vermell<128) vert=2*vermell;
		 else vert=2*(255-vermell);
		 img_RGB[3*ii+jj*3*n_r_tot]=vermell;
		 img_RGB[3*ii+1+jj*3*n_r_tot]=vert;
 		 img_RGB[3*ii+2+jj*3*n_r_tot]=blau;
		}
 return ConvertRGBToBMPBuffer ( img_BMP, img_RGB,  n_r_tot, n_z_tot);
}
// This function creates a BMP file from the array "values" of size  n_r_tot n_z_tot
// colors are graded according to Max_Value and Min_Value
// Max_Value: maximum value of the array and Min_Value: minimum value of the array

void Sample::RegBMP(double **values,char *Name,char *cami,double Max_Value,double Min_Value)
{
 char  buffer[_MAX_PATH],Nom[_MAX_PATH];
 FILE *PF;
 errno_t err;

 if(k_front_bmp>999) return;
 sprintf_s(Nom,"%s%03d.bmp",Name,k_front_bmp); 	
 strcpy_s(buffer,cami);
 strcat_s(buffer,Nom);
 err=fopen_s(&PF,buffer,"wb");
 if(err)
  {
   printf("\n\n\n\t\t\tUnable to create the file: %s\n\t\t\tProgram ended",buffer);
   return;
  }
 CreateImatge(values,Max_Value,Min_Value);
 SaveBMP ( PF,img_BMP, n_r_tot, n_z_tot, size_BMP );
 fclose(PF);
}
// Locates the position of the propagation fron if any
// Prop_Rgt=true searches a front propagating rightwards, false leftwards
// Prop_Down=true searches a front propagating downwards, false upwards

void Sample::RegPosFront(FILE *PF,double **alfa,bool Prop_Rgt,bool Prop_Down)
{
 double pos_r,width_r,pos_z,width_z;
 bool Front_found_r,Front_found_z;

 if(Reaction_exists==false) return;

 Front_found_r=false;
 Front_found_z=false;
 Front_found_r=RegPosFront_r(PF,alfa,Prop_Rgt,pos_r,width_r);
 Front_found_z=RegPosFront_z(PF,alfa,Prop_Down,pos_z,width_z);
 if(Front_found_r | Front_found_z) 
	 fprintf(PF,"%lg %lg %lg %lg %lg %lg %lg\n",elap_time,T_furnace*f_Temp,Avg_alpha,
		pos_r*f_space,width_r*f_space,pos_z*f_space,width_z*f_space);
}
// searches in the radial direction in file col_z
// Prop_Rgt=true searches a front propagating rightwards, false leftwards
bool Sample::RegPosFront_r(FILE *PF,double **alfa,bool Prop_Rgt,double &pos_r,double &width_r)
{
 int ii;
 double a_temp1,a_temp2,m_temp;
 bool Front_found,Begin_found;

 pos_r=0.0;
 width_r=0.0;
 Front_found=false;
 Begin_found=false;
 if(Prop_Rgt)  // front propagatin rightwards
	{
	 ii=n_r_sample-2;
	 do
		{
		 if(alfa[ii][col_z]>0.66) 
			{
			 m_temp=(alfa[ii][col_z]-alfa[ii+1][col_z]);
			 if(m_temp>6.0e-15)
				{
				 a_temp1=-(0.66-alfa[ii+1][col_z])/m_temp+(double)ii+1.0;
				 a_temp1*=Inc_s;
				 pos_r=a_temp1;
				 Begin_found=true;
				}
			}
		 ii--;
		}while(ii>=0 && Begin_found==false);
	  if(ii<=0) return Front_found; 
	  do
		{
		 if(alfa[ii][col_z]>0.99) 
			{
			 m_temp=(alfa[ii][col_z]-alfa[ii+1][col_z]);
			 if(m_temp>6.0e-15) a_temp2=ii;
			 else a_temp2=-(0.99-alfa[ii+1][col_z])/m_temp+(double)ii+1.0;
			 a_temp2*=Inc_s;
			 Front_found=true;
			 width_r=a_temp1-a_temp2;
			}
		 ii--;
		}while(ii>=0 && Front_found==false);
	 }
 else 
	{		 // front propagatin rightwards
	 ii=0;
	 do
		{
		 if(alfa[ii][col_z]>0.66) 
			{
			 if(ii==0) return Front_found;
			 m_temp=(alfa[ii][col_z]-alfa[ii-1][col_z]);
			 if(m_temp>6.0e-15)
				{
			     a_temp2=(0.66-alfa[ii-1][col_z])/m_temp+(double)ii-1.0;
				 a_temp2*=Inc_s;
			     pos_r=a_temp2;
				 Begin_found=true;
				}
			}
		 ii++;
		}while(ii<n_r_sample && Begin_found==false);
	 if(ii>=n_r_sample) return Front_found; 
	 do
		{
		 if(alfa[ii][col_z]>0.99) 
			{
			 m_temp=(alfa[ii][col_z]-alfa[ii-1][col_z]);
			 if(m_temp<6e-15) a_temp1=ii;
			 else a_temp1=(0.99-alfa[ii-1][col_z])/m_temp+(double)ii-1.0;
			 a_temp1*=Inc_s;
			 Front_found=true;
			 width_r=a_temp1-a_temp2;
			}
		 ii++;
		}while(ii<n_r_sample && Front_found==false);
	}
if(Front_found)
 	{
	 if(Com_Front_Created==false)
		{
		 Com_Front_Created=true;
		 puts("\n\n\tWarning!!!!!\n\n\tA combustion front has been created!!!\n\n\n");
		}
	}
 return Front_found;
}
// searches in the axial direction in column col_r
// Prop_Down=true searches a front propagating downwards, false upwards
bool Sample::RegPosFront_z(FILE *PF,double **alfa,bool Prop_Down,double &pos_z,double &width_z)
{
 int jj;
 double a_temp1,a_temp2,m_temp;
 bool Front_found,Begin_found;

 pos_z=0.0;
 width_z=0.0;
 Front_found=false;
 Begin_found=false;
 if(Prop_Down)  // front in axial direction propagating downwards
	{
	 jj=n_z_sample-2;
	 do
		{
		 if(alfa[col_r][jj]>0.66) 
			{
			 if(jj==n_z_sample-1) return Front_found;
			 m_temp=(alfa[col_r][jj]-alfa[col_r][jj+1]);
			 if(m_temp>6.0e-15)
				{
				 a_temp1=-(0.66-alfa[col_r][jj+1])/m_temp+(double)jj+1.0;
				 a_temp1*=Inc_s;
				 Begin_found=true;
				 pos_z=a_temp1;
				}
			}
		 jj--;
		}while(jj>=0 && Begin_found==false);
	 if(jj<1) return Front_found;
	 do
		{
		 if(alfa[col_r][jj]>0.99) 
			{
			 m_temp=(alfa[col_r][jj]-alfa[col_r][jj+1]);
			 if(m_temp>6.0e-15) a_temp2=jj;
			 else a_temp2=-(0.99-alfa[col_r][jj+1])/m_temp+(double)jj+1.0;
			 a_temp2*=Inc_s;
			 Front_found=true;
			 width_z=a_temp1-a_temp2;
			}
		 jj--;
		}while(jj>=0 && Front_found==false);
	}
 else 
	{// front eix y que és propaga cap amunt
	 jj=0;
	 do
		{
		 if(alfa[col_r][jj]>0.66) 
			{
			 if(jj==0) return Front_found;
			 m_temp=(alfa[col_r][jj]-alfa[col_r][jj-1]);
			 if(m_temp>6.0e-15)
				{
			     a_temp2=(0.66-alfa[col_r][jj-1])/m_temp+(double)jj-1.0;
				 a_temp2*=Inc_s;
				 Begin_found=true;
				 pos_z=a_temp2;
				}
			}
		 jj++;
		}while(jj<n_z_sample && Begin_found==false);
	 if(jj>=n_z_sample) return Front_found;
	 do
		{
		 if(alfa[col_r][jj]>0.99) 
			{
			 m_temp=(alfa[col_r][jj]-alfa[col_r][jj-1]);
			 if(m_temp<6e-15) a_temp1=jj;
			 else a_temp1=(0.99-alfa[col_r][jj-1])/m_temp+(double)jj-1.0;
			 a_temp1*=Inc_s;
			 Front_found=true;
			 width_z=a_temp1-a_temp2;
			}
		 jj++;
		}while(jj<n_z_sample && Front_found==false);
	}
 if(Front_found)
 	{
	 if(Com_Front_Created==false)
		{
		 Com_Front_Created=true;
		 puts("\n\n\tWarning!!!!!\n\n\tA combustion front has been created!!!\n\n\n");
		}
	}
 return Front_found;
}
