import numpy as np
import math
import matplotlib.pyplot as plt

#根据https://wap.sciencenet.cn/blog-531760-1053322.html
#构建Sajben扩压器上、下表面的曲线。

def fuc_xi(x,C1,C2,C3,C4):
  return C1*x*(1.0+C2*pow(x,C3))/pow((1-x),C4)
##enddef

def fuc_h(x,a):
  return a*np.cosh(x)/(a-1.0+np.cosh(x))
##enddef

alpha_c=1.4114;
l_c    =-2.598;#收缩段起点x无量纲值（除以h_th）
C1_c  =0.81;
C2_c  =1.00;
C3_c  =0.50;
C4_c  =0.60;

alpha_d=1.5;
l_d    =7.216;#扩张段终点x无量纲值（除以h_th）
C1_d  =2.25;
C2_d  =0.00;
C3_d  =0.0;
C4_d  =0.60;   

l_min =-4.1#整体起点x无量纲值（除以h_th）
l_max =8.66#整体终点x无量纲值（除以h_th）
h_th=44.0##喉道高度，mm

##收缩段上表面
start_c=0.0
stop_c=l_c
num_c=500 ##收缩段点数
#不考虑终止点（收缩段左端点），否则会后续出现分母为零。
x_c=np.linspace(start = start_c, stop = stop_c, num = num_c,endpoint = False)
x_c_b=x_c/l_c
xi_c=fuc_xi(x_c_b,C1_c,C2_c,C3_c,C4_c)
h_c=fuc_h(xi_c,alpha_c)
x1=np.append(x_c, l_min)*h_th
y1=np.append(h_c, alpha_c)*h_th
x1=x1[::-1]
y1=y1[::-1]
z1=np.zeros(y1.shape)
plt.plot(x1,y1)
plt.show()


##扩张段上表面
start_d=0.0
stop_d=l_d
num_d=500##扩张段点数
#不考虑终止点（扩张段右端点），否则会后续出现分母为零。
x_d=np.linspace(start = start_d, stop = stop_d, num = num_d,endpoint = False)
x_d_b=x_d/l_d
xi_d=fuc_xi(x_d_b,C1_d,C2_d,C3_d,C4_d)
h_d=fuc_h(xi_d,alpha_d)
x2=np.append(x_d, l_max)*h_th
y2=np.append(h_d, alpha_d)*h_th
z2=np.zeros(y2.shape)
plt.plot(x2,y2)
plt.show()


##下表面直线段
x3=np.array([l_min,l_max])*h_th
y3=np.zeros(x3.shape)
z3=y3

x=np.append(x1, x2[1:])#追加x2时，略过x2[0]，跟x1[]最后一个重复
y=np.append(y1, y2[1:])#追加y2时，略过y2[0]，跟y1[]最后一个重复
plt.plot(x,y)
plt.show()

f = open('./SajbenTopBotCurves.txt', 'w', encoding='utf-8-sig') # 打开文件， 用'w'写文件

f.write('#下表面直线段,x[mm]  y[mm]\n')
data=np.c_[x3,y3,z3]
np.savetxt(f, X=data, delimiter=' ', encoding='utf-8-sig')

f.write('\n#收缩段上表面,x[mm]  y[mm]\n')
data=np.c_[x1,y1,z1]
np.savetxt(f, X=data, delimiter=' ', encoding='utf-8-sig')

f.write('\n#扩张段上表面,x[mm]  y[mm]\n')
data=np.c_[x2,y2,z2]
np.savetxt(f, X=data, delimiter=' ', encoding='utf-8-sig')

f.close()