format compact

% Lämpötilan alue

T1 = 1650;

T2 = 1700;

% Taulukoidut arvot

% Massa (kg) [CO2 H2O N2 O2 SO2]

m = [1.13097 0.405 3.01302 0.65348 0.00036];

% Entalpia (MJ/kg) [O2 N2 CO2 H2O SO2]

h1 = [1.441 1.559 1.602 3.071 1.108];   % @ T1

h2 = [1.499 1.622 1.67 3.342 1.153];   % @ T2

% Energiatase yhtälö (vas. ja oik. puoli):

eq1 = 8.893;   % Valmiiksi laskettu

eq2 = 0;        % = sum(m_i*h_i)

% eq2 pitäisi olla lopulta yhtäsuuri kuin eq1

% Interpolaattori

i=0;

while (eq2 ~= eq1)

   i=i+1;

   avg_T = (T1+T2)/2; % Interpoloidaan T_ad

   avg_h = (h1+h2)/2; % Interpoloidaan entalpioille keskiarvot

   avg_cp = (cp1+cp2)/2;

   eq2 = m(1)*avg_h(3)+m(2)*avg_h(4)+m(3)*avg_h(2)+m(4)*avg_h(1)+m(5)*avg_h(5);

   if eq2 < eq1

       T1 = avg_T;

       h1 = avg_h;

       cp1 = avg_cp;

   else

       T2 = avg_T;

       h2 = avg_h;

       cp2 = avg_cp;

   end

   % Varmuuden vuoksi lopetus tietyn iteroinnin päästä.

   % Jos näin käy, luultavasti tapahtunut käyttäjävirhe

   % lähtöarvojen asetuksessa.

   if i>1000

       disp("TARKISTA ASETETUT ARVOT.")

       break

   end

end

% Tulosta tulokset

disp("Interpolaation kierrosten määrä: "+i)

disp("Energiatase: "+eq1+" = "+eq2)

disp("Entalpiat [O2 N2 CO2 H2O SO2] (MJ/kg): "), disp(avg_h)

disp("Lämpökapasiteetit (kJ/kgK): "), disp(avg_cp)

disp("Lämpötila T_ad: "+avg_T)