// This script is part of a calculation of the reaction rate of the chemical reaction with coal and oxygen. The calculation is done by iterating over steps. First the reaction rate of each component is calculated on each step. Using this rate the script calculates the new concentration of each component. Also the heat release of the chemical reaction is calculated.

// The reaction rate is influenced by the temperature. Due to the heat release of the reaction, the iteration needs to be done to approximate more accurate results.

//------------------------
// Input Data
//------------------------
clear;
"Beginn"
T0      =  300;         // temperature
ER      =  12682;
h       =  -2.2035D7;
cw      =  610400;
d       =  1;
cB0 =  560;
cO0     =  0.06;
cS0     =  0;
cG0     =  0;
MB      =  3948.5265;
MO      =  31.9988;
MS      =  2192.5792;
MG      =  30.7989406666667;
nuB     =  -1;
nuO     =  -89.5;
nuS     =  1;
nuG     =  150;
K1      =  2.3574D7;
K2      =  3.929D8;
abbruchT = 730;
abbruchC = 0.0001;

//----------------------------------------------------------
// Target File
//----------------------------------------------------------
pathname = get_absolute_file_path('b.sce');
getshortpathname(pathname);
tarFile2 = mopen(pathname+'b2.csv','w');

//----------------------------------------------------------
// Functions
//----------------------------------------------------------
// Reaction Rates
function RB2 = calcReactRateB2(cB2a,cO2a,T2a)
  RB2 = -cB2a * cO2a * K2 * exp(-(ER/T2a)) *60*60*24;
endfunction
function RO2 = calcReactRateO2(RB2)
  RO2 = (nuO/nuB)*(MO/MB)*RB2;
endfunction
function RS2 = calcReactRateS2(RB2)
  RS2 = (nuS/nuB)*(MS/MB)*RB2;
endfunction
function RG2 = calcReactRateG2(RB2)
  RG2 = (nuG/nuB)*(MG/MB)*RB2;
endfunction
// Concentrations
function cB2 = calcConB2(cB2a,RB2)
  cB2 = cB2a+RB2*d;
endfunction
function cO2 = calcConO2(cO2a,RO2)
  cO2 = cO2a+RO2*d;
endfunction
function cS2 = calcConS2(cS2a,RS2)
  cS2 = cS2a+RS2*d;
endfunction
function cG2 = calcConG2(cG2a,RG2)
  cG2 = cG2a+RG2*d;
endfunction
// Heat Release
function q2 = calcHeatRelease2(RB2)
  q2 = RB2*h*d;
endfunction
function Tzu2 = calcTempInc2(q2)
  Tzu2 = q2/cw;
endfunction
function T2 = calcTemp2(T2a,Tzu2)
  T2 = T2a + Tzu2;
endfunction
// Balance
function z2 = calcConBalance2(cB2,cO2,cS2,cG2)
  z2 = cB2 + cO2 + cS2 + cG2;
endfunction
function v2 = calcRateBalance2(RB2,RO2,RS2,RG2)
  v2 = RB2 + RO2 + RS2 + RG2;
endfunction

//----------------------------------------------------------
// Values in step 0
//----------------------------------------------------------
t=0;
// definieren
RB2 = 0; RO2 = 0; RS2 = 0; RG2=0;
cB2 = cB0; cO2 = cO0; cS2 = cS0; cG2 = cG0;
q2 = 0; Tzu2 = 0; T2 = T0;
z2 = calcConBalance2(cB2,cO2,cS2,cG2);
v2 = calcRateBalance2(RB2,RO2,RS2,RG2);
// schreiben
mfprintf(tarFile2,'%f\t %30.29f\t %30.29f\t %30.29f\t %30.29f\t %30.29f\t %30.29f\t %30.29f\t %30.29f\t %30.29f\t %30.29f\t %30.29f\t %30.29f\t %30.29f\n',t,RB2,RO2,RS2,RG2,cB2,cO2,cS2,cG2,q2,Tzu2,T2,z2,v2)

// iteration left out to minimize the example

//----------------------------------------------------------
// close file
//----------------------------------------------------------
mclose(tarFile2);
"Fertig!"