diff --git a/Web/coolprop/HighLevelAPI.rst b/Web/coolprop/HighLevelAPI.rst
index 8d95b994..a4b48d53 100644
--- a/Web/coolprop/HighLevelAPI.rst
+++ b/Web/coolprop/HighLevelAPI.rst
@@ -4,6 +4,9 @@
 High-Level Interface
 ********************
 
+PropsSI function
+----------------
+
 For many users, all that is needed is a simple call to the ``PropsSI`` function for pure fluids, pseudo-pure fluids and mixtures.  For humid air properties, see :ref:`Humid air properties <Humid-Air>`.  An example using ``PropsSI``:
 
 .. ipython::
@@ -24,6 +27,107 @@ All :ref:`the wrappers <wrappers>` wrap this function in exactly the same way.
 
 For pure and pseudo-pure fluids, two state points are required to fix the state.  The equations of state are based on :math:`T` and :math:`\rho` as state variables, so :math:`T, \rho` will always be the fastest inputs.  :math:`P,T` will be a bit slower (3-10 times), and then comes inputs where neither :math:`T` nor :math:`\rho` are given, like :math:`p,h`.  They will be much slower.  If speed is an issue, you can look into table-based interpolation methods using TTSE or bicubic interpolation. 
 
+PhaseSI function
+----------------
+
+It can be useful to know what the phase of a given state point is.  A high-level function called ``PhaseSI`` has been implemented to allow for access to the phase.
+
+    In [1]: import CoolProp
+    
+    In [5]: CoolProp.CoolProp.PhaseSI('Phase','P',101325,'Q',0,'Water')
+
+The phase index (as floating point number) can also be obtained using the PropsSI function. In python you would do:
+
+.. ipython::
+
+    In [1]: import CoolProp
+    
+    In [5]: CoolProp.CoolProp.PropsSI('Phase','P',101325,'Q',0,'Water')
+    
+where you can obtain the integer indices corresponding to the phase flags using the ``get_parameter_index`` function:
+
+.. ipython::
+
+    In [1]: import CoolProp
+
+    In [6]: CoolProp.CoolProp.get_parameter_index('phase_twophase')
+    
+    # Or for liquid
+    In [6]: CoolProp.CoolProp.get_parameter_index('phase_liquid')
+    
+For a given fluid, the phase can be plotted in T-p coordinates:
+
+.. plot::
+
+    import matplotlib
+    import numpy as np
+    import CoolProp as CP
+    import matplotlib.pyplot as plt
+    import scipy.interpolate
+
+    Water = CP.AbstractState("HEOS", "Water")
+    pc = Water.keyed_output(CP.iP_critical)
+    Tc = Water.keyed_output(CP.iT_critical)
+    Tmin = 200
+    Tmax = 1000
+    pmax = Water.keyed_output(CP.iP_max)
+    pt = 611.657
+    Tt = 273.16
+    fillcolor = 'g'
+
+    fig = plt.figure(figsize = (6,6))
+    ax = fig.add_subplot(111)
+    lw = 3
+
+    # --------------
+    # Melting curve
+    # --------------
+    melt_args = dict(lw = lw, solid_capstyle = 'round')
+    TT = []
+    PP = list(np.logspace(np.log10(pt), np.log10(pmax),1000))
+    for p in PP:
+        TT.append(Water.melting_line(CP.iT, CP.iP, p))
+
+    #Zone VI
+    for T in np.linspace(max(TT), 355):
+        TT.append(T)
+        theta = T/273.31
+        pi = 1-1.07476*(1-theta**4.6)
+        p = pi*632.4e6
+        PP.append(p)
+
+    plt.plot(TT,PP,'darkblue',**melt_args)
+
+    # ----------------
+    # Saturation curve
+    # ----------------
+    Ts = np.linspace(273.16, Tc, 1000)
+    ps = CP.CoolProp.PropsSI('P','T',Ts,'Q',[0]*len(Ts),'Water',[1])
+
+    # ------
+    # Labels
+    # ------
+
+    plt.plot(Ts,ps,'orange',lw = lw, solid_capstyle = 'round')
+
+    # Critical lines
+    plt.axvline(Tc, dashes = [2, 2])
+    plt.axhline(pc, dashes = [2, 2])
+
+    # Labels
+    plt.text(850, 1e8, 'supercritical',ha= 'center')
+    plt.text(850, 1e5, 'supercritical_gas', rotation = 90)
+    plt.text(450, 1e8, 'supercritical_liquid', rotation = 0, ha = 'center')
+    plt.text(350, 3e6, 'liquid', rotation = 45)
+    plt.text(450, 5e4, 'gas', rotation = 45)
+
+    plt.ylim(611,1e9)
+    plt.gca().set_yscale('log')    
+    plt.gca().set_xlim(240, 1000)
+    plt.ylabel('Pressure [Pa]')
+    plt.xlabel('Temperature [K]')
+    plt.tight_layout()
+
 Code
 ----
 .. literalinclude:: snippets/propssi.cxx
diff --git a/wrappers/Python/CoolProp/CoolProp.pxd b/wrappers/Python/CoolProp/CoolProp.pxd
index 6b2c03b5..5e7d9ae8 100644
--- a/wrappers/Python/CoolProp/CoolProp.pxd
+++ b/wrappers/Python/CoolProp/CoolProp.pxd
@@ -29,6 +29,7 @@ cdef extern from "CoolPropLib.h":
 cdef extern from "CoolProp.h" namespace "CoolProp":
     double _Props1SI "CoolProp::Props1SI"(string Ref, string Output)
     double _PropsSI "CoolProp::PropsSI"(string Output, string Name1, double Prop1, string Name2, double Prop2, string FluidName) 
+    string _PhaseSI "CoolProp::PhaseSI"(string Name1, double Prop1, string Name2, double Prop2, string FluidName) 
     vector[double] _PropsSI "CoolProp::PropsSI"(string Output, string Name1, vector[double] Prop1, string Name2, vector[double] Prop2, string FluidName, vector[double] fractions)
     vector[double] _PropsSII "CoolProp::PropsSI"(string Output, string Name1, vector[double] Prop1, string Name2, vector[double] Prop2, string FluidName)
     string _get_global_param_string "CoolProp::get_global_param_string"(string ParamName) except +
diff --git a/wrappers/Python/CoolProp/CoolProp.pyx b/wrappers/Python/CoolProp/CoolProp.pyx
index c0b1399b..9df98166 100644
--- a/wrappers/Python/CoolProp/CoolProp.pyx
+++ b/wrappers/Python/CoolProp/CoolProp.pyx
@@ -170,6 +170,14 @@ cpdef Props(in1, in2, in3 = None, in4 = None, in5 = None, in6 = None):
     else:
         return val
 
+cpdef PhaseSI(in1, in2, in3, in4, in5):
+    """
+    A Python wrapper of C++ function :cpapi:`CoolProp::PhaseSI`
+    
+    Does not support vectorization of the inputs like PropsSI
+    """
+    return _PhaseSI(in1, in2, in3, in4, in5)
+
 cpdef PropsSI(in1, in2, in3 = None, in4 = None, in5 = None, in6 = None, in7 = None):
     """
     A Python wrapper of C++ function :cpapi:`CoolProp::PropsSI` .