Revision history [back]

Hello all, This was a homework assignment for a class on thermodynamics. The MATLAB code I will post below gives the correct answer while my sagemath code that was in a worksheet gives completely different and ludicrous answers and I can't figure out what the difference is between the two. I ran the sage worksheet on cloud.sagemath.com, which is awesome by the way. I am pretty new to Sage but can definitely see the power and ease of use that it brings. The MATLAB code is split into two files, one of which is a function called by the second file.

function [ F ] = ODE( t , y )
% This function is to be used in conjunction with file Project1.m
% It uses a reduction of order to turn a 2nd order ODE into a system
% of first order ODE's.  All y(1) and y(2) are the slopes of the velocity
% the acceleration curves respectively.  Z is the height of the object,
% wn is the natural frequency of the system.  These two items are
% specifically noted here because they are functions of time in this
% system. F(1) stores values of the position, and F(2) stores the velocity.

global g m zeta V

z = V*t;

wn = sqrt((600000/(24-z)/m));

Forcing = 2*zeta*wn*(V-y(2)) + wn^2*(z-y(1));

if Forcing < 0      % check point for forcing function (no pushing of rope)
    L = 0;
else
    L = Forcing;
end

F = zeros(2,1);

F(1) = y(2);
F(2) = L - g;
%F(2) = g + 2*zeta*wn*V + wn^2*z - 2*zeta*wn*y(2) - wn^2*y(1);

end

End of first MATLAB file, ODE.m Sorry having trouble with the formatting.

% Define Constants
global g m zeta V
g    = 32.2;    % ft/s^2
%m    = 170;    % lbm
m    = 5.2795;  % slugs
zeta = 0.1;     % damping ratio
V    = 2;       % ft/s

%% Drum Calculations
GR      = 75;           % Gear Ratio
RPSin   = 1275/60;      % Input RPM from Motor (rev/sec with /60)
RPSout  = RPSin/GR;     % Output rev/sec for drum
w_motor = RPSin*2*pi;   % Omega of Motor (Rad/sec)
f_drum = RPSout;        % Natural frequency of drum (Hz)
w_drum = RPSout*2*pi;   % Natural rotational freq (rad/s)
drad    = V/w_drum;     % Drum radius (ft)  
ddiam   = drad*2;       % Drum diameter (ft)

%% ODE Solution
% Example 1
y=[0 0]; % These are the initial conditions of the system.  It begins at
% z=0 with velocity=0.  The initial jerk of the cable sets the system in 
% motion. 

tspan=[0 1]; % This is the time span that the solution will be analyzed.

options=odeset('RelTol',1e-4,'AbsTol',[1e-4, 1e-4]);
% Options defines the error that the ode45 function uses to determine the 
% time step.

[T,YV]=ode45(@ODE,tspan,y,options);
% ode45 is a 4th order Runge Kutta ode solver. It also uses a 5th order
% Runge Kutta to correct the time step and reduce error.  T is the time
% Step of the solver and YV is an array indicating the displacement
% and velocity.  Note that ode45 is calling an ODE that is in a different
% script.

Y=YV(:,1); % This is the displacements at time t of the solution.
v=YV(:,2); % This is the velocities at time t of the solution. 

% Force Time Function
for i = 1:length(T)-1
    a(i) = (v(i+1)-v(i))/(T(i+1)-T(i));
    F(i) = m*a(i)+m*g;
    Torque(i) = F(i)*drad;
    Power(i) = Torque(i) * w_motor;
end

a(length(T)) = a(length(T)-1);
F(length(T)) = F(length(T)-1);
Torque(length(T)) = Torque(length(T)-1);
Power(length(T)) = Power(length(T)-1);

%% Power Calculations
Fmax = max(F);              
Tmax = max(Torque);
Pmax = Tmax * w_drum;
Tmax_motor = Tmax/GR;
Pmax_motor = Tmax_motor*w_motor;

%% Plot results
figure(1);clf(1);
subplot(2,1,1); % Subplot has floating axes
plot(T,Y);
xlabel('Time (sec)');
ylabel('displacement (ft)');
subplot(2,1,2);
plot(T,v);
xlabel('Time (sec)');
ylabel('veolcity (ft/s)');

