function [x, sc, d] = initialize(Ths, Tcs, DTmin, n, d)

T_to_valve = Ths+DTmin;
T_to_comp = Tcs-DTmin;
Tsat_con = T_to_valve+5;  % 5 degrees subcooling
Tsat_vap = T_to_comp;
Ph = psat(Tsat_con, d);
Pl = psat(Tsat_vap, d);

% Compressor
ucom = [Pl, Ph, n, T_to_comp];
[xcom, sc.com] = initial_compressor(ucom, d);

% Condenser
ucon = [Ths, n, xcom(2), T_to_valve, Tsat_con];
[xcon, sc.con] = initial_condenser(ucon, d);

% Valve
rho_to_valve = xcon(end-1)*rhol(d)+(1-xcon(end-1))*rhov(xcon(5), Ph, d);
h_to_valve = xcon(end-1)*d.Cpl*(xcon(5)-d.Tref)+(1-xcon(end-1))*(d.Cpg*(xcon(5)-d.Tref)+d.hvap);
uval = [rho_to_valve; Ph-Pl; h_to_valve; n];
[xval, sc.val, d] = initial_valve(uval, d);

% Evaporator
uvap = [Tcs; n; xval(end); Tsat_vap];
[xvap, sc.vap] = initial_evaporator(uvap, d);

x = [xcon; xval; xvap; xcom];

function [xvap, sc] = initial_evaporator(u, d)
Tcs = u(1);
n = u(2);
hin = u(3);
Tsat = u(4);

P = psat(Tsat, d);
T = [Tsat];

h = d.Cpg*(T-d.Tref)+d.hvap;

Q = n*(h-hin);
A = Q/(d.Uvap*(Tcs-T));
V_l = 0.5*d.Vvap;
V_v = d.Vvap-V_l;
N_v = V_v*rhov(T, P, d);
N_l = V_l*rhol(d);
xvap = [T; A; Q; N_l; N_v; V_l; P];
sc = scale(xvap, d, 'vap');
u = [n; Tcs; hin];
aux.d=d;
aux.sc=sc;
aux.u=u;
res_vap = vaporizer(xvap./sc.x, aux);

function [x, sc, d] = initial_valve(u, d)
rho = u(1);
DP = u(2);
hin = u(3);
n = u(4);
z = 0.5;
d.Cv = n/(z*sqrt(DP*rho/d.MW));
hout = hin;
x = [z; n; hout];
sc = scale(x, d, 'val');
u = [DP; rho; hin];
aux.d = d;
aux.sc = sc;
aux.u = u;
res_val = valve(x./sc.x, aux);


function [x, sc] = initial_compressor(u, d)
Pl = u(1);
Ph = u(2);
n = u(3);
Tin = u(4);

Tout = Tin*(Ph/Pl)^(d.R/d.Cpg);
Ws = n*d.Cpg*(Tout-Tin);

x = [Ws; Tout];
sc = scale(x, d, 'com');
u = [Tin; n; Ph; Pl];
aux.sc = sc;
aux.d = d;
aux.u = u;
res_com = compressor(x./sc.x, aux);

function [xcon, sc] = initial_condenser(u, d)
Ths = u(1);
n = u(2);
Tin = u(3);
Tout = u(4);
Tsat = u(5);
P = psat(Tsat, d);
f = 1; % Sub-cooled
hin = d.Cpg*(Tin-d.Tref)+d.hvap;
hout = f*d.Cpl*(Tout-d.Tref)+(1-f)*(d.Cpg*(Tout-d.Tref)+d.hvap);

x0 = [[Tsat:(Tout-Tsat)/(d.ncon_sub-1):Tout]'; 1]; 

options = optimset('fsolve');
[x, fval] = fsolve(@subcooling, x0, options, d, Tsat, Tout, n, Ths);
x = [Tsat; Tsat; x(1:end-1)];
T = [Tin; x];
h = [d.Cpg*(T(1:2)-d.Tref)+d.hvap;
     f*d.Cpl*(T(3:6)-d.Tref)+(1-f)*(d.Cpg*(T(3:6)-d.Tref))];
Q = n*(h(2:6)-h(1:5));
DT = Ths-[(T(1)+T(2))/2; T(3:6)];
A = Q./(d.Ucon.*DT);
V_v = [1;0.5].*A(1:2)/sum(A)*d.Vcon;
V_l = [0.5;1].*A(2:3)/sum(A)*d.Vcon;
N_v = V_v(1:2).*rhov(T(1:2), P*ones(2,1),d);
N_l = [V_l(1);
       V_l(2)*ones(3,1)].*rhol(d);
xcon = [x; A; Q; N_l; N_v; f; P];
sc = scale(xcon, d, 'con');
u = [Tin; n; Ths];
aux.d=d;
aux.u=u;
aux.sc=sc;

res_con = condenser(xcon./sc.x, aux);
