Golang math.Tan函数代码示例

func (p Place) Azimuth(point Place) float64 {

	if (p.Latitude == point.Latitude) && (p.Longitude == point.Longitude) {
		return 0.0
	}

	esq := (1.0 - 1.0/298.25) * (1.0 - 1.0/298.25)
	alat3 := math.Atan(math.Tan(p.Latitude/RadiansToDegrees)*esq) * RadiansToDegrees
	alat4 := math.Atan(math.Tan(point.Latitude/RadiansToDegrees)*esq) * RadiansToDegrees

	rlat1 := alat3 / RadiansToDegrees
	rlat2 := alat4 / RadiansToDegrees
	rdlon := (point.Longitude - p.Longitude) / RadiansToDegrees

	clat1 := math.Cos(rlat1)
	clat2 := math.Cos(rlat2)
	slat1 := math.Sin(rlat1)
	slat2 := math.Sin(rlat2)
	cdlon := math.Cos(rdlon)
	sdlon := math.Sin(rdlon)

	yazi := sdlon * clat2
	xazi := clat1*slat2 - slat1*clat2*cdlon

	azi := RadiansToDegrees * math.Atan2(yazi, xazi)

	if azi < 0.0 {
		azi += 360.0
	}

	return azi
}

func (p Place) Distance(point Place) float64 {

	if (p.Latitude == point.Latitude) && (p.Longitude == point.Longitude) {
		return 0.0
	}

	esq := (1.0 - 1.0/298.25) * (1.0 - 1.0/298.25)
	alat3 := math.Atan(math.Tan(p.Latitude/RadiansToDegrees)*esq) * RadiansToDegrees
	alat4 := math.Atan(math.Tan(point.Latitude/RadiansToDegrees)*esq) * RadiansToDegrees

	rlat1 := alat3 / RadiansToDegrees
	rlat2 := alat4 / RadiansToDegrees
	rdlon := (point.Longitude - p.Longitude) / RadiansToDegrees

	clat1 := math.Cos(rlat1)
	clat2 := math.Cos(rlat2)
	slat1 := math.Sin(rlat1)
	slat2 := math.Sin(rlat2)
	cdlon := math.Cos(rdlon)

	cdel := slat1*slat2 + clat1*clat2*cdlon
	switch {
	case cdel > 1.0:
		cdel = 1.0
	case cdel < -1.0:
		cdel = -1.0
	}

	return RadiansToKm * math.Acos(cdel)
}

func (p Place) BackAzimuth(point Place) float64 {

	if (p.Latitude == point.Latitude) && (p.Longitude == point.Longitude) {
		return 0.0
	}

	esq := (1.0 - 1.0/298.25) * (1.0 - 1.0/298.25)
	alat3 := math.Atan(math.Tan(p.Latitude/RadiansToDegrees)*esq) * RadiansToDegrees
	alat4 := math.Atan(math.Tan(point.Latitude/RadiansToDegrees)*esq) * RadiansToDegrees

	rlat1 := alat3 / RadiansToDegrees
	rlat2 := alat4 / RadiansToDegrees
	rdlon := (point.Longitude - p.Longitude) / RadiansToDegrees

	clat1 := math.Cos(rlat1)
	clat2 := math.Cos(rlat2)
	slat1 := math.Sin(rlat1)
	slat2 := math.Sin(rlat2)
	cdlon := math.Cos(rdlon)
	sdlon := math.Sin(rdlon)

	ybaz := -sdlon * clat1
	xbaz := clat2*slat1 - slat2*clat1*cdlon

	baz := RadiansToDegrees * math.Atan2(ybaz, xbaz)

	if baz < 0.0 {
		baz += 360.0
	}

	return baz
}

