simulate_hypo_3pointsコマンドマニュアル

(The documentation of simulate_hypo_3points command)

Last Update: 2025/11/5

◆機能・用途(Purpose)

3点法震源決定をシミュレートする。
Simulate hypocenter determination by 3-point method.

●3点法の概要 (Outline of 3-point method)

3点法震源決定は3観測点でのP波・S波到着時刻の読み取り値を用いて図式的に震源決定を行う古典的な手法であり、以下の手順から成る。
The 3-point method is a classical approach to locate an earthquake graphically using picked P- and S-wave arrival times at three stations. This method is composed of the following steps.

観測点$i=1,2,3$におけるP波到着時刻の読み取り値を$t_i^p$、 S波到着時刻の読み取り値を$t_i^s$として、その差(S-P時間) \[\begin{equation} \Delta t_i\equiv t_i^s-t_i^p \label{eq.timediff.SP} \end{equation}\] を計算する。
Let $t_i^p$ and $t_i^s$ be the picked P- and S-wave arrival times at station $i=1,2,3$, respectively. Compute its difference (S-P time; Eq. \ref{eq.timediff.SP}).
S-P時間に係数$k$を掛けて震源距離 \[\begin{equation} R_i\equiv k\Delta t_i \label{eq.R} \end{equation}\] に換算する。ここで$k$は媒質のP波・S波速度によって決まる定数である。 3点法では均質媒質を仮定する。そのP波速度を$V_p$、S波速度を$V_s$として \[\begin{equation} k=\frac{V_pV_s}{V_p-V_s} \label{eq.k} \end{equation}\] である。
Multiply a factor $k$ with the S-P time to derive a hypocentral distance $R_i$ (Eq. \ref{eq.R}). The factor $k$ is determined from P- and S-wave velocities of the medium. A homogeneous medium is assumed in the 3-point method, and $k$ is determined by Eq. (\ref{eq.k}), where $V_p$ and $V_s$ are the P- and S-wave velocities, respectively.
地図上で各観測点$i$を中心とする半径$R_i$の円を描画する(図1赤・緑・青線)。 2つの円の共通弦(3通りの円の組み合わせがあるので計3本)を引くとそれらは1点で交差し、その地点が震央となる(図1黒線)。
Draw a circle of radius $R_i$ centered on each station $i$ on a map (red, green, and blue lines in Fig. 1). Draw the chords shared by two circles; because there are three combinations of two circles, a total of three chords are drawn (black line in Fig. 1). They cross at a single point, which is the epicenter.
共通弦の1つを直径とする半円を描画し、震央からその半円に垂線を引く(図1紫線)。この垂線の長さが震源の深さを表す。
Draw a half circle whose bottom line is one of the chords drawn in the previous step. Draw a line parpendicular to the chords from the epicenter to the circle (purple line in Fig. 1). Its length represents the depth of the hypocenter.

図1. 3点法の手順。後述の例題において$k=9.0$ km/sとした場合の結果を示す。
Fig. 1. Procedures of the 3-point method, from the example below with $k=9.0$ km/s.

●このプログラムの目的 (Purpose of this program)

このプログラムは3点法を教える側の立場で教材を準備する際の補助として用いることを想定している。 3点法の教材の準備に当たっては適切な観測点や$k$の選択が必要であり、そのためには様々な観測点組合せや$k$の値で3点法を繰り返し試す必要がある。その際にいちいち手で描画していたのでは効率が悪い。そこでこのプログラムでは

3観測点の座標
それらの観測点でのP波、S波到着時刻の読み取り値
仮定する$k$の値

を与えると3点法を模した計算をコンピュータ内で行い、

震源位置
描画する円や線の座標データ

をテキストファイルで出力する。これを例えば Generic Mapping Tools 等のコマンドベースの描画ツールと組合せてループ処理を行うことで観測点組合せや$k$の値を様々に変えた3点法のシミュレーションを短時間で行うことができる。
This program is assumed to be used to support preparing a lecture of the 3-point method as a teaching side. In preparing the lecture, the 3-point method has to be repeated for various trial combinations of stations and $k$ values to find out appropriate stations and $k$. Doing this work manually is inefficient. Therefore, this program simulates the 3-point method in a computer; it reads the following data:

the coordinates of three stations,
picked P- and S-wave arrival times at these stations, and
assumed $k$ value,

and outputs the following data into text files:

the hypocenter location, and
the coordinates of lines and circles to be drawn.

Combining this program with command-based plotting tools (e.g., Generic Mapping Tools) in a loop operation enables a rapid simulation of 3-point method with various combinations of stations and $k$ values.

◆ソースコード(Source code)

$YMAEDA_OPENTOOL_DIR/hypo/src/simulate_hypo_3points.c

◆使用方法(Usage)

コマンドライン引数でパラメータを指定する。パラメータの一覧を下表に示す。
Specify parameters by command-line arguments. The table below shows a list of parameters.

●「-」から始まらない引数 (Arguments not beginning with “-”)

引数
Argument

与える値
Quantity to be given

第1引数
1st argument

入力データファイル名(テキストファイルであり拡張子は自由)。
The name of an input data file (a text file with an arbitrary extension).

このファイルには使用する3観測点の座標とP波・S波到着時刻の読み取り値を以下の書式で与える。列の区切りにはタブを使用する。空行と各行の「#」から後の部分はコメントとして無視されるので自由に挿入できる。
In this file, write the coordinates of three stations and the picked P- and S-wave arrival times at these stations with the following format. Use tabs to separate the columns. Empty lines and parts after # in each line can be inserted freely because they are ignored as comments.

列 Column	値 Value
1	観測点コード。 A station code.
2	観測点の緯度(度単位または度:分:秒表記)。 The latitude of the station (in degrees or degree:minute:second format).
3	観測点の経度(度単位または度:分:秒表記)。 The longitude of the station (in degrees or degree:minute:second format).
4	観測点の標高(m)。現在のところ観測点の標高は考慮しないアルゴリズムであるが将来的な拡張性のため標高欄を設けている。 The elevation (m) of the station. This field is preserved for future extension of the program, although the current version of the algorithm does not take into account the elevation of the station.
5	P波到着時刻の読み取り値(s)。 S波とP波の到着時刻差のみが解析に用いられるので時刻原点は自由。 A picked P-wave arrival time (s). Arbitrary time origin can be used because only the time difference between S- and P-wave arrival times is used for the analysis.
6	S波到着時刻の読み取り値(s)。 A picked S-wave arrival time (s).

第2引数
2nd argument

出力ディレクトリ名。
The output directory name.

●1つの「-」から始まる引数 (Arguments beginning with a single “-”)

このコマンドでは1つの「-」から始まる引数は存在しない。
This command does not have arguments beginning with a single “-”.

●「--パラメータ名=パラメータ値」の形式の引数 (Arguments of a form “--Parameter name=Parameter Value”)

「--パラメータ名=パラメータ値」の形式の引数は自由な順番で指定できる。「-」から始まらない引数の間に挿入しても良い。相反する指定がなされた場合には後の指定が優先される。デフォルト値を持つパラメータは省略できる。
Arguments of a form “--Parameter name=Parameter Value” can be placed in an arbitrary order. They can even be inserted between arguments not beginning with “-”. In case of conflicting options being specified, the latter option has a higher priority. Parameters that have default values can be omitted.

