Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

To resolve issues starting MATLAB on Mac OS X 10.10 (Yosemite) visit: http://www.mathworks.com/matlabcentral/answers/159016

Model Migration from Code to Simulink: The Execution of Matlab Function Blocks

Asked by Niel on 17 Jul 2013

I am in the process of migrating a Matlab model to SImulink. I do this by using Matlab's function block, as can be seen below:

For this block to work, I first execute an "initialise" m-file, which sends the required values as busses to Simulink's inputs. It works quite well - but when executing the Simulink model (code given below), I receive the following error:

I can find no reason for this error, unless the solver does not execute the function from top to bottom. Any views on rectifying this? The variable given as "undefined" is quite clearly defined and is not in an _if _statement.

The code (I know it is quite hefty, but I need only general comments):

function output=calculate_autoclave_odes(t,X,data,Q_removal,flow_constants)

%This function generates the set of ordinary differential equations that %need to be solved in order to calculate the dynamic behaviour of the %autoclave in the base metal refinery (i.e. the flow rates, compositions, %and temperatures of the respective streams)

residue_solid_fraction=zeros(1,7); spent_concentration = zeros(1,7); m_Liquid_in=zeros(1,8); m_solids_in=zeros(1,8); n_phases_in=zeros(1,14); m_phases_in=zeros(1,14); m_Liquid_TK50=zeros(1,8); P_water_2nd_stage = 0; Cp_L_9=0; extent_2 = zeros(1,21); %Initialise the vector containing the extents of reaction

%Assign input data to local variables %Composition of first stage solid residue residue_solid_fraction(1) = data.residue.Cu/100; %Mass fraction Cu in first stage residue residue_solid_fraction(2) = data.residue.Ni/100; %Mass fraction Ni in first stage residue residue_solid_fraction(3) = data.residue.Fe/100; %Mass fraction Fe in first stage residue residue_solid_fraction(4) = data.residue.S/100; %Mass fraction S in first stage residue residue_solid_fraction(5) = data.residue.Rh/1000000; %Mass fraction Rh in first stage residue residue_solid_fraction(6) = data.residue.Ru/1000000; %Mass fraction Ru in first stage residue residue_solid_fraction(7) = data.residue.Ir/1000000; %Mass fraction Ir in first stage residue Cu_as_Cu9S5 = data.residue.Cu_as_Cu9S5; %Fraction of copper in first stage residue present as Cu9S5, remainder as CuS Ni_as_NiS = data.residue.Ni_as_NiS; %Fraction of nickel in first stage residue present as NiS, remainder as Ni3S4 Rh_as_Rh2S3 = data.residue.Rh_as_Rh2S3; %Fraction of rhodium in first stage residue present as Rh2S3, remainder as Rh Ru_as_RuS2 = data.residue.Ru_as_RuS2; %Fraction of ruthenium in first stage residue present as RuS2, remainder as Ru Ir_as_Ir2S3 = data.residue.Ir_as_Ir2S3; %Fraction of iridium in first stage residue present as Ir2S3, remainder as Ir Fe_as_FeOHSO4 = data.residue.Fe_as_FeOHSO4; %Fraction of iron in first stage residue present as leachable Fe(OH)SO4, remainder not leachable iron species CuS_leachable = data.residue.CuS_leachable;

    %Composition of spent electrolyte
    spent_concentration(1) = data.spent.Cu;                 %Concentration of copper in spent electrolyte, g/l
    spent_concentration(2) = data.spent.Ni;                 %Concentration of nickel in spent electrolyte, g/l
    spent_concentration(3) = data.spent.Fe/1000;            %Concentration of iron in spent electrolyte, g/l
    spent_concentration(4) = data.spent.Rh/1000;            %Concentration of rhodium in spent electrolyte, g/l
    spent_concentration(5) = data.spent.Ru/1000;            %Concentration of ruthenium in spent electrolyte, g/l
    spent_concentration(6) = data.spent.Ir/1000;            %Concentration of iridium in spent electrolyte, g/l
    spent_concentration(7) = data.spent.Sulphuric_acid;     %Concentration of sulphuric acid in spent electrolyte, g/l    
    %Component properties (densities, heat capacities, molecular weights etc.)
    density_acid = data.acid.density;                           %Density of acid, kg/l
    density_spent = data.spent.density;                         %Density of spent electrolyte, kg/l
    density_residue = data.residue.density;                     %Density of first stage residue, kg/l
    acid_concentration = data.acid.concentration;               %Mass percentage H2SO4 in acid
    MW = data.constants.MW;                                     %Data set containing the molecular weights of individual elements, g/mol
    Cp_H2SO4.a = 1000*data.constants.Cp.H2SO4_a/(2*MW.H+MW.S+4*MW.O);   %Heat capacity constant of sulphuric acid, Cp = a+bT, T in dC, Cp in kJ/kg.dC
    Cp_H2SO4.b = 1000*data.constants.Cp.H2SO4_b/(2*MW.H+MW.S+4*MW.O);   %Heat capacity constant of sulphuric acid, Cp = a+bT, T in dC, Cp in kJ/kg.dC
    Cp_H2O = 1000*data.constants.Cp.H2O/(2*MW.H+MW.O);                  %Heat capacity of water kJ/kg.dC
    Cp_steam.a = 1000*data.constants.Cp.Steam_a/(2*MW.H+MW.O);          %Heat capacity of steam, Cp = a+bT+cT^2+dT^3, T in dC, Cp in kJ/kg.dC
    Cp_steam.b = 1000*data.constants.Cp.Steam_b/(2*MW.H+MW.O);          %Heat capacity of steam, Cp = a+bT+cT^2+dT^3, T in dC, Cp in kJ/kg.dC
    Cp_steam.c = 1000*data.constants.Cp.Steam_c/(2*MW.H+MW.O);          %Heat capacity of steam, Cp = a+bT+cT^2+dT^3, T in dC, Cp in kJ/kg.dC
    Cp_steam.d = 1000*data.constants.Cp.Steam_d/(2*MW.H+MW.O);          %Heat capacity of steam, Cp = a+bT+cT^2+dT^3, T in dC, Cp in kJ/kg.dC
    Cp_O2.a = 1000*data.constants.Cp.Oxygen_a/(2*MW.O);                 %Heat capacity of oxygen, Cp = a+bT+cT^2+dT^3, T in dC, Cp in kJ/kg.dC 
    Cp_O2.b = 1000*data.constants.Cp.Oxygen_b/(2*MW.O);                 %Heat capacity of oxygen, Cp = a+bT+cT^2+dT^3, T in dC, Cp in kJ/kg.dC
    Cp_O2.c = 1000*data.constants.Cp.Oxygen_c/(2*MW.O);                 %Heat capacity of oxygen, Cp = a+bT+cT^2+dT^3, T in dC, Cp in kJ/kg.dC
    Cp_O2.d = 1000*data.constants.Cp.Oxygen_d/(2*MW.O);                 %Heat capacity of oxygen, Cp = a+bT+cT^2+dT^3, T in dC, Cp in kJ/kg.dC
    Cp_NiS = data.constants.Cp.NiS/(MW.Ni+MW.S);                        %Heat capacity of NiS, kJ/kg.K
    Cp_CuS = data.constants.Cp.CuS/(MW.Cu+MW.S);                        %Heat capacity of CuS, kJ/kg.K
    Cp_Cu9S5 = data.constants.Cp.Cu9S5/(9*MW.Cu+5*MW.S);                %Heat capacity of Cu9S5, kJ/kg.K
    Cp_Ni3S4.a = data.constants.Cp.Ni3S4_a/(3*MW.Ni+4*MW.S);            %Heat capacity of Ni3S4, Cp = a+bT+cT^2+dT^3+e/(T^2), T in K, Cp in kJ/kg.K
    Cp_Ni3S4.b = data.constants.Cp.Ni3S4_b/(3*MW.Ni+4*MW.S);            %Heat capacity of Ni3S4, Cp = a+bT+cT^2+dT^3+e/(T^2), T in K, Cp in kJ/kg.K
    Cp_Ni3S4.c = data.constants.Cp.Ni3S4_c/(3*MW.Ni+4*MW.S);            %Heat capacity of Ni3S4, Cp = a+bT+cT^2+dT^3+e/(T^2), T in K, Cp in kJ/kg.K
    Cp_Ni3S4.d = data.constants.Cp.Ni3S4_d/(3*MW.Ni+4*MW.S);            %Heat capacity of Ni3S4, Cp = a+bT+cT^2+dT^3+e/(T^2), T in K, Cp in kJ/kg.K
    Cp_Ni3S4.e = data.constants.Cp.Ni3S4_e/(3*MW.Ni+4*MW.S);            %Heat capacity of Ni3S4, Cp = a+bT+cT^2+dT^3+e/(T^2), T in K, Cp in kJ/kg.K
    %Heats of state changes (heats of reactions, heat of evaporation)
    H_Evap_H2O = data.constants.H_Evap_H2O;                             %Heat of vaporisation for water at 100dC, kJ/kg
    H_Reaction = data.constants.H_reaction;                             %Data set containing the heat of reaction for the respective reactions, kJ/mol
    %Stream information (flow rates, temperatures)
    m_S_1 = data.flow.stream_1.mass_flow_solids;                %Mass flow rate of solids in stream 1, kg/h 
    V_spent_1 = data.flow.stream_1.volume_flow_spent;           %Volumetric flow rate of spent electrolyte in stream 1, l/h  
    V_2 = data.flow.stream_2.volume_flow;                       %Volumetric flow rate of stream 2 (spent electrolyte), l/h 
    V_3 = data.flow.stream_3.volume_flow;                       %Volumetric flow rate of stream 3 (water), l/h
    V_4 = data.flow.stream_4.volume_flow;                       %Volumetric flow rate of stream 4 (acid), l/h
    m_9 = data.flow.stream_9.mass_flow;                         %Mass flow rate of stream 9, kg/h
    T_5 = data.flow.stream_5.temperature;                       %Temperature of the stream 5, dC
    T_21 = data.flow.stream_21.temperature;
    fraction_oxygen_excess = data.flow.oxygen.total_excess;     %Fraction excess oxygen fed to autoclave                
    fraction_oxygen_stream_10 = data.flow.oxygen.fraction_10;   %Fraction of total oxygen fed in stream 10                    
    fraction_oxygen_stream_11 = data.flow.oxygen.fraction_11;   %Fraction of total oxygen fed in stream 11                    
    V_18 = data.flow.stream_18.volume;                          %Volumetric flow rate of spent electrolyte in stream 18, l/h 
    V_19 = data.flow.stream_19.volume;                          %Volumetric flow rate of water in stream 19, l/h
    V_20 = data.flow.stream_20.volume;                          %Volumetric flow rate of acid in stream 20, l/h    
    Q_Removed_2 = Q_removal.compartment_2;                      %Rate at which heat is removed from autoclave compartment 2, kJ/h
    Q_Removed_3 = Q_removal.compartment_3;                      %Rate at which heat is removed from autoclave compartment 3,kJ/h
    m_13 = Q_removal.m_13;                                      %Mass flow rate of steam to compartment 3, kg/h
    %Autoclave conditions (pressure, volume, etc.)
    autoclave_P = data.equipment.autoclave_pressure;                    %Operating pressure of the autoclave, bar
    oxygen_T = data.flow.oxygen_temperature;                            %Oxygen temperature, dC
    steam_T = data.flow.steam_temperature;                              %Steam temperature, dC
    volume_1 = 1000*data.equipment.autoclave_volume.compartment_1;      %Volume of first autoclave compartment in l (second stage leach)
    volume_2 = 1000*data.equipment.autoclave_volume.compartment_2;      %Volume of second autoclave compartment in l (second stage leach) 
    volume_3 = 1000*data.equipment.autoclave_volume.compartment_3;      %Volume of third autoclave compartment in l (second stage leach)
    volume_4 = 1000*data.equipment.autoclave_volume.compartment_4;      %Volume of fourth autoclave compartment in l (third stage leach)
    agitator_power = 3600*data.equipment.agitator_power;                %Power of agitator drive, kJ/h
    autoclave_outer_diameter = data.equipment.autoclave.outer_diameter;             %Outside diameter of autoclave, m
    autoclave_surface_temperature = data.equipment.autoclave.outer_temperature;     %Temperature of autoclave outer surface, dC
    autoclave_emissivity = data.equipment.autoclave.emissivity;                     %Emissivity of autoclave outer surface  
    length_1 = data.equipment.autoclave.compartment_1_length;                       %Length of the first compartment, m
    length_2 = data.equipment.autoclave.compartment_2_length;                       %Length of the second compartment, m
    length_3 = data.equipment.autoclave.compartment_3_length;                       %Length of the third compartment, m
    length_4 = data.equipment.autoclave.compartment_4_length;                       %Length of the fourth compartment, m
    air_temperature = data.constants.air.temperature;               %Ambient temperature, dC
    air_k = data.constants.air.k;                                   %Thermal conductivity of air, W/m.K
    air_viscosity = data.constants.air.viscosity;                   %Viscosity of air, m^2/s
    air_thermal_diff = data.constants.air.thermal_diff;             %Thermal diffusivity of air, m^2/s
    air_Pr = data.constants.air.Pr;                                 %Prandtl number of air at 300K
    air_thermal_expansion = data.constants.air.thermal_expansion;   %Thermal expansion coefficient, K^-1
    %Estimate the constants to relate the outlet flow rate from a tank to the mass of slurry in the tank
    prep_tank_constant = flow_constants.prep_tank_2;
    recycle_tank_constant = flow_constants.recycle_tank;
    discharge_tank_constant = flow_constants.discharge_tank;
    thickener_constant = flow_constants.thickener;
    prep_tank_3_constant = flow_constants.prep_tank_3;
    %Calculation parameters
    max_loops = 12000;                          %Maximum number of iterative calculation loops allowed
    tolerance_percent = 0.0001;                 %Maximum percentage by which calculated temperature of composition may vary from initial guess 

%Specify step changes in process conditions at specific time instances by varying any of the operating variables above (THIS MUST BE REPLICATED IN THE FUNCTION 'calculate_dynamic_behaviour')

% if t<6; % V_4 = V_4; % elseif t<72; % V_4 = 1.5*V_4; % else % V_4 = V_4; % end;

% if t==0; % m_13 = m_13; % else % m_13 = m_13+10; % end;

for i=(1:length(X)) if (X(i)<0) X(i) =0; end; end;

%Calculate flow rates of streams fed to second stage slurry preparation tank, TK-10 m_spent_1 = V_spent_1*density_spent; %Mass flow rate of spent electrolyte in stream 1, kg/h m_2 = V_2*density_spent; %Mass flow rate of spent electrolyte in stream 2, kg/h m_3 = V_3; %Mass flow rate of water in stream 3, kg/h m_4 = V_4*density_acid; %Mass flow rate of stream 4 (acid), kg/h m_H2SO4_4 = (acid_concentration/100)*m_4; %Mass flow rate of pure H2SO4 in stream 4, kg/h

%Calculate the total flow rates of liquids and dissolved species entering the second stage slurry preparation tank, TK-10, kg/h m_Liquid_in(1) = m_3 + m_2 + m_spent_1 + m_4; %Mass balance to calculate total flow rate of aqueous part entering TK-10, kg/h m_Liquid_in(2:7) = (V_2+V_spent_1)*spent_concentration(1:6)/1000; %Mass balance to calculate flow rate of dissolved metals entering TK-10, kg/h m_Liquid_in(8) = m_H2SO4_4 + (V_2+V_spent_1)*spent_concentration(7)/1000; %Mass balance to calculate flow rate of acid entering TK-10, kg/h

%Calculate the flow rates of solid elemental species entering the second stage slurry preparation tank, TK-10, kg/h m_solids_in(1) = m_S_1; %Total solids into preparation tank, TK-10, kg/h m_solids_in(2:8) = m_solids_in(1)*residue_solid_fraction; %Calculate flow rate of individual elements;composition of solids assumed to remain unchanged in the preparation tank 400TK-010, kg/h

%Calculate the molar amounts of the solid phases entering the second stage slurry preparation tank, TK-10 (NiS, Ni3S4, Cu9S5, CuS, Fe(OH)SO4, Rh2S3, Rh, RhO2, RuS2, Ru, RuO2, Ir2S3, Ir, IrO2) n_phases_in(1) = (m_solids_in(3)/MW.Ni)*Ni_as_NiS; %Molar flow rate of NiS in stream 5, kmol/h n_phases_in(2) = ((m_solids_in(3)/MW.Ni)*(1-Ni_as_NiS))/3; %Molar flow rate of Ni3S4 in stream 5, kmol/h n_phases_in(3) = ((m_solids_in(2)/MW.Cu)*Cu_as_Cu9S5)/9; %Molar flow rate of Cu9S5 in stream 5, kmol/h n_phases_in(4) = (m_solids_in(2)/MW.Cu)*(1-Cu_as_Cu9S5); %Molar flow rate of CuS in stream 5, kmol/h n_phases_in(5) = (m_solids_in(4)/MW.Fe)*Fe_as_FeOHSO4; %Molar flow rate of FeOHSO4 in stream 5, kmol/h n_phases_in(6) = ((m_solids_in(6)/MW.Rh)*Rh_as_Rh2S3)/2; %Molar flow rate of Rh2S3 in stream 5, kmol/h n_phases_in(7) = (m_solids_in(6)/MW.Rh)*(1-Rh_as_Rh2S3); %Molar flow rate of Rh in stream 5, kmol/h n_phases_in(8) = 0; %Molar flow rate of RhO2 in stream 5, kmol/h n_phases_in(9) = ((m_solids_in(7)/MW.Ru)*Ru_as_RuS2); %Molar flow rate of RuS2 in stream 5, kmol/h n_phases_in(10) = (m_solids_in(7)/MW.Ru)*(1-Ru_as_RuS2); %Molar flow rate of Ru in stream 5, kmol/h n_phases_in(11) = 0; %Molar flow rate of RuO2 in stream 5, kmol/h n_phases_in(12) = ((m_solids_in(8)/MW.Ir)*Ir_as_Ir2S3)/2; %Molar flow rate of Ir2S3 in stream 5, kmol/h n_phases_in(13) = (m_solids_in(8)/MW.Ir)*(1-Ir_as_Ir2S3); %Molar flow rate of Ir in stream 5, kmol/h n_phases_in(14) = 0; %Molar flow rate of IrO2 in stream 5, kmol/h

%Calculate the mass amounts of the solid phases entering the second stage slurry preparation tank, TK-10 (NiS, Ni3S4, Cu9S5, CuS, Fe(OH)SO4, Rh2S3, Rh, RhO2, RuS2, Ru, RuO2, Ir2S3, Ir, IrO2) m_phases_in(1) = n_phases_in(1)*(MW.Ni+MW.S); %Mass flow rate of NiS in stream 5, kg/h m_phases_in(2) = n_phases_in(2)*(3*MW.Ni+4*MW.S); %Mass flow rate of Ni3S4 in stream 5, kg/h m_phases_in(3) = n_phases_in(3)*(9*MW.Cu+5*MW.S); %Mass flow rate of Cu9S5 in stream 5, kg/h m_phases_in(4) = n_phases_in(4)*(MW.Cu+MW.S); %Mass flow rate of CuS in stream 5, kg/h m_phases_in(5) = n_phases_in(5)*(MW.Fe+5*MW.O+MW.H+MW.S); %Mass flow rate of FeOHSO4 in stream 5, kg/h m_phases_in(6) = n_phases_in(6)*(2*MW.Rh+3*MW.S); %Mass flow rate of Rh2S3 in stream 5, kg/h m_phases_in(7) = n_phases_in(7)*(MW.Rh); %Mass flow rate of Rh in stream 5, kg/h m_phases_in(8) = n_phases_in(8)*(MW.Rh+2*MW.O); %Mass flow rate of RhO2 in stream 5, kg/h m_phases_in(9) = n_phases_in(9)*(MW.Ru+2*MW.S); %Mass flow rate of RuS2 in stream 5, kg/h m_phases_in(10) = n_phases_in(10)*(MW.Ru); %Mass flow rate of Ru in stream 5, kg/h m_phases_in(11) = n_phases_in(11)*(MW.Ru+2*MW.O); %Mass flow rate of RuO2 in stream 5, kg/h m_phases_in(12) = n_phases_in(12)*(2*MW.Ir+3*MW.S); %Mass flow rate of Ir2S3 in stream 5, kg/h m_phases_in(13) = n_phases_in(13)*(MW.Ir); %Mass flow rate of Ir in stream 5, kg/h m_phases_in(14) = n_phases_in(14)*(MW.Ir+2*MW.O); %Mass flow rate of IrO2 in stream 5, kg/h

%Calculate the total flow rates of liquids and dissolved species entering the third stage slurry preparation tank, TK-50, kg/h m_18 = V_18*density_spent; %Mass flow rate of spent electrolyte in stream 2, kg/h m_19 = V_19; %Mass flow rate of water in stream 3, kg/h m_20 = V_20*density_acid; %Mass flow rate of stream 4 (acid), kg/h m_H2SO4_20 = (acid_concentration/100)*m_20; %Mass flow rate of pure H2SO4 in stream 4, kg/h

%Calculate the total flow rates of liquids and dissolved species entering the third stage slurry preparation tank, TK-50, kg/h m_Liquid_TK50(1) = m_18 + m_19 + m_20; %Mass balance to calculate total flow rate of aqueous part entering TK-50, kg/h m_Liquid_TK50(2:7) = (V_18)*spent_concentration(1:6)/1000; %Mass balance to calculate flow rate of dissolved metals entering TK-50, kg/h m_Liquid_TK50(8) = m_H2SO4_20 + (V_18)*spent_concentration(7)/1000; %Mass balance to calculate flow rate of acid entering TK-50, kg/h

%Calculate the fraction oxygen in the vapour space of the autoclave (assume that the vapour phases in the different compartments have the same temperature, equal to the temperature in the first compartment) %Calculate water vapour pressure (total pressure equal to water vapour pressure plus oxygen partial pressure, assuming that vapour phase is saturated with water vapour) if X(73)<60 %Temperatures less than 60dC require different constants for Antoine equation P_water_2nd_stage = ((10^(8.10765-1750.286/(X(73)+235)))/760)*1.01325; %Vapour pressure of water, bar else P_water_2nd_stage = ((10^(7.96681-1668.210/(X(73)+228)))/760)*1.01325; %Vapour pressure of water,bar end; %Calculate mol fraction oxygen in vapour y_O2_2nd_stage = 1-P_water_2nd_stage/autoclave_P; %Mole fraction of oxygen in vapour of second stage leach, assuming that vapour is saturated with water vapour

%Estimate heat losses from respective autoclave compartments due to convection and radiation Ra = (9.81*air_thermal_expansion*(autoclave_surface_temperature-air_temperature)*(autoclave_outer_diameter^3))/(air_viscosity*air_thermal_diff); Nu = (0.6+(0.387*Ra^(1/6))/((1+(0.559/air_Pr)^(9/16))^(8/27)))^2; h = air_k*Nu/autoclave_outer_diameter; Q_loss = (h*22/7*autoclave_outer_diameter*(autoclave_surface_temperature-air_temperature)+autoclave_emissivity*22/7*autoclave_outer_diameter*5.67*(10^-8)*((autoclave_surface_temperature+273.15)^4-(air_temperature+273.15)^4))*3600/1000;%Heat loss per meter of autoclave, kJ/h.m

%Calculate the theoretical amount of oxygen required for complete dissolution of NiS and Cu9S5 in feed to process (assuming all copper present as Cu9S5 and all nickel present at NiS) n_Cu_solids_in = m_solids_in(2)/MW.Cu; %Molar flow rate of solid Cu in stream 5, kmol/h n_Ni_solids_in = m_solids_in(3)/MW.Ni; %Molar flow rate of solid Ni in stream 5, kmol/h theoretical_oxygen = 2*n_Ni_solids_in + 12*(n_Cu_solids_in)/9; %Theoretical amount of oxygen required for complete leaching of NiS and Cu9S5, kmol/h

%Calculate actual flow rates of the oxygen feed streams n_10 = theoretical_oxygen*fraction_oxygen_excess*fraction_oxygen_stream_10; %Flow rate of stream 10, kmol/h n_11 = theoretical_oxygen*fraction_oxygen_excess*fraction_oxygen_stream_11; %Flow rate of stream 11, kmol/h n_12 = theoretical_oxygen*fraction_oxygen_excess-n_10-n_11; %Flow rate of stream 12, kmol/h m_10 = n_10*(2*MW.O); %Oxygen mass flow rate in stream 10, kg/h m_11 = n_11*(2*MW.O); %Oxygen mass flow rate in stream 11, kg/h m_12 = n_12*(2*MW.O); %Oxygen mass flow rate in stream 12, kg/h

%Estimate the heat capacities of the liquid phase in the respective streams, kJ/kg.dC Cp_L_5 = X(24)*Cp_H2SO4.a + (1-X(24))*Cp_H2O; %Stream 5 Cp_L_7 = X(48)*Cp_H2SO4.a+(1-X(48))*Cp_H2O; %Stream 7 Cp_L_9 = X(72)*Cp_H2SO4.a+(1-X(72))*Cp_H2O; %Stream 9 Cp_L_AC1 = Cp_L_9; %Stream AC1 Cp_L_AC2 = X(96)*Cp_H2SO4.a+(1-X(96))*Cp_H2O; %Stream AC2 Cp_L_21 = X(193)*Cp_H2SO4.a+(1-X(193))*Cp_H2O; %Stream 21 Cp_L_22 = X(216)*Cp_H2SO4.a+(1-X(216))*Cp_H2O; %Stream 22 Cp_L_14 = X(120)*Cp_H2SO4.a+(1-X(120))*Cp_H2O; %Stream 14

%Normalise mass fractions of NiS, Ni3S4, CuS, and Cu9S5 in solid phase in the respective streams Normalx_NiS_S_5 = X(3)/(X(3)+X(4)+X(5)+X(6)); Normalx_Ni3S4_S_5 = X(4)/(X(3)+X(4)+X(5)+X(6)); Normalx_Cu9S5_S_5 = X(5)/(X(3)+X(4)+X(5)+X(6)); Normalx_CuS_S_5 = X(6)/(X(3)+X(4)+X(5)+X(6));

    Normalx_NiS_S_7 = X(27)/(X(27)+X(28)+X(29)+X(30));
    Normalx_Ni3S4_S_7 = X(28)/(X(27)+X(28)+X(29)+X(30));
    Normalx_Cu9S5_S_7 = X(29)/(X(27)+X(28)+X(29)+X(30));
    Normalx_CuS_S_7 = X(30)/(X(27)+X(28)+X(29)+X(30));
    Normalx_NiS_S_9 = X(51)/(X(51)+X(52)+X(53)+X(54));
    Normalx_Ni3S4_S_9 = X(52)/(X(51)+X(52)+X(53)+X(54));
    Normalx_Cu9S5_S_9 = X(53)/(X(51)+X(52)+X(53)+X(54));
    Normalx_CuS_S_9 = X(54)/(X(51)+X(52)+X(53)+X(54));
    Normalx_NiS_S_AC2 = X(75)/(X(75)+X(76)+X(77)+X(78));
    Normalx_Ni3S4_S_AC2 = X(76)/(X(75)+X(76)+X(77)+X(78));
    Normalx_Cu9S5_S_AC2 = X(77)/(X(75)+X(76)+X(77)+X(78));
    Normalx_CuS_S_AC2 = X(78)/(X(75)+X(76)+X(77)+X(78));
    Normalx_NiS_S_21 = X(172)/(X(172)+X(173)+X(174)+X(175));
    Normalx_Ni3S4_S_21 = X(173)/(X(172)+X(173)+X(174)+X(175));
    Normalx_Cu9S5_S_21 = X(174)/(X(172)+X(173)+X(174)+X(175));
    Normalx_CuS_S_21 = X(175)/(X(172)+X(173)+X(174)+X(175));
    Normalx_NiS_S_22 = X(195)/(X(195)+X(196)+X(197)+X(198));
    Normalx_Ni3S4_S_22 = X(196)/(X(195)+X(196)+X(197)+X(198));
    Normalx_Cu9S5_S_22 = X(197)/(X(195)+X(196)+X(197)+X(198));
    Normalx_CuS_S_22 = X(198)/(X(195)+X(196)+X(197)+X(198));
    Normalx_NiS_S_14 = X(99)/(X(99)+X(100)+X(101)+X(102));
    Normalx_Ni3S4_S_14 = X(100)/(X(99)+X(100)+X(101)+X(102));
    Normalx_Cu9S5_S_14 = X(101)/(X(99)+X(100)+X(101)+X(102));
    Normalx_CuS_S_14 = X(102)/(X(99)+X(100)+X(101)+X(102));

%Estimate heat capacity of the solid portion of the respective streams, kJ/kg.K Cp_S_5 = Normalx_NiS_S_5*Cp_NiS + Normalx_Cu9S5_S_5*Cp_Cu9S5 + Normalx_CuS_S_5*Cp_CuS + Normalx_Ni3S4_S_5*Cp_Ni3S4.a; Cp_S_7 = Normalx_NiS_S_7*Cp_NiS + Normalx_Cu9S5_S_7*Cp_Cu9S5 + Normalx_CuS_S_7*Cp_CuS + Normalx_Ni3S4_S_7*Cp_Ni3S4.a; Cp_S_9 = Normalx_NiS_S_9*Cp_NiS + Normalx_Cu9S5_S_9*Cp_Cu9S5 + Normalx_CuS_S_9*Cp_CuS + Normalx_Ni3S4_S_9*Cp_Ni3S4.a; Cp_S_AC1 = Cp_S_9; Cp_S_AC2 = Normalx_NiS_S_AC2*Cp_NiS + Normalx_Cu9S5_S_AC2*Cp_Cu9S5 + Normalx_CuS_S_AC2*Cp_CuS + Normalx_Ni3S4_S_AC2*Cp_Ni3S4.a; Cp_S_21 = Normalx_NiS_S_21*Cp_NiS + Normalx_Cu9S5_S_21*Cp_Cu9S5 + Normalx_CuS_S_21*Cp_CuS + Normalx_Ni3S4_S_21*Cp_Ni3S4.a; Cp_S_22 = Normalx_NiS_S_22*Cp_NiS + Normalx_Cu9S5_S_22*Cp_Cu9S5 + Normalx_CuS_S_22*Cp_CuS + Normalx_Ni3S4_S_22*Cp_Ni3S4.a; Cp_S_14 = Normalx_NiS_S_14*Cp_NiS + Normalx_Cu9S5_S_14*Cp_Cu9S5 + Normalx_CuS_S_14*Cp_CuS + Normalx_Ni3S4_S_14*Cp_Ni3S4.a;

%Calculate energy removal to cool liquid portion of stream 9 to 100dC, kJ/h Q_removed_L_9 = (m_9*X(65))*Cp_L_9*(X(73)-100);

%Calculate energy removal to cool solid portion of stream 9 to 100dC, kJ/h Q_removed_S_9 = (m_9*X(50))*Cp_S_9*(X(73)-100);

%Calculate total energy removal to cool stream 9 to 100dC, kJ/h Q_removal_9 = Q_removed_L_9+Q_removed_S_9;

%Calculate rate of water evaporation from stream 9 upon entering the flash recycle tank, kg/h %If recycle stream has a temperature less than 100dC, no evaporation will occur if (Q_removal_9>0) m_H2O_flash_evaporated = Q_removal_9/H_Evap_H2O; else m_H2O_flash_evaporated = 0; end; %If the recycle flow rate is too low (or the autoclave operating temperating temperature too high), all the water in the recycle stream will evaporate if (m_H2O_flash_evaporated>m_9*X(65)*(1-X(66)*(MW.Cu+MW.S+4*MW.O)/MW.Cu-X(67)*(MW.Ni+MW.S+4*MW.O)/MW.Ni-X(72))) m_H2O_flash_evaporated = m_9*X(65)*(1-X(66)*(MW.Cu+MW.S+4*MW.O)/MW.Cu-X(67)*(MW.Ni+MW.S+4*MW.O)/MW.Ni-X(72)); temp_flag = 1; %Raise indication flag that all the water in the recycle stream evaporates during flashing under these conditions else temp_flag = 0; end;

