Pythonで膜厚フィッテイング（実践編）

このページは準備編の続きである。準備編では薄膜における多重干渉を考慮した透過スペクトルの式を導出した。このページではこの結果に基づいてpythonのコードを実装する。まずは、準備編のおさらいから。

やりたいこと

以下の図のような状況で膜厚を算出したい。

ここでガラスの厚みは0.1mm<の十分な厚みをもつものとしよう。
（つまり、ガラスによって生じる干渉は無視できるものとする。理由は準備編を参照）
さらにガラスの屈折率と、Thin Filmの屈折率は既知であるとしよう。
さらに、測定系は以下のようなものを想定する。つまり薄膜を透過スペクトルを膜厚フィッティングに用いる。

光源は無偏光白色光源を想定している。
ここで、
$\theta_{in}$は入射媒体から入射側ガラスへの入射角、
$\theta_0$は入射側ガラスから薄膜への入射角、
$\theta_1$は薄膜から出射側ガラスへの入射角、
$\theta_2$は出射側ガラスから出射媒体への入射角、
となっている。また、
$n_{in}$は入射媒体の屈折率、
$n_0$は入射側ガラスの屈折率、
$n_1$は薄膜の屈折率、
$n_2$は出射側ガラスの屈折率、
である。
このような測定系で光源として無偏光白色光源を使用したときに得られるスペクトルは以下のようなものとなる。

さて、このスペクトルを準備編で導出した式を使ってフィッティングしよう。

フィッティング関数

準備編ではs偏光とp偏光で透過スペクトルが違うことを説明した。しかし、一方でスペクトルの変調の極大と極小の位置は両者で違いはないことも説明した。なので、どちらか一方を使ってもフィッティングすることは可能だ。（さらに言えば、測定系の光軸に対して、セルが傾いてなければ、すなわち$\theta_{in}=0$ならば、s偏光とp偏光は同じスペクトルとなるので、上のようなことを気にする必要はまったくない。）
ということで、フィッティング関数はs偏光の透過スペクトル$T_s$のみを用いて以下のような形とする。

f(\lambda)=(a_1\lambda^2+a_2\lambda+a_3)\biggl(1-\bigl(1-T_s(a_4,\lambda)\bigr)(a_5\lambda+a_6)\biggr)

ここで、$a_1$、$a_2$、$a_3$、$a_4$、$a_5$、$a_6$はフィッティング変数であり、このうち
$a_4$が求めたい膜厚に対応している。$T_s$は準備編で説明したように、$a_4$と波長$\lambda$の関数だ。
これを見てフィッティング関数は単純に$T_s$だけでも良さそうなのに、なぜこんなごちゃごちゃした式になっているのかと思われるかもしれない。だが、実際の透過スペクトルが、理想的な式とぴったり一致することはほとんどない。今回の例でいえばこんなことが起こりうる。

• 透過スペクトルの変調の振幅が$T_s$と一致しない。
• 透過スペクトルの変調の振幅が短波長なほど小さくなる。
• 多重干渉による変調の他に、透過率の変化がある

このようなことが起こりうる理由は様々だ。例えば、薄膜が可視域に吸収をもっていたり、あるいは波長分散が大きかったりすると上のようなことが起きる。また、薄膜の膜厚が場所によって違うような場合(例えば勾配があるとか)には$T_s$によって与えられる変調の振幅よりも実際のスペクトルの変調の振幅は小さくなる。これは、光源からの光束のサイズが有限であることを考えると納得がいく。つまり、この時の透過スペクトルは光束が通過する薄膜のそれぞれ膜厚の異なる各点での透過スペクトルを平均したものとして得られるため、純粋に一定の膜厚を想定した$T_s$より鈍った形となるのだろう。
とまあ、理屈をゴニョゴニョ考えてより理想的な式を立ててもよいが、そういうよくわからん要素は適当に$\lambda$の多項式にして、おパソコンに処理してもらおう、という思想でできたのが上の式だ。なのでフィッティング変数$a_1$、$a_2$、$a_3$、$a_5$、$a_6$には物理的な意味はほとんどないが、一応説明をくわえておくと、

