Havriliak-Nagami式によるフィッティング

ご案内

　誘電緩和のデータ解析で広く使われているHavriliak-Nagami式の重ね合わせによるフィッティングツールを、IgorProのプロシージャとして実現した。プロシージャファイルは以下の３つに分割した。

• Complex.proc
• FitbyManual.proc
• EpsFit.proc

IgorProでフィッティングを行うには、これら３つのファイルをあらかじめ読み込んでおく必要がある。常にこれを使うなら、IgorProフォルダのIgor Proceduresフォルダの中にファイルを３つとも入れておくと便利である。IgorProで他の解析も行うのなら、毎回この３つのファイルをダブルクリックしてIgorProに読み込ませる必要がある。

　このページは、IgorProユーザでIgorProのマニュアルをそれなりにながめたことがあって、プロシージャの書き方を知っているか、またはマニュアルを片手にこれからプロシージャを作ろうという人を対象にしている。書いてある内容がわからない時は、マニュアルを見ること。

Complex.proc

　複素数の計算をする関数の定義。IgorProのversion2.*の時に作ったので、Igor側がサポートしたため不要になったものもあるが、これまでに作ったプロシージャの互換性のために残してある。IgorProの古いバージョンでは、複素数を使うことはできたが複素数を引数にとれる関数が少なかったので、主な関数を自作した。以下の関数が含まれている。

Function/D/C cadd(z1, z2)
Function/D/C csub(z1, z2)
Function/D/C cmul(z1, z2)
Function/D/C cdiv(z1, z2)
Function/D/C cln(z)
Function/D/C clog(z)
Function/D/C cexp(z)
Function/D/C cpow(z1, z2)
Function/D/C csin(z)
Function/D/C ccos(z)
Function/D/C ctan(z)
Function/D/C csqrt(z)
Function/D/C csinh(z)
Function/D/C ccosh(z)

FitbyManual.proc

　これをコンパイルすると、IgorProのメニューバーにTDR Macroという項目が追加されて、

• Set Up Manual Fitting
• Do Fitting Automatically
• Save Current Fitting Coefficients
• Restore Current Fitting Coefficients

の４つの項目のメニューが表示されるようになる。フィッティングの手順を次に示す。

1. fittingしたい複素誘電率のスペクトルをIgorProに読み込む。周波数は対数表示になっていること。
2. Set Up Manual Fittingを選ぶ。x軸、実部、虚部のwaveがどれかきいてくるので、それぞれ（対数表示の）周波数、誘電率の実部、虚部のwaveを指定し、continueをクリックする。
3. 実部と虚部のグラフと、Havriliak-Nagami式５つのパラメータが表示されたcontrolが作られる。パラメータ入力ウィンドウは２つ出て、１つがε∞とHN式３つ、もう１つはHN式２つのパラメータが表示される。（最初HN式３つでやっていたら、５つにしろと言われたのであとから２つ追加したからこうなっている。windowを２つに分けたのは、ディスプレイが小さいとcontrolを表示しているウィンドウがうまく全部表示されないからである。普通は５つも重ね合わせることはないし、重ね合わせたとしてもよほど広い周波数範囲にわたるデータでない限りフィッテイングが収束しなくなるだけである。）
4. それぞれのパラメータに適当な値を入れて、データと計算値を合わせる。HN式５つは明らかに多すぎるので、使う分だけ値を入力する。使わない成分は、緩和強度Δεを０にして、Δε、α、β、τのholdflagをすべてON（チェックした状態）にしておく。このholdをちゃんとしておかないと、マルケート法による自動フィッティングのときにエラーになることがある。
5. データと計算値がほぼ合ったら、Do Fitting Automaticallyを選ぶ。マルケート法によるフィッティングが始まる。このとき、新たにパラメータを表示するtableができる。マルケート法による計算をしている間は、controlに表示されたパラメータは変更されないので、各パラメータがどのように変化しながら収束するか見るには、tableの方を見る。
6. マルケート法の計算が終わると、controlウィンドウのupdateボタンをクリックする。収束したときのパラメータの値がcontrolの値にコピーされる。また、残差２乗和のSigmaは数値の表示部分をマウスでクリックすると残差２乗和の値がコピーされる。
7. 結果が気に入らなければ、再度パラメータを変更して、必要ならばパラメータをholdし、自動フィットを行う。

　それなりにいいパラメータを得たのでとりあえず記録しておいて次のフィッティング結果と比べたいときは、Save Current Fitting Coefficientsを実行すると、パラメータをメモしておける。次のフィッティングをしてみて、それがよろしくないときはRestore Current Fitting Coefficientsを選ぶと記録しておいたフィッティングパラメータに戻せるので、それを初期値にして処理を再開できる。

このプロシージャファイルに含まれている関数やマクロは次の通り。

Function HN5FitManual(w, p)

パラメータを手入力したときに呼ばれるHN式の計算式。

Macro SetupFitting(xaxis, ywave_real, ywave_imag)

必要なwaveを作り、虚部と実部のグラフと、パラメータ設定用のcontrolを出す。

Macro SaveCurrentCoefs()

Save Current Fitting Coefficientsで呼ばれる

Macro RestoreCurrentCoefs()

Restore Current Fitting Coefficientsで呼ばれる

Function MakeControls()

ε∞とHN式３つのパラメータを入力するcontrolを作る

Function MakeControls_2()

HN式２つのパラメータを入力するcontrolを作る

Function SetVarProc_coefs(ctrlName,varNum,varStr,varName) : SetVariableControl

controlに表示されたパラメータの数値に値を設定する。

Function CheckProc_hold(ctrlName,checked) : CheckBoxControl

controlに表示されたholdflagの値をチェックして、FucFitに渡すhold(string)を作る。

Function SetVarProc_show(ctrlName,varNum,varStr,varName) : SetVariableControl

何もしていない。

Function ButtonProc(ctrlName) : ButtonControl

updateボタンをクリックしたときに呼ばれる。

EpsFit.proc

　Do Fitting Automaticallyで行われるマルケート法によるフィッティングの計算ルーチン。

Function/D/C HavN(omega, dele, ltau, alpha, beta)

HN式による誘電率を計算する。周波数、緩和強度、緩和時間の対数、α、βを入れると複素誘電率の値を計算して複素数で返す。

Function/D/C Fit5Relax(w, log_f)

HavN５つを足した値を返す。各パラメータはwに入れる。この形式で関数を定義すると、FuncFitに渡すことができる。/DD>

Function HN5Fit(w, p)

マルケート法で計算するときの関数の評価式の定義。この関数が返す値を最小にするようにフィッテイングが行われる。

	sigma= sqrt((dr + di)/(n-m))
	Return sigma

だと、実部と虚部について残差を計算し、和が最小にするという条件を指定したことになる。

	if(dr>di)
	return sr*sqrt(dr+di)
	else
	return si*sqrt(dr+di)
	endif

だと、実部と虚部について残差を計算し、大きい方を返すので、実部と虚部の残差を同じにするという条件を指定したのと同等である。

Macro DoFit_5Components()

適切なパラメータでFuncFitを呼び出す。

Y.Amo /