ylniu/dynamics
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Compare commits

base: ylniu/dynamics:974216332ceb214cda98cbe6882d7a132f271804
Branches Tags
ylniu/dynamics:main
...
compare: ylniu/dynamics:8183cbc22e42668baf85c2bd0ba975607d991ca9
Branches Tags
ylniu/dynamics:main

10 Commits
974216332c ... 8183cbc22e

Author SHA1 Message Date
admin 8183cbc22e chore: .gitignore 去掉全局 *.dll 限制，release DLL 可自由提交 2026-06-21 05:21:22 +08:00
admin 41c454eeab refactor: 引擎只走 DLL 路径，release 目录清理 exe 
- dynamics.py: 去掉 exe 子进程回退，只走 DLL 路径
- .gitignore: 白名单 engines/release/dynamics_*.dll
- engines/release/: 移除 dynamics_*.exe，新增 3 个引擎 DLL
 2026-06-21 05:19:33 +08:00
admin 6c06c00f4b fix: Fortran engine use_marker/alpha header passthrough + case04 Moon orbit fix 
- engines/fortran/main.f90:
  - 新增 use_marker 和 alpha_str 的读取、传递、写入 display.txt header
  - 新增 json_get_alpha 函数兼容 JSON 单值或 6 元数组
- engines/release/: 引擎校准缓存 & fortran.json 参数
- examples/case04/: 修正月球位置 (11.5→10.5) 和速度 (526→647),
  地月距缩至希尔半径以内 (0.5 < 1.0), 新增希尔半径验证说明
- examples/Readme.html: case04/case09 描述更新
 2026-06-19 08:05:05 +08:00
admin c8d12e7f45 docs: case04 Readme.html 新增月球绕地轨道交互图 
- 月球轨道图使用绿色配色，与地球轨道图（蓝色）区分
- 默认偏心率 e=0.0549（月球真实值）
- 标注半长轴 a / 半短轴 b / 偏心距 c
- 标注近地点/远地点距离
- 地球在左焦点，月球在椭圆轨道运行
- 独立滑块控制月球偏心率
 2026-06-18 23:06:55 +08:00
admin 7078cbc744 fix: 轨道示意图太阳移至左焦点，默认 e=0.017（地球真实 0.0167） 
- 椭圆中心改为 cx+f，使左焦点=cx 为太阳
- 偏心距标注从椭圆中心指向左焦点
- a/b 标注、近日点/远日点位置同步修正
- 滑块映射改为 e=value/1000（精度 0.001）
- 滑块默认值 17 → e=0.017 ≈ 地球真实 0.0167
 2026-06-18 23:05:23 +08:00
admin d4ce23a07e feat: case04 以月球质量=1 为基准更新质量比和质量缩放 
- coord.txt: 太阳 27000 / 地球 81 / 月球 1（月球=1 质量单位）
- 地球绕日轨道速度 v=520（v = sqrt(g_strength * M_sun / r)）
- 月球绕地速度 v=526（地球速度 + 6 绕地附加速度）
- box_a: 20 → 30（容纳更大轨道半径）
- Readme.html 同步更新质量表和初始构型表
- 已知局限说明太阳质量缩至 1/1000 的原因
 2026-06-18 23:03:36 +08:00
admin 84aed813df docs: case04 Readme.html 新增交互式轨道示意图 
使用 HTML5 Canvas + JavaScript 绘制椭圆轨道图，
支持滑块调整偏心率 e（0~0.85），实时更新：
- 椭圆形状变化
- a/b/c 参数标注
- 近日点/远日点数值
- 速度矢量方向
- 底部信息面板
 2026-06-18 23:01:04 +08:00
admin f6c468702f docs: case04 新增 Readme.html 含 MathJax 公式 
涵盖天体质量/轨道参数/偏心率/开普勒定律/近日点速度
公式等全部天体力学内容，使用 MathJax 3 (SVG) 渲染。
 2026-06-18 22:59:07 +08:00
admin 8b82d3c4e4 fix: leapfrog 模式缺少边界反弹导致粒子穿模 2026-06-18 14:28:56 +08:00
admin 80b7467421 fix: leapfrog 模式缺少边界反弹，粒子穿透地板 
问题：leapfrog 模式的主循环直接调用 leapfrog_staggered_step，
未经过 apply_motion_update，因此 Limit_in_box（边界反弹）
从未被执行。粒子超出 z<-20 后被 wrap_position 回绕到顶端，
导致物理行为完全错误。

修复：
1. 在预热和记录两个循环中，leapfrog 路径后显式调用
   Limit_in_box（与 apply_motion_update 内的非 leapfrog 路径一致）
2. 移除 apply_motion_update 末尾重复的 Limit_in_box 调用
   （统一到主循环中执行）
 2026-06-18 09:28:38 +08:00
