Seidelの収差論(3次収差論)の理解にむけたメモ書き③-Seidel収差をJuliaで可視化-

Julia

はじめに

本記事は、これまでの連載の続きとなります。幣サイトでは、より原理理解を目指したものとなっているので多少難しい内容となっています。こちらの記事を執筆するにあたって、光学の原理[草川1974]ならびに光学[尾崎2018]、GPT5[chatGPT]を非常に参考にした。

第1回では、瞳の定義や波面収差の基本構造を整理し、収差論の土台をまとめました。第2回では、波面収差のを分析し、その中でSeidel収差がどう構成されるのかを波面収差と光線収差の視点から解説した。

記事番号タイトル概要主な内容
1回目Seidelの収差論①-瞳面と波面収差の立式収差論の基礎として「瞳面」や「波面収差」の定式化(入射/出射瞳、NA、Fナンバーなど)を平易に解説
2回目Seidelの収差論②-光線収差とSeidel収差の各項波面収差を展開し、Seidel収差(球面収差、コマ、非点収差、像面湾曲、歪曲)の導入までを体系的に整理
3回目Seidelの収差論③-Juliaによる可視化これまで導入したSeidel収差をJuliaで数式化し、各収差(球面収差、コマ、非点収差、像面湾曲、歪曲)を可視化。理論と直感的理解をつなげる内容

今まで数式をいじってSeidel収差を理解しようとしてきたが、今のままでは直感的な理解はしづらいと考える。そこで本記事ではJuliaを使った可視化を行っていきたい。

自分の回折とは別視点での説明としてはこれらの記事も読みやすいのでおすすめしたいです。(直感的な理解ではこれらの記事のほうが適している可能性が高い)

文献内容の概要主なポイント
「収差の解説」(me‑no‑koto.com)日常・眼鏡の視点での収差解説• 収差の「実例」紹介(レンズを傾けた像のぼやけや格子の歪み、虹色の滲み)
• 球面収差、コマ収差、非点収差、像面湾曲、歪曲収差、高次収差、色収差などの種類を網羅(眼の事に関する物)
J-Stage 論文(日本視覚科学会誌 36(3), p.40)光学・理論的視点による収差の構造• 理想的な像に関する理論的条件(点→点,平面像,像の形状の保存)
• 「収差(aberrations)」の定義と基礎理論の整理 |
「球面収差とは」(LENS Review)図解・シミュレーションで直感的に説明• 球面収差を、視覚的にわかりやすく図示・シミュレーションを使って丁寧に解説(LENS Review)。球面収差以外も説明。

Seidel収差の数式

座標軸が次の時、

  • 像高ベクトル
    \( (X_0,Y_0), \quad r = \sqrt{X_0^2+Y_0^2}, \quad \theta = \arctan\frac{Y_0}{X_0} \)
  • 瞳面ベクトル
    \( (x,y), \quad \rho = \sqrt{x^2+y^2}, \quad \varphi = \arctan\frac{y}{x}
    \)
  • 内積
    \( (x,y)\cdot(X_0,Y_0) = x X_0 + y Y_0 = \rho r \cos(\varphi-\theta)
    \)

一般的にSeidel収差といわれる波面収差関数は下記の様に書ける。

\begin{equation}
\Phi^{(4)} = A\rho^4+Br \rho^3\cos(\varphi-\theta)+Cr^2\rho^2 \cos^2(\varphi-\theta)+Dr^2\rho^2+E\rho r^3 \cos(\varphi-\theta)
\end{equation}

また、これを光線収差は波面収差に対して\(x, y\)でそれぞれ偏微分すると得ることができ、\( x\)軸上\( y\)軸上以下の様に書ける。

