h = [18.8 10 7.9 5.1 2.9 2 1]; %水银测压计高度差
dp1 = (1.3332237-0.098)*h; %1cmHg=1.3332237kpa,1cmH2O=0.098kpa
dp2 = [23.485 15.647 10.070 6.494 3.5441 2.6526 1.2563]; %测压表直接读出的差压
Q1 = [10.76 8.53 6.76 5.31 3.94 3.42 1.88]; %流量，单位m^3/h
Q = Q1./3; %单位是L/s
L =2.5; %管道长度m
d = 0.0265; % 管道内径m
u = 1.0050*0.001; %20度下,水的动力粘性系数Pa*s
V = 4*Q*10^(-3)/(3.1415926*d^2); %平均速度,m/s
Re = V*d*10^(3)/u; %雷诺数
f1 = 2*10^(3) *d^(3)*1000*dp1./(L*u^(2)*Re.*Re); %(水银测压计)
f2 = 2*10^(3) *d^(3)*1000*dp2./(L*u^(2)*Re.*Re); %(测压表)
f3=0.11.*((0.015*10^(-3)/0.0265+68./Re).^0.25); %经验公式
p1 = polyfit(log10(Re), log10(f1), 2); %lgf1-lgRe拟合函数系数
p2 = polyfit(log10(Re), log10(f2), 2); %lgf2-lgRe拟合函数系数
A1 = plot(log10(Re), log10(f1),  'r.'),xlabel('lg(Re)'),ylabel('lg(f)') %水银测压计的散点图
hold on;
A2 = plot(log10(Re), log10(f2), 'b.') %测压表的散点图
T1 = plot(log10(Re), polyval(p1, log10(Re)), 'r') %水银测压计的拟合图
T2 = plot(log10(Re), polyval(p2, log10(Re)), 'b') %测压表的拟合图
T3=plot(log10(Re),log10(f3),'k');
hold off

%水平不到位，有问题告诉请我一声 xzx
%刚刚dp1改了一下
%刚刚把ln改成了lg
%Q1直接写从表上读的数，立方米每小时的那个
%之前的那个换算错误，少除了3.6
%再改是🐕
%by 车辆71spermwhale