figure(2);clf(2);
subplot(2,1,1);
plot(T,a);
ylabel('Acceleration (ft/s^2)');
xlabel('Time (sec)');
subplot(2,1,2);
plot(T,F);
xlabel('Time (sec)');
ylabel('Force (lbf)');

figure(3);clf(3);
plot(T,Power)
xlabel('Time (sec)');
ylabel('Power (lbf*ft/sec)');


%% Print Values of Interest
fprintf('Max Force in The Rope (N) %4.2f\n',Fmax);
fprintf('Required motor power (HP) %4.2f\n',Pmax*0.001818);
fprintf('Required Drum Diameter (in) %4.2f\n',ddiam);

Here is the public access link to the sagemath worksheet on the cloud.sagemath.com website, which has embedded explanatory text as to what is happening:

I can't post links yet so add preview.tinyurl.com/ before the following: SageHWhelp

Thank you in advance for the help! This is driving me crazy!

Hello all, all,

This was a homework assignment for a class on thermodynamics. The MATLAB code I will post below gives the correct answer while my sagemath code that was in a worksheet gives completely different and ludicrous answers and I can't figure out what the difference is between the two. I ran the sage worksheet on cloud.sagemath.com, http://cloud.sagemath.com/, which is awesome by the way. I am pretty new to Sage but can definitely see the power and ease of use that it brings. The MATLAB code is split into two files, one of which is a function called by the second file.

function [ F ] = ODE( t , y )
% This function is to be used in conjunction with file Project1.m
% It uses a reduction of order to turn a 2nd order ODE into a system
% of first order ODE's.  All y(1) and y(2) are the slopes of the velocity
% the acceleration curves respectively.  Z is the height of the object,
% wn is the natural frequency of the system.  These two items are
% specifically noted here because they are functions of time in this
% system. F(1) stores values of the position, and F(2) stores the velocity.

global g m zeta V

z = V*t;

wn = sqrt((600000/(24-z)/m));

Forcing = 2*zeta*wn*(V-y(2)) + wn^2*(z-y(1));

if Forcing < 0      % check point for forcing function (no pushing of rope)
    L = 0;
else
    L = Forcing;
end

F = zeros(2,1);

F(1) = y(2);
F(2) = L - g;
%F(2) = g + 2*zeta*wn*V + wn^2*z - 2*zeta*wn*y(2) - wn^2*y(1);

end

End of first MATLAB file, ODE.m Sorry having trouble with the formatting.

% Define Constants
global g m zeta V
g    = 32.2;    % ft/s^2
%m    = 170;    % lbm
m    = 5.2795;  % slugs
zeta = 0.1;     % damping ratio
V    = 2;       % ft/s

%% Drum Calculations
GR      = 75;           % Gear Ratio
RPSin   = 1275/60;      % Input RPM from Motor (rev/sec with /60)
RPSout  = RPSin/GR;     % Output rev/sec for drum
w_motor = RPSin*2*pi;   % Omega of Motor (Rad/sec)
f_drum = RPSout;        % Natural frequency of drum (Hz)
w_drum = RPSout*2*pi;   % Natural rotational freq (rad/s)
drad    = V/w_drum;     % Drum radius (ft)  
ddiam   = drad*2;       % Drum diameter (ft)

%% ODE Solution
% Example 1
y=[0 0]; % These are the initial conditions of the system.  It begins at
% z=0 with velocity=0.  The initial jerk of the cable sets the system in 
% motion. 

tspan=[0 1]; % This is the time span that the solution will be analyzed.

options=odeset('RelTol',1e-4,'AbsTol',[1e-4, 1e-4]);
% Options defines the error that the ode45 function uses to determine the 
% time step.

[T,YV]=ode45(@ODE,tspan,y,options);
% ode45 is a 4th order Runge Kutta ode solver. It also uses a 5th order
% Runge Kutta to correct the time step and reduce error.  T is the time
% Step of the solver and YV is an array indicating the displacement
% and velocity.  Note that ode45 is calling an ODE that is in a different
% script.

Y=YV(:,1); % This is the displacements at time t of the solution.
v=YV(:,2); % This is the velocities at time t of the solution. 

% Force Time Function
for i = 1:length(T)-1
    a(i) = (v(i+1)-v(i))/(T(i+1)-T(i));
    F(i) = m*a(i)+m*g;
    Torque(i) = F(i)*drad;
    Power(i) = Torque(i) * w_motor;
end

