clc
clear
syms?w?we;
w=0.4:0.05:1.6;mu=90;v=6;%kb=1;kt=1;%航速6m/s，航向90度，即橫浪，cos（90）=0
T=3;B=10;Sw=0.785;%船寬10米，吃水3米水線面系數假設為0.785
we=w+w.^2.*v/9.8;
for?i=1:24
????delta_we（i）=we（i+1）-we（i）;
end
delta_we（25）=delta_we（24）;
z=0.194375;
s_kesi=0.778./（w.^5）.*exp（-z./（w.^4））;?%含風速=15m/s即三一波高=4m在內了?%2.47./（w.^6）.*exp（-2.*9.8.^2./24.5.^2./（w.^2））;
s_alpha=w.^4.*s_kesi./（9.8^2）;%波傾角譜
%以下進行波傾角修正
%按照船寬和吃水定kbkt進而求s_alpha_v;
kb=exp（-w.^2.*T./（6.24*pi））;
kt=1-sqrt（Sw）.*B.^2.*w.^4./（1.56^2*16*pi^4）;
s_kesi_v=0.778./（w.^5）.*exp（-0.194375./（w.^4））./（1+2.*w./9.8.*v.*cos（mu*pi/180））;
s_alpha_v=w.^4.*s_kesi_v./（9.8^2）.*kb.^2.*kt.^2;

%figure;hold?on;set（gca‘XLim‘[0.22]）;
%plot（ws_kesi‘r‘ws_alpha‘g--‘）;
s_alpha_ap=sqrt（2.*s_alpha_v.*0.05）;%波傾角幅值
s_alpha_v_ap=sqrt（2.*s_alpha_v.*delta_we）;
%產生0~2*pi之間均勻分布的隨機變量epsilon
y=zeros（251）;y（1）=1977326743;b=7;M=10^10;
for?i=1:25
????y（i+1）=mod（7*y（i）M）;
????x（i）=y（i）/M;
end
epsilon=2*pi*x;

alpha=0;alpha_v=0;
t=0:0.1:30*pi;
for?i=1:25
????a=s_alpha_ap（i）.*sin（w（i）*t+epsilon（i））;
????alpha=alpha+a;
????aa=s_alpha_v_ap（i）.*sin（we（i）*t+epsilon（i））;
????alpha_v=alpha_v+aa;
end

miu_roll=0.1614;%橫搖阻尼，可以變化，0~2；
omg_roll=0.8;%橫搖固有頻率給定？根據船型計算？
A=we./omg_roll;%遭遇頻率和固有頻率之比
W_roll_2=1./（（1-A.^2）.^2+4.*（miu_roll^2）.*（A.^2））;%幅頻響應函數
s_roll=W_roll_2.*s_alpha_v;%橫搖角能譜用s_alpha_v代替s_alpha
%plot（ws_roll）

s_roll_ap=s_alpha_v_ap.*sqrt（W_roll_2）;
epsilon_roll=atan（2.*miu_roll.*A./（1-A.^2））;
hyj=0;
for?i=1:25
????r=s_roll_ap（i）.*sin（we（i）*t-epsilon_roll（i））;
????hyj=hyj+r;
end

hold?on;
figure（1）plot（t10.*alpha‘r‘t10*alpha_v‘b‘）;xlabel（‘t?（s）‘）;ylabel（‘\alpha‘）;??%繪制波傾角時域圖
legend（‘波傾角‘‘有效波傾角‘）;
%plot（t20.*alpha_v）;xlabel（‘t?（s）‘）;ylabel（‘\alpha_v‘）;
figure（2）plot（t10.*hyj）;xlabel（‘t?（s）‘）;ylabel（‘橫搖角‘）;???%繪制橫搖角時域圖

%subplot（313）;
%plot（t20.*hyj）;xlabel（‘t?（s）‘）;ylabel（‘橫搖角‘）;???%繪制橫搖角時域圖

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----

?????文件???????2075??2010-10-14?08:27??roll.m

-----------?---------??----------?-----??----

?????????????????2075????????????????????1