函数文件
function ydot=ptydfun(t,y,flag,lambda,omega)
global g
ydot=[y(2);2*omega*y(4)*sin(lambda);-2*omega*(y(6)*cos(lambda)+y(2)*sin(lambda));y(6);2*omega*y(4)*cos(lambda)-g]
主程序
clear;clc
global g;g=9.8;
omega=[0,2*pi/(60*60*24)];%不考虑地球自转角速度与考虑地球自转角速度
lambda=[60*2*pi/360, 60*2*pi/360];%北纬60度位置
figure
for i=1:2%循环解两遍微分方程
[t,y]=ode45(‘ptydfun’,[0:0.0001:6.12],[0,0,0,5,0,30],[],lambda(i),omega(i));
Z(i)=max(y(:,5))%求轨道最高点
T(i)=t(find(y(:,5)==H))%到最高点所需要的时间
Y(i)= max(y(:,3))%射程
X(i)= max(y(:,1))%偏南北射程
axis([0 0.02 0 30 0 50]);%设置坐标轴范围
title(‘考虑地球自转的斜抛运动’)；
xlabel(‘x’);
ylabel(‘y’);
zlabel(‘z’);
hold on
plot3(y(:,1), y(:,3), y(:,5))%画出轨迹
grid on
view(70,25)
end
function ydot=ptydfun(t,y,flag,lambda,omega)
global g
ydot=[y(2);2*omega*y(4)*sin(lambda);-2*omega*(y(6)*cos(lambda)+y(2)*sin(lambda));y(6);2*omega*y(4)*cos(lambda)-g]
主程序
clear;clc
global g;g=9.8;
omega=[0,2*pi/(60*60*24)];%不考虑地球自转角速度与考虑地球自转角速度
lambda=[60*2*pi/360, 60*2*pi/360];%北纬60度位置
figure
for i=1:2%循环解两遍微分方程
[t,y]=ode45(‘ptydfun’,[0:0.0001:6.12],[0,0,0,5,0,30],[],lambda(i),omega(i));
Z(i)=max(y(:,5))%求轨道最高点
T(i)=t(find(y(:,5)==H))%到最高点所需要的时间
Y(i)= max(y(:,3))%射程
X(i)= max(y(:,1))%偏南北射程
axis([0 0.02 0 30 0 50]);%设置坐标轴范围
title(‘考虑地球自转的斜抛运动’)；
xlabel(‘x’);
ylabel(‘y’);
zlabel(‘z’);
hold on
plot3(y(:,1), y(:,3), y(:,5))%画出轨迹
grid on
view(70,25)
end
