Takeoff and incident angle from TauP

TauP 1.2以上的版本可以在taup_time裡面求得該phase的takeoff angle以及incident angle，筆者目前所能下載到的最新版本為2.0 beta6，相信如果需要這些資訊的研究人員可以更輕易的使用taup。TauP的作者Philip非常貼心的把Takeoff angle以及incident angle欄位放在後面，這樣一來使用者就無需更改呼叫taup的shell script。

taup_time -ph PcP,ScP -deg 46.5 -h 122.2

Model: iasp91
Distance   Depth   Phase   Travel    Ray Param  Takeoff  Incident  Purist    Purist
  (deg)     (km)   Name    Time (s)  p (s/deg)   (deg)    (deg)   Distance   Name
-----------------------------------------------------------------------------------
   46.50   122.2   PcP      587.05     3.532     15.12    10.62    46.50   = PcP
   46.50   122.2   ScP      808.93     4.185      9.94    12.61    46.50   = ScP

請各位無痛升級到新版的TauP吧！

Share

Taup on Windows

Java在眾多程式語言裡面脫穎而出，最大的原因除了物件導向的程式設計易於維護之外，他的跨平台特性確實讓其他語言望其項背。Java透過JVM(Java虛擬機器)來執行，因此跨平台的任務就交給撰寫JVM的廠商或團體，我們只要負責寫或是用就好。Taup就是由Java語言所開發而成，理所當然他是能夠在windows上面運行的。如果你對於linux還是有抗拒但又必須用到Taup，那麼請參考這篇文章。

首先，你必須要有可以運行Java程式的環境(Java Runtime Environment，簡稱JRE)。如何知道你有沒有安裝java呢？請打開你的終端機：(哦，不！是命令提示字元....)，輸入java -version：

如果你可以看到跟筆者類似的畫面，而不是找不到這個指令，那麼就表示你的系統已經有安裝JRE(也就是大家所說的Java)，可以順利執行Taup。If not,請google一下java這個字，找到下載java的頁面。

安裝完Java後，請你再敲一次java -version，應該就可以看到跟第一張圖類似的畫面。接下來，我們需要設定環境變數讓你的Windows知道如何找到taup！
假設你下載Taup並把她解壓縮後放到C:\，則你需要設定一個名為TAUP_HOME的變數：

有了這個變數還不夠，你還需要把Taup的可執行目錄加進系統命令搜尋的PATH裡！

這個path環境變數非常重要，請小心利用分號把%TAUP_HOME%\bin加到最後面。按下確定之後還不行，請你重新開機以讓環境變數生效。(我不知道有沒有不重開機就生效的方法，如果有，可以利用留言分享給其他朋友)

上圖...有點浪費版面...
總之重新啟動系統後，TauP就可以正常執行了！

這是Taup的圖形介面，下圖是linux上執行的畫面。

各位可以發現，java程式的優勢就是設計師不必擔心不同平台有不同的顯示效果，唯一不同之處就是字型，其餘完全相同。

可以在windows上執行所有Taup的功能，或許可以對於學習Taup有些許的幫助。

後記：
筆者測試的結果是Windows上的跑Taup的速度比linux或是solaris上慢得太多，或許是因為我的windows是虛擬機器的關係；但如果你的實體windows也跑得不快，也別太意外。

Share

TauP指令筆記

TauP包含了八個指令以及一個圖形介面模式。圖形模式容易看得到圖，但是實際上，TauP強大的威力必須藉由命令列才能夠真的發揮。本文僅僅作為自己以及地震學同好的參考，若有疏漏，請不吝指正。

taup_time
主要目的：輸入深度、距離等參數，求出各種phase的走時以及ray parameter(horizontal slowness)。
可用參數：

