%升力體：通用航空運載器CAV-H
%讀取數據1 升力系數
alpha = [10 15 20];   Ma = [3.5 5 8 10 15 20 23];
alpha1 = 10:0.1:20;   Ma1 = 3.5:0.1:23; 
[Ma1, alpha1] = meshgrid(Ma1, alpha1);
CL = readmatrix('simulation.xlsx', 'Sheet', 'Sheet1', 'Range', 'B2:H4');
CL1 = interp2(Ma, alpha, CL, Ma1, alpha1, 'spline');
CD = readmatrix('simulation.xlsx', 'Sheet', 'Sheet1', 'Range', 'B7:H9');
CD1 = interp2(Ma, alpha, CD, Ma1, alpha1, 'spline');% 常量
Sm = 0.484;     % 參考面積750in2
m = 907.18;     % 質量
g = 9.80665;    % 重力加速度%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 求解區 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 初始滑翔段
% 初始值：[Vx, Vy, Vz, x, y, z]
% 初始滑翔段初始值（列向量）
y0_1 = [6500; 0; 0; 0; 0; 80000];  % 修改為6x1列向量 % Vx=6500m/s，初始高度80km
tspan1 = 0:0.02:300;
options = odeset('Events', @odeEventFun);
[t1, y1] = ode45(@(t,y) Fun1(t, y, m, Sm, CL1, CD1, g), tspan1, y0_1, options);% 滑翔段
y0_2 = y1(end, :);  % 繼承初始段終態
tspan2 = t1(end):0.02:2000;
[t2, y2] = ode45(@(t,y) Fun2(t, y, m, Sm, CL1, CD1, g), tspan2, y0_2);% 下壓段
y0_3 = y2(end, :);  % 繼承滑翔段終態
tspan3 = 2000:0.02:2300;
options1 = odeset('Events', @odeEventFun1);
[t3, y3] = ode45(@(t,y) Fun3(t, y, m, Sm, CL1, CD1, g), tspan3, y0_3, options1);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 結果可視化 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 合并軌跡數據
t = [t1; t2; t3];
y = [y1; y2; y3];
Vx = y(:,1); Vy = y(:,2); Vz = y(:,3);
x = y(:,4); y_pos = y(:,5); z = y(:,6);
V = sqrt(Vx.^2 + Vy.^2 + Vz.^2);% 繪制高度-時間曲線
figure(1)
plot(t, z