IDL中不借助STREAMLINE函数绘制二维磁场单条磁力线

需求背景 #

需要在xz平面偶极子磁场的等值线图上，叠加绘制一条从指定起点出发的磁力线，且不依赖STREAMLINE函数。现有偶极子场的计算代码如下：

N=20

z=(findgen(N)/(N-1))*10. - 5.
x=z

; 2D arrays
xx=rebin(reform(x, n, 1), n, n)
zz=rebin(reform(z, 1, n), n, n)

rr=sqrt(zz^2 + xx^2)

mu=1.
M=100.
Bx=mu/(4.*!PI)*(3.*M*zz*xx/rr^5)
Bz=mu/(4.*!PI)*(3.*M*zz^2/rr^5 - M/rr^3)

B=sqrt(Bz^2 + Bx^2)

; s=STREAMLINE(Bx, Bz, x, z, XTITLE='X', YTITLE='Z')

核心思路 #

磁力线的切线方向与该点的磁场方向完全一致，因此可以通过迭代步进的方式生成磁力线：

• 从指定起点出发，每一步沿着磁场的单位方向移动一小段距离
• 重复计算新位置的磁场方向，直到超出绘图范围或达到终止条件

完整实现代码 #

以下是在原有代码基础上，添加手动绘制磁力线的完整实现：

N=20

z=(findgen(N)/(N-1))*10. - 5.
x=z

; 2D arrays
xx=rebin(reform(x, n, 1), n, n)
zz=rebin(reform(z, 1, n), n, n)

rr=sqrt(zz^2 + xx^2)

mu=1.
M=100.
Bx=mu/(4.*!PI)*(3.*M*zz*xx/rr^5)
Bz=mu/(4.*!PI)*(3.*M*zz^2/rr^5 - M/rr^3)

B=sqrt(Bz^2 + Bx^2)

; --------------------------
; 手动绘制指定起点的磁力线
; --------------------------
; 设定初始起点（可根据需求修改）
x_start = 2.0
z_start = 3.0
; 步进距离（越小精度越高，计算量越大）
ds = 0.1
; 最大迭代步数（防止无限循环）
max_steps = 500

; 初始化存储磁力线坐标的数组
mag_x = fltarr(max_steps)
mag_z = fltarr(max_steps)
mag_x[0] = x_start
mag_z[0] = z_start

; 欧拉法迭代计算磁力线
for i=1, max_steps-1 do begin
    ; 获取当前坐标
    curr_x = mag_x[i-1]
    curr_z = mag_z[i-1]
    
    ; 计算当前点的磁场分量（直接用公式，避免网格插值误差）
    curr_rr = sqrt(curr_z^2 + curr_x^2)
    if curr_rr eq 0 then break ; 原点磁场无定义，终止迭代
    
    curr_Bx = mu/(4.*!PI)*(3.*M*curr_z*curr_x/curr_rr^5)
    curr_Bz = mu/(4.*!PI)*(3.*M*curr_z^2/curr_rr^5 - M/curr_rr^3)
    
    ; 归一化磁场方向，得到单位向量
    curr_B_mag = sqrt(curr_Bx^2 + curr_Bz^2)
    if curr_B_mag eq 0 then break ; 磁场为0时无方向，终止
    
    dx = ds * (curr_Bx / curr_B_mag)
    dz = ds * (curr_Bz / curr_B_mag)
    
    ; 更新坐标
    mag_x[i] = curr_x + dx
    mag_z[i] = curr_z + dz
    
    ; 判断是否超出绘图范围（x/z范围为-5到5）
    if (mag_x[i] lt -5. or mag_x[i] gt 5. or mag_z[i] lt -5. or mag_z[i] gt 5.) then begin
        ; 截断数组到有效步数
        mag_x = mag_x[0:i]
        mag_z = mag_z[0:i]
        break
    endif
endfor

; 绘制等值线图并叠加磁力线
p = plot(x, z, B, /contour, xtitle='X', ytitle='Z', title='偶极子磁场等值线+手动磁力线')
oplot, mag_x, mag_z, color='red', thick=2

精度优化：龙格-库塔法（RK4） #

欧拉法简单但存在累积误差，若需要更高精度，可使用4阶龙格-库塔法替代上述迭代循环：

; 替换原有的欧拉法循环为RK4实现
for i=1, max_steps-1 do begin
    curr_x = mag_x[i-1]
    curr_z = mag_z[i-1]
    
    curr_rr = sqrt(curr_z^2 + curr_x^2)
    if curr_rr eq 0 then break
    
    ; 计算k1
    k1_x = mu/(4.*!PI)*(3.*M*curr_z*curr_x/curr_rr^5)
    k1_z = mu/(4.*!PI)*(3.*M*curr_z^2/curr_rr^5 - M/curr_rr^3)
    k1_mag = sqrt(k1_x^2 + k1_z^2)
    if k1_mag eq 0 then break
    k1_x = ds * (k1_x / k1_mag)
    k1_z = ds * (k1_z / k1_mag)
    
    ; 计算k2
    x2 = curr_x + k1_x/2
    z2 = curr_z + k1_z/2
    rr2 = sqrt(z2^2 + x2^2)
    k2_x = mu/(4.*!PI)*(3.*M*z2*x2/rr2^5)
    k2_z = mu/(4.*!PI)*(3.*M*z2^2/rr2^5 - M/rr2^3)
    k2_mag = sqrt(k2_x^2 + k2_z^2)
    if k2_mag eq 0 then break
    k2_x = ds * (k2_x / k2_mag)
    k2_z = ds * (k2_z / k2_mag)
    
    ; 计算k3
    x3 = curr_x + k2_x/2
    z3 = curr_z + k2_z/2
    rr3 = sqrt(z3^2 + x3^2)
    k3_x = mu/(4.*!PI)*(3.*M*z3*x3/rr3^5)
    k3_z = mu/(4.*!PI)*(3.*M*z3^2/rr3^5 - M/rr3^3)
    k3_mag = sqrt(k3_x^2 + k3_z^2)
    if k3_mag eq 0 then break
    k3_x = ds * (k3_x / k3_mag)
    k3_z = ds * (k3_z / k3_mag)
    
    ; 计算k4
    x4 = curr_x + k3_x
    z4 = curr_z + k3_z
    rr4 = sqrt(z4^2 + x4^2)
    k4_x = mu/(4.*!PI)*(3.*M*z4*x4/rr4^5)
    k4_z = mu/(4.*!PI)*(3.*M*z4^2/rr4^5 - M/rr4^3)
    k4_mag = sqrt(k4_x^2 + k4_z^2)
    if k4_mag eq 0 then break
    k4_x = ds * (k4_x / k4_mag)
    k4_z = ds * (k4_z / k4_mag)
    
    ; 用RK4公式更新坐标
    mag_x[i] = curr_x + (k1_x + 2*k2_x + 2*k3_x + k4_x)/6
    mag_z[i] = curr_z + (k1_z + 2*k2_z + 2*k3_z + k4_z)/6
    
    ; 判断范围
    if (mag_x[i] lt -5. or mag_x[i] gt 5. or mag_z[i] lt -5. or mag_z[i] gt 5.) then begin
        mag_x = mag_x[0:i]
        mag_z = mag_z[0:i]
        break
    endif
endfor

内容的提问来源于stack exchange，提问作者freja