\begin{equation}
\begin{aligned}
\frac{\partial \Phi^{(4)}}{\partial x} &=
4A\,\rho^{2}x + B\bigl(2x\,(\rho r \cos(\varphi-\theta)) + \rho^{2}X_{0}\bigr) + 2C\,(\rho r \cos(\varphi-\theta))\,X_{0} + 2D\,r^{2}x + E\,r^{2}X_{0}\\
\frac{\partial \Phi^{(4)}}{\partial y} &=
4A\,\rho^{2}y + B\bigl(2y\,(\rho r \cos(\varphi-\theta)) + \rho^{2}Y_{0}\bigr) + 2C\,(\rho r \cos(\varphi-\theta))\,Y_{0} +2D\,r^{2}y +E\,r^{2}Y_{0}
\end{aligned}
\end{equation}

波面収差の項光線収差での依存形光線収差名
\( A \rho^4 \)\( \propto \rho^3\)球面収差:瞳座標の3乗に比例。軸上でも残る。
\( B r \rho^3 \cos(\varphi-\theta)\)\( \propto r \rho^2\cos(\varphi-\theta) \)コマ収差:像高に比例し、瞳の2乗依存。典型的なコマ尾像をつくる。
\( C r^2 \rho^2 \cos^2(\varphi-\theta) \)\( \propto r^2 \rho \cos(\varphi-\theta) \)非点収差:像高の2乗に比例、瞳半径に一次比例。メリディオナル/サジタルで焦点がずれる。
\( D r^2 \rho^2 \)\( \propto r^2 \rho \)像面湾曲:像高の2乗に比例し、像面が曲面になる。非点と対になって現れる。
\( E \rho r^3 \cos(\varphi-\theta) \)\( \propto r^3 \cos(\varphi-\theta)\)歪曲:像高の3乗に比例。像面全体の倍率変化(樽型/糸巻型)。

となる。

収差の可視化

収差の可視化を行った図を上に示す。実際のコードは後に示す。
統一的なフォーマットで収差の種類による影響の違いというのをまとめた画像はなかったので、貴重な絵が描けたと考えている。各収差の特徴を次のようになる。

(1) 球面収差 (Spherical Aberration, A)

球面収差は、入射瞳の中心光線と周辺光線で焦点位置が一致しないことに起因する。図においては、各視野点で像が同心円状に広がることが確認できる。これは球面収差が像高に依存せず、軸上像においても顕著に現れることを反映している。

(2) コマ収差 (Coma, B)

コマ収差は、軸外光線の結像位置が瞳位置に依存して異なることにより、非対称な像が生成される現象である。図に示すように、視野中心から離れるに従って星像は彗星状に尾を引き、尾の向きは視野ベクトル方向と一致する。中心像では収差が小さく、周辺において急速に顕著化する。

(3) 非点収差 (Astigmatism, C)

非点収差は、メリディオナル面とサジッタル面で焦点位置が異なることに起因する。図においては、視野方位に応じて像が放射方向あるいは円周方向に線状に伸びる様子が確認できる。この結果は、非点収差が像高依存的であり、視野周辺で顕著となることを示している。

(4) 像面湾曲 (Field Curvature, D)

像面湾曲は、結像面が平面ではなく曲面となることに起因する。図では、視野半径の増大に伴い、各セルに対応する点像が同心的に外方へずれている様子(=円の半径が広がっていくこと)が観察される。これは、平坦な像面を仮定した場合に周辺像が焦点を外れることを表している。

(5) 歪曲収差 (Distortion, E)

歪曲収差は、倍率が像高に依存することにより発生する。図では、理想的な格子配置が樽型あるいは糸巻き型に変形しているが、各点像の形状そのものは保存されている。すなわち、歪曲は像の幾何配置に影響するが、点像の広がりには寄与しない。

Juliaのコード

コードの詳細な解説については省くが、光線収差関数(ray_aberration_xy)と瞳面での像(pupil_rings)を定義してあげてあとは可視化するための細かいコードである。

using PyPlot

struct Coeffs5
    A::Float64  # spherical
    B::Float64  # coma
    C::Float64  # astigmatism
    D::Float64  # field curvature
    E::Float64  # distortion
end