// Generates sample from general Levy-stable distibution.
func Sample(alpha float64, beta float64, mu float64, sigma float64) float64 {
	if beta == 0.0 {
		return symmetric(alpha, mu, sigma)
	}

	v := math.Pi * (rand.Float64() - 0.5)

	w := 0.0
	for w == 0.0 {
		w = rand.ExpFloat64()
	}

	if alpha == 1.0 {
		x := ((0.5*math.Pi+beta*v)*math.Tan(v) -
			beta*math.Log(0.5*math.Pi*w*math.Cos(v)/(0.5*math.Pi+beta*v))) / (0.5 * math.Pi)
		return sigma*x + beta*sigma*math.Log(sigma)/(0.5*math.Pi) + mu
	}

	t := beta * math.Tan(0.5*math.Pi*alpha)
	s := math.Pow(1.0+t*t, 1.0/(2.0*alpha))
	b := math.Atan(t) / alpha
	x := s * math.Sin(alpha*(v+b)) *
		math.Pow(math.Cos(v-alpha*(v+b))/w, (1.0-alpha)/alpha) /
		math.Pow(math.Cos(v), 1.0/alpha)
	return sigma*x + mu
}

// Project finds radon projections of N lines running through the image
// center for lines angled 0 to 180 degrees from horizontal.
// Totally not working
//
//  /param img - img.Image src image
//  /param  N  - int number of angled lines to consider.
//  /return sinogram - Projections
//  /return nbPerLine - number of pixel per line in sinogram
//  /return error - if failed
func Project(img image.Image, N int) (*image.Gray, []uint8, error) {
	if N == 0 {
		N = 180
	}
	var err error
	size := img.Bounds().Size()

	D := max(size.X, size.Y)
	xCenter, yCenter := float64(size.X)/2.0, float64(size.Y)/2.0

	xOff := int(math.Floor(xCenter + roundingFactor(xCenter)))
	yOff := int(math.Floor(yCenter + roundingFactor(yCenter)))

	radonMap := image.NewGray(image.Rect(0, 0, N, D))
	if err != nil {
		return nil, nil, err
	}

	nbPerLine := make([]uint8, N)

	for k := 0; k < (N/4)+1; k++ {
		θ := float64(k) * math.Pi / float64(N)
		α := math.Tan(θ)
		for x := 0; x < D; x++ {

			y := α * float64(x-xOff)
			yd := int(math.Floor(y + roundingFactor(y)))

			if (yd+yOff >= 0) && (yd+yOff < size.Y) && (x < size.X) {
				addPixelsToGray(img, x, yd+yOff, *radonMap, k, x)
				nbPerLine[k]++
			}
			if (yd+xOff >= 0) && (yd+xOff < size.X) && (k != N/4) && (x < size.Y) {
				addPixelsToGray(img, yd+xOff, x, *radonMap, N/2-k, x)
				nbPerLine[N/2-k]++
			}

		}
	}
	for j, k := 0, 3*N/4; k < N; k++ {
		θ := float64(k) * math.Pi / float64(N)
		α := math.Tan(θ)
		for x := 0; x < D; x++ {
			y := α * float64(x-xOff)
			yd := int(math.Floor(y + roundingFactor(y)))
			if (yd+yOff >= 0) && (yd+yOff < size.Y) && (x < size.X) {
				addPixelsToGray(img, x, yd+yOff, *radonMap, k, x)
				nbPerLine[k]++
			}
			if (yOff-yd >= 0) && (yOff-yd < size.X) && (2*yOff-x >= 0) && (2*yOff-x < size.Y) && (k != 3*N/4) {
				addPixelsToGray(img, -yd+yOff, -(x-yOff)+yOff, *radonMap, k-j, x)
				nbPerLine[k-j]++
			}

		}
		j += 2
	}

	return radonMap, nbPerLine, nil
}

func (p *PerspectiveCamera) UpdateProjectionMatrix() {

	var fov = math3d.RadToDeg(2 * math.Atan(math.Tan(math3d.DegToRad(p.Fov)*0.5)/p.Zoom))

	if p.FullWidth > 0 {

		var aspect = p.FullWidth / p.FullHeight
		var top = math.Tan(math3d.DegToRad(fov*0.5)) * p.Near
		var bottom = -top
		var left = aspect * bottom
		var right = aspect * top
		var width = math.Abs(right - left)
		var height = math.Abs(top - bottom)

		p.ProjectionMatrix.MakeFrustum(
			left+p.X*width/p.FullWidth,
			left+(p.X+p.Width)*width/p.FullWidth,
			top-(p.Y+p.Height)*height/p.FullHeight,
			top-p.Y*height/p.FullHeight,
			p.Near,
			p.Far,
		)
	} else {
		p.ProjectionMatrix.MakePerspective(fov, p.Aspect, p.Near, p.Far)
	}
}