パラメータ名 Parameter name	意味 Meaning	可能なパラメータ値 Allowed parameter values	デフォルト値 Default value
k	仮定する$k$の値(km/s)。 The assumed $k$ value (km/s).	正の実数。 A positive real number.	8.0
refN	直交座標系への変換の基準点の緯度(度単位または度:分:秒表記)。 The latitude (in degrees or degree:minute:second format) of the reference point for a conversion to a cartesian coordinate system.	-90.0以上+90.0以下の実数またはその度:分:秒表記。 A real number greater than or equal to -90.0 and less than or equal to +90.0, or its degree:minute:second representation.	3観測点の緯度の平均値。 The average of the latitudes of the three stations.
refE	直交座標系への変換の基準点の経度(度単位または度:分:秒表記)。 The longitude (in degrees or degree:minute:second format) of the reference point for a conversion to a cartesian coordinate system.	-180.0以上+180.0以下の実数またはその度:分:秒表記。 A real number greater than or equal to -180.0 and less than or equal to +180.0, or its degree:minute:second representation.	3観測点の経度の平均値。 The average of the longitudes of the three stations.

◆動作(Behaviour)

第2引数で指定したディレクトリを作成し、その下に以下のファイルを出力する。
Create a directory specified by the 2nd argument and output the files listed below.

●hypo.dat

求めた震源の座標。第1行を緯度(度)、第2行を経度(度)、第3行を深さ(km)とし、それぞれ「latitude=値」「longitude=値」「depth=値」の書式で出力される。各観測点を中心とする円が交差せず震源が求まらないケース(図2)では「hypocenter not determined」と出力される。これは仮定した$k$の値が実際の速度構造から期待される値に比べて小さすぎる場合に起きる。
The coordinate of the hypocenter. The 1st to 3rd lines indicate the latitude (deg), longitude (deg), and depth (km), which are expressed in the formats “latitude=value”, “longitude=value”, and “depth=value”. If the hypocenter is not determined because of no intersection of circles centered on the stations (Fig. 2), the file consists of the texts “hypocenter not determined”; this situation occurs when the assumed $k$ is substantially smaller than the value expected from the actual velocity structure.

図2. 震源が求まらないケース。後述の例題において$k=7.0$ km/sとした場合の結果を示す。
Fig. 2. A case where the hypocenter is not determined, from the example below with $k=7.0$ km/s.

●circlei_STATION.dat

$i$番目の観測点 (観測点コード:STATION) を中心とする半径$R_i$の円の座標データ。直交座標系において観測点の真東から反時計回りに1°刻みで円周上の点の座標を計算し、それを緯度経度に戻した値を並べたものである。 1行につき円周上の1点の座標(第1列:経度、第2列:緯度、単位:度)を表す。
The coordinates of a circle of radius $R_i$ centered on $i$th station (station code: STATION). The coordinates of the circumference on a cartesian system are computed every 1° counterclockwise from the east of the station and converted to the longitudes (in the 1st column) and latitudes (in the 2nd column).

●circle_depth.dat

震源の深さを求めるための半円の座標データ。震源が求まる場合にのみ出力される。第1引数で指定した入力データファイルにおける 1行目の観測点を中心とする円と2行目の観測点を中心とする円の共通弦を直径とし、 1行目の観測点の側を向いた半円が選択される。直交座標系において反時計回りに1°刻みで円周上の点の座標を計算し、それを緯度経度に戻した値が並ぶ。 1行につき円周上の1点の座標(第1列:経度、第2列:緯度、単位:度)を表す。
The coordinates of a half circle to determine the depth of the hypocenter. This file is created only if the hypocenter is determined. The bottom line of the half circle is the chord shared by the circles centered on the 1st and 2nd stations in the input data file (1st argument), and the circular part of the half circle is oriented toward the 1st station. The coordinates of the circumference on a cartesian system are computed every 1° counterclockwise and converted to the longitudes (in the 1st column) and latitudes (in the 2nd column).

●lines.dat

観測点を中心とする円の共通弦や震源の深さを求めるための垂線の両端の座標データ。震源が求まる場合にのみ出力される。第1列は経度、第2列は緯度(単位:度)であり、以下の順番で出力される。異なる線の区切りには「>」とだけ書かれた行が出力される。
The coordinates of the end points of the chords shared by two circles and the line perpendicular to one of the chords used to determine the depth of the hypocenter. This file is created only if the hypocenter is determined. The 1st and 2nd columns of this file are the longitude and latitude (in degrees), respectively, and the order of data is as follows. Lines that consist of a string “>” are used to separate different line objects.

第1引数で指定した入力データファイルにおける 1行目の観測点を中心とする円と2行目の観測点を中心とする円の共通弦
The chord shared by the circles centered on the 1st and 2nd stations in the input data file (1st argument)
第1引数で指定した入力データファイルにおける 1行目の観測点を中心とする円と3行目の観測点を中心とする円の共通弦
The chord shared by the circles centered on the 1st and 3rd stations in the input data file (1st argument)
第1引数で指定した入力データファイルにおける 2行目の観測点を中心とする円と3行目の観測点を中心とする円の共通弦
The chord shared by the circles centered on the 2nd and 3rd stations in the input data file (1st argument)
震源の深さを求めるための垂線
The line perpendicular to one of the chords used to determine the depth of the hypocenter

◆使用例(Example)

以下の入力データファイル(sample_data.dat)を使用する。ここで[TAB]はタブを表す。この例において観測点と読み取りデータは全て架空のものである。
Use a data file (sample_data.dat) shown below, where [TAB] indicates a tab. In this example, all stations and pick data are virtual.

(sample_data.dat)

NU.STN_A[TAB]35.12345[TAB]136.78901[TAB]0.0[TAB]23.45[TAB]30.03
NU.STN_B[TAB]34.98765[TAB]137.53197[TAB]0.0[TAB]24.68[TAB]32.10
NU.STN_C[TAB]34.34343[TAB]136.54321[TAB]0.0[TAB]25.52[TAB]33.33

このデータを使用し、仮定する$k$の値を7 km/sから9 km/sまで0.1 km/s刻みで動かして3点法震源決定を行う。これはbashのループ処理を利用した以下のコマンドで実行できる。
Repeat hypocenter determination by 3-point method with assumed $k$ varying from 7 km/s to 9 km/s at an increment of 0.1 km/s. This task is realized by the command below, where loop operation is implemented with bash.

for((k10=70;k10<=90;k10++))
do
k=`echo $k10 | awk '{ printf("%.1f",$1/10.0) }'`
simulate_hypo_3points sample_data.dat result_k$k --k=$k
done

このコマンドを実行するとresult_k7.0からresult_k9.0まで21個のディレクトリが作成され、その中に仮定した各$k$に対する計算結果が格納される。
This command creates 21 directories from result_k7.0 to result_k9.0. Computation results for each assumed $k$ are stored in each directory.

例えばresult_k9.0ディレクトリの下には以下のファイルが作成される。
For example, files listed below are created in the directory result_k9.0.

(result_k9.0/hypo.dat)

latitude=34.736355
longitude=136.997982
depth=36.037779

(result_k9.0/circle1_NU.STN_A.dat)

137.438760[TAB]35.122597
137.438716[TAB]35.131914
137.438474[TAB]35.141229
137.438035[TAB]35.150538
137.437397[TAB]35.159841
…

(result_k9.0/circle2_NU.STN_B.dat)

138.263337[TAB]34.981970
138.263392[TAB]34.992475
138.263225[TAB]35.002979
138.262836[TAB]35.013479
138.262223[TAB]35.023972
…

(result_k9.0/circle3_NU.STN_C.dat)

137.307213[TAB]34.343614
137.307143[TAB]34.354674
137.306840[TAB]34.365731
137.306304[TAB]34.376783
137.305536[TAB]34.387824
…

(result_k9.0/circle_depth.dat)

137.212010[TAB]35.529398
137.202174[TAB]35.531116
137.192302[TAB]35.532693
137.182398[TAB]35.534128
137.172465[TAB]35.535420
…

