早安，苦命工程師的胡言亂語

前言:

最近開始比較有參加台灣符式語言愛好者協會 FIG Taiwan 的活動，認識了不少喜歡這個很古老語言的愛好者。就上星期，協會的理事長透過推特，傳來了一段在 GitHub 上外國同好用 eFORTH 所撰寫的太陽能熱水器的 FORTH 程式碼。

有大概喵了一下，裡面有線性內插法的程式碼。線性數值內插的程式碼對 FORTH 儀器控制而言，是非常普遍跟重要的，因為不管是量測溫度的 Thermal Couple 感測, 量測真空的真空計，... 及各式各樣的感測器，他們的特性曲線很少是線性的啦。然後你說要有一個完美的理論公式來描述它們通常也是非常困難的。所以通常是建立一個對應表，然後透過查表做線性內插的方式將所需要的數值轉換出來。

當然內插法也很多種，但線性內插法是最方便跟快速的，所以最常被使用。

因為處處都需要，之前我也常撰寫這樣的 FORTH 程式碼，所以這裡也來野人獻曝一下，為這個古老語言來一篇部落格文章分享經驗惹。最後的 Thermal Wave Sensitivity Calculator 程式，取材於 15 年前在 Palm PDA 我所寫的一個計算的 App 。這個 App 有在 Palm Gear 的社群發表過。現在，只是把 Palm GUI 的那層拿掉，抽出核心的計算指令來當範例囉。

 

線性內插

這個其實非常直覺的，就給定兩點 P1(x1,y1), P2(x2,y2)  在這兩點所形成的直線上，如果給定一點x, 然後它的 y 會是多少？

 

很簡單的數學公式  (y2 - y1)/(x2 - x1)  會是這條直線的斜率增率。所以立刻可以寫下

y = y1 + (x-x1) * (y2-y1)/(x2-x1)

這就是內插法的公式，轉成 FORTH 程式碼如下！

這裡要稍微再提一下的，在 FORTH 的儀器控制領域，我們 FORTH 程式設計師先很少人會使用浮點數及浮點運算的。大家都是直接使用固定整數運算的 (fixed point integer)。所以，大家可以順便看一下跟學習一下 FORTH 如何只使用固定整數運算來達到跟使用浮點運算一樣的結果。

這就像開手排車一樣，一開始好像覺得很難。但久了，習慣了，會覺得用純整數來做一些科學運算，也蠻好用的。

 

線性內插程式碼，這邊為了方便，定義兩個雙整數變數 Low 跟 High

Low 就儲存 P1(x1,y1) 兩個單整數值， High 就儲存 P2(x2,y2) 兩個單整數值

2variable Low
2variable High

: L2!  ( x1 y1 --)
   Low 2!  ;

: H2!  ( x2 y2 --)
   High 2!  ;

 

Slope 就計算兩點間斜率的指令，會將 Low P1(x1,y1) 跟 High P2(x2,y2) 兩點間的斜率算出來，以 dy / dx 分數的形式放到堆疊上。

很簡單呀，就 High 2@ 取出 x2, y2 ， y2 先推到返回堆疊。 Low 2@ 取出 x1, y1 ， y1 也給它推到返回堆疊。 參數堆疊兩個數字相減就是 x2-x1 ， 返回堆疊兩的數字推回來交換後相減就是 y2-y1 。  最後整理一下，數字對調就是 dy dx 的斜率惹。

: slope ( --- dy dx)
   High 2@            ( x2 y2)
   >r   Low 2@  >r    ( x2  x1  | y1 y2)
   -
   r>  r>  swap -     ( x2-x1 y2-y1)
   swap
;

 

最後，最重要的計算內插指令 inter 如下，給一個 x ，它會根據 Low, High 兩點的值來執行內插。應用上只要利用 Low 2! 跟 High 2! 存入已知兩點的數值，然後利用 inter 指令，就可以任意內插這兩點間對應任何 x 所需要的 y 值出來，使用上非常的方便。

