Matlab

Question

我必须绘制围绕中心物体运行的物体的速度矢量。这是一个开普勒语境。物体的轨迹是从经典公式（r = p /（1 + e * cos（theta））推导出的，其中e =偏心率。

我设法绘制了椭圆轨道，但现在，我想为该轨道的每个点绘制物体的速度。

为了计算速度矢量，我从经典公式（到极坐标）开始，在 2 个分量下面：

v_r = dr/dt 和 v_theta = r d(theta)/dt

为了采用时间步长 dt，我提取了与时间成正比的平均异常。

最后，我计算了这个速度向量的归一化。

clear             % clear variables

e = 0.8;    % eccentricity
a = 5;             % semi-major axis
b = a*sqrt(1-e^2); % semi-minor axis
P = 10             % Orbital period
N = 200;           % number of points defining orbit

nTerms = 10;    % number of terms to keep in infinite series defining
                % eccentric anomaly

M = linspace(0,2*pi,N);   % mean anomaly parameterizes time
                          % M varies from 0 to 2*pi over one orbit

alpha = zeros(1,N);       % preallocate space for eccentric anomaly array



%%%%%%%%%%
%%%%%%%%%%  Calculations & Plotting
%%%%%%%%%%


% Calculate eccentric anomaly at each point in orbit
for j = 1:N
    % initialize eccentric anomaly to mean anomaly
    alpha(j) = M(j);

    % include first nTerms in infinite series
    for n = 1:nTerms
        alpha(j) = alpha(j) + 2 / n * besselj(n,n*e) .* sin(n*M(j));
    end
end

% calcualte polar coordiantes (theta, r) from eccentric anomaly
theta = 2 * atan(sqrt((1+e)/(1-e)) * tan(alpha/2));
r = a * (1-e^2) ./ (1 + e*cos(theta));

% Compute cartesian coordinates with x shifted since focus
x = a*e + r.*cos(theta);
y = r.*sin(theta);
figure(1);
plot(x,y,'b-','LineWidth',2)
xlim([-1.2*a,1.2*a]);
ylim([-1.2*a,1.2*a]);
hold on;
% Plot 2 focus = foci
plot(a*e,0,'ro','MarkerSize',10,'MarkerFaceColor','r');
hold on;
plot(-a*e,0,'ro','MarkerSize',10,'MarkerFaceColor','r');

% compute velocity vectors
for i = 1:N-1
  vr(i) = (r(i+1)-r(i))/(P*(M(i+1)-M(i))/(2*pi));
  vtheta(i) = r(i)*(theta(i+1)-theta(i))/(P*(M(i+1)-M(i))/(2*pi));
  vrNorm(i) = vr(i)/norm([vr(i),vtheta(i)],1);
  vthetaNorm(i) = vtheta(i)/norm([vr(i),vtheta(i)],1);
end 

% Plot velocity vector
quiver(x(30),y(30),vrNorm(30),vthetaNorm(30),'LineWidth',2,'MaxHeadSize',1);

% Label plot with eccentricity
title(['Elliptical Orbit with e = ' sprintf('%.2f',e)]);

不幸的是，一旦执行绘图，我似乎得到了一个速度不佳的向量。例如，vrNorm 和 vthetaNorm 数组的 30th 元素：

如您所见，向量的方向错误（如果我假设从右轴取 0 的 theta 和三角函数中的正变化）。

如果有人能看到我的错误在哪里，那就太好了。

更新 1：代表椭圆轨道速度的矢量是否与椭圆曲线永久相切？

我想以正确的焦点为原点来表现它。

更新 2：

根据@MadPhysicist的解决方案，我修改了：

% compute velocity vectors
  vr(1:N-1) = (2*pi).*diff(r)./(P.*diff(M));
  vtheta(1:N-1) = (2*pi).*r(1:N-1).*diff(theta)./(P.*diff(M));

  % Plot velocity vector
  for l = 1:9    quiver(x(20*l),y(20*l),vr(20*l)*cos(vtheta(20*l)),vr(20*l)*sin(vtheta(20*l)),'LineWidth',2,'MaxHeadSize',1);
  end
  % Label plot with eccentricity
  title(['Elliptical Orbit with e = ' sprintf('%.2f',e)]);

我得到以下结果：

在轨道的某些部分，我得到了错误的方向，我不明白为什么......

Answer 1

您的代码有两个问题：

规范化不正确。 norm 计算向量的广义 p-范数，默认为欧几里得范数。它需要笛卡尔输入。将 p 设置为 1 意味着它将只是 return 向量的最大元素。在您的情况下，规范化是没有意义的。只需将 vrNorm 设置为
```
vrNorm = vr ./ max(vr)
```
您似乎将极坐标 vrNorm 和 vthetaNorm 传递给 quiver，这需要笛卡尔坐标。以矢量化的方式进行转换很容易：
```
vxNorm = vrNorm * cos(vtheta);
vyNorm = vrNorm * sin(vtheta);
```

这假设我理解你的角度是从哪里来的，并且 vtheta 是以弧度为单位的。

备注

整个循环

for i = 1:N-1
    vr(i) = (r(i+1)-r(i))/(P*(M(i+1)-M(i))/(2*pi));
    vtheta(i) = r(i)*(theta(i+1)-theta(i))/(P*(M(i+1)-M(i))/(2*pi));
    vrNorm(i) = vr(i)/norm([vr(i),vtheta(i)],1);
    vthetaNorm(i) = vtheta(i)/norm([vr(i),vtheta(i)],1);
end

可以完全向量化的方式重写：

vr = (2 * pi) .* diff(r) ./ (P .* diff(M))
vtheta = (2 * pi) .* r .* diff(theta) ./ (P .* diff(M))
vrNorm = vr ./ max(vr)
vxNorm = vrNorm * cos(vtheta);
vyNorm = vrNorm * sin(vtheta);

注2

您可以在整个数据集或子集上以矢量化方式调用 quiver：

quiver(x(20:199:20), y(20:199:20), vxNorm(20:199:20), vyNorm(20:199:20), ...)

Matlab - 在开普勒轨道上绘制卫星的速度矢量

Matlab - plot speed vector of a satellite in Keplerian orbit

vector

orbital-mechanics