4-4 Matlab 程式語言

4-4 Matlab 程式語言

在4-3 節中，我們欲求出在給定頻率下使式（4.26）中之 ∆T^*(ω, n^*) 值有最小值時之n^*值，仔細觀察式（4.26），可發現 ∆T^*(ω, n^*)值恆為 正，當∆T^*(ω, n^*)有最小值時，即是 ∆T^*(ω, n^*)→0 之時，故我們以牛 頓勘根法求解，輸入一測試起始值ntest*，由此起始值尋求方程式（4.26）

之根。

實驗進行中所撰寫的Matlab 程式語言主架構如圖 4-4 所示，主要 由參數設定、傅立葉轉換、數值計算迴圈及矩陣儲存四大結構所組 成，我們以此架構共建立兩個略有差異的程式，分為命名為Program1 及Program2，以下詳述：

Program1

1. 參數設定：

設n^*有一系列之測試值，令此集合nTest = {ntest*|ntest*

1、ntest* 2、 ntest*

3、ntest*

4……}，將 ntest*分別代入式（4.11）及（4.22）中，並以

ntest*做為計算式（4.26）之牛頓勘根法之起始值，若有一 ntest*代入使

76 Matlab 內建函數“fminsearch”以牛頓勘根法在式（4.26）中的 ∆T^*(ω, n^*)函數上尋找以 ntest*為初始值的最小值，其所尋找到之數值解根令

所求出之n^*值便較不可信，反之亦然。此指令之示意圖如圖 4-5 所 示。

整個數計算迴圈的邏輯架構圖如圖 4-6 所示，總共分成兩大部 份，虛線框外部主要為參數的輸入順序法則，虛線框內部則為主要 的邏輯判斷區。先從虛線框外看起，此部份迴圈設定主要為完成前 述之 nTest 集合及 Freq 序列依序丟入式（4.26）中進行計算；虛線 框內部則把虛線框外部的指令依序以fminsearch 指令進行最小值的 尋找，此所尋找到之最小值為 ∆p^*，其所對應的數值解根則令為 nsol*，得到這兩個值後，我們將 ∆p^*與(∆p^*)min做比較，若∆p^*較(∆p^*)min

為小則令新的(∆p^*)min 值為此次所尋找到 ∆p^*，而此 nsol*則令為

nfinal*，若∆p^*較(∆p^*)min為大，則不進行取代，即此次計算結果不予

以考慮，於是，依此規則，我們最後能找出在每一個Freq 對應序列 下所尋找到之(∆p^*)min值及nfinal*值，且此(∆p^*)min值必為此Freq 對應 序列下之最小值，而n^*即為我們所欲求之複數折射率之解。

一開始設定一遠大於∆p^*值之(∆p^*)min值，一般∆p^*值小於10，

且多小於0.1，在此處我們設計(∆p^*)min之初始值為10000。

4. 矩陣儲存：

最後可得在Freq 對應序列下所尋找到的 ntest*、n^*及(∆p^*)min三個 值，將此三個值的實、虛部各對應Freq 序列儲存成六個 2^*n 矩陣，

78

並以.dat 檔輸出。

Program2

Program2 與 Program1 的程式架構相似，主要不同處在於數值計 算迴圈在給定頻率下輸入起始值 ntest*的方式，在 Program1 中輸入一 ntest*之矩陣nTest，而在 program2 中以前一個頻率下所得到之 nfinal*為

現欲求解頻率之ntest*，以避免 Program1 中因多重解而造成之 nfinal*漂 移問題，Program2 數值計算迴圈的邏輯架構圖如圖 4-7 所示。

Program1 具在不同頻率下求得之 nfinal*各自獨立之優點，若訊號 準確度較低時，不致於因在某頻率下求得之 nfinal*失真而影響在其他 頻率下求得之 nfinal*值，但卻無法有效地避開多重解的問題；而 Program2 則具有可避開多重解和速度快的優點，但若起始之 ntest*值設 定不佳，便極有可能引起整個頻率序列下所求得之nfinal*值產生錯誤，