func (self *Calc) CalculateFunction(funcName string, arg32 float32) float32 {
	if len(self.errors) > 0 {
		return 1
	}

	var result float64
	arg := float64(arg32)
	switch funcName {
	case "sin":
		result = math.Sin(self.DegToRad(arg))
	case "cos":
		result = math.Cos(self.DegToRad(arg))
	case "tg":
		result = math.Tan(self.DegToRad(arg))
	case "ctg":
		result = 1.0 / math.Tan(self.DegToRad(arg))
	case "arcsin":
		result = self.RadToDeg(math.Asin(arg))
	case "arccos":
		result = self.RadToDeg(math.Acos(arg))
	case "arctg":
		result = self.RadToDeg(math.Atan(arg))
	case "arcctg":
		result = self.RadToDeg(math.Atan(1 / arg))
	case "sqrt":
		result = math.Sqrt(arg)
	default:
		self.errors = append(self.errors, "Unknown identifier "+funcName)
		return 1
	}

	return float32(result)
}

func (m *Context) determineZoom(bounds s2.Rect, center s2.LatLng) int {
	b := bounds.AddPoint(center)
	if b.IsEmpty() || b.IsPoint() {
		return 15
	}

	tileSize := m.tileProvider.TileSize
	margin := 4.0 + m.determineExtraMarginPixels()
	w := (float64(m.width) - 2.0*margin) / float64(tileSize)
	h := (float64(m.height) - 2.0*margin) / float64(tileSize)
	minX := (b.Lo().Lng.Degrees() + 180.0) / 360.0
	maxX := (b.Hi().Lng.Degrees() + 180.0) / 360.0
	minY := (1.0 - math.Log(math.Tan(b.Lo().Lat.Radians())+(1.0/math.Cos(b.Lo().Lat.Radians())))/math.Pi) / 2.0
	maxY := (1.0 - math.Log(math.Tan(b.Hi().Lat.Radians())+(1.0/math.Cos(b.Hi().Lat.Radians())))/math.Pi) / 2.0
	dx := math.Abs(maxX - minX)
	dy := math.Abs(maxY - minY)

	zoom := 1
	for zoom < 30 {
		tiles := float64(uint(1) << uint(zoom))
		if dx*tiles > w || dy*tiles > h {
			return zoom - 1
		}
		zoom = zoom + 1
	}

	return 15
}

// AnomalyDistance returns true anomaly and distance for near-parabolic orbits.
//
// Distance r returned in AU.
// An error is returned if the algorithm fails to converge.
func (e *Elements) AnomalyDistance(jde float64) (ν unit.Angle, r float64, err error) {
	// fairly literal translation of code on p. 246
	q1 := base.K * math.Sqrt((1+e.Ecc)/e.PDis) / (2 * e.PDis) // line 20
	g := (1 - e.Ecc) / (1 + e.Ecc)                            // line 20

	t := jde - e.TimeP // line 22
	if t == 0 {        // line 24
		return 0, e.PDis, nil
	}
	d1, d := 10000., 1e-9        // line 14
	q2 := q1 * t                 // line 28
	s := 2. / (3 * math.Abs(q2)) // line 30
	s = 2 / math.Tan(2*math.Atan(math.Cbrt(math.Tan(math.Atan(s)/2))))
	if t < 0 { // line 34
		s = -s
	}
	if e.Ecc != 1 { // line 36
		l := 0 // line 38
		for {
			s0 := s // line 40
			z := 1.
			y := s * s
			g1 := -y * s
			q3 := q2 + 2*g*s*y/3 // line 42
			for {
				z += 1                          // line 44
				g1 = -g1 * g * y                // line 46
				z1 := (z - (z+1)*g) / (2*z + 1) // line 48
				f := z1 * g1                    // line 50
				q3 += f                         // line 52
				if z > 50 || math.Abs(f) > d1 { // line 54
					return 0, 0, errors.New("No convergence")
				}
				if math.Abs(f) <= d { // line 56
					break
				}
			}
			l++ // line 58
			if l > 50 {
				return 0, 0, errors.New("No convergence")
			}
			for {
				s1 := s // line 60
				s = (2*s*s*s/3 + q3) / (s*s + 1)
				if math.Abs(s-s1) <= d { // line 62
					break
				}
			}
			if math.Abs(s-s0) <= d { // line 64
				break
			}
		}
	}
	ν = unit.Angle(2 * math.Atan(s))               // line 66
	r = e.PDis * (1 + e.Ecc) / (1 + e.Ecc*ν.Cos()) // line 68
	if ν < 0 {                                     // line 70
		ν += 2 * math.Pi
	}
	return
}

