
function [out,n5,error,n1,n2,n3,n4,n6,ItterationsDiss]=mydissfuniii(carbon,hydrogen,oxygen,temperature,fuelequiv,molar_air_stoich,pressure)
n5ntvalue=0.1;
error=1;

%TEMP CO2     H2O          CO
gibbseng=[0	-393.151	-238.921	-113.805
    100	-393.683	-236.584	-120.239
    200	-394.085	-232.766	-128.526
    298.15	-394.389	-228.582	-137.163
    300	-394.394	-228.5	-137.328
    400	-394.675	-223.901	-146.338
    500	-394.939	-219.051	-155.414
    600	-395.182	-214.007	-164.486
    700	-395.398	-208.812	-173.518
    800	-395.586	-203.496	-182.497
    900	-395.748	-198.083	-191.416
    1000	-395.886	-192.59	-200.275
    1100	-396.001	-187.033	-209.075
    1200	-396.098	-181.425	-217.819
    1300	-396.177	-175.774	-226.509
    1400	-396.24	-170.089	-235.149
    1500	-396.288	-164.376	-243.74
    1600	-396.323	-158.639	-252.284
    1700	-396.344	-152.883	-260.784
    1800	-396.353	-147.111	-269.242
    1900	-396.349	-141.325	-277.658
    2000	-396.333	-135.528	-286.034
    2100	-396.304	-129.721	-294.372
    2200	-396.262	-123.905	-302.672
    2300	-396.209	-118.082	-310.936
    2400	-396.142	-112.252	-319.165
    2500	-396.062	-106.416	-327.358
    2600	-395.969	-100.575	-335.517
    2700	-395.862	-94.729	-343.643
    2800	-395.742	-88.878	-351.736
    2900	-395.609	-83.023	-359.797
    3000	-395.461	-77.163	-367.826
    3100	-395.298	-71.298	-375.824
    3200	-395.122	-65.43	-383.792
    3300	-394.932	-59.558	-391.73
    3400	-394.726	-53.681	-399.638
    3500	-394.506	-47.801	-407.517
    3600	-394.271	-41.916	-415.366
    3700	-394.022	-36.027	-423.188
    3800	-393.756	-30.133	-430.981
    3900	-393.477	-24.236	-438.746
    4000	-393.183	-18.334	-446.485
    4100	-392.874	-12.427	-454.196
    4200	-392.55	-6.516	-461.881
    4300	-392.21	-0.6	-469.539
    4400	-391.857	5.32	-477.171
    4500	-391.488	11.245	-484.778
    4600	-391.105	17.175	-492.359
    4700	-390.706	23.111	-499.915
    4800	-390.292	29.052	-507.446
    4900	-389.862	34.998	-514.951
    5000	-389.419	40.949	-522.433
    5100	-388.959	46.906	-529.89
    5200	-388.486	52.869	-537.323
    5300	-387.996	58.838	-544.732
    5400	-387.493	64.811	-552.118
    5500	-386.974	70.791	-559.48
    5600	-386.439	76.777	-566.818
    5700	-385.89	82.769	-574.133
    5800	-385.324	88.767	-581.424
    5900	-384.745	94.77	-588.694
    6000	-384.148	100.78	-595.939];

%Gibbs free energy
gco2=interp1(gibbseng(1:end,1),gibbseng(1:end,2),temperature)*10^3;
gh2o=interp1(gibbseng(1:end,1),gibbseng(1:end,3),temperature)*10^3;
gco=interp1(gibbseng(1:end,1),gibbseng(1:end,4),temperature)*10^3;

%Equilibrium constants
kco2=exp((-(gco2-gco))*10^3/(8314.46*temperature));
kh2o=exp((-gh2o)*10^3/(8314.46*temperature));

itt = 0;

%Goodger method for calculation of dissociation
%Iterate through n5/nt guess until error goes negative or guess reaches 0
while error>0
        itt=itt+1;
    a=sqrt(n5ntvalue*pressure);
    n3=(carbon*fuelequiv)/(kco2*a+1);
    n4=((hydrogen/2)*fuelequiv)/(kh2o*a+1);
    n1=carbon*fuelequiv-n3;
    n2=(hydrogen/2)*fuelequiv-n4;
    n5=(molar_air_stoich*2+oxygen*fuelequiv-2*n1-n2-n3)/2;
    n6=molar_air_stoich*3.76;
    ntsum=n1+n2+n3+n4+n5+n6;
    error=n5ntvalue-n5/ntsum;

    n5ntvalue=n5ntvalue-0.001;
    if n5ntvalue<0
        break
    end
