diff --git a/dev/TTSE/check_TTSE.py b/dev/TTSE/check_TTSE_old.py similarity index 100% rename from dev/TTSE/check_TTSE.py rename to dev/TTSE/check_TTSE_old.py diff --git a/dev/TTSE/check_TTSE_v4.py b/dev/TTSE/check_TTSE_v4.py new file mode 100644 index 00000000..980c506a --- /dev/null +++ b/dev/TTSE/check_TTSE_v4.py @@ -0,0 +1,322 @@ +import CoolProp.CoolProp as CP +import matplotlib +matplotlib.rc('font', family='serif', serif='Times New Roman') +#from matplotlib2tikz import save as tikz_save + +import matplotlib.pyplot as plt +import matplotlib.colors as colors +import matplotlib.cm as cmx +import matplotlib.ticker +from matplotlib.patches import Ellipse +from matplotlib.transforms import ScaledTranslation +import numpy as np +import random +from numpy import linspace, meshgrid +from matplotlib.mlab import griddata +from matplotlib.gridspec import GridSpec + +# Create the colourmap +#import numpy as np +#import matplotlib.pyplot as plt +import matplotlib._cm, matplotlib.cm +specs = matplotlib._cm.cubehelix(gamma=1.4,s=0.4,r=-0.8,h=2.0) +specs_r = matplotlib.cm._reverse_cmap_spec(specs) +matplotlib.cm.register_cmap(name="jorrithelix" , data=specs) +matplotlib.cm.register_cmap(name="jorrithelix"+"_r", data=specs_r) + + + +def makeGrid(x, y, z, resX=200, resY=200): + "Convert 3 column data to matplotlib grid" + xi = linspace(min(x), max(x), resX) + yi = linspace(min(y), max(y), resY) + Z = griddata(x, y, z, xi, yi) + X, Y = meshgrid(xi, yi) + return X, Y, Z + +def getErrors(p, h, out='D', Ref=''): + "Get the relative errors from table-based interpolation" + errorTTSE = 1e3 + errorBICUBIC = 1e3 + try: + # Using the EOS + CP.disable_TTSE_LUT(Ref) + EOS = CP.PropsSI(out,'P',p,'H',h,Ref) + # Using the TTSE method + CP.enable_TTSE_LUT(Ref) + CP.set_TTSE_mode(Ref,"TTSE") + TTSE = CP.PropsSI(out,'P',p,'H',h,Ref) + # Using the Bicubic method + CP.enable_TTSE_LUT(Ref) + CP.set_TTSE_mode(Ref,"BICUBIC") + BICUBIC = CP.PropsSI(out,'P',p,'H',h,Ref) + errorTTSE = abs(TTSE /EOS-1.0)*100.0 + errorBICUBIC = abs(BICUBIC/EOS-1.0)*100.0 + except ValueError as VE: + print VE + pass + + return errorTTSE,errorBICUBIC + + + +#['YlOrRd', 'PuBuGn', 'hot', 'cubehelix', 'gnuplot', 'gnuplot2']: +for colourmap in ['jorrithelix']: + + for out in ['D']: + ## landscape figure + #fig = plt.figure(figsize=(10,5)) + #ax1 = fig.add_axes((0.08,0.1,0.32,0.83)) + #ax2 = fig.add_axes((0.50,0.1,0.32,0.83)) + #cbar_ax = fig.add_axes([0.80, 0.075, 0.05, 0.875]) + + # portrait figure + #fig = plt.figure(figsize=(5,8)) + #ax1 = plt.subplot2grid((2,8), (0,0), colspan=7) + #ax2 = plt.subplot2grid((2,8), (1,0), colspan=7) + #cbar_ax = plt.subplot2grid((2,8), (0,7), colspan=1, rowspan=2) + + #fig = plt.figure(figsize=(8,4)) + #ax1 = plt.subplot2grid((1,7), (0,0), colspan=3) + #ax2 = plt.subplot2grid((1,7), (0,3), colspan=3) + #cbar_ax = plt.subplot2grid((1,7), (0,6), colspan=1, rowspan=1) + #plt.tight_layout() + fig = plt.figure(figsize=(8,4)) + ax1 = fig.add_subplot(121) + ax2 = fig.add_subplot(122) + #cbar_ax = plt.subplot2grid((1,7), (0,6), colspan=1, rowspan=1) + #plt.tight_layout() + + #Ref = 'R245fa' + #Ref = 'Isopentane' + Ref = 'Air' + + T = np.linspace(CP.PropsSI(Ref,'Tmin')+0.1,CP.PropsSI(Ref,'Tcrit')-0.01,300) + pV = CP.PropsSI('P','T',T,'Q',1,Ref) + hL = CP.PropsSI('H','T',T,'Q',0,Ref) + hV = CP.PropsSI('H','T',T,'Q',1,Ref) + hTP= np.append(hL,[hV[::-1]]) + pTP= np.append(pV,[pV[::-1]]) + + HHH1, PPP1, EEE1 = [], [], [] + HHH2, PPP2, EEE2 = [], [], [] + + cNorm = colors.LogNorm(vmin=1e-10, vmax=1e-1) + scalarMap = cmx.ScalarMappable(norm = cNorm, cmap = plt.get_cmap(colourmap)) + + # Setting the limits for enthalpy and pressure + p_min = CP.PropsSI(Ref,'ptriple') + p_max = 60e5 + h_min = CP.PropsSI('H','T',CP.PropsSI(Ref,'Ttriple')+0.5,'Q',0,Ref) + h_max = CP.PropsSI('H','T',500+273.15,'P',p_max,Ref) + + # Creating some isotherms for better illustration of the cycle + isoT = np.array([0,100,200,300,400])+273.15 + isoP = np.logspace(np.log10(p_min),np.log10(p_max),base=10) + ones = np.ones(isoP.shape) + isoH = [ CP.PropsSI('H','T',T*ones,'P',isoP,Ref) for T in isoT ] + + + print "Lower left and upper right coordinates: ({0},{1}), ({2},{3})".format(h_min,p_min,h_max,p_max) + + CP.set_TTSESinglePhase_LUT_range(Ref,h_min,h_max*1.05,p_min,p_max*1.05) + + for a_useless_counter in range(40000): + + h = random.uniform(h_min,h_max) + p = 10**random.uniform(np.log10(p_min),np.log10(p_max)) + + try: + # Using the EOS + CP.disable_TTSE_LUT(Ref) + rhoEOS = CP.PropsSI('D','P',p,'H',h,Ref) + TEOS = CP.PropsSI('T','P',p,'H',h,Ref) + if out =='C': cpEOS = CP.PropsSI('C','P',p,'H',h,Ref) + + # Using the TTSE method + CP.enable_TTSE_LUT(Ref) + CP.set_TTSE_mode(Ref,"TTSE") + rhoTTSE = CP.PropsSI('D','P',p,'H',h,Ref) + TTTSE = CP.PropsSI('T','P',p,'H',h,Ref) + if out =='C': cpTTSE = CP.PropsSI('C','P',p,'H',h,Ref) + + # Using the Bicubic method + CP.enable_TTSE_LUT(Ref) + CP.set_TTSE_mode(Ref,"BICUBIC") + rhoBICUBIC = CP.PropsSI('D','P',p,'H',h,Ref) + TBICUBIC = CP.PropsSI('T','P',p,'H',h,Ref) + if out =='C': cpBICUBIC = CP.PropsSI('C','P',p,'H',h,Ref) + + if out == 'D': + errorTTSE = abs(rhoTTSE/rhoEOS-1)*100 + errorBICUBIC = abs(rhoBICUBIC/rhoEOS-1)*100 + elif out == 'T': + errorTTSE = abs(TTTSE/TEOS-1)*100 + errorBICUBIC = abs(TBICUBIC/TEOS-1)*100 + elif out == 'C': + errorTTSE = abs(cpTTSE/cpEOS-1)*100 + errorBICUBIC = abs(cpBICUBIC/cpEOS-1)*100 + + HHH1.append(h) + PPP1.append(p) + EEE1.append(errorTTSE) + + HHH2.append(h) + PPP2.append(p) + EEE2.append(errorBICUBIC) + + except ValueError as VE: + #print VE + pass + + HHH1 = np.array(HHH1) + PPP1 = np.array(PPP1) + SC1 = ax1.scatter(HHH1/1e3, PPP1/1e5, s=8, c=EEE1, edgecolors = 'none', cmap = plt.get_cmap(colourmap), norm = cNorm, rasterized=True) + + #X, Y, Z = makeGrid(HHH1, np.log10(PPP1), EEE1) + #SC1 = matplotlib.pyplot.contourf(X, Y, Z, + # alpha=0.75, + # norm=cNorm, + # cmap=matplotlib.pyplot.get_cmap(colourmap))#, + # #rasterized=True) + HHH2 = np.array(HHH2) + PPP2 = np.array(PPP2) + SC2 = ax2.scatter(HHH2/1e3, PPP2/1e5, s=8, c=EEE2, edgecolors = 'none', cmap = plt.get_cmap(colourmap), norm = cNorm, rasterized=True) + + + if out == 'D': + ax1.set_title('rel. density error, TTSE') + ax2.set_title('rel. density error, bicubic') + elif out == 'T': + ax1.set_title('rel. temperature error, TTSE') + ax2.set_title('rel. temperature error, bicubic') + elif out == 'C': + ax1.set_title('rel. heat capacity error, TTSE') + ax2.set_title('rel. heat capacity error, bicubic') + + for ax in [ax1, ax2]: + #h_min = np.ceil(h_min) + delta = 0.1 + delta_min = 1.0+delta + delta_max = 1.0-delta + + #ax.set_xlim(delta_min*h_min/1e3, delta_max*h_max/1e3) + #ax.set_ylim(delta_min*p_min/1e5, delta_max*p_max/1e5) + ax.set_xlim(-155, 800) + ax.set_ylim(0.025, 58) + + ax.set_yscale('log') + + #ticks = np.array([0.02,0.05,0.1,0.2,0.5,1,2,5,10,20,50]) + ticks = np.array([0.05,0.1,0.2,0.5,1,2,5,10,20,50]) + labels = [str(tick) for tick in ticks] + ax.set_yticks(ticks) + ax.set_yticklabels(labels) + ax.get_yaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter()) + + #ticks = [150,250,350,450,550] + #labels = [str(tick) for tick in ticks] + #ax.set_xticks(ticks) + #ax.set_xticklabels(labels) + + #ax.tick_params(axis='y',which='minor',left='off') + + #ax.set_xlabel('Enthalpy [kJ \cdot kg^{-1}]') + ax.set_xlabel('Specific Enthalpy [kJ$\cdot$kg$\mathdefault{^{-1}\!}$]') + ax.set_ylabel('Pressure [bar]') + + #ax.plot(hL/1e3,pV/1e5,'k',lw = 4) + #ax.plot(hV/1e3,pV/1e5,'k',lw = 4) + + ax.plot(hTP/1e3,pTP/1e5,'k',lw = 3) + + for i,T in enumerate(isoT): + ax.plot(isoH[i]/1e3,isoP/1e5,'k',lw = 1) + + + #CB = fig.colorbar(SC1) + #cbar_ax = fig.add_axes([0.80, 0.075, 0.05, 0.875]) + #CB = fig.colorbar(SC1, cax=cbar_ax) + + #CB = matplotlib.pyplot.colorbar(SC2) + #CB.solids.set_rasterized(True) + #ax2.yaxis.set_visible(False) + + #[x0,y0,width,height] + #cbar_ax = fig.add_axes([0.95, 0.00, 0.05, 1.00]) + #CB = fig.colorbar(SC2, ax=[ax1,ax2], cax=cbar_ax) + #CB.solids.set_rasterized(True) + + #from mpl_toolkits.axes_grid1 import make_axes_locatable + #divider = make_axes_locatable(ax2) + #cbar_ax = divider.append_axes("right", "5%", pad="0%") + #CB = plt.colorbar(SC2, cax=cbar_ax) + #CB.solids.set_rasterized(True) + + #CB = fig.colorbar(SC2) + #CB.solids.set_rasterized(True) + + from mpl_toolkits.axes_grid1 import make_axes_locatable + divider = make_axes_locatable(ax2) + ax_cb = divider.new_horizontal(size="5%", pad=0.05) + #fig1 = ax.get_figure() + fig.add_axes(ax_cb) + CB = fig.colorbar(SC2, cax=ax_cb) + + #aspect = 5./2. + #ax1.set_aspect(aspect) + #ax2.set_aspect(aspect) + + CB.solids.set_rasterized(True) + + if out == 'D': + CB.set_label(r'$\|\rho/\rho\mathdefault{_{EOS}-1\|\times 100}$ [%]') + elif out == 'T': + CB.set_label(r'$\|T/T\mathdefault{_{EOS}-1\|\times 100}$ [%]') + elif out == 'C': + CB.set_label(r'$\|c\mathdefault{_p}/c\mathdefault{_{p,EOS}-1\|\times 100}$ [%]') + + + # The plot is finished, now we add an ellipse + #circle=plt.Circle((5,5),.5,color='b',fill=False) + #A scale-free ellipse. + #xy - center of ellipse + #width - total length (diameter) of horizontal axis + #height - total length (diameter) of vertical axis + #angle - rotation in degrees (anti-clockwise) + p_op_min = 1e5 + p_op_max = 3e5 + h_op_min = CP.PropsSI('H','T',400+273.15,'P',p_op_max,Ref) + h_op_max = CP.PropsSI('H','T', 25+273.15,'P',p_op_max,Ref) + + p_op_cen = (p_op_min + p_op_max) / 2.0 + h_op_cen = (h_op_min + h_op_max) / 2.0 + + p_op_hei = p_op_max - p_op_min + h_op_wid = h_op_max - h_op_min + + #for ax in [ax1, ax2]: + ##x,y = 10,0 + ### use the axis scale tform to figure out how far to translate + ##circ_offset = ScaledTranslation(x,y,ax.transScale) + ### construct the composite tform + ##circ_tform = circ_offset + ax.transLimits + ax.transAxes + #ellipse = Ellipse(xy=(h_op_cen,p_op_cen), width=h_op_wid, height=p_op_hei, angle=15, color='black')#, transform=circ_tform) + #ax.add_artist(ellipse) + +# font_def = font_manager.FontProperties(family='Helvetica', style='normal', +# size=sizeOfFont, weight='normal', stretch='normal') +# +# for a in fig.axes: +# for label in [a.get_xticklabels(), a.get_yticklabels()]: +# label.set_fontproperties(ticks_font + + + #plt.savefig(out+'_'+colourmap+'_TTSE_BICUBIC.png', dpi = 300, transparent = True) + #plt.savefig(out+'_'+colourmap+'_TTSE_BICUBIC.eps') + # plt.savefig(out+'_'+colourmap+'_TTSE_BICUBIC.pdf') + plt.tight_layout() + plt.savefig('check_TTSE_'+colourmap+'.pdf' ) + #tikz_save( 'check_TTSE.tikz') + #plt.savefig(out+'_'+colourmap+'_TTSE_BICUBIC.jpg', dpi = 1200) + plt.close() \ No newline at end of file diff --git a/dev/TTSE/check_TTSE_v5.py b/dev/TTSE/check_TTSE_v5.py new file mode 100644 index 00000000..a6eb7881 --- /dev/null +++ b/dev/TTSE/check_TTSE_v5.py @@ -0,0 +1,276 @@ +import CoolProp.CoolProp as CP +import matplotlib +matplotlib.rc('font', family='serif', serif='Times New Roman') +#from matplotlib2tikz import save as tikz_save + +import matplotlib.pyplot as plt +import matplotlib.colors as colors +import matplotlib.cm as cmx +import matplotlib.ticker +from matplotlib.patches import Ellipse +from matplotlib.transforms import ScaledTranslation +import numpy as np +import random +from numpy import linspace, meshgrid +from matplotlib.mlab import griddata +from matplotlib.gridspec import GridSpec + + +def makeGrid(x, y, z, resX=200, resY=200): + "Convert 3 column data to matplotlib grid" + xi = linspace(min(x), max(x), resX) + yi = linspace(min(y), max(y), resY) + Z = griddata(x, y, z, xi, yi) + X, Y = meshgrid(xi, yi) + return X, Y, Z + +def getErrors(p, h, out='D', Ref=''): + "Get the relative errors from table-based interpolation" + errorTTSE = 1e3 + errorBICUBIC = 1e3 + try: + # Using the EOS + CP.disable_TTSE_LUT(Ref) + EOS = CP.PropsSI(out,'P',p,'H',h,Ref) + # Using the TTSE method + CP.enable_TTSE_LUT(Ref) + CP.set_TTSE_mode(Ref,"TTSE") + TTSE = CP.PropsSI(out,'P',p,'H',h,Ref) + # Using the Bicubic method + CP.enable_TTSE_LUT(Ref) + CP.set_TTSE_mode(Ref,"BICUBIC") + BICUBIC = CP.PropsSI(out,'P',p,'H',h,Ref) + errorTTSE = abs(TTSE /EOS-1.0)*100.0 + errorBICUBIC = abs(BICUBIC/EOS-1.0)*100.0 + except ValueError as VE: + print VE + pass + + return errorTTSE,errorBICUBIC + + + +#['YlOrRd', 'PuBuGn', 'hot', 'cubehelix', 'gnuplot', 'gnuplot2']: +for colourmap in ['cubehelix']: + + for out in ['D']: + ## landscape figure + #fig = plt.figure(figsize=(10,5)) + #ax1 = fig.add_axes((0.08,0.1,0.32,0.83)) + #ax2 = fig.add_axes((0.50,0.1,0.32,0.83)) + #cbar_ax = fig.add_axes([0.80, 0.075, 0.05, 0.875]) + + # portrait figure + fig = plt.figure(figsize=(5,8)) + #ax1 = fig.add_axes((0.175,0.575,0.575,0.375)) + #ax2 = fig.add_axes((0.175,0.075,0.575,0.375)) + #cbar_ax = fig.add_axes([0.80, 0.075, 0.05, 0.875]) + #ax1 = fig.add_subplot(241,colspan=3) + #ax2 = fig.add_subplot(245,colspan=3) + #cbar_ax = fig.add_subplot(244,rowspan=2) + + ax1 = plt.subplot2grid((2,8), (0,0), colspan=7) + ax2 = plt.subplot2grid((2,8), (1,0), colspan=7) + cbar_ax = plt.subplot2grid((2,8), (0,7), colspan=1, rowspan=2) + + #Ref = 'R245fa' + #Ref = 'Isopentane' + Ref = 'Air' + + T = np.linspace(CP.PropsSI(Ref,'Tmin')+0.1,CP.PropsSI(Ref,'Tcrit')-0.01,300) + pV = CP.PropsSI('P','T',T,'Q',1,Ref) + hL = CP.PropsSI('H','T',T,'Q',0,Ref) + hV = CP.PropsSI('H','T',T,'Q',1,Ref) + hTP= np.append(hL,[hV[::-1]]) + pTP= np.append(pV,[pV[::-1]]) + + HHH1, PPP1, EEE1 = [], [], [] + HHH2, PPP2, EEE2 = [], [], [] + + cNorm = colors.LogNorm(vmin=1e-10, vmax=1e-1) + scalarMap = cmx.ScalarMappable(norm = cNorm, cmap = plt.get_cmap(colourmap)) + + # Setting the limits for enthalpy and pressure + p_min = CP.PropsSI(Ref,'ptriple') + p_max = 60e5 + h_min = CP.PropsSI('H','T',CP.PropsSI(Ref,'Ttriple')+0.5,'Q',0,Ref) + h_max = CP.PropsSI('H','T',500+273.15,'P',p_max,Ref) + + # Creating some isotherms for better illustration of the cycle + isoT = np.array([0,100,200,300,400])+273.15 + isoP = np.logspace(np.log10(p_min),np.log10(p_max),base=10) + ones = np.ones(isoP.shape) + isoH = [ CP.PropsSI('H','T',T*ones,'P',isoP,Ref) for T in isoT ] + + + print "Lower left and upper right coordinates: ({0},{1}), ({2},{3})".format(h_min,p_min,h_max,p_max) + + CP.set_TTSESinglePhase_LUT_range(Ref,h_min,h_max*1.05,p_min,p_max*1.05) + + for a_useless_counter in range(40000): + + h = random.uniform(h_min,h_max) + p = 10**random.uniform(np.log10(p_min),np.log10(p_max)) + + try: + # Using the EOS + CP.disable_TTSE_LUT(Ref) + rhoEOS = CP.PropsSI('D','P',p,'H',h,Ref) + TEOS = CP.PropsSI('T','P',p,'H',h,Ref) + if out =='C': cpEOS = CP.PropsSI('C','P',p,'H',h,Ref) + + # Using the TTSE method + CP.enable_TTSE_LUT(Ref) + CP.set_TTSE_mode(Ref,"TTSE") + rhoTTSE = CP.PropsSI('D','P',p,'H',h,Ref) + TTTSE = CP.PropsSI('T','P',p,'H',h,Ref) + if out =='C': cpTTSE = CP.PropsSI('C','P',p,'H',h,Ref) + + # Using the Bicubic method + CP.enable_TTSE_LUT(Ref) + CP.set_TTSE_mode(Ref,"BICUBIC") + rhoBICUBIC = CP.PropsSI('D','P',p,'H',h,Ref) + TBICUBIC = CP.PropsSI('T','P',p,'H',h,Ref) + if out =='C': cpBICUBIC = CP.PropsSI('C','P',p,'H',h,Ref) + + if out == 'D': + errorTTSE = abs(rhoTTSE/rhoEOS-1)*100 + errorBICUBIC = abs(rhoBICUBIC/rhoEOS-1)*100 + elif out == 'T': + errorTTSE = abs(TTTSE/TEOS-1)*100 + errorBICUBIC = abs(TBICUBIC/TEOS-1)*100 + elif out == 'C': + errorTTSE = abs(cpTTSE/cpEOS-1)*100 + errorBICUBIC = abs(cpBICUBIC/cpEOS-1)*100 + + HHH1.append(h) + PPP1.append(p) + EEE1.append(errorTTSE) + + HHH2.append(h) + PPP2.append(p) + EEE2.append(errorBICUBIC) + + except ValueError as VE: + #print VE + pass + + HHH1 = np.array(HHH1) + PPP1 = np.array(PPP1) + SC1 = ax1.scatter(HHH1/1e3, PPP1/1e5, s=8, c=EEE1, edgecolors = 'none', cmap = plt.get_cmap(colourmap), norm = cNorm, rasterized=True) + + #X, Y, Z = makeGrid(HHH1, np.log10(PPP1), EEE1) + #SC1 = matplotlib.pyplot.contourf(X, Y, Z, + # alpha=0.75, + # norm=cNorm, + # cmap=matplotlib.pyplot.get_cmap(colourmap))#, + # #rasterized=True) + HHH2 = np.array(HHH2) + PPP2 = np.array(PPP2) + SC2 = ax2.scatter(HHH2/1e3, PPP2/1e5, s=8, c=EEE2, edgecolors = 'none', cmap = plt.get_cmap(colourmap), norm = cNorm, rasterized=True) + + + if out == 'D': + ax1.set_title('rel. density error, TTSE') + ax2.set_title('rel. density error, bicubic') + elif out == 'T': + ax1.set_title('rel. temperature error, TTSE') + ax2.set_title('rel. temperature error, bicubic') + elif out == 'C': + ax1.set_title('rel. heat capacity error, TTSE') + ax2.set_title('rel. heat capacity error, bicubic') + + for ax in [ax1, ax2]: + #h_min = np.ceil(h_min) + delta = 0.1 + delta_min = 1.0+delta + delta_max = 1.0-delta + + #ax.set_xlim(delta_min*h_min/1e3, delta_max*h_max/1e3) + #ax.set_ylim(delta_min*p_min/1e5, delta_max*p_max/1e5) + ax.set_xlim(-155, 800) + ax.set_ylim(0.025, 58) + + ax.set_yscale('log') + + #ticks = np.array([0.02,0.05,0.1,0.2,0.5,1,2,5,10,20,50]) + ticks = np.array([0.05,0.1,0.2,0.5,1,2,5,10,20,50]) + labels = [str(tick) for tick in ticks] + ax.set_yticks(ticks) + ax.set_yticklabels(labels) + ax.get_yaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter()) + + #ticks = [150,250,350,450,550] + #labels = [str(tick) for tick in ticks] + #ax.set_xticks(ticks) + #ax.set_xticklabels(labels) + + #ax.tick_params(axis='y',which='minor',left='off') + + #ax.set_xlabel('Enthalpy [kJ \cdot kg^{-1}]') + ax.set_xlabel('Specific Enthalpy [kJ$\cdot$kg$\mathdefault{^{-1}\!}$]') + ax.set_ylabel('Pressure [bar]') + + #ax.plot(hL/1e3,pV/1e5,'k',lw = 4) + #ax.plot(hV/1e3,pV/1e5,'k',lw = 4) + + ax.plot(hTP/1e3,pTP/1e5,'k',lw = 3) + + for i,T in enumerate(isoT): + ax.plot(isoH[i]/1e3,isoP/1e5,'k',lw = 1) + + + #CB = fig.colorbar(SC1) + #cbar_ax = fig.add_axes([0.80, 0.075, 0.05, 0.875]) + CB = fig.colorbar(SC1, cax=cbar_ax) + if out == 'D': + CB.set_label(r'$\|\rho/\rho\mathdefault{_{EOS}-1\|\times 100}$ [%]') + elif out == 'T': + CB.set_label(r'$\|T/T\mathdefault{_{EOS}-1\|\times 100}$ [%]') + elif out == 'C': + CB.set_label(r'$\|c\mathdefault{_p}/c\mathdefault{_{p,EOS}-1\|\times 100}$ [%]') + + + # The plot is finished, now we add an ellipse + #circle=plt.Circle((5,5),.5,color='b',fill=False) + #A scale-free ellipse. + #xy - center of ellipse + #width - total length (diameter) of horizontal axis + #height - total length (diameter) of vertical axis + #angle - rotation in degrees (anti-clockwise) + p_op_min = 1e5 + p_op_max = 3e5 + h_op_min = CP.PropsSI('H','T',400+273.15,'P',p_op_max,Ref) + h_op_max = CP.PropsSI('H','T', 25+273.15,'P',p_op_max,Ref) + + p_op_cen = (p_op_min + p_op_max) / 2.0 + h_op_cen = (h_op_min + h_op_max) / 2.0 + + p_op_hei = p_op_max - p_op_min + h_op_wid = h_op_max - h_op_min + + #for ax in [ax1, ax2]: + ##x,y = 10,0 + ### use the axis scale tform to figure out how far to translate + ##circ_offset = ScaledTranslation(x,y,ax.transScale) + ### construct the composite tform + ##circ_tform = circ_offset + ax.transLimits + ax.transAxes + #ellipse = Ellipse(xy=(h_op_cen,p_op_cen), width=h_op_wid, height=p_op_hei, angle=15, color='black')#, transform=circ_tform) + #ax.add_artist(ellipse) + +# font_def = font_manager.FontProperties(family='Helvetica', style='normal', +# size=sizeOfFont, weight='normal', stretch='normal') +# +# for a in fig.axes: +# for label in [a.get_xticklabels(), a.get_yticklabels()]: +# label.set_fontproperties(ticks_font + + + #plt.savefig(out+'_'+colourmap+'_TTSE_BICUBIC.png', dpi = 300, transparent = True) + #plt.savefig(out+'_'+colourmap+'_TTSE_BICUBIC.eps') + # plt.savefig(out+'_'+colourmap+'_TTSE_BICUBIC.pdf') + plt.tight_layout() + plt.savefig('check_TTSE.pdf' ) + #tikz_save( 'check_TTSE.tikz') + #plt.savefig(out+'_'+colourmap+'_TTSE_BICUBIC.jpg', dpi = 1200) + plt.close() \ No newline at end of file diff --git a/dev/Tickets/351.py b/dev/Tickets/351.py new file mode 100644 index 00000000..b06cc63f --- /dev/null +++ b/dev/Tickets/351.py @@ -0,0 +1,63 @@ +from CoolProp.CoolProp import PropsSI +from CoolProp.Plots import PropertyPlot +import matplotlib.pylab as pl +from numpy import * +import CoolProp + +eta_e_s = 0.88 +eta_p_s = 0.5 +T_max = 550 + 273.15 +p_max = 13000.e3 +p_cd = 5.e3 + +def SimpleRankineCycle(T3, p3, p1, epsilon_e, epsilon_p, fluid='water'): + h1 = PropsSI('H', 'P', p1, 'Q', 0., fluid) + s1 = PropsSI('S', 'P', p1, 'Q', 0., fluid) + T1 = PropsSI('T', 'P', p1, 'Q', 0., fluid) + + p2 = p3 + h2 = h1 + (PropsSI('H', 'P', p2, 'S', s1, fluid) - h1) / epsilon_p + s2 = PropsSI('S', 'H', h2, 'P', p2, fluid) + T2 = PropsSI('T', 'H', h2, 'P', p2, fluid) + + h3 = PropsSI('H', 'P', p3, 'T', T3, fluid) + s3 = PropsSI('S', 'H', h3, 'P', p3, fluid) + + p4 = p1 + h4 = h3 - epsilon_e * (h3 - PropsSI('H', 'P', p4, 'S', s3, fluid)) + s4 = PropsSI('S', 'H', h4, 'P', p4, fluid) + T4 = PropsSI('T', 'H', h4, 'P', p4, fluid) + + w_net = h3 - h4 + q_boiler = h3 - h2 + eta_c = w_net / q_boiler + + Ts = PropertyPlot(fluid, 'Ts', 'KSI') + Ts.set_axis_limits([0., 12., 200., 900.]) + Ts.calc_isolines(CoolProp.iP, [Ts.system.P.from_SI(p1), Ts.system.P.from_SI(p3)], num=10) + Ts.calc_isolines(CoolProp.iQ, [Ts.system.Q.from_SI(0.), Ts.system.Q.from_SI(1.)], num=11) + Ts.draw_isolines() + + states = zip(Ts.system.S.from_SI(array([s1,s2,s3,s4,s1])),Ts.system.T.from_SI(array([T1,T2,T3,T4,T1]))) + Ts.draw_process(states, iso_types=None, line_opts={'color':'red', 'lw':1.5}) + isot = [ + None, # non-isentropic pumping from 1 to 2 + CoolProp.iP, # p2=p3 + None, + CoolProp.iP, # p4=p1 + ] + Ts.draw_process(states, iso_types=isot, line_opts={'color':'green', 'lw':1.5}) + + + ax = Ts.axis + ax.text(Ts.system.S.from_SI(s1), Ts.system.T.from_SI(T1), ' 1', fontsize=10, rotation=0, color='r') + ax.text(Ts.system.S.from_SI(s2), Ts.system.T.from_SI(T2), ' 2', fontsize=10, rotation=0, color='r') + ax.text(Ts.system.S.from_SI(s3), Ts.system.T.from_SI(T3), ' 3', fontsize=10, rotation=0, color='r') + ax.text(Ts.system.S.from_SI(s4), Ts.system.T.from_SI(T4), ' 4', fontsize=10, rotation=0, color='r') + ax.text(Ts.system.S.from_SI(8e3),Ts.system.T.from_SI(850),"Efficiency: %.1f%%" %(eta_c*100.)) + ax.text(Ts.system.S.from_SI(8e3),Ts.system.T.from_SI(800),"Net work: %d kJ/kg" %(w_net/1000)) + ax.text(Ts.system.S.from_SI(8e3),Ts.system.T.from_SI(750),"Heat input: %d kJ/kg" %(q_boiler/1000)) + return Ts + +Ts = SimpleRankineCycle(T_max, p_max, p_cd, eta_e_s, eta_p_s, fluid="water") +Ts.savefig('ticket-351.pdf') \ No newline at end of file diff --git a/wrappers/Python/CoolProp/CoolProp.pyx b/wrappers/Python/CoolProp/CoolProp.pyx index cf2ede9f..4c988306 100644 --- a/wrappers/Python/CoolProp/CoolProp.pyx +++ b/wrappers/Python/CoolProp/CoolProp.pyx @@ -481,6 +481,32 @@ cpdef int get_debug_level(): # return False # + +cpdef extract_backend(string in_str): + """ + A Python wrapper of C++ function :cpapi:`CoolProp::extract_backend` . + """ + cdef string bck, fld + # Extract the backend and the fluid from the input string + _extract_backend(in_str, bck, fld) + return bck, fld + + +cpdef extract_fractions(string flds): + """ + A Python wrapper of C++ function :cpapi:`CoolProp::extract_fractions` . + """ + cdef vector[double] frcs + cdef string del_flds + # Extract the fractions + #frcs.clear() + frcs.push_back(1.0) + del_flds = _extract_fractions(flds, frcs) + # Extract the fluids + fluids = del_flds.split('&') + return fluids,frcs + + cdef toSI(constants_header.parameters key, double val): """ Convert a value in kSI system to SI system (supports a limited subset of variables) diff --git a/wrappers/Python/CoolProp/Plots/Common.py b/wrappers/Python/CoolProp/Plots/Common.py index 1376e0b8..4fb1a2db 100644 --- a/wrappers/Python/CoolProp/Plots/Common.py +++ b/wrappers/Python/CoolProp/Plots/Common.py @@ -2,258 +2,720 @@ from __future__ import print_function, unicode_literals -import matplotlib -import numpy +import matplotlib.pyplot as plt +import numpy as np import CoolProp.CoolProp as CP +from abc import ABCMeta +from CoolProp import AbstractState +from CoolProp.CoolProp import PropsSI,extract_backend,extract_fractions +import CoolProp +import warnings +from scipy.interpolate.interpolate import interp1d +from six import with_metaclass -SMALL = 1E-5 -class BasePlot(object): - #TODO: Simplify / Consolidate dictionary maps - AXIS_LABELS = {'KSI': {'T': ["Temperature", r"[K]"], - 'P': ["Pressure", r"[kPa]"], - 'S': ["Entropy", r"[kJ/kg/K]"], - 'H': ["Enthalpy", r"[kJ/kg]"], - 'U': ["Internal Energy", r"[kJ/kg]"], - 'D': ["Density", r"[kg/m$^3$]"] - }, - 'SI': {'T': ["Temperature", r"[K]"], - 'P': ["Pressure", r"[Pa]"], - 'S': ["Entropy", r"[J/kg/K]"], - 'H': ["Enthalpy", r"[J/kg]"], - 'U': ["Internal Energy", r"[J/kg]"], - 'D': ["Density", r"[kg/m$^3$]"] - } - } - COLOR_MAP = {'T': 'Darkred', - 'P': 'DarkCyan', - 'H': 'DarkGreen', - 'D': 'DarkBlue', - 'S': 'DarkOrange', - 'Q': 'black'} +def process_fluid_state(fluid_ref): + """Check input for state object or fluid string + + Parameters + ---------- + fluid_ref : str, CoolProp.AbstractState + + Returns + ------- + CoolProp.AbstractState + """ + # Process the fluid and set self._state + if isinstance(fluid_ref, basestring): + backend, fluids = extract_backend(fluid_ref) + fluids, fractions = extract_fractions(fluids) + #if backend==u'?': backend = u'HEOS' +# # TODO: Fix the backend extraction etc +# fluid_def = fluid_ref.split('::') +# if len(fluid_def)==2: +# backend = fluid_def[0] +# fluid = fluid_def[1] +# elif len(fluid_def)==1: +# backend = "HEOS" +# fluid = fluid_def[0] +# else: +# raise ValueError("This is not a valid fluid_ref string: {0:s}".format(str(fluid_ref))) + state = AbstractState(backend, '&'.join(fluids)) + #state.set_mass_fractions(fractions) + return state + elif isinstance(fluid_ref, AbstractState): + return fluid_ref + raise TypeError("Invalid fluid_ref input, expected a string or an abstract state instance.") - #: Scale factors to multiply SI units by in order to obtain kSI units - KSI_SCALE_FACTOR = {'T' : 1.0, - 'P' : 0.001, - 'H' : 0.001, - 'U' : 0.001, - 'D' : 1, - 'S' : 0.001, - 'Q' : 1.0} - - SYMBOL_MAP_KSI = {'T' : [r'$T = ', r'$ K'], - 'P' : [r'$p = ', r'$ kPa'], - 'H' : [r'$h = ', r'$ kJ/kg'], - 'U' : [r'$h = ', r'$ kJ/kg'], - 'D' : [r'$\rho = ', r'$ kg/m$^3$'], - 'S' : [r'$s = ', r'$ kJ/kg-K'], - 'Q' : [r'$x = ', r'$']} - - SYMBOL_MAP_SI = {'T' : [r'$T = ', r'$ K'], - 'P' : [r'$p = ', r'$ Pa'], - 'H' : [r'$h = ', r'$ J/kg'], - 'U' : [r'$h = ', r'$ J/kg'], - 'D' : [r'$\rho = ', r'$ kg/m$^3$'], - 'S' : [r'$s = ', r'$ J/kg-K'], - 'Q' : [r'$x = ', r'$']} - LINE_IDS = {'TS': ['P', 'D'], #'H'], - 'PH': ['S', 'T', 'D'], - 'HS': ['P'], #'T', 'D'], - 'PS': ['H', 'T', 'D'], - 'PD': ['T', 'S', 'H'], - 'TD': ['P'], #'S', 'H'], - 'PT': ['D', 'P', 'S'], - 'PU': []} +def _get_index(prop): + if isinstance(prop, basestring): + return CP.get_parameter_index(prop) + elif isinstance(prop, int): + return prop + else: + raise ValueError("Invalid input, expected a string or an int, not {0:s}.".format(str(prop))) - def __init__(self, fluid_ref, graph_type, unit_system = 'KSI', **kwargs): - if not isinstance(graph_type, str): - raise TypeError("Invalid graph_type input, expected a string") +class BaseQuantity(object): + """A very basic property that can convert an input to and from a + given unit system, note that the conversion from SI units starts + with a multiplication. If you need to remove an offset, use the + off_SI property. + Examples with temperature: + celsius = BaseQuantity(add_SI=-273.15) + fahrenheit = BaseQuantity(add_SI=32.0, mul_SI=1.8, off_SI=-273.15) + Examples with pressure: + bar = BaseQuantity(mul_SI=1e-5) + psi = BaseQuantity(mul_SI=0.000145037738) + """ + def __init__(self, add_SI=0.0, mul_SI=1.0, off_SI=0.0): + self._add_SI = add_SI + self._mul_SI = mul_SI + self._off_SI = off_SI - graph_type = graph_type.upper() - if len(graph_type) >= 2 and graph_type[1:len(graph_type)] == 'RHO': - graph_type = graph_type[0] + graph_type[1:len(graph_type)] + @property + def add_SI(self): return self._add_SI + @add_SI.setter + def add_SI(self, value): self._add_SI = value + @property + def mul_SI(self): return self._mul_SI + @mul_SI.setter + def mul_SI(self, value): self._mul_SI = value + @property + def off_SI(self): return self._off_SI + @off_SI.setter + def off_SI(self, value): self._off_SI = value + + def from_SI(self, value): return ((value+self.off_SI)*self.mul_SI)+self.add_SI + def to_SI(self, value): return (value-self.add_SI)/self.mul_SI-self.off_SI + - if graph_type.upper() not in self.LINE_IDS.keys(): - raise ValueError(''.join(["You have to specify the kind of ", - "plot, use one of", - str(self.LINE_IDS.keys())])) +class BaseDimension(BaseQuantity): + """A dimension is a class that extends the BaseQuantity and adds a label, a symbol and a unit label""" + def __init__(self, add_SI=0.0, mul_SI=1.0, off_SI=0.0, label='', symbol='', unit=''): + self._label = label + self._symbol = symbol + self._unit = unit + super(BaseDimension, self).__init__(add_SI=add_SI, mul_SI=mul_SI, off_SI=off_SI) - self.graph_drawn = False - self.fluid_ref = fluid_ref - self.graph_type = graph_type.upper() - self.unit_system = unit_system + @property + def label(self): return self._label + @label.setter + def label(self, value): self._label = value + @property + def symbol(self): return self._symbol + @symbol.setter + def symbol(self, value): self._symbol = value + @property + def unit(self): return self._unit + @unit.setter + def unit(self, value): self._unit = value - self.axis = kwargs.get('axis', None) - if self.axis is None: - self.axis = matplotlib.pyplot.gca() - def __sat_bounds(self, kind, smin=None, smax=None): - """ - Generates limits for the saturation line in either T or p determined - by 'kind'. If xmin or xmax are provided, values will be checked - against the allowable range for the EOS and an error might be - generated. +class PropertyDict(with_metaclass(ABCMeta),object): + """A collection of dimensions for all the required quantities""" + + def __init__(self): + self._D = None + self._H = None + self._P = None + self._S = None + self._T = None + self._U = None + self._Q = None + + @property + def D(self): return self._D + @D.setter + def D(self, value): self._D = value + @property + def H(self): return self._H + @H.setter + def H(self, value): self._H = value + @property + def P(self): return self._P + @P.setter + def P(self, value): self._P = value + @property + def S(self): return self._S + @S.setter + def S(self, value): self._S = value + @property + def T(self): return self._T + @T.setter + def T(self, value): self._T = value + @property + def U(self): return self._U + @U.setter + def U(self, value): self._U = value + @property + def Q(self): return self._Q + @Q.setter + def Q(self, value): self._Q = value + + @property + def dimensions(self): + return { + CoolProp.iDmass : self._D, + CoolProp.iHmass : self._H, + CoolProp.iP : self._P, + CoolProp.iSmass : self._S, + CoolProp.iT : self._T, + CoolProp.iUmass : self._U, + CoolProp.iQ : self._Q + } + + def __getitem__(self, index): + """Allow for property access via square brackets""" + idx = _get_index(index) + if idx == CoolProp.iDmass : return self.D + elif idx == CoolProp.iHmass : return self.H + elif idx == CoolProp.iP : return self.P + elif idx == CoolProp.iSmass : return self.S + elif idx == CoolProp.iT : return self.T + elif idx == CoolProp.iUmass : return self.U + elif idx == CoolProp.iQ : return self.Q + else: raise IndexError("Unknown index \"{0:s}\".".format(str(index))) + + + def __setitem__(self, index, value): + """Allow for property access via square brackets""" + idx = _get_index(index) + if idx == CoolProp.iDmass : self.D = value + elif idx == CoolProp.iHmass : self.H = value + elif idx == CoolProp.iP : self.P = value + elif idx == CoolProp.iSmass : self.S = value + elif idx == CoolProp.iT : self.T = value + elif idx == CoolProp.iUmass : self.U = value + elif idx == CoolProp.iQ : self.Q = value + else: raise IndexError("Unknown index \"{0:s}\".".format(str(index))) + + + - Returns a tuple containing (xmin, xmax) - """ - if kind == 'P': - name = 'pressure' - min_key = 'ptriple' - elif kind == 'T': - name = 'temperature' - min_key = 'Tmin' - fluid_min = CP.PropsSI(self.fluid_ref, min_key) - fluid_crit = CP.PropsSI(self.fluid_ref, ''.join([kind, 'crit'])) +class SIunits(PropertyDict): + def __init__(self): + self._D = BaseDimension(add_SI=0.0, mul_SI=1.0, off_SI=0.0, label='Density', symbol=ur'ρ', unit=ur'kg/m³') + self._H = BaseDimension(add_SI=0.0, mul_SI=1.0, off_SI=0.0, label='Specific Enthalpy', symbol=ur'h', unit=ur'J/kg') + self._P = BaseDimension(add_SI=0.0, mul_SI=1.0, off_SI=0.0, label='Pressure', symbol=ur'p', unit=ur'Pa') + self._S = BaseDimension(add_SI=0.0, mul_SI=1.0, off_SI=0.0, label='Specific Entropy', symbol=ur's', unit=ur'J/kg/K') + self._T = BaseDimension(add_SI=0.0, mul_SI=1.0, off_SI=0.0, label='Temperature', symbol=ur'T', unit=ur'K') + self._U = BaseDimension(add_SI=0.0, mul_SI=1.0, off_SI=0.0, label='Specific Internal Energy', symbol=ur'u', unit=ur'J/kg') + self._Q = BaseDimension(add_SI=0.0, mul_SI=1.0, off_SI=0.0, label='Vapour Quality', symbol=ur'x', unit=ur'') - if smin is None: - smin = fluid_min + SMALL - elif smin > fluid_crit: - raise ValueError(''.join(['Minimum ', name, - ' cannot be greater than fluid critical ', - name, '.'])) +class KSIunits(SIunits): + def __init__(self): + super(KSIunits, self).__init__() + self.H.mul_SI=1e-3 + self.H.unit=r'kJ/kg' + self.P.mul_SI=1e-3 + self.P.unit=r'kPa' + self.S.mul_SI=1e-3 + self.S.unit=r'kJ/kg/K' + self.U.mul_SI=1e-3 + self.U.unit=r'kJ/kg' - if smax is None: - smax = fluid_crit - SMALL - elif smax > fluid_crit: - raise ValueError(''.join(['Maximum ', name, - ' cannot be greater than fluid critical ', - name, '.'])) +class EURunits(KSIunits): + def __init__(self): + super(EURunits, self).__init__() + self.P.mul_SI=1e-5 + self.P.unit=r'bar' + self.T.add_SI=-273.15 + self.T.unit=ur'\u00B0C' - smin = max(smin, fluid_min + SMALL) - smax = min(smax, fluid_crit - SMALL) - return (smin, smax) +class Base2DObject(with_metaclass(ABCMeta),object): + """A container for shared settings and constants for the + isolines and the property plots.""" + + # A list of supported plot + TS = CoolProp.iT*10 + CoolProp.iSmass + PH = CoolProp.iP*10 + CoolProp.iHmass + HS = CoolProp.iHmass*10 + CoolProp.iSmass + PS = CoolProp.iP*10 + CoolProp.iSmass + PD = CoolProp.iP*10 + CoolProp.iDmass + TD = CoolProp.iT*10 + CoolProp.iDmass + PT = CoolProp.iP*10 + CoolProp.iT + PU = CoolProp.iP*10 + CoolProp.iUmass + + PLOTS = { + 'TS': TS, + 'PH': PH, + 'HS': HS, + 'PS': PS, + 'PD': PD, + 'TD': TD, + 'PT': PT, + } + + PLOTS_INV = {v: k for k, v in PLOTS.items()} + +# # A list of supported plot +# @property +# def TS(self): return type(self).TS +# @property +# def PH(self): return CoolProp.iP*10 + CoolProp.iHmass +# @property +# def HS(self): return CoolProp.iHmass*10 + CoolProp.iSmass +# @property +# def PS(self): return CoolProp.iP*10 + CoolProp.iSmass +# @property +# def PD(self): return CoolProp.iP*10 + CoolProp.iDmass +# @property +# def TD(self): return CoolProp.iT*10 + CoolProp.iDmass +# @property +# def PT(self): return CoolProp.iP*10 + CoolProp.iT +# @property +# def PU(self): return CoolProp.iP*10 + CoolProp.iUmass - def _get_fluid_data(self, req_prop, - prop1_name, prop1_vals, - prop2_name, prop2_vals): - """ - Calculates lines for constant iName (iVal) over an interval of xName - (xVal). Returns (x[],y[]) - a tuple of arrays containing the values - in x and y dimensions. - """ - if len(prop1_vals) != len(prop2_vals): - raise ValueError(''.join(['We need the same number of x value ', - 'arrays as iso quantities.'])) + def __init__(self, x_type, y_type, state=None, small=None): + self._x_index = _get_index(x_type) + self._y_index = _get_index(y_type) + if small is not None: self._small = small + else: self._small = 1e-7 + if state is not None: self.state = state + else: self._state = None - y_vals = [] - x_vals = [] + # A list of supported plot + @property + def x_index(self): return self._x_index + @property + def y_index(self): return self._y_index + @property + def state(self): return self._state + @state.setter + def state(self, value): + self._state = process_fluid_state(value) + self._T_small = self._state.trivial_keyed_output(CoolProp.iT_critical)*self._small + self._P_small = self._state.trivial_keyed_output(CoolProp.iP_critical)*self._small - # Calculate the values in SI units - for i, p1_val in enumerate(prop1_vals): - x_vals.append(prop2_vals[i]) - y_vals.append(CP.PropsSI(req_prop, - prop1_name, [p1_val]*len(prop2_vals[i]), # Convert to an iterable the same size as second input array - prop2_name, prop2_vals[i], - self.fluid_ref)) - return numpy.array([x_vals, y_vals]) + + def _get_sat_bounds(self, kind, smin=None, smax=None): + """Generates limits for the saturation line in either T or p determined + by 'kind'. If smin or smax are provided, values will be checked + against the allowable range for the EOS and a warning might be + generated. Returns a tuple containing (xmin, xmax)""" - def _get_sat_lines(self, kind='T', smin=None, - smax=None, num=500, x=[0., 1.]): - """ - Calculates bubble and dew line in the quantities for your plot. - You can specify if you need evenly spaced entries in either - pressure or temperature by supplying kind='p' and kind='T' - (default), respectively. - Limits can be set with kmin (default: minimum from EOS) and - kmax (default: critical value). - Returns lines[] - a 2D array of dicts containing 'x' and 'y' - coordinates for bubble and dew line. Additionally, the dict holds - the keys 'kmax', 'label' and 'opts', those can be used for plotting - as well. - """ - if not kind.upper() in ['T', 'P']: - raise ValueError(''.join(["Invalid input for determining the ", - "saturation lines... Expected either ", - "'T' or 'P'"])) - - smin, smax = self.__sat_bounds(kind, smin=smin, smax=smax) - sat_range = numpy.linspace(smin, smax, num) - sat_mesh = numpy.array([sat_range for i in x]) - - x_vals = sat_mesh - y_vals = sat_mesh - if self.graph_type[1] != kind: - _, x_vals = self._get_fluid_data(self.graph_type[1], - 'Q', x, - kind, sat_mesh) - - if self.graph_type[0] != kind: - _, y_vals = self._get_fluid_data(self.graph_type[0], - 'Q', x, - kind, sat_mesh) - - if self.unit_system == 'KSI': - x_vals *= self.KSI_SCALE_FACTOR[self.graph_type[1]] - y_vals *= self.KSI_SCALE_FACTOR[self.graph_type[0]] - - # Merge the two lines, capital Y holds important information. - # We merge on X values - # Every entry, eg. Xy, contains two arrays of values. - sat_lines = [] - for i in range(len(x_vals)): # two dimensions: i = {0,1} - line = {'x': x_vals[i], - 'y': y_vals[i], - 'smax': smax} - - line['label'] = self.SYMBOL_MAP_KSI['Q'][0] + str(x[i]) - line['type'] = 'Q' - line['value'] = x[i] - line['unit'] = self.SYMBOL_MAP_KSI['Q'][1] - line['opts'] = {'color': self.COLOR_MAP['Q'], - 'lw': 1.0} - - if x[i] == 0.: - line['label'] = 'bubble line' - elif x[i] == 1.: - line['label'] = 'dew line' + # TODO: REFPROP backend does not have ptriple. + T_triple = self._state.trivial_keyed_output(CoolProp.iT_triple) + try: + T_min = self._state.trivial_keyed_output(CoolProp.iT_min) + except: + T_min = T_triple + self._state.update(CoolProp.QT_INPUTS, 0, max([T_triple,T_min])+self._T_small) + kind = _get_index(kind) + if kind == CoolProp.iP: + fluid_min = self._state.keyed_output(CoolProp.iP)+self._P_small + fluid_max = self._state.trivial_keyed_output(CoolProp.iP_critical)-self._P_small + elif kind == CoolProp.iT: + fluid_min = self._state.keyed_output(CoolProp.iT)+self._T_small + fluid_max = self._state.trivial_keyed_output(CoolProp.iT_critical)-self._T_small + else: + raise ValueError("Saturation boundaries have to be defined in T or P, but not in {0:s}".format(str(kind))) + + if smin is not None: + if fluid_min < smin < fluid_max: + sat_min = smin else: - line['opts']['lw'] = 0.75 - line['opts']['alpha'] = 0.5 + warnings.warn( + "Your minimum {0:s} has been ignored, {1:f} is not between {2:f} and {3:f}".format(self.PROPERTIES[kind],smin,fluid_min,fluid_max), + UserWarning) + sat_min = fluid_min + else: + sat_min = fluid_min + + if smax is not None: + if fluid_min < smax < fluid_max: + sat_max = smax + else: + warnings.warn( + "Your maximum {0:s} has been ignored, {1:f} is not between {2:f} and {3:f}".format(self.PROPERTIES[kind],smax,fluid_min,fluid_max), + UserWarning) + sat_max = fluid_max + else: + sat_max = fluid_max - sat_lines.append(line) + return sat_min, sat_max + + + +class IsoLine(Base2DObject): + """An object that holds the functions to calculate a line of + a constant property in the dimensions of a property plot. This + class only uses SI units.""" + + # Normally we calculate a sweep in x-dimensions, but + # sometimes a sweep in y-dimensions is better. + XY_SWITCH = { + CoolProp.iDmass: { Base2DObject.TS:True , Base2DObject.PH:True , Base2DObject.HS:False, Base2DObject.PS:True , Base2DObject.PD:None , Base2DObject.TD:None , Base2DObject.PT:False}, + CoolProp.iHmass: { Base2DObject.TS:False, Base2DObject.PH:None , Base2DObject.HS:None , Base2DObject.PS:True , Base2DObject.PD:True , Base2DObject.TD:False, Base2DObject.PT:False}, + CoolProp.iP : { Base2DObject.TS:False, Base2DObject.PH:None , Base2DObject.HS:False, Base2DObject.PS:None , Base2DObject.PD:None , Base2DObject.TD:False, Base2DObject.PT:None }, + CoolProp.iSmass: { Base2DObject.TS:None , Base2DObject.PH:True , Base2DObject.HS:None , Base2DObject.PS:None , Base2DObject.PD:True , Base2DObject.TD:False, Base2DObject.PT:True }, + CoolProp.iT : { Base2DObject.TS:None , Base2DObject.PH:True , Base2DObject.HS:False, Base2DObject.PS:False, Base2DObject.PD:False, Base2DObject.TD:None , Base2DObject.PT:None }, + CoolProp.iQ : { Base2DObject.TS:True , Base2DObject.PH:True , Base2DObject.HS:True , Base2DObject.PS:True , Base2DObject.PD:True , Base2DObject.TD:True , Base2DObject.PT:False} + } + + # Abort interpolation if there are not enough + # valid entries. + VALID_REQ = 5.0/100.0 + + def __init__(self, i_index, x_index, y_index, value=0.0, state=None): + super(IsoLine, self).__init__(x_index, y_index, state) + self._i_index = _get_index(i_index) + if value is not None: self.value = value + else: self._value = None + self._x = None + self._y = None + + @property + def i_index(self): return self._i_index + @property + def value(self): return self._value + @value.setter + def value(self, value): self._value = float(value) + @property + def x(self): return self._x + @x.setter + def x(self, value): self._x = np.array(value) + @property + def y(self): return self._y + @y.setter + def y(self, value): self._y = np.array(value) + + def get_update_pair(self): + """Processes the values for the isoproperty and the graph dimensions + to figure which should be used as inputs to the state update. Returns + a tuple with the indices for the update call and the property constant. + For an isobar in a Ts-diagram it returns the default order and the + correct constant for the update pair: + get_update_pair(CoolProp.iP,CoolProp.iSmass,CoolProp.iT) -> (0,1,2,CoolProp.PSmass_INPUTS) + other values require switching and swapping. + """ + # Figure out if x or y-dimension should be used + switch = self.XY_SWITCH[self.i_index][self.y_index*10+self.x_index] + + if switch is None: + raise ValueError("This isoline cannot be calculated!") + elif switch is False: + pair, out1, _ = CP.generate_update_pair(self.i_index,0.0,self.x_index,1.0) + elif switch is True: + pair, out1, _ = CP.generate_update_pair(self.i_index,0.0,self.y_index,1.0) + else: + raise ValueError("Unknown error!") + + if out1==0.0: # Correct order + swap = False + else: # Wrong order + swap = True + + if not switch and not swap: + return 0,1,2,pair + elif switch and not swap: + return 0,2,1,pair + elif not switch and swap: + return 1,0,2,pair + elif switch and swap: + return 1,2,0,pair + else: + raise ValueError("Check the code, this should not happen!") + + def calc_sat_range(self,Trange=None,Prange=None,num=200): + if Trange is not None: + two = np.array(Trange) + one = np.resize(np.array(self.value),two.shape) + pair = CoolProp.QT_INPUTS + elif Prange is not None: + one = np.array(Prange) + two = np.resize(np.array(self.value),one.shape) + pair = CoolProp.PQ_INPUTS + else: + T_lo,T_hi = self._get_sat_bounds(CoolProp.iT) + two = np.linspace(T_lo,T_hi,num) + one = np.resize(np.array(self.value),two.shape) + pair = CoolProp.QT_INPUTS + + Tcrit = self.state.trivial_keyed_output(CoolProp.iT_critical) + Pcrit = self.state.trivial_keyed_output(CoolProp.iP_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: + warnings.warn( + "An error occurred for the critical inputs, skipping it.", + UserWarning) + xcrit = np.NaN + ycrit = np.NaN + + X = np.empty_like(one) + Y = np.empty_like(one) + + err = False + for index, _ in np.ndenumerate(one): + try: + self.state.update(pair, one[index], two[index]) + X[index] = self.state.keyed_output(self._x_index) + Y[index] = self.state.keyed_output(self._y_index) + except Exception as e: + if (pair == CoolProp.QT_INPUTS and abs(two[index]-Tcrit)<1e0) or \ + (pair == CoolProp.PQ_INPUTS and abs(one[index]-Pcrit)<1e2): + X[index] = xcrit + Y[index] = ycrit + warnings.warn( + "An error occurred for near critical inputs {0:f}, {1:f} with index {2:s}: {3:s}".format(one[index],two[index],str(index),str(e)), + UserWarning) + pass + + warnings.warn( + "An error occurred for inputs {0:f}, {1:f} with index {2:s}: {3:s}".format(one[index],two[index],str(index),str(e)), + UserWarning) + X[index] = np.NaN + Y[index] = np.NaN + err = True + self.x = X; self.y = Y + return + + def calc_range(self,xvals=None,yvals=None): + + if self.i_index == CoolProp.iQ: + warnings.warn( + "Please use \"calc_sat_range\" to calculate saturation and isoquality lines. Input ranges are discarded.", + UserWarning) + if xvals is not None: self.calc_sat_range(num=xvals.size) + elif yvals is not None: self.calc_sat_range(num=yvals.size) + else: self.calc_sat_range() + return + + ipos,xpos,ypos,pair = self.get_update_pair() + + order = [ipos,xpos,ypos] + idxs = [v for (_,v) in sorted(zip(order,[self.i_index , self.x_index , self.y_index ]))] + vals = [v for (_,v) in sorted(zip(order,[np.array(self.value), xvals , yvals ]))] + if vals[0] is None or vals[1] is None: + raise ValueError("One required input is missing, make sure to supply the correct xvals ({0:s}) or yvals ({1:s}).".format(str(xvals),str(yvals))) + + if vals[0].size > vals[1].size: + vals[1] = np.resize(vals[1],vals[0].shape) + elif vals[0].size < vals[1].size: + vals[0] = np.resize(vals[0],vals[1].shape) + + vals[2] = np.empty_like(vals[0]) + err = False + for index, _ in np.ndenumerate(vals[0]): + try: + self.state.update(pair, vals[0][index], vals[1][index]) + vals[2][index] = self.state.keyed_output(idxs[2]) + except Exception as e: + warnings.warn( + "An error occurred for inputs {0:f}, {1:f} with index {2:s}: {3:s}".format(vals[0][index],vals[1][index],str(index),str(e)), + UserWarning) + vals[2][index] = np.NaN + err = True + + for i,v in enumerate(idxs): + if v == self.x_index: self.x = vals[i] + if v == self.y_index: self.y = vals[i] + + + def sanitize_data(self): + """Fill the series via interpolation""" + validx = None; validy = None + countx = None; county = None + if self.x is not None: + validx = np.sum(np.isfinite(self.x)) + countx = float(self.x.size) + else: + raise ValueError("The x-axis is not populated, calculate values before you interpolate.") + if self.y is not None: + validy = np.sum(np.isfinite(self.y)) + county = float(self.y.size) + else: + raise ValueError("The y-axis is not populated, calculate values before you interpolate.") + + if min([validx/countx,validy/county]) < self.VALID_REQ: + warnings.warn( + "Poor data quality, there are not enough valid entries for x ({0:f}/{1:f}) or y ({2:f}/{3:f}).".format(validx,countx,validy,county), + UserWarning) + # TODO: use filter and cubic splines! + #filter = np.logical_and(np.isfinite(self.x),np.isfinite(self.y)) + if validy > validx: + y = self.y[np.isfinite(self.y)] + self.x = interp1d(self.y, self.x, kind='linear')(y) + self.y = y + else: + x = self.x[np.isfinite(self.x)] + self.y = interp1d(self.x, self.y, kind='linear')(x) + self.x = x + + + + +class BasePlot(Base2DObject): + """The base class for all plots. It can be instantiated itself, but provides many + general facilities to be used in the different plots. """ + + # Define the iteration keys + PROPERTIES = { + CoolProp.iDmass: 'density', + CoolProp.iHmass: 'specific enthalpy', + CoolProp.iP : 'pressure', + CoolProp.iSmass: 'specific entropy', + CoolProp.iT : 'temperature', + CoolProp.iUmass: 'specific internal energy' + } + + # Define the unit systems + UNIT_SYSTEMS = { + 'SI' : SIunits(), + 'KSI': KSIunits(), + 'EUR': EURunits() + } + + LINE_PROPS = { + CoolProp.iT : dict(color='Darkred' ,lw=0.25), + CoolProp.iP : dict(color='DarkCyan' ,lw=0.25), + CoolProp.iHmass: dict(color='DarkGreen' ,lw=0.25), + CoolProp.iDmass: dict(color='DarkBlue' ,lw=0.25), + CoolProp.iSmass: dict(color='DarkOrange',lw=0.25), + CoolProp.iQ : dict(color='black' ,lw=0.25) + } + + ID_FACTOR = 10.0 # Values below this number are interpreted as factors + HI_FACTOR = 2.25 # Upper default limits: HI_FACTOR*T_crit and HI_FACTOR*p_crit + LO_FACTOR = 1.01 # Lower default limits: LO_FACTOR*T_triple and LO_FACTOR*p_triple + + TP_LIMITS = { + 'NONE' : [None,None,None,None], + 'DEF' : [LO_FACTOR,HI_FACTOR,LO_FACTOR,HI_FACTOR], + 'ACHP' : [173.15,493.15,0.25e5,HI_FACTOR], + 'ORC' : [273.15,673.15,0.25e5,HI_FACTOR] + } + + def __init__(self, fluid_ref, graph_type, unit_system = 'KSI', tp_limits='DEF', **kwargs): + + # Process the graph_type and set self._x_type and self._y_type + graph_type = graph_type.upper() + graph_type = graph_type.replace(r'RHO',r'D') + if graph_type not in Base2DObject.PLOTS: + raise ValueError("Invalid graph_type input, expected a string from {0:s}".format(str(self.PLOTS))) + + # call the base class + state = process_fluid_state(fluid_ref) + Base2DObject.__init__(self, graph_type[1], graph_type[0], state, **kwargs) + + # Process the unit_system and set self._system + self.system = unit_system + # Process the plotting range based on T and p + self.limits = tp_limits + # Other properties + self.figure = kwargs.get('figure',plt.figure(tight_layout=True)) + self.axis = kwargs.get('axis', self.figure.add_subplot(111)) + self.props = kwargs.get('props', None) + + @property + def system(self): return self._system + @system.setter + def system(self, value): + value = value.upper() + if value in self.UNIT_SYSTEMS: self._system = self.UNIT_SYSTEMS[value] + else: raise ValueError("Invalid input, expected a string from {0:s}".format(str(self.UNIT_SYSTEMS.keys()))) + + @property + def limits(self): + """Returns [Tmin,Tmax,pmin,pmax] as value or factors""" + return self._limits + @limits.setter + def limits(self, value): + try: value = value.upper() + except: pass + if value in self.TP_LIMITS: self._limits = self.TP_LIMITS[value] + elif len(value)==4: self._limits = value + else: raise ValueError("Invalid input, expected a list with 4 items or a string from {0:s}".format(str(self.TP_LIMITS.keys()))) + + @property + def figure(self): return self._figure + @figure.setter + def figure(self, value): self._figure = value + + @property + def axis(self): return self._axis + @axis.setter + def axis(self, value): self._axis = value + + @property + def props(self): return self._props + @props.setter + def props(self, value): + self._props = self.LINE_PROPS.copy() + if value is not None: + self._props.update(value) + + def __sat_bounds(self, kind, smin=None, smax=None): + warnings.warn( + "You called the deprecated function \"__sat_bounds\", \ +consider replacing it with \"_get_sat_bounds\".", + DeprecationWarning) + return self._get_sat_bounds(kind, smin, smax) + + + def _get_iso_label(self, isoline, unit=True): + if self._system is not None: + dim = self._system[isoline.i_index] + return str(r"$"+dim.symbol+"="+str(dim.from_SI(isoline.value))+ "$ "+dim.unit if unit else "$").strip() + return str(isoline.value).strip() + + #def _get_phase_envelope(self): + # + #HEOS = CoolProp.AbstractState("HEOS", fluid) + #HEOS.build_phase_envelope("") + #PED = HEOS.get_phase_envelope_data() + #plt.plot(PED.T, np.log(PED.p)) + #plt.show() - return sat_lines def _plot_default_annotations(self): - def filter_fluid_ref(fluid_ref): - fluid_ref_string = fluid_ref - if fluid_ref.startswith('REFPROP-MIX'): - end = 0 - fluid_ref_string = '' - while fluid_ref.find('[', end + 1) != -1: - start = fluid_ref.find('&', end + 1) - if end == 0: - start = fluid_ref.find(':', end + 1) - end = fluid_ref.find('[', end + 1) - fluid_ref_string = ' '.join([fluid_ref_string, - fluid_ref[start+1:end], '+']) - fluid_ref_string = fluid_ref_string[0:len(fluid_ref_string)-2] - return fluid_ref_string - - if len(self.graph_type) == 2: - y_axis_id = self.graph_type[0] - x_axis_id = self.graph_type[1] - else: - y_axis_id = self.graph_type[0] - x_axis_id = self.graph_type[1:len(self.graph_type)] - - tl_str = "%s - %s Graph for %s" - if not self.axis.get_title(): - self.axis.set_title(tl_str % (self.AXIS_LABELS[self.unit_system][y_axis_id][0], - self.AXIS_LABELS[self.unit_system][x_axis_id][0], - filter_fluid_ref(self.fluid_ref))) +# def filter_fluid_ref(fluid_ref): +# fluid_ref_string = fluid_ref +# if fluid_ref.startswith('REFPROP-MIX'): +# end = 0 +# fluid_ref_string = '' +# while fluid_ref.find('[', end + 1) != -1: +# start = fluid_ref.find('&', end + 1) +# if end == 0: +# start = fluid_ref.find(':', end + 1) +# end = fluid_ref.find('[', end + 1) +# fluid_ref_string = ' '.join([fluid_ref_string, +# fluid_ref[start+1:end], '+']) +# fluid_ref_string = fluid_ref_string[0:len(fluid_ref_string)-2] +# return fluid_ref_string +# +# if len(self.graph_type) == 2: +# y_axis_id = self.graph_type[0] +# x_axis_id = self.graph_type[1] +# else: +# y_axis_id = self.graph_type[0] +# x_axis_id = self.graph_type[1:len(self.graph_type)] +# +# tl_str = "%s - %s Graph for %s" +# if not self.axis.get_title(): +# self.axis.set_title(tl_str % (self.AXIS_LABELS[self.unit_system][y_axis_id][0], +# self.AXIS_LABELS[self.unit_system][x_axis_id][0], +# filter_fluid_ref(self.fluid_ref))) + if self._x_index in [CoolProp.iDmass,CoolProp.iP]: + self.axis.set_xscale('log') + if self._y_index in [CoolProp.iDmass,CoolProp.iP]: + self.axis.set_yscale('log') + if not self.axis.get_xlabel(): - self.axis.set_xlabel(' '.join(self.AXIS_LABELS[self.unit_system][x_axis_id])) + dim = self._system[self._x_index] + self.xlabel((dim.label+u" $"+dim.symbol+u"$ / "+dim.unit).strip()) if not self.axis.get_ylabel(): - self.axis.set_ylabel(' '.join(self.AXIS_LABELS[self.unit_system][y_axis_id])) - - def _draw_graph(self): - return + dim = self._system[self._y_index] + self.ylabel((dim.label+u" $"+dim.symbol+u"$ / "+dim.unit).strip()) def title(self, title): self.axis.set_title(title) @@ -273,14 +735,297 @@ class BasePlot(object): else: self.axis.grid(kwargs) + + def set_Tp_limits(self, limits): + """Set the limits for the graphs in temperature and pressure, based on + the active units: [Tmin, Tmax, pmin, pmax]""" + dim = self._system[CoolProp.iT] + limits[0] = dim.to_SI(limits[0]) + limits[1] = dim.to_SI(limits[1]) + dim = self._system[CoolProp.iP] + limits[2] = dim.to_SI(limits[2]) + limits[3] = dim.to_SI(limits[3]) + self.limits = limits + + def get_Tp_limits(self): + """Get the limits for the graphs in temperature and pressure, based on + the active units: [Tmin, Tmax, pmin, pmax]""" + limits = self._get_Tp_limits() + dim = self._system[CoolProp.iT] + limits[0] = dim.from_SI(limits[0]) + limits[1] = dim.from_SI(limits[1]) + dim = self._system[CoolProp.iP] + limits[2] = dim.from_SI(limits[2]) + limits[3] = dim.from_SI(limits[3]) + return limits + + def _get_Tp_limits(self): + """Get the limits for the graphs in temperature and pressure, based on + SI units: [Tmin, Tmax, pmin, pmax]""" + T_lo,T_hi,P_lo,P_hi = self.limits + Ts_lo,Ts_hi = self._get_sat_bounds(CoolProp.iT) + Ps_lo,Ps_hi = self._get_sat_bounds(CoolProp.iP) + + if T_lo is None: T_lo = 0.0 + elif T_lo < self.ID_FACTOR: T_lo *= Ts_lo + if T_hi is None: T_hi = 1e6 + elif T_hi < self.ID_FACTOR: T_hi *= Ts_hi + if P_lo is None: P_lo = 0.0 + elif P_lo < self.ID_FACTOR: P_lo *= Ps_lo + if P_hi is None: P_hi = 1e10 + elif P_hi < self.ID_FACTOR: P_hi *= Ps_hi + + try: T_lo = np.nanmax([T_lo, self._state.trivial_keyed_output(CoolProp.iT_min)]) + except: pass + try: T_hi = np.nanmin([T_hi, self._state.trivial_keyed_output(CoolProp.iT_max)]) + except: pass + try: P_lo = np.nanmax([P_lo, self._state.trivial_keyed_output(CoolProp.iP_min)]) + except: pass + try: P_hi = np.nanmin([P_hi, self._state.trivial_keyed_output(CoolProp.iP_max)]) + except: pass + + return [T_lo,T_hi,P_lo,P_hi] + + def set_axis_limits(self, limits): + """Set the limits of the internal axis object based on the active units, + takes [xmin, xmax, ymin, ymax]""" self.axis.set_xlim([limits[0], limits[1]]) self.axis.set_ylim([limits[2], limits[3]]) + + def _set_axis_limits(self, limits): + """Set the limits of the internal axis object based on SI units, + takes [xmin, xmax, ymin, ymax]""" + dim = self._system[self._x_index] + self.axis.set_xlim([dim.from_SI(limits[0]), dim.from_SI(limits[1])]) + dim = self._system[self._y_index] + self.axis.set_ylim([dim.from_SI(limits[2]), dim.from_SI(limits[3])]) + + + def get_axis_limits(self,x_index=None,y_index=None): + """Returns the previously set limits or generates them and + converts the default values to the selected unit system. + Returns a list containing [xmin, xmax, ymin, ymax]""" + if x_index is None: x_index = self._x_index + if y_index is None: y_index = self._y_index + + if x_index != self.x_index or y_index != self.y_index or \ + self.axis.get_autoscalex_on() or self.axis.get_autoscaley_on(): + # One of them is not set or we work on a different set of axes + T_lo,T_hi,P_lo,P_hi = self._get_Tp_limits() + + X=[0.0]*4; Y=[0.0]*4 + i = -1 + for T in [T_lo, T_hi]: + for P in [P_lo, P_hi]: + i+=1 + try: + self._state.update(CoolProp.PT_INPUTS, P, T) + # TODO: include a check for P and T? + X[i] = self._state.keyed_output(x_index) + Y[i] = self._state.keyed_output(y_index) + except: + X[i] = np.nan; Y[i] = np.nan + # Figure out what to update + dim = self._system[x_index] + x_lim = [dim.from_SI(np.nanmin(X)),dim.from_SI(np.nanmax(X))] + dim = self._system[y_index] + y_lim = [dim.from_SI(np.nanmin(Y)),dim.from_SI(np.nanmax(Y))] + # Either update the axes limits or get them + if x_index == self._x_index: + if self.axis.get_autoscalex_on(): + self.axis.set_xlim(x_lim) + else: + x_lim = self.axis.get_xlim() + if y_index == self._y_index: + if self.axis.get_autoscaley_on(): + self.axis.set_ylim(y_lim) + else: + y_lim = self.axis.get_ylim() + else: # We only asked for the real axes limits and they are set already + x_lim = self.axis.get_xlim() + y_lim = self.axis.get_ylim() + + return [x_lim[0],x_lim[1],y_lim[0],y_lim[1]] + + def _get_axis_limits(self,x_index=None,y_index=None): + """Get the limits of the internal axis object in SI units + Returns a list containing [xmin, xmax, ymin, ymax]""" + if x_index is None: x_index = self._x_index + if y_index is None: y_index = self._y_index + limits = self.get_axis_limits(x_index,y_index) + dim = self._system[x_index] + limits[0] = dim.to_SI(limits[0]) + limits[1] = dim.to_SI(limits[1]) + dim = self._system[y_index] + limits[2] = dim.to_SI(limits[2]) + limits[3] = dim.to_SI(limits[3]) + return limits + + @staticmethod + def generate_ranges(itype, imin, imax, num): + """Generate a range for a certain property""" + if itype in [CoolProp.iP, CoolProp.iDmass]: + return np.logspace(np.log2(imin),np.log2(imax),num=num,base=2.) + return np.linspace(imin, imax, num=num) + + def _get_conversion_data(self): + [Axmin,Axmax,Aymin,Aymax] = self._get_axis_limits() + DELTAX_axis=Axmax-Axmin + DELTAY_axis=Aymax-Aymin + width=self.figure.get_figwidth() + height=self.figure.get_figheight() + pos=self.axis.get_position().get_points() + [[Fxmin,Fymin],[Fxmax,Fymax]]=pos + DELTAX_fig=width*(Fxmax-Fxmin) + DELTAY_fig=height*(Fymax-Fymin) + return [[Axmin,Axmax,Aymin,Aymax,Fxmin,Fxmax,Fymin,Fymax],[DELTAX_axis,DELTAY_axis,DELTAX_fig,DELTAY_fig]] + + def _to_pixel_coords(self,xv,yv): + [[Axmin,Axmax,Aymin,Aymax,Fxmin,Fxmax,Fymin,Fymax],[DELTAX_axis,DELTAY_axis,DELTAX_fig,DELTAY_fig]] = self._get_conversion_data() + #Convert coords to pixels + x=(xv-Axmin)/DELTAX_axis*DELTAX_fig+Fxmin + y=(yv-Aymin)/DELTAY_axis*DELTAY_fig+Fymin + return x,y + + def _to_data_coords(self,xv,yv): + [[Axmin,Axmax,Aymin,Aymax,Fxmin,Fxmax,Fymin,Fymax],[DELTAX_axis,DELTAY_axis,DELTAX_fig,DELTAY_fig]] = self._get_conversion_data() + #Convert back to measurements + x=(xv-Fxmin)/DELTAX_fig*DELTAX_axis+Axmin + y=(yv-Fymin)/DELTAY_fig*DELTAY_axis+Aymin + return x,y + + @staticmethod + 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 np.min(xv)x>xv[index+1]): # nearest above, negative inclination + if diff[index]x>xv[index]): # nearest below, negative inclination + if diff[index] (0,1,2,CoolProp.PSmass_INPUTS) + #other values require switching and swapping + #get_update_pair(CoolProp.iSmass,CoolProp.iP,CoolProp.iHmass) -> (1,0,2,CoolProp.PSmass_INPUTS) \ No newline at end of file diff --git a/wrappers/Python/CoolProp/Plots/Plots.py b/wrappers/Python/CoolProp/Plots/Plots.py index cbf7d959..80edea2a 100644 --- a/wrappers/Python/CoolProp/Plots/Plots.py +++ b/wrappers/Python/CoolProp/Plots/Plots.py @@ -7,257 +7,96 @@ 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 # +import warnings +from CoolProp.Plots.Common import IsoLine,BasePlot +import CoolProp +import sys +from CoolProp.Plots.SimpleCycles import StatePoint, StateContainer,\ + SimpleRankineCycle -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>xv[index+1]): # nearest above, negative inclination - if diff[index]x>xv[index]): # nearest below, negative inclination - if diff[index]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): +class PropertyPlot(BasePlot): + def __init__(self, fluid_name, graph_type, **kwargs): """ - Generates limits for the axes in terms of x,y defined by 'plot' - based on temperature and pressure. + Create graph for the specified fluid properties - Returns a tuple containing ((xmin, xmax), (ymin, ymax)) + 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.iS, 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 """ - # 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): + 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 @@ -283,343 +122,249 @@ class IsoLines(BasePlot): 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): + + 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. """ - 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') + + 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 = 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]) - + iso_range = self.generate_ranges(iso_type, iso_range[0], iso_range[1], 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 + 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 = numpy.logical_and(numpy.isfinite(dew.x),numpy.isfinite(dew.y)) + #dew_filter = numpy.logical_and(dew_filter,dew.x>dew.x[-1]) + stp = min([dew_filter.size,10]) + dew_filter[0:-stp] = False + bub_filter = numpy.logical_and(numpy.isfinite(bub.x),numpy.isfinite(bub.y)) + + + if self._x_index == CoolProp.iP or self._x_index == CoolProp.iDmass: + filter_x = lambda x: numpy.log10(x) + else: + filter_x = lambda x: x + if self._y_index == CoolProp.iP or self._y_index == CoolProp.iDmass: + filter_y = lambda y: numpy.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)): + f = interp1d(numpy.append(bub.x[bub_filter],dew.x[dew_filter][::-1]),numpy.append(bub.y[bub_filter],dew.y[dew_filter][::-1]),kind='cubic') + x = numpy.linspace(bub.x[bub_filter][-1], dew.x[dew_filter][-1], 11) + y = f(x) + 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(numpy.array([bub.x[bub_filter][-1], dew.x[dew_filter][-1]])),dimy.from_SI(numpy.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(numpy.append(line.x,xcrit)),dimy.from_SI(numpy.append(line.y,ycrit)),**props) + #try: + # x = numpy.append(line.x,[xcrit]) + # y = numpy.append(line.y,[ycrit]) + # fltr = numpy.logical_and(numpy.isfinite(x),numpy.isfinite(y)) + # f = interp1d(x[fltr][-3:],y[fltr][-3:],kind='linear') # could also be quadratic + # x = numpy.linspace(x[fltr][-2], x[fltr][-1], 5) + # y = f(x) + # #f = interp1d(y[fltr][-5:],x[fltr][-5:],kind='cubic') + # #y = numpy.linspace(y[fltr][-2], y[fltr][-1], 5) + # #x = f(y) + # self.axis.plot(dimx.from_SI(numpy.append(line.x,x)),dimy.from_SI(numpy.append(line.y,y)),**props) + #except: + # self.axis.plot(dimx.from_SI(numpy.append(line.x,xcrit)),dimy.from_SI(numpy.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={'color' : 'r', 'lw' : 1.5}): + """ Draw process or cycle from x and y values in axis units 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 - + 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 -------- - >>> 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) + >>> 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) + + + 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(numpy.asarray(xdata)) + ydata = dimy.from_SI(numpy.asarray(ydata)) + self.axis.plot(xdata,ydata,marker='None',linestyle=style,**line_opts) + + xdata = numpy.empty(len(points)) + ydata = numpy.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(numpy.asarray(xdata)) + ydata = dimy.from_SI(numpy.asarray(ydata)) + line_opts['label'] = '' + self.axis.plot(xdata,ydata,marker=marker,linestyle='None',**line_opts) - 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 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) - 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) +class PropsPlot(PropertyPlot): + def __init__(self, fluid_name, graph_type, units = 'KSI', reciprocal_density = False, **kwargs): + super(PropsPlot, self).__init__(fluid_name, graph_type, units=units, reciprocal_density=reciprocal_density, **kwargs) + warnings.warn("You called the deprecated class \"PropsPlot\", use \"PropertyPlot\".",DeprecationWarning) -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 +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() + diff --git a/wrappers/Python/CoolProp/Plots/SimpleCycles.py b/wrappers/Python/CoolProp/Plots/SimpleCycles.py index 621b0e93..5ba13744 100644 --- a/wrappers/Python/CoolProp/Plots/SimpleCycles.py +++ b/wrappers/Python/CoolProp/Plots/SimpleCycles.py @@ -1,7 +1,19 @@ +# -*- coding: utf-8 -*- + +from __future__ import print_function, division +from six import with_metaclass + import matplotlib,numpy +import numpy as np + +import CoolProp from CoolProp.CoolProp import PropsSI from scipy.optimize import newton +from .Common import BasePlot, process_fluid_state, PropertyDict, SIunits +import warnings +from abc import ABCMeta + def SimpleCycle(Ref,Te,Tc,DTsh,DTsc,eta_a,Ts_Ph='Ph',skipPlot=False,axis=None): """ @@ -23,6 +35,9 @@ def SimpleCycle(Ref,Te,Tc,DTsh,DTsc,eta_a,Ts_Ph='Ph',skipPlot=False,axis=None): * skipPlot : If True, won't actually plot anything, just print COP """ + + warnings.warn("This function has been deprecated. PLease consider converting it to an object inheriting from \"BaseCycle\".",DeprecationWarning) + T=numpy.zeros((6)) h=numpy.zeros_like(T) p=numpy.zeros_like(T) @@ -107,6 +122,8 @@ def TwoStage(Ref,Q,Te,Tc,DTsh,DTsc,eta_oi,f_p,Tsat_ic,DTsh_ic,Ts_Ph='Ph',prints= * skipPlot : If True, won't actually plot anything, just print COP """ + + warnings.warn("This function has been deprecated. PLease consider converting it to an object inheriting from \"BaseCycle\".",DeprecationWarning) T=numpy.zeros((8)) h=numpy.zeros_like(T) @@ -244,6 +261,8 @@ def EconomizedCycle(Ref,Qin,Te,Tc,DTsh,DTsc,eta_oi,f_p,Ti,Ts_Ph='Ts',skipPlot=Fa * skipPlot : If True, won't actually plot anything, just print COP """ + + warnings.warn("This function has been deprecated. PLease consider converting it to an object inheriting from \"BaseCycle\".",DeprecationWarning) m=1 @@ -389,34 +408,563 @@ def EconomizedCycle(Ref,Qin,Te,Tc,DTsh,DTsc,eta_oi,f_p,Ti,Ts_Ph='Ts',skipPlot=Fa print('Vdisp2: ',(mdot+mdot_inj)/(rho[1]*f*eta_v)*1e6,'cm^3') return COP -if __name__=='__main__': - from CoolProp.Plots import Ph,Ts - Ref='R290' - fig=matplotlib.pyplot.figure(figsize=(4,3)) - ax=fig.add_axes((0.15,0.15,0.8,0.8)) - Ph(Ref,Tmin=273.15-30,hbounds=[0,600],axis=ax) - COP=TwoStage('Propane',10000,273.15-5,273.15+43.3,5,7,0.7,0.3,15+273.15,3,prints = True) - matplotlib.pyplot.show() + #class SimpleCycle(object): + # """A class that calculates a simple thermodynamic cycle""" + # def __init__(self, *args, **kwargs): + # object.__init__(self, *args, **kwargs) + # (states, steps, fluid): + +# Parameters +# ---------- +# x_type : int, str +# Either a letter or an integer that specifies the property type for the x-axis +# y_type : int, str +# Either a letter or an integer that specifies the property type for the y-axis +# states : list +# A collection of state points that follows a fixed scheme defined +# in the implementing subclass. +# fluid_ref : str, CoolProp.AbstractState +# The fluid property provider, either a subclass of CoolProp.AbstractState +# or a string that can be used to generate a CoolProp.AbstractState instance +# via :func:`Common.process_fluid_state`. +# steps : int +# The number of steps used for going from one state to another +# +# for more properties, see :class:`CoolProp.Plots.Common.Base2DObject`. - Ref='R290' - fig=matplotlib.pyplot.figure(figsize=(4,3)) - ax=fig.add_axes((0.15,0.15,0.8,0.8)) - Ph(Ref,Tmin=273.15-30,hbounds=[0,600],axis=ax) - COP=SimpleCycle(Ref,273.15-5,273.15+45,5,7,0.7,Ts_Ph='Ph') - matplotlib.pyplot.show() +# # See http://stackoverflow.com/questions/1061283/lt-instead-of-cmp +# class ComparableMixin: +# """A mixin class that implements all comparing mathods except for __lt__""" +# def __eq__(self, other): +# return not selfB: return 1 + elif A>> from __future__ import print_function + >>> import CoolProp + >>> from CoolProp.PLots.SimpleCycles import StateContainer + >>> T0 = 300.000; p0 = 200000.000; h0 = 112745.749; s0 = 393.035 + >>> cycle_states = StateContainer() + >>> cycle_states[0,'H'] = h0 + >>> cycle_states[0]['S'] = s0 + >>> cycle_states[0][CoolProp.iP] = p0 + >>> cycle_states[0,CoolProp.iT] = T0 + >>> cycle_states[1,"T"] = 300.064 + >>> print(cycle_states) + Stored State Points: + state T (K) p (Pa) ρ (kg/m³) h (J/kg) s (J/kg/K) + 0 300.000 200000.000 996.601 112745.749 393.035 + 1 300.064 - - - - + + """ + + def __init__(self,unit_system=SIunits()): + self._points = {} + self._units = unit_system + + @property + def points(self): return self._points + @points.setter + def points(self, value): self._points = value + + @property + def units(self): return self._units + @units.setter + def units(self, value): self._units = value + + def get_point(self, index, SI=True): + if SI: + state = self[index] + else: + state = self[index] + for i in state: + state[i] = self.units[i].from_SI(state[i]) + return state + + def set_point(self, index, value, SI=True): + if SI: + self._points[index] = value + else: + for i in value: + self._points[index][i] = self.units[i].to_SI(value[i]) + + def _list_like(self, value): + """Try to detect a list-like structure excluding strings""" + return (not hasattr(value, "strip") and + (hasattr(value, "__getitem__") or + hasattr(value, "__iter__"))) + # return is_sequence(value) # use from pandas.core.common import is_sequence + + def __len__(self): + """Some cheating to get the correct behaviour""" + return len(self._points) + + def __iter__(self): + """Make sure we iterate in the righ order""" + for key in sorted(self._points): + yield key + + def __getitem__(self, index): + """Another tweak that changes the default access path""" + if self._list_like(index): + if len(index)==0: raise IndexError("Received empty index.") + elif len(index)==1: return self._points[index[0]] + elif len(index)==2: return self._points[index[0]][index[1]] + else: raise IndexError("Received too long index.") + return self._points[index] + + def __setitem__(self, index, value): + """Another tweak that changes the default access path""" + if self._list_like(index): + if len(index)==0: raise IndexError("Received empty index.") + elif len(index)==1: self._points[index[0]] = value + elif len(index)==2: + # safeguard against empty entries + if index[0] not in self._points: + self._points[index[0]] = StatePoint() + self._points[index[0]][index[1]] = value + else: raise IndexError("Received too long index.") + else: + self._points[index] = value + + def __str__(self): + out = "Stored State Points:\n" + keys = True + for i in self._points: + if keys: + row = ["{0:>5s}".format("state")] + for j in self._points[i]: + label = u"{0:s} ({1:s})".format(self.units[j].symbol,self.units[j].unit) + row.append(u"{0:>11s}".format(label)) + out = out + " ".join(row) + "\n" + keys = False + row = ["{0:>5s}".format(str(i))] + for j in self._points[i]: + try: + row.append(u"{0:11.3f}".format(self.units[j].from_SI(self._points[i][j]))) + except: + row.append(u"{0:>11s}".format("-")) + out = out + " ".join(row) + "\n" + return out.encode('utf8', 'replace') + + def append(self,new): + i = 0 + self.__len__() + for j in new: + self[i,j] = new[j] + return self + + def extend(self,new): + i = 0 + self.__len__() + for j in new: + for k in new[j]: + self[i,k] = new[j][k] + i = i +1 + return self + + @property + def D(self): return np.array([self._points[k].D for k in self]) + @property + def H(self): return np.array([self._points[k].H for k in self]) + @property + def P(self): return np.array([self._points[k].P for k in self]) + @property + def S(self): return np.array([self._points[k].S for k in self]) + @property + def T(self): return np.array([self._points[k].T for k in self]) + @property + def U(self): return np.array([self._points[k].U for k in self]) + @property + def Q(self): return np.array([self._points[k].Q for k in self]) + + + +class BaseCycle(BasePlot): + """A simple thermodynamic cycle, should not be used on its own.""" + + # Define the iteration keys + PROPERTIES = { + CoolProp.iDmass : 'density', + CoolProp.iHmass : 'specific enthalpy', + CoolProp.iP : 'pressure', + CoolProp.iSmass : 'specific entropy', + CoolProp.iT : 'temperature' + } + + STATECOUNT=0 + """A list of accepted numbers of states""" + + STATECHANGE=None + """A list of lists of tuples that defines how the state transitions + behave for the corresponding entry in BaseCycle.STATECOUNT""" + + def __init__(self, fluid_ref, graph_type, unit_system='EUR', **kwargs): + """Initialises a simple cycle calculator + + Parameters + ---------- + fluid_ref : str, CoolProp.AbstractState + The fluid property provider, either a subclass of CoolProp.AbstractState + or a string that can be used to generate a CoolProp.AbstractState instance + via :func:`Common.process_fluid_state`. + graph_type : string + The graph type to be plotted, like \"PH\" or \"TS\" + 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 + + for more properties, see :class:`CoolProp.Plots.Common.BasePlot`. + """ + self._cycle_states = StateContainer() + self._steps = 2 + BasePlot.__init__(self, fluid_ref, graph_type, unit_system, **kwargs) + + + @property + def cycle_states(self): return self._cycle_states + @cycle_states.setter + def cycle_states(self, value): + if len(value) != self.STATECOUNT: + raise ValueError("Your number of states ({0:d}) is not in the list of allowed state counts: {1:s}.".format(len(value),str(self.STATECOUNT))) + self._cycle_states = value + + @property + def steps(self): return self._steps + @steps.setter + def steps(self, value): self._steps = int(max([value,2])) + + @BasePlot.system.setter + def system(self, value): + if value in self.UNIT_SYSTEMS: + self._system = self.UNIT_SYSTEMS[value] + elif isinstance(value, PropertyDict): + self._system = value + else: + raise ValueError("Invalid unit_system input \"{0:s}\", expected a string from {1:s}".format(str(value),str(self.UNIT_SYSTEMS.keys()))) + self._cycle_states._units = self._system + + + def valid_states(self): + """Check the formats of BaseCycle.STATECOUNT and BaseCycle.STATECHANGE""" + if len(self.STATECHANGE) != self.STATECOUNT: + raise ValueError("Invalid number of states and or state change operations") + return True + + def fill_states(self,objs=None): + """Try to populate all fields in the state objects""" + + if objs is None: + objs = self._cycle_states + local = True + else: + local = False + + for i in objs: + full = True + for j in objs[i]: + if objs[i][j] is None: + full = False + if full: continue + if (objs[i][CoolProp.iDmass] is not None and + objs[i][CoolProp.iT] is not None): + self._state.update(CoolProp.DmassT_INPUTS, objs[i][CoolProp.iDmass], objs[i][CoolProp.iT]) + elif (objs[i][CoolProp.iP] is not None and + objs[i][CoolProp.iHmass] is not None): + self._state.update(CoolProp.HmassP_INPUTS, objs[i][CoolProp.iHmass], objs[i][CoolProp.iP]) + elif (objs[i][CoolProp.iP] is not None and + objs[i][CoolProp.iSmass] is not None): + self._state.update(CoolProp.PSmass_INPUTS, objs[i][CoolProp.iP], objs[i][CoolProp.iSmass]) + else: + warnings.warn("Please fill the state[{0:s}] manually.".format(str(i))) + continue + for j in objs[i]: + if objs[i][j] is None: + objs[i][j] = self._state.keyed_output(j) + + if local: self._cycle_states = objs + return objs + + + + def state_change(self,in1,in2,start,ty1='lin',ty2='lin'): + """Calculates a state change defined by the properties in1 and in2 + + Uses self.states[start] and self.states[start+1] (or self.states[0]) to define + the process and interpolates between the values. + + Parameters + ---------- + in1 : int + The index of the first defined property. + in2 : int + The index of the second defined property. + start : int + The index of the start state. + ty1 : str + The key that defines the type of state change for in1, lin or log. + ty2 : str + The key that defines the type of state change for in2, lin or log. + + Returns + ------- + scalar or array_like + a list of the length of self.steps+1 that describes the process. It includes start and end state. + """ + self.fill_states() + end = start + 1 + if end >= len(self.cycle_states): end -= len(self.cycle_states) + start = self.cycle_states[start] + end = self.cycle_states[end] + # + val = [] + inv = [in1,in2] + typ = [ty1,ty2] + for i,v in enumerate(inv): + if typ[i] == 'lin': + val.append(np.linspace(start[v], end[v], self.steps)) + elif typ[i] == 'log': + val.append(np.logspace(np.log10(start[v]), np.log10(end[v]), self.steps)) + else: + raise ValueError("Unknow range generator {0:s}".format(str(typ[i]))) + + sc = StateContainer(self._system) + for i,_ in enumerate(val[0]): + sc[i,inv[0]] = val[0][i] + sc[i,inv[1]] = val[1][i] + + return self.fill_states(sc) + + def get_state_change(self, index): + return self.STATECHANGE[index](self) + + def get_state_changes(self): + sc = self.get_state_change(0) + for i in range(1,self.STATECOUNT): + sc.extend(self.get_state_change(i)) + return sc + + +class BasePowerCycle(BaseCycle): + """A thermodynamic cycle for power producing processes. + + Defines the basic properties and methods to unify access to + power cycle-related quantities. + """ + + def __init__(self, fluid_ref='HEOS::Water', graph_type='TS', **kwargs): + """see :class:`CoolProp.Plots.SimpleCycles.BaseCycle` for details.""" + BaseCycle.__init__(self, fluid_ref, graph_type, **kwargs) + + def eta_carnot(self): + """Carnot efficiency + + Calculates the Carnot efficiency for the specified process, :math:`\eta_c = 1 - \frac{T_c}{T_h}`. + + Returns + ------- + float + """ + Tvector = self._cycle_states.T + return 1. - np.min(Tvector) / np.max(Tvector) + + def eta_thermal(self): + """Thermal efficiency + + The thermal efficiency for the specified process(es), :math:`\eta_{th} = \frac{\dot{W}_{exp} - \dot{W}_{pum}}{\dot{Q}_{in}}`. + + Returns + ------- + float + """ + raise NotImplementedError("Implement it in the subclass.") +class SimpleRankineCycle(BasePowerCycle): + """A simple Rankine cycle *without* regeneration""" + STATECOUNT=4 + STATECHANGE=[ + lambda inp: BaseCycle.state_change(inp,'S','P',0,ty1='log',ty2='log'), # Pumping process + lambda inp: BaseCycle.state_change(inp,'H','P',1,ty1='lin',ty2='lin'), # Heat addition + lambda inp: BaseCycle.state_change(inp,'H','P',2,ty1='log',ty2='log'), # Expansion + lambda inp: BaseCycle.state_change(inp,'H','P',3,ty1='lin',ty2='lin') # Heat removal + ] + + def __init__(self, fluid_ref='HEOS::Water', graph_type='TS', **kwargs): + """see :class:`CoolProp.Plots.SimpleCycles.BasePowerCycle` for details.""" + BasePowerCycle.__init__(self, fluid_ref, graph_type, **kwargs) + + def simple_solve(self, T0, p0, T2, p2, eta_exp, eta_pum, fluid=None, SI=True): + """" + A simple Rankine cycle calculation + + Parameters + ---------- + T0 : float + The coldest point, before the pump + p0 : float + The lowest pressure, before the pump + T2 : float + The hottest point, before the expander + p2 : float + The highest pressure, before the expander + eta_exp : float + Isentropic expander efficiency + eta_pum : float + Isentropic pump efficiency + + Examples + -------- + >>> import CoolProp + >>> from CoolProp.Plots import PropertyPlot + >>> pp = PropertyPlot('HEOS::Water', 'TS', unit_system='EUR') + >>> 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) + + """ + if fluid is not None: self.state = process_fluid_state(fluid) + if self._state is None: + raise ValueError("You have specify a fluid before you calculate.") + + cycle_states = StateContainer(unit_system=self._system) + + if not SI: + Tc = self._system[CoolProp.iT].to_SI + pc = self._system[CoolProp.iP].to_SI + T0 = Tc(T0) + p0 = pc(p0) + T2 = Tc(T2) + p2 = pc(p2) + + # Subcooled liquid + self.state.update(CoolProp.PT_INPUTS,p0,T0) + h0 = self.state.hmass() + s0 = self.state.smass() + # Just a showcase for the different accessor methods + cycle_states[0,'H'] = h0 + cycle_states[0]['S'] = s0 + cycle_states[0][CoolProp.iP] = p0 + cycle_states[0,CoolProp.iT] = T0 + + # Pressurised liquid + p1 = p2 + self.state.update(CoolProp.PSmass_INPUTS,p1,s0) + h1 = h0 + (self.state.hmass() - h0) / eta_pum + self.state.update(CoolProp.HmassP_INPUTS,h1,p1) + s1 = self.state.smass() + T1 = self.state.T() + cycle_states[1,'H'] = h1 + cycle_states[1,'S'] = s1 + cycle_states[1,'P'] = p1 + cycle_states[1,'T'] = T1 + + # Evaporated vapour + self.state.update(CoolProp.PT_INPUTS,p2,T2) + h2 = self.state.hmass() + s2 = self.state.smass() + cycle_states[2,'H'] = h2 + cycle_states[2,'S'] = s2 + cycle_states[2,'P'] = p2 + cycle_states[2,'T'] = T2 + + # Expanded gas + p3 = p0 + self.state.update(CoolProp.PSmass_INPUTS,p3,s2) + h3 = h2 - eta_exp * (h2 - self.state.hmass()) + self.state.update(CoolProp.HmassP_INPUTS,h3,p3) + s3 = self.state.smass() + T3 = self.state.T() + cycle_states[3,'H'] = h3 + cycle_states[3,'S'] = s3 + cycle_states[3,'P'] = p3 + cycle_states[3,'T'] = T3 + + w_net = h2 - h3 + q_boiler = h2 - h1 + eta_th = w_net / q_boiler + + self.cycle_states = cycle_states + self.fill_states() + + + def eta_thermal(self): + """Thermal efficiency + + The thermal efficiency for the specified process(es), :math:`\eta_{th} = \frac{\dot{W}_{exp} - \dot{W}_{pum}}{\dot{Q}_{in}}`. + + Returns + ------- + float + """ + w_net = self.cycle_states[2].H - self.cycle_states[3].H - (self.cycle_states[1].H - self.cycle_states[0].H) + q_boiler = self.cycle_states[2].H - self.cycle_states[1].H + return w_net / q_boiler + - -## for x in numpy.linspace(0,1): -## Ref='REFPROP-MIX:R152A[%g]&R32[%g]' %(x,1-x) -## COP=SimpleCycle(273.15+8,273.15+44,5,7,0.7,skipPlot=True,Ts_Ph='Ph') -## matplotlib.pyplot.show() diff --git a/wrappers/Python/CoolProp/Plots/__init__.py b/wrappers/Python/CoolProp/Plots/__init__.py index cde2a3d7..55bee4df 100644 --- a/wrappers/Python/CoolProp/Plots/__init__.py +++ b/wrappers/Python/CoolProp/Plots/__init__.py @@ -1,7 +1,9 @@ #Bring some functions into the Plots namespace for code concision from __future__ import absolute_import -from .Plots import PropsPlot, IsoLines, drawIsoLines -from .Plots import Ph, Ts, Ps, PT, Prho, Trho, hs +from .Plots import PropertyPlot,PropsPlot +from .Common import IsoLine +#from .Plots import PropsPlot, IsoLines, drawIsoLines +#from .Plots import Ph, Ts, Ps, PT, Prho, Trho, hs from .SimpleCycles import SimpleCycle, TwoStage, EconomizedCycle diff --git a/wrappers/Python/CoolProp/tests/test_Props.py b/wrappers/Python/CoolProp/tests/test_Props.py index 95f1bf53..b89fd995 100644 --- a/wrappers/Python/CoolProp/tests/test_Props.py +++ b/wrappers/Python/CoolProp/tests/test_Props.py @@ -2,7 +2,7 @@ import unittest from CoolProp.CoolProp import PropsSI import CoolProp import numpy as np - + def test_input_types(): for Fluid in ['Water']: for Tvals in [0.5*PropsSI(Fluid,'Tmin')+0.5*PropsSI(Fluid,'Tcrit'),