前幾篇我的 BLOG 文章有稍微提過， FORTH 整數運算的關鍵就是 */ 這個比例運算指令。我們利用它來依照比例來縮放一個單整數。用法為 n1 n2 n3 */ 後會得到 n4 = n1 * n2/n3 ，也就是說 將 n1 以 n2 / n3 這個分數的比例進行縮放，最後得到 n4 於堆疊上。 且 n1 * n2 的過程中，中間結果會是用雙整數來暫放，所以不用擔心最終結果會因為乘法後溢位而爆掉。

Inter 這個指令很簡單，利用 Low 2@ 把 x1, y1 取出來。 y1 先丟到返回堆疊暫存，把 x - x1 算出來，再利用 slope 指令算出 dy dx 斜率。三個數字排好好，然後利用 */ 算出 (x-x1)*dy/dx 。 最後將 y1 從返回堆疊推回來，加上剛剛的結果就完成了內插法的計算:  y = (x-x1)*dy/dx + y1

: inter ( x -- y)
   Low 2@ >r   ( x x1 | y1)
   -   
   slope       ( x-x1 dy dx | y1)
   */ 
   r> +
;

: inter ( x x1 y1 x2 y2 -- y)
   2over rot  -  >r              ( x x1 y1 x2 x1 |  y1-y2)
   swap -        >r              ( x x1 y1 | x1-x2 y1-y2)
   -rot -   r>   r>  swap        ( y1 x-x1 y1-y2 x1-x2)
   */    +
;

 

看，是不是更簡單跟精簡啊！ 

這邊我的習慣，註解裡用 | 來分隔參數堆疊跟返回堆疊。 所以像 ( x x1 y1 x2 x1 | y1-y2) 這樣的註解就代表參數堆疊目前有五個數字，分別是 x, x1, y1, x2, x1 ，而返回堆疊有一個數字為 y1-y2

照著註解所顯示的堆疊變化來看，並不難閱讀的。這段程式碼利用 2over rot -rot swap 及返回堆疊的推放 >r r> 來操控這五個數字。 最後重點就是要完成 y1, x-x1, y1-y2, x1-x2 這四個數值的計算並在堆疊裡依序排好。 後續的 */ 比例運算會算出  (x-x1)*(y1-y2)/(x1-x2)  ，最後的 + 會把 y1 加上去。 所以就完成了 y = y1 + (x-x1)*(y1-y2)/(x1-x2) 這個內插法的數值計算。

 

片段近似法

將 P1(x1,y1), P2(x2,y2) 一直延伸下去，變成很多點 P1(x1,y1), P2(x2,x2), P3(x3,y3) ... Pn(xn, yn) 這樣就叫做片段近似法 (piecewise fitting) 。就用一堆直線來逼近所要的曲線。

這樣給定一個 x ，要計算 y 的話。就是先查表，看這個 x 是落在表格裡那兩點之間，最後根據這兩點的資料作線性內插，就是所要的 y 惹。

 

程式解說

整個片段近似法就是先建立個每個線段 Pi(xi, yi) 的資料表，要內插時就是查詢這個資料表，找到對應的資料線段 P1(x1,y1), P2(x2,y2) ，找到之後再用 inter 指令執行內插，計算內插後的值。

所以主角會是內插表格的查表 LookUp 指令，及表格的結構。

表格的資料結構很簡單，就第一個cell 存表格裡面所有座標點(x,y)的組數。第二個cell 存第一組的 y，第三個 cell 存第一組的 x ，第四個 cell 存第二組的 y，第五個 cell 存第二組的 x ， ... 其他的依此類推。

所以

variable CurrentTable        \ point to the table currently using

: piecewise ( n --- )
   create  
          ,           ( element counts of this table)  
   does>  
      CurrentTable !  ( point to the table currently using)
;

piecewise 用來建造這樣的資料結構的開頭，然後用 FORTH 內建的 2, 指令來編入每個 (x,y) 的資料點。 CurrentTable 是個指標變數，所來指向目前現用的表格位置在哪裡。

