Module: GeoCalc::Calc::Intersection

Defined in:
lib/geo_calc/calc/intersection.rb

Class Method Summary collapse

Instance Method Summary collapse

Class Method Details

.intersection(p1, brng1, p2, brng2) ⇒ Object

Returns the point of intersection of two paths defined by point and bearing

see http:#williams.best.vwh.net/avform.htm#Intersection

Parameters:

  • p1: (LatLon) —

    First point

  • brng1: (Number) —

    Initial bearing from first point

  • p2: (LatLon) —

    Second point

  • brng2: (Number) —

    Initial bearing from second point



17
18
19
20
21
22
23
24
25
26
27
28
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
 
# File 'lib/geo_calc/calc/intersection.rb', line 17

def self.intersection p1, brng1, p2, brng2
  lat1 = p1.lat.to_rad
  lon1 = p1.lon.to_rad

  lat2 = p2.lat.to_rad
  lon2 = p2.lon.to_rad

  brng13 = brng1.to_rad
  brng23 = brng2.to_rad

  dlat = lat2-lat1
  dlon = lon2-lon1;

  dist12 = 2*Math.asin( Math.sqrt( Math.sin(dlat/2)*Math.sin(dlat/2) + Math.cos(lat1)*Math.cos(lat2)*Math.sin(dlon/2)*Math.sin(dlon/2) ) )
  return nil if dist12 == 0

  # initial/final bearings between points
  brng_a = begin
    Math.acos( ( Math.sin(lat2) - Math.sin(lat1)*Math.cos(dist12) ) / ( Math.sin(dist12)*Math.cos(lat1) ) )
  rescue # protect against rounding
    0
  end

  brng_b = Math.acos( ( Math.sin(lat1) - Math.sin(lat2)*Math.cos(dist12) ) / ( Math.sin(dist12)*Math.cos(lat2) ) )

  brng12, brng21 = if Math.sin(lon2-lon1) > 0
    [brng_a, 2*Math::PI - brng_b] 
  else
    [2*Math::PI - brng_a, brng_b]
  end

  alpha1 = (brng13 - brng12 + Math::PI) % (2*Math::PI) - Math::PI  # angle 2-1-3
  alpha2 = (brng21 - brng23 + Math::PI) % (2*Math::PI) - Math::PI  # angle 1-2-3

  return nil if (Math.sin(alpha1)==0 && Math.sin(alpha2)==0) # infinite intersections
  return nil if (Math.sin(alpha1)*Math.sin(alpha2) < 0)    # ambiguous intersection

  # alpha1 = Math.abs(alpha1);
  # alpha2 = Math.abs(alpha2);
  # ... Ed Williams takes abs of alpha1/alpha2, but seems to break calculation?

  alpha3 = Math.acos( -Math.cos(alpha1)*Math.cos(alpha2) + Math.sin(alpha1)*Math.sin(alpha2)*Math.cos(dist12) )

  dist13 = Math.atan2( Math.sin(dist12)*Math.sin(alpha1)*Math.sin(alpha2), Math.cos(alpha2)+Math.cos(alpha1)*Math.cos(alpha3) )

  lat3 = Math.asin( Math.sin(lat1)*Math.cos(dist13) + Math.cos(lat1)*Math.sin(dist13)*Math.cos(brng13) )

  dlon13 = Math.atan2( Math.sin(brng13)*Math.sin(dist13)*Math.cos(lat1), Math.cos(dist13)-Math.sin(lat1)*Math.sin(lat3) )

  lon3 = lon1 + dlon13;
  lon3 = (lon3 + Math::PI) % (2*Math::PI) - Math::PI  # normalise to -180..180ยบ

  GeoPoint.new lat3.to_deg, lon3.to_deg
end

Instance Method Details

#intersection(brng1, p2, brng2) ⇒ Object 



3
4
5
 
# File 'lib/geo_calc/calc/intersection.rb', line 3

def intersection brng1, p2, brng2
  GeoCalc::Calc::Intersection.intersection self, brng1, p2, brng2
end