function f = calc_f(x,Param)
% contains a system of equations that describes the steady-state operation
% of a heat exchanger network

% INPUT
% x    vector of guessed values of the root
% Param     structure of other parameters in the equations
% (both x and Param are assigned to physical variables in the code)
% 
% OUTPUT
% f    vector of the equations computed at x.  If = 0, x is root
% 
% ASSIGNMENT
% *************
% edit the file to set the appropriate variables = to elements of x
% *************
% inlet flow streams
W1 = Param.w1;  % flow in kg s-1
Cp1 = 2500;  % specific heat in J kg-1 K-1
T1 = Param.T1;   % inlet temperature in C
W2 = Param.w2;  % flow in kg s-1
Cp2 = 2500;  % specific heat in J kg-1 K-1
T2 = Param.T2;   % inlet temperature in C
W3 = Param.w3;  % flow in kg s-1
Cp3 = 2500;  % specific heat in J kg-1 K-1
T3 = Param.T3;   % inlet temperature in C

% heat exchangers
A1 = 117;  % heat transfer surface in m2
U1 = 250; % heat transfer coef in W m-2 K-1
A2 = 117;  % heat transfer surface in m2
U2 = 250; % heat transfer coef in W m-2 K-1
A3 = 117;  % heat transfer surface in m2
U3 = 250; % heat transfer coef in W m-2 K-1

% intermediate and exit temperatures
T4 = x(4);   % temperature in C
T5 = x(5);   % temperature in C
T6 = x(1);   % temperature in C
T7 = x(6);   % temperature in C
T8 = x(2);   % temperature in C
T9 = x(3);   % temperature in C

% create output vector f according to size of input matrix
% it should equal the number of equations below!
n = length(x);
f = zeros(n,1);  % creates vector f and fills with zeroes

% the network model
% exchanger E101
DTLM1 = ((T1 - T8) - (T4 - T7))/log((T1 - T8) / (T4 - T7));
if isnan(DTLM1)  % catch a uniform temperature difference
    DTLM1 = T1 - T8;
end
Q1 = U1*A1*DTLM1;
f(1) = -W1*Cp1*T1 + W1*Cp1*T4 + Q1;
f(2) = -W3*Cp3*T7 - Q1 + W3*Cp3*T8;

% exchanger E102
DTLM2 = ((T4 - T9) - (T5 - T2))/log((T4 - T9) / (T5 - T2));
if isnan(DTLM2)  % catch a uniform temperature difference
    DTLM2 = T4 - T9;
end
Q2 = U2*A2*DTLM2;
f(3) = -W1*Cp1*T4 + W1*Cp1*T5 + Q2;
f(4) = -W2*Cp2*T2 - Q2 + W2*Cp2*T9;

% exchanger E103
DTLM3 = ((T5 - T7) - (T6 - T3))/log((T5 - T7) / (T6 - T3));
if isnan(DTLM3)  % catch a uniform temperature difference
    DTLM3 = T5 - T7;
end    
Q3 = U3*A3*DTLM3;
f(5) = -W1*Cp1*T5 + W1*Cp1*T6 + Q3;
f(6) = -W3*Cp3*T3 - Q3 + W3*Cp3*T7;


% Dr. Barry S. Johnston, Copyright 2004.