就以下代码写一篇跟有心力相关的论文

%% 系统参数
N = 50000;

%% 写avi文件
elec_atom = VideoWriter('elec_atom'); % 指定视频文件名
open(elec_atom) % open开始写入


%% 物理参数
dt = 0.01;
ef = 8.854187817 * 10e-12; % 真空介电常数
e = 1.6 * 10e-19; % 电子的电荷量
k = 1 / (4 * pi * ef); % k
Vc = 20; % 光速
Rk = 10e-14;


%% 电子的参数
elec_pos = [-200, 0];
elec_speed = [0, 5];
elec_mess = 1;
elec_v = 0;
elec_m0 = 1;
elec_Q = e; % 一个负电子
elec_pos_x = zeros(1, N);
elec_pos_y = zeros(1, N);
elec_v_box = zeros(1, N);


%% 原子的参数
atom_pos = [0, 0]; % 位置
atom_speed = [0, 0]; % 速度
atom_mess = 1836; % 质量
atom_v = 0; % 速度


for i = 1:1:N
DeltaX = (atom_pos(1) - elec_pos(1)); % 计算出两个坐标的坐标差值矢量
DeltaY = (atom_pos(2) - elec_pos(2));
RadiusSquart = DeltaY ^ 2 + (DeltaX) ^ 2; % 计算两个原子的距离
DeltaPos = sqrt(RadiusSquart);
a1 = k * elec_Q / RadiusSquart / Rk / elec_mess;
sin_theta = DeltaY / DeltaPos;
cos_theta = DeltaX / DeltaPos; 
elec_speed = [elec_speed(1) + a1 * cos_theta * dt, elec_speed(2) + a1 * sin_theta * dt];
elec_pos = [elec_pos(1) + elec_speed(1) * dt, elec_pos(2) + elec_speed(2) * dt];
elec_v = sqrt(elec_speed(1) ^ 2 + elec_speed(2) ^ 2);
elec_mess = elec_m0 / (1 - abs(elec_v / Vc) ^ 2) ^ 1.5; % 质量随着电子的速度的变化而变化
%% 记录数据
elec_v_box(i) = elec_v; % 记录下每次迭代的速度
elec_pos_x(i) = elec_pos(1);
elec_pos_y(i) = elec_pos(2);
end





for i = 1:100:N 
scatter(0, 0, 'r', 'LineWidth', 3)
hold on 
xlim([-300, 300])
ylim([-300, 300])
plot(elec_pos_x(1:i), elec_pos_y(1:i), 'g', 'LineWidth', 1) % 绘制轨迹1
legend('氢原子', '电子')
Frame = getframe(gcf); % 抓取图窗
writeVideo(elec_atom, Frame); % 写入文
clf
end
close(elec_atom)

clc
clear 


%% 物理常量
G = 6.67e-11; % 万有引力常数
dt = 0.01; % 时间间隔
N = 10000; % 迭代次数
Rk = 10e-13; % 缩放


%% 地球的参数
star1_pos = [-200, 0]; % 位置
star1_mess = 200; % 质量
star1_vx = 0; % 速度x
star1_vy = 25; % 速度y
star1_v = star1_vy; % 速度
star1_pos_x = zeros(1, N);
star1_pos_y = zeros(1, N);
star1_v_box = zeros(1, N);


%% 月亮的参数
moon_pos = [-230, 0]; % 位置
moon_mess = 1; % 质量

moon_vx = 0; % 速度x
moon_vy = 35; % 速度y
moon_v = moon_vy; % 速度
moon_pos_x = zeros(1, N);
moon_pos_y = zeros(1, N);
moon_v_box = zeros(1, N);


sun_pos = [0, 0];
sun_mess = 1836;
sun_v = 0;
sun_vx = 0;
sun_vy = 0;


%% 写gif文件
ti3 = VideoWriter('3ti'); % 指定视频文件名
open(ti3) % open开始写入



%% 绘制
for i = 1 : N
DeltaX = (sun_pos(1) - star1_pos(1)); % x的距离
DeltaY = (sun_pos(2) - star1_pos(2));
RadiusSquart = DeltaY ^ 2 + (DeltaX) ^ 2;
DeltaPos = sqrt(RadiusSquart);
a1 = G * sun_mess / RadiusSquart / Rk ; % 当前位置的加速度
sin_theta = DeltaY / DeltaPos; % 夹角
cos_theta = DeltaX / DeltaPos; 
star1_vx = star1_vx + a1 * cos_theta * dt;
star1_vy = star1_vy + a1 * sin_theta * dt;
star1_v = sqrt(star1_vx^2 + star1_vy^2); 
star1_pos(1) = star1_pos(1) + star1_vx * dt;
star1_pos(2) = star1_pos(2) + star1_vy * dt;

% 记录数据
star1_v_box(i) = star1_v;
star1_pos_x(i) = star1_pos(1);
star1_pos_y(i) = star1_pos(2);