a(length(T)) = a(length(T)-1);
F(length(T)) = F(length(T)-1);
Torque(length(T)) = Torque(length(T)-1);
Power(length(T)) = Power(length(T)-1);

%% Power Calculations
Fmax = max(F);              
Tmax = max(Torque);
Pmax = Tmax * w_drum;
Tmax_motor = Tmax/GR;
Pmax_motor = Tmax_motor*w_motor;

%% Plot results
figure(1);clf(1);
subplot(2,1,1); % Subplot has floating axes
plot(T,Y);
xlabel('Time (sec)');
ylabel('displacement (ft)');
subplot(2,1,2);
plot(T,v);
xlabel('Time (sec)');
ylabel('veolcity (ft/s)');

figure(2);clf(2);
subplot(2,1,1);
plot(T,a);
ylabel('Acceleration (ft/s^2)');
xlabel('Time (sec)');
subplot(2,1,2);
plot(T,F);
xlabel('Time (sec)');
ylabel('Force (lbf)');

figure(3);clf(3);
plot(T,Power)
xlabel('Time (sec)');
ylabel('Power (lbf*ft/sec)');


%% Print Values of Interest
fprintf('Max Force in The Rope (N) %4.2f\n',Fmax);
fprintf('Required motor power (HP) %4.2f\n',Pmax*0.001818);
fprintf('Required Drum Diameter (in) %4.2f\n',ddiam);

Here is the public access link to the sagemath worksheet on the cloud.sagemath.com website, which has embedded explanatory text as to what is happening:

Thank you in advance for the help! This is driving me crazy!

function [ F ] = ODE( t , y )
% This function is to be used in conjunction with file Project1.m
% It uses a reduction of order to turn a 2nd order ODE into a system
% of first order ODE's.  All y(1) and y(2) are the slopes of the velocity
% the acceleration curves respectively.  Z is the height of the object,
% wn is the natural frequency of the system.  These two items are
% specifically noted here because they are functions of time in this
% system. F(1) stores values of the position, and F(2) stores the velocity.

global g m zeta V

z = V*t;

wn = sqrt((600000/(24-z)/m));

Forcing = 2*zeta*wn*(V-y(2)) + wn^2*(z-y(1));

if Forcing < 0      % check point for forcing function (no pushing of rope)
    L = 0;
else
    L = Forcing;
end

F = zeros(2,1);

F(1) = y(2);
F(2) = L - g;
%F(2) = g + 2*zeta*wn*V + wn^2*z - 2*zeta*wn*y(2) - wn^2*y(1);

end

% Define Constants
global g m zeta V
g    = 32.2;    % ft/s^2
%m    = 170;    % lbm
m    = 5.2795;  % slugs
zeta = 0.1;     % damping ratio
V    = 2;       % ft/s

%% Drum Calculations
GR      = 75;           % Gear Ratio
RPSin   = 1275/60;      % Input RPM from Motor (rev/sec with /60)
RPSout  = RPSin/GR;     % Output rev/sec for drum
w_motor = RPSin*2*pi;   % Omega of Motor (Rad/sec)
f_drum = RPSout;        % Natural frequency of drum (Hz)
w_drum = RPSout*2*pi;   % Natural rotational freq (rad/s)
drad    = V/w_drum;     % Drum radius (ft)  
ddiam   = drad*2;       % Drum diameter (ft)

%% ODE Solution
% Example 1
y=[0 0]; % These are the initial conditions of the system.  It begins at
% z=0 with velocity=0.  The initial jerk of the cable sets the system in 
% motion. 

tspan=[0 1]; % This is the time span that the solution will be analyzed.

options=odeset('RelTol',1e-4,'AbsTol',[1e-4, 1e-4]);
% Options defines the error that the ode45 function uses to determine the 
% time step.

[T,YV]=ode45(@ODE,tspan,y,options);
% ode45 is a 4th order Runge Kutta ode solver. It also uses a 5th order
% Runge Kutta to correct the time step and reduce error.  T is the time
% Step of the solver and YV is an array indicating the displacement
% and velocity.  Note that ode45 is calling an ODE that is in a different
% script.

Y=YV(:,1); % This is the displacements at time t of the solution.
v=YV(:,2); % This is the velocities at time t of the solution. 

