緯度経度から方位角を算出する式

Question

###前提・実現したいこと
以下のサイトの計算式を使用して緯度経度から方位角と距離を計算するプログラムをpythonを用いて作成を行いました。
CASIO2地点間の距離

python
import math
from math import sin
from math import cos
from math import tan
from math import atan2
from math import pi
from math import acos

radeg = 180/pi #radからdegに変換
r= 6378.137*10**3
x1=139.988909
y1=35.685828
x2=139.990339
y2=35.685879
deltax = x2-x1
ans= 90-atan2(sin(deltax), cos(y1)*tan(y2)-sin(y1)*cos(deltax) )*radeg #方位角計算
distance = r*acos(sin(y1)*sin(y2)+cos(y1)*cos(y2)*cos(deltax)) #距離計算

#phi= 90-math.atan2()
print('方位角:'+str(ans))
print('距離:'+str(distance/1000)+'Km')

インターネットのサイトの計算結果と一致しないため、式が間違っているのかそれとも私の実装方法が間違っているのかが知りたいです。
 また、他の言語での実装も考えてるため、申し訳ありませんがpyprojなど使用しないで純粋に数式だけを用いて実装を行うことを前提としています。
もし、CASIOの数式以外に緯度と経度から方位角と距離を求めることができる数式をご存知であればそちらでも構わないので教えてください

よろしくお願い申し上げます。

Answer 1

コードをぱぱっと置いときます。

結論から言うと、該当サイトの計算式はあっていますが、北を基準にしているところと、式にずれがあるようです。
その為、atan2の負の符号に関して条件分岐が足りないようでした。

from math import sin, cos, tan, atan2, acos, pi

def azimuth(x1, y1, x2, y2):
    # Radian角に修正
    _x1, _y1, _x2, _y2 = x1*pi/180, y1*pi/180, x2*pi/180, y2*pi/180
    Δx = _x2 - _x1
    _y = sin(Δx)
    _x = cos(_y1) * tan(_y2) - sin(_y1) * cos(Δx)
    
    psi = atan2(_y, _x) * 180 / pi
    if psi < 0:
        return 360 + atan2(_y, _x) * 180 / pi
    else:
        return atan2(_y, _x) * 180 / pi

def distance(x1, y1, x2, y2, r):
    _x1, _y1, _x2, _y2 = x1*pi/180, y1*pi/180, x2*pi/180, y2*pi/180
    Δx = _x2 - _x1
    val = sin(_y1) * sin(_y2) + cos(_y1) * cos(_y2) * cos(Δx)
    return r * acos(val)

x1 = 139.988909
y1 = 35.685828
x2 = 139.990339
y2 = 35.685879
r  = 6378.137e3

angle = azimuth(x1, y1, x2, y2)
dis   = distance(x1, y1, x2, y2, r) / 1e3 # kmに変換
print("方位角 : {0:.3f} 度".format(angle))
print("距離 : {0:.3f} km".format(dis))

# 結果
# 方位角 : 87.485 度
# 距離 : 0.129 km

これで CASIO2地点間の距離 の計算結果とはずです。

Answer 2

式が間違っていると思います。主に度とラジアンとか。

import math
from math import sin
from math import cos
from math import tan
from math import atan2
from math import acos
from math import radians
from math import degrees

r= 6378.137 #単位はkmのままでOK

# test data 1
# x1=139.74477
# y1=35.6544
# x2=39.8261
# y2=21.4225

# test data 2
x1=139.988909
y1=35.685828
x2=139.990339
y2=35.685879

# 経度・緯度を度からラジアンに変換 sin,cos,tan,etc はラジアンを使う
x1=radians(x1)
y1=radians(y1)
x2=radians(x2)
y2=radians(y2)

deltax = x2 - x1

ans = degrees(atan2(sin(deltax),(cos(y1)*tan(y2)-sin(y1)*cos(deltax))))%360 
# %300   : マイナスの場合は正の角度にする
# 90 -   : 要らないと思う なんであったんだろう？？よくわからない
# degrees :  atan2の返り値はラジアンなので度になおす

distance = r*acos(sin(y1)*sin(y2)+cos(y1)*cos(y2)*cos(deltax))

#地球の大きさとアークコサインで距離を出すから、単位は基準になった地球の大きさrに使ってた単位kmになる

print('方位角:'+str(ans))
print('距離:'+str(distance)+'km')

Comments

[takuto765]

oskbtさん
返信が遅くなり申し訳ございません。
夜分遅くにかかわらず素早いご回答ありがとうございました。
oskbtさんの回答のコードはシンプルでありベストアンサーで良かったのですが
今回他のシステムでの実装も考えている中で%360による負の値の直し方を
調べているときに負の数の剰余は言語によって答えが異なるということを
知りました。
そのためGazelleさんの答えをベストアンサーにさせていただきました。
しかしながら、%360でマイナスの場合は正の角度にするという考え方を
知ることができたため、大変勉強になりました！

今回は本当にご回答ありがとうございました！