• -ph
後面接一長串的phases list，每個phase之間以','分隔，例如PP,ScP
• -pf
後面接內含phase名稱的檔案。
• -mod[el]
接速度模型名稱，預設是iasp91,另外還可以選擇prem, qdt, ak135，以及你自訂的速度構造模型。
• -h
depth(evdp)，單位是km。
• -deg
大圓路徑的距離，單位是度。
• -km
大圓路徑的距離，單位是km。
• -sta
後面接stla及stlo。
• -evt
後面接evla及evlo。
• -rayp
只印出ray parameter。
• -time
只印出travel time。
• -o
指定output檔名。
• -debug
偵錯模式。
• -verbose
詳細輸出模式。
• -version
印出Taup版本號。
• -help
印出詳細參數資訊。

由於-o, -debug, -verbose, -version以及-help每個指令都有;且不在我們處理數據的目的之內，下面的command就不再列出。
使用範例：

在evdp=130km, gcarc=65^o，以iasp91的模型看P,PcP,ScP的走時以及ray parameter:
taup_time -h 130 -deg 65 -ph P,PcP,ScP

Model: iasp91
Distance   Depth   Phase   Travel    Ray Param   Purist    Purist
  (deg)     (km)   Name    Time (s)  p (s/deg)  Distance   Name
----------------------------------------------------------------
   65.00   130.0   P        625.99     6.471      65.00  = P
   65.00   130.0   PcP      657.70     4.134      65.00  = PcP

在evdp=130km, stla/stlo=40/120.5, evla/evlo=12/140，以prem模型看ScP的走時及ray parameter:
taup_time -mod prem -sta 40 120.5 -evt 12 140 -h 210 -ph ScP

Model: prem
Distance   Depth   Phase   Travel    Ray Param   Purist    Purist
  (deg)     (km)   Name    Time (s)  p (s/deg)  Distance   Name
----------------------------------------------------------------
   32.86   210.0   ScP      736.43     3.470      32.86  = ScP

在evdp=330km, gcarc=4350km, iasp91的模型下看PKiKP的資訊，並匯出成PKiKP.dat檔：
taup_time -evdp 330 -km 4350 -ph PKiKP -o PKiKP.dat

不會有任何畫面上的輸出，因為結果已經輸出到PKiKP.dat裡。

taup_pierce
主要目的：輸入深度、距離等參數，求出各種ray穿透到各深度的時間、距離、走時，甚至是投影到地表的經緯度位置。
可用參數：

• -ph
後面接一長串的phases list，每個phase之間以','分隔，例如PP,ScP
• -pf
後面接內含phase名稱的檔案。
• -mod[el]
接速度模型名稱，預設是iasp91,另外還可以選擇prem, qdt，以及你自訂的速度構造模型。
• -h
depth(evdp)，單位是km。
• -deg
大圓路徑的距離，單位是度。
• -km
大圓路徑的距離，單位是km。
• -sta
後面接stla及stlo。
• -evt
後面接evla及evlo。
• -az
需搭配-evt方能使用。
• -baz
需搭配-sta方能使用。
• -turn
只印出反射面上端(topside of reflection)，用於crust/mantle邊界。
• -under
只印出反射面下端(underside of reflection)，用於crust/mantle邊界。
• -rev
同時印出反射面上下端。但是可能程式有一些邏輯問題，-rev跟-turn行為完全相同。
• -pierce
接一個深度值，可以多計算該深度的到時以及ray parameter。
• -nodiscon
搭配-pierce使用，僅印出-pierce所指定深度的資訊。

使用範例：

在evdp=122.5km, gcarc=60^o時，印出ScP travel time隨距離變化的資訊：
taup_pierce -h 122.5 -deg 60 -ph ScP

> ScP at   867.49 seconds at    60.00 degrees for a    122.5 km deep source in the iasp91 model.
    0.00   122.5     0.0
    0.15   210.0    19.7
    0.53   410.0    63.2
    1.12   660.0   111.4
   13.48  2889.0   469.6
   57.47   660.0   785.6
   58.62   410.0   813.9
   59.36   210.0   838.7
   59.65   122.5   850.1
   59.92    35.0   861.6
   59.96    20.0   863.9
   60.00     0.0   867.5

此時的三個欄位分別是gcarc, depth以及travel time。

