#include "udf.h"
#include "mem.h"
#include "math.h"
#include "unsteady.h"
#include "pdf_props.h"
#include "sg_udms.h"

real 	k0p, EAp, ap, bp, k0l, EAl, al, bl, kp, kl, EAp_21, EAp_23, EAp_30, EAp_100, EAl_21, EAl_23, EAl_30, t1, t2,
		y_H2, y_O2, y_O, y_H, y_OH, y_H2O, y_HO2, y_CO, y_CO2, y_CH4, y_CH3, y_CH2O, y_HCO, y_CH3O, y_H2O2, y_N2, y_N,
		c_O2, c_H2O, c_CO2, c_N2, c_H2,	c_O, c_H, c_OH,	c_HO2, c_CO, c_CH4, c_CH3, c_CH2O, c_HCO, c_CH3O, c_H2O2, c_N, nu_H2O, Temp;

real react1, react2, react3;
int Zone_ID;
real area[ND_ND];
real area_surf;
real area_tot;

real M_H2 = 2.01588;
real M_O2 = 31.9988;
real M_O = 15.9994;
real M_H = 1.00794;
real M_OH = 17.00734;
real M_H2O = 18.01528;
real M_HO2 = 33.00674;
real M_CO = 28.0101;
real M_CO2 = 44.0095;
real M_CH4 = 16.04246;
real M_CH3 = 15.03272;
real M_CH2O = 30.02568;
real M_HCO = 29.01804;
real M_CH3O = 31.03392;
real M_H2O2 = 34.01468;
real M_N = 14.0067;
real M_N2 = 28.0134;


/* volume farction at known points */
real nu_H2O_21 = 0.171;									/* [-]*/
real nu_H2O_23 = 0.188;									/* [-]*/
real nu_H2O_30 = 0.243;									/* [-]*/
real nu_H2O_100 = 0.655;								/* [-]*/

/* Prefactors */
real k0p_21 = 989567168.3643713;						/*Prefactor of Equation*/

/* activation energy */							
real EAp_21 = 267533.7175778198;
	
/* Calculated Exponents for local atmosphere influence */
real ap_21 = 0.172724946026483;											/* [-]*/
real bp_21 = 0.787730619687136;											/* [-]*/

/* Vorerst unberücksichtigt*/
real t_ges = 9000;										/* [s]*/



DEFINE_EXECUTE_ON_LOADING(rename_UDFs, libname)
{
Set_User_Memory_Name(0,"Spec. mass gain Billet 1");
Set_User_Memory_Name(1,"Delta M Billet 1");

Set_User_Memory_Name(3, "Spec. mass gain Billet 2");
Set_User_Memory_Name(2, "Delta M Billet 2");

Set_User_Memory_Name(5, "Spec. mass gain Billet 3");
Set_User_Memory_Name(4, "Delta M Billet 3");

Set_User_Memory_Name(7, "Spec. mass gain Billet 4");
Set_User_Memory_Name(6, "Delta M Billet 4");
}
	