• $(a_1\lambda^2+a_2\lambda+a_3)$は多重干渉による変調以外のスペクトルの形の変化を記述する。例えば、薄膜が可視域に吸収特性をもっている場合などを想定している。$f(x)$が厳密に$T_s$と一致するなら、$a_1=a_2=0$、$a_3=100$だ。
• $(a_5\lambda+a_6)$は多重干渉による変調の振幅を記述する。$f(x)$が厳密に$T_s$と一致するなら、$a_5=0$、$a_6=1$だ。

ではこれを使って、単純にフィッティングしてみよう。

Trial 1 : 単純なフィッティング

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit　
import string
import codecs
import glob
import os

n_in = 1 #入射媒体の屈折率
n_out = 1 #出射媒体の屈折率
n_0 = 1.51633 #ガラス基板の屈折率(defaultはBK7の1.51633)
n_1 = 1 #サンドウィッチ媒体の屈折率(defaultは空気の1)
n_2 = 1.51633 #ガラス基板の屈折率(defaultはBK7の1.51633)
angle_in = 0 #入射角
fit_range = [500,700] #Fittingする範囲
param_list = []
d_answer = 5873 #膜厚の正解値[nm]
# fit_rangeの下限と上限に相当するインデックスを波長のデータ系列から取得する。      
def range_limiter(x,lower,upper) :
    i = 0
    j = 0
    dt = x[1] - x[0]
    if dt > 0 :
        while x[i] <= lower :
            i += 1
        while x[j] <= upper :
            j += 1
        return [i,j]
    else :
        while x[-i] <= lower :
            i += 1
        while x[-j] <= upper :
            j += 1
        return [-j,-i]

# 入射角angle_inからスネルの法則に従ってangle_0、angle_1、angle_2を計算する。
def calc_angle(angle_in, n_in, n_out, n_0, n_1, n_2) :
    angle_0 = np.arcsin((n_in / n_0) * np.sin(angle_in))
    angle_1 = np.arcsin((n_0 / n_1) * np.sin(angle_0))
    angle_2 = np.arcsin((n_1 / n_2) * np.sin(angle_1))
    return [angle_0, angle_1, angle_2]

# s偏光の振幅透過率t_sを計算する
def calc_t_for_s_polarized(angle_0, angle_1, angle_2, n_0, n_1, n_2, w, d) :
    r_1 = (n_1*np.cos(angle_1)-n_0*np.cos(angle_0)) / (n_1*np.cos(angle_1)+n_0*np.cos(angle_0))
    t_1 = 2*n_0*np.cos(angle_0) / (n_1*np.cos(angle_1)+n_0*np.cos(angle_0))
    r_2 = (n_2*np.cos(angle_2)-n_1*np.cos(angle_1)) / (n_2*np.cos(angle_2)+n_1*np.cos(angle_1))
    t_2 = 2*n_1*np.cos(angle_1) / (n_2*np.cos(angle_2)+n_1*np.cos(angle_1))
    sigma = (4*np.pi/w)*n_1*d*np.cos(angle_1)
    return (t_1*t_2*np.exp(-1j*sigma/2))/(1+r_1*r_2*np.exp(-1j*sigma))

# p偏光の振幅透過率t_pを計算する
def calc_t_for_p_polarized(angle_0, angle_1, angle_2, n_0, n_1, n_2, w, d) :
    r_1 = (n_0*np.cos(angle_1)-n_1*np.cos(angle_0)) / (n_1*np.cos(angle_0)+n_0*np.cos(angle_1))
    t_1 = 2*n_0*np.cos(angle_0) / (n_1*np.cos(angle_0)+n_0*np.cos(angle_1))
    r_2 = (n_1*np.cos(angle_2)-n_2*np.cos(angle_1)) / (n_2*np.cos(angle_1)+n_1*np.cos(angle_2))
    t_2 = 2*n_1*np.cos(angle_1) / (n_2*np.cos(angle_1)+n_1*np.cos(angle_2))
    sigma = (4*np.pi/w)*n_1*d*np.cos(angle_1)
    return (t_1*t_2*np.exp(-1j*sigma/2))/(1+r_1*r_2*np.exp(-1j*sigma))