在stla/stlo=-40/120, evla/evlo=0/140且evdp=105km的情況下，計算ScP的資訊並額外求出在440km處的travel time, ray parameter等資訊：
taup_pierce -h 105 -sta -40 120 -evt 0 140 -pierce 440 -ph ScP

> ScP at   802.18 seconds at    43.96 degrees for a    105.0 km deep source in the iasp91 model.
    0.00   105.0     0.0      0.00    140.00
    0.17   210.0    23.6     -0.15    139.94
    0.52   410.0    67.0     -0.48    139.81
    0.58   440.0    73.0     -0.53    139.78
    1.06   660.0   115.0     -0.98    139.60
   12.28  2889.0   468.3    -11.36    135.30
   41.65   660.0   721.3    -37.98    121.44
   42.58   440.0   745.7    -38.80    120.87
   42.70   410.0   749.1    -38.90    120.80
   43.37   210.0   773.6    -39.49    120.37
   43.69   105.0   787.1    -39.76    120.17
   43.89    35.0   796.3    -39.94    120.05
   43.92    20.0   798.6    -39.97    120.03
   43.96     0.0   802.2    -40.00    120.00

此時的五個欄位分別是：gcarc, depth, travel time, lat, lon。也會多算出440km處的資訊。

在stla/stlo=-40/120, evdp=105km, gcarc=1387km且baz=65^o時，求ScP的資訊：
taup_pierce -h 105 -km 1387 -sta -40 120 -baz 65 -ph ScP

> ScP at   707.47 seconds at    12.47 degrees for a    105.0 km deep source in the iasp91 model.
    0.00   105.0     0.0    -33.90    133.64
    0.06   210.0    23.3    -33.93    133.58
    0.19   410.0    66.0    -34.00    133.45
    0.39   660.0   113.1    -34.11    133.25
    4.28  2889.0   445.5    -36.16    129.21
   11.66   660.0   630.8    -39.65    120.95
   12.03   410.0   656.7    -39.81    120.53
   12.26   210.0   680.0    -39.91    120.25
   12.38   105.0   692.9    -39.96    120.12
   12.45    35.0   701.7    -39.99    120.03
   12.46    20.0   704.0    -39.99    120.02
   12.47     0.0   707.5    -40.00    120.00

靈活給定-evt/-az以及-sta/-baz，幾乎就掌握了taup_pierce這個指令。欲取得各個phase在構造邊界的反彈點位置，例如ScP反彈點投影到地表的座標，利用taup_pierce把2889.0給grep出來，取第四、第五欄位就可以了！

taup_path
主要目的：輸入深度、距離等參數，用來輸出可以用於psxy的格式，以GMT繪製出ray path。
可用參數：

• -ph
後面接一長串的phases list，每個phase之間以','分隔，例如PP,ScP
• -pf
後面接內含phase名稱的檔案。
• -mod[el]
接速度模型名稱，預設是iasp91,另外還可以選擇prem, qdt，以及你自訂的速度構造模型。
• -h
depth(evdp)，單位是km。
• -deg
大圓路徑的距離，單位是度。
• -km
大圓路徑的距離，單位是km。
• -sta
後面接stla及stlo。
• -evt
後面接evla及evlo。
• -gmt
輸出結果直接是GMT script。

使用範例：

在evdp=560km, gcarc=74^o時，印出PcP,PKiKP的ray path，需匯出成GMT script：
taup_path -h 560 -deg 74 -ph PcP,PKiKP -gmt -o 1999.009.raypath.gmt

不指定輸出檔名時，預設輸出檔名為taup_path.gmt。執行taup_path.gmt會得到這樣的圖型：

taup_curve
主要目的：輸入深度、距離等參數，用來輸出可以用於psxy的格式，以GMT繪製出travel time curve。
可用參數：