void scaleformation1(Zone_ID)
{
	Domain *d; /* Get domain */
	face_t f;
	Thread *f_thread, *t0, *c_thread; 
	cell_t c1, c0;

	area_surf = 0.;		/*initialize the surface area*/
	area_tot = 0.;		/*initialize the surface total area*/
	Temp = 0.;			/*initialize the surface temperature*/

	k0p = 0.;			/*initialize parabolic factorn in Arrhenius equation*/
	EAp = 0.;			/*initialize parabolic activation energy in Arrhenius equation*/
	ap = 0.;			/*initialize parabolic reaction exponent for O2*/
	bp = 0.;			/*initialize parabolic reaction exponent for H2O*/

	k0l = 0.;			/*initialize linear factorn in Arrhenius equation*/
	EAl = 0.;			/*initialize linear activation energy in Arrhenius equation*/
	al = 0.;			/*initialize linear reaction exponent for O2*/
	bl = 0.;			/*initialize linear reaction exponent for H2O*/

	react1 = 0.;		/*initialize reaction rate at current time */
	react2 = 0.;		/*initialize reaction rate at previous time */


	d = Get_Domain(1); 
	f_thread = Lookup_Thread(d, Zone_ID);
	
	t0 = THREAD_T0(f_thread);
	c_thread =  THREAD_T1(f_thread);

	t1 = CURRENT_TIME;									/* [s]*/
	t2 = PREVIOUS_TIME;									/* [s]*/

begin_f_loop(f, f_thread)
{
	F_AREA(area, f, f_thread);	/* get cell face */
	area_surf = NV_MAG(area);	/* get size of the cell face */
	area_tot += area_surf;		/* sum of faces */
}
end_f_loop(f, f_thread)

begin_f_loop(f,f_thread)
{ 

	c0 = F_C0(f, f_thread);
	c1 = F_C1(f, f_thread);

	/*According GRIMECH3.0*/
	y_H2 = Pdf_Yi(c0,t0,0);
	y_O2 = Pdf_Yi(c0,t0,3);
	y_O = Pdf_Yi(c0,t0,2);
	y_H = Pdf_Yi(c0,t0,1);
	y_OH = Pdf_Yi(c0,t0,4);
	y_H2O = Pdf_Yi(c0,t0,5);
	y_HO2 = Pdf_Yi(c0,t0,6);
	y_CO = Pdf_Yi(c0,t0,14);
	y_CO2 = Pdf_Yi(c0,t0,15);
	y_CH4 = Pdf_Yi(c0,t0,13);
	y_CH3 = Pdf_Yi(c0,t0,12);
	y_CH2O = Pdf_Yi(c0,t0,17);
	y_HCO = Pdf_Yi(c0,t0,16);
	y_CH3O = Pdf_Yi(c0,t0,19);
	y_H2O2 = Pdf_Yi(c0,t0,7);
	y_N2 = Pdf_Yi(c0,t0,47);
	y_N = Pdf_Yi(c0,t0,34);


	c_O2 = C_R(c0,t0)*1000*y_O2/M_O2; 
	c_H2O = C_R(c0,t0)*1000*y_H2O/M_H2O;
	c_CO2 = C_R(c0,t0)*1000*y_CO2/M_CO2;
	c_N2 = C_R(c0,t0)*1000*y_N2/M_N2;
	c_H2 = C_R(c0,t0)*1000*y_H2/M_H2;
	c_O = C_R(c0,t0)*1000*y_O/M_O;
	c_H = C_R(c0,t0)*1000*y_H/M_H;
	c_OH = C_R(c0,t0)*1000*y_OH/M_OH;
	c_HO2 = C_R(c0,t0)*1000*y_HO2/M_HO2;
	c_CO = C_R(c0,t0)*1000*y_CO/M_CO;
	c_CH4 = C_R(c0,t0)*1000*y_CH4/M_CH4;
	c_CH3 = C_R(c0,t0)*1000*y_CH3/M_CH3;
	c_CH2O = C_R(c0,t0)*1000*y_CH2O/M_CH2O;
	c_HCO = C_R(c0,t0)*1000*y_HCO/M_HCO;
	c_CH3O = C_R(c0,t0)*1000*y_CH3O/M_CH3O;
	c_H2O2 = C_R(c0,t0)*1000*y_H2O2/M_H2O2;
	c_N = C_R(c0,t0)*1000*y_N/M_N;

	nu_H2O = c_H2O / (c_O2+c_H2O+c_CO2+c_N2+c_H2+c_O+c_H+c_OH+c_HO2+c_CO+c_CH4+c_CH3+c_CH2O+c_HCO+c_CH3O+c_H2O2+c_N); 
	Temp = C_T(c1,c_thread); /* Temp  Solid */

if (Temp < 873.15) /* 1.4307 starts to form scale @ 600°C */
{
		k0l = 0.;	
		EAl = 0.;
		al = 0.;
		bl = 0.;

		k0p = 0.;
		EAp = 0.;
		ap = 0.;
		bp = 0.;
}

if (Temp > 873.15) /* Start to calc scale formation */
{
	k0l = 0.; /*Fall NGAIR */
	EAl = 0.;
	al = 0.;
	bl = 0.;

	k0p = k0p_21;
	EAp = EAp_21;
	ap = ap_21;
	bp = bp_21;
}



	react1 = pow(k0p * exp(-EAp / (UNIVERSAL_GAS_CONSTANT / 1000 * Temp)) * pow(c_O2, ap) * pow(c_H2O, bp) * t1 / t_ges, (1 / 1.7)); /* Measurements showed that its not a pure parabolic behavior */
	react2 = pow(k0p * exp(-EAp / (UNIVERSAL_GAS_CONSTANT / 1000 * Temp)) * pow(c_O2, ap) * pow(c_H2O, bp) * t2 / t_ges, (1 / 1.7));

	/*Message("c_O2 Billet 1 %f\n", c_O2);*/
	/*Message("c_H2O Billet 1 %f\n",c_H2O);*/


	C_UDMI(c0,t0,1) = react1-react2; /* delta m per m^2*/
	C_UDMI(c0,t0,0) += C_UDMI(c0,t0,1); /* Mass gain per m^2*/ 	
}
end_f_loop(f,t) 

}