(result_k9.0/lines.dat)

136.964132[TAB]34.609468
137.212010[TAB]35.529398
>
136.195281[TAB]34.908062
137.178678[TAB]34.696871
>
136.803027[TAB]34.939616
137.302481[TAB]34.416241
>
136.997982[TAB]34.736355
136.613347[TAB]34.805726

一方、result_k7.0ディレクトリの下には以下のファイルが作成される。 hypo.datにおいて「hypocenter not determined」と出力されていること、 circle_depth.datやlines.datが存在しないことから $k=7$ km/sでは震源決定ができなかったことが分かる。この場合であっても観測点を中心とする円のデータは出力される。
In the directory result_k7.0, only the files listed below are created. The content of hypo.dat (a text “hypocenter not determined”) and the absence of circle_depth.dat and lines.dat indicate that the hypocenter was not determined in case of $k=7$ km/s. Even in this case, data of circles centered on stations are created.

(result_k7.0/hypo.dat)

hypocenter not determined

(result_k7.0/circle1_NU.STN_A.dat)

137.294374[TAB]35.123087
137.294328[TAB]35.130334
137.294127[TAB]35.137579
137.293772[TAB]35.144819
137.293263[TAB]35.152054
…

(result_k7.0/circle2_NU.STN_B.dat)

138.100830[TAB]34.983613
138.100857[TAB]34.991784
138.100711[TAB]34.999953
138.100392[TAB]35.008120
138.099899[TAB]35.016281
…

(result_k7.0/circle3_NU.STN_C.dat)

137.137435[TAB]34.343985
137.137363[TAB]34.352587
137.137110[TAB]34.361187
137.136676[TAB]34.369782
137.136062[TAB]34.378369
…

コンピュータプログラムを用いる利点は結果を1つずつ見なくても済むことである。例えば仮定した$k$の値と震源の深さの関係を調べるには以下のコマンドを実行すれば良い。
An advantage of using a computer program is to avoid looking at individual results one by one. For example, the relation between assumed $k$ and estimated hypocenter depth can be derived by the command below.

for((k10=70;k10<=90;k10++))
do
    k=`echo $k10 | awk '{ printf("%.1f",$1/10.0) }'`
    depth=`grep depth result_k$k/hypo.dat | cut -f2 --delimiter='='`
    if [ "$depth" != '' ]
    then
        echo $k $depth
    fi
done

このコマンドを実行すると以下の出力が得られる。
The output of this command is as follows.

7.7 8.917243
7.8 12.927324
7.9 15.992432
8.0 18.586027
8.1 20.884441
8.2 22.976356
8.3 24.913773
8.4 26.730242
8.5 28.448878
8.6 30.086391
8.7 31.655316
8.8 33.165333
8.9 34.624097
9.0 36.037779

これにより、仮定した$k$の値と求まる震源の深さの関係を簡単に調べられるだけでなく、今回のデータでは$k\leq 7.6$ km/sを用いると震源が求まらないことも分かる。
In this way, users can not only easily survey the relation between assumed $k$ and estimated hypocenter depth but also know that the hypocenter is not determined from this data if $k\leq 7.6$ km/s is assumed.

ymaeda_opentoolsにはグラフを描画する機能は無いが、 simulate_hypo_3pointsコマンドから出力される円や線の座標データは外部のプロッティングツールを用いて簡単にプロットできる。 Generic Mapping Tools (version 6.5.0) を用いてプロットするbashスクリプトの例を以下に示す。
Although ymaeda_opentools does not include plotting tools, the data of the coordinates of circles and lines created by simulate_hypo_3points command can easily be plotted using third-party plotting tools. An example of a bash script, using Generic Mapping Tools (version 6.5.0), to plot the results is shown below.

#! /bin/bash -u

# 1つの観測点について観測点位置と観測点を中心とする円をプロットする関数
# A function to plot the location of a station and a circle centered on the station
function TF_plot_for_station () {

    # 観測点位置をプロット
    # Plot a station location
    sed -n ${line_currentStation}p $data_file | awk '{ printf("%f\t%f\n",$3,$2) }' | gmt plot $GMT_R $GMT_J $GMT_XY -St0.5 -G${color_currentStation}

    # 観測点コードをプロット
    # Plot a station code
    sed -n ${line_currentStation}p $data_file | awk '{ printf("%f\t%f\t%s\n",$3,$2-0.04,$1) }' | gmt text $GMT_R $GMT_J $GMT_XY -F+f10p,,${color_currentStation}+jCT

    # 円の座標データのファイル名を取得
    # Survey a file name for the coordinates of a circle
    station=`sed -n ${line_currentStation}p $data_file | awk '{ print $1 }'`
    inputfile=circle${line_currentStation}_${station}.dat

    # 円をプロット
    # Plot a circle
    gmt plot $inputfile $GMT_R $GMT_J $GMT_XY -W0.4,${color_currentStation}
}

# プロットする経度範囲と座標目盛間隔の設定
# Set the coordinate range and tick interval of the longitude axis
Emin=135.5
Emax=138.5
Etics=1.0

# プロットする緯度範囲と座標目盛間隔の設定
# Set the coordinate range and tick interval of the latitude axis
Nmin=33.5
Nmax=36.0
Ntics=1.0

# 紙面上での図の左下の位置(紙の左下端から、cm)
# Location of the bottom left corner of the map on a paper (measured from the bottom left corner of the paper, cm)
x0=2.0
y0=2.0

# 各観測点に対して使用する色
# The color used for each station
color_station1='255/0/0'
color_station2='0/160/0'
color_station3='0/0/255'

# 震源の深さを求めるための半円と線に対して使用する色
# The color used for a half circle and a line to determine hypocenter depth
color_depth='160/0/160'

# 共通オプション
# Common options
GMT_R="-R${Emin}/${Emax}/${Nmin}/${Nmax}"
GMT_J=`echo $Emin $Emax $Nmin $Nmax | awk '{ printf("-Jb%f/%f/%f/%f/1:2000000", ($1+$2)/2.0,($3+$4)/2.0,$3,$4) }'`
GMT_XY="-Xa${x0} -Ya${y0}"

# 仮定した各kについて地図をプロット
# Plot a map for each assumed k
for result_dir in result_k*
do
    # 計算結果のディレクトリに移動
    # Move into the result directory
    cd $result_dir

    # 入出力ファイル名を設定
    # Set input and output file names
    data_file=../sample_data.dat
    outputfile_noExt=map

    # プロットの開始
    # Start plotting
    gmt begin $outputfile_noExt ps

    # 座標軸に用いるフォントサイズの設定
    # Set font size for coordinate axis
    gmt set FONT_ANNOT_PRIMARY 10p

    # 海岸線のプロット
    # Plot the coast line
    gmt coast $GMT_R $GMT_J $GMT_XY -Bxa${Etics} -Bya${Ntics} -BWSen -S191/223/255 -G180/220/180 -W0.4,180/180/180 -Df

    # 1つ目の観測点と円のプロット
    # Plot the 1st station and circle
    line_currentStation=1
    color_currentStation=$color_station1
    TF_plot_for_station

    # 2つ目の観測点と円のプロット
    # Plot the 2nd station and circle
    line_currentStation=2
    color_currentStation=$color_station2
    TF_plot_for_station

    # 3つ目の観測点と円のプロット
    # Plot the 3rd station and circle
    line_currentStation=3
    color_currentStation=$color_station3
    TF_plot_for_station

    # 共通弦のプロット
    # Plot chords
    inputfile=lines.dat
    if [ -e $inputfile ]
    then
        sed -n '1,8p' $inputfile | gmt plot $GMT_R $GMT_J $GMT_XY -W0.4,0/0/0
    fi

    # 震源の深さを求めるための半円のプロット
    # Plot a half circle to determine hypocenter depth
    inputfile=circle_depth.dat
    if [ -e $inputfile ]
    then
        gmt plot $inputfile $GMT_R $GMT_J $GMT_XY -W0.4,${color_depth}
    fi

    # 震源の深さを求めるための垂線のプロット
    # Plot a line perpendicular to a chord to determine hypocenter depth
    inputfile=lines.dat
    if [ -e $inputfile ]
    then
        sed -n '10,11p' $inputfile | gmt plot $GMT_R $GMT_J $GMT_XY -W0.4,${color_depth}
    fi

    # プロットの終了
    # End plotting
    gmt end

    # png形式の画像ファイルに変換
    # Convert to png format
    convert -trim -density 150 $outputfile_noExt.ps $outputfile_noExt.png

    # 次のkに対する結果のプロットのため元の作業ディレクトリに戻る
    # Return back to the original working directory for the next k
    cd ..
