ECHO3D Wakefield 仿真入门：从最小算例到横向阻抗分析

Posted on In Views: Disqus:
这篇文章完整记录了使用ECHO3D进行 Wakefield 时域仿真的全过程。文章总结了物理概念、命令行流程和后处理方法，并直接附上了本文实际使用过的输入文件与Python 脚本，方便复现。

为什么使用 ECHO3D

ECHO3D 适合做束流驱动下的时域尾场分析，尤其适合研究：

  • 纵向 wake potential $W_\parallel(s)$
  • 横向 wake potential $W_\perp(s)$
  • 纵向阻抗 $Z_\parallel(\omega)$
  • 横向阻抗 $Z_\perp(\omega)$
  • loss factor
  • kick factor

简单说，ECHO3D 研究的是：

  • 束流通过结构后，激发出什么尾场
  • 这些尾场如何影响后续粒子
  • 结构在频域上有哪些显著阻抗峰

本文目标

本文做了两套最小算例：

  1. 轴上束流最小算例：ECHO3D_Minimal
    • 用于得到纵向 wake、纵向阻抗和 loss factor
  2. 偏心束流最小算例：ECHO3D_Dipole
    • 用于得到横向 wake、横向阻抗和 kick factor

同时，还编写了一套 Python 工具链用于：

  • 解析二进制 wake 结果
  • 绘制曲线图
  • 估算阻抗与指标
  • 生成对比汇总表

ECHO3D 命令行工作流

在本文环境中，ECHO3D 的命令行流程是：

1
2
3
4
Mesher.exe input.txt
InitField.exe input.txt
ECHO3D.exe input.txt
IndirectIntegration.exe input.txt

如果在 Windows + cmd 下运行，最稳妥的调用方式类似：

1
cmd /c "cd /d ECHO3D_Minimal && ""Mesher.exe"" input.txt && ""InitField.exe"" input.txt && ""ECHO3D.exe"" input.txt && ""IndirectIntegration.exe"" input.txt"

最小纵向 wake 算例：ECHO3D_Minimal

目录：

  • ECHO3D_Minimal

input.txt

这是最小主输入文件：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
%%%%%%%%%%%%%% minimal ECHO3D input %%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%% geometry %%%%%%%%%%%%%%%%%

GeometryFile = 'geometry.txt'	% [string]
Units = 'mm'	% [m/cm/mm]
BoundaryConditionsX = [0  0]	% [0-open / 1-electric / 2-PML]
BoundaryConditionsY = [0  1]	% [0-magnetic / 1-electric / 2-PML]
BoundaryConditionsZ = [0  1]	% [0-magnetic / 1-electric / 2-PML]

%%%%%%%%%%%%%% beam and field %%%%%%%%%%%%%%%%%

BunchSigma = 1	% [Units]
BunchPosition = [0  0]	% ny nz [mesh steps]
InFieldDir = '-'	% ['free'- free space / '-' -waveguide / string]

%%%%%%%%%%%%%% mesh %%%%%%%%%%%%%%%%%

TimeSteps = -1	% [integer, '-1'- by definition]
MeshLength = 200	% [integer]
dY = [0 10]	% ymin ymax [Units]
dZ = [0 10]	% zmin zmax [Units]
Steps = [0.2 0.2 0.2]	% dx dy dz [Units]
Tolerance = 0.01	% [float, relative]
PMLDepth = 0	% [integer]
PMLParameters = [1 1]	% sigma a [float, relative]

%%%%%%%%%%%%%% solver %%%%%%%%%%%%%%%%%

SolverType = 'impl'	% [expl / impl]
Conformal = 1	% [0-usc / 1-simple]
Iterations = 0	% [integer]
InitialIterations = 0	% [integer]
Damping = 0	% [float, between 0 and 0.5]
ThreadsNumber=4	% [integer]

%%%%%%%%%%%%%% monitors %%%%%%%%%%%%%%%%%

DumpMesh=1	% [0-no/1-yes]
DumpField=1	% [0-no/1-yes]

geometry.txt

这是几何定义文件：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
Background = 0 0 0 

%%%%%%%%%%% Materials %%%%%%%%%%%%%%%%%%

MaterialsNumber = 1
stepout.stl 1 1 0

%%%%%%%%%%% Meshing parts %%%%%%%%%%%%%%

MeshParts = 1
body -0.3 0.3

%%%%%%%%%%% Geometry list %%%%%%%%%%%%%%

GeometryParts = 1
body 1

这个例子使用的是最简单的圆形 step-out 结构 STL。

ECHO3D 瞬态场可视化：正确路线是 FieldMonitor

在本文前面的 wake、阻抗、loss factor 和 kick factor 分析之外，如果还想进一步观察：

  • 束流穿过结构时场分布如何随时间变化
  • 某个切面上的 Ez、Ex 如何传播
  • 如何导出可以直接展示的 GIF