# フィッティング関数
def transmittance_fit(x,a,b,c,d,e,f) :
    t_s = calc_t_for_s_polarized(angle[0], angle[1], angle[2], n_0, n_1, n_2, x, d)
    T_s = np.abs(t_s)**2
    return (a*x**2+b*x+c)*(1-(1-T_s)*(e*x+f))

# フィッティング
def simple_fit(d_ini) :
    param_ini = [0,0,100,d_ini,0,1] #初期値
    lim_idx = range_limiter(list_x, fit_range[0], fit_range[1])
    lim_x = list_x[lim_idx[0]:lim_idx[1]] #fitting範囲に基づいて、データ系列をトリミング
    lim_y = list_y[lim_idx[0]:lim_idx[1]] #fitting範囲に基づいて、データ系列をトリミング
    try :
        param, cov =  curve_fit(transmittance_fit, lim_x, lim_y, p0 = param_ini, maxfev = 10000)
    except RuntimeError :
        param = [0,0,0,0,0,0]
        pass
    return param        
       
angle = calc_angle(angle_in, n_in, n_out, n_0, n_1, n_2)  
os.chdir("透過スペクトルのファイルが格納されたディレクトリを入れてね")
file_list = sorted(glob.glob("*.txt")) #ディレクトリの.txtファイルのリストを取得する。
file_num = len(file_list)
for i in range(file_num) :
    file = file_list[i]
    file_name, ext = os.path.splitext(file)
    list = []
    list_x =[]
    list_y = []
    with codecs.open(file, "r", "utf-8", "ignore") as fileobj : #ファイルの読み込み
        for line in fileobj :
            if string.digits.find(line[0]) == -1 and line[0] != "-":
                continue
            pre_list = line.strip().split()
            list_x.append(float(pre_list[0]))
            list_y.append(float(pre_list[1]))
    
    list_x = np.array(list_x)
    list_y = np.array(list_y)
    
    d_ini_list = np.linspace(1000,10000,100) #初期値のリスト
    for j in range(len(d_ini_list)) :
        param = simple_fit(d_ini_list[j])
        param_list.append(param[3])
        print(j)

    plt.rcParams["font.size"] = 18
    fig = plt.figure(figsize=(16,10))
    ax1 = fig.add_subplot(1,1,1)
    ax1.plot(d_ini_list, np.array(param_list), "-", label="fit result")
    ax1.plot(d_ini_list, np.ones(len(d_ini_list))*d_answer, "-", label="answer")
    ax1.legend()
    ax1.set_xlabel("initial parameter [nm]")
    ax1.set_ylabel("calculated distance [nm]")
    ax1.set_ylim([np.min(d_ini_list), np.max(d_ini_list)])

このスクリプトのすることは以下の通りだ。膜厚の初期値のリストd_ini_listから初期値を取得してフィッティングを行う。そして得られた膜厚のフィッティング結果をparam_listに追加している。なんでこんなことをしているのかというと、このフィッティング方法がどれだけ有用か、つまり初期値にどれだけ依らずに正確なフィッティングをしてくれるかを確認したいからだ。ここでフィッティングするデータの膜厚の正解値はd_answerで与えてある。つまり、5873nmだ。これに対し、1000nm～10000nmの間で100個ほど初期値を選んでフィッティングしてやって、正解をだすかを確認している。ここで初期値のリストd_ini_listは以下のように設定している。

d_ini_list = np.linspace(1000,10000,100)

また、膜厚以外のフィッティング変数は$f(x)$が$T_s$と一致すると仮定して設定する。つまり、$a_1=a_2=0$、$a_3=100$、かつ$a_5=0$、$a_6=1$としよう。
さて、結果はどうだろうか…わくわk…

