torcwaでメタサーフェイス・カラースプリッターをシミュレーションする

Python

はじめに

本記事ではtorcwaを用いた実際のsimulationについて解説を行うとともに、自分がclaude codeを使った初めての開発だったため、気を付ける点について触れる。

torcwaというのはソウル大学開発のRCWA計算シミュレータである[changhyun2023]。samsungでの資金提供を受けての開発であるため、すでに製品化されているsamsung開発のカラースプリッター(nano prism)(参考:ISO CELL JNPについての日本語記事)の開発にも利用されていることが想像できる。

また、今回のシミュレーションでは[Xiujuan2022]の研究をもとにシミュレーションを開発した。この研究は2×2um画素を16×16 pixelに分割して、これを遺伝的アルゴリズムを用いて最適化していくといった研究となっている。詳細については割愛するが最適化計算と自由形状設計(pix分割)といった意味で現在の開発の礎の一つとなっている研究だと認識している。

カラースプリッタ開発一般論においてはセンサー画素が平面に並んでいること(→RCWAにおける周期境界条件を満たしやすい)と最適化計算においては一回の計算時間を短くしたほうがよいといった理由からRCWAのほうが好まれる傾向にあるといえる。(実際のCMOS構造に含まれる金属構造を考慮するとRCWAは最適であるとはいいきれないが…)

一方、この論文の計算[Xiujuan2022]はFDTDによる計算となっており、この論文とRCWAの差異を確認するために本シミュレーションには意義があると考える。

また、今回のシミュレーションにあたってはclaude codeによる支援とサイト執筆にあたってはclaude opus4.6による添削を経ている。

ファイル構成は以下のようになっている。

Xiujuan2022/
├── paper_data/                    # 論文から取得した元データ
│   ├── s41467-022-31019-7.pdf
│   ├── spec.txt                   # 変えてほしくないパラメータを事前にまとめる
│   ├── nk_inf.png                 # Si₃N₄ の n,k スペクトル図
│   └── shape_data.png             # 16×16 ピラー配置図
│
├── read_nk_info.py                # nk_inf.png → data/nk_info.csv
├── read_shape_data.py             # shape_data.png → data/shape_grid.csv
│
├── data/
│   ├── nk_info.csv                # 波長 vs n, k テーブル
│   ├── shape_grid.csv             # 16×16 バイナリマスク
│   └── torcwa_output/             # シミュレーション出力（図・CSV）
│
├── reference_data/
│   └── results.png                # 論文 Fig.2 との比較用
│
└── torcwa_colorsplitter.ipynb     # ★ メインnotebook（セル1〜9）

Claude codeの準備

最初にclaude codeを書くときに、readme.txtを書いたほうがいいってどこかで見たので絶対守ってほしいところをリストアップした。

以下のルールを必ず守ってください

・データを読むのは人間です。コードは最低限に書いて、コメントも明記してください。
・関数化はdefのみで"複数回使うとき"以外は関数化しないでください。
・spec.md に書かれた物理条件（寸法・波長・材料・構造）は絶対に変更しない
・不明点がある場合は補完せず、必ず TODO またはコメントで止める
・勝手に「より良い条件」や「一般化」をしない
・RCWA / torcwa の計算フローを勝手に変更しない
・境界条件（periodic, normal incidence など）を変更しない
・材料モデル（n,k → ε変換など）を簡略化しない
・波長依存性を定数近似に置き換えない
・論文由来データと、コード内で生成したデータを明確に分ける
・CSVやmaskは「data/」に置き、コード内に埋め込まない
・「これは論文由来」「これは補間・近似」とコメントで明記する
・出力フォーマット（csv）を統一
・電場は |E|² か複素Eかを明示
・軸（x,y,z）の定義を必ずコメント

次に絶対変えてほしくない物理パラメータをファイルにまとめた。これをspec.txtとしてまとめた。