• -ph
後面接一長串的phases list，每個phase之間以','分隔，例如PP,ScP
• -pf
後面接內含phase名稱的檔案。
• -mod[el]
接速度模型名稱，預設是iasp91,另外還可以選擇prem, qdt，以及你自訂的速度構造模型。
• -h
depth(evdp)，單位是km。
• -gmt
輸出結果直接是GMT script。
• -reddeg
輸出曲線其速度衰減單位為deg/sec。
• -redkm
輸出曲線其速度衰減單位為km/sec。

使用範例：

taup_curve -h 138 -mod prem -gmt -o curve-deg.gmt

畫出來還不是普通的醜，所以必須要自己改一下才行。

taup_setsac
主要目的：幫助標記phase名稱以及該phase的理論到時，並寫入到sac檔的headers(t0~t9以及kt0~kt9)內。
可用參數：

• -ph
後面接一長串的phases-number，每個phase之間以','分隔，例如PP-1,ScP-5
• -pf
後面接內含phase-number的檔案。
• -mod[el]
接速度模型名稱，預設是iasp91,另外還可以選擇prem, qdt，以及你自訂的速度構造模型。
• -devpkm
設置evdp單位為km(預設是meter)。這個參數很容易忘記，一定要加上去。

使用範例：

將sac檔的PcP與ScP的到時標記起來並寫入到headers裡。
taup_setsac -ph ScP-1,PcP-2 -evdpkm ST40.01.BHZ.sac

ST40.01.BHZ.sac
Swap bytes for linux: swapBytes(npts) (40743*4) + header(632) =  file length=163604

此時看一下sac檔的headers：

SAC> lh t1 kt1 t2 kt2


  FILE: ST40.01.BHZ.sac - 1
 --------------------------------------

      t1 = 8.076501e+02
     kt1 = ScP
      t2 = 5.859711e+02
     kt2 = PcP

關於這個指令，需特別注意sac headers裡的evdp跟o都必需要被設定過(也就是必須有值)；有些phase在某距離內會有多個到時，注意taup只會取第一個到時，其餘忽略。另外，這裡如果沒有加上-evdpkm，則用meter計算會天差地遠；後面的sac檔名不適用萬用字元或正規表示法，如果要三個方向的sac檔都要寫入，必須輸入三個檔名才行。

taup_table
主要目的：將速度構造模型匯出成LocSAT可以使用的travel time table格式，或是一般格式。
可用參數：

• -ph
後面接一長串的phases list，每個phase之間以','分隔，例如PP,ScP
• -pf
後面接內含phase名稱的檔案。
• -mod[el]
接速度模型名稱，預設是iasp91,另外還可以選擇prem, qdt，以及你自訂的速度構造模型。
• -header
接一個LocSAT格式的檔案，內含有depth-distance的資料檔。
• -locsat
輸出LocSAT格式的travel time table(ascii)。
• -generic
輸出一般的travel time table(ascii)。

使用範例：

將prem模型的ScP,PcP的travel time以一般格式輸出：
taup_table -model prem -ph ScP,PcP -generic -o generic_ScP_PcP.dat

輸出的結果部份如下：

prem     0.00      0.0 PcP    510.35     0.000       0.00  PcP
prem     0.00      0.0 ScP    723.01     0.000       0.00  ScP
prem     0.10      0.0 PcP    510.35     0.010       0.10  PcP
prem     0.10      0.0 ScP    723.01     0.012       0.10  ScP
prem     0.20      0.0 PcP    510.35     0.019       0.20  PcP
prem     0.20      0.0 ScP    723.01     0.025       0.20  ScP
prem     0.30      0.0 PcP    510.36     0.029       0.30  PcP
prem     0.30      0.0 ScP    723.02     0.037       0.30  ScP
prem     0.40      0.0 PcP    510.36     0.038       0.40  PcP
prem     0.40      0.0 ScP    723.02     0.050       0.40  ScP

將prem模型的ScP,PcP的travel time以LocSAT格式輸出：
taup_table -model prem -ph ScP,PcP -locsat -o locsat_ScP_PcP.dat

輸出的結果部份如下：