全然だめや！オレンジのラインが正解値で、青が初期値とフィッティング結果の関係を示すが、ほとんど初期値とフィッティング結果が変化していない。つまり、この方法でフィッティングするには、正解値にごくごく近い初期値を入れてやらねばダメで、こんなものは使い物にならない。では…どうするか。一つ方法がある。最初からナイスな初期値を見積もってやればよいのだ。ここで今まで説明していない、多重干渉のスペクトルから膜厚を見積もる別の方法を紹介しよう。
#Trail 2 : 初期値の推定
初期値の推定方法としてはピークバレー法とよばれる方法がある。この節は島津製作所が紹介している膜厚測定法を参考にした。この方法は透過スペクトルの変調の極大と極小のうち、任意の二点間の波長の差から膜厚を算出する方法だ。

このように透過率が極大、極小となる波長$\lambda_1, \lambda_2, \lambda_3, ...$がある時($\lambda_1$は任意の極大または極小から取ってきて構わない)、膜厚$d$は以下の式で記述される。

d=\frac{k-1}{4n_1\cos\theta_1\left(\frac{1}{\lambda_1}-\frac{1}{\lambda_k}\right)}

ここで$k=2,3,4...$である。リンク先の式と多少形は違うが、この理由は以下の通り。まず、リンク先の式における$\theta$はこのページにおける$\theta_0$であり、上式はこれを$\theta_1$で表現する形にしたことと、極大だけでなく極小も利用するために、リンク先の式の右辺を2で割った形になっている。
ここから言えることは、例えば横軸に$k-1$、縦軸に$\frac{1}{\lambda_1}-\frac{1}{\lambda_k}$をとってプロットすると、各点は理想的には一直線に分布するということだ。そしてその時の直線の傾きは、

slope = \frac{1}{4n_1d\cos\theta_1}

となる。実際には各点は多少ばらつくので、最小二乗法で直線を当てはめることで、傾きを得る。
さて、この方法を先程フィッティングで使用した透過スペクトルに適用してみる。すると$d$は5857nmと得られる。正解が5873nmなのでかなり近い。それほど精度が必要ない場合はこれで十分かもしれない。透過スペクトルと膜厚が5857nmの時のフィッティング関数をグラフにすると以下のようになる。

青のデータが透過スペクトルの生データで、橙がフィッティング関数を示す。かなり近い…が全体的に横にずれていて、透過率の極大、極小の絶対位置が一致していない。しかし、ここまでサポートしてやれば、後は多重干渉を考慮したフィッティングで精度を追い込めるだろう。
ここまでのコードをまとめる。ローパスのコードはこちらのページを参考にした。

# dの初期値を推定する関数
def calc_d_from_crossp(x,y):
    dm_list = [] #(k-1)のリスト
    diff_two_inv_lambda_list = [] #(1/λ_1-1/λ_k)のリスト
    dy = np.gradient(y) #yの一次微分
    fdy = low_pass_filter(x,dy) #dyにローパスをかけてノイズ除去
    lim_idx = range_limiter(x,fit_range[0],fit_range[1])
    lim_x = x[lim_idx[0]:lim_idx[1]]
    lim_fdy = fdy[lim_idx[0]:lim_idx[1]]
    zero_p_idx = zero_grad_counter(lim_fdy) #dy=0になるポイント(=極大、極小)のインデックスを取得
    for i in range(1, len(zero_p_idx)) :
        dm = i
        diff_two_inv_lambda = 1/lim_x[zero_p_idx[0]] - 1/lim_x[zero_p_idx[i]]
        dm_list.append(dm)
        diff_two_inv_lambda_list.append(diff_two_inv_lambda)
    dm_list = np.array(dm_list)
    diff_two_inv_lambda_list = np.array(diff_two_inv_lambda_list)
    slp = np.polyfit(dm_list, diff_two_inv_lambda_list, 1)[0] #ax + bでフィッティング
    d = 1/(4*n_1*slp*np.cos(angle[1]))
    return d