done

このスクリプトを実行すると各$k$に対する結果のディレクトリ(result_k∗)の下に map.pngが作成される。先に示した図1, 図2はこのうちの$k=9.0$ km/sおよび$k=7.0$ km/sに対する図である。
This script creates map.png under the result directory for each assumed $k$ (result_k∗). Figs. 1 and 2 shown above correspond to $k=9.0$ and $k_7.0$ km/s, respectively, of this output.

◆計算式(Formula)

以下で赤字が実際にコンピュータプログラム内で計算に用いている式である。
Red colors below show the equations that are actually used in the computer program.

●緯度経度と直交座標系の変換 (Conversion between geographic and cartesian coordinate systems)

このプログラムでは全ての計算を直交座標系で行う。緯度経度と直交座標の変換には関数 latlon2xy, xy2latlon を使用する。はじめに緯度経度で与えられた観測点の座標を関数latlon2xyを用いて直交座標系に変換する。次にその座標を用いて円や線の座標、震源位置を直交座標系で求め、最後にそれらの座標を関数xy2latlonにより緯度経度に変換する。
All computations in this program are conducted in the cartesian coordinate system. Functions latlon2xy and xy2latlon are used for the conversion between geographic and cartesian coordinates. First, the coordinates of stations, given by the latitudes and longitudes, and converted to cartesian coordinates using function latlon2xy. Next, the cartesian coordinates of circles, lines, and the hypocenter are derived, and then converted to the latitudes and longitudes using function xy2latlon.

●観測点を中心とする円 (A circle centered on a station)

直交座標系で表した観測点$i$の座標を$(x_i,y_i)$とする。観測点$i$を中心とする半径$R_i$の円の方程式は \[\begin{equation} (x-x_i)^2+(y-y_i)^2=R_i^2 \label{eq.circle_i} \end{equation}\] と表される。極座標系を用いて \[\begin{equation} \textcolor[named]{Red}{x=x_i+R_i\cos\theta_i} \label{eq.circle_i.x} \end{equation}\] \[\begin{equation} \textcolor[named]{Red}{y=y_i+R_i\sin\theta_i} \label{eq.circle_i.y} \end{equation}\] とおけば(\ref{eq.circle_i})式は満たされる。このプログラムでは$\theta_i$を0°から360°まで1°ずつ動かしながら (\ref{eq.circle_i.x})(\ref{eq.circle_i.y})式を繰り返し用いることで円周上の点の座標を計算する。
Let $(x_i,y_i)$ be the cartesian coordinate of $i$th station. A equation for a circle of radius $R_i$ centered on this station is expressed by Eq. (\ref{eq.circle_i}). A polar coordinate expression given by Eqs. (\ref{eq.circle_i.x}) and (\ref{eq.circle_i.y}) satisfies this equation. This program repeatedly applies Eqs. (\ref{eq.circle_i.x}) and (\ref{eq.circle_i.y}) for $\theta_i$ from 0° to 360° at an increment of 1° to compute the coordinates of the circumference.

●共通弦 (A chord shared by two circles)

観測点$j$を中心とする半径$R_j$の円の方程式は \[\begin{equation} (x-x_j)^2+(y-y_j)^2=R_j^2 \label{eq.circle_j} \end{equation}\] と表される。 (\ref{eq.circle_i})(\ref{eq.circle_j})を同時に満たす$(x,y)$が観測点$i$, $j$を中心とする半径$R_i$, $R_j$の円の共通弦の端点となる。
A circle of radius $R_j$ centered on $j$th station is expressed by Eq. (\ref{eq.circle_j}). Locations $(x,y)$ that satisfy Eqs. (\ref{eq.circle_i}) and (\ref{eq.circle_j}) are the end points of the chord shared by circles of radii $R_i$ and $R_j$ centered on $i$th and $j$th stations.

観測点$i$から共通弦の端点に向かうベクトルを観測点$i$から$j$に向かう方向の成分とそれに直交する成分に分けて \[\begin{equation} \begin{pmatrix} x-x_i \\ y-y_i \end{pmatrix} = \alpha_{ij} \begin{pmatrix} x_j-x_i \\ y_j-y_i \end{pmatrix} +\beta_{ij} \begin{pmatrix} y_j-y_i \\ -(x_j-x_i) \end{pmatrix} \label{eq.chord.vector.i} \end{equation}\] のように表す。このとき \[\begin{equation} \begin{pmatrix} x-x_j \\ y-y_j \end{pmatrix} = \begin{pmatrix} x-x_i \\ y-y_i \end{pmatrix} - \begin{pmatrix} x_j-x_i \\ y_j-y_i \end{pmatrix} = (\alpha_{ij}-1) \begin{pmatrix} x_j-x_i \\ y_j-y_i \end{pmatrix} +\beta_{ij} \begin{pmatrix} y_j-y_i \\ -(x_j-x_i) \end{pmatrix} \label{eq.chord.vector.j} \end{equation}\] である。(\ref{eq.chord.vector.i})を(\ref{eq.circle_i})に代入すると \[\begin{eqnarray} R_i^2 &=& (x-x_i)^2+(y-y_i)^2 \nonumber \\ &=& \left[\alpha_{ij}(x_j-x_i)+\beta_{ij}(y_j-y_i)\right]^2 +\left[\alpha_{ij}(y_j-y_i)-\beta_{ij}(x_j-x_i)\right]^2 \nonumber \\ &=& \alpha_{ij}^2(x_j-x_i)^2 +2\alpha_{ij}\beta_{ij}(x_j-x_i)(y_j-y_i) +\beta_{ij}^2(y_j-y_i)^2 \nonumber \\ & & +\alpha_{ij}^2(y_j-y_i)^2 -2\alpha_{ij}\beta_{ij}(x_j-x_i)(y_j-y_i) +\beta_{ij}^2(x_j-x_i)^2 \nonumber \\ &=& (\alpha_{ij}^2+\beta_{ij}^2)\left[(x_j-x_i)^2+(y_j-y_i)^2\right] \nonumber \\ &=& (\alpha_{ij}^2+\beta_{ij}^2)R_{ij}^2 \label{eq.circle_i.vector.derive1} \end{eqnarray}\] となるので \[\begin{equation} \alpha_{ij}^2+\beta_{ij}^2 =\frac{R_i^2}{R_{ij}^2} \label{eq.circle_i.vector} \end{equation}\] を得る。ここで \[\begin{equation} \textcolor[named]{Red}{R_{ij}\equiv \sqrt{(x_j-x_i)^2+(y_j-y_i)^2}} \label{eq.R_ij} \end{equation}\] は観測点$i$, $j$間の距離を表す。同様に(\ref{eq.chord.vector.j})を(\ref{eq.circle_j})に代入すると \[\begin{equation} (\alpha_{ij}-1)^2+\beta_{ij}^2 =\frac{R_j^2}{R_{ij}^2} \label{eq.circle_j.vector} \end{equation}\] を得る。(\ref{eq.circle_j.vector})$-$(\ref{eq.circle_i.vector})より \[\begin{equation} (\alpha_{ij}-1)^2-\alpha_{ij}^2 =-2\alpha_{ij}+1 =\frac{R_j^2-R_i^2}{R_{ij}^2} \label{eq.alpha.derive1} \end{equation}\] が得られ、変形して \[\begin{equation} \textcolor[named]{Red}{ \alpha_{ij}=\frac{1}{2}\left(1-\frac{R_j^2-R_i^2}{R_{ij}^2}\right) } \label{eq.alpha} \end{equation}\] を得る。
Let a vector from $i$th station to an end of the chord be expressed by a summation of parallel to perpendicular components to a vector from $i$th to $j$th stations, as Eq. (\ref{eq.chord.vector.i}). Eq. (\ref{eq.chord.vector.j}) is derived. Inserting Eq. (\ref{eq.chord.vector.i}) into (\ref{eq.circle_i}) results in (\ref{eq.circle_i.vector.derive1}), which can be arranged as (\ref{eq.circle_i.vector}), where $R_{ij}$ is defined by Eq. (\ref{eq.R_ij}). Similarly, inserting Eq. (\ref{eq.chord.vector.j}) into (\ref{eq.circle_j}) results in (\ref{eq.circle_j.vector}). Subtracting Eq. (\ref{eq.circle_i.vector}) from (\ref{eq.circle_j.vector}) results in Eq. (\ref{eq.alpha.derive1}), which can be arranged as (\ref{eq.alpha}).