n # ScP,PcP travel-time tables for prem structure. (From TauP_Table)
40     # number of depth samples
   0.00   1.00   2.00   3.00   4.00   5.00  10.00  15.00  19.00  21.00
  25.00  30.00  33.00  35.00  40.00  45.00  49.00  51.00  55.00  60.00
  70.00  80.00  90.00 100.00 120.00 140.00 160.00 180.00 200.00 220.00
 240.00 260.00 280.00 300.00 350.00 400.00 450.00 500.00 550.00 600.00
207    # number of distances
   0.00   0.10   0.20   0.30   0.40   0.50   0.60   0.70   0.80   0.90
   1.00   1.10   1.20   1.30   1.40   1.50   1.60   1.70   1.80   1.90
   2.00   2.10   2.20   2.30   2.40   2.50   2.60   2.70   2.80   2.90
   3.00   3.10   3.20   3.30   3.40   3.50   3.60   3.70   3.80   3.90
   4.00   4.10   4.20   4.30   4.40   4.50   4.60   4.70   4.80   4.90
   5.00   5.10   5.20   5.30   5.40   5.50   5.60   5.70   5.80   5.90
   6.00   6.10   6.20   6.30   6.40   6.50   6.60   6.70   6.80   6.90
   7.00   7.10   7.20   7.30   7.40   7.50   7.60   7.70   7.80   7.90

如果你想使用LocSAT，網站上提供的ftp下載處連結是錯誤的，請下載Linux版本，或是HP unix版本，或是Solaris版本。

taup_create
主要目的：自訂一個速度構造的模型。
可用參數：只有兩個參數可用

• -nd
"named discontinuities"格式的檔案，請參考$TAUP_HOME/modelFiles/prem.nd。
• -tvel
".tvel"格式的檔案，請參考$TAUP_HOME/modelFiles/iasp91.tvel。

使用範例：

修改prem速度模型後，產生一個新的model：
taup_create -nd new_prem.nd

結果會在工作目錄下產生名為new_prem.taup的一個速度模型檔，這是一個binary。

修改iasp91速度模型後，產生一個新的model：
taup_create -tvel new_iasp91.tvel

結果會在工作目錄下產生名為new_iasp91.taup的一個速度模型檔，這是一個binary。
很可惜的，taup_create並沒有-o參數，他必定會用輸入檔當成輸出的檔名。這個新的速度模型檔放在工作目錄下就可以被找到並且指定：

taup_time -mod new_prem -evdp 330 -km 4350 -ph PKiKP -o PKiKP.dat

如此一來便可以搭配自訂模型來處理問題。

taup_peek
主要目的：用於debugging先前taup_create所產生的模型細節問題。
可用參數：只有一個參數可用

• -mod
接我們所產生的model檔名。

使用範例：

taup_peek -mod new_prem

然後在互動式的問答輸入depth跟distance來作詳細的偵錯。

如果你看了這篇筆記，對你有些幫助或是有需要補充的地方，請留言給筆者，謝謝！

Share

TauP-- the travel time calculator

TauP是美國南加州大學(USC)的Lithospheric Seismology Program所開發出來的程式之一，這套程式是用來計算地震波的走時(Travel Time)。除了travel time之外，她還可以計算ray path以及地震波的穿入(pierce)及轉折點(turning)。相較於傳統的ttimes程式，她可以處理更多的速度構造(iasp91, PREM以及QDT)；並且可以分析許多虛擬波相的走時。更難能可貴的是，他是用Java所寫成的，因此不需再編譯，只要安裝JRE即可快樂使用。更詳盡的說明，可以看他的手冊，或是解壓縮後就附有使用手冊了。

TauP也採用GPL授權，他的授權聲明在這裡，如果你有趣的話，原始碼解壓縮後在src裡就有了。

目前最新版本的TauP是1.2，您可以點這裡下載，或是直接到這裡查看是否有最新的版本。我猜測作者們應該是在sparc架構上開發TauP，因為她提供了for little-endian類的處理器的版本，而little-endian恰好是x86架構處理器的排序處理方式；如果您在處理資料的時候出現NaN(Not a Number)的錯誤訊息，很可能是因為headers對不上的關係，建議拿這個版本來用；不過這裡僅提供原始碼，您如果有興趣，需要自己編譯；筆者初次嘗試的結果，Makefile並沒有寫好，自己嘗試編譯的結果也沒有成功；等有空再仔細看看好了。

