Python コンター図作成
コンター図を作成する必要があったため、そのプログラムを作成した。
作例
プログラムソース
import numpy as np import matplotlib.pyplot as plt import seaborn as sns def drawfig(xx,yy,zh,zc): q=68.64 vp1=7.0; dp1=np.sqrt(4*q/np.pi/vp1) vp2=2.5; dp2=np.sqrt(4*q/np.pi/vp2) vh1=5.0; dh1=np.sqrt(4*q/np.pi/vh1) vh2=2.0; dh2=np.sqrt(4*q/np.pi/vh2) fsz=16 xmin=3; xmax=6 ymin=4; ymax=7 plt.figure(figsize=(20,10),facecolor='w') plt.rcParams['font.size']=fsz plt.rcParams['font.family']='sans-serif' plt.subplot(121) tstr='Total Head loss Contour of Headrace and Penstock (unit: m)' plt.xlim([xmin,xmax]) plt.ylim([ymin,ymax]) plt.xlabel('Penstock diameter (m)') plt.ylabel('Headrace diameter (m)') plt.grid(color='#999999',linestyle='dashed') plt.gca().set_aspect('equal',adjustable='box') plt.title(tstr,loc='left',fontsize=fsz) plt.fill([xmin,xmax,xmax,xmin],[ymin,ymin,ymax,ymax],color='#aaaaaa',alpha=0.3) plt.fill([dp1,dp2,dp2,dp1],[dh1,dh1,dh2,dh2],color='#ffffff') co0='#555555' xs=dp1; ys=0.5*(dh1+dh2); ss='$v_P$=7.0m/s' plt.text(xs,ys,ss,color=co0,va='center',ha='left',rotation=90,fontsize=fsz) xs=dp2; ys=0.5*(dh1+dh2); ss='$v_P$=2.5m/s' plt.text(xs,ys,ss,color=co0,va='center',ha='right',rotation=90,fontsize=fsz) xs=4.25; ys=dh1; ss='$v_H$=5.0m/s' plt.text(xs,ys,ss,color=co0,va='bottom',ha='center',rotation=0,fontsize=fsz) xs=4.25; ys=dh2; ss='$v_H$=2.0m/s' plt.text(xs,ys,ss,color=co0,va='top',ha='center',rotation=0,fontsize=fsz) co1='black' co2='red' cont=plt.contour(xx,yy,zh,colors=[co1,co1,co2,co1,co1,co1,co1,co1],levels=[10,15,16.58,20,30,40,60,80]) cont.clabel(fmt='%1.2f', fontsize=fsz-1) plt.subplot(122) tstr='Total Cost Contour of Headrace and Penstock (unit: mil.NRP)' plt.xlim([xmin,xmax]) plt.ylim([ymin,ymax]) plt.xlabel('Penstock diameter (m)') plt.ylabel('Headrace diameter (m)') plt.grid(color='#999999',linestyle='dashed') plt.gca().set_aspect('equal',adjustable='box') plt.title(tstr,loc='left',fontsize=fsz) plt.fill([xmin,xmax,xmax,xmin],[ymin,ymin,ymax,ymax],color='#aaaaaa',alpha=0.3) plt.fill([dp1,dp2,dp2,dp1],[dh1,dh1,dh2,dh2],color='#ffffff') xs=dp1; ys=0.5*(dh1+dh2); ss='$v_P$=7.0m/s' plt.text(xs,ys,ss,color=co0,va='center',ha='left',rotation=90,fontsize=fsz) xs=dp2; ys=0.5*(dh1+dh2); ss='$v_P$=2.5m/s' plt.text(xs,ys,ss,color=co0,va='center',ha='right',rotation=90,fontsize=fsz) xs=4.75; ys=dh1; ss='$v_H$=5.0m/s' plt.text(xs,ys,ss,color=co0,va='bottom',ha='center',rotation=0,fontsize=fsz) xs=4.75; ys=dh2; ss='$v_H$=2.0m/s' plt.text(xs,ys,ss,color=co0,va='top',ha='center',rotation=0,fontsize=fsz) cont=plt.contour(xx,yy,zc,colors=[co1],levels=[6000,6500,7000,7500,8000,8200,8500,9000]) cont.clabel(fmt='%1.0f', fontsize=fsz-1) plt.contour(xx,yy,zh,colors=[co2],levels=[16.58]) xs=5.6; ys=5.65; ss='$h_L$=16.58m' plt.text(xs,ys,ss,color='#ff0000',va='bottom',ha='center',rotation=0,fontsize=fsz) fnameF='fig_contour.jpg' plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1) plt.show() def headrace(dd): a1=164619 # NRP/m b1=450019 # NRP/m q=68.64 al=10149.0 n=0.016 v=4*q/np.pi/dd**2 hl=124.5*n**2*al/dd**(4/3)*v**2/2/9.8 cc=(a1*dd**2/6.0**2+b1*dd/6.0)*al return hl,cc def penstock(dd): a2=143829 b2=2437167 cf=177919810 q=68.64 al=688.2 n=0.011 rho=20.0 fbr=0.4 v=4*q/np.pi/dd**2 h1=124.5*n**2*al/dd**(4/3) h2=4*(0.131+0.1632*(dd/rho)**(7/2)) h3=fbr hl=(h1+h2+h3)*v**2/2/9.8+1.991 cc=(a2*dd**2/4.0**2+b2*dd/4.0)*al+cf return hl,cc def main(): ds=0.1 d1=np.arange(4.0,7.0+ds,ds) d2=np.arange(3.0,6.0+ds,ds) hlt=np.zeros((len(d1),len(d2)),dtype=np.float64) cct=np.zeros((len(d1),len(d2)),dtype=np.float64) for i,dd1 in enumerate(d1): hl1,cc1=headrace(dd1) for j,dd2 in enumerate(d2): hl2,cc2=penstock(dd2) hlt[i,j]=(hl1+hl2)*1.1 #cct[i,j]=(cc1+cc2)/(450-hlt[i,j])*433.42/1e6 cct[i,j]=(cc1+cc2)/1e6 if 16.0<=hlt[i,j]<=16.58: print('{0:6.1f}{1:6.1f}{2:8.3f}{3:8.3f}{4:8.3f}{5:12.3f}'.format(dd1,dd2,hl1*1.1,hl2*1.1,hlt[i,j],cct[i,j])) xx, yy = np.meshgrid(d2,d1) zh=hlt zc=cct drawfig(xx,yy,zh,zc) if __name__ == '__main__': main()
以 上