# ローパスフィルタ
def low_pass_filter(x,y) :
    dt = x[1] - x[0]                  #サンプリング間隔
    fn = 1 / (2*dt)                   #ナイキスト周波数
    fp = 0.06                         #通過域端周波数[Hz]。適当に設定。
    fs = 0.12                         #阻止域端周波数[Hz]。適当に設定。
    gpass = 1                         #通過域最大損失量[dB]。適当に設定。
    gstop = 40                        #阻止域最小減衰量[dB]。適当に設定。
    Wp = fp/fn                        #正規化
    Ws = fs/fn                        #正規化
    # バターワースフィルタ
    N, Wn = signal.buttord(Wp, Ws, gpass, gstop)
    b1, a1 = signal.butter(N, Wn, "low")
    y1 = signal.filtfilt(b1, a1, y)
    return y1

#極大、極小をとるインデックスを取得
def zero_grad_counter(y) :
    zero_list = []
    for i in range(len(y)-1) :
        if np.sign(y[i])*np.sign(y[i+1]) == -1 :　#傾き0は正負が逆転する点として計算
            zero_list.append(i)
    return zero_list

#Trial 3 : 初期値を推定したフィッティング（ハイブリッド法）
ピークバレー法で推定した膜厚を初期値として、さらに多重干渉を考慮したフィッティングを行った。結果は以下のようになる。

ピークバレー法のみを用いた場合に比べて、透過スペクトルの変調の極大と極小の絶対位置も揃っている。算出された膜厚は5873nmであり、正解と一致していることが確認できた。さらにこのハイブリッド法では、初期値を入力する必要がないため、初期値によってフィッティングの成功の可否が左右することを恐れる必要もない。この方法でマニュアルで設定しないといけないのは以下の２つ。

1. フィッティング範囲
2. ローパスの透過域端周波数と阻止域端周波数

フィッティング範囲は、データの精度が高そう(ノイズが低そう)な領域を選ぶと良い。
ローパスの通過域端周波数と阻止域端周波数はデータに応じて適当なものを入れる必要があるが、阻止域端が変調の周波数以下とかよっぽど変な設定をしない限りうまくいく。これら２つも将来的にはオートにしたいところだが(ノイズレベルが低い領域を選ぶとか、変調の周波数をビークバレー法で推定した膜厚から見積もって適当なローパス設定をしてくれるとか、多分難しくはないのだろうけど。)
すべてまとめたコードは以下の通り。

import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
from scipy.optimize import curve_fit, basinhopping
import string
import codecs
import glob
import os

n_in = 1 #入射媒体の屈折率
n_out = 1 #出射媒体の屈折率
n_0 = 1.51633 #ガラス基板の屈折率(defaultはBK7の1.51633)
n_1 = 1 #サンドウィッチ媒体の屈折率(defaultは空気の1)
n_2 = 1.51633 #ガラス基板の屈折率(defaultはBK7の1.51633)
filter_use = 1 #1で使用、0でpass
fit_method = 0 #1でbasin-hopping、0でnormal-fit
angle_in = 0 #入射角
fit_range = [500,700] #Fittingする範囲
LP_lower = 0.06 #通過域端周波数[Hz]
LP_upper = 0.12 #阻止域端周波数[Hz]


# ローパスフィルター
def low_pass_filter(x,y) :
    n = len(x)
    dt = x[1] - x[0]                  #サンプリング間隔
    fn = 1 / (2*dt)                   #ナイキスト周波数
    fp = LP_lower                     #通過域端周波数[Hz]
    fs = LP_upper                     #阻止域端周波数[Hz]
    gpass = 1                         #通過域最大損失量[dB]
    gstop = 40                        #阻止域最小減衰量[dB]
    Wp = fp/fn                        #正規化
    Ws = fs/fn                        #正規化
    # バターワースフィルタ
    N, Wn = signal.buttord(Wp, Ws, gpass, gstop)
    b1, a1 = signal.butter(N, Wn, "low")
    y1 = signal.filtfilt(b1, a1, y)
    return y1

# 極大、極小のインデックスを取得
def zero_grad_counter(y) :
    zero_list = []
    for i in range(len(y)-1) :
        if np.sign(y[i])*np.sign(y[i+1]) == -1 :
            zero_list.append(i)
    return zero_list