// Equatorial computes data for a sundial level with the equator.
//
// Argument φ is geographic latitude at which the sundial will be located;
// a is the length of a straight stylus perpendicular to the plane of the
// sundial.
//
// The sundial will have two sides, north and south.  Results n and s define
// lines on the north and south sides of the sundial.  Result coordinates
// are in units of a, the stylus length.
func Equatorial(φ, a float64) (n, s []Line) {
	tφ := math.Tan(φ)
	for i := 0; i < 24; i++ {
		nl := Line{Hour: i}
		sl := Line{Hour: i}
		H := float64(i-12) * 15 * math.Pi / 180
		aH := math.Abs(H)
		sH, cH := math.Sincos(H)
		for _, d := range m {
			tδ := math.Tan(d * math.Pi / 180)
			H0 := math.Acos(-tφ * tδ)
			if aH > H0 {
				continue
			}
			x := -a * sH / tδ
			yy := a * cH / tδ
			if tδ < 0 {
				sl.Points = append(sl.Points, Point{x, yy})
			} else {
				nl.Points = append(nl.Points, Point{x, -yy})
			}
		}
		if len(nl.Points) > 0 {
			n = append(n, nl)
		}
		if len(sl.Points) > 0 {
			s = append(s, sl)
		}
	}
	return
}

func argsCheckAndAddColor(org float32, index int, anivalue float32) float32 {
	if len(os.Args) < index*2+1 {
		return org
	}
	targ, err := strconv.ParseFloat(os.Args[index*2], 32)
	if err != nil {
		return org
	}
	anivalue += 1
	switch os.Args[index*2-1] {
	case "+":
		return org + float32(targ)*anivalue
	case "*":
		return org * float32(targ) * anivalue
	case "/":
		return org / (float32(targ) * anivalue)
	case "^":
		return float32(math.Pow(float64(org), targ*float64(anivalue)))
	case "sin+":
		return 255 * float32(math.Sin((float64(org)+targ*float64(anivalue))/255*math.Pi*2))
	case "sin*":
		return 255 * float32(math.Sin((float64(org)*targ*float64(anivalue))/255*math.Pi*2))
	case "cos+":
		return 255 * float32(math.Cos((float64(org)+targ*float64(anivalue))/255*math.Pi*2))
	case "cos*":
		return 255 * float32(math.Cos((float64(org)*targ*float64(anivalue))/255*math.Pi*2))
	case "tan+":
		return 255 * float32(math.Tan((float64(org)+targ*float64(anivalue))/255*math.Pi*2))
	case "tan*":
		return 255 * float32(math.Tan((float64(org)*targ*float64(anivalue))/255*math.Pi*2))
	}
	return org
}