$\alpha_{ij}\geq R_i/R_{ij}$のとき、 (\ref{eq.circle_i.vector})式を満たす$\beta_{ij}$は存在しない。これは仮定した$k$の値が小さいために円と円が交差せず共通弦が存在しないことを表す。 $\alpha_{ij}< R_i/R_{ij}$のときは(\ref{eq.circle_i.vector})式より \[\begin{eqnarray} \textcolor[named]{Red}{\beta_{ij}} &\textcolor[named]{Red}{=}& \textcolor[named]{Red}{\pm\sqrt{\frac{R_i^2}{R_{ij}^2}-\alpha_{ij}^2}} \nonumber \\ &=& \pm\sqrt{\frac{R_i^2}{R_{ij}^2} -\frac{1}{4}\left(1-\frac{R_j^2-R_i^2}{R_{ij}^2}\right)^2} \nonumber \\ &=& \pm\sqrt{\frac{R_i^2}{R_{ij}^2} -\frac{1}{4}\left[1-2\frac{R_j^2-R_i^2}{R_{ij}^2} +\frac{(R_j^2-R_i^2)^2}{R_{ij}^4}\right]} \nonumber \\ &=& \pm\sqrt{-\frac{1}{4} +\frac{1}{2}\frac{R_i^2+R_j^2}{R_{ij}^2} -\frac{1}{4}\frac{(R_j^2-R_i^2)^2}{R_{ij}^4}} \label{eq.beta} \end{eqnarray}\] である。
If $\alpha_{ij}\geq R_i/R_{ij}$, there is no $\beta_{ij}$ that satisfies Eq. (\ref{eq.circle_i.vector}). This situation occurs because of no intersection of the circles because of too small assumption of $k$. If $\alpha_{ij}< R_i/R_{ij}$, $\beta_{ij}$ is derived from Eq. (\ref{eq.circle_i.vector}) as (\ref{eq.beta}).

(\ref{eq.beta})式の$\beta_{ij}$は$i$と$j$を入れ替えても変わらない。すなわち$\beta_{ij}$は$i$, $j$に関して対称な形になっている。一方、(\ref{eq.alpha})式の$\alpha_{ij}$は$i$, $j$に関して非対称である。これは元々の定義が(\ref{eq.chord.vector.i})式のように $i$と$j$に関して非対称な形になっていたからである。 (\ref{eq.alpha})(\ref{eq.beta})を(\ref{eq.chord.vector.i})に代入すると \[\begin{eqnarray} \begin{pmatrix} x \\ y \end{pmatrix} &=& \begin{pmatrix} x_i \\ y_i \end{pmatrix} +\alpha_{ij} \begin{pmatrix} x_j-x_i \\ y_j-y_i \end{pmatrix} +\beta_{ij} \begin{pmatrix} y_j-y_i \\ -(x_j-x_i) \end{pmatrix} \nonumber \\ &=& \begin{pmatrix} x_i \\ y_i \end{pmatrix} +\frac{1}{2}\left(1-\frac{R_j^2-R_i^2}{R_{ij}^2}\right) \begin{pmatrix} x_j-x_i \\ y_j-y_i \end{pmatrix} +\beta_{ij} \begin{pmatrix} y_j-y_i \\ -(x_j-x_i) \end{pmatrix} \nonumber \\ &=& \frac{1}{2}\left(1+\frac{R_j^2-R_i^2}{R_{ij}^2}\right) \begin{pmatrix} x_i \\ y_i \end{pmatrix} +\frac{1}{2}\left(1+\frac{R_i^2-R_j^2}{R_{ij}^2}\right) \begin{pmatrix} x_j \\ y_j \end{pmatrix} \nonumber \\ & & \pm\sqrt{-\frac{1}{4} +\frac{1}{2}\frac{R_i^2+R_j^2}{R_{ij}^2} -\frac{1}{4}\frac{(R_j^2-R_i^2)^2}{R_{ij}^4}} \begin{pmatrix} y_j-y_i \\ -(x_j-x_i) \end{pmatrix} \label{eq.chord} \end{eqnarray}\] と$i$, $j$に関して対称な形になる。
Swapping $i$ and $j$ does not change the value of $\beta_{ij}$ (Eq. \ref{eq.beta}), indicating that $\beta_{ij}$ is symmetric with respect to $i$ and $j$. On the contrary, $\alpha_{ij}$ is asymmetric with respect to $i$ and $j$ because of the asymmetric definition of Eq. (\ref{eq.chord.vector.i}). Inserting Eqs. (\ref{eq.alpha}) and (\ref{eq.beta}) into (\ref{eq.chord.vector.i}) results in (\ref{eq.chord}), which is symmetric with respect to $i$ and $j$.

(\ref{eq.chord})式の右辺第2項までの部分を \[\begin{eqnarray} \textcolor[named]{Red}{ \begin{pmatrix} x_{ij} \\ y_{ij} \end{pmatrix} } &\equiv& \frac{1}{2}\left(1+\frac{R_j^2-R_i^2}{R_{ij}^2}\right) \begin{pmatrix} x_i \\ y_i \end{pmatrix} +\frac{1}{2}\left(1+\frac{R_i^2-R_j^2}{R_{ij}^2}\right) \begin{pmatrix} x_j \\ y_j \end{pmatrix} \nonumber \\ &\textcolor[named]{red}{=}& \textcolor[named]{Red}{ \frac{1}{2} \begin{pmatrix} x_j+x_i \\ y_j+y_i \end{pmatrix} -\frac{R_j^2-R_i^2}{2R_{ij}^2} \begin{pmatrix} x_j-x_i \\ y_j-y_i \end{pmatrix} } \label{eq.xy_ij} \end{eqnarray}\] とおく。この$(x_{ij},y_{ij})$は共通弦の中点の座標であり、共通弦と観測点$i$, $j$を結ぶ線の交点の座標でもある。 $(x_{ij},y_{ij})$を用いると(\ref{eq.chord})式は \[\begin{equation} \textcolor[named]{Red}{ \begin{pmatrix} x \\ y \end{pmatrix} = \begin{pmatrix} x_{ij} \\ y_{ij} \end{pmatrix} +\beta_{ij} \begin{pmatrix} y_j-y_i \\ -(x_j-x_i) \end{pmatrix} } \label{eq.chord.simplified} \end{equation}\] という簡単な形で表現できる。
Let $(x_{ij},y_{ij})$ be the summation of the first two terms of the right hand side of Eq. (\ref{eq.chord}), i.e., (\ref{eq.xy_ij}). This coordinate corresponds to the mid point of the chord and also an intersection of the chord and a line connecting $i$th and $j$th stations. Using $(x_{ij},y_{ij})$, Eq. (\ref{eq.chord}) is simplified as (\ref{eq.chord.simplified}).