%Compartment 1 reaction rates and extents of reactions

    %Estimate volumetric flow rate of the feed to the autoclave, stream 5
        V_S_5 = X(2)*prep_tank_constant*X(1)^0.5/density_residue;                   %Volumetric flow rate of solids in stream 5, l/h
        V_L_5 = X(17)*prep_tank_constant*X(1)^0.5/density_spent;                    %Volumetric flow rate of liquid in stream 5, l/h
        V_5 = V_S_5 + V_L_5;                                                        %Total volumetric flow rate of stream 5, l/h
    %Estimate volumetric flow rate of the feed to the autoclave, stream 7
        V_S_7 = X(26)*recycle_tank_constant*X(25)^0.5/density_residue;              %Volumetric flow rate of solids in stream 7, l/h
        V_L_7 = X(41)*recycle_tank_constant*X(25)^0.5/density_spent;                %Volumetric flow rate of liquid in stream 7, l/h
        V_7 = V_S_7 + V_L_7;                                                        %Total volumetric flow rate of stream 7, l/h
    %Estimate the volumetric flow rate of stream 9
        V_S_9 = m_9*X(50)/density_residue;                                          %Volumetric flow rate of solids in stream 9, l/h
        V_L_9 = m_9*X(65)/density_spent;                                            %Volumetric flow rate of liquid in stream 9, l/h
        V_9 = V_S_9 + V_L_9;                                                        %Total volumetric flow rate of stream 9, l/h
    %Calculate the concentrations of reagents present in compartment 1 in order to calculate the reaction rates
        %Dissolved/liquid species
        liq_fractions_1.H2SO4 = X(72);       %Concentration of aqueous species in autoclave compartment 1, kg H2SO4/kg liquid
        liq_fractions_1.Rh = X(69);          %Concentration of aqueous species in autoclave compartment 1, kg Rh/kg liquid
        liq_fractions_1.Ru = X(70);          %Concentration of aqueous species in autoclave compartment 1, kg Ru/kg liquid
        liq_fractions_1.Ir = X(71);          %Concentration of aqueous species in autoclave compartment 1, kg Ir/kg liquid
        %Solid species entering compartment 1
        solids_in_1.NiS = X(26)*recycle_tank_constant*X(25)^0.5*X(27)/V_7;          %Concentration of solids entering the first autoclave compartment, kg NiS/l slurry
        solids_in_1.Ni3S4 = X(26)*recycle_tank_constant*X(25)^0.5*X(28)/V_7;        %Concentration of solids entering the first autoclave compartment, kg Ni3S4/l slurry
        solids_in_1.Cu9S5 = X(26)*recycle_tank_constant*X(25)^0.5*X(29)/V_7;        %Concentration of solids entering the first autoclave compartment, kg Cu9S5/l slurry
        solids_in_1.FeOHSO4 = X(26)*recycle_tank_constant*X(25)^0.5*X(31)/V_7;      %Concentration of solids entering the first autoclave compartment, kg FeOHSO4/l slurry
        solids_in_1.Rh2S3 = X(26)*recycle_tank_constant*X(25)^0.5*X(32)/V_7;        %Concentration of solids entering the first autoclave compartment, kg Rh2S3/l slurry
        solids_in_1.Rh = X(26)*recycle_tank_constant*X(25)^0.5*X(33)/V_7;           %Concentration of solids entering the first autoclave compartment, kg Rh/l slurry
        solids_in_1.RuS2 = X(26)*recycle_tank_constant*X(25)^0.5*X(35)/V_7;         %Concentration of solids entering the first autoclave compartment, kg RuS2/l slurry
        solids_in_1.Ru = X(26)*recycle_tank_constant*X(25)^0.5*X(36)/V_7;           %Concentration of solids entering the first autoclave compartment, kg Ru/l slurry
        solids_in_1.Ir2S3 = X(26)*recycle_tank_constant*X(25)^0.5*X(38)/V_7;        %Concentration of solids entering the first autoclave compartment, kg Ir2S3/l slurry
        solids_in_1.Ir = X(26)*recycle_tank_constant*X(25)^0.5*X(39)/V_7;           %Concentration of solids entering the first autoclave compartment, kg Ir/l slurry
        Cu9S5_fresh = X(2)*prep_tank_constant*X(1)^0.5*X(5)/V_5;                    %Concentration of Cu9S5 in the fresh feed to the flash recycle tank, kg Cu9S5/l slurry
        %Solid species leaving compartment 1
        solids_out_1.NiS = m_9*X(50)*X(51)/V_9;         %Concentration of solids leaving the first autoclave compartment, kg NiS/l slurry
        solids_out_1.Ni3S4 = m_9*X(50)*X(52)/V_9;       %Concentration of solids leaving the first autoclave compartment, kg Ni3S4/l slurry
        solids_out_1.Cu9S5 = m_9*X(50)*X(53)/V_9;       %Concentration of solids leaving the first autoclave compartment, kg Cu9S5/l slurry
        solids_out_1.CuS = m_9*X(50)*X(54)/V_9;         %Concentration of solids leaving the first autoclave compartment, kg CuS/l slurry
        solids_out_1.FeOHSO4 = m_9*X(50)*X(55)/V_9;     %Concentration of solids leaving the first autoclave compartment, kg FeOHSO4/l slurry
        solids_out_1.Rh2S3 = m_9*X(50)*X(56)/V_9;       %Concentration of solids leaving the first autoclave compartment, kg Rh2S3/l slurry
        solids_out_1.Rh = m_9*X(50)*X(57)/V_9;          %Concentration of solids leaving the first autoclave compartment, kg Rh/l slurry
        solids_out_1.RuS2 = m_9*X(50)*X(59)/V_9;        %Concentration of solids leaving the first autoclave compartment, kg RuS2/l slurry
        solids_out_1.Ru = m_9*X(50)*X(60)/V_9;          %Concentration of solids leaving the first autoclave compartment, kg Ru/l slurry
        solids_out_1.Ir2S3 = m_9*X(50)*X(62)/V_9;       %Concentration of solids leaving the first autoclave compartment, kg Ir2S3/l slurry
        solids_out_1.Ir = m_9*X(50)*X(63)/V_9;          %Concentration of solids leaving the first autoclave compartment, kg Ir/l slurry
        solids_out_1.RhO2 = m_9*X(50)*X(58)/V_9;        %Concentration of solids leaving the first autoclave compartment, kg RhO2/l slurry
        solids_out_1.RuO2 = m_9*X(50)*X(61)/V_9;        %Concentration of solids leaving the first autoclave compartment, kg RuO2/l slurry
        solids_out_1.IrO2 = m_9*X(50)*X(64)/V_9;        %Concentration of solids leaving the first autoclave compartment, kg IrO2/l slurry
        %Oxygen solubility (mol/kg) in stream 9
        Cu_molal_9 = 1000*(X(66)/(1-X(66)*(MW.Cu+MW.S+4*MW.O)/MW.Cu-X(67)*(MW.Ni+MW.S+4*MW.O)/MW.Ni-X(72)))/MW.Cu;                      %Molal concentration of copper in stream 9, mol Cu per kg solvent
        Ni_molal_9 = 1000*(X(67)/(1-X(66)*(MW.Cu+MW.S+4*MW.O)/MW.Cu-X(67)*(MW.Ni+MW.S+4*MW.O)/MW.Ni-X(72)))/MW.Ni;                      %Molal concentration of nickel in stream 9, mol Ni per kg solvent
        H2SO4_molal_9 = 1000*(X(72)/(1-X(66)*(MW.Cu+MW.S+4*MW.O)/MW.Cu-X(67)*(MW.Ni+MW.S+4*MW.O)/MW.Ni-X(72)))/(2*MW.H+MW.S+4*MW.O);    %Molal concentration of H2SO4 in stream 9, mol H2SO4 per kg solvent
        O2_solubility_9 = (2*MW.O/1000)*(calculate_oxygen_solubility(X(73),(100000*autoclave_P),Cu_molal_9,Ni_molal_9,H2SO4_molal_9));  %Oxygen solubility in stream 9, kg O2 per kg water
        molal_O2_solubility_9 = (O2_solubility_9*1000)/(2*MW.O);                                                                        %Oxygen solubility in stream 9, mol O2 per kg water
    %Calculate the reaction rates in compartment 1, mol/l.min
        r_1 = calculate_reaction_rates(X(73),liq_fractions_1,solids_in_1,solids_out_1,molal_O2_solubility_9,Cu9S5_fresh,data);
    %Calculate the amount of solid species consumed, assuming that all reagents are present in excess (no species is completely converted)
        solid_consumption_1 = zeros(1,14);                                                                      %Initialise the vector containing the rate of solid species consumption, kg/h
        solid_consumption_1(1) = volume_1*60*((MW.Ni+MW.S)/1000)*(4*r_1(1)+r_1(2));                             %Rate of NiS consumption in autoclave compartment 1, kg/h
        solid_consumption_1(2) = volume_1*60*((3*MW.Ni+4*MW.S)/1000)*(-r_1(1)+2*r_1(3)+r_1(8)+r_1(10)+r_1(12)); %Rate of Ni3S4 consumption in autoclave compartment 1, kg/h
        solid_consumption_1(3) = volume_1*60*((9*MW.Cu+5*MW.S)/1000)*(r_1(4)+3*r_1(7)+3*r_1(9)+3*r_1(11));      %Rate of Cu9S5 consumption in autoclave compartment 1, kg/h
        solid_consumption_1(4) = volume_1*60*((MW.Cu+MW.S)/1000)*(-5*r_1(4)+r_1(5));                            %Rate of CuS consumption in autoclave compartment 1, kg/h
        solid_consumption_1(5) = volume_1*60*((MW.Fe+5*MW.O+MW.H+MW.S)/1000)*(r_1(6));                          %Rate of FeOHSO4 consumption in autoclave compartment 1, kg/h
        solid_consumption_1(6) = volume_1*60*((2*MW.Rh+3*MW.S)/1000)*(r_1(13));                                 %Rate of Rh2S3 consumption in autoclave compartment 1, kg/h
        solid_consumption_1(7) = volume_1*60*((MW.Rh)/1000)*(4*r_1(14));                                        %Rate of Rh consumption in autoclave compartment 1, kg/h
        solid_consumption_1(8) = volume_1*60*((MW.Rh+2*MW.O)/1000)*(-8*r_1(7)-2*r_1(8)+2*r_1(15));              %Rate of RhO2 consumption in autoclave compartment 1, kg/h
        solid_consumption_1(9) = volume_1*60*((MW.Ru+2*MW.S)/1000)*(4*r_1(16));                                 %Rate of RuS2 consumption in autoclave compartment 1, kg/h
        solid_consumption_1(10) = volume_1*60*((MW.Ru)/1000)*(4*r_1(17));                                       %Rate of Ru consumption in autoclave compartment 1, kg/h
        solid_consumption_1(11) = volume_1*60*((MW.Ru+2*MW.O)/1000)*(-8*r_1(9)-2*r_1(10)+2*r_1(18));            %Rate of RuO2 consumption in autoclave compartment 1, kg/h
        solid_consumption_1(12) = volume_1*60*((2*MW.Ir+3*MW.S)/1000)*(r_1(19));                                %Rate of Ir2S3 consumption in autoclave compartment 1, kg/h
        solid_consumption_1(13) = volume_1*60*((MW.Ir)/1000)*(4*r_1(20));                                       %Rate of Ir consumption in autoclave compartment 1, kg/h
        solid_consumption_1(14) = volume_1*60*((MW.Ir+2*MW.O)/1000)*(-8*r_1(11)-2*r_1(12)+2*r_1(21));           %Rate of IrO2 consumption in autoclave compartment 1, kg/h
        %If solids are produced rather than consumed, check that amount produced does not exceed stochiometrically possible amount
        for i=(1:14)
            if (solid_consumption_1(i)<0)
                if (i==1 && solid_consumption_1(i)<0)
                    solid_consumption_1(i)=0; 
                elseif (i==2 && solid_consumption_1(i)<(-0.25*(X(27)*X(26)*recycle_tank_constant*X(25)^0.5)/(MW.Ni+MW.S)*(3*MW.Ni+4*MW.S)))
                    solid_consumption_1(i) = (-0.25*(X(27)*X(26)*recycle_tank_constant*X(25)^0.5)/(MW.Ni+MW.S)*(3*MW.Ni+4*MW.S));
                elseif (i==3 && solid_consumption_1(i)<0)
                    solid_consumption_1(i)=0;                
                elseif (i==4 && solid_consumption_1(i)<(-5*X(29)*X(26)*recycle_tank_constant*X(25)^0.5/(9*MW.Cu+5*MW.S)*(MW.Cu+MW.S)))
                    solid_consumption_1(i) = (-5*X(29)*X(26)*recycle_tank_constant*X(25)^0.5/(9*MW.Cu+5*MW.S)*(MW.Cu+MW.S));  
                elseif (i==5 && solid_consumption_1(i)<0)
                    solid_consumption_1(i) = 0;        
                elseif (i==6 && solid_consumption_1(i)<0)
                    solid_consumption_1(i) = 0;        
                elseif (i==7 && solid_consumption_1(i)<0)
                    solid_consumption_1(i) = 0;                
                elseif (i==8 && solid_consumption_1(i)<(-1*X(45)*X(41)*recycle_tank_constant*X(25)^0.5/(MW.Rh)*(MW.Rh+2*MW.O)))
                    solid_consumption_1(i) = (-1*X(45)*X(41)*recycle_tank_constant*X(25)^0.5/(MW.Rh)*(MW.Rh+2*MW.O));                    
                elseif (i==9 && solid_consumption_1(i)<0)
                    solid_consumption_1(i) = 0;                   
                elseif (i==10 && solid_consumption_1(i)<0)
                    solid_consumption_1(i) = 0;
                elseif (i==11 && solid_consumption_1(i)<(-1*X(46)*X(41)*recycle_tank_constant*X(25)^0.5/(MW.Ru)*(MW.Ru+2*MW.O)))
                    solid_consumption_1(i) = (-1*X(46)*X(41)*recycle_tank_constant*X(25)^0.5/(MW.Ru)*(MW.Ru+2*MW.O));
                elseif (i==12 && solid_consumption_1(i)<0)
                    solid_consumption_1(i) = 0;
                elseif (i==13 && solid_consumption_1(i)<0)
                    solid_consumption_1(i) = 0;
                elseif (i==14 && solid_consumption_1(i)<(-1*X(47)*X(41)*recycle_tank_constant*X(25)^0.5/(MW.Ir)*(MW.Ir+2*MW.O)))
                    solid_consumption_1(i) = (-1*X(47)*X(41)*recycle_tank_constant*X(25)^0.5/(MW.Ir)*(MW.Ir+2*MW.O));
                end;
            end;
        end; 
        %Check that solid consumption does not exceed the amount of solids fed to compartment 1
        for i=(1:14)
            %If solid consumption is more than the amount of solids fed,set consumption equal to feed rate (complete conversion of solid species)
            if (solid_consumption_1(i)>X(26)*recycle_tank_constant*X(25)^0.5*X(i+26))   
                solid_consumption_1(i)=X(26)*recycle_tank_constant*X(25)^0.5*X(i+26);          
            end;
            %If CuS is consumed, limit the amount of CuS leached to leachable fraction of the fresh CuS fed (empirical restriction)
            if (i==4 && solid_consumption_1(i)>(CuS_leachable*X(2)*X(6)*prep_tank_constant*X(1)^0.5))
                solid_consumption_1(i) = (CuS_leachable*X(2)*X(6)*prep_tank_constant*X(1)^0.5);
            end;
        end;
        %Calculate the extent of reaction for each reaction in compartment 1(assuming that the solid species are the limiting reactants)
        extent_1 = zeros(1,21);                                                                     %Initialise the vector containing the extents of reaction
        if (r_1(1) == 0)
            extent_1(1) = 0;                                                                        %If the reaction rate is zero, the extent of reaction will be zero
        else
            extent_1(1) = solid_consumption_1(1)/((MW.Ni+MW.S)/1000)*(r_1(1))/(4*r_1(1)+r_1(2));    %Extent of reaction 1, mol/h
        end;
        if (r_1(2) == 0)
            extent_1(2) = 0;
        else
            extent_1(2) = solid_consumption_1(1)/((MW.Ni+MW.S)/1000)*(r_1(2))/(4*r_1(1)+r_1(2));    %Extent of reaction 2, mol/h
        end;
        if (r_1(4) == 0)
            extent_1(4) = 0;
        else
            extent_1(4) = solid_consumption_1(3)/((9*MW.Cu+5*MW.S)/1000)*(r_1(4))/(r_1(4)+3*r_1(7)+3*r_1(9)+3*r_1(11)); %Extent of reaction 4, mol/h
        end;
        if (r_1(7) == 0)
            extent_1(7) = 0;
        else
            extent_1(7) = solid_consumption_1(3)/((9*MW.Cu+5*MW.S)/1000)*(r_1(7))/(r_1(4)+3*r_1(7)+3*r_1(9)+3*r_1(11)); %Extent of reaction 7, mol/h
        end;
        if (r_1(9) == 0)
            extent_1(9) = 0;
        else
            extent_1(9) = solid_consumption_1(3)/((9*MW.Cu+5*MW.S)/1000)*(r_1(9))/(r_1(4)+3*r_1(7)+3*r_1(9)+3*r_1(11)); %Extent of reaction 9, mol/h
        end;
        if (r_1(11) == 0)
            extent_1(11) = 0;
        else
            extent_1(11) = solid_consumption_1(3)/((9*MW.Cu+5*MW.S)/1000)*(r_1(11))/(r_1(4)+3*r_1(7)+3*r_1(9)+3*r_1(11));%Extent of reaction 11, mol/h
        end;
        if (r_1(13) == 0)
            extent_1(13) = 0;
        else
            extent_1(13) = solid_consumption_1(6)/((2*MW.Rh+3*MW.S)/1000);                          %Extent of reaction 13, mol/h
        end;
        if (r_1(14) == 0)
            extent_1(14) = 0;
        else
            extent_1(14) = solid_consumption_1(7)/((MW.Rh)/1000)/4;                                 %Extent of reaction 14, mol/h
        end;
        if (r_1(16) == 0)
            extent_1(16) = 0;
        else
            extent_1(16) = solid_consumption_1(9)/((MW.Ru+2*MW.S)/1000)/4;                          %Extent of reaction 16, mol/h
        end;
        if (r_1(17) == 0)
            extent_1(17) = 0;
        else
            extent_1(17) = solid_consumption_1(10)/((MW.Ru)/1000)/4;                                %Extent of reaction 17, mol/h
        end;
        if (r_1(19) == 0)
            extent_1(19) = 0;
        else
            extent_1(19) = solid_consumption_1(12)/((2*MW.Ir+3*MW.S)/1000);                         %Extent of reaction 19, mol/h
        end;
        if (r_1(20) == 0)
            extent_1(20) = 0;
        else
            extent_1(20) = solid_consumption_1(13)/((MW.Ir)/1000)/4;                                %Extent of reaction 20, mol/h
        end;
        if (r_1(6) == 0)
            extent_1(6) = 0;
        else
            extent_1(6) = solid_consumption_1(5)/((MW.Fe+5*MW.O+MW.H+MW.S)/1000);                   %Extent of reaction 6, mol/h
        end;
        if (r_1(5) == 0)
            extent_1(5) = 0;
        else
            extent_1(5) = solid_consumption_1(4)/((MW.Cu+MW.S)/1000)+5*extent_1(4);                 %Extent of reaction 5, mol/h
        end;
        if (r_1(3) == 0)
            extent_1(3) = 0;
        else
            extent_1(3) = (solid_consumption_1(2)/((3*MW.Ni+4*MW.S)/1000)+extent_1(1))*r_1(3)/(2*r_1(3)+r_1(8)+r_1(10)+r_1(12));    %Extent of reaction 3, mol/h
        end;
        if (r_1(8) == 0)
            extent_1(8) = 0;
        else
            extent_1(8) = (solid_consumption_1(2)/((3*MW.Ni+4*MW.S)/1000)+extent_1(1))*r_1(8)/(2*r_1(3)+r_1(8)+r_1(10)+r_1(12));    %Extent of reaction 8, mol/h
        end;
        if (r_1(10) == 0)
            extent_1(10) = 0;
        else
            extent_1(10) = (solid_consumption_1(2)/((3*MW.Ni+4*MW.S)/1000)+extent_1(1))*r_1(10)/(2*r_1(3)+r_1(8)+r_1(10)+r_1(12));  %Extent of reaction 10, mol/h
        end;
        if (r_1(12) == 0)
            extent_1(12) = 0;
        else
            extent_1(12) = (solid_consumption_1(2)/((3*MW.Ni+4*MW.S)/1000)+extent_1(1))*r_1(12)/(2*r_1(3)+r_1(8)+r_1(10)+r_1(12));  %Extent of reaction 12, mol/h
        end;
        if (r_1(15) == 0)
            extent_1(15) = 0;
        else
            extent_1(15) = (solid_consumption_1(8)/((MW.Rh+2*MW.O)/1000)+8*extent_1(7)+2*extent_1(8))/2;        %Extent of reaction 15, mol/h
        end;
        if (r_1(18) == 0)
            extent_1(18) = 0;
        else
            extent_1(18) = (solid_consumption_1(11)/((MW.Ru+2*MW.O)/1000)+8*extent_1(9)+2*extent_1(10))/2;      %Extent of reaction 18, mol/h
        end;
        if (r_1(21) == 0)
            extent_1(21) = 0;
        else
            extent_1(21) = (solid_consumption_1(14)/((MW.Ir+2*MW.O)/1000)+8*extent_1(11)+2*extent_1(12))/2;     %Extent of reaction 21, mol/h
        end;
        %Calculate the amount of species dissolved (generated) in compartment 1 using the extents of reaction 
        liquid_generation_1 = zeros(1,7);                                                                                           %Initialise the vector containing the liquid generation terms
        liquid_generation_1(1) = (4*extent_1(4)+extent_1(5)+27*extent_1(7)+27*extent_1(9)+27*extent_1(11))*(MW.Cu/1000);            %Rate of Cu2+ production in autoclave compartment 1, kg/h
        liquid_generation_1(2) = (extent_1(1)+extent_1(2)+6*extent_1(3)+3*extent_1(8)+3*extent_1(10)+3*extent_1(12))*(MW.Ni/1000);  %Rate of Ni2+ production in autoclave compartment 1, kg/h
        liquid_generation_1(3) = (MW.Fe/1000)*(extent_1(6));                                                                        %Rate of Fe2+/Fe3+ production in autoclave compartment 1, kg/h
        liquid_generation_1(4) = (MW.Rh/1000)*(-8*extent_1(7)-2*extent_1(8)+2*extent_1(13)+4*extent_1(14)+2*extent_1(15));          %Rate of Rh3+ production in autoclave compartment 1, kg/h
        liquid_generation_1(5) = (MW.Ru/1000)*(-8*extent_1(9)-2*extent_1(10)+4*extent_1(16)+4*extent_1(17)+2*extent_1(18));         %Rate of Ru3+ production in autoclave compartment 1, kg/h
        liquid_generation_1(6) = (MW.Ir/1000)*(-8*extent_1(11)-2*extent_1(12)+2*extent_1(19)+4*extent_1(20)+2*extent_1(21));        %Rate of Ir3+ production in autoclave compartment 1, kg/h
        liquid_generation_1(7) = ((2*MW.H+MW.S+4*MW.O)/1000)*((-2*extent_1(1)+4*extent_1(3)-8*extent_1(4)-extent_1(6)+8*extent_1(8)+8*extent_1(10)+8*extent_1(12)-12*extent_1(14)-6*extent_1(15)*4*extent_1(16)-12*extent_1(17)-6*extent_1(18)-12*extent_1(20)-6*extent_1(21))/2);%Rate of H2SO4 production in autoclave compartment 1, kg/h
        %If dissolved species are consumed rather than generated, check that the amount consumed does not exceed the available amount (only possible for Rh3+, Ru3+, and Ir3+, assuming that excess acid is available)
        for i=(1:7)
            if (liquid_generation_1(i)<(-1*X(41)*recycle_tank_constant*X(25)^0.5*X(i+41)))             %Dissolved species is the limiting reagent rather than the solid species
                liquid_generation_1(i) = (-1*X(41)*recycle_tank_constant*X(25)^0.5*X(i+41));           %Set dissolved species consumption equal to amount of dissolved species fed
                %Recalculate the extents of reaction for the reactions where the dissolved species are the limiting reagents 
                if (i==4)
                    extent_1(7) = (liquid_generation_1(i)/(MW.Rh/1000)-2*extent_1(13)-4*extent_1(14)-2*extent_1(15))*r_1(7)/(-8*r_1(7)-2*r_1(8));
                    extent_1(8) = (liquid_generation_1(i)/(MW.Rh/1000)-2*extent_1(13)-4*extent_1(14)-2*extent_1(15))*r_1(8)/(-8*r_1(7)-2*r_1(8));
                end;
                if (i==5)
                    extent_1(9) = (liquid_generation_1(i)/(MW.Ru/1000)-4*extent_1(16)-4*extent_1(17)-2*extent_1(18))*r_1(9)/(-8*r_1(9)-2*r_1(10));
                    extent_1(10) = (liquid_generation_1(i)/(MW.Ru/1000)-4*extent_1(16)-4*extent_1(17)-2*extent_1(18))*r_1(10)/(-8*r_1(9)-2*r_1(10));
                end;                
                if (i==6)
                    extent_1(11) = (liquid_generation_1(i)/(MW.Ir/1000)-2*extent_1(19)-4*extent_1(20)-2*extent_1(21))*r_1(11)/(-8*r_1(11)-2*r_1(12));
                    extent_1(12) = (liquid_generation_1(i)/(MW.Ir/1000)-2*extent_1(19)-4*extent_1(20)-2*extent_1(21))*r_1(12)/(-8*r_1(11)-2*r_1(12));
                end;
            end;
        end;
        %Recalculate the extents of reaction for the reactions dependent on the extents of reaction for reactions 7 to 12
        extent_1(3) = (solid_consumption_1(2)/((3*MW.Ni+4*MW.S)/1000)+extent_1(1)-extent_1(8)-extent_1(10)-extent_1(12))/2; %Extent of reaction 3, mol/h
        extent_1(4) = solid_consumption_1(3)/((9*MW.Cu+5*MW.S)/1000)-3*extent_1(7)-3*extent_1(9)-3*extent_1(11);            %Extent of reaction 4, mol/h
        extent_1(5) = solid_consumption_1(4)/((MW.Cu+MW.S)/1000)+5*extent_1(4);                                             %Extent of reaction 5, mol/h
        extent_1(15) = (solid_consumption_1(8)/((MW.Rh+2*MW.O)/1000)+8*extent_1(7)+2*extent_1(8))/2;                        %Extent of reaction 15, mol/h
        extent_1(18) = (solid_consumption_1(11)/((MW.Ru+2*MW.O)/1000)+8*extent_1(9)+2*extent_1(10))/2;                      %Extent of reaction 18, mol/h
        extent_1(21) = (solid_consumption_1(14)/((MW.Ir+2*MW.O)/1000)+8*extent_1(11)+2*extent_1(12))/2;                     %Extent of reaction 21, mol/h
        %Check that recalculated extents of reaction can be achieved given the reaction rates; if not, set extents equal to that determined by reaction rate and recalculate the solid consumption
        if (extent_1(3)>r_1(3)*volume_1*60)
            extent_1(3) = r_1(3)*volume_1*60;
            solid_consumption_1(2) = (2*extent_1(3)-extent_1(1)+extent_1(8)+extent_1(10)+extent_1(12))*((3*MW.Ni+4*MW.S)/1000);
        end;
        if (extent_1(4)>r_1(4)*volume_1*60)
            extent_1(4) = r_1(4)*volume_1*60;
            solid_consumption_1(3) = (extent_1(4)+3*extent_1(7)+3*extent_1(9)+3*extent_1(11))*((9*MW.Cu+5*MW.S)/1000);
        end;
        if (extent_1(5)>r_1(5)*volume_1*60)
            extent_1(5) = r_1(5)*volume_1*60;
            solid_consumption_1(4) = (extent_1(5)-5*extent_1(4))*((MW.Cu+MW.S)/1000);
        end;
        if (extent_1(15)>r_1(15)*volume_1*60)
            extent_1(15) = r_1(15)*volume_1*60;
            solid_consumption_1(8) = (2*extent_1(15)-8*extent_1(7)-2*extent_1(8))*((MW.Rh+2*MW.O)/1000);
        end;
        if (extent_1(18)>r_1(18)*volume_1*60)
            extent_1(18) = r_1(18)*volume_1*60;
            solid_consumption_1(11) = (2*extent_1(18)-8*extent_1(9)-2*extent_1(10))*((MW.Ru+2*MW.O)/1000);
        end;
        if (extent_1(21)>r_1(21)*volume_1*60)
            extent_1(21) = r_1(21)*volume_1*60;
            solid_consumption_1(14) = (2*extent_1(21)-8*extent_1(11)-2*extent_1(12))*((MW.Ir+2*MW.O)/1000);
        end;
        %Calculate the rate of oxygen consumption and water condensation in autoclave compartment 1, kg/h
        O2_consumption_1 = (0.5*extent_1(1)+2*extent_1(2)+15*extent_1(3)+2*extent_1(4)+2*extent_1(5)+38*extent_1(7)+8*extent_1(8)+38*extent_1(9)+8*extent_1(10)+38*extent_1(11)+8*extent_1(12)+6*extent_1(13)+3*extent_1(14)-0.5*extent_1(15)+15*extent_1(16)+3*extent_1(17)-0.5*extent_1(18)+6*extent_1(19)+3*extent_1(20)-0.5*extent_1(21))*(2*MW.O/1000); %Rate of oxygen consumption in compartment 1, kg/h 
        n_O2_consumption_1 = O2_consumption_1/(2*MW.O);                                 %Molar rate at which oxygen is consumed in compartment 1, kmol/h
        n_H2O_condensation_1 = (n_O2_consumption_1/y_O2_2nd_stage)*(1-y_O2_2nd_stage);  %Molar rate at which water is transferred from vapour phase to liquid phase in compartment 1, kmol/h
        H2O_condensation_1 = n_H2O_condensation_1*(2*MW.H+MW.O);                        %Mass rate at which water is transferred from vapour phase to liquid phase in compartment 1, kg/h
        %Perform overall mass balance
        m_AC1 = recycle_tank_constant*X(25)^0.5 - m_9 + O2_consumption_1 + H2O_condensation_1;  
        %Calculate the energy consumed/generated by the respective chemical reactions, kJ/h (<0 for exothermic, >0 for endothermic)
        Q_Reaction_1 = H_Reaction.*extent_1; 
        %Calculate the energy transferred to the slurry due to condensation of water from the vapour phase, kJ/h
        Q_Cond_1 = H2O_condensation_1*(H_Evap_H2O+Cp_H2O*(100-25)+Cp_steam.a*(X(73)-100));

