2 people like it.

Great Circle Distance (II)

Calculate the great-circle distance between two points on a sphere. (Not done to one-up Samual Bosch's version, but coincidentally inspired by his previous DBSCAN post.)

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 
11: 
12: 
13: 
14: 
15: 
16: 
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: 
71: 
72: 
73: 
74: 
75: 
76: 
77: 
78: 
79: 
80: 
81: 
82: 
83: 
84: 
85: 
86: 
87: 
88: 
89: 
90: 
91: 
92: 
93: 
94: 

namespace Distance

[<AutoOpen>]
module Units =

    open System

    [<Measure>] type km
    [<Measure>] type rad
    [<Measure>] type deg

    let degToRad (degrees : float<deg>) =
        degrees * Math.PI / 180.<deg/rad>

[<AutoOpen>]
module Constants =

    let earthRadius = 6371.<km>
    let marsRadius = 3397.<km>

module GreatCircle = 

    open System

    /// Calculates the great-circle distance between two Latitude/Longitude positions on a sphere of given radius.
    let DistanceBetween (radius:float<km>) lat1 long1 lat2 long2 =
        let lat1r, lat2r, long1r, long2r = lat1 |> degToRad, 
                                           lat2 |> degToRad,
                                           long1 |> degToRad,
                                           long2 |> degToRad
        let deltaLat = lat2r - lat1r
        let deltaLong = long2r - long1r

        let a = Math.Sin(deltaLat/2.<rad>) ** 2. +
                (Math.Sin(deltaLong/2.<rad>) ** 2. * Math.Cos((double)lat1r) * Math.Cos((double)lat2r))

        let c = 2. * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1.-a))

        radius * c

    /// Calculate DistanceBetween for Earth.
    let DistanceBetweenEarth = DistanceBetween earthRadius

    /// Calculate DistanceBetween for Mars.
    let DistanceBetweenMars = DistanceBetween marsRadius

module GreatCircleTests =

    open NUnit.Framework
    open FsUnit

    [<TestFixture>]
    type ``Given the DistanceBetween function for Earth``() =

        // Error margin for non-sphericality of Earth:
        let ErrorMargin = 0.003; // 0.3%

        // Travel no distance:
        [<TestCase(0., 0., 0., 0., 0.)>]
        // Travel along the equator eastwards for 90 degrees:
        [<TestCase(0., 0., 0., 90., 10018.79)>]
        // Travel along the equator westwards for 90 degrees:
        [<TestCase(0., 0., 0., -90., 10018.79)>]
        // Travel along the equator eastwards for 180 degrees:
        [<TestCase(0., 0., 0., 180., 20037.58)>]
        // Travel along the equator westwards for 180 degrees:
        [<TestCase(0., 0., 0., -180., 20037.58)>]
        // Travel along the meridian northwards 90 degrees:
        [<TestCase(0., 0., 90., 0., 10018.79)>]
        // Travel along the meridian soutwards 90 degrees:
        [<TestCase(0., 0., -90., 0., 10018.79)>]
        // Travel from Farnham to Reigate:
        [<TestCase(51.214, -0.799, 51.230, -0.188, 42.5)>]
        // Travel from London to Sidney Australia:
        [<TestCase(51.51, -0.13, -33.86, 151.21, 16998.)>]
        
        member t.``the function returns the right result``(lat1, long1, lat2, long2, expected:float<km>) =
            let actual = GreatCircle.DistanceBetweenEarth lat1 long1 lat2 long2
            let error = expected * ErrorMargin
            actual |> should (equalWithin error) expected

    [<TestFixture>]
    type ``Given the DistanceBetween function for Mars``() =

        // Error margin for non-sphericality of Mars:
        let ErrorMargin = 0.003; // 0.3%

        // Travel from Olympus Mons to Pavonis Mons:
        [<TestCase(18.65, 226.2, 1.48, 247.04, 1582.)>]
        
        member t.``the function returns the right result``(lat1, long1, lat2, long2, expected:float<km>) =
            let actual = GreatCircle.DistanceBetweenMars lat1 long1 lat2 long2
            let error = expected * ErrorMargin
            actual |> should (equalWithin error) expected
Multiple items
type AutoOpenAttribute =
  inherit Attribute
  new : unit -> AutoOpenAttribute
  new : path:string -> AutoOpenAttribute
  member Path : string

Full name: Microsoft.FSharp.Core.AutoOpenAttribute

--------------------
new : unit -> AutoOpenAttribute
new : path:string -> AutoOpenAttribute
namespace System
Multiple items
type MeasureAttribute =
  inherit Attribute
  new : unit -> MeasureAttribute

Full name: Microsoft.FSharp.Core.MeasureAttribute

--------------------
new : unit -> MeasureAttribute
[<Measure>]
type km

Full name: Distance.Units.km
[<Measure>]
type rad

Full name: Distance.Units.rad
[<Measure>]
type deg

Full name: Distance.Units.deg
val degToRad : degrees:float<deg> -> float<rad>

Full name: Distance.Units.degToRad
val degrees : float<deg>
Multiple items
val float : value:'T -> float (requires member op_Explicit)

Full name: Microsoft.FSharp.Core.Operators.float

--------------------
type float = Double

Full name: Microsoft.FSharp.Core.float

--------------------
type float<'Measure> = float