(\ref{eq.chord.simplified})式は共通弦の両端の点の座標を表す式である。これらの端点を結んだ線が共通弦であるので \[\begin{equation} \begin{pmatrix} x \\ y \end{pmatrix} = \begin{pmatrix} x_{ij} \\ y_{ij} \end{pmatrix} +\gamma_{ij} \begin{pmatrix} y_j-y_i \\ -(x_j-x_i) \end{pmatrix} \label{eq.chord.line} \end{equation}\] すなわち \[\begin{equation} x-x_{ij}=\gamma(y_j-y_i) \label{eq.chord.line.x} \end{equation}\] \[\begin{equation} y-y_{ij}=-\gamma(x_j-x_i) \label{eq.chord.line.y} \end{equation}\] が共通弦(線)を表す方程式である。2式を \[\begin{equation} \gamma=\frac{x-x_{ij}}{y_j-y_i}=-\frac{y-y_{ij}}{x_j-x_i} \label{eq.chord.gamma} \end{equation}\] と変形して$\gamma$を消去すると \[\begin{equation} (x_j-x_i)(x-x_{ij})+(y_j-y_i)(y-y_{ij})=0 \label{eq.chord.line.no_gamma} \end{equation}\] という関係式が得られる。
Eq. (\ref{eq.chord.simplified}) represents the coordinates of the end points of the chord. The chord is a line that connects these two points, and thus can be expressed as (\ref{eq.chord.line}), i.e., (\ref{eq.chord.line.x}) and (\ref{eq.chord.line.y}). Arranging these equation as (\ref{eq.chord.gamma}) to remove $\gamma$ results in (\ref{eq.chord.line.no_gamma}).

●震央 (Epicenter)

震央は3つの共通弦の交点である。実際には2つの共通弦の交点を求めれば良い。 (\ref{eq.chord.line.no_gamma})式より観測点1,2を中心とする円の共通弦は \[\begin{equation} (x_2-x_1)(x-x_{12})+(y_2-y_1)(y-y_{12})=0 \label{eq.chord.12} \end{equation}\] 観測点1,3を中心とする円の共通弦は \[\begin{equation} (x_3-x_1)(x-x_{13})+(y_3-y_1)(y-y_{13})=0 \label{eq.chord.13} \end{equation}\] であり、2式をともに満たす$(x,y)$が震央である。以下では震央を$(x_o,y_o)$と表す。 2式を連立させると \[\begin{equation} \begin{pmatrix} x_2-x_1 & y_2-y_1 \\ x_3-x_1 & y_3-y_1 \end{pmatrix} \begin{pmatrix} x_o \\ y_o \end{pmatrix} = \begin{pmatrix} (x_2-x_1)x_{12}+(y_2-y_1)y_{12} \\ (x_3-x_1)x_{13}+(y_3-y_1)y_{13} \end{pmatrix} \label{eq.chord.intersection} \end{equation}\] であるので、その解はクラメルの公式より \[\begin{equation} \textcolor[named]{Red}{x_o=\frac{\Delta_1}{\Delta}} \label{eq.hypo.x} \end{equation}\] \[\begin{equation} \textcolor[named]{Red}{y_o=\frac{\Delta_2}{\Delta}} \label{eq.hypo.y} \end{equation}\] \[\begin{eqnarray} \textcolor[named]{Red}{\Delta} &\equiv& \begin{vmatrix} x_2-x_1 & y_2-y_1 \\ x_3-x_1 & y_3-y_1 \end{vmatrix} \nonumber \\ &=& (x_2-x_1)(y_3-y_1)-(x_3-x_1)(y_2-y_1) \nonumber \\ &=& (x_2y_3-x_2y_1-x_1y_3+x_1y_1)-(x_3y_2-x_3y_1-x_1y_2+x_1y_1) \nonumber \\ &\textcolor[named]{Red}{=}& \textcolor[named]{Red}{(x_1y_2-x_2y_1)+(x_2y_3-x_3y_2)+(x_3y_1-x_1y_3)} \label{eq.Delta} \end{eqnarray}\] \[\begin{eqnarray} \Delta_1 &\equiv& \begin{vmatrix} (x_2-x_1)x_{12}+(y_2-y_1)y_{12} & y_2-y_1 \\ (x_3-x_1)x_{13}+(y_3-y_1)y_{13} & y_3-y_1 \end{vmatrix} \nonumber \\ &=& [(x_2-x_1)x_{12}+(y_2-y_1)y_{12}](y_3-y_1) -[(x_3-x_1)x_{13}+(y_3-y_1)y_{13}](y_2-y_1) \label{eq.Delta1} \end{eqnarray}\] \[\begin{eqnarray} \Delta_2 &\equiv& \begin{vmatrix} x_2-x_1 & (x_2-x_1)x_{12}+(y_2-y_1)y_{12} \\ x_3-x_1 & (x_3-x_1)x_{13}+(y_3-y_1)y_{13} \end{vmatrix} \nonumber \\ &=& (x_2-x_1)[(x_3-x_1)x_{13}+(y_3-y_1)y_{13}] -(x_3-x_1)[(x_2-x_1)x_{12}+(y_2-y_1)y_{12}] \label{eq.Delta2} \end{eqnarray}\] となる。
The epicenter is the cross point of the three chords, and actually can be obtained by solving the cross point of arbitrary two chords. From Eq. (\ref{eq.chord.line.no_gamma}), the chord shared by the circles centered on 1st and 2nd stations is expressed by Eq. (\ref{eq.chord.12}), and the chord shared by the circles centered on 1st and 3rd stations is expressed by Eq. (\ref{eq.chord.13}). The location $(x,y)$ that satisfies the two equations is the epicenter, which is expressed as $(x_o,y_o)$ in the following. Combining the two equations, we obtain (\ref{eq.chord.intersection}). Its solution is given by Eqs. (\ref{eq.hypo.x})-(\ref{eq.Delta2}) according to Cramer's formula.