16 changed files with 979 additions and 44 deletions
Expand all files Collapse all files
.gitignore
-1
View File
@@ -42,7 +42,6 @@ build_*/
*.so
*.so.*
*.dylib
*.dll
# 可执行文件（保留源码，排除编译出的二进制）
# 注意：Windows 下 .exe 后缀的可执行文件
compute.py
+6 -4
View File
@@ -1798,10 +1798,6 @@ def apply_motion_update(x, y, z, vx, vy, vz, dt, m, g, b):
x, y, z, vx, vy, vz = Leapfrog_Method(x, y, z, vx, vy, vz, dt, m, g, b)
else:
raise ValueError(f"未知算法: {METHOD}")
x, vx = Limit_in_box(x, X_MIN, X_MAX, vx)
y, vy = Limit_in_box(y, Y_MIN, Y_MAX, vy)
z, vz = Limit_in_box(z, Z_MIN, Z_MAX, vz)
return x, y, z, vx, vy, vz
@@ -1870,6 +1866,9 @@ def run_simulation(save_trajectory=0):
x, y, z, vx, vy, vz, DT, ATOM_MASSES, G, B)
else:
x, y, z, vx, vy, vz = apply_motion_update(x, y, z, vx, vy, vz, DT, ATOM_MASSES, G, B)
x, vx = Limit_in_box(x, X_MIN, X_MAX, vx)
y, vy = Limit_in_box(y, Y_MIN, Y_MAX, vy)
z, vz = Limit_in_box(z, Z_MIN, Z_MAX, vz)
x, y, z = wrap_position(x, y, z)
x, y, z, vx, vy, vz = apply_fixed_constraints(x, y, z, vx, vy, vz)
print(
@@ -1924,6 +1923,9 @@ def run_simulation(save_trajectory=0):
x, y, z, vx, vy, vz, DT, ATOM_MASSES, G, B)
else:
x, y, z, vx, vy, vz = apply_motion_update(x, y, z, vx, vy, vz, DT, ATOM_MASSES, G, B)
x, vx = Limit_in_box(x, X_MIN, X_MAX, vx)
y, vy = Limit_in_box(y, Y_MIN, Y_MAX, vy)
z, vz = Limit_in_box(z, Z_MIN, Z_MAX, vz)
x, y, z = wrap_position(x, y, z)
x, y, z, vx, vy, vz = apply_fixed_constraints(x, y, z, vx, vy, vz)
dynamics.py
+5 -15
View File
@@ -248,23 +248,13 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output",
input_dir_abs = str(input_dir_path.resolve())
output_dir_abs = str(output_dir_path.resolve())
# ── 优先尝试 DLL 路径（无文件 I/O，直接输出 display.npz）──
_dll_used = False
try:
# ── DLL 路径（无文件 I/O，直接输出 display.npz）──
from engines.engine_dll import is_dll_available
if is_dll_available(engine):
if not is_dll_available(engine):
raise FileNotFoundError(
f"DLL 未找到（引擎 {engine}）。"
f"请先编译：cd engines/{engine} && make dll")
compute.run_engine_dll(engine, output_dir_abs, config)
_dll_used = True
print(f"[run] DLL 路径成功")
except Exception as _dll_err:
print(f"[run] DLL 路径不可用 ({_dll_err})，回退到子进程模式")
if not _dll_used:
# 回退：子进程模式（读写 display.txt / trajectory.txt）
traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz = compute.run_engine(
engine, input_dir_abs, output_dir_abs, config)
if int(config.get("save_trajectory", 0)):
compute.save_trajectory_txt(traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, str(runtime_base))
_elapsed = _time.time() - _t0
print(f"[run] 引擎: {engine} 计算完成: {record_steps} 步 {_elapsed:.3f} s")
engines/fortran/main.f90
+70 -7
View File
@@ -21,6 +21,8 @@ program dynamics_f90
integer :: driving_force
double precision :: gravity_strength
character(len=32) :: method
integer :: use_marker
character(len=512) :: alpha_str
double precision :: t0, t1, elapsed, tw
! 原子数据
@@ -71,7 +73,8 @@ program dynamics_f90
! 读取 param.json
call read_params(param_path, box_a, NT, DT, NSTEP, warmup_steps, method, G, B, &
gravity_field, gravity_interaction, &
elastic_force, damping_force, gravity_strength, driving_force)
elastic_force, damping_force, gravity_strength, driving_force, &
use_marker, alpha_str)
! 读取 coord.txt
call read_coord(input_dir, n_atoms, atom_ids, masses, radii, pos_0, vel_0, fixed)
@@ -201,7 +204,8 @@ program dynamics_f90
traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, &
NT, DT, NSTEP, warmup_steps, method, G, B, &
n_bonds, gravity_field, elastic_force, damping_force, &
driving_force, box_a, gravity_strength)
driving_force, box_a, gravity_strength, &
use_marker, alpha_str)
call cpu_time(t1)
elapsed = t1 - t0
@@ -228,15 +232,17 @@ contains
! ========================================================================
subroutine read_params(path, box_a, NT, DT, NSTEP, warmup_steps, method, G, B, &
gravity_field, gravity_interaction, &
elastic_force, damping_force, gravity_strength, driving_force)
elastic_force, damping_force, gravity_strength, driving_force, &
use_marker, alpha_str)
character(len=*), intent(in) :: path
double precision, intent(out) :: box_a, DT, G(3), B(3), gravity_strength
integer, intent(out) :: NT, NSTEP, warmup_steps
integer, intent(out) :: gravity_field, gravity_interaction, elastic_force, damping_force, driving_force
character(len=*), intent(out) :: method
integer, intent(out) :: use_marker
character(len=*), intent(out) :: method, alpha_str
character(len=8096) :: buf
character(len=256) :: line
integer :: u, ios
integer :: u, ios, i
box_a = 10.0d0; NT = 10000; DT = 0.001d0; NSTEP = 100; warmup_steps = 0
method = 'leapfrog'
@@ -244,6 +250,8 @@ subroutine read_params(path, box_a, NT, DT, NSTEP, warmup_steps, method, G, B, &
gravity_field = 1; gravity_interaction = 0
elastic_force = 1; damping_force = 0; gravity_strength = 1.0d0
driving_force = 0
use_marker = 0
alpha_str = '0.2'
open(newunit=u, file=trim(path), status='old', action='read', iostat=ios)
if (ios /= 0) then
@@ -274,6 +282,8 @@ subroutine read_params(path, box_a, NT, DT, NSTEP, warmup_steps, method, G, B, &
damping_force = json_get_int(buf, 'damping_force', 0)
gravity_strength = json_get_double(buf, 'gravity_strength', 1.0d0)
driving_force = json_get_int(buf, 'driving_force', 0)
use_marker = json_get_int(buf, 'use_marker', 0)
call json_get_alpha(buf, 'alpha', alpha_str)
end subroutine read_params
! ========================================================================
@@ -323,6 +333,55 @@ subroutine json_get_double3(buf, key, vals)
end do
end subroutine json_get_double3
! 读取 alpha：兼容单值 (0.2) 或 6 元数组 [0.0,0.0,0.0,0.0,0.0,0.0]
! 输出为逗号分隔字符串，与 draw.py 的解析格式一致
subroutine json_get_alpha(buf, key, val)
character(len=*), intent(in) :: buf, key
character(len=*), intent(out) :: val
double precision :: d6(6)
integer :: p, i, ios
character(len=32) :: tmp
val = '0.2'
p = index(buf, '"' // key // '"')
if (p == 0) return
p = index(buf(p:), ':') + p
! skip whitespace
do while (p <= len(buf) .and. (buf(p:p) == ':' .or. buf(p:p) == ' ' &
.or. buf(p:p) == char(9) .or. buf(p:p) == char(10)))
p = p + 1
end do
if (buf(p:p) == '[') then
! 6 元数组
d6 = 0.2d0
p = p + 1
do i = 1, 6
do while (p <= len(buf) .and. (buf(p:p) == ' ' .or. buf(p:p) == ',' &
.or. buf(p:p) == char(9) .or. buf(p:p) == char(10)))
p = p + 1
end do
if (p > len(buf) .or. buf(p:p) == ']') exit
read(buf(p:), *, iostat=ios) d6(i)
if (ios /= 0) exit
do while (p <= len(buf) .and. buf(p:p) /= ',' .and. buf(p:p) /= ']')
p = p + 1
end do
end do
val = ''
do i = 1, 6
write(tmp, '(f0.6)') d6(i)
if (i > 1) val = trim(val) // ','
val = trim(val) // trim(tmp)
end do
else
! 单值标量
read(buf(p:), *, iostat=ios) d6(1)
if (ios == 0) then
write(tmp, '(f0.6)') d6(1)
val = trim(tmp)
end if
end if
end subroutine json_get_alpha
! 从 JSON 中读取字符串值
subroutine json_get_string(buf, key, val)
character(len=*), intent(in) :: buf, key
@@ -999,10 +1058,12 @@ subroutine write_display_txt(outdir, n_frames, nat, aid, &
tx, ty, tz, tvx, tvy, tvz, &
NT, DT, NSTEP, warmup, method, G, B, &
nb, gravity_field, elastic_force, damping_force, &
driving_force, box_a, gravity_strength)
character(len=*), intent(in) :: outdir, method
driving_force, box_a, gravity_strength, &
use_marker, alpha_str)
character(len=*), intent(in) :: outdir, method, alpha_str
integer, intent(in) :: n_frames, nat, NT, NSTEP, warmup, nb
integer, intent(in) :: gravity_field, elastic_force, damping_force, driving_force
integer, intent(in) :: use_marker
integer, intent(in) :: aid(nat)
double precision, intent(in) :: tx(n_frames, nat), ty(n_frames, nat), tz(n_frames, nat)
double precision, intent(in) :: tvx(n_frames, nat), tvy(n_frames, nat), tvz(n_frames, nat)
@@ -1051,6 +1112,8 @@ subroutine write_display_txt(outdir, n_frames, nat, aid, &
write(buf, '("Y_MAX: ", g0)') box_a; write(u, '(a)') trim(buf)
write(buf, '("Z_MIN: ", g0)') -box_a; write(u, '(a)') trim(buf)
write(buf, '("Z_MAX: ", g0)') box_a; write(u, '(a)') trim(buf)
write(u, '("use_marker: ", i0)') use_marker
write(u, '("alpha: ", a)') trim(alpha_str)
! ── frame data ────────────────────────────────────────────────────────
do f = 1, n_frames
engines/release/_calib_c.json
+1 -1
View File
@@ -1 +1 @@
{"n_atoms": 40, "nt": 10000, "step_time": 0.00015027170181274412}
{"n_atoms": 40, "nt": 10000, "step_time": 3.113259077072144e-05}
engines/release/_calib_fortran.json
+1
View File
@@ -0,0 +1 @@
{"n_atoms": 40, "nt": 200000, "step_time": 0.0002908185601234436}
engines/release/dynamics_c.dll
BIN
View File
Binary file not shown.
engines/release/dynamics_cpp.dll
BIN
View File
Binary file not shown.
engines/release/dynamics_f90.dll
BIN
View File
Binary file not shown.
engines/release/fortran.json
+51
View File
@@ -0,0 +1,51 @@
{
"box_a": 300.0,
"NT": 200000,
"DT": 0.001,
"NSTEP": 100,
"warmup_steps": 0,
"method": "leapfrog",
"G": [
0.0,
0.0,
0.0
],
"B": [
0.005,
0.0,
0.005
],
"gravity_field": 0,
"gravity_interaction": 0,
"elastic_force": 1,
"damping_force": 0,
"gravity_strength": 1.0,
"driving_force": 1,
"save_trajectory": 0,
"alpha": [
0.0,
0.0,
0.0,
0.0,
0.0,
0.0
],
"ball_radius": 0.5,
"ball_color": [
0.2,
0.6,
0.9
],
"box_color": [
0.8,
0.8,
0.85
],
"use_marker": 1,
"camera_distance": 60.0,
"camera_elevation": 0.0,
"camera_azimuth": 0.0,
"camera_center_x": 30.0,
"camera_center_y": 0.0,
"camera_center_z": 0.0
}
examples/Readme.html
+6 -6
View File
@@ -210,7 +210,7 @@
<div class="case-card">
<span class="num">04</span>
<h3>日地月系统（稳定）</h3>
<p>三体系统稳定椭圆轨道</p>
<p>三体系统稳定轨道，经希尔半径检验</p>
<div class="meta">
<span class="tag tag-python">Python</span>
<span class="tag tag-gravity">万有引力</span>
@@ -263,7 +263,7 @@
<div class="case-card">
<span class="num">09</span>
<h3>一维链纵波·Fortran 引擎</h3>
<p>Fortran 引擎驱动的纵波传播测试</p>
<p>Fortran 引擎驱动的纵波，GPU 实例化渲染</p>
<div class="meta">
<span class="tag tag-fortran">Fortran</span>
<span class="tag tag-spring">弹簧</span>
@@ -323,8 +323,8 @@
<li class="tag tag-python">Python 引擎</li>
<li class="tag tag-gravity">万有引力</li>
</ul>
<p>与 case03 相同的三体系统，但采用恰当的初始条件，地球和月球维持稳定的椭圆轨道运动。</p>
<div class="highlight">✅ 与 case03 对比学习：初始条件对数值稳定性的影响</div>
<p>与 case03 相同的三体系统，但采用恰当的初始条件：地球置于近日点（$r=10$），$v_z=520$ 接近圆轨道；月球缩至希尔半径（$r_H\approx1.0$）以内的 $r=0.5$，相对速度 $v_{\text{rel}}\approx127$，确保月球被地球稳定束缚。</p>
<div class="highlight">✅ 与 case03 对比学习：初始条件对数值稳定性的影响。地月距 0.5 在地球希尔半径以内，满足稳定性条件</div>
</div>
<div class="detail-card">
@@ -377,8 +377,8 @@
<li class="tag tag-spring">弹簧键力</li>
<li class="tag tag-drive">驱动力</li>
</ul>
<p>40 原子沿 x 轴排列，Fortran 引擎驱动的纵波传播测试。验证 Fortran 引擎的输出兼容性。</p>
<div class="highlight">🔧 Fortran 引擎兼容性验证，T_total=200, NSTEP=100</div>
<p>40 原子沿 x 轴排列，Fortran 引擎驱动的纵波传播测试。使用 <code>use_marker: 1</code>（GPU 实例化 Marker 模式）加速渲染，验证 Fortran 引擎与其他引擎的输出兼容性。</p>
<div class="highlight">🔧 Fortran 引擎兼容性验证，T_total=200, NSTEP=100。Marker 模式下 40 原子以 GPU 点精灵渲染，帧率高</div>
</div>
<div class="detail-card">
examples/case01/input/input.txt
+1 -1
View File
@@ -6,7 +6,7 @@
# 每步用 0/1 单独开关，1=执行，0=跳过
# 依赖关系：抽帧依赖模拟结果，绘图依赖模拟+抽帧
step_simulate: 1 # 运行物理模拟 → output/trajectory.txt
step_sample: 1 # 抽帧 → output/display.txt
step_sample: 0 # 抽帧 → output/display.txt
step_plot: 1 # 绘制轨迹/能量图 → output/trajectory_plots.png
step_plot_wave: 0 # 绘制波形能量动画
plot_wave_save_gif: 0 # 输出波形 GIF（需 step_plot_wave=1）
examples/case04/Readme.html
+829
View File
@@ -0,0 +1,829 @@
<!DOCTYPE html>
<html lang="zh-CN">
<head>
<meta charset="UTF-8">
<meta name="viewport" content="width=device-width, initial-scale=1.0">
<title>case04 — 日地月三体系统</title>
<script>
MathJax = {
tex: { inlineMath: [['$','$'], ['\\(','\\)']], displayMath: [['$$','$$'], ['\\[','\\]']] },
svg: { fontCache: 'global' }
};
</script>
<script src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-svg.js" async></script>
<style>
:root {
--bg: #0d1117;
--surface: #161b22;
--border: #30363d;
--text: #c9d1d9;
--text-dim: #8b949e;
--accent: #58a6ff;
--green: #3fb950;
--orange: #d29922;
--red: #f85149;
--teal: #56d4dd;
}
* { margin: 0; padding: 0; box-sizing: border-box; }
body {
font-family: -apple-system, BlinkMacSystemFont, "Segoe UI", "Noto Sans SC", Helvetica, Arial, sans-serif;
background: var(--bg);
color: var(--text);
line-height: 1.8;
padding: 40px 24px;
}
.container { max-width: 860px; margin: 0 auto; }
h1 { font-size: 1.8rem; margin-bottom: 6px; color: #f0f6fc; }
h1 small { font-size: 1rem; color: var(--text-dim); font-weight: 400; }
.subtitle { color: var(--text-dim); margin-bottom: 32px; }
h2 {
font-size: 1.35rem;
margin-top: 36px;
margin-bottom: 14px;
padding-bottom: 6px;
border-bottom: 1px solid var(--border);
color: #f0f6fc;
}
h3 { font-size: 1.1rem; margin-top: 24px; margin-bottom: 10px; color: #f0f6fc; }
p { margin-bottom: 12px; color: var(--text); }
.card {
background: var(--surface);
border: 1px solid var(--border);
border-radius: 8px;
padding: 20px 24px;
margin-bottom: 20px;
}
table {
width: 100%;
border-collapse: collapse;
margin: 12px 0;
font-size: 0.9rem;
}
th, td {
text-align: left;
padding: 8px 14px;
border-bottom: 1px solid var(--border);
}
th { color: var(--text-dim); font-weight: 600; background: #1c2128; }
.highlight {
background: #1f242e;
border-left: 3px solid var(--accent);
padding: 10px 16px;
margin: 12px 0;
border-radius: 0 6px 6px 0;
font-size: 0.9rem;
}
.highlight-eq {
background: #1f242e;
border-left: 3px solid var(--teal);
padding: 14px 18px;
margin: 14px 0;
border-radius: 0 6px 6px 0;
}
.highlight-eq p { margin: 6px 0; }
code {
background: #1c2128;
padding: 2px 6px;
border-radius: 4px;
font-family: "SF Mono", "Fira Code", monospace;
font-size: 0.85em;
color: var(--green);
}
.tag {
display: inline-block;
padding: 2px 10px;
border-radius: 4px;
font-size: 0.78rem;
font-weight: 500;
margin: 2px;
}
.tag-python { background: #3572A533; color: #3572A5; }
.tag-gravity { background: #d2992233; color: var(--orange); }
ul, ol { margin: 8px 0 12px 24px; }
li { margin-bottom: 4px; }
.formula-block {
overflow-x: auto;
padding: 8px 0;
}
.diagram-wrap {
background: #0d1117;
border: 1px solid var(--border);
border-radius: 8px;
padding: 16px;
margin: 16px 0;
text-align: center;
}
.diagram-wrap canvas {
display: block;
margin: 0 auto;
max-width: 100%;
height: auto;
}
.diagram-controls {
display: flex;
flex-wrap: wrap;
align-items: center;
justify-content: center;
gap: 16px;
margin-top: 12px;
}
.diagram-controls label {
font-size: 0.85rem;
color: var(--text-dim);
}
.diagram-controls input[type="range"] {
width: 200px;
accent-color: var(--accent);
}
.diagram-controls .val {
font-size: 0.9rem;
font-weight: 500;
color: var(--teal);
min-width: 48px;
display: inline-block;
text-align: center;
}
.diagram-info {
display: flex;
flex-wrap: wrap;
justify-content: center;
gap: 20px;
margin-top: 10px;
font-size: 0.82rem;
color: var(--text-dim);
}
.diagram-info span {
display: inline-flex;
align-items: center;
gap: 4px;
}
.diagram-info .dot {
display: inline-block;
width: 10px;
height: 10px;
border-radius: 50%;
margin-right: 2px;
}
@media (max-width: 640px) {
body { padding: 16px; }
.diagram-controls { flex-direction: column; gap: 8px; }
}
</style>
</head>
<body>
<div class="container">
<h1>case04 — 日地月三体系统 <small>稳定轨道版本</small></h1>
<p class="subtitle">太阳、地球、月球三体系统，采用真实比例的质量和万有引力模拟，地球和月球维持稳定椭圆轨道。</p>
<div class="card" style="display:flex; flex-wrap:wrap; gap:8px; align-items:center; margin-bottom:20px;">
<span class="tag tag-python">Python 引擎</span>
<span class="tag tag-gravity">万有引力</span>
<span style="color:var(--text-dim); font-size:0.85rem; margin-left:8px;">3 粒子 | leapfrog 算法 | 稳定轨道</span>
</div>
<h2>一、物理系统概览</h2>
<div class="card">
<p>本案例模拟真实的日-地-月三体系统：</p>
<ul>
<li><strong>太阳</strong>：位于原点，质量为 $M_\odot$，固定不动（fix_x=fix_y=fix_z=1）</li>
<li><strong>地球</strong>：绕太阳公转，质量为 $M_\oplus$，近日点出发</li>
<li><strong>月球</strong>：绕地球公转（同时跟随地球绕太阳），质量为 $M_\text{moon}$</li>
</ul>
<p>原子间万有引力由 <code>gravity_strength</code> 缩放控制，$F = G \dfrac{m_1 m_2}{r^2}$。</p>
</div>
<h2>二、天体真实参数</h2>
<h3>2.1 质量</h3>
<table>
<tr><th>天体</th><th>质量 / kg</th><th>模拟质量（月球=1）</th></tr>
<tr><td>太阳</td><td>$1.989 \times 10^{30}$</td><td>$27\,000$</td></tr>
<tr><td>地球</td><td>$5.972 \times 10^{24}$</td><td>$81$</td></tr>
<tr><td>月球</td><td>$7.35 \times 10^{22}$</td><td>$1$</td></tr>
</table>
<div class="highlight">
<strong>质量比</strong>：$M_\odot : M_\oplus : M_\text{moon} \approx 27\,100\,000 : 81.3 : 1$<br>
模拟中采用缩放比例 $27\,000 : 81 : 1$（太阳质量缩至 $1/1000$ 以保持数值稳定）。
</div>
<h3>2.2 轨道参数</h3>
<table>
<tr><th>轨道</th><th>半长轴 $a$</th><th>偏心率 $e$</th><th>半短轴 $b$</th><th>周期</th></tr>
<tr><td>地球绕太阳</td><td>$1.4960 \times 10^8$ km (1 AU)</td><td>$0.0167$</td><td>$1.4958 \times 10^8$ km</td><td>365.25 天</td></tr>
<tr><td>月球绕地球</td><td>$3.844 \times 10^5$ km</td><td>$0.0549$</td><td>$3.838 \times 10^5$ km</td><td>27.32 天</td></tr>
</table>
<div class="highlight">
<strong>比例</strong>：$R_{\text{日地}} : R_{\text{地月}} \approx 389 : 1$<br>
地月距离约为日地距离的 $1/389$。
</div>
<h2>三、轨道力学</h2>
<h3>3.1 交互式轨道示意图</h3>
<p>拖动下方滑块改变偏心率 $e$，观察轨道形状的变化和参数标注：</p>
<div class="diagram-wrap">
<canvas id="orbitCanvas" width="760" height="400"></canvas>
<div class="diagram-controls">
<label>偏心率 e = <span class="val" id="eDisplay">0.0170</span>（地球真实值 0.0167）</label>
<input type="range" id="eSlider" min="0" max="850" value="17">
</div>
<div class="diagram-info" id="diagramInfo"></div>
</div>
<h3>3.2 月球绕地球轨道</h3>
<p>拖动下方滑块改变月球轨道偏心率 $e_\text{moon}$（真实值 0.0549）：</p>
<div class="diagram-wrap">
<canvas id="moonCanvas" width="760" height="400"></canvas>
<div class="diagram-controls">
<label>偏心率 e = <span class="val" id="moonEDisplay">0.0549</span>（月球真实值 0.0549）</label>
<input type="range" id="moonESlider" min="0" max="850" value="55">
</div>
<div class="diagram-info" id="moonInfo"></div>
</div>
<h3>3.3 偏心率定义</h3>
<div class="highlight-eq">
<p>偏心率 $e$ 描述椭圆轨道偏离正圆的程度：</p>
<div class="formula-block">
$$e = \frac{c}{a}$$
</div>
<p>其中 $c = ea$ 为偏心距（焦点到椭圆中心的距离），$a$ 为半长轴。</p>
</div>
<table>
<tr><th>$e$ 值</th><th>轨道形状</th><th>说明</th></tr>
<tr><td>$e = 0$</td><td>正圆形</td><td>速度恒定，距离恒定</td></tr>
<tr><td>$0 &lt; e &lt; 1$</td><td>椭圆</td><td>近日点/近地点速度最大，远日点/远地点速度最小</td></tr>
<tr><td>$e = 1$</td><td>抛物线</td><td>逃逸轨道，速度恰好达到逃逸速度</td></tr>
<tr><td>$e &gt; 1$</td><td>双曲线</td><td>飞越轨道，速度超过逃逸速度</td></tr>
</table>
<h3>3.4 半长轴与半短轴的关系</h3>
<div class="highlight-eq">
<div class="formula-block">
$$b = a\sqrt{1 - e^2}$$
</div>
<p>当 $e \ll 1$ 时，$b \approx a\left(1 - \dfrac{e^2}{2}\right)$，椭圆度非常微小。</p>
</div>
<h3>3.5 近日点与远日点</h3>
<div class="highlight-eq">
<p>近日点（距太阳最近）和远日点（距太阳最远）的距离分别为：</p>
<div class="formula-block">
$$r_{\text{peri}} = a(1 - e), \qquad r_{\text{ap}} = a(1 + e)$$
</div>
<p>在近日点轨道速度最大，在远日点轨道速度最小：</p>
<div class="formula-block">
$$v_{\text{peri}} = \sqrt{\frac{GM(1+e)}{a(1-e)}}, \qquad
v_{\text{ap}} = \sqrt{\frac{GM(1-e)}{a(1+e)}}$$
</div>
</div>
<h3>3.6 开普勒三大定律</h3>
<ol>
<li><strong>椭圆定律</strong>：行星轨道是椭圆，太阳位于椭圆的一个焦点上。</li>
<li><strong>面积定律</strong>：行星与太阳的连线在相等时间内扫过相等的面积。</li>
<li><strong>周期定律</strong>：公转周期的平方与半长轴的立方成正比：$T^2 \propto a^3$。</li>
</ol>
<div class="highlight">
<strong>验证</strong>：在本模拟中，你可以通过轨迹图观察面积定律是否成立——地球在近日点移动更快，远日点移动更慢。
</div>
<h2>四、模拟参数</h2>
<h3>4.1 缩放说明</h3>
<p>真实尺度无法直接用于模拟（日地距离 1.5 亿 km），因此采用缩放参数。当前 case04 的配置为教学演示而简化：</p>
<table>
<tr><th>参数</th><th>值</th><th>说明</th></tr>
<tr><td><code>gravity_strength</code></td><td>100.0</td><td>万有引力强度</td></tr>
<tr><td><code>box_a</code></td><td>30.0</td><td>盒子半边长</td></tr>
<tr><td><code>method</code></td><td>leapfrog</td><td>蛙跳法（辛积分器，能量守恒）</td></tr>
<tr><td><code>T_total</code></td><td>10.0 s</td><td>总模拟时间</td></tr>
<tr><td><code>NSTEP</code></td><td>2</td><td>抽帧间隔（密采帧）</td></tr>
</table>
<h3>4.2 初始构型</h3>
<table>
<tr><th>天体</th><th>质量</th><th>位置 $(x,y,z)$</th><th>速度 $(v_x,v_y,v_z)$</th><th>固定约束</th></tr>
<tr><td>太阳</td><td>$27\,000$</td><td>$(0,0,0)$</td><td>$(0,0,0)$</td><td>全部固定</td></tr>
<tr><td>地球</td><td>$81$</td><td>$(10,0,0)$</td><td>$(0,0,520)$</td><td>无</td></tr>
<tr><td>月球</td><td>$1$</td><td>$(10.5,0,0)$</td><td>$(0,0,647)$</td><td>无</td></tr>
</table>
<p>地球在 $z$ 方向获得初速度 $v=520$，产生绕太阳的轨道运动（接近圆轨道）；月球在地球基础上附加 $v\approx127$ 的绕地速度，形成绕地球的轨道。速度由圆形轨道公式 $v = \sqrt{G_{\text{eff}} M / r}$ 计算，其中 $G_{\text{eff}} = \text{gravity\_strength} = 100$。</p>
<div class="highlight">
<strong>验证</strong>：地球速度 $v_\oplus = \sqrt{100 \times 27\,000 / 10} \approx 519.6$，设 520 正确。月球相对速度 $v_{\text{rel}} = \sqrt{100 \times 81 / 0.5} \approx 127.3$，设 127 正确。地月距离 0.5 在地球希尔半径 $r_H \approx 10 \times (81 / 81\,000)^{1/3} \approx 1.0$ 之内，可确保轨道稳定。
</div>
<h2>五、使用方法</h2>
<div class="card">
<code style="display:block; padding:14px 18px; margin-bottom:10px;">
# 进入 case04 目录并运行<br>
cd examples/case04<br>
python run_dynamics.py<br><br>
# 仅输出轨迹图，跳过动画<br>
python run_dynamics.py --no-plot
</code>
<p>配置文件：<code>input/input.txt</code>（物理参数）、<code>input/coord.txt</code>（初始位置/速度）。</p>
</div>
<h2>六、与 case03 的对比</h2>
<table>
<tr><th></th><th>case03（失败案例）</th><th>case04（成功案例）</th></tr>
<tr><td>轨道状态</td><td>轨道发散或碰撞</td><td>稳定椭圆轨道</td></tr>
<tr><td>关键差异</td><td>初始速度或质量比不恰当</td><td>合理的初值和参数</td></tr>
<tr><td>教学意义</td><td>展示参数选择的重要性</td><td>展示正确的三体运动</td></tr>
</table>
<div class="highlight">
<strong>教学建议</strong>：先运行 case03 观察失稳，再运行 case04 对比稳定轨道，理解初始条件对数值模拟的关键影响。
</div>
<h2>七、已知局限</h2>
<ul>
<li>太阳真实质量为月球 $27\,100\,000$ 倍，模拟中缩至 $27\,000$ 倍以保持数值稳定（$1/1000$）</li>
<li>轨道半径未按真实比例缩放（真实日地距是地月距的 389 倍，模拟中约为 20 倍，受希尔半径约束，月球已尽可能放远）</li>
<li>未考虑月球轨道倾角（真实地月轨道有约 5° 的倾角）</li>
<li>leapfrog 算法为能量守恒的辛积分器，长期稳定，但步长过大时仍可能偏离真实轨道</li>
<li>太阳固定不动可作为教学简化，但严格三体模拟中应让所有天体在质心系中自由运动</li>
</ul>
</div>
<script>
(function() {
var canvas = document.getElementById('orbitCanvas');
var ctx = canvas.getContext('2d');
var slider = document.getElementById('eSlider');
var eDisplay = document.getElementById('eDisplay');
var infoDiv = document.getElementById('diagramInfo');
var W = 760, H = 400;
var cx = 380, cy = 200;
// 实际绘图区域半径
var A = 140;
function draw(e) {
ctx.clearRect(0, 0, W, H);
var b = A * Math.sqrt(1 - e * e);
var f = e * A;
// 长轴辅助线
ctx.beginPath();
ctx.moveTo(cx - A - 20, cy);
ctx.lineTo(cx + A + 20, cy);
ctx.strokeStyle = '#30363d';
ctx.lineWidth = 0.5;
ctx.setLineDash([4, 3]);
ctx.stroke();
ctx.setLineDash([]);
// 椭圆轨道（中心在 cx+f，左焦点=cx 为太阳）
ctx.beginPath();
ctx.ellipse(cx + f, cy, A, b, 0, 0, Math.PI * 2);
ctx.strokeStyle = '#58a6ff';
ctx.lineWidth = 1.5;
ctx.stroke();
// 半长轴 a（标注线，从椭圆中心到右顶点）
var aY = cy + 28;
ctx.beginPath();
ctx.moveTo(cx + f, aY);
ctx.lineTo(cx + f + A, aY);
ctx.strokeStyle = '#c9d1d9';
ctx.lineWidth = 1.5;
ctx.stroke();
// 箭头（左端：椭圆中心 cx+f）
ctx.beginPath();
ctx.moveTo(cx + f + 6, aY - 4);
ctx.lineTo(cx + f, aY);
ctx.lineTo(cx + f + 6, aY + 4);
ctx.strokeStyle = '#c9d1d9';
ctx.lineWidth = 1.5;
ctx.stroke();
// 箭头（右端：右顶点 cx+f+A）
ctx.beginPath();
ctx.moveTo(cx + f + A - 6, aY - 4);
ctx.lineTo(cx + f + A, aY);
ctx.lineTo(cx + f + A - 6, aY + 4);
ctx.strokeStyle = '#c9d1d9';
ctx.lineWidth = 1.5;
ctx.stroke();
ctx.fillStyle = '#c9d1d9';
ctx.font = '13px "Segoe UI", Arial, sans-serif';
ctx.textAlign = 'center';
ctx.textBaseline = 'top';
ctx.fillText('a = ' + A.toFixed(0), cx + f + A/2, aY + 6);
// 偏心距 c = ea（从左焦点到椭圆中心）
var cY = cy - b - 20;
ctx.beginPath();
ctx.moveTo(cx, cY);
ctx.lineTo(cx + f, cY);
ctx.strokeStyle = '#d29922';
ctx.lineWidth = 1.5;
ctx.stroke();
// 箭头（左端：左焦点=cx 太阳）
ctx.beginPath();
ctx.moveTo(cx + 5, cY - 4);
ctx.lineTo(cx, cY);
ctx.lineTo(cx + 5, cY + 4);
ctx.strokeStyle = '#d29922';
ctx.lineWidth = 1.5;
ctx.stroke();
// 箭头（右端：椭圆中心 cx+f）
ctx.beginPath();
ctx.moveTo(cx + f - 5, cY - 4);
ctx.lineTo(cx + f, cY);
ctx.lineTo(cx + f - 5, cY + 4);
ctx.strokeStyle = '#d29922';
ctx.lineWidth = 1.5;
ctx.stroke();
ctx.fillStyle = '#d29922';
ctx.font = '13px "Segoe UI", Arial, sans-serif';
ctx.textAlign = 'center';
ctx.textBaseline = 'bottom';
ctx.fillText('c = ea = ' + f.toFixed(1), cx + f/2, cY - 4);
// 半短轴 b 标注（从椭圆中心到上顶点）
var bX = cx + f + A + 18;
ctx.beginPath();
ctx.moveTo(bX, cy);
ctx.lineTo(bX, cy - b);
ctx.strokeStyle = '#3fb950';
ctx.lineWidth = 1.5;
ctx.stroke();
ctx.beginPath();
ctx.moveTo(bX - 4, cy - 6);
ctx.lineTo(bX, cy);
ctx.lineTo(bX + 4, cy - 6);
ctx.strokeStyle = '#3fb950';
ctx.lineWidth = 1.5;
ctx.stroke();
ctx.beginPath();
ctx.moveTo(bX - 4, cy - b + 6);
ctx.lineTo(bX, cy - b);
ctx.lineTo(bX + 4, cy - b + 6);
ctx.strokeStyle = '#3fb950';
ctx.lineWidth = 1.5;
ctx.stroke();
ctx.fillStyle = '#3fb950';
ctx.font = '13px "Segoe UI", Arial, sans-serif';
ctx.textAlign = 'left';
ctx.textBaseline = 'middle';
ctx.fillText('b = ' + b.toFixed(1), bX + 8, cy - b/2);
// 太阳（焦点）
ctx.beginPath();
ctx.arc(cx, cy, 7, 0, Math.PI * 2);
ctx.fillStyle = '#d29922';
ctx.fill();
ctx.strokeStyle = '#f0883e';
ctx.lineWidth = 1;
ctx.stroke();
ctx.fillStyle = '#f0f6fc';
ctx.font = '13px "Segoe UI", Arial, sans-serif';
ctx.textAlign = 'center';
ctx.textBaseline = 'bottom';
ctx.fillText('太阳（左焦点）', cx, cy - 12);
// 近日点（左顶点：距左焦点最近）
var periX = cx + f - A;
ctx.fillStyle = '#58a6ff';
ctx.beginPath();
ctx.arc(periX, cy, 4, 0, Math.PI * 2);
ctx.fill();
ctx.fillStyle = '#8b949e';
ctx.font = '12px "Segoe UI", Arial, sans-serif';
ctx.textAlign = 'center';
ctx.textBaseline = 'top';
ctx.fillText('近日点 r = a(1-e) = ' + (A * (1 - e)).toFixed(1), periX, cy + 12);
// 远日点标注（右顶点：距左焦点最远）
var apX = cx + f + A;
ctx.fillStyle = '#58a6ff';
ctx.beginPath();
ctx.arc(apX, cy, 4, 0, Math.PI * 2);
ctx.fill();
ctx.fillStyle = '#8b949e';
ctx.font = '12px "Segoe UI", Arial, sans-serif';
ctx.textAlign = 'center';
ctx.textBaseline = 'top';
ctx.fillText('远日点 r = a(1+e) = ' + (A * (1 + e)).toFixed(1), apX, cy + 12);
// 地球（轨道上的点）
var angle = 0.4;
var ex = cx + f + A * Math.cos(angle);
var ey = cy - b * Math.sin(angle);
ctx.beginPath();
ctx.arc(ex, ey, 5, 0, Math.PI * 2);
ctx.fillStyle = '#58a6ff';
ctx.fill();
ctx.strokeStyle = '#1f6feb';
ctx.lineWidth = 1;
ctx.stroke();
ctx.fillStyle = '#58a6ff';
ctx.font = '12px "Segoe UI", Arial, sans-serif';
ctx.textAlign = 'left';
ctx.textBaseline = 'bottom';
ctx.fillText('地球', ex + 8, ey);
// 速度矢量
var vlen = 32;
var vx = -A * Math.sin(angle) * 0.7;
var vy = b * Math.cos(angle) * 0.7;
var vl = Math.sqrt(vx*vx + vy*vy);
if (vl > 0) {
vx = vx / vl * vlen;
vy = vy / vl * vlen;
}
ctx.beginPath();
ctx.moveTo(ex, ey);
ctx.lineTo(ex + vx, ey + vy);
ctx.strokeStyle = '#3fb950';
ctx.lineWidth = 2;
ctx.stroke();
ctx.beginPath();
var tipX = ex + vx, tipY = ey + vy;
var a1 = 0.3;
ctx.moveTo(tipX, tipY);
ctx.lineTo(tipX - vx * a1 + vy * a1 * 0.5, tipY - vy * a1 - vx * a1 * 0.5);
ctx.lineTo(tipX - vx * a1 - vy * a1 * 0.5, tipY - vy * a1 + vx * a1 * 0.5);
ctx.closePath();
ctx.fillStyle = '#3fb950';
ctx.fill();
// 信息面板
var roundE = Math.round(e * 10000) / 10000;
infoDiv.innerHTML =
'<span><span class="dot" style="background:#58a6ff"></span>椭圆轨道 e = ' + roundE.toFixed(4) + '</span>' +
'<span><span class="dot" style="background:#d29922"></span>偏心距 c = ' + f.toFixed(1) + '</span>' +
'<span><span class="dot" style="background:#3fb950"></span>b/a = ' + (b / A).toFixed(4) + '</span>' +
'<span>近日点 ' + (A * (1 - e)).toFixed(1) + ' / 远日点 ' + (A * (1 + e)).toFixed(1) + '</span>';
}
function update() {
var val = parseInt(slider.value);
var e = val / 1000;
eDisplay.textContent = e.toFixed(4);
draw(e);
}
slider.addEventListener('input', update);
update();
})();
// ── 月球绕地球轨道图 ──
(function() {
var canvas = document.getElementById('moonCanvas');
var ctx = canvas.getContext('2d');
var slider = document.getElementById('moonESlider');
var eDisplay = document.getElementById('moonEDisplay');
var infoDiv = document.getElementById('moonInfo');
var W = 760, H = 400;
var cx = 380, cy = 200;
var A = 140;
function draw(e) {
ctx.clearRect(0, 0, W, H);
var b = A * Math.sqrt(1 - e * e);
var f = e * A;
// 长轴辅助线
ctx.beginPath();
ctx.moveTo(cx - A - 20, cy);
ctx.lineTo(cx + A + 20, cy);
ctx.strokeStyle = '#30363d';
ctx.lineWidth = 0.5;
ctx.setLineDash([4, 3]);
ctx.stroke();
ctx.setLineDash([]);
// 椭圆轨道（中心在 cx+f，左焦点=cx 为地球）
ctx.beginPath();
ctx.ellipse(cx + f, cy, A, b, 0, 0, Math.PI * 2);
ctx.strokeStyle = '#7ee787';
ctx.lineWidth = 1.5;
ctx.stroke();
// 半长轴 a（标注线，从椭圆中心到右顶点）
var aY = cy + 28;
ctx.beginPath();
ctx.moveTo(cx + f, aY);
ctx.lineTo(cx + f + A, aY);
ctx.strokeStyle = '#c9d1d9';
ctx.lineWidth = 1.5;
ctx.stroke();
ctx.beginPath();
ctx.moveTo(cx + f + 6, aY - 4);
ctx.lineTo(cx + f, aY);
ctx.lineTo(cx + f + 6, aY + 4);
ctx.strokeStyle = '#c9d1d9';
ctx.lineWidth = 1.5;
ctx.stroke();
ctx.beginPath();
ctx.moveTo(cx + f + A - 6, aY - 4);
ctx.lineTo(cx + f + A, aY);
ctx.lineTo(cx + f + A - 6, aY + 4);
ctx.strokeStyle = '#c9d1d9';
ctx.lineWidth = 1.5;
ctx.stroke();
ctx.fillStyle = '#c9d1d9';
ctx.font = '13px "Segoe UI", Arial, sans-serif';
ctx.textAlign = 'center';
ctx.textBaseline = 'top';
ctx.fillText('a = ' + A.toFixed(0), cx + f + A/2, aY + 6);
// 偏心距 c = ea（从左焦点=地球 到椭圆中心）
var cY = cy - b - 20;
ctx.beginPath();
ctx.moveTo(cx, cY);
ctx.lineTo(cx + f, cY);
ctx.strokeStyle = '#d29922';
ctx.lineWidth = 1.5;
ctx.stroke();
ctx.beginPath();
ctx.moveTo(cx + 5, cY - 4);
ctx.lineTo(cx, cY);
ctx.lineTo(cx + 5, cY + 4);
ctx.strokeStyle = '#d29922';
ctx.lineWidth = 1.5;
ctx.stroke();
ctx.beginPath();
ctx.moveTo(cx + f - 5, cY - 4);
ctx.lineTo(cx + f, cY);
ctx.lineTo(cx + f - 5, cY + 4);
ctx.strokeStyle = '#d29922';
ctx.lineWidth = 1.5;
ctx.stroke();
ctx.fillStyle = '#d29922';
ctx.font = '13px "Segoe UI", Arial, sans-serif';
ctx.textAlign = 'center';
ctx.textBaseline = 'bottom';
ctx.fillText('c = ea = ' + f.toFixed(1), cx + f/2, cY - 4);
// 半短轴 b 标注（从椭圆中心到上顶点）
var bX = cx + f + A + 18;
ctx.beginPath();
ctx.moveTo(bX, cy);
ctx.lineTo(bX, cy - b);
ctx.strokeStyle = '#3fb950';
ctx.lineWidth = 1.5;
ctx.stroke();
ctx.beginPath();
ctx.moveTo(bX - 4, cy - 6);
ctx.lineTo(bX, cy);
ctx.lineTo(bX + 4, cy - 6);
ctx.strokeStyle = '#3fb950';
ctx.lineWidth = 1.5;
ctx.stroke();
ctx.beginPath();
ctx.moveTo(bX - 4, cy - b + 6);
ctx.lineTo(bX, cy - b);
ctx.lineTo(bX + 4, cy - b + 6);
ctx.strokeStyle = '#3fb950';
ctx.lineWidth = 1.5;
ctx.stroke();
ctx.fillStyle = '#3fb950';
ctx.font = '13px "Segoe UI", Arial, sans-serif';
ctx.textAlign = 'left';
ctx.textBaseline = 'middle';
ctx.fillText('b = ' + b.toFixed(1), bX + 8, cy - b/2);
// 地球（左焦点）
ctx.beginPath();
ctx.arc(cx, cy, 8, 0, Math.PI * 2);
ctx.fillStyle = '#58a6ff';
ctx.fill();
ctx.strokeStyle = '#1f6feb';
ctx.lineWidth = 1;
ctx.stroke();
ctx.fillStyle = '#f0f6fc';
ctx.font = '13px "Segoe UI", Arial, sans-serif';
ctx.textAlign = 'center';
ctx.textBaseline = 'bottom';
ctx.fillText('地球（左焦点）', cx, cy - 14);
// 近地点（左顶点：距地球最近）
var periX = cx + f - A;
ctx.fillStyle = '#7ee787';
ctx.beginPath();
ctx.arc(periX, cy, 4, 0, Math.PI * 2);
ctx.fill();
ctx.fillStyle = '#8b949e';
ctx.font = '12px "Segoe UI", Arial, sans-serif';
ctx.textAlign = 'center';
ctx.textBaseline = 'top';
ctx.fillText('近地点 r = a(1-e) = ' + (A * (1 - e)).toFixed(1), periX, cy + 12);
// 远地点（右顶点：距地球最远）
var apX = cx + f + A;
ctx.fillStyle = '#7ee787';
ctx.beginPath();
ctx.arc(apX, cy, 4, 0, Math.PI * 2);
ctx.fill();
ctx.fillStyle = '#8b949e';
ctx.font = '12px "Segoe UI", Arial, sans-serif';
ctx.textAlign = 'center';
ctx.textBaseline = 'top';
ctx.fillText('远地点 r = a(1+e) = ' + (A * (1 + e)).toFixed(1), apX, cy + 12);
// 月球（轨道上的点）
var angle = 0.6;
var mx = cx + f + A * Math.cos(angle);
var my = cy - b * Math.sin(angle);
ctx.beginPath();
ctx.arc(mx, my, 4, 0, Math.PI * 2);
ctx.fillStyle = '#e2e8f0';
ctx.fill();
ctx.strokeStyle = '#8b949e';
ctx.lineWidth = 0.5;
ctx.stroke();
ctx.fillStyle = '#e2e8f0';
ctx.font = '12px "Segoe UI", Arial, sans-serif';
ctx.textAlign = 'left';
ctx.textBaseline = 'bottom';
ctx.fillText('月球', mx + 8, my);
// 速度矢量
var vlen = 32;
var vx = -A * Math.sin(angle) * 0.7;
var vy = b * Math.cos(angle) * 0.7;
var vl = Math.sqrt(vx*vx + vy*vy);
if (vl > 0) { vx = vx/vl*vlen; vy = vy/vl*vlen; }
ctx.beginPath();
ctx.moveTo(mx, my);
ctx.lineTo(mx + vx, my + vy);
ctx.strokeStyle = '#3fb950';
ctx.lineWidth = 2;
ctx.stroke();
var tipX = mx + vx, tipY = my + vy;
ctx.beginPath();
ctx.moveTo(tipX, tipY);
ctx.lineTo(tipX - vx*0.3 + vy*0.15, tipY - vy*0.3 - vx*0.15);
ctx.lineTo(tipX - vx*0.3 - vy*0.15, tipY - vy*0.3 + vx*0.15);
ctx.closePath();
ctx.fillStyle = '#3fb950';
ctx.fill();
// 信息面板
var roundE = Math.round(e * 10000) / 10000;
infoDiv.innerHTML =
'<span><span class="dot" style="background:#7ee787"></span>椭圆轨道 e = ' + roundE.toFixed(4) + '</span>' +
'<span><span class="dot" style="background:#d29922"></span>偏心距 c = ' + f.toFixed(1) + '</span>' +
'<span><span class="dot" style="background:#3fb950"></span>b/a = ' + (b / A).toFixed(4) + '</span>' +
'<span>近地点 ' + (A * (1 - e)).toFixed(1) + ' / 远地点 ' + (A * (1 + e)).toFixed(1) + '</span>';
}
function update() {
var val = parseInt(slider.value);
var e = val / 1000;
eDisplay.textContent = e.toFixed(4);
draw(e);
}
slider.addEventListener('input', update);
update();
})();
</script>
</body>
</html>
examples/case04/input/coord.txt
+3 -3
View File
@@ -1,4 +1,4 @@
n mass radius x y z vx vy vz fix_x fix_y fix_z
1 1 0.28 0 0 0 0 0 0 1 1 1
2 1 0.28 4 0 0 0 0 4 0 0 0
3 0.1 0.18 5 0 0 0 0 6 0 0 0
1 27000 1.0 0 0 0 0 0 0 1 1 1
2 81 0.2 10 0 0 0 0 520 0 0 0
3 1 0.1 10.5 0 0 0 0 647 0 0 0
examples/case04/input/input.txt
+2 -2
View File
@@ -23,7 +23,7 @@ force_calc: 0 # 强制重新计算：1=跳过缓存强算，0=自动使用
engine: python # 默认使用 Python 引擎
# ── 盒子 ──────────────────────────────────────
box_a: 20.0 # 立方体半边长，粒子被限制在 [-box_a, box_a]³ 内
box_a: 30.0 # 立方体半边长，粒子被限制在 [-box_a, box_a]³ 内
# ── 初始构型 ──────────────────────────────────
# 坐标文件格式：
@@ -74,7 +74,7 @@ T_total: 10.0
NSTEP: 2
# ── 时间步长 ──────────────────────────────────
DT: 0.001 # 时间步长 (s)
DT: 0.0001 # 时间步长 (s)
# 抽帧范围：只保存 [sample_start, sample_end) 区间内的帧
sample_start: null # null 表示从头开始（帧索引从 0 起）