並且在所求之 nfinal*值未知的情形下，極難判定所求得之 nfinal*值是否 正確。前述概略性地介紹兩個程式的優缺，在實驗進行之時，我們通 常以兩個程式反覆交互檢視結果，使兩個程式各以其長補另者之短。

（a）時域訊號分佈

（b）頻域訊號分佈

圖4-1：取樣原理

Time

In ten sity Am pl it ude

Frequency

可信取樣頻率f_max 取樣頻率f_s

80

（a）THz 輻射穿透基板示意圖

（b）THz 輻射穿透薄膜/基板示意圖 圖4-2：THz 輻射穿透樣品示意圖 THz

t₁₂

t21

r₁₂

r21

r21

E_subs(t)

介質1 介質2 介質1

基板 空氣

空氣

THz

t₃₂

t₂₁ r₃₂

r21

r₂₃

E_film(t)

介質1 介質2 介質1

基板 空氣

空氣

t₁₃

介質3 薄膜 r13

r₃₁

d2

d2

d₃

（a）不考慮多重反射訊號

（b）考慮含有一個多重反射訊號

圖4-3：MgO 基板對空氣理論透射係數比 Tthe*與實驗透射係數比 Texp*

0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

0.5 0.6 0.7 0.8 0.9 1.0

Frequency (THz)

Amplitude

0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

0.5 0.6 0.7 0.8 0.9 1.0

Frequency (THz)

Amplitude

－ 理論透射係數比

˙ 實驗透射係數比

－ 理論透射係數比

˙ 實驗透射係數比

82

參數設定：

nStart nStop nStepR nStepI Time Freq 實驗訊號輸入：

E_sam(t).dat E_ref(t).dat

傅立葉轉換：

Eref(t).dat → Eref*(ω) Esam(t).dat → Esam*(ω)

數值迴圈計算：

由fminsearch 指令由牛頓勘根法尋找 n^*

矩陣結果儲存：

ntest,r*.dat ntest,I*.dat nr*.dat ni*.dat、

(Δp^*)min,r.dat (Δp^*)min,i.dat

圖 4-4：程式流程架構圖

圖4-5：fminsearch 指令示意圖 ΔT^*(ω, n^*)

n^* ntest*

1

n_test^*₂

ntest* 3

(n_sol^*₁,Δp^*₁)

(n_sol^*₂,Δp^*₂) (n_sol^*₃,Δp^*₃)

84

圖 4-6：數值計算迴圈邏輯架構圖 以fminsearch 指令尋找

∣T_the^* (ω, n^*)－T_exp^* (ω, n^*)∣最小值，

令此尋找到之最小值為Δp^*， 而使其有最小值之n_test^*值令為 n_final^*

while n_test^*實部 < nstop實部

n_test^* = n_test^* + nStepR while ntest*虛部 < nstop虛部

ntest* = nStart實部 + ntest*虛部 n_test^* = n_test^* + nStepI

for Freq = 1：FreqPoints

成立

if Δp^* < (Δp^*)_min

(Δp^*)min = Δp^* n^* = n_final^*

設定參數：nStepR、nStepI、

n_test^* = nStart、(Δp^*)_min = 10000

成立

不成立

不成立

成立

不成立

圖 4-7：Program2 數值計算迴圈邏輯架構圖 以fminsearch 指令尋找

∣T_the^* (ω, n^*)－T_exp^* (ω, n^*)∣最小值，

令此尋找到之最小值為Δp^*， 而使其有最小值之n_test^*值令為 n_final^*

n_test^* = n_final^*

for Freq = 1：FreqPoints

if Δp^* < (Δp^*)_min

(Δp^*)min = Δp^* n^* = n_final^*

設定參數：

n_test^* = nStart、(Δp^*)_min = 10000

86

第五章