那么更推荐使用 **FieldMonitor**，而不是单纯依赖 DumpField=1 导出的 OutFields/*.bin。

为什么不用 OutFields 作为主要方案

OutFields 的二进制文件当然包含场信息，但它更偏底层数据块，直接用来做瞬态动画时：

  • 文件结构不直观
  • 容易误判维度
  • 很难第一时间确认哪一维是时间、哪一维是空间

相反，FieldMonitor 在官方手册中给出了明确的定义、参数和输出格式，更适合作为“瞬态场观察”的主入口。

官方手册中 FieldMonitor 的关键定义

在 ECHO3D 章节中，FieldMonitor 的定义为：

1
FieldMonitor = { F tF x0 x1 y0 y1 z0 z1 s0 s1 N }

其中：

  • F：场分量，可取 Ex / Ey / Ez / Bx / By / Bz
  • tF：monitor 类型，可取：
    • x：x-time monitor
    • s：s-time monitor

官方手册对 ECHO3D 的说明里还明确指出：

  • 束流只能沿 x 方向传播
  • s = x - ct 是跟着 bunch 一起移动的局部纵向坐标

因此：

  • x-time monitor 更适合看某个 moving-window 里的场如何随 x 演化
  • s-time monitor 更适合看某个相对 bunch 坐标定义下的场如何变化

我在 ECHO3D_Minimal 中采用的 FieldMonitor 配置

为了得到真正可用于瞬态可视化的数据，我最终把最小算例的监视项改成了：

1
2
3
4
FieldMonitor = {'Ez' 	'x' 	20 	100 	0 	10 	0 	10 	0 	100 	1}
FieldMonitor = {'Ez' 	's' 	20 	100 	0 	10 	0 	10 	0 	100 	1}
DumpMesh=1	% [0-no/1-yes]
DumpField=1	% [0-no/1-yes]

这次和最初只做 wake 后处理的版本相比，最大的区别在于：

  • monitor 不再退化成一条线或一个点
  • y、z 两个方向都保留了非零范围
  • 导出的 Monitor_N01.txt / Monitor_N02.txt 已经是真正的大矩阵输出

新的 FieldMonitor 输出结构

重新运行 ECHO3D_Minimal 后，Results/ 目录中出现了两个关键文件：

  • Monitor_N01.txt
  • Monitor_N02.txt

Monitor_N01.txt 头部

1
2
3
4
5
6
7
8
9
10
11
12
13
field Ez
time_mode x
k_ct 105
h_ct 0.0002
ct0 0.0201
k_y 51
h_y 0.0002
y0 0.0
k_z 51
h_z 0.0002
z0 0.0001
k_s 202
h_s 0.0002

Monitor_N02.txt 头部

1
2
3
4
5
6
7
8
9
10
11
12
13
14
field Ez
time_mode s
k_ct 205
h_ct 0.0002
ct0 0.0001
k_y 51
h_y 0.0002
y0 0.0
k_z 51
h_z 0.0002
z0 0.0001
k_x 105
h_x 0.0002
x0 0.02

这两个文件不再是单列坐标，而是超大体积的 ASCII 监视矩阵。

如何正确理解 Monitor_N01

根据实际输出和官方定义，我最终确认：

  • Monitor_N01 对应 time_mode = x
  • 每一行对应一个时间/位置采样
  • 每一行后面的数据可以重构成一个：

$$
51 \times 51 \times 202
$$

的三维数据块。

我进一步固定中间 z 截面，提取出：

$$
E_z(y, s; x)
$$

于是可以得到：

  • 某个固定 z 面上的 y-s 平面图
  • 该平面随 x / 时间推进的动画

结果图：Monitor_N01 第一帧

Monitor N01 first frame

中心点时序曲线

Monitor N01 center series

瞬态 GIF

Monitor N01 official gif

这张 GIF 已经相当接近“束流穿过结构时，固定切面上的场随时间变化”的直观展示。

如何正确理解 Monitor_N02

Monitor_N02 对应：

  • time_mode = s
  • 头部给出 k_x = 105

同样固定中间 z 截面后，我提取出：

$$
E_z(y, x; s)
$$

也就是在某个局部 bunch 坐标条件下，观察 y-x 平面的场分布。

Monitor_N02 第一帧

Monitor N02 first frame

中心点时序曲线

Monitor N02 center series

由于前期有较多弱帧，我额外做了一版精简动画，并把横轴直接标成 x [m]：

x 轴显式标注动画

Monitor N02 x-axis gif

这版动画更适合直接展示“束流沿 x 方向传播时，场在 y-x 平面上的变化”。

从展示效果看：

  • Monitor_N01_official.gif 更像一个 y-s 时空切片动画
  • Monitor_N02_xaxis.gif 更适合用来直观展示束流沿 x 方向传播

Monitor_N01 与 Monitor_N02 的物理区别

可以这样理解：

Monitor_N01_official.gif

更像是在看：

$$
E_z(y,s;x)
$$

适合观察：

  • 尾场在 s 坐标中的分布
  • bunch 头尾附近的时空结构

Monitor_N02_xaxis.gif

更像是在看：

$$
E_z(y,x;s)
$$

适合观察：

  • 沿束流传播方向 x 的场传播
  • 更直观的“束流穿过结构”画面

如果你的目标是做演示或博客展示，我更推荐：

  • 主图用 Monitor_N02_xaxis.gif
  • 补充图用 Monitor_N01_official.gif

这部分使用到的脚本

为了完成上述 FieldMonitor 后处理，我新增了以下脚本：

  • plot_field_monitor_official.py
  • plot_field_monitor_n02_official.py
  • plot_field_monitor_n02_xaxis_gif.py

它们分别用于：

  • 按官方定义解析 Monitor_N01.txt
  • 按官方定义解析 Monitor_N02.txt
  • 生成横轴标为 x [m] 的 GIF

在博客中不再贴出完整脚本全文，建议直接结合仓库文件使用。如果需要展示完整代码，可以按需补充到附录。

关于瞬态可视化的一点经验总结

这次最大的教训是：

  1. DumpField=1 当然能导出底层场数据，但不一定是做动画的最佳入口
  2. 对于 ECHO3D，FieldMonitor 才是更适合做瞬态场可视化的正式接口
  3. 先看官方示例，再决定解析方式，能少走很多弯路

如果后续你想继续优化动画，我建议：

  • 自动裁剪高场区域
  • 给每一帧标注真实 x / ct / s 数值
  • 导出 MP4 版本
  • 再尝试其他切片方式，如固定 y 或固定 s

如果是第一次向外展示结果，我建议优先放这两张动画：

  • Monitor_N01_official.gif
  • Monitor_N02_xaxis.gif

解析 ECHO3D 二进制 wake：parse_wake.py

IndirectIntegration.exe 输出的是二进制文件，例如：

  • wake3D.bin
  • wake3Dindirect.bin

为了把它变成文本曲线，我写了 parse_wake.py：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
from __future__ import annotations

import struct
from pathlib import Path


def parse_wake_file(
    bin_path: Path,
    out_path: Path,
    s0_cm: float = -0.48,
    ds_cm: float = 0.02,
) -> None:
    data = bin_path.read_bytes()

    if len(data) < 44:
        raise ValueError(f"File too small: {bin_path}")

    n = struct.unpack_from('<i', data, 0)[0]
    ny = struct.unpack_from('<i', data, 4)[0]
    nz = struct.unpack_from('<i', data, 8)[0]

    values_per_sample = ny * nz
    header_size = 44
    expected_size = header_size + n * values_per_sample * 8

    if len(data) < expected_size:
        raise ValueError(
            f"Unexpected file size for {bin_path}: "
            f"got {len(data)}, expected at least {expected_size}"
        )

    lines: list[str] = []
    lines.append(f"# Parsed from {bin_path.name}")
    lines.append(f"# header: n={n}, ny={ny}, nz={nz}")
    lines.append("# Columns: index s_cm wake_avg")

    offset = header_size
    for i in range(n):
        vals = struct.unpack_from(f'<{values_per_sample}d', data, offset)
        avg = sum(vals) / values_per_sample
        s_cm = s0_cm + i * ds_cm
        lines.append(f"{i} {s_cm:.6f} {avg:.17g}")
        offset += values_per_sample * 8

    out_path.write_text("\n".join(lines) + "\n", encoding="utf-8")


if __name__ == '__main__':
    base = Path(__file__).resolve().parent
    results = base / 'Results'

    parse_wake_file(
        results / 'wake3Dindirect.bin',
        results / 'wake3Dindirect_parsed.txt',
    )

    print('Done.')
    print(f"Output: {results / 'wake3Dindirect_parsed.txt'}")

运行方式

如果 VS Code 已经绑定到 conda 环境 rl，可以直接运行：

1
python parse_wake.py

输出文件：

  • Results/wake3Dindirect_parsed.txt

如果你的系统里有多个 Python 环境，也可以把上面的命令替换成你自己的解释器路径，但在博客里推荐统一写成：

1
python parse_wake.py

绘制纵向 wake 曲线：plot_wake.py

为了直接把文本结果画成 PNG，又写了 plot_wake.py：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
from __future__ import annotations

from pathlib import Path

import matplotlib.pyplot as plt


def load_wake_data(file_path: Path) -> tuple[list[float], list[float]]:
    s_vals: list[float] = []
    wake_vals: list[float] = []

    for line in file_path.read_text(encoding='utf-8').splitlines():
        line = line.strip()
        if not line or line.startswith('#'):
            continue

        parts = line.split()
        if len(parts) < 3:
            continue

        _, s_cm, wake = parts[:3]
        s_vals.append(float(s_cm))
        wake_vals.append(float(wake))

    return s_vals, wake_vals


def main() -> None:
    base = Path(__file__).resolve().parent
    results_dir = base / 'Results'
    input_file = results_dir / 'wake3Dindirect_parsed.txt'
    output_file = results_dir / 'wake3Dindirect_plot.png'

    s_vals, wake_vals = load_wake_data(input_file)

    if not s_vals:
        raise ValueError(f'No data found in {input_file}')

    plt.figure(figsize=(9, 5.5))
    plt.plot(s_vals, wake_vals, linewidth=2, color='tab:blue', label='Wake (indirect)')
    plt.axhline(0.0, color='black', linewidth=1, linestyle='--')
    plt.axvline(0.0, color='gray', linewidth=1, linestyle=':')
    plt.xlabel('s [cm]')
    plt.ylabel('Wake')
    plt.title('ECHO3D Longitudinal Wake')
    plt.grid(True, linestyle='--', alpha=0.4)
    plt.legend()
    plt.tight_layout()
    plt.savefig(output_file, dpi=160)
    plt.close()

    print('Done.')
    print(f'Input : {input_file}')
    print(f'Output: {output_file}')


if __name__ == '__main__':
    main()

输出图：

  • Results/wake3Dindirect_plot.png

对应图像示意：

ECHO3D Longitudinal Wake

估算纵向阻抗和 loss factor：analyze_longitudinal.py

为了从纵向 wake 进一步估算：

  • 纵向阻抗 $Z_\parallel(\omega)$
  • loss factor

我写了下面这个脚本：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
from __future__ import annotations

from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt


def load_wake(file_path: Path) -> tuple[np.ndarray, np.ndarray]:
    s_vals = []
    w_vals = []
    for line in file_path.read_text(encoding='utf-8').splitlines():
        line = line.strip()
        if not line or line.startswith('#'):
            continue
        parts = line.split()
        if len(parts) < 3:
            continue
        s_vals.append(float(parts[1]))
        w_vals.append(float(parts[2]))
    return np.array(s_vals, dtype=float), np.array(w_vals, dtype=float)


def gaussian_lambda(s_cm: np.ndarray, sigma_cm: float) -> np.ndarray:
    norm = 1.0 / (np.sqrt(2.0 * np.pi) * sigma_cm)
    return norm * np.exp(-(s_cm**2) / (2.0 * sigma_cm**2))


def main() -> None:
    base = Path(__file__).resolve().parent
    results = base / 'Results'
    wake_file = results / 'wake3Dindirect_parsed.txt'

    s_cm, w = load_wake(wake_file)
    ds_cm = float(np.mean(np.diff(s_cm)))
    sigma_cm = 0.1  # BunchSigma=1 mm = 0.1 cm

    lam = gaussian_lambda(s_cm, sigma_cm)
    lam /= np.trapezoid(lam, s_cm)

    # Loss factor approximation: integral lambda(s) * W_parallel(s) ds
    loss_factor = np.trapezoid(lam * w, s_cm)

    # FFT-based longitudinal impedance approximation
    ds_m = ds_cm * 1e-2
    c = 299792458.0
    t_step = ds_m / c
    freqs_hz = np.fft.rfftfreq(len(w), d=t_step)
    z_long = np.fft.rfft(w) * ds_m

    out_txt = results / 'longitudinal_metrics.txt'
    out_txt.write_text(
        '\n'.join([
            '# Longitudinal metrics estimated from wake3Dindirect_parsed.txt',
            f'sigma_cm {sigma_cm}',
            f'ds_cm {ds_cm}',
            f'loss_factor_est {loss_factor:.17g}',
            f'wake_min {w.min():.17g}',
            f'wake_min_s_cm {s_cm[w.argmin()]:.17g}',
            f'wake_max {w.max():.17g}',
            f'wake_max_s_cm {s_cm[w.argmax()]:.17g}',
        ]) + '\n',
        encoding='utf-8'
    )

    imp_txt = results / 'Z_longitudinal_estimated.txt'
    with imp_txt.open('w', encoding='utf-8') as f:
        f.write('# f_Hz ReZ ImZ AbsZ\n')
        for ff, zz in zip(freqs_hz, z_long):
            f.write(f'{ff:.17g} {zz.real:.17g} {zz.imag:.17g} {abs(zz):.17g}\n')

    plt.figure(figsize=(9, 5.5))
    plt.plot(freqs_hz * 1e-9, np.abs(z_long), linewidth=1.5)
    plt.xlabel('f [GHz]')
    plt.ylabel('|Z_parallel| [arb.]')
    plt.title('Estimated Longitudinal Impedance')
    plt.grid(True, linestyle='--', alpha=0.4)
    plt.tight_layout()
    plt.savefig(results / 'Z_longitudinal_estimated.png', dpi=160)
    plt.close()

    print('Done.')
    print(f'Loss factor estimate: {loss_factor:.6g}')
    print(f'Output: {out_txt}')
    print(f'Output: {imp_txt}')


if __name__ == '__main__':
    main()

ECHO3D_Minimal 结果

对于轴上束流最小算例，得到：

  • loss_factor_est = -14.636628710029601
  • wake_min = -21.328521756987332
  • wake_min_s_cm = 0
  • wake_max = 3.3978856597281952
  • wake_max_s_cm = 0.4

这表明：

  • 主负峰出现在 $s=0$
  • 后续存在正向回摆和衰减振荡

对应纵向阻抗图：

Estimated Longitudinal Impedance

构造偏心束流横向算例：ECHO3D_Dipole

目录：

  • ECHO3D_Dipole

这个算例是在 ECHO3D_Minimal 基础上复制出来，并把束流改为偏心注入。

input.txt

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
%%%%%%%%%%%%%% minimal ECHO3D input %%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%% geometry %%%%%%%%%%%%%%%%%

GeometryFile = 'geometry.txt'	% [string]
Units = 'mm'	% [m/cm/mm]
BoundaryConditionsX = [0  0]	% [0-open / 1-electric / 2-PML]
BoundaryConditionsY = [0  1]	% [0-magnetic / 1-electric / 2-PML]
BoundaryConditionsZ = [0  1]	% [0-magnetic / 1-electric / 2-PML]

%%%%%%%%%%%%%% beam and field %%%%%%%%%%%%%%%%%

BunchSigma = 1	% [Units]
BunchPosition = [1  0]	% ny nz [mesh steps]
InFieldDir = '-'	% ['free'- free space / '-' -waveguide / string]

%%%%%%%%%%%%%% mesh %%%%%%%%%%%%%%%%%

TimeSteps = -1	% [integer, '-1'- by definition]
MeshLength = 200	% [integer]
dY = [0 10]	% ymin ymax [Units]
dZ = [0 10]	% zmin zmax [Units]
Steps = [0.2 0.2 0.2]	% dx dy dz [Units]
Tolerance = 0.01	% [float, relative]
PMLDepth = 0	% [integer]
PMLParameters = [1 1]	% sigma a [float, relative]

%%%%%%%%%%%%%% solver %%%%%%%%%%%%%%%%%

SolverType = 'impl'	% [expl / impl]
Conformal = 1	% [0-usc / 1-simple]
Iterations = 0	% [integer]
InitialIterations = 0	% [integer]
Damping = 0	% [float, between 0 and 0.5]
ThreadsNumber=4	% [integer]

%%%%%%%%%%%%%% monitors %%%%%%%%%%%%%%%%%

DumpMesh=0	% [0-no/1-yes]
DumpField=0	% [0-no/1-yes]
WakeMonitor= [-1 1 -1 1]	% ymin ymax zmin zmax [in mesh lines from the bunch position]

和轴上束流相比，关键改动是：

1
BunchPosition = [1 0]

这会激发横向偶极响应。

估算横向 wake、横向阻抗和 kick factor：analyze_transverse.py

为了从 3×2 的 monitor 网格里构造横向 wake，我写了如下脚本：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
from __future__ import annotations

import struct
from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt


def parse_grid(bin_path: Path, s0_cm: float = -0.48, ds_cm: float = 0.02):
    data = bin_path.read_bytes()
    n = struct.unpack_from('<i', data, 0)[0]
    ny = struct.unpack_from('<i', data, 4)[0]
    nz = struct.unpack_from('<i', data, 8)[0]
    header_size = 44
    vals_per_sample = ny * nz
    offset = header_size
    arr = np.empty((n, ny, nz), dtype=float)
    for i in range(n):
        vals = struct.unpack_from(f'<{vals_per_sample}d', data, offset)
        arr[i, :, :] = np.array(vals, dtype=float).reshape((ny, nz))
        offset += vals_per_sample * 8
    s = s0_cm + np.arange(n, dtype=float) * ds_cm
    return s, arr


def main() -> None:
    base = Path(__file__).resolve().parent
    results = base / 'Results'
    s_cm, arr = parse_grid(results / 'wake3Dindirect.bin')

    ny, nz = arr.shape[1], arr.shape[2]
    # average over z first
    yz_avg = arr.mean(axis=2)  # shape (n, ny)

    # longitudinal component: average over all monitor points
    w_long = yz_avg.mean(axis=1)

    # transverse estimate: left-right finite difference if possible
    if ny >= 3:
        left = yz_avg[:, 0]
        center = yz_avg[:, 1]
        right = yz_avg[:, 2]
        # bunch offset is 1 mesh step in y, mesh step is 0.2 mm = 0.02 cm
        dy_cm = 0.02
        w_trans = (right - left) / (2.0 * dy_cm)
    elif ny == 2:
        dy_cm = 0.02
        w_trans = (yz_avg[:, 1] - yz_avg[:, 0]) / dy_cm
    else:
        dy_cm = 0.02
        w_trans = np.zeros_like(w_long)

    sigma_cm = 0.1
    lam = (1.0 / (np.sqrt(2.0 * np.pi) * sigma_cm)) * np.exp(-(s_cm**2) / (2.0 * sigma_cm**2))
    lam /= np.trapezoid(lam, s_cm)

    kick_factor = np.trapezoid(lam * w_trans, s_cm)

    ds_m = float(np.mean(np.diff(s_cm))) * 1e-2
    c = 299792458.0
    dt = ds_m / c
    freqs_hz = np.fft.rfftfreq(len(w_trans), d=dt)
    z_trans = np.fft.rfft(w_trans) * ds_m

    out_txt = results / 'transverse_metrics.txt'
    out_txt.write_text(
        '\n'.join([
            '# Transverse metrics estimated from wake3Dindirect.bin',
            f'ny {ny}',
            f'nz {nz}',
            f'kick_factor_est {kick_factor:.17g}',
            f'wake_trans_min {w_trans.min():.17g}',
            f'wake_trans_min_s_cm {s_cm[w_trans.argmin()]:.17g}',
            f'wake_trans_max {w_trans.max():.17g}',
            f'wake_trans_max_s_cm {s_cm[w_trans.argmax()]:.17g}',
        ]) + '\n',
        encoding='utf-8'
    )

    wake_txt = results / 'wake_transverse_estimated.txt'
    with wake_txt.open('w', encoding='utf-8') as f:
        f.write('# index s_cm W_trans_est\n')
        for i, (ss, ww) in enumerate(zip(s_cm, w_trans)):
            f.write(f'{i} {ss:.6f} {ww:.17g}\n')

    z_txt = results / 'Z_transverse_estimated.txt'
    with z_txt.open('w', encoding='utf-8') as f:
        f.write('# f_Hz ReZ ImZ AbsZ\n')
        for ff, zz in zip(freqs_hz, z_trans):
            f.write(f'{ff:.17g} {zz.real:.17g} {zz.imag:.17g} {abs(zz):.17g}\n')

    plt.figure(figsize=(9, 5.5))
    plt.plot(s_cm, w_trans, linewidth=2, color='tab:red')
    plt.axhline(0.0, color='black', linewidth=1, linestyle='--')
    plt.axvline(0.0, color='gray', linewidth=1, linestyle=':')
    plt.xlabel('s [cm]')
    plt.ylabel('W_perp est.')
    plt.title('Estimated Transverse Wake')
    plt.grid(True, linestyle='--', alpha=0.4)
    plt.tight_layout()
    plt.savefig(results / 'wake_transverse_estimated.png', dpi=160)
    plt.close()

    plt.figure(figsize=(9, 5.5))
    plt.plot(freqs_hz * 1e-9, np.abs(z_trans), linewidth=1.5, color='tab:green')
    plt.xlabel('f [GHz]')
    plt.ylabel('|Z_perp| [arb.]')
    plt.title('Estimated Transverse Impedance')
    plt.grid(True, linestyle='--', alpha=0.4)
    plt.tight_layout()
    plt.savefig(results / 'Z_transverse_estimated.png', dpi=160)
    plt.close()

    print('Done.')
    print(f'Kick factor estimate: {kick_factor:.6g}')
    print(f'Output: {out_txt}')
    print(f'Output: {wake_txt}')
    print(f'Output: {z_txt}')


if __name__ == '__main__':
    main()

ECHO3D_Dipole 结果

得到：

  • kick_factor_est = 0.4717549321159662
  • wake_trans_min = -64.800902043818539
  • wake_trans_min_s_cm = -0.1
  • wake_trans_max = 72.744079625586437
  • wake_trans_max_s_cm = 0.1

这说明偏心束流已经激发出明显的横向 wake 和横向阻抗。

对应横向 wake 图：

Estimated Transverse Wake

对应横向阻抗图：

Estimated Transverse Impedance

对比两套算例：compare_echo3d_results.py

为了做总表和对比图，我写了：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
from __future__ import annotations

from pathlib import Path
import math
import matplotlib.pyplot as plt

ROOT = Path(__file__).resolve().parent
CASES = {
    'Minimal (on-axis)': ROOT / 'ECHO3D_Minimal' / 'Results',
    'Dipole (offset)': ROOT / 'ECHO3D_Dipole' / 'Results',
}


def parse_key_value_file(path: Path) -> dict[str, float | str]:
    result: dict[str, float | str] = {}
    for line in path.read_text(encoding='utf-8').splitlines():
        line = line.strip()
        if not line or line.startswith('#'):
            continue
        parts = line.split(maxsplit=1)
        if len(parts) != 2:
            continue
        key, value = parts
        try:
            result[key] = float(value)
        except ValueError:
            result[key] = value
    return result


def load_curve(path: Path, s_index: int, v_index: int):
    xs = []
    ys = []
    for line in path.read_text(encoding='utf-8').splitlines():
        line = line.strip()
        if not line or line.startswith('#'):
            continue
        parts = line.split()
        xs.append(float(parts[s_index]))
        ys.append(float(parts[v_index]))
    return xs, ys


def fmt_num(v: float | str | None) -> str:
    if v is None:
        return 'N/A'
    if isinstance(v, str):
        return v
    if math.isfinite(v):
        return f'{v:.6g}'
    return str(v)


def main() -> None:
    summary_lines = ['# ECHO3D comparison summary', '']
    rows = []

    for name, results_dir in CASES.items():
        long_metrics = parse_key_value_file(results_dir / 'longitudinal_metrics.txt')
        trans_metrics_path = results_dir / 'transverse_metrics.txt'
        trans_metrics = parse_key_value_file(trans_metrics_path) if trans_metrics_path.exists() else {}

        row = {
            'case': name,
            'loss_factor': long_metrics.get('loss_factor_est'),
            'wake_long_min': long_metrics.get('wake_min'),
            'wake_long_min_s': long_metrics.get('wake_min_s_cm'),
            'wake_long_max': long_metrics.get('wake_max'),
            'wake_long_max_s': long_metrics.get('wake_max_s_cm'),
            'kick_factor': trans_metrics.get('kick_factor_est', 0.0),
            'wake_trans_min': trans_metrics.get('wake_trans_min', 0.0),
            'wake_trans_min_s': trans_metrics.get('wake_trans_min_s_cm', 0.0),
            'wake_trans_max': trans_metrics.get('wake_trans_max', 0.0),
            'wake_trans_max_s': trans_metrics.get('wake_trans_max_s_cm', 0.0),
        }
        rows.append(row)

    summary_lines.append('| Case | loss factor | W_long min | s at min [cm] | W_long max | s at max [cm] | kick factor | W_trans min | s at trans min [cm] | W_trans max | s at trans max [cm] |')
    summary_lines.append('|---|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|')
    for r in rows:
        summary_lines.append(
            f"| {r['case']} | {fmt_num(r['loss_factor'])} | {fmt_num(r['wake_long_min'])} | {fmt_num(r['wake_long_min_s'])} | {fmt_num(r['wake_long_max'])} | {fmt_num(r['wake_long_max_s'])} | {fmt_num(r['kick_factor'])} | {fmt_num(r['wake_trans_min'])} | {fmt_num(r['wake_trans_min_s'])} | {fmt_num(r['wake_trans_max'])} | {fmt_num(r['wake_trans_max_s'])} |"
        )

    summary_path = ROOT / 'ECHO3D_summary.md'
    summary_path.write_text('\n'.join(summary_lines) + '\n', encoding='utf-8')

    # Plot longitudinal wake comparison
    plt.figure(figsize=(10, 6))
    for name, results_dir in CASES.items():
        x, y = load_curve(results_dir / 'wake3Dindirect_parsed.txt', 1, 2)
        plt.plot(x, y, linewidth=2, label=name)
    plt.axhline(0.0, color='black', linewidth=1, linestyle='--')
    plt.axvline(0.0, color='gray', linewidth=1, linestyle=':')
    plt.xlabel('s [cm]')
    plt.ylabel('W_parallel')
    plt.title('Longitudinal Wake Comparison')
    plt.grid(True, linestyle='--', alpha=0.4)
    plt.legend()
    plt.tight_layout()
    plt.savefig(ROOT / 'ECHO3D_compare_longitudinal.png', dpi=160)
    plt.close()

    # Plot transverse wake comparison where available
    plt.figure(figsize=(10, 6))
    plotted = False
    for name, results_dir in CASES.items():
        path = results_dir / 'wake_transverse_estimated.txt'
        if path.exists():
            x, y = load_curve(path, 1, 2)
            plt.plot(x, y, linewidth=2, label=name)
            plotted = True
    if plotted:
        plt.axhline(0.0, color='black', linewidth=1, linestyle='--')
        plt.axvline(0.0, color='gray', linewidth=1, linestyle=':')
        plt.xlabel('s [cm]')
        plt.ylabel('W_perp est.')
        plt.title('Transverse Wake Comparison')
        plt.grid(True, linestyle='--', alpha=0.4)
        plt.legend()
        plt.tight_layout()
        plt.savefig(ROOT / 'ECHO3D_compare_transverse.png', dpi=160)
    plt.close()

    print('Done.')
    print(f'Summary: {summary_path}')
    print(f'Figure : {ROOT / "ECHO3D_compare_longitudinal.png"}')
    print(f'Figure : {ROOT / "ECHO3D_compare_transverse.png"}')


if __name__ == '__main__':
    main()

本文最终结果总表

汇总文件：

  • ECHO3D_summary.md

内容如下：

1
2
3
4
5
6
# ECHO3D comparison summary

| Case | loss factor | W_long min | s at min [cm] | W_long max | s at max [cm] | kick factor | W_trans min | s at trans min [cm] | W_trans max | s at trans max [cm] |
|---|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:|
| Minimal (on-axis) | -14.6366 | -21.3285 | 0 | 3.39789 | 0.4 | 0 | 0 | 0 | 0 | 0 |
| Dipole (offset) | -29.2804 | -42.5884 | 0 | 6.79366 | 0.4 | 0.471755 | -64.8009 | -0.1 | 72.7441 | 0.1 |

纵向 wake 对比图：

Longitudinal Wake Comparison

横向 wake 对比图：

Transverse Wake Comparison

如何理解这些结果

纵向 wake potential $W_\parallel(s)$

  • Minimal (on-axis) 已经有明显纵向 wake
  • Dipole (offset) 的纵向负峰更深，说明偏心激励下结构响应更强

横向 wake potential $W_\perp(s)$

  • 轴上束流几乎没有横向响应
  • 偏心束流后，横向 wake 明显非零

纵向阻抗 $Z_\parallel(\omega)$

通过对 $W_\parallel(s)$ 做 FFT 得到近似阻抗谱，用于查看结构在不同频率上的纵向响应。

横向阻抗 $Z_\perp(\omega)$

通过对估算的横向 wake 做 FFT 得到近似横向阻抗谱。

loss factor

当前结果：

  • 轴上束流：-14.6366
  • 偏心束流：-29.2804

说明偏心激励时，结构导致的等效能量损失更强。

kick factor

当前结果：

  • 轴上束流：0
  • 偏心束流：0.471755

这表示偏心束流通过结构后，对后续粒子存在非零横向踢偏作用。

全部推荐运行顺序

跑最小纵向 wake 算例

1
2
3
4
5
6
7
Mesher.exe input.txt
InitField.exe input.txt
ECHO3D.exe input.txt
IndirectIntegration.exe input.txt
python parse_wake.py
python plot_wake.py
python analyze_longitudinal.py

跑偏心横向 wake 算例

1
2
3
4
5
6
7
Mesher.exe input.txt
InitField.exe input.txt
ECHO3D.exe input.txt
IndirectIntegration.exe input.txt
python parse_wake.py
python analyze_longitudinal.py
python analyze_transverse.py

生成总表和对比图

1
python compare_echo3d_results.py

输出文件总览

ECHO3D_Minimal/Results

  • wake3Dindirect_parsed.txt
  • wake3Dindirect_plot.png
  • longitudinal_metrics.txt
  • Z_longitudinal_estimated.txt
  • Z_longitudinal_estimated.png

ECHO3D_Dipole/Results

  • wake3Dindirect_parsed.txt
  • wake_transverse_estimated.txt
  • wake_transverse_estimated.png
  • transverse_metrics.txt
  • Z_transverse_estimated.txt
  • Z_transverse_estimated.png
  • longitudinal_metrics.txt
  • Z_longitudinal_estimated.txt

项目根目录

  • compare_echo3d_results.py
  • ECHO3D_summary.md
  • ECHO3D_compare_longitudinal.png
  • ECHO3D_compare_transverse.png

这套流程的局限性

需要说明的是，本文的很多后处理量是“工程上可用的近似估算”，尤其包括：

  • 二进制 wake 解析
  • 纵向阻抗 FFT
  • 横向 wake 差分构造
  • 横向阻抗估算
  • loss factor / kick factor 估算

这些方法适合作为：

  • 学习 ECHO3D 工作流
  • 快速建立物理直觉
  • 做最小验证

如果后续要进入高精度工程设计，建议进一步：

  • 查阅 ECHO3D 官方格式说明
  • 使用更严格的后处理方法
  • 对比更多 benchmark
  • 采用更真实的腔体 STL 几何

下一步建议

完成这篇文章后，下一步最值得做的事情有三种：

  1. 把当前 step-out 几何替换成自己的 pillbox / 带束管腔体 STL
  2. 进一步自动化：把 Mesher → InitField → ECHO3D → IndirectIntegration → Python 后处理 串成一键脚本
  3. 在结果图上增加：
    • 峰值标注
    • 更多物理单位说明
    • 阻抗峰识别

结语

这套工作流的核心收获是：

  • 不只是“能跑一个 ECHO3D 示例”
  • 而是已经建立起一套完整的命令行 wakefield 仿真与后处理框架

从这里出发，就可以逐步把最小算例替换成真实腔体结构，进入真正的加速器结构 wakefield 分析。