diff --git a/Simulator/Files/Thermodynamic_Functions/package.mo b/Simulator/Files/Thermodynamic_Functions/package.mo
index a582aef..cd49600 100644
--- a/Simulator/Files/Thermodynamic_Functions/package.mo
+++ b/Simulator/Files/Thermodynamic_Functions/package.mo
@@ -86,6 +86,7 @@ package Thermodynamic_Functions
     input Real P;
     input Real x;
     input Real y;
+    input Real Phase_Density;
     output Real Sliq, Svap;
   protected
     parameter Real Tref = 298.15, Pref = 101325;
@@ -103,7 +104,7 @@ package Thermodynamic_Functions
       Entr := -(T - 298.15) * (Simulator.Files.Thermodynamic_Functions.VapCpId(VapCp, T) / (2 * T) + sum(Cp[:]) + Simulator.Files.Thermodynamic_Functions.VapCpId(VapCp, 298.15) / (2 * 298.15)) / n;
     end if;
     if x > 0 and y > 0 then
-      Sliq := Entr - R * log(P / Pref) - R * log(x) - HV(HOV, Tc, T) / T;
+      Sliq := Entr - R * log(P / Pref) - R * log(x) - HV(HOV, Tc, T) / T + P/1000/Phase_Density/T;
       Svap := Entr - R * log(P / Pref) - R * log(y);
     elseif x <= 0 and y <= 0 then
       Sliq := 0;
@@ -112,7 +113,7 @@ package Thermodynamic_Functions
       Sliq := 0;
       Svap := Entr - R * log(P / Pref) - R * log(y);
     elseif y == 0 then
-      Sliq := Entr - R * log(P / Pref) - R * log(x) - HV(HOV, Tc, T) / T;
+      Sliq := Entr - R * log(P / Pref) - R * log(x) - HV(HOV, Tc, T) / T + P/1000/Phase_Density/T;
       Svap := 0;
     else
       Sliq := 0;
@@ -150,6 +151,7 @@ package Thermodynamic_Functions
 
 
 
+
 
 
   function Dens
diff --git a/Simulator/Streams/Material_Stream.mo b/Simulator/Streams/Material_Stream.mo
index 2472e28..fd7fda6 100644
--- a/Simulator/Streams/Material_Stream.mo
+++ b/Simulator/Streams/Material_Stream.mo
@@ -14,6 +14,8 @@ model Material_Stream
   Real phasMolSpHeat[3] "phase Molar Specific Heat", compMolSpHeat[3, NOC] "Component Molar Specific Heat";
   Real phasMolEnth[3] "Phase Molar Enthalpy", compMolEnth[3, NOC] "Component Molar Enthalpy";
   Real phasMolEntr[3] "Phase Molar Entropy", compMolEntr[3, NOC] "Component Molar Entropy";
+  Real Liquid_Phase_Density;
+  Real LiqDens[NOC];
   Simulator.Files.Connection.matConn inlet(connNOC = NOC) annotation(
     Placement(visible = true, transformation(origin = {-100, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0), iconTransformation(origin = {-100, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
   Simulator.Files.Connection.matConn outlet(connNOC = NOC) annotation(
@@ -74,7 +76,7 @@ equation
     compMolSpHeat[3, i] = Thermodynamic_Functions.VapCpId(comp[i].VapCp, T);
     compMolEnth[2, i] = Thermodynamic_Functions.HLiqId(comp[i].SH, comp[i].VapCp, comp[i].HOV, comp[i].Tc, T);
     compMolEnth[3, i] = Thermodynamic_Functions.HVapId(comp[i].SH, comp[i].VapCp, comp[i].HOV, comp[i].Tc, T);
-    (compMolEntr[2, i], compMolEntr[3, i]) = Thermodynamic_Functions.SId(comp[i].AS, comp[i].VapCp, comp[i].HOV, comp[i].Tb, comp[i].Tc, T, P, compMolFrac[2, i], compMolFrac[3, i]);
+    (compMolEntr[2, i], compMolEntr[3, i]) = Thermodynamic_Functions.SId(comp[i].AS, comp[i].VapCp, comp[i].HOV, comp[i].Tb, comp[i].Tc, T, P, compMolFrac[2, i], compMolFrac[3, i],Liquid_Phase_Density);
   end for;
   for i in 2:3 loop
     phasMolSpHeat[i] = sum(compMolFrac[i, :] .* compMolSpHeat[i, :]) + resMolSpHeat[i];
@@ -111,5 +113,8 @@ algorithm
   for i in 1:NOC loop
     MW[:] := MW[:] + comp[i].MW * compMolFrac[:, i];
   end for;
+   LiqDens[:] := Thermodynamic_Functions.Density_Racket(NOC, T, P, comp[:].Pc, comp[:].Tc, comp[:].Racketparam, comp[:].AF, comp[:].MW, Psat[:]);
+  
+  Liquid_Phase_Density := 1 / sum(compMasFrac[2, :] ./ LiqDens[:]) / MW[2];
 
 end Material_Stream;