29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
|
# File 'lib/terraformer/geodesic.rb', line 29
def self.compute_distance_and_bearing lat_1, lon_1, lat_2, lon_2
lat_1 = lat_1.to_rad
lat_2 = lat_2.to_rad
lon_1 = lon_1.to_rad
lon_2 = lon_2.to_rad
a = 6378137.to_d
b = 6356752.3142.to_d
f = (a - b) / a
a_sq_minus_b_sq_over_b_sq = (a * a - b * b) / (b * b)
_l = lon_2 - lon_1
_a = BigMath::ZERO
_u1 = atan((BigMath::ONE - f) * tan(lat_1))
_u2 = atan((BigMath::ONE - f) * tan(lat_2))
cos_u1 = cos(_u1)
cos_u2 = cos(_u2)
sin_u1 = sin(_u1)
sin_u2 = sin(_u2)
cos_u1_cos_u2 = cos_u1 * cos_u2
sin_u1_sin_u2 = sin_u1 * sin_u2
sigma = BigMath::ZERO
delta_sigma = BigMath::ZERO
cos_sq_alpha = BigMath::ZERO
cos2_s_m = BigMath::ZERO
cos_sigma = BigMath::ZERO
sin_sigma = BigMath::ZERO
cos_lambda = BigMath::ZERO
sin_lambda = BigMath::ZERO
_lambda = _l
MAX_ITERS.times do |n|
_lambda_orig = _lambda
cos_lambda = cos(_lambda)
sin_lambda = sin(_lambda)
t1 = cos_u2 * sin_lambda
t2 = cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos_lambda
sin_sq_sigma = t1 * t1 + t2 * t2
sin_sigma = sqrt(sin_sq_sigma)
cos_sigma = sin_u1_sin_u2 + cos_u1_cos_u2 * cos_lambda
sigma = atan2(sin_sigma, cos_sigma)
sin_alpha = (sin_sigma == BigMath::ZERO) ? BigMath::ZERO : (cos_u1_cos_u2 * sin_lambda / sin_sigma)
cos_sq_alpha = BigMath::ONE - sin_alpha * sin_alpha
cos_2_s_m = (cos_sq_alpha == BigMath::ZERO) ? BigMath::ZERO : (cos_sigma - BigMath::TWO * sin_u1_sin_u2 / cos_sq_alpha)
u_squared = cos_sq_alpha * a_sq_minus_b_sq_over_b_sq
_a = BigMath::ONE + (u_squared / 16384.0) *
(4096.0 + u_squared *
(-768.0 + u_squared * (320.0 - 175.0 * u_squared)))
_b = (u_squared / 1024.0) *
(256.0 + u_squared *
(-128.0 + u_squared * (74.0 - 47.0 * u_squared)))
_c = (f / 16.0) *
cos_sq_alpha *
(4.0 + f * (4.0 - 3.0 * cos_sq_alpha))
cos2_s_m_sq = cos2_s_m * cos2_s_m
delta_sigma = _b * sin_sigma *
(cos2_s_m + (_b / 4.0) *
(cos_sigma * (BigMath::N_ONE + BigMath::TWO * cos2_s_m_sq) -
(_b / 6.0) * cos2_s_m *
(-3.0 + 4.0 * sin_sigma * sin_sigma) *
(-3.0 + 4.0 * cos2_s_m_sq)))
_lambda = _l +
(BigMath::ONE - _c) * f * sin_alpha *
(sigma + _c * sin_sigma *
(cos2_s_m + _c * cos_sigma *
(BigMath::N_ONE + BigMath::TWO * cos2_s_m * cos2_s_m)))
delta = (_lambda - _lambda_orig) / _lambda
if delta.abs < 1.0e-12
break
end
end
{ distance: (b * _a * (sigma - delta_sigma)),
bearing: {
initial: atan2(cos_u2 * sin_lambda, cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos_lambda).to_deg,
final: atan2(cos_u1 * sin_lambda, -(sin_u1) * cos_u2 + cos_u1 * sin_u2 * cos_lambda).to_deg
}
}
end
|