Engee 文档
Notebook

我们正在模拟地球的第一颗人造卫星

在这个项目中，我们将建造1957年10月4日发射的PS-1卫星的地面轨道。

筹备及推出

让我们准备工作：打开海岸等高线图并为其设置所需的投影。

In [ ]: 
Pkg.add(["SatelliteToolbox", "Shapefile", "Proj"])

我们提前从[自然地球]存储库下载了地图（https://www.naturalearthdata.com /)。 你会发现更广泛的信息与地图工作这里。

In [ ]: 
using Shapefile, Downloads, Proj

shp_path = joinpath(@__DIR__, "map", "ne_110m_coastline.shp");
table = Shapefile.Table(shp_path);

trans = Proj.Transformation("+proj=longlat", "+proj=robin +lon_0=0 +x_0=0 +y_0=0");

现在让我们运行模型。 我们将得到几个状态向量。

image.png
In [ ]: 
engee.open("sputnik_model.engee") # 如果模型已经打开，请评论。
data = engee.run("sputnik_model")

X_icrf = collect(data["X_icrf"]).value;
V_icrf = collect(data["V_icrf"]).value;
t_utc_jd = collect(data["t_utc_jd"]).value;

df = DataFrame( time = data["X_icrf"].time, t_utc_jd = t_utc_jd, X_icrf = X_icrf, V_icrf = V_icrf);

我们已经将向量组装成一个DataFrame对象，现在我们将用其他参数填充它。

将坐标转换到其他系统

为了显示图形，将坐标转换为地心和大地测量系统通常很有趣。 请考虑以下注意事项:

-对于我们的项目，ICRF系统几乎与GCRF（ECI）系统相同，误差为一米

-观察员将位于莫斯科

-使用WGS84椭球绘制地心坐标

In [ ]: 
using DataFrames
using SatelliteToolbox
using LinearAlgebra

# 观察者在哪里
OBS_LLA = [deg2rad(55.7558), deg2rad(37.6173), 0.2e3]
eop = fetch_iers_eop()

# 添加新列
df[!, :lon] = zeros(length(df[!, :time]))
df[!, :lat] = zeros(length(df[!, :time])) 
df[!, :alt] = zeros(length(df[!, :time]))
df[!, :az] = zeros(length(df[!, :time]))
df[!, :el] = zeros(length(df[!, :time]))
df[!, :range] = zeros(length(df[!, :time]))

# 地球的参数，用于转换观察者的坐标
a = 6378137.0
lat_obs, lon_obs, alt_obs = OBS_LLA
x_obs = (a + alt_obs) * cos(lat_obs) * cos(lon_obs)
y_obs = (a + alt_obs) * cos(lat_obs) * sin(lon_obs)
z_obs = (a + alt_obs) * sin(lat_obs)

for i in 1:length(df[!, :time])
    # ICRF→ITRF转换
    sv_icrf = OrbitStateVector(df[i, :t_utc_jd], df[i, :X_icrf], df[i, :V_icrf])
    sv_ecef = sv_eci_to_ecef(sv_icrf, GCRF(), ITRF(), df[i, :t_utc_jd], eop)
    
    # 大地坐标
    lat, lon, alt = ecef_to_geodetic(sv_ecef.r)
    df[i, :lon] = rad2deg(lon)
    df[i, :lat] = rad2deg(lat)
    df[i, :alt] = alt
    
    # 地心坐标
    Δr = sv_ecef.r - [x_obs, y_obs, z_obs]
    range_val = norm(Δr)
    
    # ENU坐标（观察者的局部坐标系east-north-top）
    e = -sin(lon_obs)*Δr[1] + cos(lon_obs)*Δr[2]
    n = -sin(lat_obs)*cos(lon_obs)*Δr[1] - sin(lat_obs)*sin(lon_obs)*Δr[2] + cos(lat_obs)*Δr[3]
    u = cos(lat_obs)*cos(lon_obs)*Δr[1] + cos(lat_obs)*sin(lon_obs)*Δr[2] + sin(lat_obs)*Δr[3]
    
    az = atan(e, n) % (2π)
    el = atan(u, sqrt(e^2 + n^2))
    
    df[i, :az] = rad2deg(az)
    df[i, :el] = rad2deg(el)
    df[i, :range] = range_val
end

航天器的地球轨迹图

