Skip to content

Instantly share code, notes, and snippets.

@rddaz2013

rddaz2013/gist:f7efd141bfc9f2ce024c

Created February 4, 2016 12:45
Show Gist options
  • Save rddaz2013/f7efd141bfc9f2ce024c to your computer and use it in GitHub Desktop.
Save rddaz2013/f7efd141bfc9f2ce024c to your computer and use it in GitHub Desktop.
Download ZIP
coolprop 5.1.1 old plot.py
Raw
gistfile1.txt
# -*- coding: utf-8 -*-
from __future__ import print_function, absolute_import
import numpy, matplotlib, matplotlib.pyplot, math, re
from scipy.interpolate import interp1d
import CoolProp.CoolProp as CP
from .Common import BasePlot
from scipy import interpolate
from scipy.spatial.kdtree import KDTree
class IsoLine(object):
def __init__(self):
self.DEBUG = False
# direct geometry
self.X = None #
self.Y = None #
self.type = None #
self.value = None #
self.unit = None #
self.opts = None #
def InlineLabel(xv,yv,x = None, y= None, axis = None, fig = None):
"""
This will give the coordinates and rotation required to align a label with
a line on a plot
"""
def ToPixelCoords(xv,yv,axis,fig):
[Axmin,Axmax]=axis.get_xlim()
[Aymin,Aymax]=axis.get_ylim()
DELTAX_axis=Axmax-Axmin
DELTAY_axis=Aymax-Aymin
width=fig.get_figwidth()
height=fig.get_figheight()
pos=axis.get_position().get_points()
[[Fxmin,Fymin],[Fxmax,Fymax]]=pos
DELTAX_fig=width*(Fxmax-Fxmin)
DELTAY_fig=height*(Fymax-Fymin)
#Convert coords to pixels
x=(xv-Axmin)/DELTAX_axis*DELTAX_fig+Fxmin
y=(yv-Aymin)/DELTAY_axis*DELTAY_fig+Fymin
return x,y
def ToDataCoords(xv,yv,axis,fig):
[Axmin,Axmax]=axis.get_xlim()
[Aymin,Aymax]=axis.get_ylim()
DELTAX_axis=Axmax-Axmin
DELTAY_axis=Aymax-Aymin
width=fig.get_figwidth()
height=fig.get_figheight()
pos=axis.get_position().get_points()
[[Fxmin,Fymin],[Fxmax,Fymax]]=pos
DELTAX_fig=(Fxmax-Fxmin)*width
DELTAY_fig=(Fymax-Fymin)*height
#Convert back to measurements
x=(xv-Fxmin)/DELTAX_fig*DELTAX_axis+Axmin
y=(yv-Fymin)/DELTAY_fig*DELTAY_axis+Aymin
return x,y
def get_x_y_dydx(xv,yv,x):
"""Get x and y coordinates and the linear interpolation derivative"""
# Old implementation:
##Get the rotation angle
#f = interp1d(xv, yv)
#y = f(x)
#h = 0.00001*x
#dy_dx = (f(x+h)-f(x-h))/(2*h)
#return x,y,dy_dx
if len(xv)==len(yv)>1: # assure same length
if len(xv)==len(yv)==2: # only two points
if numpy.min(xv)<x<numpy.max(xv):
dx = xv[1] - xv[0]
dy = yv[1] - yv[0]
dydx = dy/dx
y = yv[0] + dydx * (x-xv[0])
return x,y,dydx
else:
raise ValueError("Your coordinate has to be between the input values.")
else:
limit = 1e-10 # avoid hitting a point directly
diff = numpy.array(xv)-x # get differences
index = numpy.argmin(diff*diff) # nearest neighbour
if (xv[index]<x<xv[index+1] # nearest below, positive inclination
or xv[index]>x>xv[index+1]): # nearest above, negative inclination
if diff[index]<limit:
index = [index-1,index+1]
else:
index = [index, index+1]
elif (xv[index-1]<x<xv[index] # nearest above, positive inclination
or xv[index-1]>x>xv[index]): # nearest below, negative inclination
if diff[index]<limit:
index = [index-1,index+1]
else:
index = [index-1,index]
xvnew = xv[index]
yvnew = yv[index]
return get_x_y_dydx(xvnew,yvnew,x) # Allow for a single recursion
else:
raise ValueError("You have to provide the same amount of x- and y-pairs with at least two entries each.")
if axis is None:
axis=matplotlib.pyplot.gca()
if fig is None:
fig=matplotlib.pyplot.gcf()
if y is None and x is not None:
trash=0
(xv,yv)=ToPixelCoords(xv,yv,axis,fig)
#x is provided but y isn't
(x,trash)=ToPixelCoords(x,trash,axis,fig)
#Get the rotation angle and y-value
x,y,dy_dx = get_x_y_dydx(xv,yv,x)
rot = numpy.arctan(dy_dx)/numpy.pi*180.
elif x is None and y is not None:
#y is provided, but x isn't
_xv = xv[::-1]
_yv = yv[::-1]
#Find x by interpolation
x = interp1d(yv, xv)(y)
trash=0
(xv,yv)=ToPixelCoords(xv,yv,axis,fig)
(x,trash)=ToPixelCoords(x,trash,axis,fig)
#Get the rotation angle and y-value
x,y,dy_dx = get_x_y_dydx(xv,yv,x)
rot = numpy.arctan(dy_dx)/numpy.pi*180.
(x,y)=ToDataCoords(x,y,axis,fig)
return (x,y,rot)
def drawLines(Ref,lines,axis,plt_kwargs=None):
"""
Just an internal method to systematically plot values from
the generated 'line' dicts, method is able to cover the whole
saturation curve. Closes the gap at the critical point and
adds a marker between the two last points of bubble and
dew line if they reach up to critical point.
Returns an array of line objects that can be used to change
the colour or style afterwards.
"""
if not plt_kwargs is None:
for line in lines:
line['opts'] = plt_kwargs
plottedLines = []
if len(lines)==2 and (
'q' in str(lines[0]['type']).lower() and 'q' in str(lines[1]['type']).lower()
) and (
( 0 == lines[0]['value'] and 1 == lines[1]['type'] ) or ( 1 == lines[0]['value'] and 0 == lines[1]['type'] ) ):
# We plot the saturation curve
bubble = lines[0]
dew = lines[1]
line, = axis.plot(bubble['x'],bubble['y'],**bubble['opts'])
plottedLines.extend([line])
line, = axis.plot(dew['x'], dew['y'], **dew['opts'])
plottedLines.extend([line])
# Do we need to test if this is T or p?
Tmax = min(bubble['kmax'],dew['kmax'])
if Tmax>CP.PropsSI(Ref,'Tcrit')-2e-5:
axis.plot(numpy.r_[bubble['x'][-1],dew['x'][-1]],numpy.r_[bubble['y'][-1],dew['y'][-1]],**bubble['opts'])
#axis.plot((bubble['x'][-1]+dew['x'][-1])/2.,(bubble['y'][-1]+dew['y'][-1])/2.,'o',color='Tomato')
else:
for line in lines:
line, = axis.plot(line['x'],line['y'],**line['opts'])
plottedLines.extend([line])
return plottedLines
class IsoLines(BasePlot):
def __init__(self, fluid_ref, graph_type, iso_type, unit_system='SI', **kwargs):
BasePlot.__init__(self, fluid_ref, graph_type, unit_system=unit_system,**kwargs)
if not isinstance(iso_type, str):
raise TypeError("Invalid iso_type input, expected a string")
iso_type = iso_type.upper()
if iso_type not in self.COLOR_MAP.keys() and iso_type != 'Q':
raise ValueError('This kind of isoline is not supported for a ' \
+ str(graph_type) + \
' plot. Please choose from '\
+ str(self.COLOR_MAP.keys()) + ' or Q.')
self.iso_type = iso_type
def __set_axis_limits(self, swap_xy):
"""
Generates limits for the axes in terms of x,y defined by 'plot'
based on temperature and pressure.
Returns a tuple containing ((xmin, xmax), (ymin, ymax))
"""
# Get current axis limits, be sure to set those before drawing isolines
# if no limits are set, use triple point and critical conditions
X = [CP.PropsSI(self.graph_type[1],
'T', 1.5*CP.PropsSI(self.fluid_ref, 'Tcrit'),
'P', CP.PropsSI(self.fluid_ref, 'ptriple'),
self.fluid_ref),
CP.PropsSI(self.graph_type[1],
'T', 1.1*CP.PropsSI(self.fluid_ref, 'Tmin'),
'P', 1.5*CP.PropsSI(self.fluid_ref, 'pcrit'),
self.fluid_ref),
CP.PropsSI(self.graph_type[1],
'T', 1.5*CP.PropsSI(self.fluid_ref, 'Tcrit'),
'P', 1.5*CP.PropsSI(self.fluid_ref, 'pcrit'),
self.fluid_ref),
CP.PropsSI(self.graph_type[1],
'T', 1.1*CP.PropsSI(self.fluid_ref, 'Tmin'),
'P', CP.PropsSI(self.fluid_ref, 'ptriple'),
self.fluid_ref)]
Y = [CP.PropsSI(self.graph_type[0],
'T', 1.5*CP.PropsSI(self.fluid_ref, 'Tcrit'),
'P', CP.PropsSI(self.fluid_ref, 'ptriple'),
self.fluid_ref),
CP.PropsSI(self.graph_type[0],
'T', 1.1*CP.PropsSI(self.fluid_ref, 'Tmin') ,
'P', 1.5*CP.PropsSI(self.fluid_ref, 'pcrit'),
self.fluid_ref),
CP.PropsSI(self.graph_type[0],
'T', 1.1*CP.PropsSI(self.fluid_ref, 'Tcrit'),
'P', 1.5*CP.PropsSI(self.fluid_ref, 'pcrit'),
self.fluid_ref),
CP.PropsSI(self.graph_type[0],
'T', 1.5*CP.PropsSI(self.fluid_ref, 'Tmin') ,
'P', CP.PropsSI(self.fluid_ref, 'ptriple'),
self.fluid_ref)]
limits = [[min(X), max(X)], [min(Y), max(Y)]]
if not self.axis.get_autoscalex_on():
limits[0][0] = max([limits[0][0], min(self.axis.get_xlim())])
limits[0][1] = min([limits[0][1], max(self.axis.get_xlim())])
limits[1][0] = max([limits[1][0], min(self.axis.get_ylim())])
limits[1][1] = min([limits[1][1], max(self.axis.get_ylim())])
# Limits correction in case of KSI unit_system
if self.unit_system == 'KSI':
limits[0] = [l*self.KSI_SCALE_FACTOR[self.graph_type[1]] for l in limits[0]]
limits[1] = [l*self.KSI_SCALE_FACTOR[self.graph_type[0]] for l in limits[1]]
self.axis.set_xlim(limits[0])
self.axis.set_ylim(limits[1])
return limits
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 = numpy.unique(numpy.sort(numpy.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 = ( numpy.around(numpy.log10(numpy.abs(inVal))) * -1) + digits + 1
val = numpy.where(val < lim, val, lim)
val = numpy.where(val >-lim, val, -lim)
output = numpy.zeros(inVal.shape)
for i in range(len(inVal)):
output[i] = numpy.around(inVal[i],decimals=int(val[i]))
output = numpy.unique(output)
return output
def get_isolines(self, iso_range=[], num=None, rounding=False):
"""
This is the core method to obtain lines in the dimensions defined
by 'plot' that describe the behaviour of fluid 'Ref'. The constant
value is determined by 'iName' and has the values of 'iValues'.
'iValues' is an array-like object holding at least one element. Lines
are calculated for every entry in 'iValues'. If the input 'num' is
larger than the amount of entries in 'iValues', an internally defined
pattern is used to calculate an appropriate line spacing between the maximum
and minimum values provided in 'iValues'.
Returns lines[num] - an array of dicts containing 'x' and 'y'
coordinates for bubble and dew line. Additionally, the dict holds
the keys 'label' and 'opts', those can be used for plotting as well.
"""
if iso_range is None or (len(iso_range) == 1 and num != 1):
raise ValueError('Automatic interval detection for isoline \
boundaries is not supported yet, use the \
iso_range=[min, max] parameter.')
if len(iso_range) == 2 and num is None:
raise ValueError('Please specify the number of isoline you want \
e.g. num=10')
iso_range = numpy.sort(numpy.unique(iso_range))
def generate_ranges(xmin, xmax, num):
if self.iso_type in ['P', 'D']:
return numpy.logspace(math.log(xmin, 2.),
math.log(xmax, 2.),
num=num,
base=2.)
return numpy.linspace(xmin, xmax, num=num)
# Generate iso ranges
if len(iso_range) == 2:
iso_range = generate_ranges(iso_range[0], iso_range[1], num)
#iso_range = plotRound(iso_range)
#else:
# TODO: Automatic interval detection
# iVal = [CP.PropsSI(iName,'T',T_c[i],'D',rho_c[i],Ref) for i in range(len(T_c))]
# iVal = patterns[iName]([numpy.min(iVal),numpy.max(iVal),num])
if rounding:
iso_range = self.__plotRound(iso_range)
switch_xy_map = {'D': ['TS', 'PH', 'PS'],
'S': ['PH', 'PD', 'PT'],
'T': ['PH', 'PS'],
'H': ['PD']}
#TS: TD is defined, SD is not
#PH: PD is defined, HD is not
#PS: PD is defined, SD is not
#PH: PS is more stable than HS
#PD: PS is defined, DS is not
#PT: PS is defined, TS is not
#PH: PT is defined, HT is not
#PS: PT is defined, ST is not
#PD: PH is defined, DH is not
iso_error_map = {'TD': ['S', 'H'],
'HS': ['T', 'D'],}
switch_xy = False
if self.iso_type in ['D', 'S', 'T', 'H']:
if self.graph_type in switch_xy_map[self.iso_type]:
switch_xy = True
if self.graph_type in ['TD', 'HS']:
if self.iso_type in iso_error_map[self.graph_type]:
raise ValueError('You should not reach this point!')
axis_limits = self.__set_axis_limits(switch_xy)
req_prop = self.graph_type[0]
prop2_name = self.graph_type[1]
if switch_xy:
axis_limits.reverse()
req_prop = self.graph_type[1]
prop2_name = self.graph_type[0]
# Calculate the points
if self.iso_type == 'Q':
lines = self._get_sat_lines(x=iso_range)
return lines
# TODO: Determine saturation state if two phase region present
x_range = numpy.linspace(axis_limits[0][0], axis_limits[0][1], 1000.)
x_mesh = [x_range for i in iso_range]
plot_data = self._get_fluid_data(req_prop,
self.iso_type, iso_range,
prop2_name, x_mesh)
if switch_xy:
plot_data = plot_data[::-1]
lines = []
for j in range(len(plot_data[0])):
line = {
'x': plot_data[0][j],
'y': plot_data[1][j],
# TODO
'label': "", #_getIsoLineLabel(self.iso_type, iso_range[j]),
'type': self.iso_type,
'opts': {'color': self.COLOR_MAP[self.iso_type], 'lw':0.75, 'alpha':0.5 }
}
lines.append(line)
return lines
def draw_isolines(self, iso_range, num=None, rounding=False):
"""
Draw lines with constant values of type 'which' in terms of x and y as
defined by 'plot'. 'iMin' and 'iMax' are minimum and maximum value between
which 'num' get drawn.
There should also be helpful error messages...
"""
if iso_range is None or (len(iso_range) == 1 and num != 1):
raise ValueError('Automatic interval detection for isoline \
boundaries is not supported yet, use the \
iso_range=[min, max] parameter.')
if len(iso_range) == 2 and num is None:
raise ValueError('Please specify the number of isoline you want \
e.g. num=10')
if self.iso_type == 'all':
raise ValueError('Plotting all lines automatically is not \
supported, yet..')
if self.iso_type != 'all':
lines = self.get_isolines(iso_range, num, rounding)
drawn_lines = drawLines(self.fluid_ref, lines, self.axis)
self._plot_default_annotations()
return drawn_lines
#else:
# # TODO: assign limits to values automatically
# ll = _getIsoLineIds(plot)
# if not len(ll)==len(iValues):
# raise ValueError('Please provide a properly sized array of bounds.')
# for c,l in enumerate(ll):
# drawIsoLines(Ref, plot, l, iValues=iValues[c], num=num, axis=axis, fig=fig)
class PropsPlot(BasePlot):
def __init__(self, fluid_name, graph_type, units = 'KSI', reciprocal_density = False, **kwargs):
"""
Create graph for the specified fluid properties
Parameters
----------
fluid_ref : string
The name of the fluid to be plotted
graph_type : string
The graph type to be plotted
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
units : string, ['KSI','SI']
Select the units used for the plotting. 'KSI' is kPa, kJ, K; 'SI' is Pa, J, K
reciprocal_density : bool
If True, 1/rho will be plotted instead of rho
Examples
--------
>>> from CoolProp.Plots import PropsPlot
>>> plt = PropsPlot('Water', 'Ph')
>>> plt.show()
>>> plt = PropsPlot('n-Pentane', 'Ts')
>>> plt.set_axis_limits([-0.5, 1.5, 300, 530])
>>> plt.draw_isolines('Q', [0.1, 0.9])
>>> plt.draw_isolines('P', [100, 2000])
>>> plt.draw_isolines('D', [2, 600])
>>> plt.show()
.. note::
See the online documentation for a list of the available fluids and
graph types
"""
BasePlot.__init__(self, fluid_name, graph_type, unit_system=units, **kwargs)
self.smin = kwargs.get('smin', None)
self.smax = kwargs.get('smax', None)
self._draw_graph()
def __draw_region_lines(self):
lines = self._get_sat_lines(kind='T',
smin=self.smin,
smax=self.smax)
drawLines(self.fluid_ref, lines, self.axis)
def _draw_graph(self):
self.__draw_region_lines()
self._plot_default_annotations()
def draw_isolines(self, iso_type, iso_range, num=10, rounding=False):
iso_lines = IsoLines(self.fluid_ref,
self.graph_type,
iso_type, unit_system = self.unit_system,
axis=self.axis)
iso_lines.draw_isolines(iso_range, num, rounding)
def Ts(Ref, Tmin=None, Tmax=None, show=False, axis=None, *args, **kwargs):
"""
Deprecated. Use :py:func:`CoolProps.Plots.PropsPlot`
"""
plt = PropsPlot(Ref, 'Ts', smin=Tmin, smax=Tmax, axis=axis, *args, **kwargs)
plt._draw_graph()
if show:
plt.show()
else:
plt._draw_graph()
return plt.axis
def Ph(Ref, Tmin=None, Tmax=None, show=False, axis=None, *args, **kwargs):
"""
Deprecated. Use :py:func:`CoolProps.Plots.PropsPlot`
"""
plt = PropsPlot(Ref, 'Ph', smin=Tmin, smax=Tmax, axis=axis, *args, **kwargs)
if show:
plt.show()
else:
plt._draw_graph()
return plt.axis
def Ps(Ref, Tmin=None, Tmax=None, show=False, axis=None, *args, **kwargs):
"""
Deprecated. Use :py:func:`CoolProps.Plots.PropsPlot`
"""
plt = PropsPlot(Ref, 'Ps', smin=Tmin, smax=Tmax, axis=axis, *args, **kwargs)
if show:
plt.show()
else:
plt._draw_graph()
return plt.axis
def PT(Ref, Tmin=None, Tmax=None, show=False, axis=None, *args, **kwargs):
"""
Deprecated. Use :py:func:`CoolProps.Plots.PropsPlot`
"""
plt = PropsPlot(Ref, 'PT', smin=Tmin, smax=Tmax, axis=axis, *args, **kwargs)
if show:
plt.show()
else:
plt._draw_graph()
return plt.axis
def Prho(Ref, Tmin=None, Tmax=None, show=False, axis=None, *args, **kwargs):
"""
"""
plt = PropsPlot(Ref, 'PD', smin=Tmin, smax=Tmax, axis=axis, *args, **kwargs)
if show:
plt.show()
else:
plt._draw_graph()
return plt.axis
def Trho(Ref, Tmin=None, Tmax=None, show=False, axis=None, *args, **kwargs):
"""
Deprecated. Use :py:func:`CoolProps.Plots.PropsPlot`
"""
plt = PropsPlot(Ref, 'TD', smin=Tmin, smax=Tmax, axis=axis, *args, **kwargs)
if show:
plt.show()
else:
plt._draw_graph()
return plt.axis
def hs(Ref, Tmin=None, Tmax=None, show=False, axis=None, *args, **kwargs):
"""
Deprecated. Use :py:func:`CoolProps.Plots.PropsPlot`
"""
plt = PropsPlot(Ref, 'hs', smin=Tmin, smax=Tmax, axis=axis, *args, **kwargs)
if show:
plt.show()
else:
plt._draw_graph()
return plt.axis
def drawIsoLines(Ref, plot, which, iValues=[], num=0, show=False, axis=None):
"""
Draw lines with constant values of type 'which' in terms of x and y as
defined by 'plot'. 'iMin' and 'iMax' are minimum and maximum value
between which 'num' get drawn.
:Note:
:func:`CoolProps.Plots.drawIsoLines` will be depreciated in future
releases and replaced with :func:`CoolProps.Plots.IsoLines`
Parameters
----------
Ref : str
The given reference fluid
plot : str
The plot type used
which : str
The iso line type
iValues : list
The list of constant iso line values
num : int, Optional
The number of iso lines
(Default: 0 - Use iValues list only)
show : bool, Optional
Show the current plot
(Default: False)
axis : :func:`matplotlib.pyplot.gca()`, Optional
The current axis system to be plotted to.
(Default: create a new axis system)
Examples
--------
>>> from matplotlib import pyplot
>>> from CoolProp.Plots import Ts, drawIsoLines
>>>
>>> Ref = 'n-Pentane'
>>> ax = Ts(Ref)
>>> ax.set_xlim([-0.5, 1.5])
>>> ax.set_ylim([300, 530])
>>> quality = drawIsoLines(Ref, 'Ts', 'Q', [0.3, 0.5, 0.7, 0.8], axis=ax)
>>> isobars = drawIsoLines(Ref, 'Ts', 'P', [100, 2000], num=5, axis=ax)
>>> isochores = drawIsoLines(Ref, 'Ts', 'D', [2, 600], num=7, axis=ax)
>>> pyplot.show()
"""
isolines = IsoLines(Ref, plot, which, axis=axis)
lines = isolines.draw_isolines(iValues, num)
if show:
isolines.show()
return lines
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment