format long
syms x z w
L=10;B=1.0;T=0.625;g=9.81;p=1000;        %L、B、T為主尺度；g為重力加速度；p為密度；分別采用國際單位
fx=eval(diff('(B/2)*(1-(2*x/L)^2)*(1-(z/T)^2)','x'));   %船體表面方程f對(duì)x的偏導(dǎo)
fz=eval(diff('(B/2)*(1-(2*x/L)^2)*(1-(z/T)^2)','z'));   %船體表面方程f對(duì)z的偏導(dǎo)
S1=inline('sqrt(1+fx.^2+fz.^2)');            %積分S的被積函數(shù)表達(dá)式
S=2*dblquad(S1,-L/2,L/2,-T,0);           %船體靜浮狀態(tài)下的濕表面積S
i=0;
%計(jì)算不同航速下的Fn、Cw
for c=1.5:0.02:6;                        %c為船速 1.5m/s～6m/s
K=g/(c^2);                               %K為波數(shù)
i=i+1;
G1=(1-(2*x/L)^2)*(1-(z/T)^2)*(exp(K*(sec(w))^2*z))*cos(K*sec(w)*x);   %積分G的被積函數(shù)表達(dá)式
G2=int(G1,x,-L/2,L/2);                   %計(jì)算G的內(nèi)層積分x
G3=int(G2,z,-T,0);                       %計(jì)算G的外層積分z 
G=B/2*G3;                                %求G
R1=inline('(G^2)*(sec(w))^5');             %積分R的被積函數(shù)表達(dá)式
R=(4*p*g*(K^3)/pi)*quadl((G^2)*(sec(w))^5,0,pi/2);     %求R
Cw(i)=0.5*R/(p*c^2*S);                   %不同航速下的興波阻力系數(shù)
Fn(i)=c/sqrt(g*L);                       %不同航速下的弗氏數(shù)
end
plot(Fn,Cw);                             %繪出Fn-Cw曲線
xlabel('Fn'),ylabel('Cw');
title('興波阻力系數(shù)隨航速的變化');