語法像這樣，例如有 10 個元素的表格

10 piecewise myTable

  11 22 2,

333 444 2,

...

555 777 2,

 

這裡比較手動一點啦，自己要去數一下元素的數目然後放上來

piecewise 一開始會用 create 建個 myTable詞頭， 逗點會把你給的元素的數目 10 編入這個詞資料區的第一個位置。接下來每個 2, 會將每個堆疊上的兩個數字 x y 以 y x 的順序依序的編入字典。要取出資料，只要利用 2@ 就可以取出來惹。

最後 does> 定義一下這個 Table 被執行時要幹嘛啊？就簡單啊，把目前這個 Table 的位址存入 CurrentTable 以供後續使用。這樣會非常方便的，當你有很多 Table 的時候，例如 myTable1 ， mayTable2 ， myTable3 ... 

例如只要簡單執行一下 myTable3 ，就會自動選擇 MyTable3 將它的位址存入 CurrentTable 這個指標變數。

要使用 myTable2，只要再簡單的執行一下 myTable2 ，就會自動切到 myTable2 去使用惹！

 

解說 LookUp 前，先定義一下 查表過程假如產生錯誤時，會回傳的 error code,

我們定義當 error code > 0 代表有 error 產生

0 代表沒 error順利找到對應的 P1(x1,y1), P2(x2,y2) 。 1 代表 x 值落在表格之外，比最低的值還低。而 2 代表 x 值落在表格之外，比表格所定義的最大值還大。

0 constant NO_ERROR
1 constant UNDER_LOW_LIMIT
2 constant OVER_HIGH_LIMIT

