Matlab - график вектора скорости спутника на кеплеровской орбите - PullRequest
Новая
       40

Matlab - график вектора скорости спутника на кеплеровской орбите

0 голосов
/ 20 января 2019

Я должен построить вектор скорости объекта, вращающегося вокруг центрального тела.Это кеплеровский контекст.Траектория объекта выводится из классической формулы (r = p / (1 + e * cos (theta)) с e = эксцентриситетом.

Мне удается построить эллиптическую орбиту, но сейчас я бы хотелнанесите для каждой точки этой орбиты скорость скорости объекта.

Чтобы вычислить вектор скорости, я начну с классических формул (в полярные координаты), ниже двух компонентов:

v_r = dr / dt и v_theta = rd (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)]);

К сожалению, после выполнения графика кажется, что я получил плохой вектор для скорости. Вот, например, элемент 30th из vrNorm и vthetaNorm массивы: bad direction

Как вы можете видеть, вектор имеет неправильное направление (если я предполагаю взять 0 для тета от правой оси и положительное отклонение, как в тригонометрии).

Если бы кто-то мог видетьЭто моя ошибка, это было бы неплохо.

ОБНОВЛЕНИЕ 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)]);

Я получаю следующий результат:

enter image description here

На некоторых участках орбиты я получаю неправильные указания и не понимаю, почему ...

1 Ответ

0 голосов
/ 20 января 2019

С вашим кодом есть две проблемы:

  1. Нормализация выполнена неправильно. norm вычисляет обобщенную p-норму для вектора, которая по умолчанию соответствует евклидовой норме. Он ожидает декартовых входов. Установка p в 1 означает, что он просто вернет самый большой элемент вашего вектора. В вашем случае нормализация бессмысленна. Просто установите vrNorm как 

    vrNorm = vr ./ max(vr)

  2. Похоже, что вы передаете в полярных координатах 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), ...)
...