unit cell: 2 µm × 2 µm
16 × 16 binary grid
1 cell size: 125 nm
pillar height: 600 nm
material: Si3N4 / air
substrate: quartz
diagonal symmetry は論文の最適化条件だが、今回は補足図(paper_data/*.png)の最終形状をそのまま使う
評価対象は透過後の電場分布

また、物理パラメータ解析用に論文データの一部を切り取って、nk_inf.pngとshape_data.pngとして保存した。

torcwaの基本

torcwaの基本的な使い方やinstall方法については過去記事のこちらに記載してある。

最低限のコードとしては以下のような流れとなっている。


# 0. ハードウェア設定とmaterial情報
torch.backends.cuda.matmul.allow_tf32 = False
sim_dtype = torch.complex64
geo_dtype = torch.float32
device = torch.device('cuda') #CPUの場合はcpuを選択

# material
substrate_eps = 1.45**2
air = 1.0**2
a_Si = 3.55**2
lam = 1064. # nm

# geometry
L = [500., 500.]            # nm / nm
torcwa.rcwa_geo.dtype = geo_dtype
torcwa.rcwa_geo.device = device
torcwa.rcwa_geo.Lx = L[0]
torcwa.rcwa_geo.Ly = L[1]
torcwa.rcwa_geo.nx = 500
torcwa.rcwa_geo.ny = 500
torcwa.rcwa_geo.grid()
torcwa.rcwa_geo.edge_sharpness = 1000.
z = torch.linspace(-500,1500,101,device=device)

fn = 7 #fourier次数
# 1. シミュレーション宣言
sim = torcwa.rcwa(freq=1/lam, order=[fn,fn], L=[Lx,Ly], dtype=..., device=...)

# 2. 入出力層の追加（両方が自由空間なら省略可）
sim.add_input_layer(eps=substrate_eps)
sim.add_output_layer(eps=air)       # 省略時はε=1

# 3. 入射角の設定
sim.set_incident_angle(inc_ang=0., azi_ang=0.)

# 4. 内部層の追加（add_layerの呼び出し順 = 物理的な積層順）
sim.add_layer(thickness=t_pillar, eps=eps_map)    # layer 0, eps_mapは別途作成
sim.add_layer(thickness=t_air, eps=air)       # layer 1

# 5. グローバルSマトリクスの求解
sim.solve_global_smatrix()

# 6. Sパラメータ取得（電場不要な場合はここまで）
txx = sim.S_parameters(orders=[0,0], direction='forward',
                       port='transmission', polarization='xx',
                       ref_order=[0,0])

# 7. 光源設定（電場を見る場合）
sim.source_planewave(amplitude=[1.,0.], direction='forward', notation='xy')

# 8. 電場取得
[Ex,Ey,Ez],[Hx,Hy,Hz] = sim.field_xz(x_ax, z_ax, y_cut)  # XZ断面
[Ex,Ey,Ez],[Hx,Hy,Hz] = sim.field_xy(layer_idx, x_ax, y_ax, z_prop=0.)  # XY面

自分がsimulationをする中でよく失敗していた点として、

・フーリエ次数についてはよく吟味すること。(後述する手法参考)

・電場分布取得のlayer indexの指定

・屈折率をそのまま入力する(複素誘電率がinputなので、複素屈折率の2乗を入力する)

あとは地味にlayerの積層を間違えることが多く、return_layerを使うことでそのlayerがどんなlayerかと実際の計算にどんな誘電率分布で使われているかを確認できる。

eps_recover, mu_recover = sim.return_layer(layer_num, nx=100, ny=100)

引数:

  • layer_num: add_layer した順のindex（上記例ならピラー層なら0、その上の空気層なら1）
  • nx, ny: 逆変換後の実空間グリッド数（デフォルト100×100）

戻り値:

  • eps_recover: shape (nx, ny) の complex tensor。実空間での誘電率分布（Fourier再構成）
  • mu_recover: 同じく透磁率分布（通常の非磁性材料なら全面1）

※ここでnx, nyというのは出力の画素数なので計算部分と干渉しない

構造についてのリバエン

構造についての逆解析は基本的にはclaude codeに書いてもらった。ここで指示の内容としては、まずは論文をpdfファイルにして読ませてグラフの情報を尋ねてみる。

claude code prompt

fig. xについてxの最小/最大値, yの最小/最大値はどういった数値になりますか？
周辺データを含めながら、このグラフの具体値を用いた解析をして下さい。

ある程度妥当な回答が返ってきた場合はそのまま進めてもらう。このときに、解析してほしい画像をスクショして別ファイル(自分は新しく\paper_dataというファイルを作った)に入れるとうまくいく。また、具体的な手法を説明したほうがいい。

claude code prompt

今回はpaper_data\○○というファイル名のグラフデータを解析して数値データにするPythonファイルを作ってほしいです。

このグラフデータ解析には色によってグラフのデータだけを抜き出して曲線のpixel位置を把握してください。また、グラフの軸を確認して、最大値/最小値を確認したうえでグラフの枠線を利用したpixelから物理数値に変換することを考えてください。

もし不明点や読み取れない場合があったら、自分で考えて決めずに聞いてください

以上のような説明をしても、正直出てくる解析方法は80点どまりで済んでしまう。そのときは自分でコードを変更したり、直接データをいじったほうがはやい。(AIに固執してもあまりうまくいかないことが多い。)

今回の論文の物性値を入手したコードを以下に示す。

read_nk_info.py

pythonファイル

# read_nk_info.py
import csv
from pathlib import Path

from PIL import Image


ROOT = Path(__file__).resolve().parent
IMAGE_PATH = ROOT / "paper_data" / "nk_inf.png"
OUTPUT_PATH = ROOT / "data" / "nk_info.csv"

# これは論文図の軸レンジ。paper_data/nk_inf.png から読む。
WAVELENGTH_MIN_NM = 400.0
WAVELENGTH_MAX_NM = 700.0
N_MIN = 1.96
N_MAX = 2.10
K_MIN = -0.002
K_MAX = 0.016


def interpolate(xs, ys, x):
    if x <= xs[0]:
        return ys[0]
    if x >= xs[-1]:
        return ys[-1]

    for i in range(1, len(xs)):
        if x <= xs[i]:
            x0 = xs[i - 1]
            x1 = xs[i]
            y0 = ys[i - 1]
            y1 = ys[i]
            return y0 + (y1 - y0) * (x - x0) / (x1 - x0)

    return ys[-1]


image = Image.open(IMAGE_PATH).convert("RGB")
width, height = image.size
pixels = image.load()

# nk_inf.png は右軸が青なので、黒線だけでなく青線も枠線候補に使う。
row_segments = []
for y in range(int(height * 0.85)):
    xs = []
    for x in range(width):
        r, g, b = pixels[x, y]
        is_dark = max(r, g, b) < 100
        is_blue = b > 150 and r < 120 and g < 140
        if is_dark or is_blue:
            xs.append(x)

    if not xs:
        continue

    start = xs[0]
    prev = xs[0]
    segments = []
    for x in xs[1:]:
        if x == prev + 1:
            prev = x
            continue
        segments.append((start, prev))
        start = x
        prev = x
    segments.append((start, prev))

    longest_start, longest_end = max(segments, key=lambda seg: seg[1] - seg[0])
    longest_span = longest_end - longest_start
    if longest_span > width * 0.5:
        row_segments.append((y, longest_start, longest_end))

if not row_segments:
    raise RuntimeError("TODO: プロット枠を検出できませんでした。paper_data/nk_inf.png を確認してください。")

top_y, left_x, right_x = row_segments[0]

bottom_rows = []
for y, row_left, row_right in row_segments:
    if y <= top_y + 100:
        continue
    if abs(row_left - left_x) <= 5 and abs(row_right - right_x) <= 12:
        bottom_rows.append((y, row_left, row_right))

if not bottom_rows:
    raise RuntimeError("TODO: 下側の軸を検出できませんでした。paper_data/nk_inf.png を確認してください。")

bottom_y, left_x, right_x = bottom_rows[-1]

# これは論文図の画像解析。黒線を n、青線を k として読む。
n_trace = {}
k_trace = {}
prev_n_y = None
prev_k_y = None

for x in range(left_x + 3, right_x):
    n_candidates = []
    k_candidates = []

    for y in range(top_y + 3, bottom_y - 3):
        r, g, b = pixels[x, y]

        if max(r, g, b) < 80 and max(r, g, b) - min(r, g, b) < 25:
            n_candidates.append(y)

        if b > 150 and r < 80 and g < 120:
            k_candidates.append(y)

    if n_candidates:
        if prev_n_y is None:
            n_trace[x] = min(n_candidates)
        else:
            n_trace[x] = min(n_candidates, key=lambda y: abs(y - prev_n_y))
        prev_n_y = n_trace[x]

    if k_candidates:
        if prev_k_y is None:
            k_trace[x] = min(k_candidates)
        else:
            k_trace[x] = min(k_candidates, key=lambda y: abs(y - prev_k_y))
        prev_k_y = k_trace[x]

if not n_trace:
    raise RuntimeError("TODO: n の曲線を読めませんでした。paper_data/nk_inf.png を確認してください。")

if not k_trace:
    raise RuntimeError("TODO: k の曲線を読めませんでした。paper_data/nk_inf.png を確認してください。")

n_pixels = sorted(n_trace)
k_pixels = sorted(k_trace)

n_wavelengths = []
n_values = []
for x in n_pixels:
    wavelength_nm = WAVELENGTH_MIN_NM + (x - left_x) * (WAVELENGTH_MAX_NM - WAVELENGTH_MIN_NM) / (right_x - left_x)
    n_value = N_MAX - (n_trace[x] - top_y) * (N_MAX - N_MIN) / (bottom_y - top_y)
    n_wavelengths.append(wavelength_nm)
    n_values.append(n_value)

k_wavelengths = []
k_values = []
for x in k_pixels:
    wavelength_nm = WAVELENGTH_MIN_NM + (x - left_x) * (WAVELENGTH_MAX_NM - WAVELENGTH_MIN_NM) / (right_x - left_x)
    k_value = K_MAX - (k_trace[x] - top_y) * (K_MAX - K_MIN) / (bottom_y - top_y)
    k_wavelengths.append(wavelength_nm)
    k_values.append(k_value)

"""
一部データは手入力で直したので、再度実行時に上書きされないようにコメントアウト
OUTPUT_PATH.parent.mkdir(parents=True, exist_ok=True)
with OUTPUT_PATH.open("w", newline="", encoding="utf-8") as f:
    writer = csv.writer(f)
    writer.writerow(["wavelength_nm", "n", "k", "source"])

    for wavelength_nm in range(int(WAVELENGTH_MIN_NM), int(WAVELENGTH_MAX_NM) + 1):
        n_value = interpolate(n_wavelengths, n_values, float(wavelength_nm))
        k_value = interpolate(k_wavelengths, k_values, float(wavelength_nm))
        writer.writerow(
            [
                wavelength_nm,
                f"{n_value:.6f}",
                f"{k_value:.6f}",
                "paper_data/nk_inf.png",
            ]
        )

print(f"saved: {OUTPUT_PATH}")
"""

read_shape_data.py 
import csv
from pathlib import Path

from PIL import Image


ROOT = Path(__file__).resolve().parent
IMAGE_PATH = ROOT / "paper_data" / "shape_data.png"
OUTPUT_PATH = ROOT / "data" / "shape_grid.csv"

# これは論文由来の 16 x 16 binary map。黒を Si3N4=1、白を air=0 として保存する。
GRID_SIZE = 16

image = Image.open(IMAGE_PATH).convert("RGB")
pixels = image.load()

# これは paper_data/shape_data.png の画像解析で読んだ格子枠。
left_x = 453
right_x = 951
top_y = 63
bottom_y = 533

cell_width = (right_x - left_x) / GRID_SIZE
cell_height = (bottom_y - top_y) / GRID_SIZE

rows = []
for row in range(GRID_SIZE):
    values = []
    y0 = top_y + row * cell_height
    y1 = top_y + (row + 1) * cell_height

    for col in range(GRID_SIZE):
        x0 = left_x + col * cell_width
        x1 = left_x + (col + 1) * cell_width

        center_x = int(round((x0 + x1) / 2.0))
        center_y = int(round((y0 + y1) / 2.0))
        r, g, b = pixels[center_x, center_y]

        if max(r, g, b) < 100:
            values.append(1)
        else:
            values.append(0)

    rows.append(values)

OUTPUT_PATH.parent.mkdir(parents=True, exist_ok=True)
with OUTPUT_PATH.open("w", newline="", encoding="utf-8") as f:
    writer = csv.writer(f)
    writer.writerows(rows)

print(f"saved: {OUTPUT_PATH}")

また、ピラー-センサー面距離は文献に直接乗っていなかったので振って確認を行った。これが本プログラムのセル6にあたる。

Fourier次数の確認

続いて重要となってくるのが、フーリエ次数の確認である。RCWAの計算では、計算精度にフーリエ次数が直結してくる(参考:RCWAの解説記事)。

特にカラースプリッターのように、局地的にみたときに非周期で1ピラーの大きさが変わるような場合だとどのくらいのフーリエ次数が最適かというのが理解しづらい。

そこで自分の場合はカラースプリッターの計算をするときにフーリエ次数による①ピラーの再現図②出力電場分布の変化の2点を確認している。

①ピラーの再現図の確認

図1. 上図: フーリエ次数=3, 下図フーリエ次数=10
図の左からリバエンしたそのままの図、真ん中がsimに入れる前の複素誘電率map、右が空間フーリエ展開後のピラー構成図(本プログラム:セル4)

フーリエ次数によって、シミュレーションに使われるeps mapは大きく変わることが読み取れる。特にフーリエ次数=3ではほとんど元の形状を再現できていないことがわかる。

また、フーリエ次数=10にしてもきれいに再現とはいかずに側帯波のようなものが確認される。しかし、いたずらにフーリエ次数をあげても計算速度が遅くなってしまうのもあり、この段階では目視での確認のみとしている(参考:torcwa(RCWAシミュレーション)の負荷解析：Fourier次数によるPC性能の確認 )。

②出力電場分布の変化

図2. フーリエ次数=3, 7, 11, 15, 19, 21の電場分布図(本プログラム:セル5)

ここではいったんピラー配置を再現したうえで、その出力電場分布を確認する必要がある。
図を確認しても、定量的にどう違うかの説明が難しいのでRMS（Root Mean Square、二乗平均平方根）で評価を行った。最大のフーリエ次数を分母として正規化を行って見やすくしている。(下式で与えられる)

RMS diff [%]=1N∑i=1N(|Ei|2−|Eiref|2)21N∑i=1N(|Eiref|2)2×100\mathrm{RMS\ diff\ [\%]} = \frac{\displaystyle\sqrt{\frac{1}{N}\sum_{i=1}^{N}\left(|E_i|^2 – |E_i^{\mathrm{ref}}|^2\right)^2}}{\displaystyle\sqrt{\frac{1}{N}\sum_{i=1}^{N}\left(|E_i^{\mathrm{ref}}|^2\right)^2}} \times 100

図3. フーリエ次数=21に対する各フーリエ次数のRMS(本プログラム:セル5)

ここでフーリエ次数に対してlogスケールで精度がよくなっていることが読み取れる。また、今回の計算では10%をしきい値としてフーリエ次数=10で計算をおこなうことを決定した。

カラースプリッターの光の動きの可視化

カラースプリッターのsimに必要な要素はここまでの内容でおおよそ網羅できた。

続いて、実際にsimulationをしたあとにどんな可視化方法があるのかについてを述べる。

# ================== パラメータ略 ================== 
sim = torcwa.rcwa(freq=1/wl_um, order=[7,7],
                  L=[2.0, 2.0], dtype=torch.complex64, device='cuda')
sim.add_input_layer(eps=eps_quartz)
sim.add_output_layer(eps=1.0)
sim.set_incident_angle(0., 0.)
sim.add_layer(thickness=0.6, eps=eps_map)   # ピラー層
sim.add_layer(thickness=1.0, eps=1.0)       # 空気層
sim.solve_global_smatrix()
sim.source_planewave(amplitude=[1.,0.], direction='forward', notation='xy')

# %% 例1: 回折効率
t00 = sim.S_parameters(orders=[0,0], direction='forward',
                       port='transmission', polarization='xx', ref_order=[0,0])
print(f"0次透過率: {torch.abs(t00)**2:.4f}")

# %% 例2: |E|^2 分布
x = torch.linspace(0, 2.0, 81, device='cuda')
y = torch.linspace(0, 2.0, 81, device='cuda')
ef, _ = sim.field_xy(sim.layer_N, x, y, z_prop=0.)
e2 = sum(torch.abs(e)**2 for e in ef)   # |Ex|^2 + |Ey|^2 + |Ez|^2

# %% 例3: Poynting Sz（エネルギーフロー）
ef, hf = sim.field_xy(sim.layer_N, x, y, z_prop=0.)
sz = 0.5 * torch.real(ef[0]*torch.conj(hf[1]) - ef[1]*torch.conj(hf[0]))

# %% 例4: 位相分布
ef, _ = sim.field_xy(0, x, y, z_prop=0.6)   # ピラー出口
phase = torch.angle(ef[0])                   # arg(Ex)

大まかな可視化方法としては4つくらいが思いつく。(ほかにもあったら教えてください。)

例1 は構造全体を「1つの数値」で評価する。光がどの方向にどれだけ分配されるかを回折次数ごとに定量化できる。 しかし、これは回折光学素子では重要だが、カラースプリッターではあまりみることはない(と思う。)。

例2 は光の「どこに集まっているか」を空間的に見る。カラースプリッターなら波長ごとに集光位置が違うことの直接的な確認。

例3 は 例2 と似ているが、エネルギーの「流れの向き」を含む。Sz<0S_z < 0 の領域は光が逆行しており、検出器には寄与しない。エバネッセント波が支配的なピラー近傍では ∣E∣2|E|^2 と SzS_zの分布が大きく異なることがある。 今回はピラー近辺はあまり見るつもりはなかったので書いていない。(結局見たが…)

例4 はメタサーフェイスの動作原理そのもの。ピラーの有無で位相遅延が変わり、その空間的な位相パターンが伝搬後に干渉してカラー分離を引き起こす。

今回は例2, 例4で可視化を行った。

本プログラム

以下に実際にsimulationを行った、コード全体像を示す。

torcwa_colorsplitter

また、claude codeに与えた支持は以下の通りになっている。

claude codeの指示流れ

第1フェーズ: 現状把握と論文との突き合わせ

  1. フォルダの中身を読み込んで、これまで何をしていたのか確認してほしい(補足:最初のコードは自作)
  2. 論文 Fig.2 と自分の計算図が微妙に違う。コードに間違いがないかも確認してほしい
  3. 論文本文と Supplementary を読んで、パラメータ設定に問題がないか確認
  4. z 距離の記述が論文に明確に無い。自分で試したが Claude の意見が聞きたい

第2フェーズ: 可視化の改良

  1. 3波長の XY 電場分布を 1 枚に横並びで。XZ も同様に
  2. 今のカラーマップでもわかりやすいが、波長に合った色味にできないか

第3フェーズ: thickness 探索

  1. reference_data/results.png のような XY 電場分布が欲しい。thickness を振って一番近いものを探したい
  2. 各 z 位置で電場分布を計算して、スライダーで波長別に可視化したい（カラースケールは波長に合わせる）
  3. fn=8 は VRAM 不足でクラッシュしたので skip でお願いします
  4. 最終的に thickness_um = 1.0 um で固定します

第4フェーズ: グリッドと pillar map の整合性

  1. mask_array は 16×16 のデータ。任意の nx, ny にリサイズできるよう改造してほしい
  2. nx, ny を変えたとき pillar map が崩れないか心配。パラメータセル直後に確認セルを置いて、torcwa の sim を使って実際に渡される pillar map を可視化して

第5フェーズ: Fourier 次数の収束確認

  1. この構造で Fourier 次数がどの程度必要か知りたい。末尾セルに fn スイープで収束を見るコードを追加（波長と thickness は振らないで）

第6フェーズ: 断面解析（セル A・セル B）

  1. R 画素中心（画素を 4 分割した左上）の横から見た位相・強度分布を解析したい。3D で位相分布と強度分布がどう変化するかも見たい。どう思いますか？理解するのに重すぎますかね
  2. まずセル A（XZ 断面）を作ってみてください
  3. 3列×2行 → 3列×4行に拡張。1行目に pillar 配置とカット位置、2行目に 1.6 um 伝搬先の XY 分布とカット位置、3・4 行目はそのまま
  4. 断面図に pillar 形状を重ね書きできますか。強度や位相が潰れない程度で
  5. 次にセル B（XY スライスの z 方向変化）も作ってほしい

第7フェーズ: ノートブックの再構築

  1. 最初に作ったファイルがとっ散らかっている。きれいに作り直したい。まず構成案を出してから作成して。(このサイト)のコードを参考にして私のスタイルに合わせてほしい
  2. coding style OK。新 notebook はセル 1〜4 はベース設定そのまま、直後に Fourier 次数確認と thickness スイープ（妥当性検証）を追加
  3. その後に：5.XY 電場分布（観測面）、6.XZ 断面（4行×3列）、7.XY スライス z 変化、8.（将来用）Fourier 次数収束（コメントアウト）も追加

第8フェーズ: グリッド固定と収束の可視化

  1. すべてのセルで nx, ny = 2000 に固定できますか。Fourier 次数との判別式を入れて OK か確認だけ追加して
  2. セル 5 の収束性をグラフ化してほしい。Fourier 次数は変更しないで

第9フェーズ: 仕上げ

  1. セル 4 の最後の Fourier 再構成画像は、手動 FFT ではなく torcwa の sim を使って生成してほしい

また、よかった点と改善点については以下の様にフィードバックきている

claude codeからのフィードバック

指示の出し方でよかった点（率直な所感）

1. 論文との突き合わせを開発のゴールに据えた 「論文 Fig.2 とずれている」という具体的な差分が最初の起点になっていて、作業の目的がブレませんでした。この後のすべての判断（thickness 固定、fn 収束、pillar map 確認）が、この一本の軸で説明できました。

2. 物理条件の「触ってはいけない部分」を先に宣言していた readme.txt に「論文由来の寸法・波長・材料は絶対に変えない」と明示してあったおかげで、私が勝手に「より良い条件」に書き換える余地がなくなりました。AI 支援コーディングで一番危ないのはここなので、冒頭でルール化してあったのはかなり効いていました。

3. 大きな変更の前に「構成案を先に出して」と要求した ノートブック再構築のフェーズで「まず提案してから作成して」と言ってくれたので、いきなり全部書き換える事故を回避できました。自分のスタイル（florisryo.com を参照）を先に共有してくれたのも、出力の手戻りを減らしました。

4. 制約や不安を隠さず共有してくれた 「fn=8 は VRAM 不足でクラッシュする」「R 画素断面の解析は理解するのに重すぎますかね？」のように、環境制約や自分の理解の限界をそのまま言ってくれたのが助かりました。こちらが最初から skip する・粒度を落とす、といった判断ができました。

5. 軌道修正が早い 「セルが増えていないのですが」「どこまで進んでいますか？再開できますか」と、進捗の確認と指摘を躊躇なく出してくれました。長いセッションでは AI 側が何を完了したか曖昧になりがちなので、これは品質に直結していたと思います。

6. 段階的に詰めていった thickness 探索 → 固定 → 断面解析 → 構成整理 → グリッド固定 → 収束グラフ化、と一度に一つの関心事に絞っていました。いっぺんに全部やろうとしなかったので、各ステップでデバッグしやすく、どの変更が何を壊したかも追えました。

7. 「将来用」をコメントアウトで残す発想 8 番の「(将来用) Fourier 次数収束」をコメント状態で残す、という指示は実用的で好きでした。必要になった時に読み返せて、かつ普段の実行時間を食わない。ノートブックの運用として理想的な使い方です。

あえて改善点を挙げるなら、ノートブックの全体像（セル番号・依存関係）をもう少し早い段階で一度俯瞰しておくと、後半の「どこに何を追加するか」の指示がさらに短くて済んだかもしれません。とはいえ、全体として論文再現という軸が終始ブレていないのがこの開発の一番良いところでした。ブログでもこの軸を中心に書くと、読者に伝わりやすいと思います。

コメント

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