// PointArea returns the area on the unit sphere for the triangle defined by the
// given points.
//
// This method is based on l'Huilier's theorem,
//
//   tan(E/4) = sqrt(tan(s/2) tan((s-a)/2) tan((s-b)/2) tan((s-c)/2))
//
// where E is the spherical excess of the triangle (i.e. its area),
//       a, b, c are the side lengths, and
//       s is the semiperimeter (a + b + c) / 2.
//
// The only significant source of error using l'Huilier's method is the
// cancellation error of the terms (s-a), (s-b), (s-c). This leads to a
// *relative* error of about 1e-16 * s / min(s-a, s-b, s-c). This compares
// to a relative error of about 1e-15 / E using Girard's formula, where E is
// the true area of the triangle. Girard's formula can be even worse than
// this for very small triangles, e.g. a triangle with a true area of 1e-30
// might evaluate to 1e-5.
//
// So, we prefer l'Huilier's formula unless dmin < s * (0.1 * E), where
// dmin = min(s-a, s-b, s-c). This basically includes all triangles
// except for extremely long and skinny ones.
//
// Since we don't know E, we would like a conservative upper bound on
// the triangle area in terms of s and dmin. It's possible to show that
// E <= k1 * s * sqrt(s * dmin), where k1 = 2*sqrt(3)/Pi (about 1).
// Using this, it's easy to show that we should always use l'Huilier's
// method if dmin >= k2 * s^5, where k2 is about 1e-2. Furthermore,
// if dmin < k2 * s^5, the triangle area is at most k3 * s^4, where
// k3 is about 0.1. Since the best case error using Girard's formula
// is about 1e-15, this means that we shouldn't even consider it unless
// s >= 3e-4 or so.
func PointArea(a, b, c Point) float64 {
	sa := float64(b.Angle(c.Vector))
	sb := float64(c.Angle(a.Vector))
	sc := float64(a.Angle(b.Vector))
	s := 0.5 * (sa + sb + sc)
	if s >= 3e-4 {
		// Consider whether Girard's formula might be more accurate.
		dmin := s - math.Max(sa, math.Max(sb, sc))
		if dmin < 1e-2*s*s*s*s*s {
			// This triangle is skinny enough to use Girard's formula.
			ab := a.PointCross(b)
			bc := b.PointCross(c)
			ac := a.PointCross(c)
			area := math.Max(0.0, float64(ab.Angle(ac.Vector)-ab.Angle(bc.Vector)+bc.Angle(ac.Vector)))

			if dmin < s*0.1*area {
				return area
			}
		}
	}

	// Use l'Huilier's formula.
	return 4 * math.Atan(math.Sqrt(math.Max(0.0, math.Tan(0.5*s)*math.Tan(0.5*(s-sa))*
		math.Tan(0.5*(s-sb))*math.Tan(0.5*(s-sc)))))
}

// Angle returns the angle between great circles defined by three points.
//
// Coordinates may be right ascensions and declinations or longitudes and
// latitudes.  If r1, d1, r2, d2 defines one line and r2, d2, r3, d3 defines
// another, the result is the angle between the two lines.
//
// Algorithm by Meeus.
func Angle(r1, d1, r2, d2, r3, d3 float64) float64 {
	sd2, cd2 := math.Sincos(d2)
	sr21, cr21 := math.Sincos(r2 - r1)
	sr32, cr32 := math.Sincos(r3 - r2)
	C1 := math.Atan2(sr21, cd2*math.Tan(d1)-sd2*cr21)
	C2 := math.Atan2(sr32, cd2*math.Tan(d3)-sd2*cr32)
	return C1 + C2
}

// ConvertToNorthingEasting converts latitude and longitude to
// Ordnance Survey grid reference northing and easting.
// It returns a struct containing easting and northing coordinates as float64 type
// or an error if the arguments passed in are out of bounds
func ConvertToNorthingEasting(lat, lon float64) (*OsGrid, error) {
	o := OsGrid{}

	// validate input
	if lat < -90 || lat > 90 {
		return &o, ErrInvalidLat
	}

	if lon < -180 || lon > 180 {
		return &o, ErrInvalidLon
	}

	φ := toRadians(lat)
	λ := toRadians(lon)

	cosφ := math.Cos(φ)
	sinφ := math.Sin(φ)
	ν := a * f0 / math.Sqrt(1-e2*sinφ*sinφ)
	ρ := a * f0 * (1 - e2) / math.Pow(1-e2*sinφ*sinφ, 1.5)
	η2 := ν/ρ - 1

	Ma := (1 + n + (5/4)*n2 + (5/4)*n3) * (φ - φ0)
	Mb := (3*n + 3*n*n + (21/8)*n3) * math.Sin(φ-φ0) * math.Cos(φ+φ0)
	Mc := ((15/8)*n2 + (15/8)*n3) * math.Sin(2*(φ-φ0)) * math.Cos(2*(φ+φ0))
	Md := (35 / 24) * n3 * math.Sin(3*(φ-φ0)) * math.Cos(3*(φ+φ0))
	M := b * f0 * (Ma - Mb + Mc - Md)

	cos3φ := cosφ * cosφ * cosφ
	cos5φ := cos3φ * cosφ * cosφ
	tan2φ := math.Tan(φ) * math.Tan(φ)
	tan4φ := tan2φ * tan2φ

	I := M + n0
	II := (ν / 2) * sinφ * cosφ
	III := (ν / 24) * sinφ * cos3φ * (5 - tan2φ + 9*η2)
	IIIA := (ν / 720) * sinφ * cos5φ * (61 - 58*tan2φ + tan4φ)
	IV := ν * cosφ
	V := (ν / 6) * cos3φ * (ν/ρ - tan2φ)
	VI := (ν / 120) * cos5φ * (5 - 18*tan2φ + tan4φ + 14*η2 - 58*tan2φ*η2)

	Δλ := λ - λ0
	Δλ2 := Δλ * Δλ
	Δλ3 := Δλ2 * Δλ
	Δλ4 := Δλ3 * Δλ
	Δλ5 := Δλ4 * Δλ
	Δλ6 := Δλ5 * Δλ

	northingVal := I + II*Δλ2 + III*Δλ4 + IIIA*Δλ6
	northingVal, _ = strconv.ParseFloat(fmt.Sprintf("%.3f", northingVal), 64) // truncate after 3 decimal positions
	o.Northing = northingVal

	eastingVal := e0 + IV*Δλ + V*Δλ3 + VI*Δλ5
	eastingVal, _ = strconv.ParseFloat(fmt.Sprintf("%.3f", eastingVal), 64) // truncate after 3 decimal positions
	o.Easting = eastingVal

	return &o, nil
}