Full name: Microsoft.FSharp.Core.float<_>
type Math =
  static val PI : float
  static val E : float
  static member Abs : value:sbyte -> sbyte + 6 overloads
  static member Acos : d:float -> float
  static member Asin : d:float -> float
  static member Atan : d:float -> float
  static member Atan2 : y:float * x:float -> float
  static member BigMul : a:int * b:int -> int64
  static member Ceiling : d:decimal -> decimal + 1 overload
  static member Cos : d:float -> float
  ...

Full name: System.Math
field Math.PI = 3.14159265359
val earthRadius : float<km>

Full name: Distance.Constants.earthRadius
val marsRadius : float<km>

Full name: Distance.Constants.marsRadius
val DistanceBetween : radius:float<km> -> lat1:float<deg> -> long1:float<deg> -> lat2:float<deg> -> long2:float<deg> -> float<km>

Full name: Distance.GreatCircle.DistanceBetween


 Calculates the great-circle distance between two Latitude/Longitude positions on a sphere of given radius.
val radius : float<km>
val lat1 : float<deg>
val long1 : float<deg>
val lat2 : float<deg>
val long2 : float<deg>
val lat1r : float<rad>
val lat2r : float<rad>
val long1r : float<rad>
val long2r : float<rad>
val deltaLat : float<rad>
val deltaLong : float<rad>
val a : float
Math.Sin(a: float) : float
Math.Cos(d: float) : float
Multiple items
val double : value:'T -> float (requires member op_Explicit)

Full name: Microsoft.FSharp.Core.ExtraTopLevelOperators.double

--------------------
type double = Double

Full name: Microsoft.FSharp.Core.double
val c : float
Math.Atan2(y: float, x: float) : float
Math.Sqrt(d: float) : float
val DistanceBetweenEarth : (float<deg> -> float<deg> -> float<deg> -> float<deg> -> float<km>)

Full name: Distance.GreatCircle.DistanceBetweenEarth


 Calculate DistanceBetween for Earth.
val DistanceBetweenMars : (float<deg> -> float<deg> -> float<deg> -> float<deg> -> float<km>)

Full name: Distance.GreatCircle.DistanceBetweenMars


 Calculate DistanceBetween for Mars.
module GreatCircleTests

from Distance
namespace NUnit
namespace NUnit.Framework
namespace FsUnit
Multiple items
type TestFixtureAttribute =
  inherit Attribute
  new : unit -> TestFixtureAttribute + 1 overload
  member Arguments : obj[]
  member Categories : IList
  member Category : string with get, set
  member Description : string with get, set
  member Ignore : bool with get, set
  member IgnoreReason : string with get, set
  member TypeArgs : Type[] with get, set

Full name: NUnit.Framework.TestFixtureAttribute

--------------------
TestFixtureAttribute() : unit
TestFixtureAttribute([<System.ParamArray>] arguments: obj []) : unit
val ErrorMargin : float
Multiple items
type TestCaseAttribute =
  inherit Attribute
  new : [<ParamArray>] arguments:obj[] -> TestCaseAttribute + 3 overloads
  member Arguments : obj[]
  member Categories : IList
  member Category : string with get, set
  member Description : string with get, set
  member ExpectedException : Type with get, set
  member ExpectedExceptionName : string with get, set
  member ExpectedMessage : string with get, set
  member ExpectedResult : obj with get, set
  member Explicit : bool with get, set
  ...

Full name: NUnit.Framework.TestCaseAttribute

--------------------
TestCaseAttribute([<System.ParamArray>] arguments: obj []) : unit
TestCaseAttribute(arg: obj) : unit
TestCaseAttribute(arg1: obj, arg2: obj) : unit
TestCaseAttribute(arg1: obj, arg2: obj, arg3: obj) : unit
val t : Given the DistanceBetween function for Earth
member Given the DistanceBetween function for Earth.( the function returns the right result ) : lat1:float<deg> * long1:float<deg> * lat2:float<deg> * long2:float<deg> * expected:float<km> -> unit

Full name: Distance.GreatCircleTests.Given the DistanceBetween function for Earth.( the function returns the right result )
val expected : float<km>
Multiple items
val float : value:'T -> float (requires member op_Explicit)

Full name: Microsoft.FSharp.Core.Operators.float

--------------------
type float = System.Double

Full name: Microsoft.FSharp.Core.float

--------------------
type float<'Measure> = float

Full name: Microsoft.FSharp.Core.float<_>
val actual : float<km>
module GreatCircle

from Distance
val error : float<km>
val should : f:('a -> #Constraints.Constraint) -> x:'a -> y:obj -> unit

Full name: FsUnit.TopLevelOperators.should
val equalWithin : tolerance:'a -> x:'b -> Constraints.EqualConstraint

Full name: FsUnit.TopLevelOperators.equalWithin
val t : Given the DistanceBetween function for Mars
member Given the DistanceBetween function for Mars.( the function returns the right result ) : lat1:float<deg> * long1:float<deg> * lat2:float<deg> * long2:float<deg> * expected:float<km> -> unit

Full name: Distance.GreatCircleTests.Given the DistanceBetween function for Mars.( the function returns the right result )
Raw view Test code New version

More information

Link:http://fssnip.net/jA
Posted:11 years ago
Author:Kit Eason
Tags: circle distance , sphere , nunit , fsunit

This web site is created using F# and Suave web server. It is hosted on Azure and the source code is on GitHub. Contributions are welcome!

The first version of fssnip.net has been created by @tomaspetricek back in 2010. This web site is a new, open-source and contribution-friendly version.