# 光線収差(∂Φ/∂X, ∂Φ/∂Y)on 直交座標
function ray_aberration_xy(x, y, X0, Y0, s::Float64, c::Coeffs5)
    P  = x^2 + y^2       # 平面上の点の距離
    S  = x*X0 + y*Y0     # 像面上の点の距離
    R2 = X0^2 + Y0^2     # 視野中心からの距離
    dΦdX = 4c.A*P*x + c.B*(2x*S + P*X0) + 2c.C*S*X0 + 2c.D*R2*x + c.E*R2*X0
    dΦdY = 4c.A*P*y + c.B*(2y*S + P*Y0) + 2c.C*S*Y0 + 2c.D*R2*y + c.E*R2*Y0
    return s*dΦdX, s*dΦdY
end

# 同心円の瞳サンプリング
function pupil_rings(; ρmax=1.0, Nr=6, Mmax=120)
    xs = Float64[]; ys = Float64[]; ridx = Int[]
    for (k, ρ) in enumerate(range(0, ρmax, length=Nr)[2:end]) # 0を除く
        M = max(12, round(Int, (ρ/ρmax)*Mmax))
        ϕ = range(0, , length=M+1)[1:end-1]
        append!(xs, ρ .* cos.(ϕ))
        append!(ys, ρ .* sin.(ϕ))
        append!(ridx, fill(k, length(ϕ)))
    end
    return xs, ys, ridx
end

# n×n の格子中心
function grid_centers(n::Int; sep::Float64=20.0)
    offs = range(-(n-1)/2*sep, (n-1)/2*sep, length=n)
    centers = Tuple{Float64,Float64}[]
    for cy in offs, cx in offs
        push!(centers, (cx, cy))
    end
    centers
end

# 1セルぶんを、セル中心に応じた視野(X0,Y0)で再計算して配置
function make_aberration(cx, cy; c, s=1.0, field_scale=0.05,
                             ρmax=1.0, Nr=7, Mmax=160)
    # 視野はセル中心に比例して与える(スケールはお好みで)
    X0, Y0 = field_scale*cx, field_scale*cy

    # サンプル点を瞳面上に作る
    xv, yv, _ = pupil_rings(ρmax=ρmax, Nr=Nr, Mmax=Mmax)

    Δx = similar(xv); Δy = similar(yv)
    for i in eachindex(xv)
        Δx[i], Δy[i] = ray_aberration_xy(xv[i], yv[i], X0, Y0, s, c)
    end
    # セル中心へ平行移動して返す
    return cx .+ Δx, cy .+ Δy
end

# まとめて作る
function make_grid_points(n::Int; sep=20.0, kwargs...)
    Xs = Float64[]; Ys = Float64[]
    for (cx, cy) in grid_centers(n; sep=sep)
        x1, y1 = make_aberration(cx, cy; kwargs...)
        append!(Xs, x1); append!(Ys, y1)
    end
    return Xs, Ys
end

# 参照リング(収差ゼロ)
function plot_reference_rings!(ax; center=(0.0,0.0), ρmax=1.0, Nr=5, Mmax=120,
                               color=(0.75,0.75,0.75), alpha=0.6)
    cx, cy = center
    xr, yr, _ = pupil_rings(ρmax=ρmax, Nr=Nr, Mmax=Mmax)
    ax.scatter(cx .+ xr, cy .+ yr; s = 5.0, color=color, alpha=alpha, linewidths=0)
end

# 1パネル描画
function plot_aberration_panel!(ax; name, c::Coeffs5, n=5, sep=22.0, field_scale=0.06,
                                ρmax=1.0, Nr=7, Mmax=180, t=1.0, show_reference=true)
    if show_reference
        for ctr in grid_centers(n; sep=sep)
            plot_reference_rings!(ax; center=ctr, ρmax=ρmax, Nr=Nr, Mmax=Mmax)
        end
    end
    Xs, Ys = make_grid_points(n; sep=sep, c=c, field_scale=field_scale,
                              ρmax=ρmax, Nr=Nr, Mmax=Mmax, s=t)
    ax.scatter(Xs, Ys; s = 7.0, color="red", linewidths=0)

    ax.axis("equal")            # ← ここをオブジェクトメソッドで
    ax.grid(true)
    ax.set_title(name)
    ax.set_xlabel("x axis")
    ax.set_ylabel("y axis")
