%升力体：通用航空运载器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;     % 参考面积750in²
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