%% 月亮和地球
DeltaX = (star1_pos(1) - moon_pos(1)); % x的距离
DeltaY = (star1_pos(2) - moon_pos(2));
RadiusSquart = DeltaY ^ 2 + (DeltaX) ^ 2;
DeltaPos = sqrt(RadiusSquart);
a2 = G * star1_mess / RadiusSquart / Rk; % 当前位置的加速度
sin_theta = DeltaY / DeltaPos; % 夹角
cos_theta = DeltaX / DeltaPos; 
moon_vx = moon_vx + a2 * cos_theta * dt;
moon_vy = moon_vy + a2 * sin_theta * dt;
moon_v = sqrt(moon_vx^2 + moon_vy^2); 
moon_pos(1) = moon_pos(1) + moon_vx * dt;
moon_pos(2) = moon_pos(2) + moon_vy * dt;

% 记录数据
moon_v_box(i) = moon_v;
moon_pos_x(i) = moon_pos(1);
moon_pos_y(i) = moon_pos(2);

end 


% 绘制动态图
DrawSun() % 绘制太阳
hold on

% for i = 1:10:N
% plot(star1_pos_x(1:i), star1_pos_y(1:i), 'g', 'LineWidth', 4) % 绘制轨迹
% frame = getframe(gcf); % 获取当前的关键帧
% I = frame2im(frame); % 将关键帧转化成图像
% [I, map] = rgb2ind(I, 256);
% imwrite(I, map, 'test.gif', 'gif', 'LoopCount', inf, 'DelayTime', 0.2); %将关键帧保存到gif中
% end
disp(1)

for i = 1:50:N
% scatter(star1_pos_x(i), star1_pos_y(i), 'g', 'LineWidth', 1) % 绘制轨迹
% scatter(moon_pos_x(i), moon_pos_y(i), 'b', 'LineWidth', 1) % 绘制轨迹
plot(star1_pos_x(1:i), star1_pos_y(1:i), 'g', 'LineWidth', 3) % 绘制轨迹1
plot(moon_pos_x(1:i), moon_pos_y(1:i), 'b', 'LineWidth', 2) % 绘制轨迹
Frame = getframe(gcf); % 抓取图窗
writeVideo(ti3, Frame); % 写入文
end
close(ti3)


%% 绘图
function out = DrawSun()
out = 0;
scatter(0, 0, 'r', 'LineWidth', 1)
xlim([-300, 300])
ylim([-300, 300])
hold on 
end

clf;
clear;


% 设定模拟参数
m = 0.2; % 小球质量
k = 0.1; % 弹簧劲度系数
r = 1; % 弹簧原长
repulsion = 0.1; % 斥力系数，用于控制小球之间的斥力大小


% 初始化小球位置和速度
x1 = -1; % 小球1的初始位置
x2 = 1; % 小球2的初始位置
v1 = -0.2; % 小球1的初始速度
v2 = -0.1; % 小球2的初始速度


% 设置动画参数
dt = 0.1; % 时间步长
t = 0; % 初始时间
duration = 15; % 模拟时长


% 创建图形窗口
% figure('Position', [100, 100, 800, 400]);


% 开始模拟循环
while t < duration
% 计算弹簧力
force = k * (x2 - x1 - r) - repulsion / ((x2 - x1)^2 + 0.01); % 斥力项用于防止小球碰撞
% 计算加速度
a1 = force / m;
a2 = -force / m;
% 更新速度和位置
v1 = v1 + a1 * dt;
v2 = v2 + a2 * dt;
x1 = x1 + v1 * dt;
x2 = x2 + v2 * dt;
% 绘制小球和弹簧
% plot([x1, x2], [0, 0], 'k', 'LineWidth', 2); % 绘制弹簧
draw_tan(x1, x2)
hold on;
plot(x1, 0, 'ro', 'MarkerSize', 20, 'MarkerFaceColor', 'r'); % 绘制小球1
plot(x2, 0, 'bo', 'MarkerSize', 20, 'MarkerFaceColor', 'b'); % 绘制小球2
hold off;
% 设置坐标轴范围
xlim([-3, 3]);
ylim([-1, 1]);
% 设置图形标题和时间显示
title(['Time: ', num2str(t)]);
% 更新图形窗口
drawnow;
% 增加时间步长
t = t + dt;
end


function out = draw_tan(x1, x2)
out = 0;
n = 15;
tempx = linspace(x1, x2, n);
tempy = zeros(1, length(tempx));
n = length(tempx);
for i = 1:1:n
if mod(i, 2) == 0
tempy(i) = 0.2;
else
tempy(i) = -0.2;
end 
end
tempy(1) = 0;
tempy(2) = 0;
tempy(3) = 0;
tempy(end) = 0;
tempy(n - 1) = 0;
tempy(n - 2) = 0;
plot(tempx, tempy, 'g','LineWidth', 2)
hold on 
end