Prev: format_ticks.m
Next: plotting
From: Gavin on 15 Jun 2010 07:08
I'm trying to model a biomass gasifier using chemical equilibrium. I have a system of equations and unknowns that I would like to solve with Matlab using the "fsolve" command. I have two m-files which I use. The first m-file is named "gasf_funcB.m" where I create my system of equations to be used by fsolve. The second m-file is named "gasf_solB.m" where I state my inputs (total moles of C, H, O, N into the system) and initial guess for the outputs (mole fractions) obtained by fsolve. I am trying to solve the system for the unknowns which are the compounds created from the reaction. I am getting complex answers from using fsolve. Please help and let me know what I'm doing wrong. All m-files are shown below. My output to the workspace from running gasf_solB.m is also shown below.

---------------------------------------- Here is my gasf_funcB.m file --------------------------------------
% gasifier functions

function F = gasf_funcB(x,nC,nH,nO,nN,T)

% Overall reaction is
% CHON + air -> CO + CO2 + H2O + H2 + H + N + N2 + NO + O2 + O + OH
yCO = x(1);
yCO2 = x(2);
yH2O = x(3);
yH2 = x(4);
yH = x(5);
yN = x(6);
yN2 = x(7);
yNO = x(8);
yO2 = x(9);
yO = x(10);
yOH = x(11);

% total number of moles
nt = yCO+yCO2+yH2O+yH2+yH+yN+yN2+yNO+yO2+yO+yOH;

% determine equilibrium constant at specified temperature
[K] = K_est(T);
K1 = K(1);
K2 = K(2);
K3 = K(3);
K4 = K(4);
K5 = K(5);
K6 = K(6);
K7 = K(7);
K8 = K(8);

% balance equations to solve
F(1) = yCO+yCO2-nC/nt; % C balance
F(2) = 2*yH2O+2*yH2+yH+yOH-nH/nt; % H balance
F(3) = yCO+2*yCO2+yH2O+2*yO2+yO+yNO+yOH-nO/nt; % O balance
F(4) = 2*yN2+yN+yNO-nN/nt; % N balance

% equilibrium reactions to solve where 1) refers to K1 etc...
% assume P = 1 atm therefore P/Pref = 1
% 1) CO2 <-> CO + 0.5O2
% 2) CO2 + H2 <-> CO + H2O
% 3) N2 <-> 2N
% 4) H2 <-> 2H
% 5) O2 <-> 2O
% 6) 0.5O2 + 0.5N2 <-> NO
% 7) H2O <-> H2 + 0.5O2
% 8) H2O <-> OH + 0.5H2
F(5) = (yCO*yO2^0.5)/yCO2*(1/nt)^0.5-K1;
F(6) = (yCO*yH2O)/(yCO2*yH2)-K2;
F(7) = yN^2/yN2*(1/nt)-K3;
F(8) = yH^2/yH2*(1/nt)-K4;
F(9) = yO^2/yO2*(1/nt)-K5;
F(10) = yNO/(yO2^0.5*yN2^0.5)-K6;
F(11) = (yO2^0.5*yH2)/yH2O*(1/nt)^0.5-K7;
F(12) = (yOH*yH2^0.5)/yH2O*(1/nt)^0.5-K8;

end

% estimate of equilibrium constant (K) at temperature T in Kelvin for
% various reactions

function [K] = K_est(T)

% database of temperatures in Kelvin
Tdb = [298 500 1000 1200 1400 1600 1700 1800 1900 2000 2100 2200 2300 ...
2400 2500 2600 2700 2800 2900 3000];