観測点1, 2, 3に関して対称な形にするために (\ref{eq.Delta1})(\ref{eq.Delta2})を更に変形する。まず(\ref{eq.xy_ij})式より \[\begin{eqnarray} (x_2-x_1)x_{12}+(y_2-y_1)y_{12} &=& \frac{1}{2}(x_2-x_1)(x_2+x_1) -\frac{R_2^2-R_1^2}{2R_{12}^2}(x_2-x_1)^2 \nonumber \\ & & +\frac{1}{2}(y_2-y_1)(y_2+y_1) -\frac{R_2^2-R_1^2}{2R_{12}^2}(y_2-y_1)^2 \nonumber \\ &=& \frac{1}{2}(x_2^2-x_1^2+y_2^2-y_1^2) \nonumber \\ & & -\frac{R_2^2-R_1^2}{2R_{12}^2}\left[(x_2-x_1)^2+(y_2-y_1)^2\right] \nonumber \\ &=& \frac{1}{2}(x_2^2+y_2^2)-\frac{1}{2}(x_1^2+y_1^2) -\frac{R_2^2-R_1^2}{2R_{12}^2}R_{12}^2 \nonumber \\ &=& \frac{1}{2}\left(x_2^2+y_2^2-R_2^2\right) -\frac{1}{2}\left(x_1^2+y_1^2-R_1^2\right) \nonumber \\ &=& \delta_2-\delta_1 \label{eq.Delta.term12} \end{eqnarray}\] 同様に \[\begin{equation} (x_3-x_1)x_{13}+(y_3-y_1)y_{13}=\delta_3-\delta_1 \label{eq.Delta.term13} \end{equation}\] である。ここで \[\begin{equation} \textcolor[named]{Red}{ \delta_i\equiv \frac{1}{2}\left(x_i^2+y_i^2-R_i^2\right) } \label{eq.delta_i} \end{equation}\] とおいた。(\ref{eq.Delta.term12})(\ref{eq.Delta.term13})を (\ref{eq.Delta1})(\ref{eq.Delta2})に代入すると \[\begin{eqnarray} \textcolor[named]{Red}{\Delta_1} &=& (\delta_2-\delta_1)(y_3-y_1) -(\delta_3-\delta_1)(y_2-y_1) \nonumber \\ &=& (\delta_2y_3-\delta_2y_1-\delta_1y_3+\delta_1y_1) -(\delta_3y_2-\delta_3y_1-\delta_1y_2+\delta_1y_1) \nonumber \\ &\textcolor[named]{Red}{=}& \textcolor[named]{Red}{ (\delta_1y_2-\delta_2y_1) +(\delta_2y_3-\delta_3y_2) +(\delta_3y_1-\delta_1y_3) } \label{eq.Delta1.symmetric} \end{eqnarray}\] \[\begin{eqnarray} \textcolor[named]{Red}{\Delta_2} &=& (x_2-x_1)(\delta_3-\delta_1)-(x_3-x_1)(\delta_2-\delta_1) \nonumber \\ &=& (x_2\delta_3-x_2\delta_1-x_1\delta_3+x_1\delta_1) -(x_3\delta_2-x_3\delta_1-x_1\delta_2+x_1\delta_1) \nonumber \\ &\textcolor[named]{Red}{=}& \textcolor[named]{Red}{ (x_1\delta_2-x_2\delta_1) +(x_2\delta_3-x_3\delta_2) +(x_3\delta_1-x_1\delta_2) } \label{eq.Delta2.symmetric} \end{eqnarray}\] となって対称的な形の解が得られる。
We further arrange Eqs. (\ref{eq.Delta1}) and (\ref{eq.Delta2}) to make them symmetric with respect to stations 1, 2, and 3. Eq. (\ref{eq.xy_ij}) gives (\ref{eq.Delta.term12}) and (\ref{eq.Delta.term13}), where $\delta_i$ is defined by Eq. (\ref{eq.delta_i}). Inserting Eqs. (\ref{eq.Delta.term12}) and (\ref{eq.Delta.term13}) into (\ref{eq.Delta1}) and (\ref{eq.Delta2}) results in (\ref{eq.Delta1.symmetric}) and (\ref{eq.Delta2.symmetric}), which are symmetric with respect to stations 1, 2, and 3.

●震源の深さを求める半円 (A half circle to determine the depth of the hypocenter)

(\ref{eq.chord.simplified})式より、観測点1, 2を中心とする円の共通弦を直径とする半円は中心点が$(x_{12},y_{12})$、半径が$\beta_{12}R_{12}$である。このような半円は2つあり、そのうち観測点1の側を向いた半円を求めれば良い。
Eq. (\ref{eq.chord.simplified}) indicates that $(x_{12},y_{12})$ and $\beta_{12}R_{12}$ are the center and radius, respectively, of the half circle whose bottom line is the chord shared by the circles centered on 1st and 2nd stations. There are two half circles that satisfy these requirements, and choose the one that is oriented to the 1st station.

観測点1の側を向くという条件の数式による表現を得るためにまずは観測点2から観測点1に向かうベクトル \[\begin{equation} \myvector{v_{2\rightarrow 1}}\equiv \begin{pmatrix} x_1-x_2 \\ y_1-y_2 \end{pmatrix} \label{eq.v21} \end{equation}\] を考える。$\myvector{v_{2\rightarrow 1}}$の方向を $x$軸から反時計回りに測った角度を$\theta_{2\rightarrow 1}$とすると \[\begin{equation} \tan\theta_{2\rightarrow 1}=\frac{y_1-y_2}{x_1-x_2} \label{eq.theta21} \end{equation}\] である。観測点1から観測点2に向かうベクトルの角度$\theta_{1\rightarrow 2}$も (\ref{eq.theta21})を満たしてしまうが、C言語プログラムにおいては関数 atan2(y1-y2,x1-x2) を用いることで両者は区別され、$\theta_{2\rightarrow 1}$が正しく求まる。
Now, let us consider a mathematical expression for the condition that the half circle orients to the 1st station. Consider a vector from the 2nd to 1st station (Eq. \ref{eq.v21}) and let $\theta_{2\rightarrow 1}$ be the direction of this vector measured counterclockwise from the $x$-axis. It satisfies Eq. (\ref{eq.theta21}). The angle $\theta_{1\rightarrow 2}$ of a vector from the 1st to 2nd stations also satisfies Eq. (\ref{eq.theta21}), but they can be distinguished in C language using a function atan2(y1-y2,x1-x2), which correctly returns $\theta_{2\rightarrow 1}$.

この$\theta_{2\rightarrow 1}$から±90°の範囲が観測点1の側を向いた半円であるので \[\begin{equation} \textcolor[named]{Red}{ x=x_{12}+\beta_{12}R_{12}\cos(\theta_{2\rightarrow 1}+\theta) } \label{eq.circle_for_depth.x} \end{equation}\] \[\begin{equation} \textcolor[named]{Red}{ y=y_{12}+\beta_{12}R_{12}\sin(\theta_{2\rightarrow 1}+\theta) } \label{eq.circle_for_depth.y} \end{equation}\] として$\theta$を-90°から+90-°まで1°ずつ動かしながら (\ref{eq.circle_for_depth.x})(\ref{eq.circle_for_depth.y})式を繰り返し用いれば半円の円周上の点の座標リストが得られる。
The half circle oriented to the 1st station is obtained by using the range of the angle within ±90° from this $\theta_{2\rightarrow 1}$. Therefore, the coordinates of the circumference of the half circle are computed by repeatedly applying Eqs. (\ref{eq.circle_for_depth.x}) and (\ref{eq.circle_for_depth.y}) for $\theta$ from -90° to +90° with an increment of 1°.

●震源の深さを求める垂線 (A line perpendicular to a chord to determine the depth of the hypocenter)