% Force Time Function
for i = 1:length(T)-1
    a(i) = (v(i+1)-v(i))/(T(i+1)-T(i));
    F(i) = m*a(i)+m*g;
    Torque(i) = F(i)*drad;
    Power(i) = Torque(i) * w_motor;
end

a(length(T)) = a(length(T)-1);
F(length(T)) = F(length(T)-1);
Torque(length(T)) = Torque(length(T)-1);
Power(length(T)) = Power(length(T)-1);

%% Power Calculations
Fmax = max(F);              
Tmax = max(Torque);
Pmax = Tmax * w_drum;
Tmax_motor = Tmax/GR;
Pmax_motor = Tmax_motor*w_motor;

%% Plot results
figure(1);clf(1);
subplot(2,1,1); % Subplot has floating axes
plot(T,Y);
xlabel('Time (sec)');
ylabel('displacement (ft)');
subplot(2,1,2);
plot(T,v);
xlabel('Time (sec)');
ylabel('veolcity (ft/s)');

figure(2);clf(2);
subplot(2,1,1);
plot(T,a);
ylabel('Acceleration (ft/s^2)');
xlabel('Time (sec)');
subplot(2,1,2);
plot(T,F);
xlabel('Time (sec)');
ylabel('Force (lbf)');

figure(3);clf(3);
plot(T,Power)
xlabel('Time (sec)');
ylabel('Power (lbf*ft/sec)');


%% Print Values of Interest
fprintf('Max Force in The Rope (N) %4.2f\n',Fmax);
fprintf('Required motor power (HP) %4.2f\n',Pmax*0.001818);
fprintf('Required Drum Diameter (in) %4.2f\n',ddiam);

Here is the public access link to the sagemath worksheet on the cloud.sagemath.com website, which has embedded explanatory text as to what is happening:

https://cloud.sagemath.com/projects/059085c6-509b-45b7-b9aa-f5f07c7e0194/files/EMEC%20425%20Thermal%20Systems/HW3/Thermal%20Systems%20HW%203.sagews

also at

http://preview.tinyurl.com/SageHWhelp

Thank you in advance for the help! This is driving me crazy!

f1(x1,x2,x3)= (yC1^vC1*yD1^vD1) / (yA1^vA1)*(P/P0)^(vC1+vD1-vA1)-K1
f2(x1,x2,x3)= (yC2^vC2) / (yA2^vA2*yB2^vB2)*(P/P0)^(vC2-vB2-vA2)-K2
f3(x1,x2,x3)= (yC3^vC3) / (yA3^vA3*yB3^vB3)*(P/P0)^(vC3-vB3-vA3)-K3
solutions = solve([f1, f2,f3],x1,x2,x3,solution_dict=True)

Here is the public access link to the sagemath worksheet on the cloud.sagemath.com website, which has embedded explanatory text as to what is happening:

https://cloud.sagemath.com/projects/059085c6-509b-45b7-b9aa-f5f07c7e0194/files/EMEC%20425%20Thermal%20Systems/HW3/Thermal%20Systems%20HW%203.sagews

also at

http://preview.tinyurl.com/SageHWhelp

Thank you in advance for the help! This is driving me crazy!

EDIT: I guess I should add some more explanatory text here about the sage worksheet. The answers in that are all correct up until the very last cell that contains the following:

f1(x1,x2,x3)= (yC1^vC1*yD1^vD1) / (yA1^vA1)*(P/P0)^(vC1+vD1-vA1)-K1
f2(x1,x2,x3)= (yC2^vC2) / (yA2^vA2*yB2^vB2)*(P/P0)^(vC2-vB2-vA2)-K2
f3(x1,x2,x3)= (yC3^vC3) / (yA3^vA3*yB3^vB3)*(P/P0)^(vC3-vB3-vA3)-K3
solutions = solve([f1, f2,f3],x1,x2,x3,solution_dict=True)

EDIT again: Wow that is embarrassing, I put up the wrong MATLAB code for comparison. It also appears that I have lost the MATLAB code to compare with. I guess the question is still there, though: Is sage not capable of solving these equations? I can tell you that I compared each one of the final equations f1 through f3 and they are all correct going into the solver as far as values go. It's just that the solver spits back different values than I expect. I had no less than 4 others check the numbers going in, including my professor, although I will say that none of us are very familiar with sage or python.

function HW3_d()
%% Embedded Function

function F = HW3_d(x)