% database values for log base 10 equilibrium
% these correspond to equilibrium reactions stated previously
log10K1 = [-45.066 -25.025 -10.221 -7.764 -6.014 -4.706 -4.169 -3.693...
-3.267 -2.884 -2.539 -2.226 -1.940 -1.679 -1.440 -1.219 -1.015...
-0.825 -0.649 -0.485];
log10K2 = [-5.018 -2.139 -0.159 0.135 0.333 0.474 0.530 0.577 0.619...
0.656 0.688 0.716 0.742 0.764 0.784 0.802 0.818 0.833 0.846...
0.858];
log10K3 = [-159.6 -92.672 -43.056 -34.754 -28.812 -24.350 -22.512...
-20.874 -19.410 -18.092 -16.898 -15.810 -14.818 -13.908...
-13.070 -12.298 -11.580 -10.914 -10.294 -9.716];
log10K4 = [-71.224 -40.316 -17.292 -13.414 -10.63 -8.532 -7.666 -6.896...
-6.204 -5.58 -5.016 -4.502 -4.032 -3.6 -3.202 -2.836 -2.494...
-2.178 -1.882 -1.606];
log10K5 = [-81.208 -45.88 -19.614 -15.208 -12.054 -9.684 -8.706 -7.836...
-7.058 -6.356 -5.72 -5.142 -4.614 -4.13 -3.684 -3.272 -2.892...
-2.536 -2.206 -1.898];
log10K6 = [-15.171 -8.783 -4.062 -3.275 -2.712 -2.29 -2.116 -1.962...
-1.823 -1.699 -1.586 -1.484 -1.391 -1.305 -1.227 -1.154...
-1.087 -1.025 -0.967 -0.913];
log10K7 = [-40.048 -22.886 -10.062 -7.899 -6.347 -5.18 -4.699 -4.27...
-3.886 -3.54 -3.227 -2.942 -2.682 -2.443 -2.224 -2.021...
-1.833 -1.658 -1.495 -1.343];
log10K8 = [-46.054 -26.130 -11.280 -8.811 -7.021 -5.677 -5.124 -4.613...
-4.19 -3.776 -3.434 -3.091 -2.809 -2.52 -2.27 -2.038 -1.823...
-1.624 -1.438 -1.265];

% interpolate equilibrium constant at T
K1_interp = interp1(Tdb,log10K1,T,'spline','extrap');
K1 = 10^K1_interp;
K2_interp = interp1(Tdb,log10K2,T,'spline','extrap');
K2 = 10^K2_interp;
K3_interp = interp1(Tdb,log10K3,T,'spline','extrap');
K3 = 10^K3_interp;
K4_interp = interp1(Tdb,log10K4,T,'spline','extrap');
K4 = 10^K4_interp;
K5_interp = interp1(Tdb,log10K5,T,'spline','extrap');
K5 = 10^K5_interp;
K6_interp = interp1(Tdb,log10K6,T,'spline','extrap');
K6 = 10^K6_interp;
K7_interp = interp1(Tdb,log10K7,T,'spline','extrap');
K7 = 10^K7_interp;
K8_interp = interp1(Tdb,log10K8,T,'spline','extrap');
K8 = 10^K8_interp;

K = [K1 K2 K3 K4 K5 K6 K7 K8];

end

----------------------------------------------------------------------------------------------------------------------

--------------------------------------------- Here is my gasf_solB.m file -----------------------------------
% gasifier equilibrium solution
clear
clc

% Inputs to system, total moles of C, H, O, N coming from CHON and air
nC = 2; % C total moles
nH = 2; % H total moles
nO = 4; % O total moles
nN = 1; % N total moles
T = 1000; % equilibrium temperature, K

% initial guess of outputs
x0 = [.1 .1 .1 .1 .1 .1 .1 .1 .1 .1 .1];

options = optimset('Display','iter','MaxIter',10000,'MaxFunEvals',10000);
[x,exitflag] = fsolve(@(x) gasf_funcB(x,nC,nH,nO,nN,T),x0,options);

yCO = x(1)
yCO2 = x(2)
yH2O = x(3)
yH2 = x(4)
yH = x(5)
yN = x(6)
yN2 = x(7)
yNO = x(8)
yO2 = x(9)
yO = x(10)
yOH = x(11)
nt = yCO+yCO2+yH2O+yH2+yH+yN+yN2+yNO+yO2+yO+yOH
T

----------------------------------------------------------------------------------------------------------------------

-------------------------------- Output from running gasf_solB.m file ---------------------------------
Warning: Trust-region-dogleg algorithm of FSOLVE cannot handle non-square systems; using
Levenberg-Marquardt algorithm instead.
> In fsolve at 324
In gasf_solB at 16