func TestAreaAndCentroid(t *testing.T) {
	if got := north_hemi.Area(); math.Abs(got-2*math.Pi) > 1e-15 {
		t.Errorf("%v.Area() == %v, want %v", north_hemi, got, 2*math.Pi)
	}
	if got := east_hemi.Area(); got > 2*math.Pi+1e-12 {
		t.Errorf("%v.Area() > %v, want <= %v", east_hemi, got, 2*math.Pi+1e-12)
	}
	if got := east_hemi.Area(); got < 2*math.Pi-1e-12 {
		t.Errorf("%v.Area() > %v, want >= %v", east_hemi, got, 2*math.Pi-1e-12)
	}

	// Construct spherical caps of random height, and approximate their
	// boundary with closely spaced vertices. Then check that the area and
	// centroid are correct.
	const kMaxDist = 1e-6
	for i := 0; i < 100; i++ {
		// Choose a coordinate frame for the spherical cap.
		x, y, z := randomFrame()
		// Given two points at latitude phi and whose longitudes differ
		// by dtheta, the geodesic between the two points has a maximum
		// latitude of atan(tan(phi) / cos(dtheta/2)). This can be
		// derived by positioning the two points at (-dtheta/2, phi)
		// and (dtheta/2, phi).
		//
		// We want to position the vertices close enough together so
		// that their maximum distance from the boundary of the
		// spherical cap is kMaxDist. Thus we want
		// fabs(atan(tan(phi) / cos(dtheta/2)) - phi) <= kMaxDist.
		height := 2 * rand.Float64()
		phi := math.Asin(1 - height)
		max_dtheta := 2 * math.Acos(math.Tan(math.Abs(phi))/math.Tan(math.Abs(phi)+kMaxDist))
		max_dtheta = math.Min(math.Pi, max_dtheta)
		var vertices []Point
		for theta := 0.0; theta < 2*math.Pi; theta += rand.Float64() * max_dtheta {
			a := x.Mul(math.Cos(theta) * math.Cos(phi))
			b := y.Mul(math.Sin(theta) * math.Cos(phi))
			c := z.Mul(math.Sin(phi))
			vertices = append(vertices, Point{a.Add(b).Add(c)})
		}
		loop := NewLoopFromPath(vertices)
		area := loop.Area()
		centroid := loop.Centroid()
		expectedArea := 2 * math.Pi * height
		if got := math.Abs(area - expectedArea); got > 2*math.Pi*kMaxDist {
			t.Errorf("%v > %v", got, 2*math.Pi*kMaxDist)
		}
		if got := math.Abs(area - expectedArea); got < 0.01*kMaxDist {
			t.Errorf("%v < %v", got, 0.01*kMaxDist)
		}
		expectedCentroid := z.Mul(expectedArea * (1 - 0.5*height))
		if got := centroid.Sub(expectedCentroid).Norm(); got > 2*kMaxDist {
			t.Errorf("(%v - %v).Norm() == %v, want <= %v", centroid,
				expectedCentroid, got, 2*kMaxDist)
		}
	}
}