安裝TauP on Linux/MacOS/Unix

假設我們想把TauP安裝在/opt(當然你也可以依照linux的習慣放在/usr/local底下)：

wget -c http://www.seis.sc.edu/downloads/TauP/TauP-2.0.tgz
sudo tar zxvf TauP-2.0.tgz -C /opt

接下來設定環境變數：新增兩行到/etc/profile

export TAUP_HOME=/opt/TauP-1.2
export PATH=$PATH:$TAUP_HOME/bin

完成之後重新登入，輸入taup看看有沒有結果。

如果看不到圖形跑出來，表示你沒有安裝Java runtime environment。通知偉大的系統管理員安裝他吧！

地科系學生常遇到的 Q&A

1. 為什麼系統檔案裡有許多地方有.profile這個檔案？該修改哪個才好呢？
系統裡一共會有三個地方會有profile這個名稱：/etc/profile, /etc/skel/.profile以及~/.profile(如果你不知道~是什麼意思，麻煩自己google)

/etc/profile
系統啟動時所需要的所有環境變數，所有使用者登入時第一個讀取的檔案，包含root也是。

/etc/skel/.profile
系統新增使用者時，會複製一份這個檔案到使用者的家目錄。這個檔案與系統管理員比較有關。

~/.profile
個人登入時讀取的第二個檔案。如果變數名稱與/etc/profile相同，以後讀的檔案為主。

所以如果是系統管理員，而所有的使用者都是需要用到這個程式，請寫在/etc/profile；如果你的能力勝過你的系統管理員，是自己安裝來自己用的只是沒root權限，請修改~/.profile。
2. 為什麼我的系統裡沒有~/.profile也沒有/etc/profile這個檔案？
因為你的linux是Redhat系列(包含CentOS, ScientificLinux, Fedora, OracleLinux)，請改用/etc/bash_profile以及~/.bash_profile
3. 為什麼我修改完之後，一定需要重新登出？有沒有不重新登出的辦法？
修改完檔案之後，系統並不會自己去讀取環境變數，因此重新登入會強迫使用者重新讀取/etc/profile以及~/.profile這兩個檔案，也就是取得我們新賦予的環境變數。
然而，手動讀取這些環境變數檔也是可以的：

source /etc/profile

或是

source ~/.profile

你會發現這樣一樣有效。這是個『有限』的有效，因為你的桌面環境是前一次登入所讀取的，你桌面的終端機的變數是後來再source而來的，因此你在桌面再開啟的終端機仍然是繼承先前桌面的環境變數，仍然是舊的變數設定。因此強烈建議不熟悉的使用者，重新登入。
4. 我懷疑我的電腦沒安裝java，該怎麼確定？
1. 問你的系統管理員。
2. 假設你的系統管理員不知道，可以自救：

java -version

如果出現command not found，表示沒有安裝。如果出現類似

java version "1.6.0_24"
Java(TM) SE Runtime Environment (build 1.6.0_24-b07)
Java HotSpot(TM) Server VM (build 19.1-b02, mixed mode)

表示你的系統有安裝Java Runtime Environment。
如果出現command not found的話，

aptitude install sun-java6-jre

在Debian/Ubuntu下安裝，或

yum install java-1.6.0-openjdk

在CentOS/Fedora下安裝，或

zypper install java-1_6_0-openjdk

在OpenSuSE下安裝。
5. TauP有沒有文字介面能快速計算地震學的一些參數呢？
當然有，詳細指令請參考TauP指令筆記。TauP在1.2之後新增一個taup_console的指令，這是在taup裡新增一個shell來撰寫Jython語言的一個環境，與計算較無關連。你應該閱讀TauP說明檔來彌補本文的不足。

Share