%% 恒力始终垂直于速度的运动轨迹模拟
clear; clc; close all;
hold on;
grid on;
axis equal;
%% 参数设置
m = 2.0;            % 质量 (kg)
v0 = 5.0;           % 初始速度大小 (m/s)
f = 2.0;            % 恒力大小 (N)
theta0 = 0;         % 初始速度方向 (弧度，0表示+x方向)
t_total = 100;       % 总时间 (s)
dt = 0.0001;          % 时间步长 (s)

% 计算理论值
r_theory = m * v0^2 / f;           % 理论半径
omega_theory = f / (m * v0);       % 理论角速度
T_theory = 2*pi / omega_theory;    % 理论周期

fprintf('理论计算值：\n');
fprintf('  半径 r = %.3f m\n', r_theory);
fprintf('  角速度 ω = %.3f rad/s\n', omega_theory);
fprintf('  周期 T = %.3f s\n', T_theory);

%% 数值模拟（欧拉法）
% 初始化数组
n_steps = round(t_total / dt);
t = linspace(0, t_total, n_steps);
x = zeros(1, n_steps);
y = zeros(1, n_steps);
vx = zeros(1, n_steps);
vy = zeros(1, n_steps);

% 初始条件
x(1) = 0;                       % 初始位置
y(1) = 0;
vx(1) = v0 * cos(theta0);       % 初始速度
vy(1) = v0 * sin(theta0);

% 数值积分
for i = 1:n_steps-1
    % 当前速度大小和方向
    v = sqrt(vx(i)^2 + vy(i)^2);
    
    % 计算垂直于速度方向的单位向量
    % 垂直于(vx, vy)的向量可以是(-vy, vx)或(vy, -vx)
    % 这里选择使质点顺时针旋转的方向
    nx = -vy(i) / v;
    ny = vx(i) / v;
    
    % 计算加速度（力垂直于速度）
    ax = (f / m) * nx;
    ay = (f / m) * ny;

    % 更新速度和位置
    vx(i+1) = vx(i) + ax * dt;
    vy(i+1) = vy(i) + ay * dt;
    
    x(i+1) = x(i) + vx(i) * dt;
    y(i+1) = y(i) + vy(i) * dt;
end



plot(x, y, 'b-', 'LineWidth', 1.5);