// Position precesses equatorial coordinates from one epoch to another,
// including proper motions.
//
// If proper motions are not to be considered or are not applicable, pass 0, 0
// for mα, mδ
//
// Both eqFrom and eqTo must be non-nil, although they may point to the same
// struct.  EqTo is returned for convenience.
func Position(eqFrom, eqTo *coord.Equatorial, epochFrom, epochTo float64, mα base.HourAngle, mδ base.Angle) *coord.Equatorial {
	p := NewPrecessor(epochFrom, epochTo)
	t := epochTo - epochFrom
	eqTo.RA = eqFrom.RA + mα.Rad()*t
	eqTo.Dec = eqFrom.Dec + mδ.Rad()*t
	return p.Precess(eqTo, eqTo)
}

// ApproxPosition uses ApproxAnnualPrecession to compute a simple and quick
// precession while still considering proper motion.
//
// Both eqFrom and eqTo must be non-nil, although they may point to the same
// struct.  EqTo is returned for convenience.
func ApproxPosition(eqFrom, eqTo *coord.Equatorial, epochFrom, epochTo float64, mα base.HourAngle, mδ base.Angle) *coord.Equatorial {
	Δα, Δδ := ApproxAnnualPrecession(eqFrom, epochFrom, epochTo)
	dy := epochTo - epochFrom
	eqTo.RA = eqFrom.RA + (Δα+mα).Rad()*dy
	eqTo.Dec = eqFrom.Dec + (Δδ+mδ).Rad()*dy
	return eqTo
}

// ApproxPosition uses ApproxAnnualPrecession to compute a simple and quick
// precession while still considering proper motion.
//
// Both eqFrom and eqTo must be non-nil, although they may point to the same
// struct.  EqTo is returned for convenience.
func ApproxPosition(eqFrom, eqTo *coord.Equatorial, epochFrom, epochTo float64, mα unit.HourAngle, mδ unit.Angle) *coord.Equatorial {
	Δα, Δδ := ApproxAnnualPrecession(eqFrom, epochFrom, epochTo)
	dy := epochTo - epochFrom
	eqTo.RA = eqFrom.RA.Add((Δα + mα).Mul(dy))
	eqTo.Dec = eqFrom.Dec + (Δδ + mδ).Mul(dy)
	return eqTo
}

// Position precesses equatorial coordinates from one epoch to another,
// including proper motions.
//
// If proper motions are not to be considered or are not applicable, pass 0, 0
// for mα, mδ
//
// Both eqFrom and eqTo must be non-nil, although they may point to the same
// struct.  EqTo is returned for convenience.
func Position(eqFrom, eqTo *coord.Equatorial, epochFrom, epochTo float64, mα unit.HourAngle, mδ unit.Angle) *coord.Equatorial {
	p := NewPrecessor(epochFrom, epochTo)
	t := epochTo - epochFrom
	eqTo.RA = unit.RAFromRad(eqFrom.RA.Rad() + mα.Rad()*t)
	eqTo.Dec = eqFrom.Dec + mδ*unit.Angle(t)
	return p.Precess(eqTo, eqTo)
}

// Position computes the apparent position of an object.
//
// Position is computed for equatorial coordinates in eqFrom, considering
// proper motion, precession, nutation, and aberration.  Result is in
// eqTo.  EqFrom and eqTo must be non-nil, but may point to the same struct.
func Position(eqFrom, eqTo *coord.Equatorial, epochFrom, epochTo float64, mα sexa.HourAngle, mδ sexa.Angle) *coord.Equatorial {
	precess.Position(eqFrom, eqTo, epochFrom, epochTo, mα, mδ)
	jd := base.JulianYearToJDE(epochTo)
	Δα1, Δδ1 := Nutation(eqTo.RA, eqTo.Dec, jd)
	Δα2, Δδ2 := Aberration(eqTo.RA, eqTo.Dec, jd)
	eqTo.RA += Δα1 + Δα2
	eqTo.Dec += Δδ1 + Δδ2
	return eqTo
}

// Precess precesses coordinates eqFrom, leaving result in eqTo.
//
// The same struct may be used for eqFrom and eqTo.
// EqTo is returned for convenience.
func (p *Precessor) Precess(eqFrom, eqTo *coord.Equatorial) *coord.Equatorial {
	// (21.4) p. 134
	sδ, cδ := math.Sincos(eqFrom.Dec)
	sαζ, cαζ := math.Sincos(eqFrom.RA + p.ζ)
	A := cδ * sαζ
	B := p.cθ*cδ*cαζ - p.sθ*sδ
	C := p.sθ*cδ*cαζ + p.cθ*sδ
	eqTo.RA = math.Atan2(A, B) + p.z
	if C < base.CosSmallAngle {
		eqTo.Dec = math.Asin(C)
	} else {
		eqTo.Dec = math.Acos(math.Hypot(A, B)) // near pole
	}
	return eqTo
}

