function [c,ceq] = constraints(u)

% Mass flows (m) in kg/s
% Temperatures in Celcius
% mCp's in W/K
% UA's in W/K
% n = number of shells

global T0;
global m0;
global T1; global Th1; global T2; global Th2;

%% Parameters for stream A

m0A     = m0*u(1);
cp0A    = 2570.944432;
wcA     = m0A*cp0A;

nA      = 4;
whA     = 48187.09;
Th1inA  = 300.416;
UAA     = 131025.48;

[T1A Th1A] = tempCalc(T0,Th1inA,wcA,whA,UAA,nA);

%% Parameters for stream B

m0B     = m0*u(2);

cp0B1   = 2445.743628;
wc1B    = m0B*cp0B1;

wh1B    = 92541.46;
nB1     = 2;
UA1B    = 102709.10;

cp0B2    = 2618.432483;
wc2B    = m0B*cp0B2;
nB2     = 2;
wh2B    = 129499.93;
UA2B    = 88644.82;

Th2inB  = 270.254;

tparB = [T0 wc1B wh1B UA1B nB1 Th2inB wc2B wh2B UA2B nB2 1];

%Using fsolve to find Th1in=Th2
Th1inB  = fsolve(@(x) FindTh1in(x,tparB),232,optimset('Display','off'));

% Temperatures calculated in the function FindTh1in
T1B = T1; Th1B = Th1; T2B = T2; Th2B = Th2; 

%% Parameters for stream C

m0C     = m0*u(3);

cp0C1   = 2407.952102;
wc1C    = m0C*cp0C1;

wh1C    = 35816.4165;
UA1C    = 84642.10;
nC1     = 2;

cp0C2   = 2584.355153;
wc2C    = m0C*cp0C2;
wh2C    = 38747.8648;
UA2C    = 133605.46;
nC2     = 3;

Th2inC  = 245.375;

tparC = [T0 wc1C wh1C UA1C nC1 Th2inC wc2C wh2C UA2C nC2 2];

%Using fsolve to find Th1in=Th2
Th1inC  = fsolve(@(x) FindTh1in(x,tparC),175.5,optimset('Display','off'));

% Temperatures calculated in the function FindTh1in
T1C = T1; Th1C = Th1; T2C = T2; Th2C = Th2; 

%% Parameters for stream D

m0D     = m0*u(4);
cp0D1   = 2483.068751;
wcD     = m0D*cp0D1;

n1D     = 4;
wh1D    = 118477.94;
Th1inD  = 226.181;
UA1D    = 132831.92;

[T1D Th1D] = tempCalc(T0,Th1inD,wcD,wh1D,UA1D,n1D);

cp0D2   = 2649.127849;
wc2D    = m0D*cp0D2;
n2D     = 1;
wh2D    = 33935.43;
Th2inD  = 273.813;
UA2D = 41563.90;

[T2D Th2D] = tempCalc(T1D,Th2inD,wc2D,wh2D,UA2D,n2D);

%% Parameters for stream E

m0E     = m0*u(5);
cp0E    = 2527.081131;
wcE     = m0E*cp0E;

nE      = 4;
whE     = 79890.13;
Th1inE  = 256.408;
UAE     = 190981.94;

[T1E Th1E] = tempCalc(T0,Th1inE,wcE,whE,UAE,nE);

%% Parameters for stream F

m0F     = m0*u(6);
cp0F1   = 2382.242363;
wcF     = m0F*cp0F1;

n1F     = 1;
wh1F    = 47248.68;
Th1inF  = 203.201;
UA1F    = 49437.00;

[T1F Th1F] = tempCalc(T0,Th1inF,wcF,wh1F,UA1F,n1F);

cp0F2   = 2535.118619;
wc2F    = m0F*cp0F2; 
n2F     = 4;
wh2F    = 175073.77;
Th2inF  = 248.07;
UA2F    = 224053.26;

[T2F Th2F] = tempCalc(T1F,Th2inF,wc2F,wh2F,UA2F,n2F);


%% Saving temperatures
global Tcold
global Thot

Tcold=[T1A T1B T2B T1C T2C T1D T2D T1E T1F T2F];
Thot=[Th1A Th1B Th2B Th1C Th2C Th1D Th2D Th1E Th1F Th2F];

%% Controlled variables
fA = (T1A-T0).^2./(Th1inA-T0);
fB = ((Th2inB-T2B)./(Th1inB-T0)-1).*(T1B-T0).^2./(Th2inB-T1B)+(T2B-T0).^2./(Th2inB-T1B);
fC = ((Th2inC-T2C)./(Th1inC-T0)-1).*(T1C-T0).^2./(Th2inC-T1C)+(T2C-T0).^2./(Th2inC-T1C);
fD = ((Th2inD-T2D)./(Th1inD-T0)-1).*(T1D-T0).^2./(Th2inD-T1D)+(T2D-T0).^2./(Th2inD-T1D);
fE = (T1E-T0).^2./(Th1inE-T0);
fF = ((Th2inF-T2F)./(Th1inF-T0)-1).*(T1F-T0).^2./(Th2inF-T1F)+(T2F-T0).^2./(Th2inF-T1F);

cA = fA - fF;
cB = fB - fF;
cC = fC - fF;
cD = fD - fF;
cE = fE - fF;

%% Constraints
c=[]; %No inequality constraints

% Equality constraints = 0
ceq=[cA
     cB
     cC
     cD
     cE
     sum(u)-1];
 