end
%---------------------------------------------
%If guess reached 0, reset guess to 0.001 and iterate down in 10^-3 steps
if n5ntvalue<0
    n5ntvalue=0.001;
    while error>0
            itt=itt+1;
        a=sqrt(n5ntvalue*pressure);
        n3=(carbon*fuelequiv)/(kco2*a+1);
        n4=((hydrogen/2)*fuelequiv)/(kh2o*a+1);
        n1=carbon*fuelequiv-n3;
        n2=(hydrogen/2)*fuelequiv-n4;
        n5=(molar_air_stoich*2+oxygen*fuelequiv-2*n1-n2-n3)/2;
        n6=molar_air_stoich*3.76;
        ntsum=n1+n2+n3+n4+n5+n6;
        error=n5ntvalue-n5/ntsum;

        n5ntvalue=n5ntvalue-0.000001;
        if n5ntvalue<0
            break
        end
    end
    %if error went negative, iterate guess in +10^-6 increments until error is
    %positive
else
    while error<0
            itt=itt+1;
        a=sqrt(n5ntvalue*pressure);
        n3=(carbon*fuelequiv)/(kco2*a+1);
        n4=((hydrogen/2)*fuelequiv)/(kh2o*a+1);
        n1=carbon*fuelequiv-n3;
        n2=(hydrogen/2)*fuelequiv-n4;
        n5=(molar_air_stoich*2+oxygen*fuelequiv-2*n1-n2-n3)/2;
        n6=molar_air_stoich*3.76;
        ntsum=n1+n2+n3+n4+n5+n6;
        error=n5ntvalue-n5/ntsum;

        n5ntvalue=n5ntvalue+0.000001; %-6
    end
end
%----------------------------------------------------
%repeat with smaller increments
%now positive error
if n5ntvalue<0
    n5ntvalue=0.000001;
    while error>0
            itt=itt+1;
        a=sqrt(n5ntvalue*pressure);
        n3=(carbon*fuelequiv)/(kco2*a+1);
        n4=((hydrogen/2)*fuelequiv)/(kh2o*a+1);
        n1=carbon*fuelequiv-n3;
        n2=(hydrogen/2)*fuelequiv-n4;
        n5=(molar_air_stoich*2+oxygen*fuelequiv-2*n1-n2-n3)/2;
        n6=molar_air_stoich*3.76;
        ntsum=n1+n2+n3+n4+n5+n6;
        error=n5ntvalue-n5/ntsum;

        n5ntvalue=n5ntvalue-0.000000001;
        if n5ntvalue<0
            break
        end
    end

else
    while error>0
            itt=itt+1;
        a=sqrt(n5ntvalue*pressure);
        n3=(carbon*fuelequiv)/(kco2*a+1);
        n4=((hydrogen/2)*fuelequiv)/(kh2o*a+1);
        n1=carbon*fuelequiv-n3;
        n2=(hydrogen/2)*fuelequiv-n4;
        n5=(molar_air_stoich*2+oxygen*fuelequiv-2*n1-n2-n3)/2;
        n6=molar_air_stoich*3.76;
        ntsum=n1+n2+n3+n4+n5+n6;
        error=n5ntvalue-n5/ntsum;

        n5ntvalue=n5ntvalue-0.000000001; %-9
    end
end
%--------------------------------------------------
%error now negative
if n5ntvalue<0
    n5ntvalue=0.000000001;
    while error>0
            itt=itt+1;
        a=sqrt(n5ntvalue*pressure);
        n3=(carbon*fuelequiv)/(kco2*a+1);
        n4=((hydrogen/2)*fuelequiv)/(kh2o*a+1);
        n1=carbon*fuelequiv-n3;
        n2=(hydrogen/2)*fuelequiv-n4;
        n5=(molar_air_stoich*2+oxygen*fuelequiv-2*n1-n2-n3)/2;
        n6=molar_air_stoich*3.76;
        ntsum=n1+n2+n3+n4+n5+n6;
        error=n5ntvalue-n5/ntsum;

        n5ntvalue=n5ntvalue-0.000000000001;
        if n5ntvalue<0
            break
        end
    end

else
    while error<0
            itt=itt+1;
        a=sqrt(n5ntvalue*pressure);
        n3=(carbon*fuelequiv)/(kco2*a+1);
        n4=((hydrogen/2)*fuelequiv)/(kh2o*a+1);
        n1=carbon*fuelequiv-n3;
        n2=(hydrogen/2)*fuelequiv-n4;
        n5=(molar_air_stoich*2+oxygen*fuelequiv-2*n1-n2-n3)/2;
        n6=molar_air_stoich*3.76;
        ntsum=n1+n2+n3+n4+n5+n6;
        error=n5ntvalue-n5/ntsum;

        n5ntvalue=n5ntvalue+0.000000000001; %-12
    end
end
% %-------------------------------------