end

# ===== 5収差を同条件で並べる =====
fig, axs = subplots(2, 3, figsize=(12, 8))
axs = collect(axs)  # 2x3 -> ベクタ化

# 係数(まずは見やすい強さの例)
A = 1.0
B = 0.8
C = 0.8
D = 0.8
E = 0.8
c_spherical   = Coeffs5( A  ,  0.0,   0.0,   0.0,   0.0)  # A>0: 同心状にぼける
c_coma        = Coeffs5( 0.0,  B  ,   0.0,   0.0,   0.0)  # B>0: 視野方向へコメット尾
c_astig       = Coeffs5( 0.0,  0.0,   C  ,   0.0,   0.0)  # C>0: 線状伸長(方位依存)
c_field_curv  = Coeffs5( 0.0,  0.0,   0.0,   D  ,   0.0)  # D>0: 視野^2に比例して全体がずれる
c_distortion  = Coeffs5( 0.0,  0.0,   0.0,   0.0,   E  )  # E>0: 樽/糸巻き

# 共通パラメータ(ここを変えれば全パネル揃って動きます)
n           = 5       # グリッドの星数(n×n)
sep         = 22.0    # 星間隔
field_scale = 0.06    # セル中心 -> 視野スケール
ρmax, Nr, Mmax = 1.0, 7, 180
s           = 1.0     # 画像スケール
show_ref    = true

plot_aberration_panel!(axs[1]; name="Spherical (A = $A)",      c=c_spherical,
    n=n, sep=sep, field_scale=field_scale, ρmax=ρmax, Nr=Nr, Mmax=Mmax, t=s, show_reference=show_ref)

plot_aberration_panel!(axs[2]; name="Coma (B = $B)",           c=c_coma,
    n=n, sep=sep, field_scale=field_scale, ρmax=ρmax, Nr=Nr, Mmax=Mmax, t=s, show_reference=show_ref)

plot_aberration_panel!(axs[3]; name="Astigmatism (C = $C)",    c=c_astig,
    n=n, sep=sep, field_scale=field_scale, ρmax=ρmax, Nr=Nr, Mmax=Mmax, t=s, show_reference=show_ref)

plot_aberration_panel!(axs[4]; name="Field Curvature (D = $D)",c=c_field_curv,
    n=n, sep=sep, field_scale=field_scale, ρmax=ρmax, Nr=Nr, Mmax=Mmax, t=s, show_reference=show_ref)

plot_aberration_panel!(axs[5]; name="Distortion (E = $E)",     c=c_distortion,
    n=n, sep=sep, field_scale=field_scale, ρmax=ρmax, Nr=Nr, Mmax=Mmax, t=s, show_reference=show_ref)

# 6枚目(空き)は凡例的な注釈にしておく
axs[6].axis("off")  
axs[6].text(-0.1, 0.85,
    "gray: reference (no aberration)",
    transform=axs[6].transAxes,
    fontsize=15,
    color="gray")   # ← gray色

# 次に "blue: aberration" の行
axs[6].text(-0.1, 0.75,
    "red: aberration",
    transform=axs[6].transAxes,
    fontsize=15,
    color="red")   # ← blue色

# 残りの説明はまとめて黒などデフォルトで
axs[6].text(-0.1, 0.35,
    "Spherical (A): $A \nComa (B): $B \nAstigmatism (C): $C \nField Curvature (D): $D \nDistortion (E): $E",
    transform=axs[6].transAxes,
    fontsize=12,
    color="black")


tight_layout()
display(gcf())

コメント

タイトルとURLをコピーしました