DEFINE_EXECUTE_AT_END(verzunderung1)
{

	/* Aufruf der Funktion scaleformation für die Zone_ID's, an denen die Verzunderung berechnet werden soll. Beliebig oft aufrufbar! */
	/* Immer die Zone_ID auswaehlen, die in den Gasraum zeigt! */

	scaleformation1(34); /* Oberfläche Billet 1*/
}
/*----------------------------------------------------------BILLET 2 ---------------------------------------------------------------------------------------*/
void scaleformation2(Zone_ID)
{
	Domain* d; /* Get domain */
	face_t f;
	Thread* f_thread, * t0, * c_thread;
	cell_t c1, c0;

	area_surf = 0.;		/*initialize the surface area*/
	area_tot = 0.;		/*initialize the surface total area*/
	Temp = 0.;			/*initialize the surface temperature*/

	k0p = 0.;			/*initialize parabolic factorn in Arrhenius equation*/
	EAp = 0.;			/*initialize parabolic activation energy in Arrhenius equation*/
	ap = 0.;			/*initialize parabolic reaction exponent for O2*/
	bp = 0.;			/*initialize parabolic reaction exponent for H2O*/

	k0l = 0.;			/*initialize linear factorn in Arrhenius equation*/
	EAl = 0.;			/*initialize linear activation energy in Arrhenius equation*/
	al = 0.;			/*initialize linear reaction exponent for O2*/
	bl = 0.;			/*initialize linear reaction exponent for H2O*/

	react1 = 0.;		/*initialize reaction rate at current time */
	react2 = 0.;		/*initialize reaction rate at previous time */


	d = Get_Domain(1);
	f_thread = Lookup_Thread(d, Zone_ID);

	t0 = THREAD_T0(f_thread);
	c_thread = THREAD_T1(f_thread);

	t1 = CURRENT_TIME;									/* [s]*/
	t2 = PREVIOUS_TIME;									/* [s]*/

	begin_f_loop(f, f_thread)
	{
		F_AREA(area, f, f_thread);	/* get cell face */
		area_surf = NV_MAG(area);	/* get size of the cell face */
		area_tot += area_surf;		/* sum of faces */
	}
	end_f_loop(f, f_thread)

		begin_f_loop(f, f_thread)
	{

		c0 = F_C0(f, f_thread);
		c1 = F_C1(f, f_thread);

		/*According GRIMECH3.0*/
		y_H2 = Pdf_Yi(c0, t0, 0);
		y_O2 = Pdf_Yi(c0, t0, 3);
		y_O = Pdf_Yi(c0, t0, 2);
		y_H = Pdf_Yi(c0, t0, 1);
		y_OH = Pdf_Yi(c0, t0, 4);
		y_H2O = Pdf_Yi(c0, t0, 5);
		y_HO2 = Pdf_Yi(c0, t0, 6);
		y_CO = Pdf_Yi(c0, t0, 14);
		y_CO2 = Pdf_Yi(c0, t0, 15);
		y_CH4 = Pdf_Yi(c0, t0, 13);
		y_CH3 = Pdf_Yi(c0, t0, 12);
		y_CH2O = Pdf_Yi(c0, t0, 17);
		y_HCO = Pdf_Yi(c0, t0, 16);
		y_CH3O = Pdf_Yi(c0, t0, 19);
		y_H2O2 = Pdf_Yi(c0, t0, 7);
		y_N2 = Pdf_Yi(c0, t0, 47);
		y_N = Pdf_Yi(c0, t0, 34);


		c_O2 = C_R(c0, t0) * 1000 * y_O2 / M_O2;
		c_H2O = C_R(c0, t0) * 1000 * y_H2O / M_H2O;
		c_CO2 = C_R(c0, t0) * 1000 * y_CO2 / M_CO2;
		c_N2 = C_R(c0, t0) * 1000 * y_N2 / M_N2;
		c_H2 = C_R(c0, t0) * 1000 * y_H2 / M_H2;
		c_O = C_R(c0, t0) * 1000 * y_O / M_O;
		c_H = C_R(c0, t0) * 1000 * y_H / M_H;
		c_OH = C_R(c0, t0) * 1000 * y_OH / M_OH;
		c_HO2 = C_R(c0, t0) * 1000 * y_HO2 / M_HO2;
		c_CO = C_R(c0, t0) * 1000 * y_CO / M_CO;
		c_CH4 = C_R(c0, t0) * 1000 * y_CH4 / M_CH4;
		c_CH3 = C_R(c0, t0) * 1000 * y_CH3 / M_CH3;
		c_CH2O = C_R(c0, t0) * 1000 * y_CH2O / M_CH2O;
		c_HCO = C_R(c0, t0) * 1000 * y_HCO / M_HCO;
		c_CH3O = C_R(c0, t0) * 1000 * y_CH3O / M_CH3O;
		c_H2O2 = C_R(c0, t0) * 1000 * y_H2O2 / M_H2O2;
		c_N = C_R(c0, t0) * 1000 * y_N / M_N;

		nu_H2O = c_H2O / (c_O2 + c_H2O + c_CO2 + c_N2 + c_H2 + c_O + c_H + c_OH + c_HO2 + c_CO + c_CH4 + c_CH3 + c_CH2O + c_HCO + c_CH3O + c_H2O2 + c_N);
		Temp = C_T(c1, c_thread); /* Temp  Solid */

		if (Temp < 873.15) /* 1.4307 starts to form scale @ 600°C */
		{
			k0l = 0.;
			EAl = 0.;
			al = 0.;
			bl = 0.;

			k0p = 0.;
			EAp = 0.;
			ap = 0.;
			bp = 0.;
		}

		if (Temp > 873.15) /* Start to calc scale formation */
		{
			k0l = 0.; /*Fall NGAIR */
			EAl = 0.;
			al = 0.;
			bl = 0.;

			k0p = k0p_21;
			EAp = EAp_21;
			ap = ap_21;
			bp = bp_21;
		}



		react1 = pow(k0p * exp(-EAp / (UNIVERSAL_GAS_CONSTANT / 1000 * Temp)) * pow(c_O2, ap) * pow(c_H2O, bp) * t1 / t_ges, (1 / 1.7)); /* Measurements showed that its not a pure parabolic behavior */
		react2 = pow(k0p * exp(-EAp / (UNIVERSAL_GAS_CONSTANT / 1000 * Temp)) * pow(c_O2, ap) * pow(c_H2O, bp) * t2 / t_ges, (1 / 1.7));

		/*Message("c_O2 Billet 1 %f\n", c_O2);*/
		/*Message("c_H2O Billet 1 %f\n",c_H2O);*/


		C_UDMI(c0, t0, 2) = react1 - react2; /* delta m per m^2*/
		C_UDMI(c0, t0, 3) += C_UDMI(c0, t0, 2); /* Mass gain per m^2*/
	}
	end_f_loop(f, t)

}