# fit_rangeの下限と上限に相当するインデックスを波長のデータ系列から取得する。         
def range_limiter(x,lower,upper) :
    i = 0
    j = 0
    dt = x[1] - x[0]
    if dt > 0 :
        while x[i] <= lower :
            i += 1
        while x[j] <= upper :
            j += 1
        return [i,j]
    else :
        while x[-i] <= lower :
            i += 1
        while x[-j] <= upper :
            j += 1
        return [-j,-i]

#ピークバレー法で初期値を推定する。
def calc_d_from_crossp(x,y):
    dm_list = []
    diff_two_inv_lambda_list = []
    dy = np.gradient(y)
    if filter_use == 1:
        fdy = low_pass_filter(x,dy)
    else :
        fdy = dy
    lim_idx = range_limiter(x,fit_range[0],fit_range[1])
    lim_x = x[lim_idx[0]:lim_idx[1]]
    lim_fdy = fdy[lim_idx[0]:lim_idx[1]]
    zero_p_idx = zero_grad_counter(lim_fdy)
    for i in range(1, len(zero_p_idx)) :
        dm = i
        diff_two_inv_lambda = 1/lim_x[zero_p_idx[0]] - 1/lim_x[zero_p_idx[i]]
        dm_list.append(dm)
        diff_two_inv_lambda_list.append(diff_two_inv_lambda)
    dm_list = np.array(dm_list)
    diff_two_inv_lambda_list = np.array(diff_two_inv_lambda_list)
    slp = np.polyfit(dm_list, diff_two_inv_lambda_list, 1)[0]
    d = 1/(4*n_1*slp*np.cos(angle[1]))
    return d
    
# 入射角angle_inからスネルの法則に従ってangle_0、angle_1、angle_2を計算する。
def calc_angle(angle_in, n_in, n_out, n_0, n_1, n_2) :
    angle_0 = np.arcsin((n_in / n_0) * np.sin(angle_in))
    angle_1 = np.arcsin((n_0 / n_1) * np.sin(angle_0))
    angle_2 = np.arcsin((n_1 / n_2) * np.sin(angle_1))
    return [angle_0, angle_1, angle_2]

# s偏光の振幅透過率t_sを計算する
def calc_t_for_s_polarized(angle_0, angle_1, angle_2, n_0, n_1, n_2, w, d) :
    r_1 = (n_1*np.cos(angle_1)-n_0*np.cos(angle_0)) / (n_1*np.cos(angle_1)+n_0*np.cos(angle_0))
    t_1 = 2*n_0*np.cos(angle_0) / (n_1*np.cos(angle_1)+n_0*np.cos(angle_0))
    r_2 = (n_2*np.cos(angle_2)-n_1*np.cos(angle_1)) / (n_2*np.cos(angle_2)+n_1*np.cos(angle_1))
    t_2 = 2*n_1*np.cos(angle_1) / (n_2*np.cos(angle_2)+n_1*np.cos(angle_1))
    sigma = (4*np.pi/w)*n_1*d*np.cos(angle_1)
    return (t_1*t_2*np.exp(-1j*sigma/2))/(1+r_1*r_2*np.exp(-1j*sigma))

# p偏光の振幅透過率t_pを計算する
def calc_t_for_p_polarized(angle_0, angle_1, angle_2, n_0, n_1, n_2, w, d) :
    r_1 = (n_0*np.cos(angle_1)-n_1*np.cos(angle_0)) / (n_1*np.cos(angle_0)+n_0*np.cos(angle_1))
    t_1 = 2*n_0*np.cos(angle_0) / (n_1*np.cos(angle_0)+n_0*np.cos(angle_1))
    r_2 = (n_1*np.cos(angle_2)-n_2*np.cos(angle_1)) / (n_2*np.cos(angle_1)+n_1*np.cos(angle_2))
    t_2 = 2*n_1*np.cos(angle_1) / (n_2*np.cos(angle_1)+n_1*np.cos(angle_2))
    sigma = (4*np.pi/w)*n_1*d*np.cos(angle_1)
    return (t_1*t_2*np.exp(-1j*sigma/2))/(1+r_1*r_2*np.exp(-1j*sigma))