%% Given Conditions
P0=100;     % Ambient pressure (kPa)
T0=298.15;  % Ambient temperature (K)
P=500;      % Products pressure (kPa)
T=1600;     % Products temperature (K)


%% Combustion Reaction: CH8H18 + 37.5*(02 + 3.76*N2) -> 8*CO2 + 9*H2O + 141*N2 + 25*02
nO2=25+x(1)-x(2)-2*x(3); 
nN2=141-x(2)-x(3); 
nH2O=9; 
nCO2=8-2*x(1); 
% constants are initial concentrations


%% Dissociation Reactions 

% Reaction #1: 2CO2 -> 2CO + O2
% Reaction #2: N2 + 02 -> 2NO
% Reaction #3: N2 + 202 -> 2NO2

% Number of Moles
nCO = 2*x(1);
nNO = 2*x(2);
nNO2 = 2*x(3);
ntot = nCO2 + nH2O + nN2 + nO2 + nCO + nNO + nNO2;


%% Equilibrium constants at 1600 K from Table A.11
K1=exp(-21.656); K2=exp(-10.547); K3=exp(-20.126);


%% Mole Fractions
yA1=nCO2/ntot; yB1=0; yC1=nCO/ntot; yD1=nO2/ntot;
yA2=nN2/ntot; yB2=nO2/ntot; yC2=nNO/ntot; yD2=0;
yA3=nN2/ntot; yB3=nO2/ntot; yC3=nNO2/ntot; yD3=0;


%% Stoichiometric Coefficients
vA1=2; vB1=0; vC1=2; vD1=1;
vA2=1; vB2=1; vC2=2; vD2=0;
vA3=1; vB3=2; vC3=2; vD3=0;


%% Equilibrium Constants
F=[(yC1^vC1*yD1^vD1)/(yA1^vA1*yB1^vB1)*(P/P0)^(vC1+vD1-vA1-vB1)-K1;
    (yC2^vC2*yD2^vD2)/(yA2^vA2*yB2^vB2)*(P/P0)^(vC2+vD2-vA2-vB2)-K2;
    (yC3^vC3*yD3^vD3)/(yA3^vA3*yB3^vB3)*(P/P0)^(vC3+vD3-vA3-vB3)-K3];
end

x=fsolve(@(x) HW3_d(x),[0.3;0.3;0.3])


%% Combustion Reaction: CH8H18 + 37.5*(02 + 3.76*N2) -> 8*CO2 + 9*H2O + 141*N2 + 25*02
    % Number of moles
nH2O = 9
nCO2 = 8 - 2*x(1)
nN2 = 141 - x(2) - x(3)
nO2 = 25 + x(1) - x(2) - 2*x(3)


%% Dissociation Reactions 

% Reaction #1: 2CO2 -> 2CO + O2
% Reaction #2: N2 + 02 -> 2NO
% Reaction #3: N2 + 202 -> 2NO2

% Number of Moles
nCO = 2*x(1)
nNO = 2*x(2)
nNO2 = 2*x(3)
ntot = nCO2 + nH2O + nN2 + nO2 + nCO + nNO + nNO2;


%% Mole Fractions
yH2O=nH2O/ntot
yO2=nO2/ntot
yN2=nN2/ntot
yCO2=nCO2/ntot
yCO=nCO/ntot
yNO=nNO/ntot
yNO2=nNO2/ntot
ytot = yH2O + yO2 + yN2 + yCO2 + yCO + yNO + yNO2
end

Here is the public access link to the sagemath worksheet on the cloud.sagemath.com website, which has embedded explanatory text as to what is happening:

https://cloud.sagemath.com/projects/059085c6-509b-45b7-b9aa-f5f07c7e0194/files/EMEC%20425%20Thermal%20Systems/HW3/Thermal%20Systems%20HW%203.sagews

also at

http://preview.tinyurl.com/SageHWhelp

Thank you in advance for the help! This is driving me crazy!

EDIT: I guess I should add some more explanatory text here about the sage worksheet. The answers in that are all correct up until the very last cell that contains the following:

f1(x1,x2,x3)= (yC1^vC1*yD1^vD1) / (yA1^vA1)*(P/P0)^(vC1+vD1-vA1)-K1
f2(x1,x2,x3)= (yC2^vC2) / (yA2^vA2*yB2^vB2)*(P/P0)^(vC2-vB2-vA2)-K2
f3(x1,x2,x3)= (yC3^vC3) / (yA3^vA3*yB3^vB3)*(P/P0)^(vC3-vB3-vA3)-K3
solutions = solve([f1, f2,f3],x1,x2,x3,solution_dict=True)