DEFINE_EXECUTE_AT_END(verzunderung2)
{

	/* Aufruf der Funktion scaleformation für die Zone_ID's, an denen die Verzunderung berechnet werden soll. Beliebig oft aufrufbar! */
	/* Immer die Zone_ID auswaehlen, die in den Gasraum zeigt! */

	scaleformation2(33); /* Oberfläche Billet 2*/
}

/*----------------------------------------------------------BILLET 3 ---------------------------------------------------------------------------------------*/
void scaleformation3(Zone_ID)
{
	Domain* d; /* Get domain */
	face_t f;
	Thread* f_thread, * t0, * c_thread;
	cell_t c1, c0;

	area_surf = 0.;		/*initialize the surface area*/
	area_tot = 0.;		/*initialize the surface total area*/
	Temp = 0.;			/*initialize the surface temperature*/

	k0p = 0.;			/*initialize parabolic factorn in Arrhenius equation*/
	EAp = 0.;			/*initialize parabolic activation energy in Arrhenius equation*/
	ap = 0.;			/*initialize parabolic reaction exponent for O2*/
	bp = 0.;			/*initialize parabolic reaction exponent for H2O*/

	k0l = 0.;			/*initialize linear factorn in Arrhenius equation*/
	EAl = 0.;			/*initialize linear activation energy in Arrhenius equation*/
	al = 0.;			/*initialize linear reaction exponent for O2*/
	bl = 0.;			/*initialize linear reaction exponent for H2O*/

	react1 = 0.;		/*initialize reaction rate at current time */
	react2 = 0.;		/*initialize reaction rate at previous time */


	d = Get_Domain(1);
	f_thread = Lookup_Thread(d, Zone_ID);

	t0 = THREAD_T0(f_thread);
	c_thread = THREAD_T1(f_thread);

	t1 = CURRENT_TIME;									/* [s]*/
	t2 = PREVIOUS_TIME;									/* [s]*/

	begin_f_loop(f, f_thread)
	{
		F_AREA(area, f, f_thread);	/* get cell face */
		area_surf = NV_MAG(area);	/* get size of the cell face */
		area_tot += area_surf;		/* sum of faces */
	}
	end_f_loop(f, f_thread)

		begin_f_loop(f, f_thread)
	{

		c0 = F_C0(f, f_thread);
		c1 = F_C1(f, f_thread);

		/*According GRIMECH3.0*/
		y_H2 = Pdf_Yi(c0, t0, 0);
		y_O2 = Pdf_Yi(c0, t0, 3);
		y_O = Pdf_Yi(c0, t0, 2);
		y_H = Pdf_Yi(c0, t0, 1);
		y_OH = Pdf_Yi(c0, t0, 4);
		y_H2O = Pdf_Yi(c0, t0, 5);
		y_HO2 = Pdf_Yi(c0, t0, 6);
		y_CO = Pdf_Yi(c0, t0, 14);
		y_CO2 = Pdf_Yi(c0, t0, 15);
		y_CH4 = Pdf_Yi(c0, t0, 13);
		y_CH3 = Pdf_Yi(c0, t0, 12);
		y_CH2O = Pdf_Yi(c0, t0, 17);
		y_HCO = Pdf_Yi(c0, t0, 16);
		y_CH3O = Pdf_Yi(c0, t0, 19);
		y_H2O2 = Pdf_Yi(c0, t0, 7);
		y_N2 = Pdf_Yi(c0, t0, 47);
		y_N = Pdf_Yi(c0, t0, 34);


		c_O2 = C_R(c0, t0) * 1000 * y_O2 / M_O2;
		c_H2O = C_R(c0, t0) * 1000 * y_H2O / M_H2O;
		c_CO2 = C_R(c0, t0) * 1000 * y_CO2 / M_CO2;
		c_N2 = C_R(c0, t0) * 1000 * y_N2 / M_N2;
		c_H2 = C_R(c0, t0) * 1000 * y_H2 / M_H2;
		c_O = C_R(c0, t0) * 1000 * y_O / M_O;
		c_H = C_R(c0, t0) * 1000 * y_H / M_H;
		c_OH = C_R(c0, t0) * 1000 * y_OH / M_OH;
		c_HO2 = C_R(c0, t0) * 1000 * y_HO2 / M_HO2;
		c_CO = C_R(c0, t0) * 1000 * y_CO / M_CO;
		c_CH4 = C_R(c0, t0) * 1000 * y_CH4 / M_CH4;
		c_CH3 = C_R(c0, t0) * 1000 * y_CH3 / M_CH3;
		c_CH2O = C_R(c0, t0) * 1000 * y_CH2O / M_CH2O;
		c_HCO = C_R(c0, t0) * 1000 * y_HCO / M_HCO;
		c_CH3O = C_R(c0, t0) * 1000 * y_CH3O / M_CH3O;
		c_H2O2 = C_R(c0, t0) * 1000 * y_H2O2 / M_H2O2;
		c_N = C_R(c0, t0) * 1000 * y_N / M_N;

		nu_H2O = c_H2O / (c_O2 + c_H2O + c_CO2 + c_N2 + c_H2 + c_O + c_H + c_OH + c_HO2 + c_CO + c_CH4 + c_CH3 + c_CH2O + c_HCO + c_CH3O + c_H2O2 + c_N);
		Temp = C_T(c1, c_thread); /* Temp  Solid */

		if (Temp < 873.15) /* 1.4307 starts to form scale @ 600°C */
		{
			k0l = 0.;
			EAl = 0.;
			al = 0.;
			bl = 0.;

			k0p = 0.;
			EAp = 0.;
			ap = 0.;
			bp = 0.;
		}

		if (Temp > 873.15) /* Start to calc scale formation */
		{
			k0l = 0.; /*Fall NGAIR */
			EAl = 0.;
			al = 0.;
			bl = 0.;

			k0p = k0p_21;
			EAp = EAp_21;
			ap = ap_21;
			bp = bp_21;
		}



		react1 = pow(k0p * exp(-EAp / (UNIVERSAL_GAS_CONSTANT / 1000 * Temp)) * pow(c_O2, ap) * pow(c_H2O, bp) * t1 / t_ges, (1 / 1.7)); /* Measurements showed that its not a pure parabolic behavior */
		react2 = pow(k0p * exp(-EAp / (UNIVERSAL_GAS_CONSTANT / 1000 * Temp)) * pow(c_O2, ap) * pow(c_H2O, bp) * t2 / t_ges, (1 / 1.7));

		/*Message("c_O2 Billet 1 %f\n", c_O2);*/
		/*Message("c_H2O Billet 1 %f\n",c_H2O);*/


		C_UDMI(c0, t0, 4) = react1 - react2; /* delta m per m^2*/
		C_UDMI(c0, t0, 5) += C_UDMI(c0, t0, 4); /* Mass gain per m^2*/
	}
	end_f_loop(f, t)

}

