ここのコードを改変した。
butterfly-effect.hatenablog.com
clear;cla global A; %有次元安定微係数 Xu = -0.01; Zu = -0.1; Mu = 0.001; Xa = 30; Za = -200; Ma = -4; Xq = 0.3; Zq = -5; Mq = -1; Yb = -45; Lb_= -2; Nb_= 1; Yp = 0.5; Lp_= -1; Np_= -0.1; Yr = 3; Lr_= 0.2; Nr_=-0.2; %その他のパラメタ W0 = 0; U0 = 100; theta0 = 0.05; g = 9.8; %重力加速度 %縦のシステム A_lat =[Xu, Xa, -W0, -g*cos(theta0); Zu/U0, Za/U0, (U0+Zq)/U0, -g*sin(theta0)/U0; Mu, Ma, Mq, 0; 0, 0, 1, 0]; %横・方向のシステム A_lon = [Yb, (W0+Yp), -(U0-Yr), g*cos(theta0), 0; Lb_, Lp_, Lr_, 0, 0; Nb_, Np_, Nr_, 0, 0; 0, 1, tan(theta0), 0, 0; 0, 0, sec(theta0), 0, 0]; %対角ブロックとしてシステムを結合する A = blkdiag(A_lat,A_lon); %計算条件の設定 endurance = 100;%飛行時間[sec] step = 10;%1.0[sec]あたりの時間ステップ数 t = linspace(0,endurance,endurance*step); %初期値 x0 = [u,alpha,q,theta, beta,p,r,phi,psi] x0_lat = [10;0.1;0.4;0.2]; %縦の初期値 x0_lon = [0.0;0.6;0.4;0.2;0.2]; %横・方向の初期値 x0 = vertcat(x0_lat,x0_lon); %運動方程式を解く [t,x]=ode45(@dynamical_system,t,x0); % x = [u;α;q;θ; β;p;r;φ;ψ] u=x(:,1); alpha=x(:,2); beta=x(:,5); v=u.*beta; w=u.*alpha; phi=x(:,8); theta=x(:,4); psi=x(:,9); dat=[t,u,v,w,phi,theta,psi]; %Matlab Aerospace toolbox によるシミュレーション h = Aero.Animation; h.createBody('pa24-250_orange.ac','Ac3d'); h.Bodies{1}.TimeSeriesSource = dat; h.Camera.PositionFcn = @staticCameraPosition; h.Figure.Position(1) = h.Figure.Position(1) + 572/2; h.show(); h.play(); %運動方程式 function dx = dynamical_system(t,x) % x = [u;α;q;θ; β;p;r;φ;ψ] global A; dx = A*x; end
