まず具体的に計算をするために(3)式を(1)(2)式により書き換える。
(1)式により、
\(\displaystyle \frac{\partial Vx}{\partial x}＝-cosθ, \) \(\displaystyle \frac{\partial Vx}{\partial y} = -sinθ, \)\(\displaystyle \frac{\partial Vx}{\partial θ}=-(Xi-x)sinθ+(Yi-y)cosθ=-a, \)\(\displaystyle \frac{\partial Vx}{\partial D}=-K \) --- (5)
(2)式により、
\(\displaystyle \frac{\partial Vy}{\partial x}＝sinθ, \) \(\displaystyle \frac{\partial Vy}{\partial y} = -cosθ, \)\(\displaystyle \frac{\partial Vy}{\partial θ}=-(Yi-y)sinθ-(Xi-x)cosθ=-b, \)\(\displaystyle \frac{\partial Vy}{\partial D}=-K \) ---- (6)
(5)式により(3)式は、
\(Vx = Vx0(x0,y0,θ0,D0) -cosθΔx -sinθΔy -aΔθ-kΔD \)--- (7)
(6)式により(3)式は、
\(Vy= Vy0(x0,y0,θ0,D0) +sinθΔx -cosθΔy -bΔθ-kΔD \)--- (8)
(7)式より正規方程式を作る。(東西方向の条坊に適用する)
\(\displaystyle \frac{\partial ΣVx^2}{\partial x} =\)
\(\displaystyle \frac{\partial Σ(Vx0(x0,y0,θ0,D0) -cosθΔx -sinθΔy -aΔθ-KΔD)^2}{\partial x} \)
\( = Σ- 2cosθ(Vx0(x0,y0,θ0,D0) -cosθΔx -sinθΔy -aΔθ-KΔD )= 0 \) --- (9)
同様に
\(\displaystyle \frac{\partial ΣVx^2}{\partial y} \)
\( = Σ- 2sinθ(Vx0(x0,y0,θ0,D0) -cosθΔx -sinθΔy -aΔθ-KΔD )= 0 \) --- (10)
\(\displaystyle \frac{\partial ΣVx^2}{\partial θ}\)
\( = Σ- 2a(Vx0(x0,y0,θ0,D0) -cosθΔx -sinθΔy -aΔθ-KΔD )= 0 \) --- (11)
\(\displaystyle \frac{\partial ΣVx^2}{\partial D} \)
\( = Σ- 2K(Vx0(x0,y0,θ0,D0) -cosθΔx -sinθΔy -aΔθ-KΔD )= 0 \) --- (12)
同様に(8)式により正規方程式を作る(南北方向の条坊に適用する)
\(\displaystyle \frac{\partial ΣVy^2}{\partial x} \)
\( = Σ- 2sinθ(Vy0(x0,y0,θ0,D0) +sinθΔx -cosθΔy -bΔθ-KΔD )= 0 \) --- (13)
\(\displaystyle \frac{\partial ΣVy^2}{\partial y} \)
\( = Σ- 2cosθ(Vy0(x0,y0,θ0,D0) +sinθΔx -cosθΔy -bΔθ-KΔD )= 0 \) --- (14)
\(\displaystyle \frac{\partial ΣVy^2}{\partial θ}\)
\( = Σ- 2b(Vy0(x0,y0,θ0,D0) +sinθΔx -cosθΔy -bΔθ-KΔD)= 0 \) --- (15)
\(\displaystyle \frac{\partial ΣVy^2}{\partial D} \)
\( = Σ- 2K(Vy0(x0,y0,θ0,D0) +sinθΔx -cosθΔy -bΔθ-KΔD)= 0 \) --- (16)
なお"-2"は行列全てにかけてあるので、正規方程式の行列プログラムでは省略した。
東西方向の条坊データは(9)〜(12)式で行列の係数に足し込み、南北方向は(13)〜(16)式で足し込む。
これは、東西方向、南北方向ともに、４変数は同じという考え方。
例えば(9)〜(12)式は以下の行列式を意味する。
\(
\left(\begin{array}{ccc}
Σcosθcosθ & Σcosθsinθ & Σcosθa & ΣcosθK \\
Σsinθcosθ & Σsinθsinθ & Σsinθa & ΣsinθK \\
Σacosθ & Σasinθ & Σaa & ΣaK \\
ΣKcosθ & ΣKsinθ & ΣKa & ΣKK
\end{array}\right)
\left(\begin{array}{cc}
Δx \\
Δy \\
Δθ \\
ΔD
\end{array}\right)
=\left(\begin{array}{cc}
Σcosθ Vx0(x0,y0,θ0,D0) \\
ΣsinθVx0(x0,y0,θ0,D0) \\
ΣaVx0(x0,y0,θ0,D0) \\
ΣKVx0(x0,y0,θ0,D0)
\end{array}\right)
\)
これは4元1次の方程式なので、この行列式を解くことでΔx,Δy,Δθ,ΔDが求まる。
以上の内容で作ったプログラムと出力結果が以下。
Cプログラム(テキスト文)
プログラム出力文(テキスト文)
【計算結果の比較】
【計算結果の考察】
上記論文のオリジナル計算では約19回ループして計算結果に収斂しているが、今回の作成したプログラムではループの3回目でほぼ収斂した。違いはオリジナルにはプログラムが添付されていないので比較はできない。1980年代のPCライブラリソフトの精度によるものかもしれない。しかし、本体の最終結果はほぼ同じ結果である。標準偏差も総合は同じであるが、個別の偏差値がプログラムでは２割ほど大き目にでている。これは計算方法の違いと思われる。添付のプログラムの標準偏差の個別計算は長谷川一郎「天体起動論」(1983)p.224-225の方式で計算している。(求める項目の右辺を1とし他の項目の右辺を0として解Aをもとめる。全体の標準偏差をσ、行列の解をAとすると、個別項目の標準偏差は\( σ \sqrt{A} \)で計算できるとする。）