%Compartment 2 reaction rates and extents of reactions

        %Estimate volumetric flow rate of the feed to the second compartment of the autoclave, stream AC1
        V_S_AC1 = X(50)*m_AC1/density_residue;                  %Volumetric flow rate of solids in stream AC1, l/h
        V_L_AC1 = X(65)*m_AC1/density_spent;                    %Volumetric flow rate of liquid in stream AC1, l/h
        V_AC1 = V_S_AC1 + V_L_AC1;                              %Total volumetric flow rate of stream AC1, l/h
        %Calculate concentrations of the reagents in autoclave compartment 2
        %Dissolved species
        liq_fractions_2.H2SO4 = X(96);  %Concentration of aqueous species in autoclave compartment 2, kg H2SO4/kg liquid
        liq_fractions_2.Rh = X(93);     %Concentration of aqueous species in autoclave compartment 2, kg Rh/kg liquid
        liq_fractions_2.Ru = X(94);     %Concentration of aqueous species in autoclave compartment 2, kg Ru/kg liquid
        liq_fractions_2.Ir = X(95);     %Concentration of aqueous species in autoclave compartment 2, kg Ir/kg liquid
        %Solid species entering
        solids_in_2.NiS = X(50)*m_AC1*X(51)/V_AC1;        %Concentration of solids entering the second autoclave compartment, kg NiS/l slurry
        solids_in_2.Ni3S4 = X(50)*m_AC1*X(52)/V_AC1;      %Concentration of solids entering the second autoclave compartment, kg Ni3S4/l slurry
        solids_in_2.Cu9S5 = X(50)*m_AC1*X(53)/V_AC1;      %Concentration of solids entering the second autoclave compartment, kg Cu9S5/l slurry
        solids_in_2.FeOHSO4 = X(50)*m_AC1*X(55)/V_AC1;    %Concentration of solids entering the second autoclave compartment, kg FeOHSO4/l slurry
        solids_in_2.Rh2S3 = X(50)*m_AC1*X(56)/V_AC1;      %Concentration of solids entering the second autoclave compartment, kg Rh2S3/l slurry
        solids_in_2.Rh = X(50)*m_AC1*X(57)/V_AC1;         %Concentration of solids entering the second autoclave compartment, kg Rh/l slurry
        solids_in_2.RuS2 = X(50)*m_AC1*X(59)/V_AC1;       %Concentration of solids entering the second autoclave compartment, kg RuS2/l slurry
        solids_in_2.Ru = X(50)*m_AC1*X(60)/V_AC1;         %Concentration of solids entering the second autoclave compartment, kg Ru/l slurry
        solids_in_2.Ir2S3 = X(50)*m_AC1*X(62)/V_AC1;      %Concentration of solids entering the second autoclave compartment, kg Ir2S3/l slurry
        solids_in_2.Ir = X(50)*m_AC1*X(63)/V_AC1;         %Concentration of solids entering the second autoclave compartment, kg Ir/l slurry
        %Estimate the mass (kg/h) and volumetric (l/h) flow rates of stream AC2
        m_AC2_guess = m_AC1;
        error_comp2 = 1;
        count_comp2 = 1;
        while (error_comp2 == 1)
        error_comp2 = 0;
        V_S_AC2 = X(74)*m_AC2_guess/density_residue;
        V_L_AC2 = X(89)*m_AC2_guess/density_spent;
        V_AC2 = V_S_AC2+V_L_AC2;            
        %Solid species present/leaving
        solids_out_2.NiS = m_AC2_guess*X(74)*X(75)/V_AC2;       %Concentration of solids leaving the second autoclave compartment, kg NiS/l slurry
        solids_out_2.Ni3S4 = m_AC2_guess*X(74)*X(76)/V_AC2;     %Concentration of solids leaving the secondt autoclave compartment, kg Ni3S4/l slurry
        solids_out_2.Cu9S5 = m_AC2_guess*X(74)*X(77)/V_AC2;     %Concentration of solids leaving the second autoclave compartment, kg Cu9S5/l slurry
        solids_out_2.CuS = m_AC2_guess*X(74)*X(78)/V_AC2;       %Concentration of solids leaving the second autoclave compartment, kg CuS/l slurry
        solids_out_2.FeOHSO4 = m_AC2_guess*X(74)*X(79)/V_AC2;   %Concentration of solids leaving the second autoclave compartment, kg FeOHSO4/l slurry
        solids_out_2.Rh2S3 = m_AC2_guess*X(74)*X(80)/V_AC2;     %Concentration of solids leaving the second autoclave compartment, kg Rh2S3/l slurry
        solids_out_2.Rh = m_AC2_guess*X(74)*X(81)/V_AC2;        %Concentration of solids leaving the second autoclave compartment, kg Rh/l slurry
        solids_out_2.RuS2 = m_AC2_guess*X(74)*X(83)/V_AC2;      %Concentration of solids leaving the second autoclave compartment, kg RuS2/l slurry
        solids_out_2.Ru = m_AC2_guess*X(74)*X(84)/V_AC2;        %Concentration of solids leaving the second autoclave compartment, kg Ru/l slurry
        solids_out_2.Ir2S3 = m_AC2_guess*X(74)*X(86)/V_AC2;     %Concentration of solids leaving the second autoclave compartment, kg Ir2S3/l slurry
        solids_out_2.Ir = m_AC2_guess*X(74)*X(87)/V_AC2;        %Concentration of solids leaving the second autoclave compartment, kg Ir/l slurry
        solids_out_2.RhO2 = m_AC2_guess*X(74)*X(82)/V_AC2;      %Concentration of solids leaving the second autoclave compartment, kg RhO2/l slurry
        solids_out_2.RuO2 = m_AC2_guess*X(74)*X(85)/V_AC2;      %Concentration of solids leaving the second autoclave compartment, kg RuO2/l slurry
        solids_out_2.IrO2 = m_AC2_guess*X(74)*X(88)/V_AC2;      %Concentration of solids leaving the second autoclave compartment, kg IrO2/l slurry
        %Oxygen solubility (mol/kg) in stream AC2
        Cu_molal_AC2 = 1000*(X(90)/(1-X(90)*(MW.Cu+MW.S+4*MW.O)/MW.Cu-X(91)*(MW.Ni+MW.S+4*MW.O)/MW.Ni-X(96)))/MW.Cu;                    %Molal concentration of copper in stream AC2, mol Cu per kg solvent
        Ni_molal_AC2 = 1000*(X(91)/(1-X(90)*(MW.Cu+MW.S+4*MW.O)/MW.Cu-X(91)*(MW.Ni+MW.S+4*MW.O)/MW.Ni-X(96)))/MW.Ni;                    %Molal concentration of nickel in stream AC2, mol Ni per kg solvent
        H2SO4_molal_AC2 = 1000*(X(96)/(1-X(90)*(MW.Cu+MW.S+4*MW.O)/MW.Cu-X(91)*(MW.Ni+MW.S+4*MW.O)/MW.Ni-X(96)))/(2*MW.H+MW.S+4*MW.O);  %Molal concentration of H2SO4 in stream AC2, mol H2SO4 per kg solvent
        O2_solubility_AC2 = (2*MW.O/1000)*(calculate_oxygen_solubility(X(97),(100000*autoclave_P),Cu_molal_AC2,Ni_molal_AC2,H2SO4_molal_AC2));                                              %Oxygen solubility in stream AC2, kg O2 per kg water
        molal_O2_solubility_AC2 = (O2_solubility_AC2*1000)/(2*MW.O);
        %Calculate the reaction rates in compartment 2, mol/l.min
        r_2 = calculate_reaction_rates(X(97),liq_fractions_2,solids_in_2,solids_out_2,molal_O2_solubility_AC2,Cu9S5_fresh,data);%Calculate the reaction rates of the respective reactions, mol/l.min
        %Calculate the amount of solid species consumed in compartment 2, assuming that all reagents are present in excess (no species is completely converted)
        solid_consumption_2 = zeros(1,14);%Initialise the vector containing the rate of solid species consumption, kg/h
        solid_consumption_2(1) = volume_2*60*((MW.Ni+MW.S)/1000)*(4*r_2(1)+r_2(2));                             %Rate of NiS consumption in autoclave compartment 2, kg/h
        solid_consumption_2(2) = volume_2*60*((3*MW.Ni+4*MW.S)/1000)*(-r_2(1)+2*r_2(3)+r_2(8)+r_2(10)+r_2(12)); %Rate of Ni3S4 consumption in autoclave compartment 2, kg/h
        solid_consumption_2(3) = volume_2*60*((9*MW.Cu+5*MW.S)/1000)*(r_2(4)+3*r_2(7)+3*r_2(9)+3*r_2(11));      %Rate of Cu9S5 consumption in autoclave compartment 2, kg/h
        solid_consumption_2(4) = volume_2*60*((MW.Cu+MW.S)/1000)*(-5*r_2(4)+r_2(5));                            %Rate of CuS consumption in autoclave compartment 2, kg/h
        solid_consumption_2(5) = volume_2*60*((MW.Fe+5*MW.O+MW.H+MW.S)/1000)*(r_2(6));                          %Rate of FeOHSO4 consumption in autoclave compartment 2, kg/h
        solid_consumption_2(6) = volume_2*60*((2*MW.Rh+3*MW.S)/1000)*(r_2(13));                                 %Rate of Rh2S3 consumption in autoclave compartment 2, kg/h
        solid_consumption_2(7) = volume_2*60*((MW.Rh)/1000)*(4*r_2(14));                                        %Rate of Rh consumption in autoclave compartment 2, kg/h
        solid_consumption_2(8) = volume_2*60*((MW.Rh+2*MW.O)/1000)*(-8*r_2(7)-2*r_2(8)+2*r_2(15));              %Rate of RhO2 consumption in autoclave compartment 2, kg/h
        solid_consumption_2(9) = volume_2*60*((MW.Ru+2*MW.S)/1000)*(4*r_2(16));                                 %Rate of RuS2 consumption in autoclave compartment 2, kg/h
        solid_consumption_2(10) = volume_2*60*((MW.Ru)/1000)*(4*r_2(17));                                       %Rate of Ru consumption in autoclave compartment 2, kg/h
        solid_consumption_2(11) = volume_2*60*((MW.Ru+2*MW.O)/1000)*(-8*r_2(9)-2*r_2(10)+2*r_2(18));            %Rate of RuO2 consumption in autoclave compartment 2, kg/h
        solid_consumption_2(12) = volume_2*60*((2*MW.Ir+3*MW.S)/1000)*(r_2(19));                                %Rate of Ir2S3 consumption in autoclave compartment 2, kg/h
        solid_consumption_2(13) = volume_2*60*((MW.Ir)/1000)*(4*r_2(20));                                       %Rate of Ir consumption in autoclave compartment 2, kg/h
        solid_consumption_2(14) = volume_2*60*((MW.Ir+2*MW.O)/1000)*(-8*r_2(11)-2*r_2(12)+2*r_2(21));           %Rate of IrO2 consumption in autoclave compartment 2, kg/h
        %If solids are produced rather than consumed, check that amount produced does not exceed stochiometrically possible amount
        for i=(1:14)
            if (solid_consumption_2(i)<0)
                if (i==1 && solid_consumption_2(i)<0)
                    solid_consumption_2(i)=0; 
                elseif (i==2 && solid_consumption_2(i)<(-0.25*(X(51)*X(50)*m_AC1)/(MW.Ni+MW.S)*(3*MW.Ni+4*MW.S)))
                    solid_consumption_2(i) = (-0.25*(X(51)*X(50)*m_AC1)/(MW.Ni+MW.S)*(3*MW.Ni+4*MW.S));
                elseif (i==3 && solid_consumption_2(i)<0)
                    solid_consumption_2(i)=0;                
                elseif (i==4 && solid_consumption_2(i)<(-5*X(53)*X(50)*m_AC1/(9*MW.Cu+5*MW.S)*(MW.Cu+MW.S)))
                    solid_consumption_2(i) = (-5*X(53)*X(50)*m_AC1/(9*MW.Cu+5*MW.S)*(MW.Cu+MW.S));  
                elseif (i==5 && solid_consumption_2(i)<0)
                    solid_consumption_2(i) = 0;        
                elseif (i==6 && solid_consumption_2(i)<0)
                    solid_consumption_2(i) = 0;        
                elseif (i==7 && solid_consumption_2(i)<0)
                    solid_consumption_2(i) = 0;                
                elseif (i==8 && solid_consumption_2(i)<(-1*X(69)*X(65)*m_AC1/(MW.Rh)*(MW.Rh+2*MW.O)))
                    solid_consumption_2(i) = (-1*X(69)*X(65)*m_AC1/(MW.Rh)*(MW.Rh+2*MW.O));                    
                elseif (i==9 && solid_consumption_2(i)<0)
                    solid_consumption_2(i) = 0;                   
                elseif (i==10 && solid_consumption_2(i)<0)
                    solid_consumption_2(i) = 0;
                elseif (i==11 && solid_consumption_2(i)<(-1*X(70)*X(65)*m_AC1/(MW.Ru)*(MW.Ru+2*MW.O)))
                    solid_consumption_2(i) = (-1*X(70)*X(65)*m_AC1/(MW.Ru)*(MW.Ru+2*MW.O));
                elseif (i==12 && solid_consumption_2(i)<0)
                    solid_consumption_2(i) = 0;
                elseif (i==13 && solid_consumption_2(i)<0)
                    solid_consumption_2(i) = 0;
                elseif (i==14 && solid_consumption_2(i)<(-1*X(71)*X(65)*m_AC1/(MW.Ir)*(MW.Ir+2*MW.O)))
                    solid_consumption_2(i) = (-1*X(71)*X(65)*
    

4 Comments

Niel on 18 Jul 2013

Hi. Thanks for the reply.

I specifically checked for this - and no... The variable is assigned another value normally. When I comment out the variable in question, another variable elsewhere in the code gives the same error.

That is why I believe that the problem lies somewhere with how Simulink executes the function block?

Kaustubha Govind on 18 Jul 2013

It's unlikely this has anything to do with how Simulink is executing, because these are compile-time errors (not run-time errors) being thrown by the block. You don't have to hit "Play" to see them, just run Ctrl+D to compile the model without running it. Also, when I do a text search, I don't see any instance of 'H2O_condensation_2' in your code.

Niel on 18 Jul 2013

I see - thank you. It seems that not all the code was copied successfully (only about half). Below is the rest of it (it starts with last section of the previous part of the code):

          %If solids are produced rather than consumed, check that amount produced does not exceed stochiometrically possible amount
          for i=(1:14)
              if (solid_consumption_2(i)<0)
                  if (i==1 && solid_consumption_2(i)<0)
                      solid_consumption_2(i)=0; 
                  elseif (i==2 && solid_consumption_2(i)<(-0.25*(X(51)*X(50)*m_AC1)/(MW.Ni+MW.S)*(3*MW.Ni+4*MW.S)))
                      solid_consumption_2(i) = (-0.25*(X(51)*X(50)*m_AC1)/(MW.Ni+MW.S)*(3*MW.Ni+4*MW.S));
                  elseif (i==3 && solid_consumption_2(i)<0)
                      solid_consumption_2(i)=0;                
                  elseif (i==4 && solid_consumption_2(i)<(-5*X(53)*X(50)*m_AC1/(9*MW.Cu+5*MW.S)*(MW.Cu+MW.S)))
                      solid_consumption_2(i) = (-5*X(53)*X(50)*m_AC1/(9*MW.Cu+5*MW.S)*(MW.Cu+MW.S));  
                  elseif (i==5 && solid_consumption_2(i)<0)
                      solid_consumption_2(i) = 0;        
                  elseif (i==6 && solid_consumption_2(i)<0)
                      solid_consumption_2(i) = 0;        
                  elseif (i==7 && solid_consumption_2(i)<0)
                      solid_consumption_2(i) = 0;                
                  elseif (i==8 && solid_consumption_2(i)<(-1*X(69)*X(65)*m_AC1/(MW.Rh)*(MW.Rh+2*MW.O)))
                      solid_consumption_2(i) = (-1*X(69)*X(65)*m_AC1/(MW.Rh)*(MW.Rh+2*MW.O));                    
                  elseif (i==9 && solid_consumption_2(i)<0)
                      solid_consumption_2(i) = 0;                   
                  elseif (i==10 && solid_consumption_2(i)<0)
                      solid_consumption_2(i) = 0;
                  elseif (i==11 && solid_consumption_2(i)<(-1*X(70)*X(65)*m_AC1/(MW.Ru)*(MW.Ru+2*MW.O)))
                      solid_consumption_2(i) = (-1*X(70)*X(65)*m_AC1/(MW.Ru)*(MW.Ru+2*MW.O));
                  elseif (i==12 && solid_consumption_2(i)<0)
                      solid_consumption_2(i) = 0;
                  elseif (i==13 && solid_consumption_2(i)<0)
                      solid_consumption_2(i) = 0;
                  elseif (i==14 && solid_consumption_2(i)<(-1*X(71)*X(65)*m_AC1*/(MW.Ir)*(MW.Ir+2*MW.O)))
                      solid_consumption_2(i) = (-1*X(71)*X(65)*m_AC1/(MW.Ir)*(MW.Ir+2*MW.O));
                  end;
              end;
          end;  
          %Check that solid consumption does not exceed the amount of solids fed
          for i=(1:14)
              %If complete conversion(consumption) of solid species, set the consumption equal to the amount fed
              if (solid_consumption_2(i)>X(50)*m_AC1*X(i+50))    
                  solid_consumption_2(i)=X(50)*m_AC1*X(i+50);          
              end;
              %Limit the amount of CuS leached(consumed) to leachable fraction of the fresh CuS fed (empirical restriction)
              if (i==4 && (solid_consumption_1(i)+solid_consumption_2(i))>(CuS_leachable*X(2)*X(6)*prep_tank_constant*(X(1)^0.5)))
                  solid_consumption_2(i) = (CuS_leachable*X(2)*X(6)*prep_tank_constant*(X(1)^0.5))-solid_consumption_1(i);
              end;
          end;
          %Calculate the extent of reaction for each reaction (assuming that the solid species are limiting reactants)
  %         extent_2 = zeros(1,21);                                                                     %Initialise the vector containing the extents of reaction
          if (r_2(1) == 0)
              extent_2(1) = 0;                                                                        %If reaction rate is zero, extent of reaction must be zero
          else
              extent_2(1) = solid_consumption_2(1)/((MW.Ni+MW.S)/1000)*(r_2(1))/(4*r_2(1)+r_2(2));    %Extent of reaction 1, mol/h
          end;
          if (r_2(2) == 0)
              extent_2(2) = 0;
          else
              extent_2(2) = solid_consumption_2(1)/((MW.Ni+MW.S)/1000)*(r_2(2))/(4*r_2(1)+r_2(2));    %Extent of reaction 2, mol/h
          end;
          if (r_2(4) == 0)
              extent_2(4) = 0;
          else
              extent_2(4) = solid_consumption_2(3)/((9*MW.Cu+5*MW.S)/1000)*(r_2(4))/(r_2(4)+3*r_2(7)+3*r_2(9)+3*r_2(11)); %Extent of reaction 4, mol/h
          end;
          if (r_2(7) == 0)
              extent_2(7) = 0;
          else
              extent_2(7) = solid_consumption_2(3)/((9*MW.Cu+5*MW.S)/1000)*(r_2(7))/(r_2(4)+3*r_2(7)+3*r_2(9)+3*r_2(11)); %Extent of reaction 7, mol/h
          end;
          if (r_2(9) == 0)
              extent_2(9) = 0;
          else
              extent_2(9) = solid_consumption_2(3)/((9*MW.Cu+5*MW.S)/1000)*(r_2(9))/(r_2(4)+3*r_2(7)+3*r_2(9)+3*r_2(11)); %Extent of reaction 9, mol/h
          end;
          if (r_2(11) == 0)
              extent_2(11) = 0;
          else
              extent_2(11) = solid_consumption_2(3)/((9*MW.Cu+5*MW.S)/1000)*(r_2(11))/(r_2(4)+3*r_2(7)+3*r_2(9)+3*r_2(11));%Extent of reaction 11, mol/h
          end;
          if (r_2(13) == 0)
              extent_2(13) = 0;
          else
              extent_2(13) = solid_consumption_2(6)/((2*MW.Rh+3*MW.S)/1000);                          %Extent of reaction 13, mol/h
          end;
          if (r_2(14) == 0)
              extent_2(14) = 0;
          else
              extent_2(14) = solid_consumption_2(7)/((MW.Rh)/1000)/4;                                 %Extent of reaction 14, mol/h
          end;
          if (r_2(16) == 0)
              extent_2(16) = 0;
          else
              extent_2(16) = solid_consumption_2(9)/((MW.Ru+2*MW.S)/1000)/4;                          %Extent of reaction 16, mol/h
          end;
          if (r_2(17) == 0)
              extent_2(17) = 0;
          else
              extent_2(17) = solid_consumption_2(10)/((MW.Ru)/1000)/4;                                %Extent of reaction 17, mol/h
          end;
          if (r_2(19) == 0)
              extent_2(19) = 0;
          else
              extent_2(19) = solid_consumption_2(12)/((2*MW.Ir+3*MW.S)/1000);                         %Extent of reaction 19, mol/h
          end;
          if (r_2(20) == 0)
              extent_2(20) = 0;
          else
              extent_2(20) = solid_consumption_2(13)/((MW.Ir)/1000)/4;                                %Extent of reaction 20, mol/h
          end;
          if (r_2(6) == 0)
              extent_2(6) = 0;
          else
              extent_2(6) = solid_consumption_2(5)/((MW.Fe+5*MW.O+MW.H+MW.S)/1000);                   %Extent of reaction 6, mol/h
          end;
          if (r_2(5) == 0)
              extent_2(5) = 0;
          else
              extent_2(5) = solid_consumption_2(4)/((MW.Cu+MW.S)/1000)+5*extent_2(4);                 %Extent of reaction 5, mol/h
          end;
          if (r_2(3) == 0)
              extent_2(3) = 0;
          else
              extent_2(3) = (solid_consumption_2(2)/((3*MW.Ni+4*MW.S)/1000)+extent_2(1))*r_2(3)/(2*r_2(3)+r_2(8)+r_2(10)+r_2(12));    %Extent of reaction 3, mol/h
          end;
          if (r_2(8) == 0)
              extent_2(8) = 0;
          else
              extent_2(8) = (solid_consumption_2(2)/((3*MW.Ni+4*MW.S)/1000)+extent_2(1))*r_2(8)/(2*r_2(3)+r_2(8)+r_2(10)+r_2(12));    %Extent of reaction 8, mol/h
          end;
          if (r_2(10) == 0)
              extent_2(10) = 0;
          else
              extent_2(10) = (solid_consumption_2(2)/((3*MW.Ni+4*MW.S)/1000)+extent_2(1))*r_2(10)/(2*r_2(3)+r_2(8)+r_2(10)+r_2(12));  %Extent of reaction 10, mol/h
          end;
          if (r_2(12) == 0)
              extent_2(12) = 0;
          else
              extent_2(12) = (solid_consumption_2(2)/((3*MW.Ni+4*MW.S)/1000)+extent_2(1))*r_2(12)/(2*r_2(3)+r_2(8)+r_2(10)+r_2(12));  %Extent of reaction 12, mol/h
          end;
          if (r_2(15) == 0)
              extent_2(15) = 0;
          else
              extent_2(15) = (solid_consumption_2(8)/((MW.Rh+2*MW.O)/1000)+8*extent_2(7)+2*extent_2(8))/2;        %Extent of reaction 15, mol/h
          end;
          if (r_2(18) == 0)
              extent_2(18) = 0;
          else
              extent_2(18) = (solid_consumption_2(11)/((MW.Ru+2*MW.O)/1000)+8*extent_2(9)+2*extent_2(10))/2;      %Extent of reaction 18, mol/h
          end;
          if (r_2(21) == 0)
              extent_2(21) = 0;
          else
              extent_2(21) = (solid_consumption_2(14)/((MW.Ir+2*MW.O)/1000)+8*extent_2(11)+2*extent_2(12))/2;     %Extent of reaction 21, mol/h
          end;
          %Calculate the amount of species dissolved (generated) using the extents of reaction 
          liquid_generation_2 = zeros(1,7);                                                                                           %Initialise the vector containing the liquid generation terms
          liquid_generation_2(1) = (4*extent_2(4)+extent_2(5)+27*extent_2(7)+27*extent_2(9)+27*extent_2(11))*(MW.Cu/1000);            %Rate of Cu2+ production in autoclave compartment 2, kg/h
          liquid_generation_2(2) = (extent_2(1)+extent_2(2)+6*extent_2(3)+3*extent_2(8)+3*extent_2(10)+3*extent_2(12))*(MW.Ni/1000);  %Rate of Ni2+ production in autoclave compartment 2, kg/h
          liquid_generation_2(3) = (MW.Fe/1000)*(extent_2(6));                                                                        %Rate of Fe2+/Fe3+ production in autoclave compartment 2, kg/h
          liquid_generation_2(4) = (MW.Rh/1000)*(-8*extent_2(7)-2*extent_2(8)+2*extent_2(13)+4*extent_2(14)+2*extent_2(15));          %Rate of Rh3+ production in autoclave compartment 2, kg/h
          liquid_generation_2(5) = (MW.Ru/1000)*(-8*extent_2(9)-2*extent_2(10)+4*extent_2(16)+4*extent_2(17)+2*extent_2(18));         %Rate of Ru3+ production in autoclave compartment 2, kg/h
          liquid_generation_2(6) = (MW.Ir/1000)*(-8*extent_2(11)-2*extent_2(12)+2*extent_2(19)+4*extent_2(20)+2*extent_2(21));        %Rate of Ir3+ production in autoclave compartment 2, kg/h
          liquid_generation_2(7) = ((2*MW.H+MW.S+4*MW.O)/1000)*((-2*extent_2(1)+4*extent_2(3)-8*extent_2(4)-extent_2(6)+8*extent_2(8)+8*extent_2(10)+8*extent_2(12)-12*extent_2(14)-6*extent_2(15)*4*extent_2(16)-12*extent_2(17)-6*extent_2(18)-12*extent_2(20)-6*extent_2(21))/2);%Rate of H2SO4 production in autoclave compartment 2, kg/h
          %If dissolved species are consumed rather than generated, check that the amount consumed does not exceed the available amount
          for i=(1:7)
              if (liquid_generation_2(i)<(-1*X(65)*m_AC1*X(i+65)))           %Dissolved species is the limiting reagent rather than the solid species
                  liquid_generation_2(i) = -1*X(65)*m_AC1*X(i+65);           %Set dissolved species consumption equal to amount of dissolved species fed
                  %Recalculate the extents of reaction for the reactions where the dissolved species are the limiting reagents
                  if (i==4)
                      extent_2(7) = (liquid_generation_2(i)/(MW.Rh/1000)-2*extent_2(13)-4*extent_2(14)-2*extent_2(15))*r_2(7)/(-8*r_2(7)-2*r_2(8));
                      extent_2(8) = (liquid_generation_2(i)/(MW.Rh/1000)-2*extent_2(13)-4*extent_2(14)-2*extent_2(15))*r_2(8)/(-8*r_2(7)-2*r_2(8));
                  end;
                  if (i==5)
                      extent_2(9) = (liquid_generation_2(i)/(MW.Ru/1000)-4*extent_2(16)-4*extent_2(17)-2*extent_2(18))*r_2(9)/(-8*r_2(9)-2*r_2(10));
                      extent_2(10) = (liquid_generation_2(i)/(MW.Ru/1000)-4*extent_2(16)-4*extent_2(17)-2*extent_2(18))*r_2(10)/(-8*r_2(9)-2*r_2(10));
                  end;                
                  if (i==6)
                      extent_2(11) = (liquid_generation_2(i)/(MW.Ir/1000)-2*extent_2(19)-4*extent_2(20)-2*extent_2(21))*r_2(11)/(-8*r_2(11)-2*r_2(12));
                      extent_2(12) = (liquid_generation_2(i)/(MW.Ir/1000)-2*extent_2(19)-4*extent_2(20)-2*extent_2(21))*r_2(12)/(-8*r_2(11)-2*r_2(12));
                  end;
              end;
          end;
          %Recalculate the extents of reaction for the reactions dependent on the extents of reaction for reactions 7 to 12
          extent_2(3) = (solid_consumption_2(2)/((3*MW.Ni+4*MW.S)/1000)+extent_2(1)-extent_2(8)-extent_2(10)-extent_2(12))/2; %Extent of reaction 3, mol/h
          extent_2(4) = solid_consumption_2(3)/((9*MW.Cu+5*MW.S)/1000)-3*extent_2(7)-3*extent_2(9)-3*extent_2(11);            %Extent of reaction 4, mol/h
          extent_2(5) = solid_consumption_2(4)/((MW.Cu+MW.S)/1000)+5*extent_2(4);                                             %Extent of reaction 5, mol/h
          extent_2(15) = (solid_consumption_2(8)/((MW.Rh+2*MW.O)/1000)+8*extent_2(7)+2*extent_2(8))/2;                        %Extent of reaction 15, mol/h
          extent_2(18) = (solid_consumption_2(11)/((MW.Ru+2*MW.O)/1000)+8*extent_2(9)+2*extent_2(10))/2;                      %Extent of reaction 18, mol/h
          extent_2(21) = (solid_consumption_2(14)/((MW.Ir+2*MW.O)/1000)+8*extent_2(11)+2*extent_2(12))/2;                     %Extent of reaction 21, mol/h
          %Check that recalculated extents of reaction can be achieved given the reaction rates; if not, set extents equal to that determined by reaction rate
          if (extent_2(3)>r_2(3)*volume_2*60)
              extent_2(3) = r_2(3)*volume_2*60;
              solid_consumption_2(2) = (2*extent_2(3)-extent_2(1)+extent_2(8)+extent_2(10)+extent_2(12))*((3*MW.Ni+4*MW.S)/1000);
          end;
          if (extent_2(4)>r_2(4)*volume_2*60)
              extent_2(4) = r_2(4)*volume_2*60;
              solid_consumption_2(3) = (extent_2(4)+3*extent_2(7)+3*extent_2(9)+3*extent_2(11))*((9*MW.Cu+5*MW.S)/1000);
          end;
          if (extent_2(5)>r_2(5)*volume_2*60)
              extent_2(5) = r_2(5)*volume_2*60;
              solid_consumption_2(4) = (extent_2(5)-5*extent_2(4))*((MW.Cu+MW.S)/1000);
          end;
          if (extent_2(15)>r_2(15)*volume_2*60)
              extent_2(15) = r_2(15)*volume_2*60;
              solid_consumption_2(8) = (2*extent_2(15)-8*extent_2(7)-2*extent_2(8))*((MW.Rh+2*MW.O)/1000);
          end;
          if (extent_2(18)>r_2(18)*volume_2*60)
              extent_2(18) = r_2(18)*volume_2*60;
              solid_consumption_2(11) = (2*extent_2(18)-8*extent_2(9)-2*extent_2(10))*((MW.Ru+2*MW.O)/1000);
          end;
          if (extent_2(21)>r_2(21)*volume_2*60)
              extent_2(21) = r_2(21)*volume_2*60;
              solid_consumption_2(14) = (2*extent_2(21)-8*extent_2(11)-2*extent_2(12))*((MW.Ir+2*MW.O)/1000);
          end;
          %Calculate the rate of oxygen consumption in autoclave compartment 2, kg/h
          O2_consumption_2 = (0.5*extent_2(1)+2*extent_2(2)+15*extent_2(3)+2*extent_2(4)+2*extent_2(5)+38*extent_2(7)+8*extent_2(8)+38*extent_2(9)+8*extent_2(10)+38*extent_2(11)+8*extent_2(12)+6*extent_2(13)+3*extent_2(14)-0.5*extent_2(15)+15*extent_2(16)+3*extent_2(17)-0.5*extent_2(18)+6*extent_2(19)+3*extent_2(20)-0.5*extent_2(21))*(2*MW.O/1000);
          O2_consumption_AC4 = O2_consumption_2-m_10;                                         %Rate at which oxygen in vapour space is consumed (if <0, more oxygen in AC3 than AC4)
          n_O2_consumption_AC4 = O2_consumption_AC4/(2*MW.O);                                 %Molar rate at which oxygen is consumed from stream AC4, kmol/h
          n_H2O_condensation_2 = (n_O2_consumption_AC4/y_O2_2nd_stage)*(1-y_O2_2nd_stage);    %Molar rate at which water is transferred from vapour phase to liquid phase in compartment 1, kmol/h (if negative, water evaporates)
          H2O_condensation_2 = n_H2O_condensation_2*(2*MW.H+MW.O);                            %Mass rate at which water is transferred from vapour phase to liquid phase in compartment 1, kg/h
          %Perform overall mass balance
          m_AC2 = m_AC1+O2_consumption_2+H2O_condensation_2;  
          variable_change_comp2 = abs(m_AC2_guess-m_AC2);
          %Increment loop counter
          count_comp2 = count_comp2+1;
          if ((variable_change_comp2>tolerance_percent)&&(count_comp2<max_loops))
              error_comp2 = 1;
          end;
          %Recalculate the initial guesses based on the calculated values
          m_AC2_guess = m_AC2_guess+0.2*(m_AC2-m_AC2_guess);
          end;
          %Calculate the energy consumed/generated by the respective chemical reactions, kJ/h (<0 for exothermic, >0 for endothermic)
          Q_Reaction_2 = H_Reaction.*extent_2; 
          Q_Cond_2 = H2O_condensation_2*(H_Evap_H2O+Cp_H2O*(100-25)+Cp_steam.a*(X(97)-100));
          %Calculate properties of the oxygen fed to autoclave compartment 2 (stream 10) 
          Cp_10 = Cp_O2.a;        %Oxygen heat capacity, kJ/kg.dC
          T_10 = oxygen_T;        %Oxygen temperature, dC
%Compartment 3 reaction rates and extents of reactions
          %Calculate the reagent concentrations in autoclave compartment 3 in order to estimate the reaction rates
          %Dissolved species
          liq_fractions_3.H2SO4 = X(120);          %Concentration of aqueous species in autoclave compartment 3, kg H2SO4/kg liquid
          liq_fractions_3.Rh = X(117);             %Concentration of aqueous species in autoclave compartment 3, kg Rh/kg liquid
          liq_fractions_3.Ru = X(118);             %Concentration of aqueous species in autoclave compartment 3, kg Ru/kg liquid
          liq_fractions_3.Ir = X(119);             %Concentration of aqueous species in autoclave compartment 3, kg Ir/kg liquid
          %Solid species entering
          solids_in_3.NiS = X(74)*m_AC2*X(75)/V_AC2;          %Concentration of solids entering the third autoclave compartment, kg NiS/l slurry
          solids_in_3.Ni3S4 = X(74)*m_AC2*X(76)/V_AC2;        %Concentration of solids entering the third autoclave compartment, kg Ni3S4/l slurry
          solids_in_3.Cu9S5 = X(74)*m_AC2*X(77)/V_AC2;        %Concentration of solids entering the third autoclave compartment, kg Cu9S5/l slurry
          solids_in_3.FeOHSO4 = X(74)*m_AC2*X(79)/V_AC2;      %Concentration of solids entering the third autoclave compartment, kg FeOHSO4/l slurry
          solids_in_3.Rh2S3 = X(74)*m_AC2*X(80)/V_AC2;        %Concentration of solids entering the third autoclave compartment, kg Rh2S3/l slurry
          solids_in_3.Rh = X(74)*m_AC2*X(81)/V_AC2;           %Concentration of solids entering the third autoclave compartment, kg Rh/l slurry
          solids_in_3.RuS2 = X(74)*m_AC2*X(83)/V_AC2;         %Concentration of solids entering the third autoclave compartment, kg RuS2/l slurry
          solids_in_3.Ru = X(74)*m_AC2*X(84)/V_AC2;           %Concentration of solids entering the third autoclave compartment, kg Ru/l slurry
          solids_in_3.Ir2S3 = X(74)*m_AC2*X(86)/V_AC2;        %Concentration of solids entering the third autoclave compartment, kg Ir2S3/l slurry
          solids_in_3.Ir = X(74)*m_AC2*X(87)/V_AC2;           %Concentration of solids entering the third autoclave compartment, kg Ir/l slurry
          m_14_guess = m_AC2;
          error_comp3 = 1;
          count_comp3 = 1;
          while (error_comp3 == 1)
          error_comp3 = 0;
          %Estimate the volumetric flow rate of stream 14
          V_S_14 = m_14_guess*X(98)/density_residue;          %Volumetric flow rate of solids in stream 14, l/h
          V_L_14 = m_14_guess*X(113)/density_spent;           %Volumetric flow rate of liquid in stream 14, l/h
          V_14 = V_S_14 + V_L_14;                             %Total volumetric flow rate of stream 14, l/h
          %Solids in compartment 3, leaving the compartment
          solids_out_3.NiS = m_14_guess*X(98)*X(99)/V_14;        %Concentration of solids leaving the third autoclave compartment, kg NiS/l slurry
          solids_out_3.Ni3S4 = m_14_guess*X(98)*X(100)/V_14;     %Concentration of solids leaving the third autoclave compartment, kg Ni3S4/l slurry
          solids_out_3.Cu9S5 = m_14_guess*X(98)*X(101)/V_14;     %Concentration of solids leaving the third autoclave compartment, kg Cu9S5/l slurry
          solids_out_3.CuS = m_14_guess*X(98)*X(102)/V_14;       %Concentration of solids leaving the third autoclave compartment, kg CuS/l slurry
          solids_out_3.FeOHSO4 = m_14_guess*X(98)*X(103)/V_14;   %Concentration of solids leaving the third autoclave compartment, kg FeOHSO4/l slurry
          solids_out_3.Rh2S3 = m_14_guess*X(98)*X(104)/V_14;     %Concentration of solids leaving the third autoclave compartment, kg Rh2S3/l slurry
          solids_out_3.Rh = m_14_guess*X(98)*X(105)/V_14;        %Concentration of solids leaving the third autoclave compartment, kg Rh/l slurry
          solids_out_3.RuS2 = m_14_guess*X(98)*X(107)/V_14;      %Concentration of solids leaving the third autoclave compartment, kg RuS2/l slurry
          solids_out_3.Ru = m_14_guess*X(98)*X(108)/V_14;        %Concentration of solids leaving the third autoclave compartment, kg Ru/l slurry
          solids_out_3.Ir2S3 = m_14_guess*X(98)*X(110)/V_14;     %Concentration of solids leaving the third autoclave compartment, kg Ir2S3/l slurry
          solids_out_3.Ir = m_14_guess*X(98)*X(111)/V_14;        %Concentration of solids leaving the third autoclave compartment, kg Ir/l slurry
          solids_out_3.RhO2 = m_14_guess*X(98)*X(106)/V_14;      %Concentration of solids leaving the third autoclave compartment, kg RhO2/l slurry
          solids_out_3.RuO2 = m_14_guess*X(98)*X(109)/V_14;      %Concentration of solids leaving the third autoclave compartment, kg RuO2/l slurry
          solids_out_3.IrO2 = m_14_guess*X(98)*X(112)/V_14;      %Concentration of solids leaving the third autoclave compartment, kg IrO2/l slurry
          %Oxygen solubility (mol/kg) in stream 14
          Cu_molal_14 = 1000*(X(114)/(1-X(114)*(MW.Cu+MW.S+4*MW.O)/MW.Cu-X(115)*(MW.Ni+MW.S+4*MW.O)/MW.Ni-X(120)))/MW.Cu; %Molal concentration of copper in stream 14, mol Cu per kg solvent
          Ni_molal_14 = 1000*(X(115)/(1-X(114)*(MW.Cu+MW.S+4*MW.O)/MW.Cu-X(115)*(MW.Ni+MW.S+4*MW.O)/MW.Ni-X(120)))/MW.Ni; %Molal concentration of nickel in stream 14, mol Ni per kg solvent
          H2SO4_molal_14 = 1000*(X(120)/(1-X(114)*(MW.Cu+MW.S+4*MW.O)/MW.Cu-X(115)*(MW.Ni+MW.S+4*MW.O)/MW.Ni-X(120)))/(2*MW.H+MW.S+4*MW.O); %Molal concentration of H2SO4 in stream 14, mol H2SO4 per kg solvent
          O2_solubility_14 = (2*MW.O/1000)*(calculate_oxygen_solubility(X(121),(100000*autoclave_P),Cu_molal_14,Ni_molal_14,H2SO4_molal_14)); %Oxygen solubility in stream 14, kg O2 per kg water
          molal_O2_solubility_14 = (O2_solubility_14*1000)/(2*MW.O);
          %Calculate the reaction rates for the respective reactions in autoclave compartment 3
          r_3 = calculate_reaction_rates(X(121),liq_fractions_3,solids_in_3,solids_out_3,molal_O2_solubility_14,Cu9S5_fresh,data);%Calculate the reaction rates of the respective reactions, mol/l.min
          %Calculate the amount of solid species consumed, assuming that all reagents are present in excess (no species is completely converted)
          solid_consumption_3 = zeros(1,14);                                                                      %Initialise the vector containing the rate of solid species consumption, kg/h
          solid_consumption_3(1) = volume_3*60*((MW.Ni+MW.S)/1000)*(4*r_3(1)+r_3(2));                             %Rate of NiS consumption in autoclave compartment 3, kg/h
          solid_consumption_3(2) = volume_3*60*((3*MW.Ni+4*MW.S)/1000)*(-r_3(1)+2*r_3(3)+r_3(8)+r_3(10)+r_3(12)); %Rate of Ni3S4 consumption in autoclave compartment 3, kg/h
          solid_consumption_3(3) = volume_3*60*((9*MW.Cu+5*MW.S)/1000)*(r_3(4)+3*r_3(7)+3*r_3(9)+3*r_3(11));      %Rate of Cu9S5 consumption in autoclave compartment 3, kg/h
          solid_consumption_3(4) = volume_3*60*((MW.Cu+MW.S)/1000)*(-5*r_3(4)+r_3(5));                            %Rate of CuS consumption in autoclave compartment 3, kg/h
          solid_consumption_3(5) = volume_3*60*((MW.Fe+5*MW.O+MW.H+MW.S)/1000)*(r_3(6));                          %Rate of FeOHSO4 consumption in autoclave compartment 3, kg/h
          solid_consumption_3(6) = volume_3*60*((2*MW.Rh+3*MW.S)/1000)*(r_3(13));                                 %Rate of Rh2S3 consumption in autoclave compartment 3, kg/h
          solid_consumption_3(7) = volume_3*60*((MW.Rh)/1000)*(4*r_3(14));                                        %Rate of Rh consumption in autoclave compartment 3, kg/h
          solid_consumption_3(8) = volume_3*60*((MW.Rh+2*MW.O)/1000)*(-8*r_3(7)-2*r_3(8)+2*r_3(15));              %Rate of RhO2 consumption in autoclave compartment 3, kg/h
          solid_consumption_3(9) = volume_3*60*((MW.Ru+2*MW.S)/1000)*(4*r_3(16));                                 %Rate of RuS2 consumption in autoclave compartment 3, kg/h
          solid_consumption_3(10) = volume_3*60*((MW.Ru)/1000)*(4*r_3(17));                                       %Rate of Ru consumption in autoclave compartment 3, kg/h
          solid_consumption_3(11) = volume_3*60*((MW.Ru+2*MW.O)/1000)*(-8*r_3(9)-2*r_3(10)+2*r_3(18));            %Rate of RuO2 consumption in autoclave compartment 3, kg/h
          solid_consumption_3(12) = volume_3*60*((2*MW.Ir+3*MW.S)/1000)*(r_3(19));                                %Rate of Ir2S3 consumption in autoclave compartment 3, kg/h
          solid_consumption_3(13) = volume_3*60*((MW.Ir)/1000)*(4*r_3(20));                                       %Rate of Ir consumption in autoclave compartment 3, kg/h
          solid_consumption_3(14) = volume_3*60*((MW.Ir+2*MW.O)/1000)*(-8*r_3(11)-2*r_3(12)+2*r_3(21));           %Rate of IrO2 consumption in autoclave compartment 3, kg/h
          %If solids are produced rather than consumed, check that amount produced does not exceed stochiometrically possible amount
          for i=(1:14)
              if (solid_consumption_3(i)<0)
                  if (i==1 && solid_consumption_3(i)<0)
                      solid_consumption_3(i)=0; 
                  elseif (i==2 && solid_consumption_3(i)<(-0.25*(X(75)*X(74)*m_AC2)/(MW.Ni+MW.S)*(3*MW.Ni+4*MW.S)))
                      solid_consumption_3(i) = (-0.25*(X(75)*X(74)*m_AC2)/(MW.Ni+MW.S)*(3*MW.Ni+4*MW.S));
                  elseif (i==3 && solid_consumption_3(i)<0)
                      solid_consumption_3(i)=0;                
                  elseif (i==4 && solid_consumption_3(i)<(-5*X(77)*X(74)*m_AC2/(9*MW.Cu+5*MW.S)*(MW.Cu+MW.S)))
                      solid_consumption_3(i) = (-5*X(77)*X(74)*m_AC2/(9*MW.Cu+5*MW.S)*(MW.Cu+MW.S));  
                  elseif (i==5 && solid_consumption_3(i)<0)
                      solid_consumption_3(i) = 0;        
                  elseif (i==6 && solid_consumption_3(i)<0)
                      solid_consumption_3(i) = 0;        
                  elseif (i==7 && solid_consumption_3(i)<0)
                      solid_consumption_3(i) = 0;                
                  elseif (i==8 && solid_consumption_3(i)<(-1*X(93)*X(89)*m_AC2/(MW.Rh)*(MW.Rh+2*MW.O)))
                      solid_consumption_3(i) = (-1*X(93)*X(89)*m_AC2/(MW.Rh)*(MW.Rh+2*MW.O));                    
                  elseif (i==9 && solid_consumption_3(i)<0)
                      solid_consumption_3(i) = 0;                   
                  elseif (i==10 && solid_consumption_3(i)<0)
                      solid_consumption_3(i) = 0;
                  elseif (i==11 && solid_consumption_3(i)<(-1*X(94)*X(89)*m_AC2/(MW.Ru)*(MW.Ru+2*MW.O)))
                      solid_consumption_3(i) = (-1*X(94)*X(89)*m_AC2/(MW.Ru)*(MW.Ru+2*MW.O));
                  elseif (i==12 && solid_consumption_3(i)<0)
                      solid_consumption_3(i) = 0;
                  elseif (i==13 && solid_consumption_3(i)<0)
                      solid_consumption_3(i) = 0;
                  elseif (i==14 && solid_consumption_3(i)<(-1*X(95)*X(89)*m_AC2/(MW.Ir)*(MW.Ir+2*MW.O)))
                      solid_consumption_3(i) = (-1*X(95)*X(89)*m_AC2/(MW.Ir)*(MW.Ir+2*MW.O));
                  end;
              end;
          end; 
          %Check that solid consumption does not exceed the amount of solids fed
          for i=(1:14)
              %If consumption exceeds the amount fed, set consumption equal to maximum possible amount (complete conversion)
              if (solid_consumption_3(i)>X(74)*m_AC2*X(i+74))
                  solid_consumption_3(i)=X(74)*m_AC2*X(i+74);          
              end;
              %Limit the amount of CuS leached(consumed) to leachable fraction of the fresh CuS fed (empirical restriction)
              if (i==4 && (solid_consumption_1(i)+solid_consumption_2(i)+solid_consumption_3(i))>(CuS_leachable*X(2)*X(6)*prep_tank_constant*X(1)^0.5))
                  solid_consumption_3(i) = (CuS_leachable*X(2)*X(6)*prep_tank_constant*X(1)^0.5)-solid_consumption_1(i)-solid_consumption_2(i);
              end;
          end;
          %Calculate the extent of reaction for each reaction (assuming that the solid species are limiting reactants)
          extent_3 = zeros(1,21);                                                                     %Initialise the vector containing the extents of reaction
          if (r_3(1) == 0)
              extent_3(1) = 0;
          else
              extent_3(1) = solid_consumption_3(1)/((MW.Ni+MW.S)/1000)*(r_3(1))/(4*r_3(1)+r_3(2));    %Extent of reaction 1, mol/h
          end;
          if (r_3(2) == 0)
              extent_3(2) = 0;
          else
              extent_3(2) = solid_consumption_3(1)/((MW.Ni+MW.S)/1000)*(r_3(2))/(4*r_3(1)+r_3(2));    %Extent of reaction 2, mol/h
          end;
          if (r_3(4) == 0)
              extent_3(4) = 0;
          else
              extent_3(4) = solid_consumption_3(3)/((9*MW.Cu+5*MW.S)/1000)*(r_3(4))/(r_3(4)+3*r_3(7)+3*r_3(9)+3*r_3(11)); %Extent of reaction 4, mol/h
          end;
          if (r_3(7) == 0)
              extent_3(7) = 0;
          else
              extent_3(7) = solid_consumption_3(3)/((9*MW.Cu+5*MW.S)/1000)*(r_3(7))/(r_3(4)+3*r_3(7)+3*r_3(9)+3*r_3(11)); %Extent of reaction 7, mol/h
          end;
          if (r_3(9) == 0)
              extent_3(9) = 0;
          else
              extent_3(9) = solid_consumption_3(3)/((9*MW.Cu+5*MW.S)/1000)*(r_3(9))/(r_3(4)+3*r_3(7)+3*r_3(9)+3*r_3(11)); %Extent of reaction 9, mol/h
          end;
          if (r_3(11) == 0)
              extent_3(11) = 0;
          else
              extent_3(11) = solid_consumption_3(3)/((9*MW.Cu+5*MW.S)/1000)*(r_3(11))/(r_3(4)+3*r_3(7)+3*r_3(9)+3*r_3(11));%Extent of reaction 11, mol/h
          end;
          if (r_3(13) == 0)
              extent_3(13) = 0;
          else
              extent_3(13) = solid_consumption_3(6)/((2*MW.Rh+3*MW.S)/1000);                          %Extent of reaction 13, mol/h
          end;
          if (r_3(14) == 0)
              extent_3(14) = 0;
          else
              extent_3(14) = solid_consumption_3(7)/((MW.Rh)/1000)/4;                                 %Extent of reaction 14, mol/h
          end;
          if (r_3(16) == 0)
              extent_3(16) = 0;
          else
              extent_3(16) = solid_consumption_3(9)/((MW.Ru+2*MW.S)/1000)/4;                          %Extent of reaction 16, mol/h
          end;
          if (r_3(17) == 0)
              extent_3(17) = 0;
          else
              extent_3(17) = solid_consumption_3(10)/((MW.Ru)/1000)/4;                                %Extent of reaction 17, mol/h
          end;
          if (r_3(19) == 0)
              extent_3(19) = 0;
          else
              extent_3(19) = solid_consumption_3(12)/((2*MW.Ir+3*MW.S)/1000);                         %Extent of reaction 19, mol/h
          end;
          if (r_3(20) == 0)
              extent_3(20) = 0;
          else
              extent_3(20) = solid_consumption_3(13)/((MW.Ir)/1000)/4;                                %Extent of reaction 20, mol/h
          end;
          if (r_3(6) == 0)
              extent_3(6) = 0;
          else
              extent_3(6) = solid_consumption_3(5)/((MW.Fe+5*MW.O+MW.H+MW.S)/1000);                   %Extent of reaction 6, mol/h
          end;
          if (r_3(5) == 0)
              extent_3(5) = 0;
          else
              extent_3(5) = solid_consumption_3(4)/((MW.Cu+MW.S)/1000)+5*extent_3(4);                 %Extent of reaction 5, mol/h
          end;
          if (r_3(3) == 0)
              extent_3(3) = 0;
          else
              extent_3(3) = (solid_consumption_3(2)/((3*MW.Ni+4*MW.S)/1000)+extent_3(1))*r_3(3)/(2*r_3(3)+r_3(8)+r_3(10)+r_3(12));    %Extent of reaction 3, mol/h
          end;
          if (r_3(8) == 0)
              extent_3(8) = 0;
          else
              extent_3(8) = (solid_consumption_3(2)/((3*MW.Ni+4*MW.S)/1000)+extent_3(1))*r_3(8)/(2*r_3(3)+r_3(8)+r_3(10)+r_3(12));    %Extent of reaction 8, mol/h
          end;
          if (r_3(10) == 0)
              extent_3(10) = 0;
          else
              extent_3(10) = (solid_consumption_3(2)/((3*MW.Ni+4*MW.S)/1000)+extent_3(1))*r_3(10)/(2*r_3(3)+r_3(8)+r_3(10)+r_3(12));  %Extent of reaction 10, mol/h
          end;
          if (r_3(12) == 0)
              extent_3(12) = 0;
          else
              extent_3(12) = (solid_consumption_3(2)/((3*MW.Ni+4*MW.S)/1000)+extent_3(1))*r_3(12)/(2*r_3(3)+r_3(8)+r_3(10)+r_3(12));  %Extent of reaction 12, mol/h
          end;
          if (r_3(15) == 0)
              extent_3(15) = 0;
          else
              extent_3(15) = (solid_consumption_3(8)/((MW.Rh+2*MW.O)/1000)+8*extent_3(7)+2*extent_3(8))/2;        %Extent of reaction 15, mol/h
          end;
          if (r_3(18) == 0)
              extent_3(18) = 0;
          else
              extent_3(18) = (solid_consumption_3(11)/((MW.Ru+2*MW.O)/1000)+8*extent_3(9)+2*extent_3(10))/2;      %Extent of reaction 18, mol/h
          end;
          if (r_3(21) == 0)
              extent_3(21) = 0;
          else
              extent_3(21) = (solid_consumption_3(14)/((MW.Ir+2*MW.O)/1000)+8*extent_3(11)+2*extent_3(12))/2;     %Extent of reaction 21, mol/h
          end;
          %Calculate the amount of species dissolved (generated) using the extents of reaction 
          liquid_generation_3 = zeros(1,7);                                                                                           %Initialise the vector containing the liquid generation terms
          liquid_generation_3(1) = (4*extent_3(4)+extent_3(5)+27*extent_3(7)+27*extent_3(9)+27*extent_3(11))*(MW.Cu/1000);            %Rate of Cu2+ production in autoclave compartment 3, kg/h
          liquid_generation_3(2) = (extent_3(1)+extent_3(2)+6*extent_3(3)+3*extent_3(8)+3*extent_3(10)+3*extent_3(12))*(MW.Ni/1000);  %Rate of Ni2+ production in autoclave compartment 3, kg/h
          liquid_generation_3(3) = (MW.Fe/1000)*(extent_3(6));                                                                        %Rate of Fe2+/Fe3+ production in autoclave compartment 3, kg/h
          liquid_generation_3(4) = (MW.Rh/1000)*(-8*extent_3(7)-2*extent_3(8)+2*extent_3(13)+4*extent_3(14)+2*extent_3(15));          %Rate of Rh3+ production in autoclave compartment 3, kg/h
          liquid_generation_3(5) = (MW.Ru/1000)*(-8*extent_3(9)-2*extent_3(10)+4*extent_3(16)+4*extent_3(17)+2*extent_3(18));         %Rate of Ru3+ production in autoclave compartment 3, kg/h
          liquid_generation_3(6) = (MW.Ir/1000)*(-8*extent_3(11)-2*extent_3(12)+2*extent_3(19)+4*extent_3(20)+2*extent_3(21));        %Rate of Ir3+ production in autoclave compartment 3, kg/h
          liquid_generation_3(7) = ((2*MW.H+MW.S+4*MW.O)/1000)*((-2*extent_3(1)+4*extent_3(3)-8*extent_3(4)-extent_3(6)+8*extent_3(8)+8*extent_3(10)+8*extent_3(12)-12*extent_3(14)-6*extent_3(15)*4*extent_3(16)-12*extent_3(17)-6*extent_3(18)-12*extent_3(20)-6*extent_3(21))/2);%Rate of H2SO4 production in autoclave compartment 3, kg/h
          %If dissolved species are consumed rather than generated, check that the amount consumed does not exceed the available amount
          for i=(1:7)
              if (liquid_generation_3(i)<(-1*X(89)*m_AC2*X(i+89)))           %Dissolved species is the limiting reagent rather than the solid species
                  liquid_generation_3(i) = -1*X(89)*m_AC2*X(i+89);           %Set dissolved species consumption equal to amount of dissolved species fed
                  %Recalculate the extents of reaction for the reactions where the dissolved species are the limiting reagents
                  if (i==4)
                      extent_3(7) = (liquid_generation_3(i)/(MW.Rh/1000)-2*extent_3(13)-4*extent_3(14)-2*extent_3(15))*r_3(7)/(-8*r_3(7)-2*r_3(8));
                      extent_3(8) = (liquid_generation_3(i)/(MW.Rh/1000)-2*extent_3(13)-4*extent_3(14)-2*extent_3(15))*r_3(8)/(-8*r_3(7)-2*r_3(8));
                  end;                
                  if (i==5)
                      extent_3(9) = (liquid_generation_3(i)/(MW.Ru/1000)-4*extent_3(16)-4*extent_3(17)-2*extent_3(18))*r_3(9)/(-8*r_3(9)-2*r_3(10));
                      extent_3(10) = (liquid_generation_3(i)/(MW.Ru/1000)-4*extent_3(16)-4*extent_3(17)-2*extent_3(18))*r_3(10)/(-8*r_3(9)-2*r_3(10));
                  end;                
                  if (i==6)
                      extent_3(11) = (liquid_generation_3(i)/(MW.Ir/1000)-2*extent_3(19)-4*extent_3(20)-2*extent_3(21))*r_3(11)/(-8*r_3(11)-2*r_3(12));
                      extent_3(12) = (liquid_generation_3(i)/(MW.Ir/1000)-2*extent_3(19)-4*extent_3(20)-2*extent_3(21))*r_3(12)/(-8*r_3(11)-2*r_3(12));
                  end;
              end;
          end;
          %Recalculate the extents of reaction dependent on the extents of reactions 7 to 12
          extent_3(3) = (solid_consumption_3(2)/((3*MW.Ni+4*MW.S)/1000)+extent_3(1)-extent_3(8)-extent_3(10)-extent_3(12))/2; %Extent of reaction 3, mol/h
          extent_3(4) = solid_consumption_3(3)/((9*MW.Cu+5*MW.S)/1000)-3*extent_3(7)-3*extent_3(9)-3*extent_3(11);            %Extent of reaction 4, mol/h
          extent_3(5) = solid_consumption_3(4)/((MW.Cu+MW.S)/1000)+5*extent_3(4);                                             %Extent of reaction 5, mol/h
          extent_3(15) = (solid_consumption_3(8)/((MW.Rh+2*MW.O)/1000)+8*extent_3(7)+2*extent_3(8))/2;                        %Extent of reaction 15, mol/h
          extent_3(18) = (solid_consumption_3(11)/((MW.Ru+2*MW.O)/1000)+8*extent_3(9)+2*extent_3(10))/2;                      %Extent of reaction 18, mol/h
          extent_3(21) = (solid_consumption_3(14)/((MW.Ir+2*MW.O)/1000)+8*extent_3(11)+2*extent_3(12))/2;                     %Extent of reaction 21, mol/h
          %Check that recalculated extents of reaction can be achieved given the reaction rates; if not, set extents equal to that determined by reaction rate
          if (extent_3(3)>r_3(3)*volume_3*60)
              extent_3(3) = r_3(3)*volume_3*60;
              solid_consumption_3(2) = (2*extent_3(3)-extent_3(1)+extent_3(8)+extent_3(10)+extent_3(12))*((3*MW.Ni+4*MW.S)/1000);
          end;
          if (extent_3(4)>r_3(4)*volume_3*60)
              extent_3(4) = r_3(4)*volume_3*60;
              solid_consumption_3(3) = (extent_3(4)+3*extent_3(7)+3*extent_3(9)+3*extent_3(11))*((9*MW.Cu+5*MW.S)/1000);
          end;
          if (extent_3(5)>r_3(5)*volume_3*60)
              extent_3(5) = r_3(5)*volume_3*60;
              solid_consumption_3(4) = (extent_3(5)-5*extent_3(4))*((MW.Cu+MW.S)/1000);
          end;
          if (extent_3(15)>r_3(15)*volume_3*60)
              extent_3(15) = r_3(15)*volume_3*60;
              solid_consumption_3(8) = (2*extent_3(15)-8*extent_3(7)-2*extent_3(8))*((MW.Rh+2*MW.O)/1000);
          end;
          if (extent_3(18)>r_3(18)*volume_3*60)
              extent_3(18) = r_3(18)*volume_3*60;
              solid_consumption_3(11) = (2*extent_3(18)-8*extent_3(9)-2*extent_3(10))*((MW.Ru+2*MW.O)/1000);
          end;
          if (extent_3(21)>r_3(21)*volume_3*60)
              extent_3(21) = r_3(21)*volume_3*60;
              solid_consumption_3(14) = (2*extent_3(21)-8*extent_3(11)-2*extent_3(12))*((MW.Ir+2*MW.O)/1000);
          end;
          %Calculate the rate of oxygen consumption in autoclave compartment 3, kg/h
          O2_consumption_3 = (0.5*extent_3(1)+2*extent_3(2)+15*extent_3(3)+2*extent_3(4)+2*extent_3(5)+38*extent_3(7)+8*extent_3(8)+38*extent_3(9)+8*extent_3(10)+38*extent_3(11)+8*extent_3(12)+6*extent_3(13)+3*extent_3(14)-0.5*extent_3(15)+15*extent_3(16)+3*extent_3(17)-0.5*extent_3(18)+6*extent_3(19)+3*extent_3(20)-0.5*extent_3(21))*(2*MW.O/1000); %Rate of oxygen consumption in compartment 3, kg/h 
          O2_consumption_AC5 = O2_consumption_3-m_11;                                     %Rate at which oxygen in vapour space is consumed (if <0, more oxygen in AC3 than AC4)
          n_O2_consumption_AC5 = O2_consumption_AC5/(2*MW.O);                             %Molar rate at which oxygen is consumed from stream AC4, kmol/h
          n_H2O_condensation_3 = (n_O2_consumption_AC5/y_O2_2nd_stage)*(1-y_O2_2nd_stage);%Molar rate at which water is transferred from vapour phase to liquid phase in compartment 1, kmol/h (if negative, water evaporates)
          H2O_condensation_3 = n_H2O_condensation_3*(2*MW.H+MW.O);                        %Mass rate at which water is transferred from vapour phase to liquid phase in compartment 1, kg/h
          %Perform overall mass balance
          m_14 = m_AC2+O2_consumption_3+H2O_condensation_3;  
          %Compare initial guess with the calculated value
          variable_change_comp3 = abs(m_14_guess-m_14);
          %Increment loop counter
          count_comp3 = count_comp3+1;
          %Reset the error flag to one if the convergence criteria has not been met
          if ((variable_change_comp3>tolerance_percent)&&(count_comp3<max_loops))
              error_comp3 = 1;
          end;
          %Recalculate the initial guesses based on the calculated values
          m_14_guess = m_14_guess+0.2*(m_14-m_14_guess);
          end;
          %Calculate the energy consumed/generated by the respective chemical reactions, kJ/h (<0 for exothermic, >0 for endothermic)
          Q_Reaction_3 = H_Reaction.*extent_3; 
          Q_Cond_3 = H2O_condensation_3*(H_Evap_H2O+Cp_H2O*(100-25)+Cp_steam.a*(X(121)-100));
          %Calculate the properties of the oxygen fed to autoclave compartment 3
          Cp_11 = Cp_O2.a;        %Heat capacity of oxygen, kJ/kg.dC
          T_11 = oxygen_T;        %Temperature of oxygen, dC
%Compartment 4 reaction rates and extents of reactions
          m_22_guess = prep_tank_3_constant*(X(170)^0.5);
          %Estimate volumetric flow rate of the feed to the fourth compartment of the autoclave, stream 21
          V_S_21 = X(171)*prep_tank_3_constant*(X(170)^0.5)/density_residue;              %Volumetric flow rate of solids in stream 21, l/h
          V_L_21 = X(186)*prep_tank_3_constant*(X(170)^0.5)/density_spent;                %Volumetric flow rate of liquid in stream 21, l/h
          V_21 = V_S_21 + V_L_21;                       %Total volumetric flow rate of stream 21, l/h
          %Calculate the concentration of reagens present in autoclave compartment 4
          %Dissolved species
          liq_fractions_4.H2SO4 = X(216);          %Concentration of aqueous species in autoclave compartment 4, kg H2SO4/kg liquid
          liq_fractions_4.Rh = X(213);             %Concentration of aqueous species in autoclave compartment 4, kg Rh/kg liquid
          liq_fractions_4.Ru = X(214);             %Concentration of aqueous species in autoclave compartment 4, kg Ru/kg liquid
          liq_fractions_4.Ir = X(215);             %Concentration of aqueous species in autoclave compartment 4, kg Ir/kg liquid
          %Solid species entering
          solids_in_4.NiS = X(171)*prep_tank_3_constant*(X(170)^0.5).*X(172)/V_21;        %Concentration of solids entering the fourth autoclave compartment, kg NiS/l slurry
          solids_in_4.Ni3S4 = X(171)*prep_tank_3_constant*(X(170)^0.5).*X(173)/V_21;      %Concentration of solids entering the fourth autoclave compartment, kg Ni3S4/l slurry
          solids_in_4.Cu9S5 = X(171)*prep_tank_3_constant*(X(170)^0.5).*X(174)/V_21;      %Concentration of solids entering the fourth autoclave compartment, kg Cu9S5/l slurry
          solids_in_4.FeOHSO4 = X(171)*prep_tank_3_constant*(X(170)^0.5).*X(176)/V_21;    %Concentration of solids entering the fourth autoclave compartment, kg FeOHSO4/l slurry
          solids_in_4.Rh2S3 = X(171)*prep_tank_3_constant*(X(170)^0.5).*X(177)/V_21;      %Concentration of solids entering the fourth autoclave compartment, kg Rh2S3/l slurry
          solids_in_4.Rh = X(171)*prep_tank_3_constant*(X(170)^0.5).*X(178)/V_21;         %Concentration of solids entering the fourth autoclave compartment, kg Rh/l slurry
          solids_in_4.RuS2 = X(171)*prep_tank_3_constant*(X(170)^0.5).*X(180)/V_21;       %Concentration of solids entering the fourth autoclave compartment, kg RuS2/l slurry
          solids_in_4.Ru = X(171)*prep_tank_3_constant*(X(170)^0.5).*X(181)/V_21;         %Concentration of solids entering the fourth autoclave compartment, kg Ru/l slurry
          solids_in_4.Ir2S3 = X(171)*prep_tank_3_constant*(X(170)^0.5).*X(183)/V_21;      %Concentration of solids entering the fourth autoclave compartment, kg Ir2S3/l slurry
          solids_in_4.Ir = X(171)*prep_tank_3_constant*(X(170)^0.5).*X(184)/V_21;         %Concentration of solids entering the fourth autoclave compartment, kg Ir/l slurry
          %Calculate the properties of oxygen fed to compartment 4
          Cp_12 = Cp_O2.a;                %Heat capacity of oxygen, kJ/kg.dC
          T_12 = oxygen_T;                %Temperature of oxygen, dC
          %Specify properties of the steam used to heat the content of compartment 4
          T_13 = steam_T;                 %Temperature of stream 13, dC
          error_comp4 = 1;
          count_comp4 = 1;
          while (error_comp4 == 1)
          error_comp4 = 0;
          %Estimate the volumetric flow rate of stream 22
          V_S_22 = m_22_guess*X(194)/density_residue;        %Volumetric flow rate of solids in stream 22, l/h
          V_L_22 = m_22_guess*X(209)/density_spent;          %Volumetric flow rate of liquid in stream 22, l/h
          V_22 = V_S_22 + V_L_22;                            %Total volumetric flow rate of stream 22, l/h
          %Solid species present/leaving
          solids_out_4.NiS = m_22_guess*X(194).*X(195)/V_22;       %Concentration of solids leaving the fourth autoclave compartment, kg NiS/l slurry
          solids_out_4.Ni3S4 = m_22_guess*X(194).*X(196)/V_22;     %Concentration of solids leaving the fourth autoclave compartment, kg Ni3S4/l slurry
          solids_out_4.Cu9S5 = m_22_guess*X(194).*X(197)/V_22;     %Concentration of solids leaving the fourth autoclave compartment, kg Cu9S5/l slurry
          solids_out_4.CuS = m_22_guess*X(194).*X(198)/V_22;       %Concentration of solids leaving the fourth autoclave compartment, kg CuS/l slurry
          solids_out_4.FeOHSO4 = m_22_guess*X(194).*X(199)/V_22;   %Concentration of solids leaving the fourth autoclave compartment, kg FeOHSO4/l slurry
          solids_out_4.Rh2S3 = m_22_guess*X(194).*X(200)/V_22;     %Concentration of solids leaving the fourth autoclave compartment, kg Rh2S3/l slurry
          solids_out_4.Rh = m_22_guess*X(194).*X(201)/V_22;        %Concentration of solids leaving the fourth autoclave compartment, kg Rh/l slurry
          solids_out_4.RuS2 = m_22_guess*X(194).*X(203)/V_22;      %Concentration of solids leaving the fourth autoclave compartment, kg RuS2/l slurry
          solids_out_4.Ru = m_22_guess*X(194).*X(204)/V_22;        %Concentration of solids leaving the fourth autoclave compartment, kg Ru/l slurry
          solids_out_4.Ir2S3 = m_22_guess*X(194).*X(206)/V_22;     %Concentration of solids leaving the fourth autoclave compartment, kg Ir2S3/l slurry
          solids_out_4.Ir = m_22_guess*X(194).*X(207)/V_22;        %Concentration of solids leaving the fourth autoclave compartment, kg Ir/l slurry
          solids_out_4.RhO2 = m_22_guess*X(194).*X(202)/V_22;      %Concentration of solids leaving the fourth autoclave compartment, kg RhO2/l slurry
          solids_out_4.RuO2 = m_22_guess*X(194).*X(205)/V_22;      %Concentration of solids leaving the fourth autoclave compartment, kg RuO2/l slurry
          solids_out_4.IrO2 = m_22_guess*X(194).*X(208)/V_22;      %Concentration of solids leaving the fourth autoclave compartment, kg IrO2/l slurry
          %Oxygen solubility
          Cu_molal_22 = 1000*(X(210)/(1-X(210)*(MW.Cu+MW.S+4*MW.O)/MW.Cu-X(211)*(MW.Ni+MW.S+4*MW.O)/MW.Ni-X(216)))/MW.Cu;                     %Molal concentration of copper in stream 22, mol Cu per kg solvent
          Ni_molal_22 = 1000*(X(211)/(1-X(210)*(MW.Cu+MW.S+4*MW.O)/MW.Cu-X(211)*(MW.Ni+MW.S+4*MW.O)/MW.Ni-X(216)))/MW.Ni;                     %Molal concentration of nickel in stream 22, mol Ni per kg solvent
          H2SO4_molal_22 = 1000*(X(216)/(1-X(210)*(MW.Cu+MW.S+4*MW.O)/MW.Cu-X(211)*(MW.Ni+MW.S+4*MW.O)/MW.Ni-X(216)))/(2*MW.H+MW.S+4*MW.O);   %Molal concentration of H2SO4 in stream 22, mol H2SO4 per kg solvent
          O2_solubility_22 = (2*MW.O/1000)*(calculate_oxygen_solubility(X(217),(100000*autoclave_P),Cu_molal_22,Ni_molal_22,H2SO4_molal_22));                                         %Oxygen solubility in stream 22, kg O2 per kg water
          molal_O2_solubility_22 = (O2_solubility_22*1000)/(2*MW.O);
          %Calculate the reaction rates in autoclave compartment 4, mol/l.min
          r_4 = calculate_reaction_rates(X(217),liq_fractions_4,solids_in_4,solids_out_4,molal_O2_solubility_22,Cu9S5_fresh,data);
          %Calculate the amount of solid species consumed, assuming that all reagents are present in excess (no species is completely converted)
          solid_consumption_4 = zeros(1,14);                                                                      %Initialise the vector containing the rate of solid species consumption, kg/h
          solid_consumption_4(1) = volume_4*60*((MW.Ni+MW.S)/1000)*(4*r_4(1)+r_4(2));                             %Rate of NiS consumption in autoclave compartment 4, kg/h
          solid_consumption_4(2) = volume_4*60*((3*MW.Ni+4*MW.S)/1000)*(-r_4(1)+2*r_4(3)+r_4(8)+r_4(10)+r_4(12)); %Rate of Ni3S4 consumption in autoclave compartment 4, kg/h
          solid_consumption_4(3) = volume_4*60*((9*MW.Cu+5*MW.S)/1000)*(r_4(4)+3*r_4(7)+3*r_4(9)+3*r_4(11));      %Rate of Cu9S5 consumption in autoclave compartment 4, kg/h
          solid_consumption_4(4) = volume_4*60*((MW.Cu+MW.S)/1000)*(-5*r_4(4)+r_4(5));                            %Rate of CuS consumption in autoclave compartment 4, kg/h
          solid_consumption_4(5) = volume_4*60*((MW.Fe+5*MW.O+MW.H+MW.S)/1000)*(r_4(6));                          %Rate of FeOHSO4 consumption in autoclave compartment 4, kg/h
          solid_consumption_4(6) = volume_4*60*((2*MW.Rh+3*MW.S)/1000)*(r_4(13));                                 %Rate of Rh2S3 consumption in autoclave compartment 4, kg/h
          solid_consumption_4(7) = volume_4*60*((MW.Rh)/1000)*(4*r_4(14));                                        %Rate of Rh consumption in autoclave compartment 4, kg/h
          solid_consumption_4(8) = volume_4*60*((MW.Rh+2*MW.O)/1000)*(-8*r_4(7)-2*r_4(8)+2*r_4(15));              %Rate of RhO2 consumption in autoclave compartment 4, kg/h
          solid_consumption_4(9) = volume_4*60*((MW.Ru+2*MW.S)/1000)*(4*r_4(16));                                 %Rate of RuS2 consumption in autoclave compartment 4, kg/h
          solid_consumption_4(10) = volume_4*60*((MW.Ru)/1000)*(4*r_4(17));                                       %Rate of Ru consumption in autoclave compartment 4, kg/h
          solid_consumption_4(11) = volume_4*60*((MW.Ru+2*MW.O)/1000)*(-8*r_4(9)-2*r_4(10)+2*r_4(18));            %Rate of RuO2 consumption in autoclave compartment 4, kg/h
          solid_consumption_4(12) = volume_4*60*((2*MW.Ir+3*MW.S)/1000)*(r_4(19));                                %Rate of Ir2S3 consumption in autoclave compartment 4, kg/h
          solid_consumption_4(13) = volume_4*60*((MW.Ir)/1000)*(4*r_4(20));                                       %Rate of Ir consumption in autoclave compartment 4, kg/h
          solid_consumption_4(14) = volume_4*60*((MW.Ir+2*MW.O)/1000)*(-8*r_4(11)-2*r_4(12)+2*r_4(21));           %Rate of IrO2 consumption in autoclave compartment 4, kg/h
          %If solids are produced rather than consumed, check that amount produced does not exceed stochiometrically possible amount
          for i=(1:14)
              if (solid_consumption_4(i)<0)
                  if (i==1 && solid_consumption_4(i)<0)
                      solid_consumption_4(i)=0; 
                  elseif (i==2 && solid_consumption_4(i)<(-0.25*(X(172)*X(171)*prep_tank_3_constant*(X(170)^0.5))/(MW.Ni+MW.S)*(3*MW.Ni+4*MW.S)))
                      solid_consumption_4(i) = (-0.25*(X(172)*X(171)*X(171)*prep_tank_3_constant*(X(170)^0.5))/(MW.Ni+MW.S)*(3*MW.Ni+4*MW.S));
                  elseif (i==3 && solid_consumption_4(i)<0)
                      solid_consumption_4(i)=0;                
                  elseif (i==4 && solid_consumption_4(i)<(-5*X(174)*X(171)*prep_tank_3_constant*(X(170)^0.5)/(9*MW.Cu+5*MW.S)*(MW.Cu+MW.S)))
                      solid_consumption_4(i) = (-5*X(174)*X(171)*prep_tank_3_constant*(X(170)^0.5)/(9*MW.Cu+5*MW.S)*(MW.Cu+MW.S));  
                  elseif (i==5 && solid_consumption_4(i)<0)
                      solid_consumption_4(i) = 0;        
                  elseif (i==6 && solid_consumption_4(i)<0)
                      solid_consumption_4(i) = 0;        
                  elseif (i==7 && solid_consumption_4(i)<0)
                      solid_consumption_4(i) = 0;                
                  elseif (i==8 && solid_consumption_4(i)<(-1*X(190)*X(186)*prep_tank_3_constant*(X(170)^0.5)/(MW.Rh)*(MW.Rh+2*MW.O)))
                      solid_consumption_4(i) = (-1*X(190)*X(186)*prep_tank_3_constant*(X(170)^0.5)/(MW.Rh)*(MW.Rh+2*MW.O));                    
                  elseif (i==9 && solid_consumption_4(i)<0)
                      solid_consumption_4(i) = 0;                   
                  elseif (i==10 && solid_consumption_4(i)<0)
                      solid_consumption_4(i) = 0;
                  elseif (i==11 && solid_consumption_4(i)<(-1*X(191)*X(186)*prep_tank_3_constant*(X(170)^0.5)/(MW.Ru)*(MW.Ru+2*MW.O)))
                      solid_consumption_4(i) = (-1*X(191)*X(186)*prep_tank_3_constant*(X(170)^0.5)/(MW.Ru)*(MW.Ru+2*MW.O));
                  elseif (i==12 && solid_consumption_4(i)<0)
                      solid_consumption_4(i) = 0;
                  elseif (i==13 && solid_consumption_4(i)<0)
                      solid_consumption_4(i) = 0;
                  elseif (i==14 && solid_consumption_4(i)<(-1*X(192)*X(186)*prep_tank_3_constant*(X(170)^0.5)/(MW.Ir)*(MW.Ir+2*MW.O)))
                      solid_consumption_4(i) = (-1*X(192)*X(186)*prep_tank_3_constant*(X(170)^0.5)/(MW.Ir)*(MW.Ir+2*MW.O));
                  end;
              end;
          end; 
          %Check that solid consumption does not exceed the amount of solids fed
          for i=(1:14)
               %Complete conversion(consumption) of solid species
              if (solid_consumption_4(i)>X(171)*prep_tank_3_constant*(X(170)^0.5)*X(i+171))   
                  solid_consumption_4(i)=X(171)*prep_tank_3_constant*(X(170)^0.5)*X(i+171);          
              end;
              %Limit the amount of CuS leached(consumed) to leachable fraction of the fresh CuS fed (empirical restriction)
              if (i==4 && (solid_consumption_1(i)+solid_consumption_2(i)+solid_consumption_3(i)+solid_consumption_4(i))>(CuS_leachable*X(2)*X(6)*prep_tank_constant*(X(1)^0.5)))
                  solid_consumption_4(i) = (CuS_leachable*X(2)*X(6)*prep_tank_constant*(X(1)^0.5))-solid_consumption_1(i)-solid_consumption_2(i)-solid_consumption_3(i);
              end;
          end;
          %Calculate the extent of reaction for each reaction (assuming that the solid species are limiting reactants)
          extent_4 = zeros(1,21);                                                                     %Initialise the vector containing the extents of reaction
          if (r_4(1) == 0)
              extent_4(1) = 0;
          else
              extent_4(1) = solid_consumption_4(1)/((MW.Ni+MW.S)/1000)*(r_4(1))/(4*r_4(1)+r_4(2));    %Extent of reaction 1, mol/h
          end;
          if (r_4(2) == 0)
              extent_4(2) = 0;
          else
              extent_4(2) = solid_consumption_4(1)/((MW.Ni+MW.S)/1000)*(r_4(2))/(4*r_4(1)+r_4(2));    %Extent of reaction 2, mol/h
          end;
          if (r_4(4) == 0)
              extent_4(4) = 0;
          else
              extent_4(4) = solid_consumption_4(3)/((9*MW.Cu+5*MW.S)/1000)*(r_4(4))/(r_4(4)+3*r_4(7)+3*r_4(9)+3*r_4(11)); %Extent of reaction 4, mol/h
          end;
          if (r_4(7) == 0)
              extent_4(7) = 0;
          else
              extent_4(7) = solid_consumption_4(3)/((9*MW.Cu+5*MW.S)/1000)*(r_4(7))/(r_4(4)+3*r_4(7)+3*r_4(9)+3*r_4(11)); %Extent of reaction 7, mol/h
          end;
          if (r_4(9) == 0)
              extent_4(9) = 0;
          else
              extent_4(9) = solid_consumption_4(3)/((9*MW.Cu+5*MW.S)/1000)*(r_4(9))/(r_4(4)+3*r_4(7)+3*r_4(9)+3*r_4(11)); %Extent of reaction 9, mol/h
          end;
          if (r_4(11) == 0)
              extent_4(11) = 0;
          else
              extent_4(11) = solid_consumption_4(3)/((9*MW.Cu+5*MW.S)/1000)*(r_4(11))/(r_4(4)+3*r_4(7)+3*r_4(9)+3*r_4(11));%Extent of reaction 11, mol/h
          end;
          if (r_4(13) == 0)
              extent_4(13) = 0;
          else
              extent_4(13) = solid_consumption_4(6)/((2*MW.Rh+3*MW.S)/1000);                          %Extent of reaction 13, mol/h
          end;
          if (r_4(14) == 0)
              extent_4(14) = 0;
          else
              extent_4(14) = solid_consumption_4(7)/((MW.Rh)/1000)/4;                                 %Extent of reaction 14, mol/h
          end;
          if (r_4(16) == 0)
              extent_4(16) = 0;
          else
              extent_4(16) = solid_consumption_4(9)/((MW.Ru+2*MW.S)/1000)/4;                          %Extent of reaction 16, mol/h
          end;
          if (r_4(17) == 0)
              extent_4(17) = 0;
          else
              extent_4(17) = solid_consumption_4(10)/((MW.Ru)/1000)/4;                                %Extent of reaction 17, mol/h
          end;
          if (r_4(19) == 0)
              extent_4(19) = 0;
          else
              extent_4(19) = solid_consumption_4(12)/((2*MW.Ir+3*MW.S)/1000);                         %Extent of reaction 19, mol/h
          end;
          if (r_4(20) == 0)
              extent_4(20) = 0;
          else
              extent_4(20) = solid_consumption_4(13)/((MW.Ir)/1000)/4;                                %Extent of reaction 20, mol/h
          end;
          if (r_4(6) == 0)
              extent_4(6) = 0;
          else
              extent_4(6) = solid_consumption_4(5)/((MW.Fe+5*MW.O+MW.H+MW.S)/1000);                   %Extent of reaction 6, mol/h
          end;
          if (r_4(5) == 0)
              extent_4(5) = 0;
          else
              extent_4(5) = solid_consumption_4(4)/((MW.Cu+MW.S)/1000)+5*extent_4(4);                 %Extent of reaction 5, mol/h
          end;
          if (r_4(3) == 0)
              extent_4(3) = 0;
          else
              extent_4(3) = (solid_consumption_4(2)/((3*MW.Ni+4*MW.S)/1000)+extent_4(1))*r_4(3)/(2*r_4(3)+r_4(8)+r_4(10)+r_4(12));    %Extent of reaction 3, mol/h
          end;
          if (r_4(8) == 0)
              extent_4(8) = 0;
          else
              extent_4(8) = (solid_consumption_4(2)/((3*MW.Ni+4*MW.S)/1000)+extent_4(1))*r_4(8)/(2*r_4(3)+r_4(8)+r_4(10)+r_4(12));    %Extent of reaction 8, mol/h
          end;
          if (r_4(10) == 0)
              extent_4(10) = 0;
          else
              extent_4(10) = (solid_consumption_4(2)/((3*MW.Ni+4*MW.S)/1000)+extent_4(1))*r_4(10)/(2*r_4(3)+r_4(8)+r_4(10)+r_4(12));  %Extent of reaction 10, mol/h
          end;
          if (r_4(12) == 0)
              extent_4(12) = 0;
          else
              extent_4(12) = (solid_consumption_4(2)/((3*MW.Ni+4*MW.S)/1000)+extent_4(1))*r_4(12)/(2*r_4(3)+r_4(8)+r_4(10)+r_4(12));  %Extent of reaction 12, mol/h
          end;
          if (r_4(15) == 0)
              extent_4(15) = 0;
          else
              extent_4(15) = (solid_consumption_4(8)/((MW.Rh+2*MW.O)/1000)+8*extent_4(7)+2*extent_4(8))/2;        %Extent of reaction 15, mol/h
          end;
          if (r_4(18) == 0)
              extent_4(18) = 0;
          else
              extent_4(18) = (solid_consumption_4(11)/((MW.Ru+2*MW.O)/1000)+8*extent_4(9)+2*extent_4(10))/2;      %Extent of reaction 18, mol/h
          end;
          if (r_4(21) == 0)
              extent_4(21) = 0;
          else
              extent_4(21) = (solid_consumption_4(14)/((MW.Ir+2*MW.O)/1000)+8*extent_4(11)+2*extent_4(12))/2;     %Extent of reaction 21, mol/h
          end;
          %Calculate the amount of species dissolved (generated) using the extents of reaction 
          liquid_generation_4 = zeros(1,7);                                                                                           %Initialise the vector containing the liquid generation terms
          liquid_generation_4(1) = (4*extent_4(4)+extent_4(5)+27*extent_4(7)+27*extent_4(9)+27*extent_4(11))*(MW.Cu/1000);            %Rate of Cu2+ production in autoclave compartment 4, kg/h
          liquid_generation_4(2) = (extent_4(1)+extent_4(2)+6*extent_4(3)+3*extent_4(8)+3*extent_4(10)+3*extent_4(12))*(MW.Ni/1000);  %Rate of Ni2+ production in autoclave compartment 4, kg/h
          liquid_generation_4(3) = (MW.Fe/1000)*(extent_4(6));                                                                        %Rate of Fe2+/Fe3+ production in autoclave compartment 4, kg/h
          liquid_generation_4(4) = (MW.Rh/1000)*(-8*extent_4(7)-2*extent_4(8)+2*extent_4(13)+4*extent_4(14)+2*extent_4(15));          %Rate of Rh3+ production in autoclave compartment 4, kg/h
          liquid_generation_4(5) = (MW.Ru/1000)*(-8*extent_4(9)-2*extent_4(10)+4*extent_4(16)+4*extent_4(17)+2*extent_4(18));         %Rate of Ru3+ production in autoclave compartment 4, kg/h
          liquid_generation_4(6) = (MW.Ir/1000)*(-8*extent_4(11)-2*extent_4(12)+2*extent_4(19)+4*extent_4(20)+2*extent_4(21));        %Rate of Ir3+ production in autoclave compartment 4, kg/h
          liquid_generation_4(7) = ((2*MW.H+MW.S+4*MW.O)/1000)*((-2*extent_4(1)+4*extent_4(3)-8*extent_4(4)-extent_4(6)+8*extent_4(8)+8*extent_4(10)+8*extent_4(12)-12*extent_4(14)-6*extent_4(15)*4*extent_4(16)-12*extent_4(17)-6*extent_4(18)-12*extent_4(20)-6*extent_4(21))/2);%Rate of H2SO4 production in autoclave compartment 4, kg/h
          %If dissolved species are consumed rather than generated, check that the amount consumed does not exceed the available amount
          for i=(1:7)
              if (liquid_generation_4(i)<(-1*X(186)*prep_tank_3_constant*(X(170)^0.5)*X(i+186)))           %Dissolved species is the limiting reagent rather than the solid species
                  liquid_generation_4(i) = -1*X(186)*prep_tank_3_constant*(X(170)^0.5)*X(i+186);           %Set dissolved species consumption equal to amount of dissolved species fed
                  %Recalculate the extents of reaction for the reactions where the dissolved species are the limiting reagents
                  if (i==4)
                      extent_4(7) = (liquid_generation_4(i)/(MW.Rh/1000)-2*extent_4(13)-4*extent_4(14)-2*extent_4(15))*r_4(7)/(-8*r_4(7)-2*r_4(8));
                      extent_4(8) = (liquid_generation_4(i)/(MW.Rh/1000)-2*extent_4(13)-4*extent_4(14)-2*extent_4(15))*r_4(8)/(-8*r_4(7)-2*r_4(8));
                  end;
                  if (i==5)
                      extent_4(9) = (liquid_generation_4(i)/(MW.Ru/1000)-4*extent_4(16)-4*extent_4(17)-2*exten
Niel

Products

1 Answer

Answer by Kaustubha Govind on 19 Jul 2013
Accepted answer

Niel: Don't mean to discourage you, but please post your comments on this thread instead of sending personal messages in reply (we highly discourage contacting contributors personally on this forum since our participation is voluntary and not part of our job responsibilities. Thanks or understanding!)

---Message from Niel follows---

Thank you for pointing out that the errors I received were "compile time" errors. By pressing Cntl+D, I indeed received the same errors. The reason I thought that the problem lies with Simulink is that the same code ran perfectly in Matlab's normal code format. I did make a few changes to the code in order to incorporate it into a Matlab function block. Why would the same piece of code suddenly given errors in Simulink?

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

When your code runs in the MATLAB Function block, it actually generates C code from the MATLAB code and generates a MEX-file from it. It is this MEX-file that is used for execution of the block. The compilation process for generating the C code is a lot stricter than how the MATLAB interpeter works, which is why you're seeing the failure only when using the code in the block. You should however see the same compilation error if you attempt to generate code from the function using MATLAB Coder (outside of Simulink).

Now, coming back to the original question - it looks like your code is too long to read legibly even if you paste it here completely. Are you absolutely sure that the variable that the error complains about is assigned outside of the "if" section as well? To be safe, you could try assigning it to a value at the top of your function, and then let it get reassigned later.

3 Comments

Niel on 21 Jul 2013

Apologies for the personal message. I am new to this forum and I thought it was some kind of Matlab Central PM, and not an email...

I found that the variable in question is indeed in a very big "while" loop.

Kaustubha Govind on 22 Jul 2013

Usually, searching the documentation is a good start. For example, I searched for the error you ran into and found this page: http://www.mathworks.com/help/simulink/ug/best-practices-for-defining-variables-for-c-c-code-generation.html

In rare cases, when you don't find help in the documentation, this forum is a good place to turn to.

Niel on 22 Jul 2013

I appreciate your help - thank you.

Kaustubha Govind

Contact us