%% 带阻力（f=kv）的平抛运动模拟
close
clc
clear
hold on
grid on
axis equal;


%% 参数设置
m = 1;           % 质量 (kg)
k = 0.98;          % 阻力系数 (kg/s)
g = 9.8;           % 重力加速度 (m/s^2)
v0 = 10;           % 初始速度 (m/s)
h0 = 20;          % 初始高度 (m)



% 初始位置
x0 = 0;
y0 = h0;

%% 时间设置
dt = 0.0001;         % 时间步长 (s)
t_max = 15;        % 最大模拟时间 (s)
t = 0:dt:t_max;    % 时间向量
N = length(t);

%% 初始化数组
x = [];
y = [];
vy = -v0; %垂直方向做匀速运动，整个过程速度不变
v_tilt = [];
a_tilt = []; %倾斜方向


% 设置初始条件
x(1) = x0;
y(1) = y0;

v_tilt(1) =sqrt(2) * v0 ;

%% 数值求解（欧拉法）
s_theta=atan(-vy/v0);
for i = 1:N-1
    
    
    % 计算加速度（考虑阻力）
    a_tilt(i)=-(k/m) * v_tilt(i);
    
    % 更新速度和位置
    v_tilt(i+1) = v_tilt(i) + a_tilt(i) * dt;
    s_tilt=v_tilt(i)*dt;%倾斜方向运行的距离，需要分析成水平与垂直方向

    x(i+1) = x(i) + s_tilt * sqrt(1/2);
    y(i+1) = y(i) + vy*dt + s_tilt * sqrt(1/2);


   if length(x) >=2
    %45度退出
    x_end=x(end)- x(end-1);
    y_end=y(end)- y(end-1);
    d=abs(atand (y_end/x_end));
    if d >=45
        fprintf('到达45:\n');
         v_a=[0,vy];
        v_b=[v_tilt(i) * cos(s_theta),v_tilt(i) * sin(s_theta)];
        v_c=v_a+v_b;
        fprintf('当前速度:%.2f\n',norm(v_c));
        % fprintf('当前速度:%.2f\n',sqrt(vy^2 +v_tilt(i)^2 -
        % 2*vy*v_tilt(i)*cos(3*pi/4) )); %有问题不能按余弦定理计算
        break;
    end
   end

end

plot(x, y, 'b-', 'LineWidth', 1);
fprintf("计算速度:%.2f\n",sqrt(1/2)*v0);


%% 无阻力
x_n=[x0];
y_n=[y0];
vx_n=v0;
vy_n=0;

i=0;
while y_n(end)>=0
 i=i+1;
 vx_n=v0;
 vy_n=vy_n - g*dt;
 x_n(i+1)=x_n(i)+vx_n*dt;
 y_n(i+1)=y_n(i)+vy_n*dt;

 if(y_n(end)<= y(end))
    break;
 end
end
plot(x_n, y_n, 'r-', 'LineWidth', 1);

%% 带阻力（f=kv）的平抛运动模拟
close
clc
clear
hold on
grid on
axis equal;


%% 参数设置
m = 1;           % 质量 (kg)
k = 0.98;          % 阻力系数 (kg/s)
g = 9.8;           % 重力加速度 (m/s^2)
v0 = 10;           % 初始速度 (m/s)
theta = 0;         % 平抛，角度为0度
h0 = 20;          % 初始高度 (m)

% 将角度转换为弧度
theta_rad = deg2rad(theta);

% 初始速度分量
vx0 = v0 * cos(theta_rad);
vy0 = v0 * sin(theta_rad);

% 初始位置
x0 = 0;
y0 = h0;

%% 时间设置
dt = 0.0001;         % 时间步长 (s)
t_max = 15;        % 最大模拟时间 (s)
t = 0:dt:t_max;    % 时间向量
N = length(t);

%% 初始化数组
x = [];
y = [];
vx = [];
vy = [];
ax = [];
ay = [];

% 设置初始条件
x(1) = x0;
y(1) = y0;
vx(1) = vx0;
vy(1) = vy0;

%% 数值求解（欧拉法）
for i = 1:N-1

    
    % 计算加速度（考虑阻力）
    ax(i) = -(k/m) * vx(i);           % x方向：只有阻力
    ay(i) = -g - (k/m) * vy(i);       % y方向：重力 + 阻力,vy内容是负的，所以这里取负号
    
    % 更新速度和位置
    vx(i+1) = vx(i) + ax(i) * dt;
    vy(i+1) = vy(i) + ay(i) * dt;
    x(i+1) = x(i) + vx(i) * dt;
    y(i+1) = y(i) + vy(i) * dt;


   if length(x) >=2
    %45度退出
    x_end=x(end)- x(end-1);
    y_end=y(end)- y(end-1);
    d=abs(atand (y_end/x_end));
    if d >=45
        fprintf('到达45:\n');
        fprintf('当前速度:%.2f\n',sqrt(vx(i)^2 +vy(i)^2));
         break;
    end
   end

end

plot(x, y, 'b-', 'LineWidth', 1);
fprintf("计算速度:%.2f\n",sqrt(1/2)*v0);


%% 无阻力
x_n=[x0];
y_n=[y0];
vx_n=v0;
vy_n=0;

i=0;
while y_n(end)>=0
 i=i+1;
 vx_n=v0;
 vy_n=vy_n - g*dt;
 x_n(i+1)=x_n(i)+vx_n*dt;
 y_n(i+1)=y_n(i)+vy_n*dt;

 if(y_n(end)<= y(end))
    break;
 end
end
plot(x_n, y_n, 'r-', 'LineWidth', 1);