Delphi で２地点間の距離と方位角を求める

# はじめに
[keisan.casio.jp](https://keisan.casio.jp/) に２地点間の距離と方位角を求める記事がありまして、
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/21785/f29e6811-3539-aa09-111c-e678e3d134d3.png)
これを 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.
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/21785/1df07a6c-adef-db08-25fc-6e07574f02ed.png)
**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 の算術関数](https://qiita.com/ht_deko/items/15b27e00fac9d723d9f0#%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.
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/21785/cbcb62d0-a857-d593-df7c-99c4206e7cf5.png)
ちょっと長くなりましたね。コードの 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)なのですが、
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/21785/621135ac-a9f7-7b26-3425-6a71f3cd967d.png)
方位角は、**世界地図を見た時の直感とは異なる方向を指しています**よね。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/21785/062bcb27-f61e-faab-f831-a6c639063277.png)


[^1]: [System.Math](http://docwiki.embarcadero.com/Libraries/ja/System.Math) には他にも見慣れない関数が...いつの間に？
[^2]: `ArcTan2()` がアルゴリズム的にこれでいいのかは要検証。Delphi の `System.Math` にある `ArcTan2()` はルックアップテーブルを使って値を求めている。恐らくこちらの方が高速に動作するアルゴリズムなのだろう。