EDIT again: Wow that is embarrassing, I put up the wrong MATLAB code for comparison. It also appears that I have lost the MATLAB code to compare with. I guess the question is still there, though: Is sage not capable of solving these equations? I can tell you that I compared each one of the final equations f1 through f3 and they are all correct going into the solver as far as values go. It's just that the solver spits back different values than I expect. I had no less than 4 others check the numbers going in, including my professor, although I will say that none of us are very familiar with sage or python.

EDIT again again: Christ Almighty. I found the MATLAB code. I really need to reorganize my files.

function HW3_d()
%% Embedded Function

function F = HW3_d(x)


%% Given Conditions
P0=100;     % Ambient pressure (kPa)
T0=298.15;  % Ambient temperature (K)
P=500;      % Products pressure (kPa)
T=1600;     % Products temperature (K)


%% Combustion Reaction: CH8H18 + 37.5*(02 + 3.76*N2) -> 8*CO2 + 9*H2O + 141*N2 + 25*02
nO2=25+x(1)-x(2)-2*x(3); 
nN2=141-x(2)-x(3); 
nH2O=9; 
nCO2=8-2*x(1); 
% constants are initial concentrations


%% Dissociation Reactions 

% Reaction #1: 2CO2 -> 2CO + O2
% Reaction #2: N2 + 02 -> 2NO
% Reaction #3: N2 + 202 -> 2NO2

% Number of Moles
nCO = 2*x(1);
nNO = 2*x(2);
nNO2 = 2*x(3);
ntot = nCO2 + nH2O + nN2 + nO2 + nCO + nNO + nNO2;


%% Equilibrium constants at 1600 K from Table A.11
K1=exp(-21.656); K2=exp(-10.547); K3=exp(-20.126);


%% Mole Fractions
yA1=nCO2/ntot; yB1=0; yC1=nCO/ntot; yD1=nO2/ntot;
yA2=nN2/ntot; yB2=nO2/ntot; yC2=nNO/ntot; yD2=0;
yA3=nN2/ntot; yB3=nO2/ntot; yC3=nNO2/ntot; yD3=0;


%% Stoichiometric Coefficients
vA1=2; vB1=0; vC1=2; vD1=1;
vA2=1; vB2=1; vC2=2; vD2=0;
vA3=1; vB3=2; vC3=2; vD3=0;


%% Equilibrium Constants
F=[(yC1^vC1*yD1^vD1)/(yA1^vA1*yB1^vB1)*(P/P0)^(vC1+vD1-vA1-vB1)-K1;
    (yC2^vC2*yD2^vD2)/(yA2^vA2*yB2^vB2)*(P/P0)^(vC2+vD2-vA2-vB2)-K2;
    (yC3^vC3*yD3^vD3)/(yA3^vA3*yB3^vB3)*(P/P0)^(vC3+vD3-vA3-vB3)-K3];
end

x=fsolve(@(x) HW3_d(x),[0.3;0.3;0.3])


%% Combustion Reaction: CH8H18 + 37.5*(02 + 3.76*N2) -> 8*CO2 + 9*H2O + 141*N2 + 25*02
    % Number of moles
nH2O = 9
nCO2 = 8 - 2*x(1)
nN2 = 141 - x(2) - x(3)
nO2 = 25 + x(1) - x(2) - 2*x(3)


%% Dissociation Reactions 

% Reaction #1: 2CO2 -> 2CO + O2
% Reaction #2: N2 + 02 -> 2NO
% Reaction #3: N2 + 202 -> 2NO2

% Number of Moles
nCO = 2*x(1)
nNO = 2*x(2)
nNO2 = 2*x(3)
ntot = nCO2 + nH2O + nN2 + nO2 + nCO + nNO + nNO2;


%% Mole Fractions
yH2O=nH2O/ntot
yO2=nO2/ntot
yN2=nN2/ntot
yCO2=nCO2/ntot
yCO=nCO/ntot
yNO=nNO/ntot
yNO2=nNO2/ntot
ytot = yH2O + yO2 + yN2 + yCO2 + yCO + yNO + yNO2
end