/* 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?

*/
769364
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
}
898824
func (p *MercatorProjection) FromLatLng(ll LatLng) r2.Point {
	// This formula is more accurate near zero than the log(tan()) version.
	// Note that latitudes of +/- 90 degrees yield "y" values of +/- infinity.
	sinPhi := math.Sin(float64(ll.Lat))
	y := 0.5 * math.Log((1+sinPhi)/(1-sinPhi))
	return r2.Point{p.fromRadians * float64(ll.Lng), p.fromRadians * y}
}
601919
func ParseGeometry(gi interface{}) (geom Geometry, err error) {
	switch g := gi.(type) {
	case map[string]interface{}:
		coord := g["coordinates"]
		switch typ := g["type"]; typ {
		case "Point":
			return parsePoint(coord)
		case "LineString":
			return parseLineString(coord)
		case "MultiPoint":
			return parseMultiPoint(coord)
		case "MultiLineString":
			return parseMultiLineString(coord)
		case "Polygon":
			return parsePolygon(coord)
		case "MultiPolygon":
			return parseMultiPolygon(coord)
		case "GeometryCollection":
			return parseGeometryCollection(g["geometries"])
		default:
			err = errors.New(fmt.Sprintf("ParseError: Unknown geometry type %s", typ))
			break
		}
	}
	return geom, err
}
997321
func intersectsLatEdge(a, b Point, lat s1.Angle, lng s1.Interval) bool {
	// Unfortunately, lines of constant latitude are curves on
	// the sphere. They can intersect a straight edge in 0, 1, or 2 points.

	// First, compute the normal to the plane AB that points vaguely north.
	z := Point{a.PointCross(b).Normalize()}
	if z.Z < 0 {
		z = Point{z.Mul(-1)}
	}

	// Extend this to an orthonormal frame (x,y,z) where x is the direction
	// where the great circle through AB achieves its maximium latitude.
	y := Point{z.PointCross(PointFromCoords(0, 0, 1)).Normalize()}
	x := y.Cross(z.Vector)

	// Compute the angle "theta" from the x-axis (in the x-y plane defined
	// above) where the great circle intersects the given line of latitude.
	sinLat := math.Sin(float64(lat))
	if math.Abs(sinLat) >= x.Z {
		// The great circle does not reach the given latitude.
		return false
	}

	cosTheta := sinLat / x.Z
	sinTheta := math.Sqrt(1 - cosTheta*cosTheta)
	theta := math.Atan2(sinTheta, cosTheta)

	// The candidate intersection points are located +/- theta in the x-y
	// plane. For an intersection to be valid, we need to check that the
	// intersection point is contained in the interior of the edge AB and
	// also that it is contained within the given longitude interval "lng".

	// Compute the range of theta values spanned by the edge AB.
	abTheta := s1.IntervalFromPointPair(
		math.Atan2(a.Dot(y.Vector), a.Dot(x)),
		math.Atan2(b.Dot(y.Vector), b.Dot(x)))

	if abTheta.Contains(theta) {
		// Check if the intersection point is also in the given lng interval.
		isect := x.Mul(cosTheta).Add(y.Mul(sinTheta))
		if lng.Contains(math.Atan2(isect.Y, isect.X)) {
			return true
		}
	}

	if abTheta.Contains(-theta) {
		// Check if the other intersection point is also in the given lng interval.
		isect := x.Mul(cosTheta).Sub(y.Mul(sinTheta))
		if lng.Contains(math.Atan2(isect.Y, isect.X)) {
			return true
		}
	}
	return false
}
238422
func (s GeoDistanceSort) Point(lat, lon float64) GeoDistanceSort {
	s.points = append(s.points, GeoPointFromLatLon(lat, lon))
	return s
}
519843
func (p *Point) GeoDistanceFrom(point *Point, haversine ...bool) float64 {
	dLat := deg2rad(point.Lat() - p.Lat())
	dLng := math.Abs(deg2rad(point.Lng() - p.Lng()))

	if yesHaversine(haversine) {
		// yes trig functions
		dLat2Sin := math.Sin(dLat / 2)
		dLng2Sin := math.Sin(dLng / 2)
		a := dLat2Sin*dLat2Sin + math.Cos(deg2rad(p.Lat()))*math.Cos(deg2rad(point.Lat()))*dLng2Sin*dLng2Sin

		return 2.0 * EarthRadius * math.Atan2(math.Sqrt(a), math.Sqrt(1-a))
	}

	if dLng > math.Pi {
		dLng = 2*math.Pi - dLng
	}

	// fast way using pythagorean theorem on an equirectangular projection
	x := dLng * math.Cos(deg2rad((p.Lat()+point.Lat())/2.0))
	return math.Sqrt(dLat*dLat+x*x) * EarthRadius
}
166132
func InterpolateAtDistance(ax s1.Angle, a, b Point) Point {
	aRad := ax.Radians()

	// Use PointCross to compute the tangent vector at A towards B. The
	// result is always perpendicular to A, even if A=B or A=-B, but it is not
	// necessarily unit length. (We effectively normalize it below.)
	normal := a.PointCross(b)
	tangent := normal.Vector.Cross(a.Vector)

	// Now compute the appropriate linear combination of A and "tangent". With
	// infinite precision the result would always be unit length, but we
	// normalize it anyway to ensure that the error is within acceptable bounds.
	// (Otherwise errors can build up when the result of one interpolation is
	// fed into another interpolation.)
	return Point{(a.Mul(math.Cos(aRad)).Add(tangent.Mul(math.Sin(aRad) / tangent.Norm()))).Normalize()}
}
916872
func locator(p Point, precision int) (string, error) {
	if math.IsNaN(p.Latitude) {
		return "", errors.New("maidenhead: latitude is not a digit")
	}
	if math.IsInf(p.Latitude, 0) {
		return "", errors.New("maidenhead: latitude is infinite")
	}
	if math.IsNaN(p.Longitude) {
		return "", errors.New("maidenhead: longitude is not a digit")
	}
	if math.IsInf(p.Longitude, 0) {
		return "", errors.New("maidenhead: longitude is infinite")
	}
	if math.Abs(p.Latitude) == 90.0 {
		return "", errors.New("maidenhead: grid square invalid at poles")
	} else if math.Abs(p.Latitude) > 90.0 {
		return "", fmt.Errorf("maidenhead: invalid latitude %.04f", p.Latitude)
	} else if math.Abs(p.Longitude) > 180.0 {
		return "", fmt.Errorf("maidenhead: invalid longitude %.05f", p.Longitude)
	}

	var (
		lat = p.Latitude + 90.0
		lng = p.Longitude + 180.0
		loc string
	)

	lat = lat/10.0 + 0.0000001
	lng = lng/20.0 + 0.0000001
	loc = loc + string(upper[int(lng)]) + string(upper[int(lat)])
	if precision == 1 {
		return loc, nil
	}
	lat = 10 * (lat - math.Floor(lat))
	lng = 10 * (lng - math.Floor(lng))
	loc = loc + fmt.Sprintf("%d%d", int(lng)%10, int(lat)%10)
	if precision == 2 {
		return loc, nil
	}
	lat = 24 * (lat - math.Floor(lat))
	lng = 24 * (lng - math.Floor(lng))
	loc = loc + string(upper[int(lng)]) + string(upper[int(lat)])
	if precision == 3 {
		return loc, nil
	}
	lat = 10 * (lat - math.Floor(lat))
	lng = 10 * (lng - math.Floor(lng))
	loc = loc + fmt.Sprintf("%d%d", int(lng)%10, int(lat)%10)
	if precision == 4 {
		return loc, nil
	}
	lat = 24 * (lat - math.Floor(lat))
	lng = 24 * (lng - math.Floor(lng))
	loc = loc + string(lower[int(lng)]) + string(lower[int(lat)])
	return loc, nil
}
1237738
func EncodeMultiPolygon(mp *geom.MultiPolygon) kml.Element {
	polygons := make([]kml.Element, mp.NumPolygons())
	flatCoords := mp.FlatCoords()
	endss := mp.Endss()
	stride := mp.Stride()
	d := dim(mp.Layout())
	offset := 0
	for i, ends := range endss {
		boundaries := make([]kml.Element, len(ends))
		for j, end := range ends {
			linearRing := kml.LinearRing(kml.CoordinatesFlat(flatCoords, offset, end, stride, d))
			if j == 0 {
				boundaries[j] = kml.OuterBoundaryIs(linearRing)
			} else {
				boundaries[j] = kml.InnerBoundaryIs(linearRing)
			}
			offset = end
		}
		polygons[i] = kml.Polygon(boundaries...)
	}
	return kml.MultiGeometry(polygons...)
}
685648
func TriMa(real []float64, timePeriod int32) []float64 {
	var outBegIdx C.int
	var outNBElement C.int
	outReal := make([]float64, len(real))
	C.TA_TRIMA(0, C.int(len(real)-1), (*C.double)(unsafe.Pointer(&real[0])), C.int(timePeriod), &outBegIdx, &outNBElement, (*C.double)(unsafe.Pointer(&outReal[0])))
	return outReal[:outNBElement]
}
830651
func (i Address) Latitude(v reflect.Value) (interface{}, error) {
	kind := v.Kind()
	val := i.latitute()
	if kind == reflect.Float32 {
		return val, nil
	}
	return float64(val), nil
}
1033210
func (r Rect) Lo() LatLng {
	return LatLng{s1.Angle(r.Lat.Lo) * s1.Radian, s1.Angle(r.Lng.Lo) * s1.Radian}
}
1293070
func Atan2(x, y float32) float32 {
	return float32(math.Atan2(float64(x), float64(y)))
}
605256
func CartesianToCylindical(coord Vec3) (rho, phi, z float32) {
	rho = float32(math.Hypot(float64(coord[0]), float64(coord[1])))

	phi = float32(math.Atan2(float64(coord[1]), float64(coord[0])))

	z = coord[2]

	return
}
947459
func (r Rect) Vertex(i int) LatLng {
	var lat, lng float64

	switch i {
	case 0:
		lat = r.Lat.Lo
		lng = r.Lng.Lo
	case 1:
		lat = r.Lat.Lo
		lng = r.Lng.Hi
	case 2:
		lat = r.Lat.Hi
		lng = r.Lng.Hi
	case 3:
		lat = r.Lat.Hi
		lng = r.Lng.Lo
	}
	return LatLng{s1.Angle(lat) * s1.Radian, s1.Angle(lng) * s1.Radian}
}
1346847
func (p *PlateCarreeProjection) ToLatLng(pt r2.Point) LatLng {
	return LatLng{
		Lat: s1.Angle(p.toRadians * pt.Y),
		Lng: s1.Angle(p.toRadians * math.Remainder(pt.X, p.xWrap)),
	}
}
1254444
func (c cachedGeocoder) ReverseGeocode(lat, lng float64) (*geo.Address, error) {
	// Check if we've cached this response
	locKey := fmt.Sprintf("geo.Location{%f,%f}", lat, lng)
	if cachedAddr, found := c.Cache.Get(locKey); found {
		return cachedAddr.(*geo.Address), nil
	}

	if addr, err := c.Geocoder.ReverseGeocode(lat, lng); err != nil {
		return nil, err
	} else {
		c.Cache.Set(locKey, addr, 0)
		return addr, nil
	}
}