Расчет угла точки впереди между спутниками LEO-GEO (сценарий ISL)

Итак, у меня есть два эталонных TLE от 2 разных спутников в LEO-GEO, которые выглядят следующим образом: (1-й TLE для LEO, 2-й TLE для GEO)

1 44072U 19015A   19265.80540496 -.00000053  00000-0  00000+0 0  9990
2 44072  97.8892 339.4753 0001195  83.2985 276.8367 14.83660044 27382

1 44476U 19049B   19263.72236756 +.00000078 +00000-0 +00000-0 0  9992
2 44476 000.0697 100.7846 0001501 038.3605 175.5638 01.00275593000497

Я использовал SGP4 Orbit Propagator and Integrated (период анализа с 10:00 20 сентября 2019 г. по 10:00 21 сентября 2019 г.) в Matlab и получил вектор состояния орбиты обоих спутников в декартовой координате. А также с помощью этого ответа я рассчитал угол точки вперед и доплеровский сдвиг. И у меня есть это:

введите описание изображения здесь

Но я не уверен, правильно это или неправильно из-за изменения угла?

ОБНОВЛЕНИЕ Я использую λ "=" 1550 н м для расчета доплеровского сдвига. Так что сюжет Δ ф против Т я м е . Я также добавляю свой код в MATLAB; (где оба файла .mat являются вектором состояния rx ry rz vx vy vz)

clc
clear all
close all
format long g
t = 1:86401;
% LEO SATELLITE
load ('LEOPriPosVel.mat')
r1_x = LEOPriPosVel(:,1);   % Inertial Cartesian Coordinate Position X-axis of LEO Sat
r1_y = LEOPriPosVel(:,2);   % Inertial Cartesian Coordinate Position Y-axis of LEO Sat
r1_z = LEOPriPosVel(:,3);   % Inertial Cartesian Coordinate Position Z-axis of LEO Sat
v1_x = LEOPriPosVel(:,4);   % Inertial Cartesian Coordinate Velocity X-axis of LEO Sat
v1_y = LEOPriPosVel(:,5);   % Inertial Cartesian Coordinate Velocity Y-axis of LEO Sat
v1_z = LEOPriPosVel(:,6);   % Inertial Cartesian Coordinate Velocity Z-axis of LEO Sat
%GEO SATELLITE
load ('GEOIn39PosVel.mat')
r2_x = GEOIn39PosVel(:,1);   % Inertial Cartesian Coordinate Position X-axis of GEO Sat
r2_y = GEOIn39PosVel(:,2);   % Inertial Cartesian Coordinate Position Y-axis of GEO Sat
r2_z = GEOIn39PosVel(:,3);   % Inertial Cartesian Coordinate Position Z-axis of GEO Sat
v2_x = GEOIn39PosVel(:,4);   % Inertial Cartesian Coordinate Velocity X-axis of GEO Sat
v2_y = GEOIn39PosVel(:,5);   % Inertial Cartesian Coordinate Velocity Y-axis of GEO Sat
v2_z = GEOIn39PosVel(:,6);   % Inertial Cartesian Coordinate Velocity Z-axis of GEO Sat
for i = 1:86401
    r(i,1) = r1_x(i) - r2_x(i);
    r(i,2) = r1_y(i) - r2_y(i);
    r(i,3) = r1_z(i) - r2_z(i);
    v(i,1) = v1_x(i) - v2_x(i);
    v(i,2) = v1_y(i) - v2_y(i);
    v(i,3) = v1_z(i) - v2_z(i);
    modr12(i) = sqrt((r(i,1)*r(i,1)) + (r(i,2)*r(i,2)) + (r(i,3)*r(i,3)));
    modv12(i) = sqrt((v(i,1)*v(i,1)) + (v(i,2)*v(i,2)) + (v(i,3)*v(i,3)));
    unitvecR(i,1) = r(i,1)/modr12(i);
    unitvecR(i,2) = r(i,2)/modr12(i);
    unitvecR(i,3) = r(i,3)/modr12(i);
    crossVR (i,1) = v(i,2)*unitvecR(i,3) - v(i,3)*unitvecR(i,2);
    crossVR (i,2) = -(v(i,1)*unitvecR(i,3) - v(i,3)*unitvecR(i,1));
    crossVR (i,3) = v(i,1)*unitvecR(i,2) - v(i,2)*unitvecR(i,1);
    dotVR12 (i) = -(v(i,1)*unitvecR(i,1) + v(i,2)*unitvecR(i,2) + v(i,3)*unitvecR(i,3));
    modcrossVR12 (i) = sqrt((crossVR (i,1)*crossVR (i,1)) + (crossVR (i,2)*crossVR (i,2)) + (crossVR (i,3)*crossVR (i,3)));
end
modr = modr12';
modv = modv12';
modcrossVR = modcrossVR12';
dotVR = dotVR12';
for i = 1:86401
    c =  299792.458;
    lambda = 1.55e-9; 
    PAA12(i) = 2*modcrossVR(i)/c;
    CF12(i) = dotVR(i)/lambda;