DEFINE_EXECUTE_AT_END(verzunderung3)
{

	/* Aufruf der Funktion scaleformation für die Zone_ID's, an denen die Verzunderung berechnet werden soll. Beliebig oft aufrufbar! */
	/* Immer die Zone_ID auswaehlen, die in den Gasraum zeigt! */

	scaleformation3(32); /* Oberfläche Billet 3*/
}

/*----------------------------------------------------------BILLET 4 ---------------------------------------------------------------------------------------*/
void scaleformation4(Zone_ID)
{
	Domain* d; /* Get domain */
	face_t f;
	Thread* f_thread, * t0, * c_thread;
	cell_t c1, c0;

	area_surf = 0.;		/*initialize the surface area*/
	area_tot = 0.;		/*initialize the surface total area*/
	Temp = 0.;			/*initialize the surface temperature*/

	k0p = 0.;			/*initialize parabolic factorn in Arrhenius equation*/
	EAp = 0.;			/*initialize parabolic activation energy in Arrhenius equation*/
	ap = 0.;			/*initialize parabolic reaction exponent for O2*/
	bp = 0.;			/*initialize parabolic reaction exponent for H2O*/

	k0l = 0.;			/*initialize linear factorn in Arrhenius equation*/
	EAl = 0.;			/*initialize linear activation energy in Arrhenius equation*/
	al = 0.;			/*initialize linear reaction exponent for O2*/
	bl = 0.;			/*initialize linear reaction exponent for H2O*/

	react1 = 0.;		/*initialize reaction rate at current time */
	react2 = 0.;		/*initialize reaction rate at previous time */


	d = Get_Domain(1);
	f_thread = Lookup_Thread(d, Zone_ID);

	t0 = THREAD_T0(f_thread);
	c_thread = THREAD_T1(f_thread);

	t1 = CURRENT_TIME;									/* [s]*/
	t2 = PREVIOUS_TIME;									/* [s]*/

	begin_f_loop(f, f_thread)
	{
		F_AREA(area, f, f_thread);	/* get cell face */
		area_surf = NV_MAG(area);	/* get size of the cell face */
		area_tot += area_surf;		/* sum of faces */
	}
	end_f_loop(f, f_thread)

		begin_f_loop(f, f_thread)
	{

		c0 = F_C0(f, f_thread);
		c1 = F_C1(f, f_thread);

		/*According GRIMECH3.0*/
		y_H2 = Pdf_Yi(c0, t0, 0);
		y_O2 = Pdf_Yi(c0, t0, 3);
		y_O = Pdf_Yi(c0, t0, 2);
		y_H = Pdf_Yi(c0, t0, 1);
		y_OH = Pdf_Yi(c0, t0, 4);
		y_H2O = Pdf_Yi(c0, t0, 5);
		y_HO2 = Pdf_Yi(c0, t0, 6);
		y_CO = Pdf_Yi(c0, t0, 14);
		y_CO2 = Pdf_Yi(c0, t0, 15);
		y_CH4 = Pdf_Yi(c0, t0, 13);
		y_CH3 = Pdf_Yi(c0, t0, 12);
		y_CH2O = Pdf_Yi(c0, t0, 17);
		y_HCO = Pdf_Yi(c0, t0, 16);
		y_CH3O = Pdf_Yi(c0, t0, 19);
		y_H2O2 = Pdf_Yi(c0, t0, 7);
		y_N2 = Pdf_Yi(c0, t0, 47);
		y_N = Pdf_Yi(c0, t0, 34);


		c_O2 = C_R(c0, t0) * 1000 * y_O2 / M_O2;
		c_H2O = C_R(c0, t0) * 1000 * y_H2O / M_H2O;
		c_CO2 = C_R(c0, t0) * 1000 * y_CO2 / M_CO2;
		c_N2 = C_R(c0, t0) * 1000 * y_N2 / M_N2;
		c_H2 = C_R(c0, t0) * 1000 * y_H2 / M_H2;
		c_O = C_R(c0, t0) * 1000 * y_O / M_O;
		c_H = C_R(c0, t0) * 1000 * y_H / M_H;
		c_OH = C_R(c0, t0) * 1000 * y_OH / M_OH;
		c_HO2 = C_R(c0, t0) * 1000 * y_HO2 / M_HO2;
		c_CO = C_R(c0, t0) * 1000 * y_CO / M_CO;
		c_CH4 = C_R(c0, t0) * 1000 * y_CH4 / M_CH4;
		c_CH3 = C_R(c0, t0) * 1000 * y_CH3 / M_CH3;
		c_CH2O = C_R(c0, t0) * 1000 * y_CH2O / M_CH2O;
		c_HCO = C_R(c0, t0) * 1000 * y_HCO / M_HCO;
		c_CH3O = C_R(c0, t0) * 1000 * y_CH3O / M_CH3O;
		c_H2O2 = C_R(c0, t0) * 1000 * y_H2O2 / M_H2O2;
		c_N = C_R(c0, t0) * 1000 * y_N / M_N;

		nu_H2O = c_H2O / (c_O2 + c_H2O + c_CO2 + c_N2 + c_H2 + c_O + c_H + c_OH + c_HO2 + c_CO + c_CH4 + c_CH3 + c_CH2O + c_HCO + c_CH3O + c_H2O2 + c_N);
		Temp = C_T(c1, c_thread); /* Temp  Solid */

		if (Temp < 873.15) /* 1.4307 starts to form scale @ 600°C */
		{
			k0l = 0.;
			EAl = 0.;
			al = 0.;
			bl = 0.;

			k0p = 0.;
			EAp = 0.;
			ap = 0.;
			bp = 0.;
		}

		if (Temp > 873.15) /* Start to calc scale formation */
		{
			k0l = 0.; /*Fall NGAIR */
			EAl = 0.;
			al = 0.;
			bl = 0.;

			k0p = k0p_21;
			EAp = EAp_21;
			ap = ap_21;
			bp = bp_21;
		}



		react1 = pow(k0p * exp(-EAp / (UNIVERSAL_GAS_CONSTANT / 1000 * Temp)) * pow(c_O2, ap) * pow(c_H2O, bp) * t1 / t_ges, (1 / 1.7)); /* Measurements showed that its not a pure parabolic behavior */
		react2 = pow(k0p * exp(-EAp / (UNIVERSAL_GAS_CONSTANT / 1000 * Temp)) * pow(c_O2, ap) * pow(c_H2O, bp) * t2 / t_ges, (1 / 1.7));

		/*Message("c_O2 Billet 1 %f\n", c_O2);*/
		/*Message("c_H2O Billet 1 %f\n",c_H2O);*/


		C_UDMI(c0, t0, 6) = react1 - react2; /* delta m per m^2*/
		C_UDMI(c0, t0, 7) += C_UDMI(c0, t0, 6); /* Mass gain per m^2*/
	}
	end_f_loop(f, t)

}

DEFINE_EXECUTE_AT_END(verzunderung4)
{

	/* Aufruf der Funktion scaleformation für die Zone_ID's, an denen die Verzunderung berechnet werden soll. Beliebig oft aufrufbar! */
	/* Immer die Zone_ID auswaehlen, die in den Gasraum zeigt! */

	scaleformation4(35); /* Oberfläche Billet 4*/
}