震央から上記の半円に垂線を引く。その一方の端は震央$(x_o,y_o)$である。他方の端の座標を求める。垂線は$\myvector{v_{2\rightarrow 1}}$に平行であるのでその向きは$\theta_{2\rightarrow 1}$方向であり、垂線上の任意の点の座標は \[\begin{equation} \textcolor[named]{Red}{x=x_o+\epsilon\cos\theta_{2\rightarrow 1}} \label{eq.line_for_depth.x} \end{equation}\] \[\begin{equation} \textcolor[named]{Red}{y=y_o+\epsilon\sin\theta_{2\rightarrow 1}} \label{eq.line_for_depth.y} \end{equation}\] と表される。一方、円周上の点では \[\begin{equation} (x-x_{12})^2+(y-y_{12})^2=\beta_{12}^2R_{12}^2 \label{eq.circle_for_depth} \end{equation}\] が成り立つ。 (\ref{eq.line_for_depth.x})(\ref{eq.line_for_depth.y})式と (\ref{eq.circle_for_depth})式を同時に満たす$(x,y)$が円周と垂線の交点である。その座標を求めるため、 (\ref{eq.line_for_depth.x})(\ref{eq.line_for_depth.y})を (\ref{eq.circle_for_depth})に代入すると \[\begin{eqnarray} \beta_{12}^2R_{12}^2 &=& (x-x_{12})^2+(y-y_{12})^2 \nonumber \\ &=& (x_o-x_{12}+\epsilon\cos\theta_{2\rightarrow 1})^2 +(y_o-y_{12}+\epsilon\sin\theta_{2\rightarrow 1})^2 \nonumber \\ &=& (x_o-x_{12})^2 +2(x_o-x_{12})\epsilon\cos\theta_{2\rightarrow 1} +\epsilon^2\cos^2\theta_{2\rightarrow 1} \nonumber \\ & & +(y_o-y_{12})^2 +2(y_o-y_{12})\epsilon\sin\theta_{2\rightarrow 1} +\epsilon^2\sin^2\theta_{2\rightarrow 1} \nonumber \\ &=& \epsilon^2 +2[(x_o-x_{12})\cos\theta_{2\rightarrow 1} +(y_o-y_{12})\sin\theta_{2\rightarrow 1}]\epsilon +(x_o-x_{12})^2+(y_o-y_{12})^2 \label{eq.line_for_depth.end.derive1} \end{eqnarray}\] 変形して \[\begin{eqnarray} \epsilon^2 +2[(x_o-x_{12})\cos\theta_{2\rightarrow 1} +(y_o-y_{12})\sin\theta_{2\rightarrow 1}]\epsilon & & \nonumber \\ +(x_o-x_{12})^2+(y_o-y_{12})^2-\beta_{12}^2R_{12}^2 &=& 0 \label{eq.line_for_depth.end.derive2} \end{eqnarray}\] を得る。
Draw a line perpendicular to the chord from the epicenter to the circle derived above. Its one end is the epicenter $(x_o,y_o)$. Let us derive the coordinate of the other end. Because the line is parallel to $\myvector{v_{2\rightarrow 1}}$ and thus oriented to $\theta_{2\rightarrow 1}$, arbitrary points on the line are expressed by Eqs. (\ref{eq.line_for_depth.x}) and (\ref{eq.line_for_depth.y}). Arbitrary points on the circumference of the circle are expressed by Eq. (\ref{eq.circle_for_depth}). The intersect of the line and circumference satisfies all these equations, and can be derived by inserting Eqs. (\ref{eq.line_for_depth.x}) and (\ref{eq.line_for_depth.y}) into (\ref{eq.circle_for_depth}), resulting in Eq. (\ref{eq.line_for_depth.end.derive1}), which is arranged as (\ref{eq.line_for_depth.end.derive2}).

(\ref{eq.line_for_depth.end.derive2})式の解を求める前に左辺第2項について考える。 \[\begin{equation} (x_o-x_{12})\cos\theta_{2\rightarrow 1} +(y_o-y_{12})\sin\theta_{2\rightarrow 1} = \begin{pmatrix} x_o-x_{12} \\ y_o-y_{12} \end{pmatrix} \cdot \begin{pmatrix} \cos\theta_{2\rightarrow 1} \\ \sin\theta_{2\rightarrow 1} \end{pmatrix} \label{eq.line_for_depth.end.derive2.term2} \end{equation}\] と変形してみると(\ref{eq.line_for_depth.end.derive2})式の左辺第2項の係数はベクトル$(x_o-x_{12},y_o-y_{12})$と $(\cos\theta_{2\rightarrow 1},\sin\theta_{2\rightarrow 1})$の内積の2倍であることが分かる。 $(x_o,y_o)$は震央、$(x_{12},y_{12})$は共通弦の中点であるのでこれらはともに共通弦上に位置し、 $(x_o-x_{12},y_o-y_{12})$は共通弦に平行なベクトルである。一方、$(\cos\theta_{2\rightarrow 1},\sin\theta_{2\rightarrow 1})$は共通弦に垂直な方向を向いたベクトルである。したがって両者の内積は0になる。このことを用いると (\ref{eq.line_for_depth.end.derive2})式の左辺第2項が0であることが分かるので (\ref{eq.line_for_depth.end.derive2})式は \[\begin{equation} \epsilon^2+(x_o-x_{12})^2+(y_o-y_{12})^2-\beta_{12}^2R_{12}^2=0 \label{eq.line_for_depth.end.derive3} \end{equation}\] と簡単になり、その解は \[\begin{equation} \epsilon=\pm\sqrt{\beta_{12}^2R_{12}^2-(x_o-x_{12})^2-(y_o-y_{12})^2} \label{eq.line_for_depth.end.derive4} \end{equation}\] と求まる。
Before solving Eq. (\ref{eq.line_for_depth.end.derive2}), let us consider the 2nd term of the left hand side of this equation. Eq. (\ref{eq.line_for_depth.end.derive2.term2}) shows that the coefficient of the the 2nd term of the left hand side of Eq. (\ref{eq.line_for_depth.end.derive2}) is twice the inner product of vectors $(x_o-x_{12},y_o-y_{12})$ and $(\cos\theta_{2\rightarrow 1},\sin\theta_{2\rightarrow 1})$. Remind that $(x_o,y_o)$ is the epicenter and $(x_{12},y_{12})$ is the mid point of the chord. Both of them are on the chord, indicating that the vector $(x_o-x_{12},y_o-y_{12})$ is parallel to the chord. The vector $(\cos\theta_{2\rightarrow 1},\sin\theta_{2\rightarrow 1})$ is perpendicular to the chord. Therefore, the inner product is zero, which means that the 2nd term of the left hand side of Eq. (\ref{eq.line_for_depth.end.derive2}) is zero. Using this conclusion, Eq. (\ref{eq.line_for_depth.end.derive2}) is simplified as (\ref{eq.line_for_depth.end.derive3}), and its solution is given by Eq. (\ref{eq.line_for_depth.end.derive4}).

$\epsilon$の候補が2つ出てきたのはこれまでの導出過程において「半円」であることを用いておらず垂線と「円」の交点を求めたためである。前節で考えた半円は観測点1の側を向いているので (\ref{eq.line_for_depth.x})(\ref{eq.line_for_depth.y})式において $\epsilon>0$のものを考えれば良く、 \[\begin{equation} \textcolor[named]{Red}{ \epsilon=\sqrt{\beta_{12}^2R_{12}^2-(x_o-x_{12})^2-(y_o-y_{12})^2} } \label{eq.epsilon} \end{equation}\] が得られる。この$\epsilon$を (\ref{eq.line_for_depth.x})(\ref{eq.line_for_depth.y})式に代入すれば震源の深さを求める垂線の震央ではない方の端点が求まる。
Two candidates of $\epsilon$ were derived because the “half” circle has not been taken into account in the derivation above; we solved the intersection of the line and circle. The half circle considered in the previous section is oriented to the 1st station, which introduces a constraint $\epsilon>0$ in Eqs. (\ref{eq.line_for_depth.x}) and (\ref{eq.line_for_depth.y}). Therefore, we obtain Eq. (\ref{eq.epsilon}). Inserting this $\epsilon$ into Eqs. (\ref{eq.line_for_depth.x}) and (\ref{eq.line_for_depth.y}), we have the coordinate of the end point of the line to the opposite side of the epicenter used to determine the depth of the hypocenter.

simulate_hypo_3pointsコマンド マニュアル