summaryrefslogtreecommitdiff
path: root/vendor/github.com/StefanSchroeder
diff options
context:
space:
mode:
Diffstat (limited to 'vendor/github.com/StefanSchroeder')
-rwxr-xr-xvendor/github.com/StefanSchroeder/Golang-Ellipsoid/LICENSE7
-rwxr-xr-xvendor/github.com/StefanSchroeder/Golang-Ellipsoid/ellipsoid/ellipsoid.go713
2 files changed, 720 insertions, 0 deletions
diff --git a/vendor/github.com/StefanSchroeder/Golang-Ellipsoid/LICENSE b/vendor/github.com/StefanSchroeder/Golang-Ellipsoid/LICENSE
new file mode 100755
index 0000000..f7f56d0
--- /dev/null
+++ b/vendor/github.com/StefanSchroeder/Golang-Ellipsoid/LICENSE
@@ -0,0 +1,7 @@
+Copyright (c) 2011 Stefan Schroeder
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS ORIMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS INTHE SOFTWARE.
diff --git a/vendor/github.com/StefanSchroeder/Golang-Ellipsoid/ellipsoid/ellipsoid.go b/vendor/github.com/StefanSchroeder/Golang-Ellipsoid/ellipsoid/ellipsoid.go
new file mode 100755
index 0000000..91874ef
--- /dev/null
+++ b/vendor/github.com/StefanSchroeder/Golang-Ellipsoid/ellipsoid/ellipsoid.go
@@ -0,0 +1,713 @@
+package ellipsoid
+
+// Written in Go by Stefan Schroeder, New York, 2013
+// Version 1.0 based on Geo::Ellipsoid Version 1.12.
+// Version 1.1 Added ECEF functions.
+// Version 1.2 Replaced Fabs with Abs.
+// Version 1.3 Added Displacement function
+
+/*
+
+SYNOPSIS
+
+See hello-world.go example.
+
+DESCRIPTION
+
+ellipsoid performs geometrical calculations on the surface of
+an ellipsoid. An ellipsoid is a three-dimension object formed from
+the rotation of an ellipse about one of its axes. The approximate
+shape of the earth is an ellipsoid, so ellipsoid can accurately
+calculate distance and bearing between two widely-separated locations
+on the earth's surface.
+
+The shape of an ellipsoid is defined by the lengths of its
+semi-major and semi-minor axes. The shape may also be specifed by
+the flattening ratio f as:
+
+ f = ( semi-major - semi-minor ) / semi-major
+
+which, since f is a small number, is normally given as the reciprocal
+of the flattening 1/f.
+
+The shape of the earth has been surveyed and estimated differently
+at different times over the years. The two most common sets of values
+used to describe the size and shape of the earth in the United States
+are 'NAD27', dating from 1927, and 'WGS84', from 1984. United States
+Geological Survey topographical maps, for example, use one or the
+other of these values, and commonly-available Global Positioning
+System (GPS) units can be set to use one or the other.
+See "DEFINED ELLIPSOIDS" below for the ellipsoid survey values
+that may be selected for use by ellipsoid.
+
+*/
+
+import "math"
+import "fmt"
+
+const (
+ pi = math.Pi
+ twopi = math.Pi * 2.0
+ maxLoopCount = 20
+ eps = 1.0e-23
+ debug = false
+ // Meter is one of the output/input units.
+ Meter = 0 // 1.0 meter
+ // Foot is one of the output/input units.
+ Foot = 1 // 0.3048 meter are a foot
+ // Kilometer is one of the output/input units.
+ Kilometer = 2 // 1000.0 meter are a kilometer
+ // Mile is one of the output/input units.
+ Mile = 3 // 1609.344 meter are a mile
+ // Nm (nautical mile) is one of the output/input units.
+ Nm = 4 // 1852.0 meter are a nautical mile,
+ // Degrees is one of the possible angle units for input/output.
+ Degrees = iota
+ // Radians is one of the possible angle units for input/output.
+ Radians = iota
+ // LongitudeIsSymmetric determines that the output longitude shall be symmetric.
+ LongitudeIsSymmetric = true
+ // LongitudeNotSymmetric determines that the output longitude shall not be symmetric.
+ LongitudeNotSymmetric = false
+ // BearingIsSymmetric determines that the output bearing shall be symmetric.
+ BearingIsSymmetric = true
+ // BearingNotSymmetric determines that the output bearing shall not be symmetric.
+ BearingNotSymmetric = false
+)
+
+// Ellipsoid is the main object to store information about one ellispoid.
+type Ellipsoid struct {
+ Ellipse ellipse
+ Units int
+ DistanceUnits int
+ LongitudeSymmetric bool
+ BearingSymmetry bool
+ DistanceFactor float64
+ // Having the DistanceFactor AND the DistanceUnits in this struct is redundant
+ // but it looks nicer in the code.
+}
+
+type ellipse struct {
+ Equatorial float64
+ InvFlattening float64
+}
+
+// Location is one coordinate in LLA.
+type Location struct {
+ Lat float64
+ Lon float64
+ Ele float64
+}
+
+func deg2rad(d float64) (r float64) {
+ return d * pi / 180.0
+}
+func rad2deg(d float64) (r float64) {
+ return d * 180.0 / pi
+}
+
+/* Init
+
+The Init constructor must be called with a list of parameters to set
+the value of the ellipsoid to be used, the value of the units to be
+used for angles and distances, and whether or not the output range
+of longitudes and bearing angles should be symmetric around zero
+or always greater than zero. There is no default constructor, all
+arguments are required; they may not be abbreviated.
+
+Example:
+
+ geo := ellipsoid.Init(
+ "WGS84", // for possible values see below.
+ ellipsoid.Degrees, // possible values: Degrees or Radians
+ ellipsoid.Meter, // possible values: Meter, Kilometer,
+ // Foot, Nm, Mile
+ ellipsoid.LongitudeIsSymmetric, // possible values
+ // LongitudeIsSymmetric or
+ // LongitudeNotSymmetric
+ ellipsoid.BearingIsSymmetric // possible
+ // values BearingIsSymmetric or
+ // BearingNotSymmetric
+ )
+
+*/
+func Init(name string, units int, distUnits int, longSym bool, bearSym bool) (e Ellipsoid) {
+ m := map[string]ellipse{
+ "AIRY": {6377563.396, 299.3249646},
+ "AIRY-MODIFIED": {6377340.189, 299.3249646},
+ "AUSTRALIAN": {6378160.0, 298.25},
+ "BESSEL-1841": {6377397.155, 299.1528128},
+ "BESSEL-1841-NAMIBIA": {6377483.865, 299.152813},
+ "CLARKE-1866": {6378206.400, 294.978698},
+ "CLARKE-1880": {6378249.145, 293.465},
+ "EVEREST-1830": {6377276.345, 300.8017},
+ "EVEREST-1948": {6377304.063, 300.8017},
+ "EVEREST-SABAH-SARAWAK": {6377298.556, 300.801700},
+ "EVEREST-1956": {6377301.243, 300.801700},
+ "EVEREST-1969": {6377295.664, 300.801700},
+ "FISHER-1960": {6378166.0, 298.3},
+ "FISCHER-1960-MODIFIED": {6378155.000, 298.300000},
+ "FISHER-1968": {6378150.0, 298.3},
+ "GRS80": {6378137.0, 298.25722210088},
+ "HELMERT-1906": {6378200.000, 298.300000},
+ "HOUGH-1956": {6378270.0, 297.0},
+ "HAYFORD": {6378388.0, 297.0},
+ "IAU76": {6378140.0, 298.257},
+ "INTERNATIONAL": {6378388.000, 297.000000},
+ "KRASSOVSKY-1938": {6378245.0, 298.3},
+ "NAD27": {6378206.4, 294.9786982138},
+ "NWL-9D": {6378145.0, 298.25},
+ "SGS85": {6378136.000, 298.257000},
+ "SOUTHAMERICAN-1969": {6378160.0, 298.25},
+ "SOVIET-1985": {6378136.0, 298.257},
+ "WGS60": {6378165.000, 298.300000},
+ "WGS66": {6378145.000, 298.250000},
+ "WGS72": {6378135.0, 298.26},
+ "WGS84": {6378137.0, 298.257223563},
+ }
+
+ e2, ok := m[name]
+ if !ok {
+ fmt.Printf("ellipsoid.go: Warning: Invalid ellipse type '%v'\n", name)
+ }
+
+ // m ft km mi nm
+ conversion := []float64{1.0, 0.3048, 1000.0, 1609.344, 1852.0}
+ ellipsoid := Ellipsoid{e2, units, distUnits, longSym, bearSym, conversion[distUnits]}
+ return ellipsoid
+}
+
+/* Intermediate
+
+Takes two coordinates with longitude and latitude; and a step count and
+returns range and bearing and an array with the lons and lats of intermediate
+points on a straight line (whatever that is on an ellipsoid), INCLUDING the
+start and the endpoint.
+
+So if you put in point1 and point2 with step count 4, the output will be
+(you make 4 hops, right?)
+
+ point1
+ i1
+ i2
+ i3
+ point2
+
+Each point is two float64 values, lat and lon, thus you have an array
+with 4*2 + 2 = 5*2 cells.
+
+steps shall not be 0.
+
+I havent tested the upper limit for steps.
+
+*/
+func (ellipsoid Ellipsoid) Intermediate(lat1, lon1, lat2, lon2 float64, steps int) (distance, bearing float64, arr []float64) {
+ if steps == 0 {
+ return
+ }
+ r, phi := ellipsoid.To(lat1, lon1, lat2, lon2)
+ v := make([]float64, steps*2+2)
+ for i := 0; i <= steps; i++ {
+ a, b := ellipsoid.At(lat1, lon1, r*float64(i)/float64(steps), phi)
+ v[i*2], v[i*2+1] = a, b
+ }
+ arr = v
+ return r, phi, arr
+
+}
+
+/* To returns range, bearing between two specified locations.
+
+ dist, theta = geo.To( lat1, lon1, lat2, lon2 )
+
+*/
+func (ellipsoid Ellipsoid) To(lat1, lon1, lat2, lon2 float64) (distance, bearing float64) {
+
+ if ellipsoid.Units == Degrees {
+ lat1 = deg2rad(lat1)
+ lon1 = deg2rad(lon1)
+ lat2 = deg2rad(lat2)
+ lon2 = deg2rad(lon2)
+ }
+
+ distance, bearing = ellipsoid.calculateBearing(lat1, lon1, lat2, lon2)
+ if ellipsoid.Units == Degrees {
+ bearing = rad2deg(bearing)
+ }
+
+ distance /= ellipsoid.DistanceFactor
+
+ return
+}
+
+/* At returns the list latitude,longitude in degrees or radians that is a
+specified range and bearing from a given location.
+
+ lat2, lon2 = geo.At( lat1, lon1, range, bearing )
+
+*/
+func (ellipsoid Ellipsoid) At(lat1, lon1, distance, bearing float64) (lat2, lon2 float64) {
+
+ if ellipsoid.Units == Degrees {
+ lat1 = deg2rad(lat1)
+ lon1 = deg2rad(lon1)
+ bearing = deg2rad(bearing)
+ }
+
+ lat2, lon2 = ellipsoid.calculateTargetlocation(lat1, lon1, distance, bearing)
+
+ if ellipsoid.LongitudeSymmetric == LongitudeIsSymmetric {
+ if lon2 > pi {
+ lon2 -= twopi
+ }
+ }
+ if ellipsoid.LongitudeSymmetric == LongitudeNotSymmetric {
+ if lon2 < 0.0 {
+ lon2 += twopi
+ }
+ }
+
+ if ellipsoid.Units == Degrees {
+ lat2 = rad2deg(lat2)
+ lon2 = rad2deg(lon2)
+ }
+
+ return
+}
+
+/* Displacement returns the (x,y) displacement in distance units between the two specified
+locations.
+
+ x, y = geo.Displacement( lat1, lon1, lat2, lon2 )
+
+NOTE: The x and y displacements are only approximations and only valid
+between two locations that are fairly near to each other. Beyond 10 kilometers
+or more, the concept of X and Y on a curved surface loses its meaning.
+
+*/
+func (ellipsoid Ellipsoid) Displacement(lat1, lon1, lat2, lon2 float64) (x, y float64) {
+ // FIXME: Normalize!!! before use.
+ r, bearing := ellipsoid.To(lat1, lon1, lat2, lon2)
+
+ if ellipsoid.Units == Degrees {
+ bearing = deg2rad(bearing)
+ }
+
+ x = r * math.Sin(bearing)
+ y = r * math.Cos(bearing)
+ return x, y
+}
+
+/* Location returns the list (latitude,longitude) of a location at a given (x,y)
+displacement from a given location.
+
+ lat2, lon2 = geo.Location( lat1, lon1, x, y )
+
+The note from Displacement applies.
+
+*/
+func (ellipsoid Ellipsoid) Location(lat1, lon1, x, y float64) (lat, lon float64) {
+ degreesPerRadian := 180.0 / math.Pi
+
+ range1 := math.Sqrt(x*x + y*y)
+ bearing1 := math.Atan2(x, y)
+
+ if ellipsoid.Units == Degrees {
+ bearing1 *= degreesPerRadian
+ }
+
+ return ellipsoid.At(lat1, lon1, range1, bearing1)
+}
+
+func (ellipsoid Ellipsoid) calculateTargetlocation(lat1, lon1, distance, bearing float64) (lat2, lon2 float64) {
+
+ if debug == true {
+ fmt.Printf("_forward(lat1=%v,lon1=%v,range=%v,bearing=%v)\n", lat1, lon1, distance, bearing)
+ }
+
+ eps := 0.5e-13
+
+ a := ellipsoid.Ellipse.Equatorial
+ f := 1.0 / ellipsoid.Ellipse.InvFlattening
+ r := 1.0 - f
+
+ clat1 := math.Cos(lat1)
+ if clat1 == 0 {
+ fmt.Printf("WARNING: Division by Zero in ellipsoid.go.\n")
+ return 0.0, 0.0
+ }
+ tu := r * math.Sin(lat1) / clat1
+ faz := bearing
+
+ s := ellipsoid.DistanceFactor * distance
+
+ sf := math.Sin(faz)
+ cf := math.Cos(faz)
+
+ baz := 0.0
+ if cf != 0.0 {
+ baz = 2.0 * math.Atan2(tu, cf)
+ }
+
+ cu := 1.0 / math.Sqrt(1.0+tu*tu)
+ su := tu * cu
+ sa := cu * sf
+ c2a := 1.0 - (sa * sa)
+ x := 1.0 + math.Sqrt((((1.0/(r*r))-1.0)*c2a)+1.0)
+ x = (x - 2.0) / x
+ c := 1.0 - x
+ c = (((x * x) / 4.0) + 1.0) / c
+ d := x * ((0.375 * x * x) - 1.0)
+ tu = ((s / r) / a) / c
+ y := tu
+
+ if debug == true {
+ fmt.Printf("r=%.8f, tu=%.8f, faz=%.8f\n", r, tu, faz)
+ fmt.Printf("baz=%.8f, sf=%.8f, cf=%.8f\n", baz, sf, cf)
+ fmt.Printf("cu=%.8f, su=%.8f, sa=%.8f\n", cu, su, sa)
+ fmt.Printf("x=%.8f, c=%.8f, y=%.8f\n", x, c, y)
+ }
+
+ var cy, cz, e, sy float64
+ for true {
+ sy = math.Sin(y)
+ cy = math.Cos(y)
+ cz = math.Cos(baz + y)
+ e = (2.0 * cz * cz) - 1.0
+ c = y
+ x = e * cy
+ y = (2.0 * e) - 1.0
+ y = (((((((((sy * sy * 4.0) - 3.0) * y * cz * d) / 6.0) + x) * d) / 4.0) - cz) * sy * d) + tu
+
+ if math.Abs(y-c) <= eps {
+ break
+ }
+ }
+ baz = (cu * cy * cf) - (su * sy)
+ c = r * math.Sqrt((sa*sa)+(baz*baz))
+ d = su*cy + cu*sy*cf
+ lat2 = math.Atan2(d, c)
+ c = cu*cy - su*sy*cf
+ x = math.Atan2(sy*sf, c)
+ c = (((((-3.0 * c2a) + 4.0) * f) + 4.0) * c2a * f) / 16.0
+ d = ((((e * cy * c) + cz) * sy * c) + y) * sa
+ lon2 = lon1 + x - (1.0-c)*d*f
+
+ if debug == true {
+ fmt.Printf("returns(lat2=%v,lon2=%v)\n", lat2, lon2)
+ }
+ return lat2, lon2
+}
+
+func (ellipsoid Ellipsoid) calculateBearing(lat1, lon1, lat2, lon2 float64) (distance, bearing float64) {
+ a := ellipsoid.Ellipse.Equatorial
+ f := 1 / ellipsoid.Ellipse.InvFlattening
+
+ if lon1 < 0 {
+ lon1 += twopi
+ }
+ if lon2 < 0 {
+ lon2 += twopi
+ }
+
+ r := 1.0 - f
+ clat1 := math.Cos(lat1)
+ if clat1 == 0 {
+ fmt.Printf("WARNING: Division by Zero in ellipsoid.go.\n")
+ return 0.0, 0.0
+ }
+ clat2 := math.Cos(lat2)
+ if clat2 == 0 {
+ fmt.Printf("WARNING: Division by Zero in ellipsoid.go.\n")
+ return 0.0, 0.0
+ }
+ tu1 := r * math.Sin(lat1) / clat1
+ tu2 := r * math.Sin(lat2) / clat2
+ cu1 := 1.0 / (math.Sqrt((tu1 * tu1) + 1.0))
+ su1 := cu1 * tu1
+ cu2 := 1.0 / (math.Sqrt((tu2 * tu2) + 1.0))
+ s := cu1 * cu2
+ baz := s * tu2
+ faz := baz * tu1
+ dlon := lon2 - lon1
+
+ if debug == true {
+ fmt.Printf("a=%v, f=%v\n", a, f)
+ fmt.Printf("lat1=%v, lon1=%v\n", lat1, lon1)
+ fmt.Printf("lat2=%v, lon2=%v\n", lat2, lon2)
+
+ fmt.Printf("r=%v, tu1=%v, tu2=%v\n", r, tu1, tu2)
+ fmt.Printf("faz=%.8f, dlon=%.8f, su1=%v\n", faz, dlon, su1)
+ }
+
+ x := dlon
+ cnt := 0
+
+ var c2a, c, cx, cy, cz, d, del, e, sx, sy, y float64
+ // This originally was a do-while loop. Exit condition is at end of loop.
+ for true {
+ if debug == true {
+ fmt.Printf(" x=%.8f\n", x)
+ }
+ sx = math.Sin(x)
+ cx = math.Cos(x)
+ tu1 = cu2 * sx
+ tu2 = baz - (su1 * cu2 * cx)
+
+ if debug == true {
+ fmt.Printf(" sx=%.8f, cx=%.8f, tu1=%.8f, tu2=%.8f\n", sx, cx, tu1, tu2)
+ }
+
+ sy = math.Sqrt(tu1*tu1 + tu2*tu2)
+ cy = s*cx + faz
+ y = math.Atan2(sy, cy)
+ var sa float64
+ if sy == 0.0 {
+ sa = 1.0
+ } else {
+ sa = (s * sx) / sy
+ }
+
+ if debug == true {
+ fmt.Printf(" sy=%.8f, cy=%.8f, y=%.8f, sa=%.8f\n", sy, cy, y, sa)
+ }
+
+ c2a = 1.0 - (sa * sa)
+ cz = faz + faz
+ if c2a > 0.0 {
+ cz = ((-cz) / c2a) + cy
+ }
+ e = (2.0 * cz * cz) - 1.0
+ c = (((((-3.0 * c2a) + 4.0) * f) + 4.0) * c2a * f) / 16.0
+ d = x
+ x = ((e*cy*c+cz)*sy*c + y) * sa
+ x = (1.0-c)*x*f + dlon
+ del = d - x
+
+ if debug == true {
+ fmt.Printf(" c2a=%.8f, cz=%.8f\n", c2a, cz)
+ fmt.Printf(" e=%.8f, d=%.8f\n", e, d)
+ fmt.Printf(" (d-x)=%.8g\n", del)
+ }
+ if math.Abs(del) <= eps {
+ break
+ }
+ cnt++
+ if cnt > maxLoopCount {
+ break
+ }
+
+ }
+
+ faz = math.Atan2(tu1, tu2)
+ baz = math.Atan2(cu1*sx, (baz*cx-su1*cu2)) + pi
+ x = math.Sqrt(((1.0/(r*r))-1.0)*c2a+1.0) + 1.0
+ x = (x - 2.0) / x
+ c = 1.0 - x
+ c = ((x*x)/4.0 + 1.0) / c
+ d = ((0.375 * x * x) - 1.0) * x
+ x = e * cy
+
+ if debug == true {
+ fmt.Printf("e=%.8f, cy=%.8f, x=%.8f\n", e, cy, x)
+ fmt.Printf("sy=%.8f, c=%.8f, d=%.8f\n", sy, c, d)
+ fmt.Printf("cz=%.8f, a=%.8f, r=%.8f\n", cz, a, r)
+ }
+
+ s = 1.0 - e - e
+ s = ((((((((sy * sy * 4.0) - 3.0) * s * cz * d / 6.0) - x) * d / 4.0) + cz) * sy * d) + y) * c * a * r
+
+ if debug == true {
+ fmt.Printf("s=%.8f\n", s)
+ }
+
+ // adjust azimuth to (0,360) or (-180,180) as specified
+ if ellipsoid.BearingSymmetry == BearingIsSymmetric {
+ if faz < -(pi) {
+ faz += twopi
+ }
+ if faz >= pi {
+ faz -= twopi
+ }
+ } else {
+ if faz < 0 {
+ faz += twopi
+ }
+ if faz >= twopi {
+ faz -= twopi
+ }
+ }
+
+ distance, bearing = s, faz
+ return
+}
+
+/* ToLLA takes three cartesian coordinates x, y, z and returns
+the latitude, longitude, elevation list.
+
+FIXME: This algorithm cannot handle x==0, although this is a valid value.
+WARNING: I put in an if condition to catch this. Is it still necessary?
+
+*/
+func (ellipsoid Ellipsoid) ToLLA(x, y, z float64) (lat1, lon1, alt1 float64) {
+
+ if x == 0 {
+ fmt.Printf("FATAL: Caught x==0 (div by zero).\n")
+ return
+ }
+
+ a := ellipsoid.Ellipse.Equatorial
+ f := 1 / ellipsoid.Ellipse.InvFlattening
+
+ b := a * (1.0 - f)
+ e := math.Sqrt((a*a - b*b) / (a * a))
+ e2 := math.Sqrt((a*a - b*b) / (b * b)) // e'
+ esq := e * e // e squared
+ e2sq := e2 * e2 // e' squared
+ p := math.Sqrt(x*x + y*y)
+
+ theta := math.Atan2(z*a, p*b)
+ stheta3 := math.Sin(theta) * math.Sin(theta) * math.Sin(theta)
+ ctheta3 := math.Cos(theta) * math.Cos(theta) * math.Cos(theta)
+
+ lon1 = math.Atan2(y, x)
+ phi := math.Atan2(z+e2sq*b*stheta3, p-esq*a*ctheta3)
+ lat1 = phi
+
+ sphisq := math.Sin(phi) * math.Sin(phi)
+ N := a / (math.Sqrt(1 - esq*sphisq))
+ alt1 = p/math.Cos(phi) - N
+
+ if ellipsoid.LongitudeSymmetric == LongitudeIsSymmetric {
+ if lon1 > pi {
+ lon1 -= twopi
+ }
+ }
+ if ellipsoid.LongitudeSymmetric == LongitudeNotSymmetric {
+ if lon1 < 0.0 {
+ lon1 += twopi
+ }
+ }
+
+ if ellipsoid.Units == Degrees {
+ lat1 = rad2deg(lat1)
+ lon1 = rad2deg(lon1)
+ }
+ return lat1, lon1, alt1
+}
+
+/* ToECEF takes the latitude, longitude, elevation list and
+ returns three cartesian coordinates x, y, z */
+func (ellipsoid Ellipsoid) ToECEF(lat1, lon1, alt1 float64) (x, y, z float64) {
+ a := ellipsoid.Ellipse.Equatorial
+ f := 1 / ellipsoid.Ellipse.InvFlattening
+
+ b := a * (1.0 - f)
+ e := math.Sqrt((a*a - b*b) / (a * a))
+ esq := e * e // e squared
+
+ if ellipsoid.Units == Degrees {
+ lat1 = deg2rad(lat1)
+ lon1 = deg2rad(lon1)
+ }
+
+ h := alt1 // renamed for convenience
+ phi := lat1
+ lambda := lon1
+
+ cphi := math.Cos(phi)
+ sphi := math.Sin(phi)
+ sphisq := sphi * sphi
+ clam := math.Cos(lambda)
+ slam := math.Sin(lambda)
+ N := a / (math.Sqrt(1 - esq*sphisq))
+ x = (N + h) * cphi * clam
+ y = (N + h) * cphi * slam
+ z = ((b*b*N)/(a*a) + h) * sphi
+
+ return x, y, z
+}
+
+/*
+ DEFINED ELLIPSOIDS
+
+The following ellipsoids are defined in Geo::Ellipsoid, with the
+semi-major axis in meters and the reciprocal flattening as shown.
+
+
+ Ellipsoid Semi-Major Axis (m.) 1/Flattening
+ --------- ------------------- ---------------
+ AIRY 6377563.396 299.3249646
+ AIRY-MODIFIED 6377340.189 299.3249646
+ AUSTRALIAN 6378160.0 298.25
+ BESSEL-1841 6377397.155 299.1528128
+ CLARKE-1880 6378249.145 293.465
+ EVEREST-1830 6377276.345 290.8017
+ EVEREST-MODIFIED 6377304.063 290.8017
+ FISHER-1960 6378166.0 298.3
+ FISHER-1968 6378150.0 298.3
+ GRS80 6378137.0 298.25722210088
+ HOUGH-1956 6378270.0 297.0
+ HAYFORD 6378388.0 297.0
+ IAU76 6378140.0 298.257
+ KRASSOVSKY-1938 6378245.0 298.3
+ NAD27 6378206.4 294.9786982138
+ NWL-9D 6378145.0 298.25
+ SOUTHAMERICAN-1969 6378160.0 298.25
+ SOVIET-1985 6378136.0 298.257
+ WGS72 6378135.0 298.26
+ WGS84 6378137.0 298.257223563
+
+plus a few more ...
+
+ LIMITATIONS
+
+The methods should not be used on points which are too near the poles
+(above or below 89 degrees), and should not be used on points which
+are antipodal, i.e., exactly on opposite sides of the ellipsoid. The
+methods will not return valid results in these cases.
+
+The Go-version does not support all features of the Perl module. If you
+need advanced features, like defining your own ellipses at runtime,
+calculating x,y dislocations, etc, please refer to the package on CPAN
+
+Geo::Ellipsoid
+
+http://search.cpan.org/~jgibson/Geo-Ellipsoid-1.12/lib/Geo/Ellipsoid.pm
+
+FIXME: Add more checks for div by 0.
+
+ ACKNOWLEDGEMENTS
+
+The conversion algorithms used here are Perl translations of Fortran
+routines written by LCDR L. Pfeifer NGS Rockville MD that implement
+T. Vincenty's Modified Rainsford's method with Helmert's elliptical
+terms as published in "Direct and Inverse Solutions of Ellipsoid on
+the Ellipsoid with Application of Nested Equations", T. Vincenty,
+Survey Review, April 1975.
+
+The Fortran source code files inverse.for and forward.for
+may be obtained from
+
+ ftp://ftp.ngs.noaa.gov/pub/pcsoft/for_inv.3d/source/
+
+ AUTHOR
+
+Jim Gibson, <Jim@Gibson.org> (Perl version)
+Stefan Schroeder <ondekoza@gmail.com> (Port from Perl to Golang)
+
+ BUGS
+
+See LIMITATIONS, above.
+
+Please report any bugs or feature requests to
+the author.
+
+COPYRIGHT & LICENSE
+
+Copyright 2005-2008 Jim Gibson, all rights reserved.
+
+This program is free software; you can redistribute it and/or modify it
+under the same terms as Perl.
+
+*/