// ProperMotion3D takes the 3D equatorial coordinates of an object
// at one epoch and computes its coordinates at a new epoch, considering
// proper motion and radial velocity.
//
// Radial distance (r) must be in parsecs, radial velocitiy (mr) in
// parsecs per year.
//
// Both eqFrom and eqTo must be non-nil, although they may point to the same
// struct.  EqTo is returned for convenience.
func ProperMotion3D(eqFrom, eqTo *coord.Equatorial, epochFrom, epochTo, r, mr float64, mα base.HourAngle, mδ base.Angle) *coord.Equatorial {
	sα, cα := math.Sincos(eqFrom.RA)
	sδ, cδ := math.Sincos(eqFrom.Dec)
	x := r * cδ * cα
	y := r * cδ * sα
	z := r * sδ
	mrr := mr / r
	zmδ := z * mδ.Rad()
	mx := x*mrr - zmδ*cα - y*mα.Rad()
	my := y*mrr - zmδ*sα + x*mα.Rad()
	mz := z*mrr + r*mδ.Rad()*cδ
	t := epochTo - epochFrom
	xp := x + t*mx
	yp := y + t*my
	zp := z + t*mz
	eqTo.RA = math.Atan2(yp, xp)
	eqTo.Dec = math.Atan2(zp, math.Hypot(xp, yp))
	return eqTo
}

// PositionRonVondrak computes the apparent position of an object using
// the Ron-Vondrák expression for aberration.
//
// Position is computed for equatorial coordinates in eqFrom, considering
// proper motion, aberration, precession, and nutation.  Result is in
// eqTo.  EqFrom and eqTo must be non-nil, but may point to the same struct.
//
// Note the Ron-Vondrák expression is only valid for the epoch J2000.
// EqFrom must be coordinates at epoch J2000.
func PositionRonVondrak(eqFrom, eqTo *coord.Equatorial, epochTo float64, mα sexa.HourAngle, mδ sexa.Angle) *coord.Equatorial {
	t := epochTo - 2000
	eqTo.RA = eqFrom.RA + mα.Rad()*t
	eqTo.Dec = eqFrom.Dec + mδ.Rad()*t
	jd := base.JulianYearToJDE(epochTo)
	Δα, Δδ := AberrationRonVondrak(eqTo.RA, eqTo.Dec, jd)
	eqTo.RA += Δα
	eqTo.Dec += Δδ
	precess.Position(eqTo, eqTo, 2000, epochTo, 0, 0)
	Δα1, Δδ1 := Nutation(eqTo.RA, eqTo.Dec, jd)
	eqTo.RA += Δα1
	eqTo.Dec += Δδ1
	return eqTo
}

// PositionRonVondrak computes the apparent position of an object using
// the Ron-Vondrák expression for aberration.
//
// Position is computed for equatorial coordinates in eqFrom, considering
// proper motion, aberration, precession, and nutation.  Result is in
// eqTo.  EqFrom and eqTo must be non-nil, but may point to the same struct.
//
// Note the Ron-Vondrák expression is only valid for the epoch J2000.
// EqFrom must be coordinates at epoch J2000.
func PositionRonVondrak(eqFrom, eqTo *coord.Equatorial, epochTo float64, mα unit.HourAngle, mδ unit.Angle) *coord.Equatorial {
	t := epochTo - 2000
	eqTo.RA = eqFrom.RA.Add(mα.Mul(t))
	eqTo.Dec = eqFrom.Dec + mδ.Mul(t)
	jd := base.JulianYearToJDE(epochTo)
	Δα, Δδ := AberrationRonVondrak(eqTo.RA, eqTo.Dec, jd)
	eqTo.RA = eqTo.RA.Add(Δα)
	eqTo.Dec += Δδ
	precess.Position(eqTo, eqTo, 2000, epochTo, 0, 0)
	Δα1, Δδ1 := Nutation(eqTo.RA, eqTo.Dec, jd)
	eqTo.RA = eqTo.RA.Add(Δα1)
	eqTo.Dec += Δδ1
	return eqTo
}