Plots.py

# -*- coding: utf-8 -*-
from __future__ import print_function, division, absolute_import

import numpy as np
import warnings

import CoolProp
from CoolProp.Plots.Common import IsoLine, BasePlot, interpolate_values_1d
from CoolProp.Plots.SimpleCycles import StateContainer


class PropertyPlot(BasePlot):
    def __init__(self, fluid_name, graph_type, **kwargs):
        """
        Create graph for the specified fluid properties

        Parameters
        ----------
        fluid_name : string or AbstractState
            The name of the fluid to be plotted or a state instance
        graph_type : string
            The graph type to be plotted, like \"PH\" or \"TS\"
        axis : :func:`matplotlib.pyplot.gca()`, Optional
            The current axis system to be plotted to.
            Default: create a new axis system
        fig : :func:`matplotlib.pyplot.figure()`, Optional
            The current figure to be plotted to.
            Default: create a new figure
        unit_system : string, ['EUR','KSI','SI']
            Select the units used for the plotting.  'EUR' is bar, kJ, C; 'KSI' is kPa, kJ, K; 'SI' is Pa, J, K
        tp_limits : string, ['NONE','DEF','ACHP','ORC']
            Select the limits in T and p.
        reciprocal_density : bool
            NOT IMPLEMENTED: If True, 1/rho will be plotted instead of rho

        Examples
        --------
        >>> from CoolProp.Plots import PropertyPlot
        >>> plot = PropertyPlot('HEOS::Water', 'TS')
        >>> plot.calc_isolines()
        >>> plot.show()

        >>> import CoolProp
        >>> from CoolProp.Plots import PropertyPlot
        >>> plot = PropertyPlot('HEOS::R134a', 'PH', unit_system='EUR', tp_limits='ACHP')
        >>> plot.calc_isolines(CoolProp.iQ, num=11)
        >>> plot.calc_isolines(CoolProp.iT, num=25)
        >>> plot.calc_isolines(CoolProp.iSmass, num=15)
        >>> plot.show()

        >>> import CoolProp
        >>> from CoolProp.Plots import PropertyPlot
        >>> plot = PropertyPlot('HEOS::R245fa', 'TS', unit_system='EUR', tp_limits='ORC')
        >>> plot.calc_isolines(CoolProp.iQ, num=11)
        >>> plot.calc_isolines(CoolProp.iP, iso_range=[1,50], num=10, rounding=True)
        >>> plot.draw()
        >>> plot.isolines.clear()
        >>> plot.props[CoolProp.iP]['color'] = 'green'
        >>> plot.props[CoolProp.iP]['lw'] = '0.5'
        >>> plot.calc_isolines(CoolProp.iP, iso_range=[1,50], num=10, rounding=False)
        >>> plot.show()

        .. note::

            See the online documentation for a list of the available fluids and
            graph types
        """
        super(PropertyPlot, self).__init__(fluid_name, graph_type, **kwargs)
        self._isolines = {}
        #self._plines = {}
        #self._ppoints = {}
        self.get_axis_limits()
        self._plot_default_annotations()

    @property
    def isolines(self): return self._isolines
    # @property
    #def plines(self): return self._plines
    # @property
    #def ppoints(self): return self._ppoints

    def show(self):
        self.draw()
        super(PropertyPlot, self).show()

    def savefig(self, *args, **kwargs):
        self.draw()
        super(PropertyPlot, self).savefig(*args, **kwargs)

    def _plotRound(self, values):
        """
        A function round an array-like object while maintaining the
        amount of entries. This is needed for the isolines since we
        want the labels to look pretty (=rounding), but we do not
        know the spacing of the lines. A fixed number of digits after
        rounding might lead to reduced array size.
        """
        inVal = np.unique(np.sort(np.array(values)))
        output = inVal[1:] * 0.0
        digits = -1
        limit = 10
        lim = inVal * 0.0 + 10
        # remove less from the numbers until same length,
        # more than 10 significant digits does not really
        # make sense, does it?
        while len(inVal) > len(output) and digits < limit:
            digits += 1
            val = (np.around(np.log10(np.abs(inVal))) * -1) + digits + 1
            val = np.where(val < lim, val, lim)
            val = np.where(val > -lim, val, -lim)
            output = np.zeros(inVal.shape)
            for i in range(len(inVal)):
                output[i] = np.around(inVal[i], decimals=int(val[i]))
            output = np.unique(output)
        return output

    def calc_isolines(self, iso_type=None, iso_range=None, num=15, rounding=False, points=250):
        """Calculate lines with constant values of type 'iso_type' in terms of x and y as
        defined by the plot object. 'iso_range' either is a collection of values or
        simply the minimum and maximum value between which 'num' lines get calculated.
        The 'rounding' parameter can be used to generate prettier labels if needed.
        """

        if iso_type is None or iso_type == 'all':
            for i_type in IsoLine.XY_SWITCH:
                if IsoLine.XY_SWITCH[i_type].get(self.y_index * 10 + self.x_index, None) is not None:
                    self.calc_isolines(i_type, None, num, rounding, points)
            return

        if iso_range is None:
            if iso_type is CoolProp.iQ:
                iso_range = [0.0, 1.0]
            else:
                limits = self.get_axis_limits(iso_type, CoolProp.iT)
                iso_range = [limits[0], limits[1]]

        if len(iso_range) <= 1 and num != 1:
            raise ValueError('You have to provide two values for the iso_range, {0} is not valid.'.format(iso_range))

        if len(iso_range) == 2 and (num is None or num < 2):
            raise ValueError('Please specify the number of isoline you want e.g. num=10.')

        iso_range = np.sort(np.unique(iso_range))
        # Generate iso ranges
        if len(iso_range) == 2:
            iso_range = self.generate_ranges(iso_type, iso_range[0], iso_range[1], num)
        if rounding:
            iso_range = self._plotRound(iso_range)

        # Limits are already in SI units
        limits = self._get_axis_limits()

        ixrange = self.generate_ranges(self._x_index, limits[0], limits[1], points)
        iyrange = self.generate_ranges(self._y_index, limits[2], limits[3], points)

        dim = self._system[iso_type]

        lines = self.isolines.get(iso_type, [])
        for i in range(num):
            lines.append(IsoLine(iso_type, self._x_index, self._y_index, value=dim.to_SI(iso_range[i]), state=self._state))
            lines[-1].calc_range(ixrange, iyrange)
            lines[-1].sanitize_data()
        self.isolines[iso_type] = lines
        return

    def draw_isolines(self):
        dimx = self._system[self._x_index]
        dimy = self._system[self._y_index]

        sat_props = self.props[CoolProp.iQ].copy()
        if 'lw' in sat_props: sat_props['lw'] *= 2.0
        else: sat_props['lw'] = 1.0
        if 'alpha' in sat_props: min([sat_props['alpha'] * 2.0, 1.0])
        else: sat_props['alpha'] = 1.0

        for i in self.isolines:
            props = self.props[i]
            dew = None; bub = None
            xcrit = None; ycrit = None
            if i == CoolProp.iQ:
                for line in self.isolines[i]:
                    if line.value == 0.0: bub = line
                    elif line.value == 1.0: dew = line
                if dew is not None and bub is not None:
                    xmin, xmax, ymin, ymax = self.get_axis_limits()
                    xmin = dimx.to_SI(xmin)
                    xmax = dimx.to_SI(xmax)
                    ymin = dimy.to_SI(ymin)
                    ymax = dimy.to_SI(ymax)
                    dx = xmax - xmin
                    dy = ymax - ymin
                    dew_filter = np.logical_and(np.isfinite(dew.x), np.isfinite(dew.y))
                    #dew_filter = np.logical_and(dew_filter,dew.x>dew.x[-1])
                    stp = min([dew_filter.size, 10])
                    dew_filter[0:-stp] = False
                    bub_filter = np.logical_and(np.isfinite(bub.x), np.isfinite(bub.y))

                    if self._x_index == CoolProp.iP or self._x_index == CoolProp.iDmass:
                        filter_x = lambda x: np.log10(x)
                    else:
                        filter_x = lambda x: x
                    if self._y_index == CoolProp.iP or self._y_index == CoolProp.iDmass:
                        filter_y = lambda y: np.log10(y)
                    else:
                        filter_y = lambda y: y

                    if (  # (filter_x(dew.x[dew_filter][-1])-filter_x(bub.x[bub_filter][-1])) > 0.010*filter_x(dx) and
                        (filter_x(dew.x[dew_filter][-1]) - filter_x(bub.x[bub_filter][-1])) < 0.050 * filter_x(dx) or
                        (filter_y(dew.y[dew_filter][-1]) - filter_y(bub.y[bub_filter][-1])) < 0.010 * filter_y(dy)):
                        x = np.linspace(bub.x[bub_filter][-1], dew.x[dew_filter][-1], 11)
                        y = interpolate_values_1d(
                          np.append(bub.x[bub_filter], dew.x[dew_filter][::-1]),
                          np.append(bub.y[bub_filter], dew.y[dew_filter][::-1]),
                          x_points=x,
                          kind='cubic')
                        self.axis.plot(dimx.from_SI(x), dimy.from_SI(y), **sat_props)
                        warnings.warn("Detected an incomplete phase envelope, fixing it numerically.")
                        xcrit = x[5]; ycrit = y[5]
                        #Tcrit = self.state.trivial_keyed_output(CoolProp.iT_critical)
                        #Dcrit = self.state.trivial_keyed_output(CoolProp.irhomass_critical)
                        # try:
                        #    self.state.update(CoolProp.DmassT_INPUTS, Dcrit, Tcrit)
                        #    xcrit = self.state.keyed_output(self._x_index)
                        #    ycrit = self.state.keyed_output(self._y_index)
                        # except:
                        #    xcrit = x[5]; ycrit = y[5]
                        #    pass
                        #self.axis.plot(dimx.from_SI(np.array([bub.x[bub_filter][-1], dew.x[dew_filter][-1]])),dimy.from_SI(np.array([bub.y[bub_filter][-1], dew.y[dew_filter][-1]])),'o')
            for line in self.isolines[i]:
                if line.i_index == CoolProp.iQ:
                    if line.value == 0.0 or line.value == 1.0:
                        self.axis.plot(dimx.from_SI(line.x), dimy.from_SI(line.y), **sat_props)
                    else:
                        if xcrit is not None and ycrit is not None:
                            self.axis.plot(dimx.from_SI(np.append(line.x, xcrit)), dimy.from_SI(np.append(line.y, ycrit)), **props)
                            # try:
                            #    x = np.append(line.x,[xcrit])
                            #    y = np.append(line.y,[ycrit])
                            #    fltr = np.logical_and(np.isfinite(x),np.isfinite(y))
                            #    f = interp1d(x[fltr][-3:],y[fltr][-3:],kind='linear') # could also be quadratic
                            #    x = np.linspace(x[fltr][-2], x[fltr][-1], 5)
                            #    y = f(x)
                            #    #f = interp1d(y[fltr][-5:],x[fltr][-5:],kind='cubic')
                            #    #y = np.linspace(y[fltr][-2], y[fltr][-1], 5)
                            #    #x = f(y)
                            #    self.axis.plot(dimx.from_SI(np.append(line.x,x)),dimy.from_SI(np.append(line.y,y)),**props)
                            # except:
                            #    self.axis.plot(dimx.from_SI(np.append(line.x,xcrit)),dimy.from_SI(np.append(line.y,ycrit)),**props)
                            #    pass
                else:
                    self.axis.plot(dimx.from_SI(line.x), dimy.from_SI(line.y), **props)

    def draw(self):
        self.get_axis_limits()
        self.draw_isolines()

    # def label_isolines(self, dx=0.075, dy=0.100):
    #    [xmin, xmax, ymin, ymax] = self.get_axis_limits()
    #    for i in self.isolines:
    #         for line in self.isolines[i]:
    #             if self.get_x_y_dydx(xv, yv, x)

    def draw_process(self, statecontainer, points=None, line_opts=None):
        """ Draw process or cycle from x and y values in axis units

        Parameters
        ----------
        statecontainer : CoolProp.Plots.SimpleCycles.StateContainer()
            A state container object that contains all the information required to draw the process.
            Note that points that appear several times get added to a special of highlighted points.
        line_opts : dict
            Line options (please see :func:`matplotlib.pyplot.plot`), optional
            Use this parameter to pass a label for the legend.

        Examples
        --------
        >>> import CoolProp
        >>> from CoolProp.Plots import PropertyPlot
        >>> pp = PropertyPlot('HEOS::Water', 'TS', unit_system='EUR')
        >>> pp.calc_isolines(CoolProp.iP        )
        >>> pp.calc_isolines(CoolProp.iHmass    )
        >>> pp.calc_isolines(CoolProp.iQ, num=11)
        >>> cycle = SimpleRankineCycle('HEOS::Water', 'TS', unit_system='EUR')
        >>> T0 = 300
        >>> pp.state.update(CoolProp.QT_INPUTS,0.0,T0+15)
        >>> p0 = pp.state.keyed_output(CoolProp.iP)
        >>> T2 = 700
        >>> pp.state.update(CoolProp.QT_INPUTS,1.0,T2-150)
        >>> p2 = pp.state.keyed_output(CoolProp.iP)
        >>> cycle.simple_solve(T0, p0, T2, p2, 0.7, 0.8, SI=True)
        >>> cycle.steps = 50
        >>> sc = cycle.get_state_changes()
        >>> pp.draw_process(sc)
        >>> # The same calculation can be carried out in another unit system:
        >>> cycle.simple_solve(T0-273.15-10, p0/1e5, T2-273.15+50, p2/1e5-5, 0.7, 0.8, SI=False)
        >>> sc2 = cycle.get_state_changes()
        >>> pp.draw_process(sc2, line_opts={'color':'blue', 'lw':1.5})
        >>> pp.show()

        """
        warnings.warn("You called the function \"draw_process\", which is not tested.", UserWarning)

        # Default values
        line_opts = line_opts or {'color': 'r', 'lw': 1.5}

        dimx = self.system[self.x_index]
        dimy = self.system[self.y_index]

        marker = line_opts.pop('marker', 'o')
        style = line_opts.pop('linestyle', 'solid')
        style = line_opts.pop('ls', style)

        if points is None: points = StateContainer()

        xdata = []
        ydata = []
        old = statecontainer[len(statecontainer) - 1]
        for i in statecontainer:
            point = statecontainer[i]
            if point == old:
                points.append(point)
                old = point
                continue
            xdata.append(point[self.x_index])
            ydata.append(point[self.y_index])
            old = point
        xdata = dimx.from_SI(np.asarray(xdata))
        ydata = dimy.from_SI(np.asarray(ydata))
        self.axis.plot(xdata, ydata, marker='None', linestyle=style, **line_opts)

        xdata = np.empty(len(points))
        ydata = np.empty(len(points))
        for i in points:
            point = points[i]
            xdata[i] = point[self.x_index]
            ydata[i] = point[self.y_index]
        xdata = dimx.from_SI(np.asarray(xdata))
        ydata = dimy.from_SI(np.asarray(ydata))
        line_opts['label'] = ''
        self.axis.plot(xdata, ydata, marker=marker, linestyle='None', **line_opts)


def InlineLabel(xv, yv, x=None, y=None, axis=None, fig=None):
    warnings.warn("You called the deprecated function \"InlineLabel\", use \"BasePlot.inline_label\".", DeprecationWarning)
    plot = PropertyPlot("water", "TS", figure=fig, axis=axis)
    return plot.inline_label(xv, yv, x, y)


class PropsPlot(PropertyPlot):
    def __init__(self, fluid_name, graph_type, units='KSI', reciprocal_density=False, **kwargs):
        super(PropsPlot, self).__init__(fluid_name, graph_type, unit_system=units, reciprocal_density=reciprocal_density, **kwargs)
        warnings.warn("You called the deprecated class \"PropsPlot\", use \"PropertyPlot\".", DeprecationWarning)


if __name__ == "__main__":
    plot = PropertyPlot('HEOS::n-Pentane', 'PD', unit_system='EUR')  # , reciprocal_density=True)
    plot.calc_isolines(CoolProp.iT)
    plot.calc_isolines(CoolProp.iQ, num=11)
    # plot.calc_isolines(CoolProp.iSmass)
    # plot.calc_isolines(CoolProp.iHmass)
    plot.show()