我找到的-船舶横摇纵摇角仿真程序

acqua_23

后面关于船舶横摇模型求解、船舶纵摇模型求解部分请帮忙注释一下，
[Yt,Xt]=lsim(AAt,BBt,CCt,DDt,at,t);
[Yf,Xf]=lsim(AAf,BBf,CCf,DDf,af,t);
怎么就简单的两个语句就实现了横摇角和纵摇角的仿真曲线了呢？
船舶的传递函数是个2阶震荡环节，下面的程序是如何实现离散化的？
用其他的离散化方法怎么样？谁会教我一下！小生感激不尽……
clear;
%%%%%船舶设计参数.
L=56.5;              %船长L（m）
B=10.4;              %船宽B（m）  
d=3.1;               %吃水d（m）
D=1323.2;            %排水D（m）
Cw=0.686;            %设计水线面积系数
Cb=0.559;            %方形系数
%%%%%船舶横摇参数.
g=9.8;               %重力加速度g
c=0.35;              %计算横摇惯性半径系数
Lou=c*B;             %横摇惯性半径（m）
Jt=D/g*Lou^2;        %横摇总转动惯量Jθ+△Jθ
dJt=0.25*Jt;         %横摇附加转动惯量 △Jθ（kg*m*s^2）
h=0.94;              %横稳性中心高 h（m）
wt=(D*h/Jt)^1/2;     %横摇固有角频率ωθ
mt=0.16/2;           %横摇阻尼因子μθ
%%%%%船舶纵摇参数.
Jf=0.07*D/g*Cw*L^2;  %纵摇总转动惯量Jφ
dJf=Jf;              %纵摇附加转动惯量△Jφ
R=(Cw*L)^2/(14*Cb*d);%纵稳心半径（m）
H=R;                 %纵稳性中心高（m）
wf=(D*H/(Jf+dJf))^1/2;%纵摇固有角频率ωφ
mf=0.36;             %纵摇阻尼因子μφ
%%%%%随机海浪参数.
Ut=11;               %风速（m/s）（22Kn）
Lm=128;              %海浪波长（m）
h1=3.3;              %有义波高（m）
T=8.6;               %波浪固有周期（s）
PI=3.1415926;
w=2*PI/T;            %波圆频率ω
v=Lm/T;              %波相速度（m/s）
k=2*PI/Lm;           %波数
hw=0.17*Lm^(3/4);    %波高（m）
waveA=hw/2;          %波幅（m）
Xt=exp(-k*d/2);      %有限吃水修正系数Xθt
Xb=1-sqrt(Cw)*(B/Lm)^2;%有限船宽修正系数Xθb
x=Xt*Xb;             %波倾修正系数Xθ
a0=k*waveA;          %最大波倾角a0
ae0=x*a0;            %最大有效波倾角ae0
%%%%%随机海浪的波倾角仿真.
beta=PI/6;           %遭遇波倾角β（浪向角）
V=5;                 %船舶航速V（m/s）
dw=0.08;             %仿真频率增量△ω
wi=0.24:0.08:2.4;    %各仿真点频率
a=0.779;             %海浪谱系数A
b=3.116/h1^2;        %海浪谱系数B  
for i=1:28
    we(i)=wi(i)-wi(i)^2/g*V*cos(beta);%各仿真点遭遇频率ωe
    Swi(i)=a/wi(i)^5*exp(-b/wi(i)^4); %单参数海浪谱S（ω）
    Sawi(i)=wi(i)^4/g^2*Swi(i);       %波倾角能量谱密度Sa（ω）
    Sadw(i)=sqrt(2*Sawi(i)*dw);      
end

MinT=0;             %仿真采样的最小时间
MaxT=MinT+600;      %仿真采样的最大时间
t=MinT:0.2:MaxT;    %时间t（s）
at=0;               %横摇波倾角
af=0;               %纵摇波倾角
for i=1:28
    e(i)=rand;      %产生随机数εi
    at=at+(Sadw(i)*(cos(we(i)*t+e(i))))*sin(beta);
    af=af+(Sadw(i)*(cos(we(i)*t+e(i))))*cos(beta);
end

%%%%%船舶横摇模型求解.
%状态方程:[DX1;DX2]=[0 1;-a2 -a1]*[X1;X2]+[b1;b2]*u,即;X'=AX+BU
%输出方程;y=[1 0]*[X1;X2],即:Y=CX+DU
%X1,X2为状态变量,输入信号 u=a(t),输出信号 y=θ
a1t=2*mt*wt;
a2t=wt^2;
b0t=wt^2*dJt/(D*h);
b1t=a1t*(1-b0t);
b2t=wt^2-a1t^2*(1-b0t)-wt^2*b0t;
AAt=[0 1;-a2t -a1t];
BBt=[b1t;b2t];
CCt=[1 0];
DDt=b0t;
[Yt,Xt]=lsim(AAt,BBt,CCt,DDt,at,t);%得到船舶摇荡的横摇角信号
%%%%%船舶纵摇模型
a1f=2*mf*wf;
a2f=wf^2;
b0f=wf^2*dJf/(D*H);
b1f=a1f*(1-b0f);
b2f=wf^2-a1f^2*(1-b0f)-wf^2*b0f;
AAf=[0 1;-a2f -a1f];
BBf=[b1f;b2f];
CCf=[1 0];
DDf=b0f;
[Yf,Xf]=lsim(AAf,BBf,CCf,DDf,af,t);%得到船舶摇荡的纵摇角信号
%%%%%船舶横摇角.纵摇角仿真曲线
figure(1)
plot(t,Yt,'-r',t,0)
xlabel('时间（秒）','FontWeight','bold')
ylabel('船舶横摇角','FontWeight','bold')
grid on
figure(2)
plot(t,Yf,'-k',t,0)
xlabel('时间（秒）','FontWeight','bold')
ylabel('船舶纵摇角','FontWeight','bold')
grid on

czg0915

有没有舵角的仿真程序啊！！急用！！

stream

好东西，先支持。