: lookup ( x addr -- addr' 0, or error-code)  
   dup @ 0
   do                         ( x addr )
      [ 2 cells ] literal +   ( x addr')      
      over over @ <           ( x addr'  x<@addr') 
      if   i 0=   ( first element?)
           if     2drop       UNDER_LOW_LIMIT 
           else   swap drop   NO_ERROR   then
           unloop  exit
      then
   loop                         
   2drop  OVER_HIGH_LIMIT
;

 

 

離子佈值機和Thermal Wave 量測

這裡大概提一下，在半導體製程裡面，要形成 n-type 或 p-type 的半導體，是透過參雜 doping 的這個製程來完成的。將5價的雜質(P, As, ..)原子參雜進去4價的矽晶圓中就會產生多餘的自由電子流動形成n-type半導體。 將 3價的雜質(B, Al, Ga, In)原子參雜進去4價的矽晶圓中就會多餘的電洞流動形成 p-type 半導體。參雜製程一直是整個半導體製程的主角之一。

這參雜製程中最有名的莫過於離子佈值機惹。就把我們需要參雜的雜質氣體透過直流電漿解離成電漿態，透過負電壓電極的吸引將這些被解離的雜質正離子給拉出來並加速形成具有特定能量的離子束，再透過內建大型質譜儀的純化，再經過二次加速提升到更高的加速能量及內部離子束掃描系統及計量系統的協同合作，最後就可以很精準的將所需要的特定離子以很精確的能量跟數目均勻地打到整面晶圓的特定深度。這樣就完成了很有名的離子佈值製程。

打進去的雜質原子，當然不可能立刻變成 n-type 或 p-type 半導體。要形成共價鍵結才會有電子電洞的載子形成啦。所以離子佈值製程完畢後，晶圓還要進入高溫爐管或快速回火爐進行高溫回火的製程，高溫能夠提供晶格間的電子跟原子有足夠的動能，讓這些雜質原子可以透過晶格原子間共價鍵結的形成，穩定形成對應的載子跟我們所需的半導體，也同時修復離子佈值過程中因為離子撞擊所造成的晶格缺陷損傷。

那，怎麼監控離子佈值這個製程呢?? 特別是中電流離子佈值機 (佈值的劑量約在 1E10 - 1E13 ions/cm^2 之間)，因為佈值的劑量非常的淡，所以高溫回火完裡面的電子電洞載子少，所以會有極高的面電阻，很容易被回火製程干擾而造成難以判斷是否為離子佈值製程的問題還是回火製程的問題。所以一般來說，要監控低劑量的離子佈值製程都會採用一種很特別的光學監控儀器，稱作 Thermal Wave 量測。

Thermal Wave 量測的儀器是一種利用光學來做量測的非接觸性量測，大概的原理是這樣的。儀器裏頭有兩個雷射源，第一個雷射源是功率比較強的 Ar 雷射，第二個雷射源是功率弱的 He-Ne 雷射。一開始 Ar 雷射會以脈波的方式打入晶圓表面，吸收雷射的能量後會在離子佈值那層造成聲子的晶格熱振動。這時候第二道 He-Ne 雷射會透過雷射照射晶格跟反射至光學感測器的方式，將這些微小的振動訊號給讀取取回來。廠商根據這些訊號會算出一個他們定義的測量值，稱之為 Thermal Wave Unit，或簡寫為 TWU

TWU 可以看作是晶格破壞程度度量的一種量測，離子佈值製程中，當佈值能量越強，佈值計量越多，佈值離子電流越大，佈值離子的原子質量越重，破壞的程度就越強， 量測到的 TWU 數值就愈大。 TWU量測所呈現的晶格破壞程度可以反應離子佈值製程完畢的第一手製程結果，因為量測非常快速，所以被用在這個整個製程的 SPC 監控中。

Thermal Wave 量測原理示意圖

Thermal Wave 根據雷射訊號計算 TWU 的示意圖

 

 

Thermal Wave 靈敏度及SPC 監控

利用 Thermal Wave 來監控離子佈值的製程結果，畢盡這還是一種間接監控的方式。整個製程的結果最重要的還是實際所打進去的離子劑量 (稱之為 Dose 量，單位為 ions/cm^2) 是否均勻跟準確。TWU 畢竟不是我們所真正關心的量啊! 真正的製程結果還是要回歸到 Dose 量。 所以這裡有個重要的量，稱之為靈敏度 Thermal Wave Sensitivity.

它的定義是 Sensitiivity = TWU 變化量 / Dose 變化量

所以就是當單位 Dose 變化時，所能造成的 Thermal Wave TWU 變化的一個比值時。亦即當Thermal Wave 越靈敏，一單位的 Dose 變化時會造成很大的 TWU 變化。這時候 Sensitivity 的數值越大。

為什麼會有這個 Sensitivity 數值呢？？？ 原來，TWU 跟 Dose 量的比例關係並不是恆定的，在不同的離子佈值製程的劑量，能量，離子束電流大小，離子束的寬度，離子的種類... 這兩個比例關係是會改變的。這個比例關係就是 Sensitivity。

製程工程師在設計離子佈值機製程 SPC 監控程式時，往往設計好一個離子佈值監控製程程式，就要實際以不同的離子劑量來試打模擬預估一下這個離子佈值監控製程程式實際可能造成的 TWU 值變化如何，是否足夠靈敏到足以預防當機台實際 Dose 已經改變了，能夠被偵測出來，預先停機，以避免進一步的製程問題而造成大量的產品損失。

所以，假如能夠事先預知在各種情況下的 TWU Sensitivity 的話，會對製程工程師有很大的幫助。

 

美商瓦里安 Varian

在離子佈值的領域中，其中總部在美國波士頓的美國廠商瓦里安 Varian 算是個中的翹楚，因為百分之 90% 以上的商用離子佈值機都是他們所製造銷售的！ (當然，他們現在已經跟有名的應用材料 Applied Material 合併了!)

他們的製程部門曾經在 IEEE 上發表了一篇跟 Thermal Wave Sensitivity 有關的論文。這篇論文說，雖然 Thermal Wave 受很多剛剛提到的參數所影響，但根據他們製程實驗室對 Thermal Wave 十年以上的 SPC 監控資料，他們最後歸納出一個大概的 TWU Sensitiivity vs TWU 的經驗曲線。 (下圖)

 

他們發現， TWU 的 Sensitivity 數值大概會跟 TWU 的絕對值有一個如上圖的經驗曲線關係。也就是雖然 Thermal Wave Sensitivity 受很多參數影響，但是只要能知道 TWU 的數值，就可以透過這個經驗曲線大概預估 Sensitivity，最後可以把大概可能對應的 Dose Error 數值估計出來。

 

Thermal Wave 靈敏度 FORTH 程式碼

這裡就要用我們前面發展的片段近似法程式碼，根據 Varian 所公開的 Thermal Wave Sensitivity 曲線，來撰寫一個 Thermal Wave Sensitivity 計算的程式碼惹，來幫助我們在給定一個 TWU 數值下，能快速的算出 Sensitivity 的數值，以及對應的 Dose Error

例如，我們只要輸入 500 TWU Sensitivity 它就立刻回報此處的 Sensitivity 為 21.77%

然後再輸入 18 TWU DoseError 它就立刻回報，這樣的 Sensitivity 下， +/- 18 TWU 的變化會造成 16.53% 的 DoseError

 

這段程式如下,

先用 piecewise 建個 5個元素的內插表 DoseSens ，用來呈現 Varian 的經驗公式。

這裡 x 值是 TWU， 一位小數點， 所以 28500 實際代表是 2850.0  TWU

這裡 y 值是 Sensitivity ，四位小數點，所以 5000 實際代表是 0.5

再舉例，表裡的 (6250, 1000) 實際代表是 (TWU = 625.0 , Sensitivity = 0.1)

 

5 piecewise DoseSens

     0  6500 2,
  2000  5000 2,
  6500  1000 2,
 22200  7000 2,
 28500 10000 2,

 

定義兩個變數，做結果的暫存，給後面計算 DoseError 時使用

variable Sens
variable TWVal

 

Sensitivity 的計算其實超簡單的，就呼叫我們前面的內插法指令做內插啊，結果就出來惹！ 

: Sensitivity ( TWU --- Sens)   DoseSens pinter ;

 

立刻將結果存入變數，提供後續 DoseError 計算使用，順便印出結果惹！

: ThermalWave ( TWU ---)
   10 * dup TWVal !  Sensitivity  dup Sens ! ."  --> Sensitivity = " .% cr
;

 

Dose Error 的計算，就利用剛剛的結果，再加減乘除一下就是結果惹！  XD

因為 TWU 的改變是 +/- 都算在內的，所以是所給的參數的兩倍。  TWU% 將 TWU的改變乘上2，並適當的調整你看不到，隱形小數點的位數。 XDD

Dose% 會立刻將 Dose 的誤差算出來，然後因為是 +/- 都算在內的，所以最終結果要除2 (就是要乘上0.5) 並適當的調整你看不到，隱形小數點的位數。 XDDD

最後就是兩位小數的結果，然後列印完收工去！ 

 

: TWU% ( DTWU TWU -- TWU%)  20000 swap */ ;

: Dose% ( TWU% Sensitivity -- Dose%) 5000 swap */ ;

: DoseError ( DTWU ---)
   10 * TWVal @ TWU%  Sens @  Dose%  ."  -->  Dose Error = " .% cr
;

 

最後的最後，列印兩位小數的指令

: .% ( data --)  0  <# 37 hold # #  46 hold  #S    #> type ;

 

 

 

xxx

 

最後，原始程式碼列表

 

\
\ Piecewise Fitting
\      2019.07.21   Frank Lin
\

\
\  two points P1(x1,y1), P2(x2,y2) form a line
\  interpolate a new point (x,y) on this line between these two points.
\
\  given x,  got y   by  y = y1 + (x-x1)*(y2-y1)/(x2-x1)
\

 

\
\ Table LookUp
\

 

0 constant NO_ERROR
1 constant UNDER_LOW_LIMIT
2 constant OVER_HIGH_LIMIT

: lookup ( x addr -- addr' 0, or error-code)   \ when error, no result. only error code
   dup @ 0
   do                         ( x addr )
      [ 2 cells ] literal +   ( x addr')           \ step to next element
      over over @ <           ( x addr'  x<@addr') \ found?
      if   i 0=   ( first element?)
           if     2drop       UNDER_LOW_LIMIT      \ error, UNDER_LOW_LIMIT
           else   swap drop   NO_ERROR   then
           unloop  exit
      then
   loop                                            \ error, OVER_HIGH_LIMIT
   2drop  OVER_HIGH_LIMIT
;

 

\
\ Piecewise Table Creation
\

variable CurrentTable   \ point to the table currently using

: piecewise ( n --- )
   create  
          ,                ( element counts of this table)  
   does>  
      CurrentTable !       ( select and point to the table currently using)
;

 

false   ( true = compile method #1, false = method #2)
[if]

\
\ method #1 , more readable version
\

\
\ linear interpolate
\

2variable Low
2variable High

: L2!  ( x1 y1 --)
   Low 2!  ;

: H2!  ( x2 y2 --)
   High 2!  ;

: slope ( --- dy dx)
   High 2@            ( x2 y2)
   >r   Low 2@  >r    ( x2  x1  | y1 y2)
   -
   r>  r>  swap -     ( x2-x1 y2-y1)
   swap
;

: inter ( x -- y)
   Low 2@ >r   ( x x1 | y1)
   -   
   slope       ( x-x1 dy dx | y1)
   */ 
   r> +
;

 

\
\ piecewise fitting
\

: pinter ( x --- y)      \ piecewise inter
   dup   CurrentTable @    
   lookup                         ( x addr' 0, or x err)
   abort" Out Of Range..."
   [ 1 cells ] literal -  dup 2@ H2!
   [ 2 cells ] literal -      2@ L2!
   inter
;

 

[else]

 

\
\ method #2 , less readable but more compact by stack operations
\

: inter ( x x1 y1 x2 y2 -- y)
   2over rot  -  >r              ( x x1 y1 x2 x1 |  y1-y2)
   swap -        >r              ( x x1 y1 | x1-x2 y1-y2)
   -rot -   r>   r>  swap        ( y1 x-x1 y1-y2 x1-x2)
   */    +
;

 

\
\ piecewise fitting
\

: pinter ( x --- y)      \ piecewise inter
   dup   CurrentTable @    
   lookup                           ( x addr' 0, or x err)
   abort" Out Of Range..."
   
       [ 3 cells ] literal - 
       dup                     2@   ( x low.addr x1 y1)
   rot [ 2 cells ] literal +   2@   ( x x1 y1 x2 y2)
   
   inter
;

[then]

 

\
\  Example code
\
\     Thermal Wave Sensitivity Calculator
\

\ piecewise fitting table for thermal wave sensitivity

5 piecewise DoseSens

     0  6500 2,
  2000  5000 2,
  6500  1000 2,
 22200  7000 2,
 28500 10000 2,
 

: .% ( data --)  0  <# 37 hold # #  46 hold  #S    #> type ;

variable Sens
variable TWVal

: TWU ;  ( Dummy word)

: Sensitivity ( TWU --- Sens)   DoseSens pinter ;

 

: TWU% ( DTWU TWU -- TWU%)  20000 swap */ ;

: Dose% ( TWU% Sensitivity -- Dose%) 5000 swap */ ;

 

: ThermalWave ( TWU ---)
   10 * dup TWVal !  Sensitivity  dup Sens ! ."  --> Sensitivity = " .% cr
;

: DoseError ( DTWU ---)
   10 * TWVal @ TWU%  Sens @  Dose%  ."  -->  Dose Error = " .% cr
;

\
\ Example
\

DoseSens

cr
500 TWU ThermalWave 

\ Sensitivity = 21.77% at 500 TWU

18 TWU DoseError

\ Dose Error = 16.53%  at 500 TWU, +/-18 TWU change will cause 16.53% dose error