获得的大地坐标可以简单地绘制，但为了将它们绘制为线，您首先需要将行星周围的每个圆划分为单独的段，以免在图上观察到不必要的跳跃（经度180和-

In [ ]: 
ground_track_df = select(df, [:lon, :lat, :time])

# 轨道的地面路线
ground_track_lons = df.lon
ground_track_lats = df.lat

# 我们把它分成几段
jumps = findall(i -> abs(ground_track_lons[i+1] - ground_track_lons[i]) > 180, 1:length(ground_track_lons)-1)
segments = []
start_idx = 1

for jump_idx in jumps
    push!(segments, (ground_track_lons[start_idx:jump_idx], ground_track_lats[start_idx:jump_idx]))
    start_idx = jump_idx + 1
end

# 添加最后一个段
push!(segments, (ground_track_lons[start_idx:end], ground_track_lats[start_idx:end]))

# 我们分别绘制每个段
plot(title="轨道的地面路线")
for (seg_lons, seg_lats) in segments
    xy_transf = trans.(seg_lons, seg_lats)
    xx,yy = getindex.(xy_transf, 1), getindex.(xy_transf, 2)
    plot!(xx, yy, legend=false, linecolor=:blue)
end
plot!()

for shape in table
    x = [point.x for point in shape.geometry.points]
    y = [point.y for point in shape.geometry.points]
    xy_transf = trans.(x, y)
    xx,yy = getindex.(xy_transf, 1), getindex.(xy_transf, 2)
    plot!(xx, yy, c=1, lw=0.5, fill=:lightblue)
end
plot!()
Out[0]:

地球表面卫星可观测性图

为了绘制卫星从观察者所在点的可观测性，我们计算了地心方位角和赤纬。 当卫星在地平线上方时，我们将只分离轨迹的那些部分。

In [ ]: 
# 地平线以上的点
skyplot_mask = df.el .> 0

skyplot_az = df.az # df.az[skyplot_mask]#这个数据已经可以推导出来了
skyplot_el = df.el # df．el[skyplot_mask]#但是让我们把它们分成段

true_indices = findall(skyplot_mask)    # 我们把它分成连续的段
gaps = findall(diff(true_indices) .> 1) # 我们发现索引中的空白
boundaries = vcat([1], gaps .+ 1, [length(true_indices) + 1]) # 段边界

segments_az = Vector{Vector{eltype(skyplot_az)}}()
segments_el = Vector{Vector{eltype(skyplot_el)}}()
for i in 1:length(boundaries)-1
    start_idx = boundaries[i]
    end_idx = boundaries[i+1] - 1
    segment_indices = true_indices[start_idx:end_idx]
    push!(segments_az, skyplot_az[segment_indices])
    push!(segments_el, skyplot_el[segment_indices])
end

在一个单独的单元格中，我们将显示一个半球并在其上绘制所有线段，以便看到卫星从地球表面看到的轨迹。

In [ ]: 
# 创造一个半球
n = 30
u = range(0, 2π; length=n)
v = range(0, π/2; length=n)
x_sphere = cos.(u) .* sin.(v)'
y_sphere = sin.(u) .* sin.(v)'
z_sphere = ones(n) .* cos.(v)'

# 我们画一个半球
surface(x_sphere, y_sphere, z_sphere, alpha=0.4, legend=false)
xlabel!("东"); ylabel!("北区"); zlabel!("顶部")

# 绘制轨迹段
for (azimuths,elevations) in zip(segments_az, segments_el)
    # 转换为笛卡尔坐标
    az_rad = deg2rad.(azimuths)
    el_rad = deg2rad.(elevations)
    x = cos.(el_rad) .* cos.(az_rad)
    y = cos.(el_rad) .* sin.(az_rad)
    z = sin.(el_rad)

    plot!(x, y, z, lw=4)
end
plot!(legend=false, title="卫星能见度图(莫斯科)")
Out[0]:

我们可以看到卫星何时以及如何从莫斯科的一个点观察到。 由于卫星的状态矢量还包括速度，我们可以计算多普勒频移并准确指定其信号的声音，这使我们取得了新的成就。

结论

我们在60,000秒的观测中获得了PS-1卫星的坐标（真正的卫星在轨道上花费了大约800万秒），并构建了几个有用的图表来解释这个航天器的轨迹。