%uncomment the rest if higher accuracy is desired. (increases run time)

% %error now positive
% if n5ntvalue<0
% n5ntvalue=0.000000000001;
% while error>0
% a=sqrt(n5ntvalue*pressure);
% n3=(carbon*fuelequiv)/(kco2*a+1);
% n4=((hydrogen/2)*fuelequiv)/(kh2o*a+1);
% n1=carbon*fuelequiv-n3;
% n2=(hydrogen/2)*fuelequiv-n4;
% n5=(molar_air_stoich*2+oxygen*fuelequiv-2*n1-n2-n3)/2;
% n6=molar_air_stoich*3.76;
% ntsum=n1+n2+n3+n4+n5+n6;
% error=n5ntvalue-n5/ntsum;
%
% n5ntvalue=n5ntvalue-0.000000000000001;
% if n5ntvalue<0
% break
% end
% end
%
% else
% while error>0
% a=sqrt(n5ntvalue*pressure);
% n3=(carbon*fuelequiv)/(kco2*a+1);
% n4=((hydrogen/2)*fuelequiv)/(kh2o*a+1);
% n1=carbon*fuelequiv-n3;
% n2=(hydrogen/2)*fuelequiv-n4;
% n5=(molar_air_stoich*2+oxygen*fuelequiv-2*n1-n2-n3)/2;
% n6=molar_air_stoich*3.76;
% ntsum=n1+n2+n3+n4+n5+n6;
% error=n5ntvalue-n5/ntsum;
%
% n5ntvalue=n5ntvalue-0.000000000000001; %-15
% end
% end
% %---------------------------
% %error now negative
% if n5ntvalue<0
% n5ntvalue=0.000000000000001;
% while error>0
% a=sqrt(n5ntvalue*pressure);
% n3=(carbon*fuelequiv)/(kco2*a+1);
% n4=((hydrogen/2)*fuelequiv)/(kh2o*a+1);
% n1=carbon*fuelequiv-n3;
% n2=(hydrogen/2)*fuelequiv-n4;
% n5=(molar_air_stoich*2+oxygen*fuelequiv-2*n1-n2-n3)/2;
% n6=molar_air_stoich*3.76;
% ntsum=n1+n2+n3+n4+n5+n6;
% error=n5ntvalue-n5/ntsum;
%
% n5ntvalue=n5ntvalue-0.000000000000000001;
% if n5ntvalue<0
% break
% end
% end
%
% else
% while error<0
% a=sqrt(n5ntvalue*pressure);
% n3=(carbon*fuelequiv)/(kco2*a+1);
% n4=((hydrogen/2)*fuelequiv)/(kh2o*a+1);
% n1=carbon*fuelequiv-n3;
% n2=(hydrogen/2)*fuelequiv-n4;
% n5=(molar_air_stoich*2+oxygen*fuelequiv-2*n1-n2-n3)/2;
% n6=molar_air_stoich*3.76;
% ntsum=n1+n2+n3+n4+n5+n6;
% error=n5ntvalue-n5/ntsum;
%
% n5ntvalue=n5ntvalue+0.000000000000000001; %-18
% end
% end
% %-------------------------
% %error now positive
% if n5ntvalue<0
% n5ntvalue=0.000000000000000001;
% while error>0
% a=sqrt(n5ntvalue*pressure);
% n3=(carbon*fuelequiv)/(kco2*a+1);
% n4=((hydrogen/2)*fuelequiv)/(kh2o*a+1);
% n1=carbon*fuelequiv-n3;
% n2=(hydrogen/2)*fuelequiv-n4;
% n5=(molar_air_stoich*2+oxygen*fuelequiv-2*n1-n2-n3)/2;
% n6=molar_air_stoich*3.76;
% ntsum=n1+n2+n3+n4+n5+n6;
% error=n5ntvalue-n5/ntsum;
%
% n5ntvalue=n5ntvalue-0.000000000000000000001;
% if n5ntvalue<0
% break
% end
% end
%
% else
% while error>0
% a=sqrt(n5ntvalue*pressure);
% n3=(carbon*fuelequiv)/(kco2*a+1);
% n4=((hydrogen/2)*fuelequiv)/(kh2o*a+1);
% n1=carbon*fuelequiv-n3;
% n2=(hydrogen/2)*fuelequiv-n4;
% n5=(molar_air_stoich*2+oxygen*fueq-2*n1-n2-n3)/2;
% n6=molar_air_stoich*3.76;
% ntsum=n1+n2+n3+n4+n5+n6;
% error=n5ntvalue-n5/ntsum;
%
% n5ntvalue=n5ntvalue-0.000000000000000000001; %-21
% end
% end
% %-------------------------
% %error now negative
out=n5ntvalue;
ItterationsDiss = itt;
end








