damyarou

python, GMT などのプログラム

Python コンター図作成

コンター図を作成する必要があったため、そのプログラムを作成した。

作例

f:id:damyarou:20200110160608j:plain

プログラムソース

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()

以 上