end
denomin = denom';
PAA = PAA12';
CF = CF12';
figure (1)
subplot(2,1,1)
plot (t,PAA)
title('Changes in Point Ahead Angle $(rad/s)$ LEO-GEO','Interpreter','latex')
xlabel('Time (sec)','Interpreter','latex')
ylabel('Angle (rad)','Interpreter','latex')
subplot(2,1,2)
plot (t,CF)
title('Changes in Frequency $(Hz/s)$ LEO-GEO','Interpreter','latex')
xlabel('Time (sec)','Interpreter','latex')
ylabel('Frequency (Hz)','Interpreter','latex')

Ответы (1)

Частичный ответ. Вот что у меня есть до сих пор. Я использую Python вместо Matlab и «прокручиваю свои собственные» точечные произведения, но эти графики очень похожи на ваши графики! Я думаю, что в выражении для угла может отсутствовать «2», но теперь, когда вы упомянули, что используете свет с длиной волны 1550 нм, мы, похоже, согласны с размером доплеровского сдвига, хотя разница в знаках все еще существует.

Взгляни.

Тем временем я проведу более тщательный численный анализ нескольких конкретных моментов.

введите описание изображения здесь

введите описание изображения здесь

Это Python 3, использующий пакет Skyfield .

TLEs = """1 44072U 19015A   19265.80540496 -.00000053  00000-0  00000+0 0  9990
2 44072  97.8892 339.4753 0001195  83.2985 276.8367 14.83660044 27382
1 44476U 19049B   19263.72236756 +.00000078 +00000-0 +00000-0 0  9992
2 44476 000.0697 100.7846 0001501 038.3605 175.5638 01.00275593000497"""

import numpy as np
import matplotlib.pyplot as plt
from skyfield.api import Topos, Loader, EarthSatellite
from mpl_toolkits.mplot3d import Axes3D

load  = Loader('~/Documents/fishing/SkyData')  # single instance for big files
ts    = load.timescale()
de421 = load('de421.bsp')
earth = de421['earth']

minutes = np.arange(24*60 + 1)
seconds = 60. * minutes
times   = ts.utc(2019, 9, 20, 10, minutes) # starts 09-Sep-2019 10:00 UTC

L0, L1, L2, L3 = TLEs.splitlines()

LEO = EarthSatellite(L0, L1)
GEO = EarthSatellite(L2, L3)

LEOposns = LEO.at(times).position.km   # kilometers
GEOposns = GEO.at(times).position.km

LEOvels  = LEO.at(times).velocity.km_per_s
GEOvels  = GEO.at(times).velocity.km_per_s

if True:
    for i, positions in enumerate((LEOposns, GEOposns)):
        plt.subplot(2, 1, i+1)
        for component in positions:
            plt.plot(seconds, component)
    plt.show()

r    = LEOposns - GEOposns
rhat = r / np.sqrt((r**2).sum(axis=0))

clight = 2.9979E+05  # km/sec
lam    = 1550E-12    # km (1550 nanometers expressed in kilometers)
f      = clight / lam

df_f = -((LEOvels - GEOvels) * rhat).sum(axis=0) / clight
df   = df_f * f
cross = np.cross( (LEOvels - GEOvels).T, rhat.T).T
angle = 2 * np.sqrt((cross**2).sum(axis=0)) / clight

if True:
    fig = plt.figure()
    ax  = fig.add_subplot(2, 1, 1)
    ax.ticklabel_format(style='sci',scilimits=(-3,4),axis='both')
    ax.plot(seconds, angle)
    ax.set_title('Lookahead angle (rads)', fontsize=16)
    ax  = fig.add_subplot(2, 1, 2)
    ax.ticklabel_format(style='sci',scilimits=(-3,4),axis='both')
    ax.plot(seconds, df)
    ax.set_title('Doppler shift Hz (@1550 nm)', fontsize=16)
    plt.show()

if True:
    fig = plt.figure()
    ax  = fig.add_subplot(1, 1, 1, projection='3d')
    x, y, z = LEOposns
    print(x.max())
    ax.plot(x, y, z)
    x, y, z = GEOposns
    print(x.max())
    ax.plot(x, y, z)
    ax.set_xlim(-42000, 42000)
    ax.set_ylim(-42000, 42000)
    ax.set_zlim(-42000, 42000)
    plt.show()
+1 отлично. Ну, я взял Δ ф "=" р е л а т я в е р а д я а л В е л о с я т у / λ , где λ "=" 1550 н м
Можете ли вы также показать мне в том же коде, как мы можем получить вектор орбитального состояния в python из этих TLE?
Я обновил его, с кодами, может быть вы можете посмотреть.
@JOY Хорошо, я сделал обновление, используя 1550 нм, и соглашение улучшилось. Что хорошо в Python, так это то, что вы можете найти пакеты, которые делают практически все, что вам нравится. Здесь я использую пакет Skyfield , и я также добавлю примечание в ответ. Python бесплатный, с открытым исходным кодом и очень, очень интересный. Это все, что я могу сделать сегодня, я проверю завтра. Спасибо за матлаб!
+1 @uhoh Раньше я использовал Python, но потом перешел на Matlab, так что я потерял свой путь в Python, но все же я знаю кое-что о Python, но не хардкорно. Я начну использовать пакет Skyfield в python3. Спасибо за информацию.
Расчет угла точки впереди между спутниками LEO-GEO (сценарий ISL)
СпросиСетьПолитика конфиденциальности1y, 0v