// toCartesian converts lat/lon coordinates to cartesian x/y/z coordinates.
// It returns a vector representing a lat/lon point on x-y-z axes in metres
// from the earth centre.
497135
func toCartesian(coords Coordinates, datum Datum) Vector3d {

	aa := datum.a
	bb := datum.b

	φ := toRadians(coords.Lat)
	λ := toRadians(coords.Lon)
	sinφ := math.Sin(φ)
	cosφ := math.Cos(φ)
	sinλ := math.Sin(λ)
	cosλ := math.Cos(λ)

	eSq := (aa*aa - bb*bb) / (aa * aa) // make a and b const datum properties
	ν := aa / math.Sqrt(1-eSq*sinφ*sinφ)
	x := ν * cosφ * cosλ
	y := ν * cosφ * sinλ
	z := (1 - eSq) * ν * sinφ

	// apply Helmert transform
	// move to another function
	rx := toRadians((datum.rx * -1) / 3600) // normalise seconds to radians
	ry := toRadians((datum.ry * -1) / 3600) // normalise seconds to radians
	rz := toRadians((datum.rz * -1) / 3600) // normalise seconds to radians
	s1 := (datum.s*-1)/1e6 + 1              // normalise ppm to (s+1)

	// apply transform
	x2 := (datum.tx * -1) + x*s1 - y*rz + z*ry
	y2 := (datum.ty * -1) + x*rz + y*s1 - z*rx
	z2 := (datum.tz * -1) - x*ry + y*rx + z*s1

	// instantiate new vector
	vect := Vector3d{x: x2, y: y2, z: z2}

	return vect
}
883095
func (t T) Circle(center kml.Coordinate, radius, maxErr float64) []kml.Coordinate {
	numVertices := int(math.Ceil(math.Pi / math.Acos((radius-maxErr)/(radius+maxErr))))
	cs := make([]kml.Coordinate, numVertices+1)
	for i := 0; i < numVertices; i++ {
		cs[i] = t.Offset(center, radius, 360*float64(i)/float64(numVertices))
	}
	cs[numVertices] = cs[0]
	return cs
}
502711
func (g Geometry) MarshalJSON() ([]byte, error) {
	if g.Coordinates == nil && len(g.Geometries) == 0 {
		return []byte(`null`), nil
	}

	ng := &jsonGeometryMarshall{}
	switch g := g.Coordinates.(type) {
	case orb.Ring:
		ng.Coordinates = orb.Polygon{g}
	case orb.Bound:
		ng.Coordinates = g.ToPolygon()
	case orb.Collection:
		ng.Geometries = make([]*Geometry, 0, len(g))
		for _, c := range g {
			ng.Geometries = append(ng.Geometries, NewGeometry(c))
		}
		ng.Type = g.GeoJSONType()
	default:
		ng.Coordinates = g
	}

	if ng.Coordinates != nil {
		ng.Type = ng.Coordinates.GeoJSONType()
	}

	if len(g.Geometries) > 0 {
		ng.Geometries = g.Geometries
		ng.Type = orb.Collection{}.GeoJSONType()
	}
	return json.Marshal(ng)
}
700992
func (sc SpaceComponent) Corners() (points [4]engo.Point) {
	points[0].X = sc.Position.X
	points[0].Y = sc.Position.Y

	sin, cos := math.Sincos(sc.Rotation * math.Pi / 180)

	points[1].X = points[0].X + sc.Width*cos
	points[1].Y = points[0].Y + sc.Width*sin

	points[2].X = points[0].X - sc.Height*sin
	points[2].Y = points[0].Y + sc.Height*cos

	points[3].X = points[0].X + sc.Width*cos - sc.Height*sin
	points[3].Y = points[0].Y + sc.Height*cos + sc.Width*sin

	return
}
734904
func Geometry(g orb.Geometry, z maptile.Zoom) maptile.Set {
	if g == nil {
		return nil
	}

	switch g := g.(type) {
	case orb.Point:
		return Point(g, z)
	case orb.MultiPoint:
		return MultiPoint(g, z)
	case orb.LineString:
		return LineString(g, z)
	case orb.MultiLineString:
		return MultiLineString(g, z)
	case orb.Ring:
		return Ring(g, z)
	case orb.Polygon:
		return Polygon(g, z)
	case orb.MultiPolygon:
		return MultiPolygon(g, z)
	case orb.Collection:
		return Collection(g, z)
	case orb.Bound:
		return Bound(g, z)
	}

	panic(fmt.Sprintf("geometry type not supported: %T", g))
}
1277399
func QuatRotate(angle float32, axis Vec3) Quat {
	// angle = (float32(math.Pi) * angle) / 180.0

	c, s := float32(math.Cos(float64(angle/2))), float32(math.Sin(float64(angle/2)))

	return Quat{c, axis.Mul(s)}
}
78166
func (c Client) Cos(x, y Expr) (float32, error) {
	req := cosQuery{A: x, B: y}
	body, err := c.fetch(req, "cos")
	if err != nil {
		return 0.0, err
	}

	var data cosResponse
	err = json.Unmarshal(body, &data)
	if err != nil {
		return 0.0, fmt.Errorf("error unmarshalling result: %v", err)
	}

	return data.Value, nil
}
107430
func Cos(z, x *decimal.Big) *decimal.Big {
	if z.CheckNaNs(x, nil) {
		return z
	}

	if x.IsInf(0) {
		z.Context.Conditions |= decimal.InvalidOperation
		return z.SetNaN(false)
	}

	ctx := decimal.Context{Precision: precision(z) + defaultExtraPrecision}

	negXSq, halved, ok := prepCosine(z, x, ctx)
	if !ok {
		z.Context.Conditions |= decimal.InvalidOperation
		return z.SetNaN(false)
	}

	ctx.Precision += halved
	z.Copy(BinarySplitDynamic(ctx,
		func(_ uint64) *decimal.Big { return one },
		getCosineP(negXSq),
		func(_ uint64) *decimal.Big { return one },
		getCosineQ(ctx),
	))

	// now undo the half angle bit
	for i := 0; i < halved; i++ {
		ctx.Mul(z, z, z)
		ctx.Mul(z, z, two)
		ctx.Sub(z, z, one)
	}
	ctx.Precision -= defaultExtraPrecision + halved
	return ctx.Round(z)
}
760743
func (v PreciseVector) Cross(ov PreciseVector) PreciseVector {
	return PreciseVector{
		X: precSub(precMul(v.Y, ov.Z), precMul(v.Z, ov.Y)),
		Y: precSub(precMul(v.Z, ov.X), precMul(v.X, ov.Z)),
		Z: precSub(precMul(v.X, ov.Y), precMul(v.Y, ov.X)),
	}
}
340009
func (g *geom1) Coords() []Coord {
	return inflate1(g.flatCoords, 0, len(g.flatCoords), g.stride)
}
1033210
func (r Rect) Lo() LatLng {
	return LatLng{s1.Angle(r.Lat.Lo) * s1.Radian, s1.Angle(r.Lng.Lo) * s1.Radian}
}
67896
func Midpoint(p, p2 orb.Point) orb.Point {
	dLon := deg2rad(p2[0] - p[0])

	aLatRad := deg2rad(p[1])
	bLatRad := deg2rad(p2[1])

	x := math.Cos(bLatRad) * math.Cos(dLon)
	y := math.Cos(bLatRad) * math.Sin(dLon)

	r := orb.Point{
		deg2rad(p[0]) + math.Atan2(y, math.Cos(aLatRad)+x),
		math.Atan2(math.Sin(aLatRad)+math.Sin(bLatRad), math.Sqrt((math.Cos(aLatRad)+x)*(math.Cos(aLatRad)+x)+y*y)),
	}

	// convert back to degrees
	r[0] = rad2deg(r[0])
	r[1] = rad2deg(r[1])

	return r
}
33026
func (mat *T) Quaternion() quaternion.T {
	tr := mat.Trace()

	s := math.Sqrt(tr + 1)
	w := s * 0.5
	s = 0.5 / s

	q := quaternion.T{
		(mat[1][2] - mat[2][1]) * s,
		(mat[2][0] - mat[0][2]) * s,
		(mat[0][1] - mat[1][0]) * s,
		w,
	}
	return q.Normalized()
}
953351
func Rotate(p, axis Point, angle s1.Angle) Point {
	// Let M be the plane through P that is perpendicular to axis, and let
	// center be the point where M intersects axis. We construct a
	// right-handed orthogonal frame (dx, dy, center) such that dx is the
	// vector from center to P, and dy has the same length as dx. The
	// result can then be expressed as (cos(angle)*dx + sin(angle)*dy + center).
	center := axis.Mul(p.Dot(axis.Vector))
	dx := p.Sub(center)
	dy := axis.Cross(p.Vector)
	// Mathematically the result is unit length, but normalization is necessary
	// to ensure that numerical errors don't accumulate.
	return Point{dx.Mul(math.Cos(angle.Radians())).Add(dy.Mul(math.Sin(angle.Radians()))).Add(center).Normalize()}
}
244535
func (r Region) MidIn(point int) bool {
	return point > r.From && point < r.To
}
834562
func Distance(p1, p2 orb.Point) float64 {
	dLat := deg2rad(p1[1] - p2[1])
	dLon := deg2rad(p1[0] - p2[0])

	dLon = math.Abs(dLon)
	if dLon > math.Pi {
		dLon = 2*math.Pi - dLon
	}

	// fast way using pythagorean theorem on an equirectangular projection
	x := dLon * math.Cos(deg2rad((p1[1]+p2[1])/2.0))
	return math.Sqrt(dLat*dLat+x*x) * orb.EarthRadius
}
221238
func lagrange(pts []point, x float64) (y float64) {
	y = 0.0
	for j := range pts {
		l := 1.0
		for m := range pts {
			if j == m {
				continue
			}
			l *= (x - pts[m].X) / (pts[j].X - pts[m].X)
		}
		y += pts[j].Y * l
	}
	return y
}