First-Order Norm of
Iteration Func-count Residual optimality Lambda step
0 12 13.2407 25.4 0.01
1 24 0.924578 4.48 0.001 0.478901
2 36 0.0103891 0.336 0.0001 0.288275
3 48 0.00146625 0.142 1e-05 0.103138
4 65 0.00115348 0.146 1 0.0106871
5 77 0.000951684 0.143 0.1 0.00882184
6 90 0.000838769 0.148 1 0.00729546
7 102 0.000739851 0.145 0.1 0.00620988
8 115 0.00067335 0.148 1 0.00564338
9 127 0.000609141 0.146 0.1 0.00504862
10 140 0.000564525 0.148 1 0.00464933
143 1806 3.28161e-06 1.65 10000 4.43205e-06
144 1819 2.42909e-06 1.18 100000 1.00569e-06
145 1831 2.26005e-06 0.0609 10000 1.69937e-07
Optimization terminated: the magnitude of the search direction is less than options.TolX
yCO =
0.32134 - 1.4945e-06i
yCO2 =
0.74179 + 1.223e-06i
yH2O =
0.32147 + 1.2992e-06i
yH2 =
0.20083 - 5.1165e-07i
yH =
0.018548 + 2.2209e-07i
yN =
0.022879 + 1.2783e-07i
yN2 =
0.25434 - 8.5025e-07i
yNO =
-5.2241e-09 + 3.6983e-08i
yO2 =
-5.1624e-07 + 1.7704e-08i
yO =
-1.6125e-05 - 4.0885e-07i
yOH =
-2.9233e-06 - 1.6344e-07i
nt =
1.8812 - 5.0192e-07i
T =
1000
----------------------------------------------------------------------------------------------------------------------
From: Torsten Hennig on 15 Jun 2010 03:33
> I'm trying to model a biomass gasifier using chemical
> equilibrium. I have a system of equations and
> unknowns that I would like to solve with Matlab using
> the "fsolve" command. I have two m-files which I
> use. The first m-file is named "gasf_funcB.m" where
> I create my system of equations to be used by fsolve.
> The second m-file is named "gasf_solB.m" where I
> I state my inputs (total moles of C, H, O, N into the
> system) and initial guess for the outputs (mole
> fractions) obtained by fsolve. I am trying to solve
> the system for the unknowns which are the compounds
> created from the reaction. I am getting complex
> answers from using fsolve. Please help and let me
> know what I'm doing wrong. All m-files are shown
> below. My output to the workspace from running
> gasf_solB.m is also shown below.
>
> ---------------------------------------- Here is my
> gasf_funcB.m file
> --------------------------------------
> % gasifier functions
>
> function F = gasf_funcB(x,nC,nH,nO,nN,T)
>
> % Overall reaction is
> % CHON + air -> CO + CO2 + H2O + H2 + H + N + N2 + NO
> + O2 + O + OH
> yCO = x(1);
> yCO2 = x(2);
> yH2O = x(3);
> yH2 = x(4);
> yH = x(5);
> yN = x(6);
> yN2 = x(7);
> yNO = x(8);
> yO2 = x(9);
> yO = x(10);
> yOH = x(11);
>
> % total number of moles
> nt = yCO+yCO2+yH2O+yH2+yH+yN+yN2+yNO+yO2+yO+yOH;
>
> % determine equilibrium constant at specified
> temperature
> [K] = K_est(T);
> K1 = K(1);
> K2 = K(2);
> K3 = K(3);
> K4 = K(4);
> K5 = K(5);
> K6 = K(6);
> K7 = K(7);
> K8 = K(8);
>
> % balance equations to solve
> F(1) = yCO+yCO2-nC/nt;
> % C
> % C
> % C balance
> F(2) = 2*yH2O+2*yH2+yH+yOH-nH/nt;
> % H balance
> F(3) = yCO+2*yCO2+yH2O+2*yO2+yO+yNO+yOH-nO/nt; % O
> balance
> F(4) = 2*yN2+yN+yNO-nN/nt;
> % N balance
>
> % equilibrium reactions to solve where 1) refers to
> K1 etc...
> % assume P = 1 atm therefore P/Pref = 1
> % 1) CO2 <-> CO + 0.5O2
> % 2) CO2 + H2 <-> CO + H2O
> % 3) N2 <-> 2N
> % 4) H2 <-> 2H
> % 5) O2 <-> 2O
> % 6) 0.5O2 + 0.5N2 <-> NO
> % 7) H2O <-> H2 + 0.5O2
> % 8) H2O <-> OH + 0.5H2
> F(5) = (yCO*yO2^0.5)/yCO2*(1/nt)^0.5-K1;
> F(6) = (yCO*yH2O)/(yCO2*yH2)-K2;
> F(7) = yN^2/yN2*(1/nt)-K3;
> F(8) = yH^2/yH2*(1/nt)-K4;
> F(9) = yO^2/yO2*(1/nt)-K5;
> F(10) = yNO/(yO2^0.5*yN2^0.5)-K6;
> F(11) = (yO2^0.5*yH2)/yH2O*(1/nt)^0.5-K7;
> F(12) = (yOH*yH2^0.5)/yH2O*(1/nt)^0.5-K8;
>
> end
>
> % estimate of equilibrium constant (K) at temperature
> T in Kelvin for
> % various reactions
>
> function [K] = K_est(T)
>
> % database of temperatures in Kelvin
> Tdb = [298 500 1000 1200 1400 1600 1700 1800 1900
> 2000 2100 2200 2300 ...
> 2400 2500 2600 2700 2800 2900 3000];
>
> % database values for log base 10 equilibrium
> % these correspond to equilibrium reactions stated
> previously
> log10K1 = [-45.066 -25.025 -10.221 -7.764 -6.014
> -4.706 -4.169 -3.693...
> -3.267 -2.884 -2.539 -2.226 -1.940 -1.679
> 940 -1.679 -1.440 -1.219 -1.015...
> -0.825 -0.649 -0.485];
> log10K2 = [-5.018 -2.139 -0.159 0.135 0.333 0.474
> 0.530 0.577 0.619...
> 0.656 0.688 0.716 0.742 0.764 0.784 0.802
> .784 0.802 0.818 0.833 0.846...
> 0.858];
> log10K3 = [-159.6 -92.672 -43.056 -34.754 -28.812
> -24.350 -22.512...
> -20.874 -19.410 -18.092 -16.898 -15.810
> 98 -15.810 -14.818 -13.908...
> -13.070 -12.298 -11.580 -10.914 -10.294
> 14 -10.294 -9.716];
> log10K4 = [-71.224 -40.316 -17.292 -13.414 -10.63
> -8.532 -7.666 -6.896...
> -6.204 -5.58 -5.016 -4.502 -4.032 -3.6
> 4.032 -3.6 -3.202 -2.836 -2.494...
> -2.178 -1.882 -1.606];
> log10K5 = [-81.208 -45.88 -19.614 -15.208 -12.054
> -9.684 -8.706 -7.836...
> -7.058 -6.356 -5.72 -5.142 -4.614 -4.13
> .614 -4.13 -3.684 -3.272 -2.892...
> -2.536 -2.206 -1.898];
> log10K6 = [-15.171 -8.783 -4.062 -3.275 -2.712 -2.29
> -2.116 -1.962...
> -1.823 -1.699 -1.586 -1.484 -1.391 -1.305
> .391 -1.305 -1.227 -1.154...
> -1.087 -1.025 -0.967 -0.913];
> log10K7 = [-40.048 -22.886 -10.062 -7.899 -6.347
> -5.18 -4.699 -4.27...
> -3.886 -3.54 -3.227 -2.942 -2.682 -2.443
> .682 -2.443 -2.224 -2.021...
> -1.833 -1.658 -1.495 -1.343];
> log10K8 = [-46.054 -26.130 -11.280 -8.811 -7.021
> -5.677 -5.124 -4.613...
> -4.19 -3.776 -3.434 -3.091 -2.809 -2.52
> 2.809 -2.52 -2.27 -2.038 -1.823...
> -1.624 -1.438 -1.265];
>
> % interpolate equilibrium constant at T
> K1_interp = interp1(Tdb,log10K1,T,'spline','extrap');
> K1 = 10^K1_interp;
> K2_interp = interp1(Tdb,log10K2,T,'spline','extrap');
> K2 = 10^K2_interp;
> K3_interp = interp1(Tdb,log10K3,T,'spline','extrap');
> K3 = 10^K3_interp;
> K4_interp = interp1(Tdb,log10K4,T,'spline','extrap');
> K4 = 10^K4_interp;
> K5_interp = interp1(Tdb,log10K5,T,'spline','extrap');
> K5 = 10^K5_interp;
> K6_interp = interp1(Tdb,log10K6,T,'spline','extrap');
> K6 = 10^K6_interp;
> K7_interp = interp1(Tdb,log10K7,T,'spline','extrap');
> K7 = 10^K7_interp;
> K8_interp = interp1(Tdb,log10K8,T,'spline','extrap');
> K8 = 10^K8_interp;
>
> K = [K1 K2 K3 K4 K5 K6 K7 K8];
>
> end
>
> ------------------------------------------------------
> ------------------------------------------------------
> ----------
>
> --------------------------------------------- Here is
> my gasf_solB.m file
> -----------------------------------
> % gasifier equilibrium solution
> clear
> clc
>
> % Inputs to system, total moles of C, H, O, N coming
> from CHON and air
> nC = 2; % C total moles
> nH = 2; % H total moles
> nO = 4; % O total moles
> nN = 1; % N total moles
> T = 1000; % equilibrium temperature, K
>
> % initial guess of outputs
> x0 = [.1 .1 .1 .1 .1 .1 .1 .1 .1 .1 .1];
>
> options =
> optimset('Display','iter','MaxIter',10000,'MaxFunEvals
> ',10000);
> [x,exitflag] = fsolve(@(x)
> gasf_funcB(x,nC,nH,nO,nN,T),x0,options);
>
> yCO = x(1)
> yCO2 = x(2)
> yH2O = x(3)
> yH2 = x(4)
> yH = x(5)
> yN = x(6)
> yN2 = x(7)
> yNO = x(8)
> yO2 = x(9)
> yO = x(10)
> yOH = x(11)
> nt = yCO+yCO2+yH2O+yH2+yH+yN+yN2+yNO+yO2+yO+yOH
> T
>
> ------------------------------------------------------
> ------------------------------------------------------
> ----------
>
> -------------------------------- Output from running
> gasf_solB.m file ---------------------------------
> Warning: Trust-region-dogleg algorithm of FSOLVE
> cannot handle non-square systems; using
> Levenberg-Marquardt algorithm instead.
> > In fsolve at 324
> In gasf_solB at 16
>
> First-Order
> First-Order
> First-Order Norm
> Norm of
> Iteration Func-count Residual optimality
> y Lambda step
> 0 12 13.2407 25.4
> 25.4 0.01
> 1 24 0.924578 4.48
> 4.48 0.001 0.478901
> 2 36 0.0103891 0.336
> 0.336 0.0001 0.288275
> 3 48 0.00146625 0.142
> 0.142 1e-05 0.103138
> 4 65 0.00115348 0.146
> 0.146 1 0.0106871
> 5 77 0.000951684 0.143
> 0.143 0.1 0.00882184
> 6 90 0.000838769 0.148
> 0.148 1 0.00729546
> 7 102 0.000739851 0.145
> 0.145 0.1 0.00620988
> 8 115 0.00067335 0.148
> 0.148 1 0.00564338
> 9 127 0.000609141 0.146
> 0.146 0.1 0.00504862
> 10 140 0.000564525 0.148
> .148 1 0.00464933
> 143 1806 3.28161e-06 1.65
> .65 10000 4.43205e-06
> 144 1819 2.42909e-06 1.18
> .18 100000 1.00569e-06
> 145 1831 2.26005e-06 0.0609
> 609 10000 1.69937e-07
> Optimization terminated: the magnitude of the search
> direction is less than options.TolX
> yCO =
> 0.32134 - 1.4945e-06i
> yCO2 =
> 0.74179 + 1.223e-06i
> yH2O =
> 0.32147 + 1.2992e-06i
> yH2 =
> 0.20083 - 5.1165e-07i
> yH =
> 0.018548 + 2.2209e-07i
> yN =
> 0.022879 + 1.2783e-07i
> yN2 =
> 0.25434 - 8.5025e-07i
> yNO =
> -5.2241e-09 + 3.6983e-08i
> yO2 =
> -5.1624e-07 + 1.7704e-08i
> yO =
> -1.6125e-05 - 4.0885e-07i
> yOH =
> -2.9233e-06 - 1.6344e-07i
> nt =
> 1.8812 - 5.0192e-07i
> T =
> 1000
> ------------------------------------------------------
> ------------------------------------------------------
> ----------

Hi Gavin,

during iteration with fsolve some of your mole fractions
seem to become negative - and taking the square root
gives complex values.
Whereever there is a square root of a mole fraction y,
you should substitute y := y~^2 which makes the
square root disappear.

Best wishes
Torsten.
From: Gavin on 15 Jun 2010 09:45
Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote in message <963770834.348103.1276601652274.JavaMail.root(a)gallium.mathforum.org>...
> >
> > ------------------------------------------------------
> > ----------
>
> Hi Gavin,
>
> during iteration with fsolve some of your mole fractions
> seem to become negative - and taking the square root
> gives complex values.
> Whereever there is a square root of a mole fraction y,
> you should substitute y := y~^2 which makes the
> square root disappear.
>
> Best wishes
> Torsten.

Torsten,

I have changed the following equations to the following to prevent complex numbers:

F(5) = (yCO*(yO2^2)^0.5)/yCO2*(1/nt)^0.5-K1;
F(6) = (yCO*yH2O)/(yCO2*yH2)-K2;
F(7) = yN^2/yN2*(1/nt)-K3;
F(8) = yH^2/yH2*(1/nt)-K4;
F(9) = yO^2/yO2*(1/nt)-K5;
F(10) = yNO/((yO2^2)^0.5*yN2)^0.5-K6;
F(11) = ((yO2^2)^0.5*yH2)/yH2O*(1/nt)^0.5-K7;
F(12) = (yOH*(yH2^2)^0.5)/yH2O*(1/nt)^0.5-K8;

The complex answers are gone but I'm getting negative values for some of my answers which is not correct:
yO2 =
-2.2808e-10
yO =
-9.1383e-11

Is there a way to tell fsolve that a term cannot be negative when evaluating a system of equations?
From: Torsten Hennig on 15 Jun 2010 06:17
> Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote
> in message
> <963770834.348103.1276601652274.JavaMail.root(a)gallium.
> mathforum.org>...
> > >
> > >
> ------------------------------------------------------
> > > ----------
> >
> > Hi Gavin,
> >
> > during iteration with fsolve some of your mole
> fractions
> > seem to become negative - and taking the square
> root
> > gives complex values.
> > Whereever there is a square root of a mole fraction
> y,
> > you should substitute y := y~^2 which makes the
> > square root disappear.
> >
> > Best wishes
> > Torsten.
>
> Torsten,
>
> I have changed the following equations to the
> following to prevent complex numbers:
>
> F(5) = (yCO*(yO2^2)^0.5)/yCO2*(1/nt)^0.5-K1;
> F(6) = (yCO*yH2O)/(yCO2*yH2)-K2;
> F(7) = yN^2/yN2*(1/nt)-K3;
> F(8) = yH^2/yH2*(1/nt)-K4;
> F(9) = yO^2/yO2*(1/nt)-K5;
> F(10) = yNO/((yO2^2)^0.5*yN2)^0.5-K6;
> F(11) = ((yO2^2)^0.5*yH2)/yH2O*(1/nt)^0.5-K7;
> F(12) = (yOH*(yH2^2)^0.5)/yH2O*(1/nt)^0.5-K8;
>
> The complex answers are gone but I'm getting negative
> values for some of my answers which is not correct:
> yO2 =
> -2.2808e-10
> yO =
> -9.1383e-11
>
> Is there a way to tell fsolve that a term cannot be
> negative when evaluating a system of equations?

Hi Gavin,

when returning from fsolve, did you retransform your
solution, i.e. squared the mole fraction for
which you inserted the square in the algebraic
equations ?

Best wishes
Torsten.
From: Gavin on 15 Jun 2010 11:45
Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote in message <800314764.348942.1276611471963.JavaMail.root(a)gallium.mathforum.org>...
> >
> > Is there a way to tell fsolve that a term cannot be
> > negative when evaluating a system of equations?
>
> Hi Gavin,
>
> when returning from fsolve, did you retransform your
> solution, i.e. squared the mole fraction for
> which you inserted the square in the algebraic
> equations ?
>
> Best wishes
> Torsten.

I'm not sure what you mean by "retransform your solution when returning to fsolve". Also, another way of dealing with the complex numbers was to make sure the mole fraction could not be less than 0. I added the following to the gasf_funcB.m file:

if yCO<=0; yCO=1e-20; end
if yCO2<=0; yCO2=1e-20; end
if yH2O<=0; yH2O=1e-20; end
if yH2<=0; yH2=1e-20; end
if yH<=0; yH=1e-20; end
if yN<=0; yN=1e-20; end
if yN2<=0; yN2=1e-20; end
if yNO<=0; yNO=1e-20; end
if yO2<=0; yO2=1e-20; end
if yO<=0; yO=1e-20; end
if yOH<=0; yOH=1e-20; end
 |  Next  |  Last
Pages: 1 2 3
Prev: format_ticks.m
Next: plotting