ArduinoとGPS（正確にはGNSS)受信モジュール、OLED（有機ELディスプレイ）を使って目的地までの直線距離を表示するアイテムを作ってみました。ところが 目的地に近づくと表示がおかしくなる 事象が発生。当初はきっと原因は桁数まるめの影響ですぐに解決するだろうとタカをくくっていたのですが、調べれば調べるほどに謎が深まっています。どうにも行き詰まりつつあるので詳しい方のお知恵をお借りしたく記事にすることにしました。

テスト中の様子。シャッタースピードの関係で上部の文字が消えてしまっていますが、OLEDの最上段が目的地までの直線距離、中段が現在地の標高、下段が受信している衛星数のバーグラフ表示です。

＜事象＞

１）目的地の5kmくらい手前まで近づくと距離表示がおかしくなる。

　　（具体的には3.114km、2.202km、0.000kmをランダムに表示するようになる。）

２）目的地から大きく離れているときには正確な距離が表示される。

３）GPSから受信した位置情報（緯度経度）を元に手計算すると、目的地付近でもちゃんと距離が出る。

４）GPSから受信した位置情報のバラツキは大きくても20m程度。分解能は1.5m程度。

５）GPSからの受信の代わりに位置情報を手入力した値を用いた場合も、現在地と目的地が近いと表示がおかしくなる。

上記３）、４）、５）から、 GPSの位置精度のバラツキが原因ではない ことが分かりました。

手入力した固定値の位置情報でも表示がおかしい ことから、計算過程に何らかの問題があると思われます。計算式は下記のものを用いています。

直線距離 ＝ 地球の半径 × Arccos(sin(y1) × sin(y2) + cos(y1) × cos(y2) × cos(x2-x1))

地球の半径：6378.137km

x1：現在地の経度

y1：現在地の緯度

x2：目的地の緯度

y2：目的地の経度

計算式や地球の半径の値はネットから引用しています。Google Mapで適当な2地点の緯度経度を調べて上記の式で計算してみると正確に距離が算出されるので、計算式は問題ないと思います・・・少なくとも手計算では。

ところがこれをArduinoに書き込んで計算させてみると、直線距離が近いときに手計算とアンマッチが生じてしまうのです。

式のどの部分でアンマッチが生じているかを検証したのが下のグラフです。

横軸は式の中の各項、縦軸はArduinoと手計算との差異のパーセンテージです。直線距離0kmの場合から710kmの場合まで13パターンをプロットしています。これを見るとArccosのカッコの中までは差異が±0.00001%程度に収まっており問題ないように思えるのですが・・・

その先、 Arccosを計算すると途端にブッ飛びます 。−100%になっているのは、Arduinoでは計算結果がゼロになってしまっているものです。

※横軸の”acos()”はArccos(sin(y1) × sin(y2) + cos(y1) × cos(y2) × cos(x2-x1))を、”kyori”は目的地までの直線距離を表しています。

ちなみにArduinoと手計算での直線距離の算出結果は以下の通りです。目的地までの直線距離が28kmでようやく差異50mとなります。せめてこのくらいの誤差で収まってほしいのですが・・・。目的地の5km以上手前からあてにならないのでは役に立ちません。まぁそもそも「直線距離を表示」って時点で実用性などほとんどなく、スマホやクルマのナビの方が便利なのは分かりきっていますが（笑。

＜疑問点＞

・Arccosの段階で途端に誤差が大きくなるのは何故か？

・手入力した固定値をそのままserial.printで表示させるとなぜか少し値がズレるのは何故か？

　たとえば135.001000，35.001000と入力　→　表示は135.001007080,35.000999451となる。

　ただし整数で入力したものはズレない（例：135.000000）。

以下に検証用のスケッチを記します。

これのx2,y2の値を色々と変えて各項の値や直線距離の計算値を吐き出させたのが上記グラフです。

疑問点についての解答やヒント、ここ確認してみたら？などアドバイスいただけると嬉しいです。

#include <stdio.h>

#include <math.h>

double EARTH_RAD = 6378.137;

void setup() {

 Serial.begin(57600);

 while (!Serial) {

 ;

 }

} 

void loop() {

    double x1 = 135.000000 * (M_PI/180);  //現在地の経度

    double y1 = 35.000000 * (M_PI/180);   //現在地の緯度

    double x2 = 135.100000 * (M_PI/180);  //目的地の経度

    double y2 = 35.100000 * (M_PI/180);   //目的地の緯度

    double kyori = EARTH_RAD * acos(sin(y1) * sin(y2) + cos(y1) * cos(y2) * cos(x2 - x1));  //目的地までの直線距離

    

    Serial.print(kyori,3);

    }