# Delphi で２地点間の距離と方位角を求める --- tags: Delphi programming Pascal objectpascal created_at: 2020-12-19 updated_at: 2020-12-20 --- # はじめに [keisan.casio.jp](https://keisan.casio.jp/) に２地点間の距離と方位角を求める記事がありまして、  これを Delphi で書いてみようという趣旨となります。 - [２地点間の距離と方位角 (keisan.casio.jp)](https://keisan.casio.jp/exec/system/1257670779) なお、地球を楕円体として正確に計算するアルゴリズムは国土地理院のサイトに掲載されていますので、興味のある方はそちらも参照してみてください。 **See also:** - [距離と方位角の計算 (国土地理院)](https://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/bl2stf.html) # コード ## Delphi の持つ数値計算ルーチン 数値計算ルーチンは割と豊富に揃っていますので、アルゴリズムの記述に困る事はないと思います。 - [数値計算ルーチン (DocWiki)](http://docwiki.embarcadero.com/RADStudio/ja/%E6%95%B0%E5%80%A4%E8%A8%88%E7%AE%97%E3%83%AB%E3%83%BC%E3%83%81%E3%83%B3) ## Delphi で書いてみたコード 先のサイトにある計算式をコードにしてみるとこんな感じでしょうか？ ```pascal:ll2da.dpr program ll2da; {$APPTYPE CONSOLE} uses System.Math; const r = 6378.137; var x1, y1, x2, y2: Double; begin Write('Longitude [A]: '); Readln(x1); x1 := DegToRad(x1); Write('Latitude [A]: '); Readln(y1); y1 := DegToRad(y1); Write('Longitude [B]: '); Readln(x2); x2 := DegToRad(x2); Write('Latitude [B]: '); Readln(y2); y2 := DegToRad(y2); Writeln; var dx := x2 - x1; var dst := r * ArcCos(Sin(y1) * Sin(y2) + Cos(y1) * Cos(y2) * Cos(dx)); var phi := FMod(360 + 90 - RadToDeg(ArcTan2(Cos(y1) * Tan(y2) - Sin(y1) * Cos(dx), Sin(dx))), 360); Writeln('Distance:', dst:10:2); Writeln('Azimuth :', phi:10:2); end. ```  **XE8 以降には実数の剰余を求める事ができる FMod() なんていう関数があります [^1]。**`FMod()` を使わず、従来の数値計算ルーチンだけで記述すると次のようになるかと思います。 ```pascal var phi := 360 + 90 - RadToDeg(ArcTan2(Cos(y1) * Tan(y2) - Sin(y1) * Cos(dx), Sin(dx))); phi := Trunc(phi) mod 360 + Frac(phi); ``` 剰余を求める部分は次のようにも書けるのですが、 ```pascal phi := phi - Trunc(phi / 360) * 360; ``` 実数除算する事に (なんとなく) 抵抗があり、先のコードのように整数部と小数部を分けて計算する方が個人的には好きです。 **See also:** - [System.Math.FMod (DocWiki)](http://docwiki.embarcadero.com/Libraries/ja/System.Math.FMod) ## 標準 Pascal で書いてみたコード [標準 Pascal の算術関数](./15b27e00fac9d723d9f0.md#%E7%AE%97%E8%A1%93%E9%96%A2%E6%95%B0-arithmetic-functions)には `Tan()` も `ArcCos()` もないけど、Pascal 50 周年だから書いてみるしかないね！(？) ```pascal:ll2da.pas program ll2da(Input, Output); const Pi = 3.1415926535; r = 6378.137; var x1, y1, x2, y2, dx, dst, phi: Real; function ArcTan2(Y, X: Real): Real; Forward; function ArcCos(X: Real): Real; begin ArcCos := ArcTan2(Sqrt((1 + X) * (1 - X)), X); end; { ArcCos } function ArcTan2; var v: Real; begin if X = 0 then begin if Y > 0 then ArcTan2 := Pi / 2 else if Y < 0 then ArcTan2 := -Pi / 2 else ArcTan2 := 0; end else begin v := ArcTan(Y / X); if X > 0 then ArcTan2 := v else begin if Y < 0 then ArcTan2 := v - Pi else ArcTan2 := v + Pi; end; end; end; { ArcTan2 } function DegToRad(Degrees: Real): Real; begin DegToRad := Degrees * Pi / 180; end; { DegToRad } function Frac(X: Real): Real; begin Frac := X - Trunc(X); end; { Frac } function RadToDeg(Radians: Real): Real; begin RadToDeg := Radians * 180 / Pi; end; { RadToDeg } function Tan(X: Real): Real; begin Tan := Sin(X) / Cos(X); end; { Tan } begin Write('Longitude [A]: '); Readln(x1); x1 := DegToRad(x1); Write('Latitude [A]: '); Readln(y1); y1 := DegToRad(y1); Write('Longitude [B]: '); Readln(x2); x2 := DegToRad(x2); Write('Latitude [B]: '); Readln(y2); y2 := DegToRad(y2); Writeln; dx := x2 - x1; dst := r * ArcCos(Sin(y1) * Sin(y2) + Cos(y1) * Cos(y2) * Cos(dx)); phi := 360 + 90 - RadToDeg(ArcTan2(Cos(y1) * Tan(y2) - Sin(y1) * Cos(dx), Sin(dx))); phi := Trunc(phi) mod 360 + Frac(phi); Writeln('Distance:', dst:10:2); Writeln('Azimuth :', phi:10:2); end. ```  ちょっと長くなりましたね。コードの 4 割位を `ArcTan2()` [^2] 関数が占めています (w # おわりに 参考にしたサイトの入力ボックスにデフォルトで設定されている２地点は[東京とメッカ](https://www.google.co.jp/maps/dir/%E6%9D%B1%E4%BA%AC%E3%80%81%E6%9D%B1%E4%BA%AC%E9%83%BD/%E3%83%A1%E3%83%83%E3%82%AB,+%E3%82%B5%E3%82%A6%E3%82%B8%E3%82%A2%E3%83%A9%E3%83%93%E3%82%A2/@27.3613539,53.7129842,3z/data=!4m14!4m13!1m5!1m1!1s0x60188b857628235d:0xcdd8aef709a2b520!2m2!1d139.7690174!2d35.6803997!1m5!1m1!1s0x15c21b4ced818775:0x98ab2469cf70c9ce!2m2!1d39.8579118!2d21.3890824!3e4?hl=ja)なのですが、  方位角は、**世界地図を見た時の直感とは異なる方向を指しています**よね。  [^1]: [System.Math](http://docwiki.embarcadero.com/Libraries/ja/System.Math) には他にも見慣れない関数が...いつの間に？ [^2]: `ArcTan2()` がアルゴリズム的にこれでいいのかは要検証。Delphi の `System.Math` にある `ArcTan2()` はルックアップテーブルを使って値を求めている。恐らくこちらの方が高速に動作するアルゴリズムなのだろう。