# フィッティング関数
def transmittance_fit(x,a,b,c,d,e,f) :
    t_s = calc_t_for_s_polarized(angle[0], angle[1], angle[2], n_0, n_1, n_2, x, d)
    T_s = np.abs(t_s)**2
    return (a*x**2+b*x+c)*(1-(1-T_s)*(e*x+f))

# フィッティング
def normal_fit() :
    param_ini = [0,0,0,d_ave,0,1]
    lim_idx = range_limiter(list_x, fit_range[0], fit_range[1])
    lim_x = list_x[lim_idx[0]:lim_idx[1]]
    lim_y = list_y[lim_idx[0]:lim_idx[1]]
    try :
        param, cov =  curve_fit(transmittance_fit, lim_x, lim_y, p0 = param_ini, maxfev = 10000)
    except RuntimeError :
        pass
    return param

# フィッティング(basin hopping method)
def basin_hopping_fit() :
    param_ini = [0,0,0,d_ave,0,1]
    lim_idx = range_limiter(list_x, fit_range[0], fit_range[1])
    lim_x = list_x[lim_idx[0]:lim_idx[1]]
    lim_y = list_y[lim_idx[0]:lim_idx[1]]    
    def calc_cost(param, x, y) :
        return ((y - transmittance_fit(x, param[0], param[1], param[2], param[3], param[4], param[5]))**2).sum()
    minimizer_kwargs = {"args":(np.array(lim_x), np.array(lim_y))}
    try :
        result = basinhopping(calc_cost, param_ini, stepsize=2.,minimizer_kwargs=minimizer_kwargs)
    except RuntimeWarning or RuntimeError  :
        pass
    
    return result.x
   
       
# Main       
angle = calc_angle(angle_in, n_in, n_out, n_0, n_1, n_2)  
os.chdir("透過スペクトルのファイルが格納されたディレクトリを入れてね")
file_list = sorted(glob.glob("*.txt"))
file_num = len(file_list)
for i in range(file_num) :
    file = file_list[i]
    file_name, ext = os.path.splitext(file)
    list = []
    list_x =[]
    list_y = []
    with codecs.open(file, "r", "utf-8", "ignore") as fileobj :
        for line in fileobj :
            if string.digits.find(line[0]) == -1 and line[0] != "-":
                continue
            pre_list = line.strip().split()
            list_x.append(float(pre_list[0]))
            list_y.append(float(pre_list[1]))
    
    list_x = np.array(list_x)
    list_y = np.array(list_y)
    
    d_ave = calc_d_from_crossp(list_x, list_y)
    if fit_method == 0 :
        param = normal_fit()
    else :
        param = basin_hopping_fit()
    plt.rcParams["font.size"] = 18
    fig = plt.figure(figsize=(16,10))
    ax1 = fig.add_subplot(1, 1, 1)
    ax1.plot(list_x, list_y, "-", label="raw")
    ax1.plot(list_x, transmittance_fit(list_x, param[0], param[1], param[2], param[3], param[4], param[5]), label="fit")
    ax1.set_xlim([400,800])
    ax1.set_ylim([0,100])
    ax1.set_xlabel("Wavelength [nm]")
    ax1.set_ylabel("Transmittance [%]")
    ax1.legend()
    ax1.text(500, 50,f" thickness={param[3]} \n a={param[0]} \n b={param[1]} \n c={param[2]} \n e={param[4]} \n f={param[5]}")
    plt.savefig(f"{file_name}.png")
    plt.close()

おわりに

薄膜を通過した透過スペクトルを用いて薄膜の膜厚を推定するスクリプトを作成した。ピークバレー法による膜厚の見積もりと、それを初期値に用いて多重干渉を考慮したフィッティングを行うハイブリッド法により、初期値の入力を必要とせずに精度良好な膜厚を得ることが可能となった。
一方、フィッティング範囲とローパス設定はマニュアルで行う必要が残されており、いずれこれらのオート化も取り組んでみたい。