请求协助排查TE/TM波相移绘图代码返回NaN值的问题

咱们来一步步拆解你遇到的问题，核心是公式适用场景错误+代码执行顺序问题，另外你还混淆了TE/TM波的相移逻辑，下面逐一解决：

1. 为什么会出现NaN？ #

你的代码里计算np.sqrt(np.sin(theta*np.pi/180)**2-n**2)时，n=1.5/1.33≈1.1278，而sinθ的最大值是1，所以sin²θ - n²永远是负数，根号负数直接返回NaN——这是直接原因。

根本问题是：你用了全反射时的相移公式，但你的场景是光从折射率1.33的光疏介质入射到1.5的光密介质，这种情况不会发生全反射（全反射仅当光从光密到光疏时才会出现），所以这个公式完全不适用。

2. 正确的TE/TM波相移规律（光疏→光密介质） #

针对你的场景（n₁=1.33，n₂=1.5）：

• TE波（垂直偏振，s波）：无论入射角是0°到90°，反射波都会有180°（π弧度）的相移（光疏到光密反射时，电场分量直接反向）。
• TM波（平行偏振，p波）：
  • 入射角小于布儒斯特角时，反射波有180°相移；
  • 入射角大于布儒斯特角时，反射波相移为0°；
  • 布儒斯特角θ_B满足tanθ_B = n₂/n₁，计算得θ_B≈48.4°。

3. 修正后的代码 #

先解决执行顺序问题（把x的定义放到函数调用之前），再根据正确的相移逻辑编写代码：

import numpy as np

def phase_te(theta):
    # TE波在光疏→光密介质时，反射相移始终为180°
    return np.full_like(theta, 180.0)

def phase_tm(theta):
    n1 = 1.33
    n2 = 1.5
    theta_rad = np.radians(theta)
    # 计算布儒斯特角
    theta_b = np.degrees(np.arctan(n2 / n1))
    # 初始化相移数组
    shift = np.full_like(theta, 0.0)
    # 入射角小于布儒斯特角时，相移180°
    mask = theta < theta_b
    shift[mask] = 180.0
    return shift

# 生成入射角数组（0°到90°即可，180°和0°对称，超过90°属于反向入射）
x = np.linspace(0, 90, 30)
print("TE波相移：", phase_te(x))
print("TM波相移：", phase_tm(x))

4. 如果实际是光密→光疏场景（比如1.5→1.33） #

如果你的需求其实是光从1.5到1.33（光密到光疏），这时候会发生全反射，临界角θ_c≈62.5°。此时全反射时的TE/TM相移公式如下，修正后的代码：

import numpy as np

def phase_te_total_reflection(theta):
    n1 = 1.5
    n2 = 1.33
    n = n2 / n1
    theta_rad = np.radians(theta)
    # 计算临界角
    theta_c = np.degrees(np.arcsin(n))
    shift = np.zeros_like(theta)
    # 处理全反射区域（入射角>临界角）
    mask_total = theta > theta_c
    sqrt_term = np.sqrt(np.sin(theta_rad[mask_total])**2 - n**2)
    shift[mask_total] = 2 * np.degrees(np.arctan(sqrt_term / np.cos(theta_rad[mask_total])))
    # 非全反射区域，TE波相移180°
    shift[~mask_total] = 180.0
    return shift

def phase_tm_total_reflection(theta):
    n1 = 1.5
    n2 = 1.33
    n = n2 / n1
    theta_rad = np.radians(theta)
    theta_c = np.degrees(np.arcsin(n))
    theta_b = np.degrees(np.arctan(n2 / n1))
    shift = np.zeros_like(theta)
    # 全反射区域的TM相移
    mask_total = theta > theta_c
    sqrt_term = np.sqrt(np.sin(theta_rad[mask_total])**2 - n**2)
    shift[mask_total] = 2 * np.degrees(np.arctan((sqrt_term * n**2) / np.cos(theta_rad[mask_total])))
    # 非全反射区域：小于布儒斯特角相移180°，大于则0°
    mask_b = theta < theta_b
    shift[mask_b & ~mask_total] = 180.0
    return shift

x = np.linspace(0, 90, 30)
print("TE波全反射相移：", phase_te_total_reflection(x))
print("TM波全反射相移：", phase_tm_total_reflection(x))

总结 #

• 先明确入射/出射介质的折射率大小，判断是否会发生全反射；
• 针对不同场景选择对应的相移公式；
• 代码执行顺序要保证变量定义在调用之前。

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