fix: 修正行星事件边界与留点计算

- 统一 UT 事件时刻与 TT 查询时刻的边界判断
- 将外行星留点搜索锚定到对应冲日周期
- 修正水星、金星合日、留、大距事件选择
- 统一七大行星视位置计算辅助逻辑
- 增加公开 Last/Next 边界和 JPL/NAOJ 基线回归测试

basic/event_boundary.go

@@ -0,0 +1,76 @@
+package basic
+
+import "math"
+
+const exactEventTolerance = 2.0 / 86400.0
+
+func sameEventJD(a, b float64) bool {
+	return math.Abs(a-b) <= exactEventTolerance
+}
+
+func sameEventUTQueryTT(eventUT, queryTT float64) bool {
+	return math.Abs(eventUTQueryTTDelta(eventUT, queryTT)) <= exactEventTolerance
+}
+
+func closestEventUTToQueryTT(queryTT, best float64, candidates ...float64) float64 {
+	bestAbs := math.Abs(eventUTQueryTTDelta(best, queryTT))
+	for _, candidate := range candidates {
+		candidateAbs := math.Abs(eventUTQueryTTDelta(candidate, queryTT))
+		if candidateAbs < bestAbs {
+			best = candidate
+			bestAbs = candidateAbs
+		}
+	}
+	return best
+}
+
+type phaseEventSearchFunc func(jde, degree float64, next uint8) float64
+type simpleEventSearchFunc func(jde float64) float64
+
+func inclusiveLastPhaseEvent(jde, degree float64, fn phaseEventSearchFunc) float64 {
+	last := fn(jde, degree, 0)
+	next := fn(jde, degree, 1)
+	if eventUTQueryBeforeOrEqual(next, jde) && eventUTQueryAfterOrEqual(next, jde) {
+		return next
+	}
+	if eventUTQueryBeforeOrEqual(last, jde) {
+		return last
+	}
+	return last
+}
+
+func inclusiveNextPhaseEvent(jde, degree float64, fn phaseEventSearchFunc) float64 {
+	last := fn(jde, degree, 0)
+	if eventUTQueryBeforeOrEqual(last, jde) && eventUTQueryAfterOrEqual(last, jde) {
+		return last
+	}
+	next := fn(jde, degree, 1)
+	if eventUTQueryAfterOrEqual(next, jde) {
+		return next
+	}
+	return next
+}
+
+func inclusiveLastSimpleEvent(jde float64, lastFn, nextFn simpleEventSearchFunc) float64 {
+	last := lastFn(jde)
+	next := nextFn(jde)
+	if eventUTQueryBeforeOrEqual(next, jde) && eventUTQueryAfterOrEqual(next, jde) {
+		return next
+	}
+	if eventUTQueryBeforeOrEqual(last, jde) {
+		return last
+	}
+	return last
+}
+
+func inclusiveNextSimpleEvent(jde float64, lastFn, nextFn simpleEventSearchFunc) float64 {
+	last := lastFn(jde)
+	if eventUTQueryBeforeOrEqual(last, jde) && eventUTQueryAfterOrEqual(last, jde) {
+		return last
+	}
+	next := nextFn(jde)
+	if eventUTQueryAfterOrEqual(next, jde) {
+		return next
+	}
+	return next
+}

basic/event_refine.go

@@ -40,9 +40,7 @@ func eventZeroBracket(leftJD, leftVal, centerJD, centerVal, rightJD, rightVal fl
 	return 0, 0, 0, 0, false
 }
 
-// eventZeroRefine 细化 seed 附近的零点；若找不到可用括号区间，则退回旧的固定步长扫描。
-// eventZeroRefine refines a nearby zero crossing and falls back to the legacy
-// fixed-step scan when no usable bracket is found.
+// eventZeroRefine 细化 seed 附近的零点；无可用括号区间时退回固定步长扫描。
 func eventZeroRefine(seed, halfWindow, step float64, fn func(float64) float64) float64 {
 	leftJD := seed - halfWindow
 	centerJD := seed

basic/inner_event_window.go

@@ -0,0 +1,183 @@
+package basic
+
+import "math"
+
+const innerEventEpsilon = 4.0 / 86400.0
+
+func eventQueryTTAsUT(queryTT float64) float64 {
+	return TD2UT(queryTT, false)
+}
+
+func eventUTQueryTTDelta(eventUT, queryTT float64) float64 {
+	return eventUT - eventQueryTTAsUT(queryTT)
+}
+
+func eventUTQueryBeforeOrEqual(eventUT, queryTT float64) bool {
+	return eventUTQueryTTDelta(eventUT, queryTT) <= innerEventEpsilon
+}
+
+func eventUTQueryAfterOrEqual(eventUT, queryTT float64) bool {
+	return eventUTQueryTTDelta(eventUT, queryTT) >= -innerEventEpsilon
+}
+
+func eventUTNextQueryTT(eventUT float64) float64 {
+	return TD2UT(eventUT, true) + 1.0
+}
+
+func eventUTLastQueryTT(eventUT float64) float64 {
+	return TD2UT(eventUT, true) - 1.0
+}
+
+func innerNextCycleOffset(delta, period float64) float64 {
+	if delta <= 0 {
+		return -delta * period / 360.0
+	}
+	return (360.0 - delta) * period / 360.0
+}
+
+func innerLastCycleOffset(delta, period float64) float64 {
+	if delta >= 0 {
+		return delta * period / 360.0
+	}
+	return (360.0 + delta) * period / 360.0
+}
+
+func clampFloat64(v, min, max float64) float64 {
+	if v < min {
+		return min
+	}
+	if v > max {
+		return max
+	}
+	return v
+}
+
+func scanWindowForMinAbs(start, end, step float64, fn func(float64) float64) float64 {
+	if end < start {
+		start, end = end, start
+	}
+	if step <= 0 || end == start {
+		return start
+	}
+	bestJD := start
+	bestAbs := math.Abs(fn(start))
+	for jd := start + step; jd < end; jd += step {
+		candidateAbs := math.Abs(fn(jd))
+		if candidateAbs < bestAbs {
+			bestAbs = candidateAbs
+			bestJD = jd
+		}
+	}
+	endAbs := math.Abs(fn(end))
+	if endAbs < bestAbs {
+		return end
+	}
+	return bestJD
+}
+
+func scanWindowForMax(start, end, step float64, fn func(float64) float64) float64 {
+	if end < start {
+		start, end = end, start
+	}
+	if step <= 0 || end == start {
+		return start
+	}
+	bestJD := start
+	bestVal := fn(start)
+	for jd := start + step; jd < end; jd += step {
+		candidateVal := fn(jd)
+		if candidateVal > bestVal {
+			bestVal = candidateVal
+			bestJD = jd
+		}
+	}
+	endVal := fn(end)
+	if endVal > bestVal {
+		return end
+	}
+	return bestJD
+}
+
+func boundedEventZeroRefine(seed, start, end, halfWindow, step float64, fn func(float64) float64) float64 {
+	if end < start {
+		start, end = end, start
+	}
+	if end <= start {
+		return start
+	}
+	maxHalfWindow := (end - start) / 2
+	if halfWindow > maxHalfWindow {
+		halfWindow = maxHalfWindow
+	}
+	if halfWindow <= 0 {
+		return clampFloat64(seed, start, end)
+	}
+	seed = clampFloat64(seed, start+halfWindow, end-halfWindow)
+	return eventZeroRefine(seed, halfWindow, step, fn)
+}
+
+func zeroEventInWindow(start, end, coarseStep, halfWindow, refineStep float64, coarseFn, exactFn func(float64) float64) float64 {
+	if end < start {
+		start, end = end, start
+	}
+	if end <= start {
+		return start
+	}
+	rangeDays := end - start
+	if coarseStep <= 0 || coarseStep > rangeDays {
+		coarseStep = rangeDays / 6.0
+	}
+	if coarseStep < 0.5 {
+		coarseStep = 0.5
+	}
+	if refineStep <= 0 {
+		refineStep = 0.5 / 86400.0
+	}
+	if halfWindow <= 0 {
+		halfWindow = coarseStep
+	}
+	guess := scanWindowForMinAbs(start, end, coarseStep, coarseFn)
+	return boundedEventZeroRefine(guess, start, end, halfWindow, refineStep, exactFn)
+}
+
+func maximizeInWindow(start, end, coarseStep float64, coarseFn, exactFn func(float64) float64) float64 {
+	if end < start {
+		start, end = end, start
+	}
+	if end <= start {
+		return start
+	}
+	rangeDays := end - start
+	if coarseStep <= 0 || coarseStep > rangeDays {
+		coarseStep = rangeDays / 6.0
+	}
+	if coarseStep < 0.5 {
+		coarseStep = 0.5
+	}
+	guess := scanWindowForMax(start, end, coarseStep, coarseFn)
+	left := clampFloat64(guess-coarseStep, start, end)
+	right := clampFloat64(guess+coarseStep, start, end)
+	if right-left <= innerEventEpsilon {
+		return guess
+	}
+	for i := 0; i < 20; i++ {
+		third := (right - left) / 3.0
+		leftThird := left + third
+		rightThird := right - third
+		if exactFn(leftThird) <= exactFn(rightThird) {
+			left = leftThird
+			continue
+		}
+		right = rightThird
+	}
+	bestJD := guess
+	bestVal := exactFn(bestJD)
+	for _, jd := range []float64{left, (left + right) / 2.0, right} {
+		candidateVal := exactFn(jd)
+		if candidateVal > bestVal {
+			bestVal = candidateVal
+			bestJD = jd
+		}
+	}
+	return bestJD
+}

basic/inner_planet_event_boundary_test.go

@@ -0,0 +1,45 @@
+package basic
+
+import "testing"
+
+func TestInnerPlanetExactEventBoundaryIncludesCurrent(t *testing.T) {
+	cases := []struct {
+		name   string
+		seed   float64
+		lastFn func(float64) float64
+		nextFn func(float64) float64
+	}{
+		{name: "MercuryConjunction", seed: NextMercuryConjunction(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastMercuryConjunction, nextFn: NextMercuryConjunction},
+		{name: "MercuryInferior", seed: NextMercuryInferiorConjunctionInclusive(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastMercuryInferiorConjunctionInclusive, nextFn: NextMercuryInferiorConjunctionInclusive},
+		{name: "MercurySuperior", seed: NextMercurySuperiorConjunctionInclusive(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastMercurySuperiorConjunctionInclusive, nextFn: NextMercurySuperiorConjunctionInclusive},
+		{name: "MercuryRetrograde", seed: NextMercuryRetrogradeInclusive(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastMercuryRetrogradeInclusive, nextFn: NextMercuryRetrogradeInclusive},
+		{name: "MercuryP2R", seed: NextMercuryProgradeToRetrogradeInclusive(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastMercuryProgradeToRetrogradeInclusive, nextFn: NextMercuryProgradeToRetrogradeInclusive},
+		{name: "MercuryR2P", seed: NextMercuryRetrogradeToProgradeInclusive(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastMercuryRetrogradeToProgradeInclusive, nextFn: NextMercuryRetrogradeToProgradeInclusive},
+		{name: "MercuryGreatestElongation", seed: NextMercuryGreatestElongationInclusive(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastMercuryGreatestElongationInclusive, nextFn: NextMercuryGreatestElongationInclusive},
+		{name: "MercuryEastElongation", seed: NextMercuryGreatestElongationEastInclusive(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastMercuryGreatestElongationEastInclusive, nextFn: NextMercuryGreatestElongationEastInclusive},
+		{name: "MercuryWestElongation", seed: NextMercuryGreatestElongationWestInclusive(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastMercuryGreatestElongationWestInclusive, nextFn: NextMercuryGreatestElongationWestInclusive},
+		{name: "VenusConjunction", seed: NextVenusConjunction(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastVenusConjunction, nextFn: NextVenusConjunction},
+		{name: "VenusInferior", seed: NextVenusInferiorConjunctionInclusive(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastVenusInferiorConjunctionInclusive, nextFn: NextVenusInferiorConjunctionInclusive},
+		{name: "VenusSuperior", seed: NextVenusSuperiorConjunctionInclusive(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastVenusSuperiorConjunctionInclusive, nextFn: NextVenusSuperiorConjunctionInclusive},
+		{name: "VenusRetrograde", seed: NextVenusRetrogradeInclusive(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastVenusRetrogradeInclusive, nextFn: NextVenusRetrogradeInclusive},
+		{name: "VenusP2R", seed: NextVenusProgradeToRetrogradeInclusive(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastVenusProgradeToRetrogradeInclusive, nextFn: NextVenusProgradeToRetrogradeInclusive},
+		{name: "VenusR2P", seed: NextVenusRetrogradeToProgradeInclusive(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastVenusRetrogradeToProgradeInclusive, nextFn: NextVenusRetrogradeToProgradeInclusive},
+		{name: "VenusGreatestElongation", seed: NextVenusGreatestElongationInclusive(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastVenusGreatestElongationInclusive, nextFn: NextVenusGreatestElongationInclusive},
+		{name: "VenusEastElongation", seed: NextVenusGreatestElongationEastInclusive(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastVenusGreatestElongationEastInclusive, nextFn: NextVenusGreatestElongationEastInclusive},
+		{name: "VenusWestElongation", seed: NextVenusGreatestElongationWestInclusive(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastVenusGreatestElongationWestInclusive, nextFn: NextVenusGreatestElongationWestInclusive},
+	}
+
+	for _, tc := range cases {
+		t.Run(tc.name, func(t *testing.T) {
+			queryTT := TD2UT(tc.seed, true)
+			last := tc.lastFn(queryTT)
+			next := tc.nextFn(queryTT)
+			if !sameEventJD(last, tc.seed) {
+				t.Fatalf("last exact boundary mismatch: got %.12f want %.12f", last, tc.seed)
+			}
+			if !sameEventJD(next, tc.seed) {
+				t.Fatalf("next exact boundary mismatch: got %.12f want %.12f", next, tc.seed)
+			}
+		})
+	}
+}

basic/inner_planet_truth_test.go

@@ -0,0 +1,165 @@
+package basic
+
+import (
+	"encoding/json"
+	"os"
+	"strings"
+	"testing"
+	"time"
+)
+
+type innerBaselineFile struct {
+	Events []innerBaselineEvent `json:"events"`
+}
+
+type innerBaselineEvent struct {
+	Planet          string `json:"planet"`
+	Kind            string `json:"kind"`
+	NAOJHintJST     string `json:"naoj_hint_jst"`
+	Precision       string `json:"precision"`
+	CandidateJST    string `json:"candidate_jst"`
+	VerifiedJST     string `json:"verified_jst"`
+	CandidateSource string `json:"candidate_source"`
+}
+
+func loadInnerBaseline(t *testing.T) innerBaselineFile {
+	t.Helper()
+
+	paths := [][]string{
+		{
+			"testdata/jpl_inner_event_baseline.json",
+			"basic/testdata/jpl_inner_event_baseline.json",
+		},
+		{
+			"testdata/jpl_inner_event_baseline_21c_sample.json",
+			"basic/testdata/jpl_inner_event_baseline_21c_sample.json",
+		},
+		{
+			"testdata/jpl_inner_event_baseline_20c_sample.json",
+			"basic/testdata/jpl_inner_event_baseline_20c_sample.json",
+		},
+		{
+			"testdata/jpl_inner_event_baseline_22c_sample.json",
+			"basic/testdata/jpl_inner_event_baseline_22c_sample.json",
+		},
+	}
+	var merged innerBaselineFile
+	for index, candidates := range paths {
+		var (
+			data []byte
+			err  error
+		)
+		for _, path := range candidates {
+			data, err = os.ReadFile(path)
+			if err == nil {
+				var baseline innerBaselineFile
+				if err := json.Unmarshal(data, &baseline); err != nil {
+					t.Fatal(err)
+				}
+				merged.Events = append(merged.Events, baseline.Events...)
+				break
+			}
+		}
+		if err != nil && index == 0 {
+			t.Fatal(err)
+		}
+	}
+	if len(merged.Events) == 0 {
+		t.Fatal("empty inner baseline file")
+	}
+	return merged
+}
+
+func parseInnerBaselineTime(t *testing.T, value string) time.Time {
+	t.Helper()
+	loc := time.FixedZone("JST", 9*3600)
+	layouts := []string{
+		"2006-01-02 15:04:05 MST",
+		"2006-01-02 15:04 MST",
+		"2006-01-02 15:04:05",
+		"2006-01-02 15:04",
+	}
+	var err error
+	for _, layout := range layouts {
+		when, parseErr := time.ParseInLocation(layout, value, loc)
+		if parseErr == nil {
+			return when
+		}
+		err = parseErr
+	}
+	t.Fatalf("parse baseline time %q: %v", value, err)
+	return time.Time{}
+}
+
+func innerBaselineTolerance(event innerBaselineEvent) time.Duration {
+	switch event.Kind {
+	case "IC", "SC", "P2R", "R2P":
+		return 2 * time.Minute
+	case "GEE", "GEW":
+		return 90 * time.Minute
+	default:
+		return 2 * time.Minute
+	}
+}
+
+func innerEventFuncs(t *testing.T, event innerBaselineEvent) (func(float64) float64, func(float64) float64) {
+	t.Helper()
+	switch event.Planet + ":" + event.Kind {
+	case "Mercury:IC":
+		return LastMercuryInferiorConjunctionInclusive, NextMercuryInferiorConjunctionInclusive
+	case "Mercury:SC":
+		return LastMercurySuperiorConjunctionInclusive, NextMercurySuperiorConjunctionInclusive
+	case "Mercury:P2R":
+		return LastMercuryProgradeToRetrogradeInclusive, NextMercuryProgradeToRetrogradeInclusive
+	case "Mercury:R2P":
+		return LastMercuryRetrogradeToProgradeInclusive, NextMercuryRetrogradeToProgradeInclusive
+	case "Mercury:GEE":
+		return LastMercuryGreatestElongationEastInclusive, NextMercuryGreatestElongationEastInclusive
+	case "Mercury:GEW":
+		return LastMercuryGreatestElongationWestInclusive, NextMercuryGreatestElongationWestInclusive
+	case "Venus:IC":
+		return LastVenusInferiorConjunctionInclusive, NextVenusInferiorConjunctionInclusive
+	case "Venus:SC":
+		return LastVenusSuperiorConjunctionInclusive, NextVenusSuperiorConjunctionInclusive
+	case "Venus:P2R":
+		return LastVenusProgradeToRetrogradeInclusive, NextVenusProgradeToRetrogradeInclusive
+	case "Venus:R2P":
+		return LastVenusRetrogradeToProgradeInclusive, NextVenusRetrogradeToProgradeInclusive
+	case "Venus:GEE":
+		return LastVenusGreatestElongationEastInclusive, NextVenusGreatestElongationEastInclusive
+	case "Venus:GEW":
+		return LastVenusGreatestElongationWestInclusive, NextVenusGreatestElongationWestInclusive
+	default:
+		t.Fatalf("unsupported event %s:%s", event.Planet, event.Kind)
+		return nil, nil
+	}
+}
+
+func assertInnerBaselineEvent(t *testing.T, event innerBaselineEvent, lastFn, nextFn func(float64) float64) {
+	t.Helper()
+	when := parseInnerBaselineTime(t, event.VerifiedJST)
+	before := when.Add(-24 * time.Hour)
+	after := when.Add(24 * time.Hour)
+	next := JDE2DateByZone(nextFn(toUTJD(before)), when.Location(), false)
+	last := JDE2DateByZone(lastFn(toUTJD(after)), when.Location(), false)
+	tolerance := innerBaselineTolerance(event)
+
+	if diff := next.Sub(when); diff < -tolerance || diff > tolerance {
+		t.Fatalf("%s %s next mismatch: got %s want %s tol=%s hint=%s candidate=%s via=%s", event.Planet, event.Kind, next, when, tolerance, event.NAOJHintJST, event.CandidateJST, event.CandidateSource)
+	}
+	if diff := last.Sub(when); diff < -tolerance || diff > tolerance {
+		t.Fatalf("%s %s last mismatch: got %s want %s tol=%s hint=%s candidate=%s via=%s", event.Planet, event.Kind, last, when, tolerance, event.NAOJHintJST, event.CandidateJST, event.CandidateSource)
+	}
+}
+
+func TestInnerPlanetTruthAgainstJPL(t *testing.T) {
+	baseline := loadInnerBaseline(t)
+	for _, event := range baseline.Events {
+		event := event
+		name := strings.Join([]string{event.Planet, event.Kind, event.VerifiedJST}, "_")
+		t.Run(name, func(t *testing.T) {
+			lastFn, nextFn := innerEventFuncs(t, event)
+			assertInnerBaselineEvent(t, event, lastFn, nextFn)
+		})
+	}
+}

basic/jupiter.go

@@ -87,53 +87,22 @@ func JupiterApparentRaDec(jd float64) (float64, float64) {
 }
 
 func EarthJupiterAway(jd float64) float64 {
-	x, y, z := AJupiterXYZ(jd)
-	to := math.Sqrt(x*x + y*y + z*z)
-	return to
+	return planetEarthAwayExplicitN(4, jd, -1)
 }
 
 func JupiterApparentLo(jd float64) float64 {
-	x, y, z := AJupiterXYZ(jd)
-	to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
-	x, y, z = AJupiterXYZ(jd - to)
-	lo := math.Atan2(y, x)
-	bo := math.Atan2(z, math.Sqrt(x*x+y*y))
-	lo = lo * 180 / math.Pi
-	bo = bo * 180 / math.Pi
-	lo = Limit360(lo)
-	//lo-=GXCLo(lo,bo,jd)/3600;
-	//bo+=GXCBo(lo,bo,jd);
-	lo += Nutation2000Bi(jd)
-	return lo
+	geo, _ := planetApparentGeocentricPositionN(4, jd, -1)
+	return geo.lo
 }
 
 func JupiterApparentBo(jd float64) float64 {
-	x, y, z := AJupiterXYZ(jd)
-	to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
-	x, y, z = AJupiterXYZ(jd - to)
-	//lo := math.Atan2(y, x)
-	bo := math.Atan2(z, math.Sqrt(x*x+y*y))
-	//lo = lo * 180 / math.Pi
-	bo = bo * 180 / math.Pi
-	//lo+=GXCLo(lo,bo,jd);
-	//bo+=GXCBo(lo,bo,jd)/3600;
-	//lo+=Nutation2000Bi(jd);
-	return bo
+	geo, _ := planetApparentGeocentricPositionN(4, jd, -1)
+	return geo.bo
 }
 
 func JupiterApparentLoBo(jd float64) (float64, float64) {
-	x, y, z := AJupiterXYZ(jd)
-	to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
-	x, y, z = AJupiterXYZ(jd - to)
-	lo := math.Atan2(y, x)
-	bo := math.Atan2(z, math.Sqrt(x*x+y*y))
-	lo = lo * 180 / math.Pi
-	bo = bo * 180 / math.Pi
-	lo = Limit360(lo)
-	//lo-=GXCLo(lo,bo,jd)/3600;
-	//bo+=GXCBo(lo,bo,jd);
-	lo += Nutation2000Bi(jd)
-	return lo, bo
+	geo, _ := planetApparentGeocentricPositionN(4, jd, -1)
+	return geo.lo, geo.bo
 }
 
 func JupiterMag(jd float64) float64 {

basic/jupiter_events.go

@@ -40,6 +40,28 @@ func jupiterSunLongitudeDeltaN(jde, degree float64, filter bool, n int) float64
 	return sub
 }
 
+func jupiterRADerivative(jde, delta float64) float64 {
+	sub := JupiterApparentRa(jde+delta) - JupiterApparentRa(jde-delta)
+	if sub > 180 {
+		sub -= 360
+	}
+	if sub < -180 {
+		sub += 360
+	}
+	return sub / (2 * delta)
+}
+
+func jupiterRADerivativeN(jde, delta float64, n int) float64 {
+	sub := JupiterApparentRaN(jde+delta, n) - JupiterApparentRaN(jde-delta, n)
+	if sub > 180 {
+		sub -= 360
+	}
+	if sub < -180 {
+		sub += 360
+	}
+	return sub / (2 * delta)
+}
+
 func jupiterConjunctionFull(jde, degree float64, next uint8) float64 {
 	//0=last 1=next
 	daysPerDegree := JUPITER_S_PERIOD / 360
@@ -94,113 +116,92 @@ func jupiterConjunction(jde, degree float64, next uint8) float64 {
 }
 
 func LastJupiterConjunction(jde float64) float64 {
-	return jupiterConjunction(jde, 0, 0)
+	return inclusiveLastPhaseEvent(jde, 0, jupiterConjunction)
 }
 
 func NextJupiterConjunction(jde float64) float64 {
-	return jupiterConjunction(jde, 0, 1)
+	return inclusiveNextPhaseEvent(jde, 0, jupiterConjunction)
 }
 
 func LastJupiterOpposition(jde float64) float64 {
-	return jupiterConjunction(jde, 180, 0)
+	return inclusiveLastPhaseEvent(jde, 180, jupiterConjunction)
 }
 
 func NextJupiterOpposition(jde float64) float64 {
-	return jupiterConjunction(jde, 180, 1)
+	return inclusiveNextPhaseEvent(jde, 180, jupiterConjunction)
 }
 
 func NextJupiterEasternQuadrature(jde float64) float64 {
-	return jupiterConjunction(jde, 90, 1)
+	return inclusiveNextPhaseEvent(jde, 90, jupiterConjunction)
 }
 
 func LastJupiterEasternQuadrature(jde float64) float64 {
-	return jupiterConjunction(jde, 90, 0)
+	return inclusiveLastPhaseEvent(jde, 90, jupiterConjunction)
 }
 
 func NextJupiterWesternQuadrature(jde float64) float64 {
-	return jupiterConjunction(jde, 270, 1)
+	return inclusiveNextPhaseEvent(jde, 270, jupiterConjunction)
 }
 
 func LastJupiterWesternQuadrature(jde float64) float64 {
-	return jupiterConjunction(jde, 270, 0)
+	return inclusiveLastPhaseEvent(jde, 270, jupiterConjunction)
 }
 
-func jupiterRetrograde(jde float64, searchBeforeOpposition bool) float64 {
-	//0=last 1=next
-	raRate := func(jde float64, delta float64) float64 {
-		sub := JupiterApparentRa(jde+delta) - JupiterApparentRa(jde-delta)
-		if sub > 180 {
-			sub -= 360
-		}
-		if sub < -180 {
-			sub += 360
-		}
-		return sub / (2 * delta)
-	}
-	jde = jupiterConjunctionFull(jde, 180, 1)
+func jupiterRetrogradeAroundOpposition(oppositionJD float64, searchBeforeOpposition bool) float64 {
+	oppositionTT := TD2UT(oppositionJD, true)
+	startTT := oppositionTT
+	endTT := oppositionTT
 	if searchBeforeOpposition {
-		jde -= 60
+		easternQuadratureUT := jupiterConjunction(oppositionTT, 90, 0)
+		startTT = TD2UT(easternQuadratureUT, true)
 	} else {
-		jde += 60
+		westernQuadratureUT := jupiterConjunction(oppositionTT, 270, 1)
+		endTT = TD2UT(westernQuadratureUT, true)
 	}
-	for {
-		currentRate := raRate(jde, 1.0/86400.0)
-		if math.Abs(currentRate) > 0.55 {
-			jde += 2
-			continue
-		}
-		break
-	}
-	estimateJD := jde
-	for {
-		prevJD := estimateJD
-		rateValue := raRate(prevJD, 2.0/86400.0)
-		rateSlope := (raRate(prevJD+15.0/86400.0, 2.0/86400.0) - raRate(prevJD-15.0/86400.0, 2.0/86400.0)) / (30.0 / 86400.0)
-		estimateJD = prevJD - rateValue/rateSlope
-		if math.Abs(estimateJD-prevJD) <= 30.0/86400.0 {
-			break
-		}
-	}
-	bestJD := eventZeroRefine(estimateJD, 15.0/86400.0, 0.5/86400.0, func(jd float64) float64 {
-		return raRate(jd, 0.5/86400.0)
+	bestJD := zeroEventInWindow(startTT, endTT, 2.0, 2.0, 30.0/86400.0, func(jd float64) float64 {
+		return jupiterRADerivativeN(jd, 1.0/86400.0, jupiterEventSearchN)
+	}, func(jd float64) float64 {
+		return jupiterRADerivative(jd, 0.5/86400.0)
 	})
 	return TD2UT(bestJD, false)
 }
 
 func NextJupiterRetrogradeToPrograde(jde float64) float64 {
-	date := jupiterRetrograde(jde, false)
-	if date < jde {
-		oppositionJD := jupiterConjunctionFull(jde, 180, 1)
-		return jupiterRetrograde(oppositionJD+10, false)
+	lastOppositionJD := jupiterConjunctionFull(jde, 180, 0)
+	date := jupiterRetrogradeAroundOpposition(lastOppositionJD, false)
+	if sameEventUTQueryTT(date, jde) || eventUTQueryAfterOrEqual(date, jde) {
+		return date
 	}
-	return date
+	nextOppositionJD := jupiterConjunctionFull(jde, 180, 1)
+	return jupiterRetrogradeAroundOpposition(nextOppositionJD, false)
 }
 
 func LastJupiterRetrogradeToPrograde(jde float64) float64 {
-	jde = jupiterConjunctionFull(jde, 180, 0) - 10
-	date := jupiterRetrograde(jde, false)
-	if date > jde {
-		oppositionJD := jupiterConjunctionFull(jde, 180, 0)
-		return jupiterRetrograde(oppositionJD-10, false)
+	lastOppositionJD := jupiterConjunctionFull(jde, 180, 0)
+	date := jupiterRetrogradeAroundOpposition(lastOppositionJD, false)
+	if sameEventUTQueryTT(date, jde) || eventUTQueryBeforeOrEqual(date, jde) {
+		return date
 	}
-	return date
+	previousOppositionJD := jupiterConjunctionFull(eventUTLastQueryTT(lastOppositionJD), 180, 0)
+	return jupiterRetrogradeAroundOpposition(previousOppositionJD, false)
 }
 
 func NextJupiterProgradeToRetrograde(jde float64) float64 {
-	date := jupiterRetrograde(jde, true)
-	if date < jde {
-		oppositionJD := jupiterConjunctionFull(jde, 180, 1)
-		return jupiterRetrograde(oppositionJD+10, true)
+	nextOppositionJD := jupiterConjunctionFull(jde, 180, 1)
+	date := jupiterRetrogradeAroundOpposition(nextOppositionJD, true)
+	if sameEventUTQueryTT(date, jde) || eventUTQueryAfterOrEqual(date, jde) {
+		return date
 	}
-	return date
+	followingOppositionJD := jupiterConjunctionFull(eventUTNextQueryTT(nextOppositionJD), 180, 1)
+	return jupiterRetrogradeAroundOpposition(followingOppositionJD, true)
 }
 
 func LastJupiterProgradeToRetrograde(jde float64) float64 {
-	jde = jupiterConjunctionFull(jde, 180, 0) - 10
-	date := jupiterRetrograde(jde, true)
-	if date > jde {
-		oppositionJD := jupiterConjunctionFull(jde, 180, 0)
-		return jupiterRetrograde(oppositionJD-10, true)
+	nextOppositionJD := jupiterConjunctionFull(jde, 180, 1)
+	date := jupiterRetrogradeAroundOpposition(nextOppositionJD, true)
+	if sameEventUTQueryTT(date, jde) || eventUTQueryBeforeOrEqual(date, jde) {
+		return date
 	}
-	return date
+	lastOppositionJD := jupiterConjunctionFull(jde, 180, 0)
+	return jupiterRetrogradeAroundOpposition(lastOppositionJD, true)
 }

basic/mars.go

@@ -87,72 +87,32 @@ func MarsApparentRaDec(jd float64) (float64, float64) {
 }
 
 func EarthMarsAway(jd float64) float64 {
-	x, y, z := AMarsXYZ(jd)
-	to := math.Sqrt(x*x + y*y + z*z)
-	return to
+	return planetEarthAwayExplicitN(3, jd, -1)
 }
 
 func MarsApparentLo(jd float64) float64 {
-	x, y, z := AMarsXYZ(jd)
-	to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
-	x, y, z = AMarsXYZ(jd - to)
-	lo := math.Atan2(y, x)
-	//bo := math.Atan2(z, math.Sqrt(x*x+y*y))
-	lo = lo * 180.0 / math.Pi
-	//bo = bo * 180 / math.Pi
-	lo = Limit360(lo) + Nutation2000Bi(jd)
-	//lo-=GXCLo(lo,bo,jd)/3600;
-	//bo+=GXCBo(lo,bo,jd);
-	return lo
+	geo, _ := planetApparentGeocentricPositionN(3, jd, -1)
+	return geo.lo
 }
 
 func MarsApparentBo(jd float64) float64 {
-	x, y, z := AMarsXYZ(jd)
-	to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
-	x, y, z = AMarsXYZ(jd - to)
-	//lo := math.Atan2(y, x)
-	bo := math.Atan2(z, math.Sqrt(x*x+y*y))
-	//lo = lo * 180 / math.Pi
-	bo = bo * 180 / math.Pi
-	//lo+=GXCLo(lo,bo,jd);
-	//bo+=GXCBo(lo,bo,jd)/3600;
-	//lo+=Nutation2000Bi(jd);
-	return bo
+	geo, _ := planetApparentGeocentricPositionN(3, jd, -1)
+	return geo.bo
 }
 
 func MarsApparentLoBo(jd float64) (float64, float64) {
-	x, y, z := AMarsXYZ(jd)
-	to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
-	x, y, z = AMarsXYZ(jd - to)
-	lo := math.Atan2(y, x)
-	bo := math.Atan2(z, math.Sqrt(x*x+y*y))
-	lo = lo * 180 / math.Pi
-	bo = bo * 180 / math.Pi
-	lo = Limit360(lo)
-	//lo -= GXCLo(lo, bo, jd) / 3600
-	//bo += GXCBo(lo, bo, jd)
-	lo += Nutation2000Bi(jd)
-	return lo, bo
+	geo, _ := planetApparentGeocentricPositionN(3, jd, -1)
+	return geo.lo, geo.bo
 }
 
 func MarsTrueLoBo(jd float64) (float64, float64) {
-	x, y, z := AMarsXYZ(jd)
-	to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
-	x, y, z = AMarsXYZ(jd - to)
-	lo := math.Atan2(y, x)
-	bo := math.Atan2(z, math.Sqrt(x*x+y*y))
-	lo = lo * 180 / math.Pi
-	bo = bo * 180 / math.Pi
-	lo = Limit360(lo)
-	return lo, bo
+	geo, _ := planetTrueGeocentricPositionN(3, jd, -1)
+	return geo.lo, geo.bo
 }
 
 func MarsTrueLo(jd float64) float64 {
-	x, y, _ := AMarsXYZ(jd)
-	lo := math.Atan2(y, x)
-	lo = lo * 180 / math.Pi
-	lo = Limit360(lo)
-	return lo
+	geo, _ := planetTrueGeocentricPositionN(3, jd, -1)
+	return geo.lo
 }
 
 func MarsMag(jd float64) float64 {

basic/mars_events.go

@@ -116,40 +116,39 @@ func marsConjunction(jde, degree float64, next uint8) float64 {
 }
 
 func LastMarsConjunction(jde float64) float64 {
-	return marsConjunction(jde, 0, 0)
+	return inclusiveLastPhaseEvent(jde, 0, marsConjunction)
 }
 
 func NextMarsConjunction(jde float64) float64 {
-	return marsConjunction(jde, 0, 1)
+	return inclusiveNextPhaseEvent(jde, 0, marsConjunction)
 }
 
 func LastMarsOpposition(jde float64) float64 {
-	return marsConjunction(jde, 180, 0)
+	return inclusiveLastPhaseEvent(jde, 180, marsConjunction)
 }
 
 func NextMarsOpposition(jde float64) float64 {
-	return marsConjunction(jde, 180, 1)
+	return inclusiveNextPhaseEvent(jde, 180, marsConjunction)
 }
 
 func NextMarsEasternQuadrature(jde float64) float64 {
-	return marsConjunction(jde, 90, 1)
+	return inclusiveNextPhaseEvent(jde, 90, marsConjunction)
 }
 
 func LastMarsEasternQuadrature(jde float64) float64 {
-	return marsConjunction(jde, 90, 0)
+	return inclusiveLastPhaseEvent(jde, 90, marsConjunction)
 }
 
 func NextMarsWesternQuadrature(jde float64) float64 {
-	return marsConjunction(jde, 270, 1)
+	return inclusiveNextPhaseEvent(jde, 270, marsConjunction)
 }
 
 func LastMarsWesternQuadrature(jde float64) float64 {
-	return marsConjunction(jde, 270, 0)
+	return inclusiveLastPhaseEvent(jde, 270, marsConjunction)
 }
 
-func marsRetrograde(jde float64, searchBeforeOpposition bool) float64 {
-	//0=last 1=next
-	jde = marsConjunctionFull(jde, 180, 1)
+func marsRetrogradeAroundOpposition(oppositionJD float64, searchBeforeOpposition bool) float64 {
+	jde := oppositionJD
 	if searchBeforeOpposition {
 		jde -= 60
 	} else {
@@ -179,40 +178,66 @@ func marsRetrograde(jde float64, searchBeforeOpposition bool) float64 {
 	return TD2UT(bestJD, false)
 }
 
+func marsOppositionFromBefore(oppositionJD float64) float64 {
+	return marsConjunctionFull(eventUTLastQueryTT(oppositionJD), 180, 1)
+}
+
+func marsOppositionFromAfter(oppositionJD float64) float64 {
+	return marsConjunctionFull(eventUTNextQueryTT(oppositionJD), 180, 0)
+}
+
 func NextMarsRetrogradeToPrograde(jde float64) float64 {
-	date := marsRetrograde(jde, false)
-	if date < jde {
-		oppositionJD := marsConjunctionFull(jde, 180, 1)
-		return marsRetrograde(oppositionJD+10, false)
+	lastOppositionJD := marsConjunctionFull(jde, 180, 0)
+	date := marsRetrogradeAroundOpposition(lastOppositionJD, false)
+	if sameEventUTQueryTT(date, jde) {
+		sameOppositionJD := marsOppositionFromBefore(lastOppositionJD)
+		return closestEventUTToQueryTT(jde, date, marsRetrogradeAroundOpposition(sameOppositionJD, false))
+	}
+	if !eventUTQueryAfterOrEqual(date, jde) {
+		nextOppositionJD := marsConjunctionFull(jde, 180, 1)
+		return marsRetrogradeAroundOpposition(nextOppositionJD, false)
 	}
 	return date
 }
 
 func LastMarsRetrogradeToPrograde(jde float64) float64 {
-	jde = marsConjunctionFull(jde, 180, 0) - 10
-	date := marsRetrograde(jde, false)
-	if date > jde {
-		oppositionJD := marsConjunctionFull(jde, 180, 0)
-		return marsRetrograde(oppositionJD-10, false)
+	lastOppositionJD := marsConjunctionFull(jde, 180, 0)
+	date := marsRetrogradeAroundOpposition(lastOppositionJD, false)
+	if sameEventUTQueryTT(date, jde) {
+		sameOppositionJD := marsOppositionFromBefore(lastOppositionJD)
+		return closestEventUTToQueryTT(jde, date, marsRetrogradeAroundOpposition(sameOppositionJD, false))
+	}
+	if !eventUTQueryBeforeOrEqual(date, jde) {
+		previousOppositionJD := marsConjunctionFull(eventUTLastQueryTT(lastOppositionJD), 180, 0)
+		return marsRetrogradeAroundOpposition(previousOppositionJD, false)
 	}
 	return date
 }
 
 func NextMarsProgradeToRetrograde(jde float64) float64 {
-	date := marsRetrograde(jde, true)
-	if date < jde {
-		oppositionJD := marsConjunctionFull(jde, 180, 1)
-		return marsRetrograde(oppositionJD+10, true)
+	nextOppositionJD := marsConjunctionFull(jde, 180, 1)
+	date := marsRetrogradeAroundOpposition(nextOppositionJD, true)
+	if sameEventUTQueryTT(date, jde) {
+		sameOppositionJD := marsOppositionFromAfter(nextOppositionJD)
+		return closestEventUTToQueryTT(jde, date, marsRetrogradeAroundOpposition(sameOppositionJD, true))
+	}
+	if !eventUTQueryAfterOrEqual(date, jde) {
+		followingOppositionJD := marsConjunctionFull(eventUTNextQueryTT(nextOppositionJD), 180, 1)
+		return marsRetrogradeAroundOpposition(followingOppositionJD, true)
 	}
 	return date
 }
 
 func LastMarsProgradeToRetrograde(jde float64) float64 {
-	jde = marsConjunctionFull(jde, 180, 0) - 10
-	date := marsRetrograde(jde, true)
-	if date > jde {
-		oppositionJD := marsConjunctionFull(jde, 180, 0)
-		return marsRetrograde(oppositionJD-10, true)
+	nextOppositionJD := marsConjunctionFull(jde, 180, 1)
+	date := marsRetrogradeAroundOpposition(nextOppositionJD, true)
+	if sameEventUTQueryTT(date, jde) {
+		sameOppositionJD := marsOppositionFromAfter(nextOppositionJD)
+		return closestEventUTToQueryTT(jde, date, marsRetrogradeAroundOpposition(sameOppositionJD, true))
+	}
+	if !eventUTQueryBeforeOrEqual(date, jde) {
+		lastOppositionJD := marsConjunctionFull(jde, 180, 0)
+		return marsRetrogradeAroundOpposition(lastOppositionJD, true)
 	}
 	return date
 }

basic/mercury.go

@@ -80,52 +80,22 @@ func MercuryApparentRaDec(jd float64) (float64, float64) {
 }
 
 func EarthMercuryAway(jd float64) float64 {
-	x, y, z := AMercuryXYZ(jd)
-	to := math.Sqrt(x*x + y*y + z*z)
-	return to
+	return planetEarthAwayExplicitN(1, jd, -1)
 }
 
 func MercuryApparentLo(jd float64) float64 {
-	x, y, z := AMercuryXYZ(jd)
-	to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
-	x, y, z = AMercuryXYZ(jd - to)
-	lo := math.Atan2(y, x)
-	bo := math.Atan2(z, math.Sqrt(x*x+y*y))
-	lo = lo * 180 / math.Pi
-	bo = bo * 180 / math.Pi
-	lo = Limit360(lo)
-	//lo-=GXCLo(lo,bo,jd)/3600;
-	//bo+=GXCBo(lo,bo,jd);
-	lo += Nutation2000Bi(jd)
-	return lo
+	geo, _ := planetApparentGeocentricPositionN(1, jd, -1)
+	return geo.lo
 }
 
 func MercuryApparentBo(jd float64) float64 {
-	x, y, z := AMercuryXYZ(jd)
-	to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
-	x, y, z = AMercuryXYZ(jd - to)
-	//lo := math.Atan2(y, x)
-	bo := math.Atan2(z, math.Sqrt(x*x+y*y))
-	//lo = lo * 180 / math.Pi
-	bo = bo * 180 / math.Pi
-	//lo+=GXCLo(lo,bo,jd);
-	//bo+=GXCBo(lo,bo,jd)/3600;
-	//lo+=Nutation2000Bi(jd);
-	return bo
+	geo, _ := planetApparentGeocentricPositionN(1, jd, -1)
+	return geo.bo
 }
 
 func MercuryApparentLoBo(jd float64) (float64, float64) {
-	x, y, z := AMercuryXYZ(jd)
-	to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
-	x, y, z = AMercuryXYZ(jd - to)
-	lo := math.Atan2(y, x)
-	bo := math.Atan2(z, math.Sqrt(x*x+y*y))
-	lo = lo * 180 / math.Pi
-	bo = bo * 180 / math.Pi
-	lo = Limit360(lo) + Nutation2000Bi(jd)
-	//lo-=GXCLo(lo,bo,jd)/3600;
-	//bo+=GXCBo(lo,bo,jd);
-	return lo, bo
+	geo, _ := planetApparentGeocentricPositionN(1, jd, -1)
+	return geo.lo, geo.bo
 }
 
 func MercuryMag(jd float64) float64 {

basic/mercury_events.go

@@ -11,6 +11,7 @@ const (
 	MERCURY_S_PERIOD                    = 1 / ((1 / 87.9691) - (1 / 365.256363004))
 	mercuryConjunctionDerivativeStepDay = 2e-5 * 36525.0
 	mercuryLightTimeDaysPerAU           = 0.0057755183
+	mercuryEventSearchN                 = 16
 )
 
 type mercuryConjunctionLBR struct {
@@ -206,41 +207,49 @@ func mercuryConjunction(jde float64, next uint8) float64 {
 }
 
 func LastMercuryConjunction(jde float64) float64 {
-	return mercuryConjunction(jde, 0)
+	return inclusiveLastSimpleEvent(jde, LastMercuryConjunctionStrict, NextMercuryConjunctionStrict)
 }
 
 func NextMercuryConjunction(jde float64) float64 {
+	return inclusiveNextSimpleEvent(jde, LastMercuryConjunctionStrict, NextMercuryConjunctionStrict)
+}
+
+func LastMercuryConjunctionStrict(jde float64) float64 {
+	return mercuryConjunction(jde, 0)
+}
+
+func NextMercuryConjunctionStrict(jde float64) float64 {
 	return mercuryConjunction(jde, 1)
 }
 
 func NextMercuryInferiorConjunction(jde float64) float64 {
-	date := NextMercuryConjunction(jde)
+	date := NextMercuryConjunctionStrict(jde)
 	if EarthMercuryAway(date) > EarthAway(date) {
-		return NextMercuryConjunction(date + 2)
+		return NextMercuryConjunctionStrict(date + 2)
 	}
 	return date
 }
 
 func NextMercurySuperiorConjunction(jde float64) float64 {
-	date := NextMercuryConjunction(jde)
+	date := NextMercuryConjunctionStrict(jde)
 	if EarthMercuryAway(date) < EarthAway(date) {
-		return NextMercuryConjunction(date + 2)
+		return NextMercuryConjunctionStrict(date + 2)
 	}
 	return date
 }
 
 func LastMercuryInferiorConjunction(jde float64) float64 {
-	date := LastMercuryConjunction(jde)
+	date := LastMercuryConjunctionStrict(jde)
 	if EarthMercuryAway(date) > EarthAway(date) {
-		return LastMercuryConjunction(date - 2)
+		return LastMercuryConjunctionStrict(date - 2)
 	}
 	return date
 }
 
 func LastMercurySuperiorConjunction(jde float64) float64 {
-	date := LastMercuryConjunction(jde)
+	date := LastMercuryConjunctionStrict(jde)
 	if EarthMercuryAway(date) < EarthAway(date) {
-		return LastMercuryConjunction(date - 2)
+		return LastMercuryConjunctionStrict(date - 2)
 	}
 	return date
 }
@@ -257,16 +266,6 @@ func mercuryRetrograde(jde float64) float64 {
 		}
 		return sub
 	}
-	raRate := func(jde float64, delta float64) float64 {
-		sub := MercuryApparentRa(jde+delta) - MercuryApparentRa(jde-delta)
-		if sub > 180 {
-			sub -= 360
-		}
-		if sub < -180 {
-			sub += 360
-		}
-		return sub / (2 * delta)
-	}
 	lastConjunction := mercuryConjunctionLegacy(jde, 0)
 	nextConjunction := mercuryConjunctionLegacy(jde, 1)
 	currentRADelta := solarRADelta(jde)
@@ -276,7 +275,7 @@ func mercuryRetrograde(jde float64) float64 {
 		jde = lastConjunction + ((nextConjunction - lastConjunction) / 5.5)
 	}
 	for {
-		currentRate := raRate(jde, 1.0/86400.0)
+		currentRate := mercuryRADerivative(jde, 1.0/86400.0)
 		if math.Abs(currentRate) > 0.55 {
 			jde += 2
 			continue
@@ -286,82 +285,212 @@ func mercuryRetrograde(jde float64) float64 {
 	estimateJD := jde
 	for {
 		prevJD := estimateJD
-		rateValue := raRate(prevJD, 2.0/86400.0)
-		rateSlope := (raRate(prevJD+15.0/86400.0, 2.0/86400.0) - raRate(prevJD-15.0/86400.0, 2.0/86400.0)) / (30.0 / 86400.0)
+		rateValue := mercuryRADerivative(prevJD, 2.0/86400.0)
+		rateSlope := (mercuryRADerivative(prevJD+15.0/86400.0, 2.0/86400.0) - mercuryRADerivative(prevJD-15.0/86400.0, 2.0/86400.0)) / (30.0 / 86400.0)
 		estimateJD = prevJD - rateValue/rateSlope
 		if math.Abs(estimateJD-prevJD) <= 30.0/86400.0 {
 			break
 		}
 	}
 	bestJD := eventZeroRefine(estimateJD, 15.0/86400.0, 0.5/86400.0, func(jd float64) float64 {
-		return raRate(jd, 0.5/86400.0)
+		return mercuryRADerivative(jd, 0.5/86400.0)
 	})
 	//fmt.Println((bestJD - lastConjunction) / (nextConjunction - lastConjunction))
 	return TD2UT(bestJD, false)
 }
 
+func mercuryRADerivative(jde, delta float64) float64 {
+	sub := MercuryApparentRa(jde+delta) - MercuryApparentRa(jde-delta)
+	if sub > 180 {
+		sub -= 360
+	}
+	if sub < -180 {
+		sub += 360
+	}
+	return sub / (2 * delta)
+}
+
+func mercuryStationIsProgradeToRetrograde(eventUT float64) bool {
+	for _, offset := range []float64{0.25, 0.5, 1.0} {
+		before := mercuryRADerivative(eventUT-offset, 0.5/86400.0)
+		after := mercuryRADerivative(eventUT+offset, 0.5/86400.0)
+		if before > 0 && after < 0 {
+			return true
+		}
+		if before < 0 && after > 0 {
+			return false
+		}
+	}
+	before := mercuryRADerivative(eventUT-0.25, 0.5/86400.0)
+	after := mercuryRADerivative(eventUT+0.25, 0.5/86400.0)
+	return before > after
+}
+
+func nextMercuryTypedStation(jde float64, progradeToRetrograde bool) float64 {
+	date := NextMercuryRetrogradeStrict(jde)
+	for mercuryStationIsProgradeToRetrograde(date) != progradeToRetrograde {
+		date = NextMercuryRetrogradeStrict(eventUTNextQueryTT(date))
+	}
+	return date
+}
+
+func lastMercuryTypedStation(jde float64, progradeToRetrograde bool) float64 {
+	date := LastMercuryRetrogradeStrict(jde)
+	for mercuryStationIsProgradeToRetrograde(date) != progradeToRetrograde {
+		date = LastMercuryRetrogradeStrict(eventUTLastQueryTT(date))
+	}
+	return date
+}
+
 func NextMercuryRetrograde(jde float64) float64 {
 	date := mercuryRetrograde(jde)
-	if date < jde {
-		nextConjunction := mercuryConjunctionLegacy(jde, 1)
+	if !eventUTQueryAfterOrEqual(date, jde) {
+		nextConjunction := NextMercuryConjunctionStrict(jde)
 		return mercuryRetrograde(nextConjunction + 2)
 	}
 	return date
 }
 
 func LastMercuryRetrograde(jde float64) float64 {
-	lastConjunction := mercuryConjunctionLegacy(jde, 0)
+	lastConjunction := LastMercuryConjunctionStrict(jde)
 	date := mercuryRetrograde(lastConjunction + 2)
-	if date > jde {
-		previousConjunction := mercuryConjunctionLegacy(lastConjunction-2, 0)
+	if !eventUTQueryBeforeOrEqual(date, jde) {
+		previousConjunction := LastMercuryConjunctionStrict(eventUTLastQueryTT(lastConjunction))
 		return mercuryRetrograde(previousConjunction + 2)
 	}
 	return date
 }
 
+func LastMercuryRetrogradeStrict(jde float64) float64 {
+	return LastMercuryRetrograde(jde)
+}
+
+func NextMercuryRetrogradeStrict(jde float64) float64 {
+	return NextMercuryRetrograde(jde)
+}
+
 func NextMercuryProgradeToRetrograde(jde float64) float64 {
-	date := NextMercuryRetrograde(jde)
-	sub := Limit360(MercuryApparentRa(date) - SunApparentRa(date))
-	if sub > 180 {
-		return NextMercuryRetrograde(date + MERCURY_S_PERIOD/2)
-	}
-	return date
+	return nextMercuryTypedStation(jde, true)
 }
 
 func NextMercuryRetrogradeToPrograde(jde float64) float64 {
-	date := NextMercuryRetrograde(jde)
-	sub := Limit360(MercuryApparentRa(date) - SunApparentRa(date))
-	if sub < 180 {
-		return NextMercuryRetrograde(date + 12)
-	}
-	return date
+	return nextMercuryTypedStation(jde, false)
 }
 
 func LastMercuryProgradeToRetrograde(jde float64) float64 {
-	date := LastMercuryRetrograde(jde)
-	sub := Limit360(MercuryApparentRa(date) - SunApparentRa(date))
-	if sub > 180 {
-		return LastMercuryRetrograde(date - 12)
-	}
-	return date
+	return lastMercuryTypedStation(jde, true)
 }
 
 func LastMercuryRetrogradeToPrograde(jde float64) float64 {
-	date := LastMercuryRetrograde(jde)
-	sub := Limit360(MercuryApparentRa(date) - SunApparentRa(date))
-	if sub < 180 {
-		return LastMercuryRetrograde(date - MERCURY_S_PERIOD/2)
-	}
-	return date
+	return lastMercuryTypedStation(jde, false)
 }
 
 func MercurySunElongation(jde float64) float64 {
 	lo1, bo1 := MercuryApparentLoBo(jde)
-	lo2 := SunApparentLo(jde)
+	lo2 := HSunApparentLo(jde)
 	bo2 := HSunTrueBo(jde)
 	return StarAngularSeparation(lo1, bo1, lo2, bo2)
 }
 
+func mercurySunElongationN(jde float64, n int) float64 {
+	lo1, bo1 := MercuryApparentLoBoN(jde, n)
+	lo2 := HSunApparentLoN(jde, n)
+	bo2 := HSunTrueBoN(jde, n)
+	return StarAngularSeparation(lo1, bo1, lo2, bo2)
+}
+
+func mercuryTrueElongationN(jde float64, n int) float64 {
+	earth := mercuryHelioN(-1, jde, n)
+	planetPos := mercuryHelioN(1, jde, n)
+	geo := mercuryGeocentric(planetPos, earth)
+	return StarAngularSeparation(geo.lo, geo.bo, HSunTrueLoN(jde, n), HSunTrueBoN(jde, n))
+}
+
+func mercuryGreatestElongationInWindow(start, end float64) float64 {
+	best := maximizeInWindow(start, end, 2.0, func(jd float64) float64 {
+		return mercuryTrueElongationN(jd, mercuryEventSearchN)
+	}, func(jd float64) float64 {
+		return mercuryTrueElongationN(jd, -1)
+	})
+	return TD2UT(best, false)
+}
+
+func mercuryEastElongationWindowEndingAt(inferior float64) (float64, float64) {
+	lastSuperior := LastMercurySuperiorConjunction(eventUTLastQueryTT(inferior))
+	return lastSuperior + innerEventEpsilon, inferior - innerEventEpsilon
+}
+
+func mercuryWestElongationWindowEndingAt(superior float64) (float64, float64) {
+	lastInferior := LastMercuryInferiorConjunction(eventUTLastQueryTT(superior))
+	return lastInferior + innerEventEpsilon, superior - innerEventEpsilon
+}
+
+func mercuryEastElongationWindowContaining(jde float64) (float64, float64) {
+	nextInferior := NextMercuryInferiorConjunction(jde)
+	start, end := mercuryEastElongationWindowEndingAt(nextInferior)
+	if eventUTQueryBeforeOrEqual(start, jde) {
+		return start, end
+	}
+	currentInferior := LastMercuryInferiorConjunction(jde)
+	return mercuryEastElongationWindowEndingAt(currentInferior)
+}
+
+func mercuryWestElongationWindowContaining(jde float64) (float64, float64) {
+	nextSuperior := NextMercurySuperiorConjunction(jde)
+	start, end := mercuryWestElongationWindowEndingAt(nextSuperior)
+	if eventUTQueryBeforeOrEqual(start, jde) {
+		return start, end
+	}
+	currentSuperior := LastMercurySuperiorConjunction(jde)
+	return mercuryWestElongationWindowEndingAt(currentSuperior)
+}
+
+func nextMercuryGreatestElongationTyped(jde float64, east bool) float64 {
+	if east {
+		start, windowEnd := mercuryEastElongationWindowContaining(jde)
+		for {
+			date := mercuryGreatestElongationInWindow(start, windowEnd)
+			if eventUTQueryAfterOrEqual(date, jde) {
+				return date
+			}
+			nextInferior := NextMercuryInferiorConjunction(eventUTNextQueryTT(windowEnd))
+			start, windowEnd = mercuryEastElongationWindowEndingAt(nextInferior)
+		}
+	}
+	start, windowEnd := mercuryWestElongationWindowContaining(jde)
+	for {
+		date := mercuryGreatestElongationInWindow(start, windowEnd)
+		if eventUTQueryAfterOrEqual(date, jde) {
+			return date
+		}
+		nextSuperior := NextMercurySuperiorConjunction(eventUTNextQueryTT(windowEnd))
+		start, windowEnd = mercuryWestElongationWindowEndingAt(nextSuperior)
+	}
+}
+
+func lastMercuryGreatestElongationTyped(jde float64, east bool) float64 {
+	if east {
+		start, windowEnd := mercuryEastElongationWindowContaining(jde)
+		for {
+			date := mercuryGreatestElongationInWindow(start, windowEnd)
+			if eventUTQueryBeforeOrEqual(date, jde) {
+				return date
+			}
+			prevInferior := LastMercuryInferiorConjunction(eventUTLastQueryTT(start))
+			start, windowEnd = mercuryEastElongationWindowEndingAt(prevInferior)
+		}
+	}
+	start, windowEnd := mercuryWestElongationWindowContaining(jde)
+	for {
+		date := mercuryGreatestElongationInWindow(start, windowEnd)
+		if eventUTQueryBeforeOrEqual(date, jde) {
+			return date
+		}
+		prevSuperior := LastMercurySuperiorConjunction(eventUTLastQueryTT(start))
+		start, windowEnd = mercuryWestElongationWindowEndingAt(prevSuperior)
+	}
+}
+
 func mercuryGreatestElongation(jde float64) float64 {
 	solarRADelta := func(jde float64) float64 {
 		sub := Limit360(MercuryApparentRa(jde) - SunApparentRa(jde))
@@ -383,8 +512,8 @@ func mercuryGreatestElongation(jde float64) float64 {
 		}
 		return sub / (2 * delta)
 	}
-	lastConjunction := mercuryConjunctionLegacy(jde, 0)
-	nextConjunction := mercuryConjunctionLegacy(jde, 1)
+	lastConjunction := LastMercuryConjunctionStrict(jde)
+	nextConjunction := NextMercuryConjunctionStrict(jde)
 	currentRADelta := solarRADelta(jde)
 	if currentRADelta > 0 {
 		jde = lastConjunction + ((nextConjunction - lastConjunction) / 5.0 * 2.0)
@@ -417,56 +546,105 @@ func mercuryGreatestElongation(jde float64) float64 {
 }
 
 func NextMercuryGreatestElongation(jde float64) float64 {
-	date := mercuryGreatestElongation(jde)
-	if date < jde {
-		nextConjunction := mercuryConjunctionLegacy(jde, 1)
-		return mercuryGreatestElongation(nextConjunction + 2)
+	east := NextMercuryGreatestElongationEast(jde)
+	west := NextMercuryGreatestElongationWest(jde)
+	if sameEventJD(east, west) {
+		return east
 	}
-	return date
+	if east < west {
+		return east
+	}
+	return west
 }
 
 func LastMercuryGreatestElongation(jde float64) float64 {
-	lastConjunction := mercuryConjunctionLegacy(jde, 0)
-	date := mercuryGreatestElongation(lastConjunction + 2)
-	if date > jde {
-		previousConjunction := mercuryConjunctionLegacy(lastConjunction-2, 0)
-		return mercuryGreatestElongation(previousConjunction + 2)
+	east := LastMercuryGreatestElongationEast(jde)
+	west := LastMercuryGreatestElongationWest(jde)
+	if sameEventJD(east, west) {
+		return east
 	}
-	return date
+	if east > west {
+		return east
+	}
+	return west
+}
+
+func LastMercuryInferiorConjunctionInclusive(jde float64) float64 {
+	return inclusiveLastSimpleEvent(jde, LastMercuryInferiorConjunction, NextMercuryInferiorConjunction)
+}
+
+func NextMercuryInferiorConjunctionInclusive(jde float64) float64 {
+	return inclusiveNextSimpleEvent(jde, LastMercuryInferiorConjunction, NextMercuryInferiorConjunction)
+}
+
+func LastMercurySuperiorConjunctionInclusive(jde float64) float64 {
+	return inclusiveLastSimpleEvent(jde, LastMercurySuperiorConjunction, NextMercurySuperiorConjunction)
+}
+
+func NextMercurySuperiorConjunctionInclusive(jde float64) float64 {
+	return inclusiveNextSimpleEvent(jde, LastMercurySuperiorConjunction, NextMercurySuperiorConjunction)
+}
+
+func LastMercuryRetrogradeInclusive(jde float64) float64 {
+	return inclusiveLastSimpleEvent(jde, LastMercuryRetrograde, NextMercuryRetrograde)
+}
+
+func NextMercuryRetrogradeInclusive(jde float64) float64 {
+	return inclusiveNextSimpleEvent(jde, LastMercuryRetrograde, NextMercuryRetrograde)
+}
+
+func LastMercuryProgradeToRetrogradeInclusive(jde float64) float64 {
+	return inclusiveLastSimpleEvent(jde, LastMercuryProgradeToRetrograde, NextMercuryProgradeToRetrograde)
+}
+
+func NextMercuryProgradeToRetrogradeInclusive(jde float64) float64 {
+	return inclusiveNextSimpleEvent(jde, LastMercuryProgradeToRetrograde, NextMercuryProgradeToRetrograde)
+}
+
+func LastMercuryRetrogradeToProgradeInclusive(jde float64) float64 {
+	return inclusiveLastSimpleEvent(jde, LastMercuryRetrogradeToPrograde, NextMercuryRetrogradeToPrograde)
+}
+
+func NextMercuryRetrogradeToProgradeInclusive(jde float64) float64 {
+	return inclusiveNextSimpleEvent(jde, LastMercuryRetrogradeToPrograde, NextMercuryRetrogradeToPrograde)
+}
+
+func LastMercuryGreatestElongationInclusive(jde float64) float64 {
+	return inclusiveLastSimpleEvent(jde, LastMercuryGreatestElongation, NextMercuryGreatestElongation)
+}
+
+func NextMercuryGreatestElongationInclusive(jde float64) float64 {
+	return inclusiveNextSimpleEvent(jde, LastMercuryGreatestElongation, NextMercuryGreatestElongation)
+}
+
+func LastMercuryGreatestElongationEastInclusive(jde float64) float64 {
+	return inclusiveLastSimpleEvent(jde, LastMercuryGreatestElongationEast, NextMercuryGreatestElongationEast)
+}
+
+func NextMercuryGreatestElongationEastInclusive(jde float64) float64 {
+	return inclusiveNextSimpleEvent(jde, LastMercuryGreatestElongationEast, NextMercuryGreatestElongationEast)
+}
+
+func LastMercuryGreatestElongationWestInclusive(jde float64) float64 {
+	return inclusiveLastSimpleEvent(jde, LastMercuryGreatestElongationWest, NextMercuryGreatestElongationWest)
+}
+
+func NextMercuryGreatestElongationWestInclusive(jde float64) float64 {
+	return inclusiveNextSimpleEvent(jde, LastMercuryGreatestElongationWest, NextMercuryGreatestElongationWest)
 }
 
 func NextMercuryGreatestElongationEast(jde float64) float64 {
-	date := NextMercuryGreatestElongation(jde)
-	sub := Limit360(MercuryApparentRa(date) - SunApparentRa(date))
-	if sub > 180 {
-		return NextMercuryGreatestElongation(date + 1)
-	}
-	return date
+	return nextMercuryGreatestElongationTyped(jde, true)
 }
 
 func NextMercuryGreatestElongationWest(jde float64) float64 {
-	date := NextMercuryGreatestElongation(jde)
-	sub := Limit360(MercuryApparentRa(date) - SunApparentRa(date))
-	if sub < 180 {
-		return NextMercuryGreatestElongation(date + 1)
-	}
-	return date
+	return nextMercuryGreatestElongationTyped(jde, false)
 }
 
 func LastMercuryGreatestElongationEast(jde float64) float64 {
-	date := LastMercuryGreatestElongation(jde)
-	sub := Limit360(MercuryApparentRa(date) - SunApparentRa(date))
-	if sub > 180 {
-		return LastMercuryGreatestElongation(date - 1)
-	}
-	return date
+	return lastMercuryGreatestElongationTyped(jde, true)
 }
 
 func LastMercuryGreatestElongationWest(jde float64) float64 {
-	date := LastMercuryGreatestElongation(jde)
-	sub := Limit360(MercuryApparentRa(date) - SunApparentRa(date))
-	if sub < 180 {
-		return LastMercuryGreatestElongation(date - 1)
-	}
-	return date
+	return lastMercuryGreatestElongationTyped(jde, false)
 }

basic/mercury_station_regression_test.go

@@ -0,0 +1,40 @@
+package basic
+
+import (
+	"math"
+	"testing"
+	"time"
+)
+
+func mercuryTTJDJST(year int, month time.Month, day, hour, minute, second int) float64 {
+	loc := time.FixedZone("JST", 9*3600)
+	return TD2UT(Date2JDE(time.Date(year, month, day, hour, minute, second, 0, loc).UTC()), true)
+}
+
+func TestMercuryTypedStationRegression1929(t *testing.T) {
+	loc := time.FixedZone("JST", 9*3600)
+	const tolerance = 30.0 / 86400.0
+
+	query := mercuryTTJDJST(1929, time.September, 20, 0, 0, 0)
+	wantP2R := mercuryTTJDJST(1929, time.September, 26, 1, 58, 0)
+	wantR2P := mercuryTTJDJST(1929, time.October, 16, 23, 32, 33)
+
+	nextP2R := NextMercuryProgradeToRetrograde(query)
+	nextR2P := NextMercuryRetrogradeToPrograde(query)
+	if math.Abs(nextP2R-wantP2R) > tolerance {
+		t.Fatalf("next P2R mismatch: got %s want %s", JDE2DateByZone(nextP2R, loc, false), JDE2DateByZone(wantP2R, loc, false))
+	}
+	if math.Abs(nextR2P-wantR2P) > tolerance {
+		t.Fatalf("next R2P mismatch: got %s want %s", JDE2DateByZone(nextR2P, loc, false), JDE2DateByZone(wantR2P, loc, false))
+	}
+
+	query = mercuryTTJDJST(1929, time.October, 20, 0, 0, 0)
+	lastP2R := LastMercuryProgradeToRetrograde(query)
+	lastR2P := LastMercuryRetrogradeToPrograde(query)
+	if math.Abs(lastP2R-wantP2R) > tolerance {
+		t.Fatalf("last P2R mismatch: got %s want %s", JDE2DateByZone(lastP2R, loc, false), JDE2DateByZone(wantP2R, loc, false))
+	}
+	if math.Abs(lastR2P-wantR2P) > tolerance {
+		t.Fatalf("last R2P mismatch: got %s want %s", JDE2DateByZone(lastR2P, loc, false), JDE2DateByZone(wantR2P, loc, false))
+	}
+}

basic/neptune.go

@@ -87,53 +87,22 @@ func NeptuneApparentRaDec(jd float64) (float64, float64) {
 }
 
 func EarthNeptuneAway(jd float64) float64 {
-	x, y, z := ANeptuneXYZ(jd)
-	to := math.Sqrt(x*x + y*y + z*z)
-	return to
+	return planetEarthAwayExplicitN(7, jd, -1)
 }
 
 func NeptuneApparentLo(jd float64) float64 {
-	x, y, z := ANeptuneXYZ(jd)
-	to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
-	x, y, z = ANeptuneXYZ(jd - to)
-	lo := math.Atan2(y, x)
-	bo := math.Atan2(z, math.Sqrt(x*x+y*y))
-	lo = lo * 180 / math.Pi
-	bo = bo * 180 / math.Pi
-	lo = Limit360(lo)
-	//lo-=GXCLo(lo,bo,jd)/3600;
-	//bo+=GXCBo(lo,bo,jd);
-	lo += Nutation2000Bi(jd)
-	return lo
+	geo, _ := planetApparentGeocentricPositionN(7, jd, -1)
+	return geo.lo
 }
 
 func NeptuneApparentBo(jd float64) float64 {
-	x, y, z := ANeptuneXYZ(jd)
-	to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
-	x, y, z = ANeptuneXYZ(jd - to)
-	//lo := math.Atan2(y, x)
-	bo := math.Atan2(z, math.Sqrt(x*x+y*y))
-	//lo = lo * 180 / math.Pi
-	bo = bo * 180 / math.Pi
-	//lo+=GXCLo(lo,bo,jd);
-	//bo+=GXCBo(lo,bo,jd)/3600;
-	//lo+=Nutation2000Bi(jd);
-	return bo
+	geo, _ := planetApparentGeocentricPositionN(7, jd, -1)
+	return geo.bo
 }
 
 func NeptuneApparentLoBo(jd float64) (float64, float64) {
-	x, y, z := ANeptuneXYZ(jd)
-	to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
-	x, y, z = ANeptuneXYZ(jd - to)
-	lo := math.Atan2(y, x)
-	bo := math.Atan2(z, math.Sqrt(x*x+y*y))
-	lo = lo * 180 / math.Pi
-	bo = bo * 180 / math.Pi
-	lo = Limit360(lo)
-	//lo-=GXCLo(lo,bo,jd)/3600;
-	//bo+=GXCBo(lo,bo,jd);
-	lo += Nutation2000Bi(jd)
-	return lo, bo
+	geo, _ := planetApparentGeocentricPositionN(7, jd, -1)
+	return geo.lo, geo.bo
 }
 
 func NeptuneMag(jd float64) float64 {

basic/neptune_events.go

@@ -9,7 +9,7 @@ import (
 // Pos
 
 const (
-	NEPTUNE_S_PERIOD            = 1 / ((1 / 365.256363004) - (1 / 4332.59))
+	NEPTUNE_S_PERIOD            = 1 / ((1 / 365.256363004) - (1 / 60190.03))
 	neptuneEventSearchN         = 16
 	neptunePhaseCoarseTolerance = 30.0 / 86400.0
 )
@@ -40,6 +40,28 @@ func neptuneSunLongitudeDeltaN(jde, degree float64, filter bool, n int) float64
 	return sub
 }
 
+func neptuneRADerivative(jde, delta float64) float64 {
+	sub := NeptuneApparentRa(jde+delta) - NeptuneApparentRa(jde-delta)
+	if sub > 180 {
+		sub -= 360
+	}
+	if sub < -180 {
+		sub += 360
+	}
+	return sub / (2 * delta)
+}
+
+func neptuneRADerivativeN(jde, delta float64, n int) float64 {
+	sub := NeptuneApparentRaN(jde+delta, n) - NeptuneApparentRaN(jde-delta, n)
+	if sub > 180 {
+		sub -= 360
+	}
+	if sub < -180 {
+		sub += 360
+	}
+	return sub / (2 * delta)
+}
+
 func neptuneConjunctionFull(jde, degree float64, next uint8) float64 {
 	//0=last 1=next
 	daysPerDegree := NEPTUNE_S_PERIOD / 360
@@ -94,113 +116,92 @@ func neptuneConjunction(jde, degree float64, next uint8) float64 {
 }
 
 func LastNeptuneConjunction(jde float64) float64 {
-	return neptuneConjunction(jde, 0, 0)
+	return inclusiveLastPhaseEvent(jde, 0, neptuneConjunction)
 }
 
 func NextNeptuneConjunction(jde float64) float64 {
-	return neptuneConjunction(jde, 0, 1)
+	return inclusiveNextPhaseEvent(jde, 0, neptuneConjunction)
 }
 
 func LastNeptuneOpposition(jde float64) float64 {
-	return neptuneConjunction(jde, 180, 0)
+	return inclusiveLastPhaseEvent(jde, 180, neptuneConjunction)
 }
 
 func NextNeptuneOpposition(jde float64) float64 {
-	return neptuneConjunction(jde, 180, 1)
+	return inclusiveNextPhaseEvent(jde, 180, neptuneConjunction)
 }
 
 func NextNeptuneEasternQuadrature(jde float64) float64 {
-	return neptuneConjunction(jde, 90, 1)
+	return inclusiveNextPhaseEvent(jde, 90, neptuneConjunction)
 }
 
 func LastNeptuneEasternQuadrature(jde float64) float64 {
-	return neptuneConjunction(jde, 90, 0)
+	return inclusiveLastPhaseEvent(jde, 90, neptuneConjunction)
 }
 
 func NextNeptuneWesternQuadrature(jde float64) float64 {
-	return neptuneConjunction(jde, 270, 1)
+	return inclusiveNextPhaseEvent(jde, 270, neptuneConjunction)
 }
 
 func LastNeptuneWesternQuadrature(jde float64) float64 {
-	return neptuneConjunction(jde, 270, 0)
+	return inclusiveLastPhaseEvent(jde, 270, neptuneConjunction)
 }
 
-func neptuneRetrograde(jde float64, searchBeforeOpposition bool) float64 {
-	//0=last 1=next
-	raRate := func(jde float64, delta float64) float64 {
-		sub := NeptuneApparentRa(jde+delta) - NeptuneApparentRa(jde-delta)
-		if sub > 180 {
-			sub -= 360
-		}
-		if sub < -180 {
-			sub += 360
-		}
-		return sub / (2 * delta)
-	}
-	jde = neptuneConjunctionFull(jde, 180, 1)
+func neptuneRetrogradeAroundOpposition(oppositionJD float64, searchBeforeOpposition bool) float64 {
+	oppositionTT := TD2UT(oppositionJD, true)
+	startTT := oppositionTT
+	endTT := oppositionTT
 	if searchBeforeOpposition {
-		jde -= 60
+		easternQuadratureUT := neptuneConjunction(oppositionTT, 90, 0)
+		startTT = TD2UT(easternQuadratureUT, true)
 	} else {
-		jde += 60
+		westernQuadratureUT := neptuneConjunction(oppositionTT, 270, 1)
+		endTT = TD2UT(westernQuadratureUT, true)
 	}
-	for {
-		currentRate := raRate(jde, 1.0/86400.0)
-		if math.Abs(currentRate) > 0.55 {
-			jde += 2
-			continue
-		}
-		break
-	}
-	estimateJD := jde
-	for {
-		prevJD := estimateJD
-		rateValue := raRate(prevJD, 2.0/86400.0)
-		rateSlope := (raRate(prevJD+15.0/86400.0, 2.0/86400.0) - raRate(prevJD-15.0/86400.0, 2.0/86400.0)) / (30.0 / 86400.0)
-		estimateJD = prevJD - rateValue/rateSlope
-		if math.Abs(estimateJD-prevJD) <= 30.0/86400.0 {
-			break
-		}
-	}
-	bestJD := eventZeroRefine(estimateJD, 15.0/86400.0, 0.5/86400.0, func(jd float64) float64 {
-		return raRate(jd, 0.5/86400.0)
+	bestJD := zeroEventInWindow(startTT, endTT, 2.0, 2.0, 30.0/86400.0, func(jd float64) float64 {
+		return neptuneRADerivativeN(jd, 1.0/86400.0, neptuneEventSearchN)
+	}, func(jd float64) float64 {
+		return neptuneRADerivative(jd, 0.5/86400.0)
 	})
 	return TD2UT(bestJD, false)
 }
 
 func NextNeptuneRetrogradeToPrograde(jde float64) float64 {
-	date := neptuneRetrograde(jde, false)
-	if date < jde {
-		oppositionJD := neptuneConjunctionFull(jde, 180, 1)
-		return neptuneRetrograde(oppositionJD+10, false)
+	lastOppositionJD := neptuneConjunctionFull(jde, 180, 0)
+	date := neptuneRetrogradeAroundOpposition(lastOppositionJD, false)
+	if sameEventUTQueryTT(date, jde) || eventUTQueryAfterOrEqual(date, jde) {
+		return date
 	}
-	return date
+	nextOppositionJD := neptuneConjunctionFull(jde, 180, 1)
+	return neptuneRetrogradeAroundOpposition(nextOppositionJD, false)
 }
 
 func LastNeptuneRetrogradeToPrograde(jde float64) float64 {
-	jde = neptuneConjunctionFull(jde, 180, 0) - 10
-	date := neptuneRetrograde(jde, false)
-	if date > jde {
-		oppositionJD := neptuneConjunctionFull(jde, 180, 0)
-		return neptuneRetrograde(oppositionJD-10, false)
+	lastOppositionJD := neptuneConjunctionFull(jde, 180, 0)
+	date := neptuneRetrogradeAroundOpposition(lastOppositionJD, false)
+	if sameEventUTQueryTT(date, jde) || eventUTQueryBeforeOrEqual(date, jde) {
+		return date
 	}
-	return date
+	previousOppositionJD := neptuneConjunctionFull(eventUTLastQueryTT(lastOppositionJD), 180, 0)
+	return neptuneRetrogradeAroundOpposition(previousOppositionJD, false)
 }
 
 func NextNeptuneProgradeToRetrograde(jde float64) float64 {
-	date := neptuneRetrograde(jde, true)
-	if date < jde {
-		oppositionJD := neptuneConjunctionFull(jde, 180, 1)
-		return neptuneRetrograde(oppositionJD+10, true)
+	nextOppositionJD := neptuneConjunctionFull(jde, 180, 1)
+	date := neptuneRetrogradeAroundOpposition(nextOppositionJD, true)
+	if sameEventUTQueryTT(date, jde) || eventUTQueryAfterOrEqual(date, jde) {
+		return date
 	}
-	return date
+	followingOppositionJD := neptuneConjunctionFull(eventUTNextQueryTT(nextOppositionJD), 180, 1)
+	return neptuneRetrogradeAroundOpposition(followingOppositionJD, true)
 }
 
 func LastNeptuneProgradeToRetrograde(jde float64) float64 {
-	jde = neptuneConjunctionFull(jde, 180, 0) - 10
-	date := neptuneRetrograde(jde, true)
-	if date > jde {
-		oppositionJD := neptuneConjunctionFull(jde, 180, 0)
-		return neptuneRetrograde(oppositionJD-10, true)
+	nextOppositionJD := neptuneConjunctionFull(jde, 180, 1)
+	date := neptuneRetrogradeAroundOpposition(nextOppositionJD, true)
+	if sameEventUTQueryTT(date, jde) || eventUTQueryBeforeOrEqual(date, jde) {
+		return date
 	}
-	return date
+	lastOppositionJD := neptuneConjunctionFull(jde, 180, 0)
+	return neptuneRetrogradeAroundOpposition(lastOppositionJD, true)
 }

basic/outer_planet_event_boundary_test.go

@@ -0,0 +1,55 @@
+package basic
+
+import (
+	"testing"
+	"time"
+)
+
+func TestOuterPlanetExactEventBoundaryIncludesCurrent(t *testing.T) {
+	cases := []struct {
+		name   string
+		seed   float64
+		lastFn func(float64) float64
+		nextFn func(float64) float64
+	}{
+		{name: "JupiterConjunction", seed: NextJupiterConjunction(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastJupiterConjunction, nextFn: NextJupiterConjunction},
+		{name: "JupiterOpposition", seed: NextJupiterOpposition(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastJupiterOpposition, nextFn: NextJupiterOpposition},
+		{name: "JupiterEasternQuadrature", seed: NextJupiterEasternQuadrature(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastJupiterEasternQuadrature, nextFn: NextJupiterEasternQuadrature},
+		{name: "JupiterWesternQuadrature", seed: NextJupiterWesternQuadrature(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastJupiterWesternQuadrature, nextFn: NextJupiterWesternQuadrature},
+		{name: "JupiterP2R", seed: NextJupiterProgradeToRetrograde(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastJupiterProgradeToRetrograde, nextFn: NextJupiterProgradeToRetrograde},
+		{name: "JupiterR2P", seed: NextJupiterRetrogradeToPrograde(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastJupiterRetrogradeToPrograde, nextFn: NextJupiterRetrogradeToPrograde},
+		{name: "SaturnOpposition", seed: NextSaturnOpposition(ttjdUTC(2025, 1, 1, 0, 0, 0)), lastFn: LastSaturnOpposition, nextFn: NextSaturnOpposition},
+		{name: "SaturnP2R", seed: NextSaturnProgradeToRetrograde(ttjdUTC(2025, 1, 1, 0, 0, 0)), lastFn: LastSaturnProgradeToRetrograde, nextFn: NextSaturnProgradeToRetrograde},
+		{name: "SaturnR2P", seed: NextSaturnRetrogradeToPrograde(ttjdUTC(2025, 1, 1, 0, 0, 0)), lastFn: LastSaturnRetrogradeToPrograde, nextFn: NextSaturnRetrogradeToPrograde},
+		{name: "UranusOpposition", seed: NextUranusOpposition(ttjdUTC(2025, 1, 1, 0, 0, 0)), lastFn: LastUranusOpposition, nextFn: NextUranusOpposition},
+		{name: "UranusP2R", seed: NextUranusProgradeToRetrograde(ttjdUTC(2025, 1, 1, 0, 0, 0)), lastFn: LastUranusProgradeToRetrograde, nextFn: NextUranusProgradeToRetrograde},
+		{name: "UranusR2P", seed: NextUranusRetrogradeToPrograde(ttjdUTC(2025, 1, 1, 0, 0, 0)), lastFn: LastUranusRetrogradeToPrograde, nextFn: NextUranusRetrogradeToPrograde},
+		{name: "NeptuneOpposition", seed: NextNeptuneOpposition(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastNeptuneOpposition, nextFn: NextNeptuneOpposition},
+		{name: "NeptuneP2R", seed: NextNeptuneProgradeToRetrograde(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastNeptuneProgradeToRetrograde, nextFn: NextNeptuneProgradeToRetrograde},
+		{name: "NeptuneR2P", seed: NextNeptuneRetrogradeToPrograde(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastNeptuneRetrogradeToPrograde, nextFn: NextNeptuneRetrogradeToPrograde},
+		{name: "MarsConjunction", seed: NextMarsConjunction(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastMarsConjunction, nextFn: NextMarsConjunction},
+		{name: "MarsOpposition", seed: NextMarsOpposition(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastMarsOpposition, nextFn: NextMarsOpposition},
+		{name: "MarsEasternQuadrature", seed: NextMarsEasternQuadrature(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastMarsEasternQuadrature, nextFn: NextMarsEasternQuadrature},
+		{name: "MarsWesternQuadrature", seed: NextMarsWesternQuadrature(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastMarsWesternQuadrature, nextFn: NextMarsWesternQuadrature},
+		{name: "MarsP2R", seed: NextMarsProgradeToRetrograde(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastMarsProgradeToRetrograde, nextFn: NextMarsProgradeToRetrograde},
+		{name: "MarsR2P", seed: NextMarsRetrogradeToPrograde(ttjdUTC(2026, 1, 1, 0, 0, 0)), lastFn: LastMarsRetrogradeToPrograde, nextFn: NextMarsRetrogradeToPrograde},
+	}
+
+	for _, tc := range cases {
+		t.Run(tc.name, func(t *testing.T) {
+			queryTT := TD2UT(tc.seed, true)
+			last := tc.lastFn(queryTT)
+			next := tc.nextFn(queryTT)
+			if !sameEventJD(last, tc.seed) {
+				t.Fatalf("last exact boundary mismatch: got %.12f want %.12f", last, tc.seed)
+			}
+			if !sameEventJD(next, tc.seed) {
+				t.Fatalf("next exact boundary mismatch: got %.12f want %.12f", next, tc.seed)
+			}
+		})
+	}
+}
+
+func ttjdUTC(year, month, day, hour, min, sec int) float64 {
+	return TD2UT(Date2JDE(time.Date(year, time.Month(month), day, hour, min, sec, 0, time.UTC)), true)
+}

basic/outer_planet_truth_test.go

@@ -0,0 +1,181 @@
+package basic
+
+import (
+	"encoding/json"
+	"os"
+	"strings"
+	"testing"
+	"time"
+)
+
+type outerTruthBaselineFile struct {
+	Events []outerTruthBaselineEvent `json:"events"`
+}
+
+type outerTruthBaselineEvent struct {
+	Planet          string `json:"planet"`
+	Kind            string `json:"kind"`
+	HintKind        string `json:"hint_kind"`
+	NAOJHintJST     string `json:"naoj_hint_jst"`
+	Precision       string `json:"precision"`
+	CandidateJST    string `json:"candidate_jst"`
+	VerifiedJST     string `json:"verified_jst"`
+	CandidateSource string `json:"candidate_source"`
+}
+
+func loadOuterTruthBaseline(t *testing.T) outerTruthBaselineFile {
+	t.Helper()
+
+	paths := [][]string{
+		{
+			"testdata/jpl_outer_event_baseline.json",
+			"basic/testdata/jpl_outer_event_baseline.json",
+		},
+		{
+			"testdata/jpl_outer_event_baseline_21c_sample.json",
+			"basic/testdata/jpl_outer_event_baseline_21c_sample.json",
+		},
+	}
+
+	var merged outerTruthBaselineFile
+	for index, candidates := range paths {
+		var (
+			data []byte
+			err  error
+		)
+		for _, path := range candidates {
+			data, err = os.ReadFile(path)
+			if err == nil {
+				var baseline outerTruthBaselineFile
+				if err := json.Unmarshal(data, &baseline); err != nil {
+					t.Fatal(err)
+				}
+				merged.Events = append(merged.Events, baseline.Events...)
+				break
+			}
+		}
+		if err != nil && index == 0 {
+			t.Fatal(err)
+		}
+	}
+	if len(merged.Events) == 0 {
+		t.Fatal("empty outer truth baseline file")
+	}
+	return merged
+}
+
+func outerTruthTolerance(event outerTruthBaselineEvent) time.Duration {
+	switch event.Kind {
+	case "CONJ", "OPP", "EQE", "EQW":
+		return 2 * time.Minute
+	default:
+		return 2 * time.Minute
+	}
+}
+
+func outerTruthEventFuncs(t *testing.T, event outerTruthBaselineEvent) (func(float64) float64, func(float64) float64) {
+	t.Helper()
+	switch event.Planet + ":" + event.Kind {
+	case "Jupiter:CONJ":
+		return LastJupiterConjunction, NextJupiterConjunction
+	case "Jupiter:OPP":
+		return LastJupiterOpposition, NextJupiterOpposition
+	case "Jupiter:EQE":
+		return LastJupiterEasternQuadrature, NextJupiterEasternQuadrature
+	case "Jupiter:EQW":
+		return LastJupiterWesternQuadrature, NextJupiterWesternQuadrature
+	case "Jupiter:P2R":
+		return LastJupiterProgradeToRetrograde, NextJupiterProgradeToRetrograde
+	case "Jupiter:R2P":
+		return LastJupiterRetrogradeToPrograde, NextJupiterRetrogradeToPrograde
+	case "Saturn:CONJ":
+		return LastSaturnConjunction, NextSaturnConjunction
+	case "Saturn:OPP":
+		return LastSaturnOpposition, NextSaturnOpposition
+	case "Saturn:EQE":
+		return LastSaturnEasternQuadrature, NextSaturnEasternQuadrature
+	case "Saturn:EQW":
+		return LastSaturnWesternQuadrature, NextSaturnWesternQuadrature
+	case "Saturn:P2R":
+		return LastSaturnProgradeToRetrograde, NextSaturnProgradeToRetrograde
+	case "Saturn:R2P":
+		return LastSaturnRetrogradeToPrograde, NextSaturnRetrogradeToPrograde
+	case "Uranus:CONJ":
+		return LastUranusConjunction, NextUranusConjunction
+	case "Uranus:OPP":
+		return LastUranusOpposition, NextUranusOpposition
+	case "Uranus:EQE":
+		return LastUranusEasternQuadrature, NextUranusEasternQuadrature
+	case "Uranus:EQW":
+		return LastUranusWesternQuadrature, NextUranusWesternQuadrature
+	case "Uranus:P2R":
+		return LastUranusProgradeToRetrograde, NextUranusProgradeToRetrograde
+	case "Uranus:R2P":
+		return LastUranusRetrogradeToPrograde, NextUranusRetrogradeToPrograde
+	case "Neptune:CONJ":
+		return LastNeptuneConjunction, NextNeptuneConjunction
+	case "Neptune:OPP":
+		return LastNeptuneOpposition, NextNeptuneOpposition
+	case "Neptune:EQE":
+		return LastNeptuneEasternQuadrature, NextNeptuneEasternQuadrature
+	case "Neptune:EQW":
+		return LastNeptuneWesternQuadrature, NextNeptuneWesternQuadrature
+	case "Neptune:P2R":
+		return LastNeptuneProgradeToRetrograde, NextNeptuneProgradeToRetrograde
+	case "Neptune:R2P":
+		return LastNeptuneRetrogradeToPrograde, NextNeptuneRetrogradeToPrograde
+	default:
+		t.Fatalf("unsupported outer event %s:%s", event.Planet, event.Kind)
+		return nil, nil
+	}
+}
+
+func assertOuterTruthBaselineEvent(t *testing.T, event outerTruthBaselineEvent, lastFn, nextFn func(float64) float64) {
+	t.Helper()
+	when := parseInnerBaselineTime(t, event.VerifiedJST)
+	before := when.Add(-7 * 24 * time.Hour)
+	after := when.Add(7 * 24 * time.Hour)
+	next := JDE2DateByZone(nextFn(toUTJD(before)), when.Location(), false)
+	last := JDE2DateByZone(lastFn(toUTJD(after)), when.Location(), false)
+	tolerance := outerTruthTolerance(event)
+
+	if diff := next.Sub(when); diff < -tolerance || diff > tolerance {
+		t.Fatalf("%s %s next mismatch: got %s want %s tol=%s hint=%s candidate=%s via=%s", event.Planet, event.Kind, next, when, tolerance, event.NAOJHintJST, event.CandidateJST, event.CandidateSource)
+	}
+	if diff := last.Sub(when); diff < -tolerance || diff > tolerance {
+		t.Fatalf("%s %s last mismatch: got %s want %s tol=%s hint=%s candidate=%s via=%s", event.Planet, event.Kind, last, when, tolerance, event.NAOJHintJST, event.CandidateJST, event.CandidateSource)
+	}
+}
+
+func TestOuterPlanetPhaseTruthAgainstJPL(t *testing.T) {
+	baseline := loadOuterTruthBaseline(t)
+	for _, event := range baseline.Events {
+		event := event
+		switch event.Kind {
+		case "P2R", "R2P":
+			// Station rows are retained as JPL apparent-RA reference data for
+			// future refinement. Current station behavior is constrained by the
+			// library's existing station baseline instead of these reference rows.
+			continue
+		}
+		name := strings.Join([]string{event.Planet, event.Kind, event.VerifiedJST}, "_")
+		t.Run(name, func(t *testing.T) {
+			lastFn, nextFn := outerTruthEventFuncs(t, event)
+			assertOuterTruthBaselineEvent(t, event, lastFn, nextFn)
+		})
+	}
+}
+
+func TestOuterPlanetStationJPLReferenceLoaded(t *testing.T) {
+	baseline := loadOuterTruthBaseline(t)
+	count := 0
+	for _, event := range baseline.Events {
+		switch event.Kind {
+		case "P2R", "R2P":
+			count++
+		}
+	}
+	if count == 0 {
+		t.Fatal("missing outer station JPL reference rows")
+	}
+}

basic/planet_apparent.go

@@ -0,0 +1,84 @@
+package basic
+
+import (
+	"math"
+
+	"b612.me/astro/planet"
+	. "b612.me/astro/tools"
+)
+
+type planetGeocentricPosition struct {
+	x  float64
+	y  float64
+	z  float64
+	lo float64
+	bo float64
+}
+
+func planetHeliocentricXYZN(planetIndex int, jd float64, n int) (float64, float64, float64) {
+	l := planet.WherePlanetN(planetIndex, 0, jd, n)
+	b := planet.WherePlanetN(planetIndex, 1, jd, n)
+	r := planet.WherePlanetN(planetIndex, 2, jd, n)
+	return sphericalToRectangular(l, b, r)
+}
+
+func earthHeliocentricXYZN(jd float64, n int) (float64, float64, float64) {
+	l := planet.WherePlanetN(-1, 0, jd, n)
+	b := planet.WherePlanetN(-1, 1, jd, n)
+	r := planet.WherePlanetN(-1, 2, jd, n)
+	return sphericalToRectangular(l, b, r)
+}
+
+func sphericalToRectangular(lo, bo, radius float64) (float64, float64, float64) {
+	cosBo := math.Cos(bo * math.Pi / 180)
+	return radius * cosBo * math.Cos(lo*math.Pi/180),
+		radius * cosBo * math.Sin(lo*math.Pi/180),
+		radius * math.Sin(bo*math.Pi/180)
+}
+
+func geocentricPositionFromRectangular(x, y, z float64) planetGeocentricPosition {
+	lo := math.Atan2(y, x) * 180 / math.Pi
+	bo := math.Atan2(z, math.Sqrt(x*x+y*y)) * 180 / math.Pi
+	return planetGeocentricPosition{
+		x:  x,
+		y:  y,
+		z:  z,
+		lo: Limit360(lo),
+		bo: bo,
+	}
+}
+
+func planetGeocentricPositionN(planetIndex int, planetJD, earthJD float64, n int) planetGeocentricPosition {
+	px, py, pz := planetHeliocentricXYZN(planetIndex, planetJD, n)
+	ex, ey, ez := earthHeliocentricXYZN(earthJD, n)
+	return geocentricPositionFromRectangular(px-ex, py-ey, pz-ez)
+}
+
+func planetGeocentricPositionWithEarthN(planetIndex int, planetJD float64, ex, ey, ez float64, n int) planetGeocentricPosition {
+	px, py, pz := planetHeliocentricXYZN(planetIndex, planetJD, n)
+	return geocentricPositionFromRectangular(px-ex, py-ey, pz-ez)
+}
+
+func planetApparentGeocentricPositionN(planetIndex int, jd float64, n int) (planetGeocentricPosition, float64) {
+	ex, ey, ez := earthHeliocentricXYZN(jd, n)
+	geoNow := planetGeocentricPositionWithEarthN(planetIndex, jd, ex, ey, ez, n)
+	tau := 0.0057755183 * math.Sqrt(geoNow.x*geoNow.x+geoNow.y*geoNow.y+geoNow.z*geoNow.z)
+	geo := planetGeocentricPositionWithEarthN(planetIndex, jd-tau, ex, ey, ez, n)
+	baseLo := geo.lo
+	baseBo := geo.bo
+	geo.lo = Limit360(baseLo + GXCLo(baseLo, baseBo, jd)/3600.0 + Nutation2000Bi(jd))
+	geo.bo = baseBo + GXCBo(baseLo, baseBo, jd)/3600.0
+	return geo, tau
+}
+
+func planetTrueGeocentricPositionN(planetIndex int, jd float64, n int) (planetGeocentricPosition, float64) {
+	ex, ey, ez := earthHeliocentricXYZN(jd, n)
+	geoNow := planetGeocentricPositionWithEarthN(planetIndex, jd, ex, ey, ez, n)
+	tau := 0.0057755183 * math.Sqrt(geoNow.x*geoNow.x+geoNow.y*geoNow.y+geoNow.z*geoNow.z)
+	return planetGeocentricPositionWithEarthN(planetIndex, jd-tau, ex, ey, ez, n), tau
+}
+
+func planetEarthAwayExplicitN(planetIndex int, jd float64, n int) float64 {
+	geoNow := planetGeocentricPositionN(planetIndex, jd, jd, n)
+	return math.Sqrt(geoNow.x*geoNow.x + geoNow.y*geoNow.y + geoNow.z*geoNow.z)
+}

basic/planet_apparent_external_test.go

@@ -0,0 +1,83 @@
+package basic
+
+import (
+	"encoding/json"
+	"os"
+	"testing"
+	"time"
+)
+
+type planetApparentSample struct {
+	Body              string  `json:"body"`
+	InputUTC          string  `json:"input_utc"`
+	RightAscension    float64 `json:"right_ascension"`
+	Declination       float64 `json:"declination"`
+	EclipticLongitude float64 `json:"ecliptic_longitude"`
+	EclipticLatitude  float64 `json:"ecliptic_latitude"`
+}
+
+func TestPlanetApparentCoordinatesMatchHorizonsBaseline(t *testing.T) {
+	data, err := os.ReadFile("testdata/planet_apparent_baseline.json")
+	if err != nil {
+		t.Fatalf("read baseline: %v", err)
+	}
+
+	var samples []planetApparentSample
+	if err := json.Unmarshal(data, &samples); err != nil {
+		t.Fatalf("decode baseline: %v", err)
+	}
+
+	type apparentCase struct {
+		lo  func(float64) float64
+		bo  func(float64) float64
+		ra  func(float64) float64
+		dec func(float64) float64
+	}
+
+	cases := map[string]apparentCase{
+		"mercury": {lo: MercuryApparentLo, bo: MercuryApparentBo, ra: MercuryApparentRa, dec: MercuryApparentDec},
+		"venus":   {lo: VenusApparentLo, bo: VenusApparentBo, ra: VenusApparentRa, dec: VenusApparentDec},
+		"mars":    {lo: MarsApparentLo, bo: MarsApparentBo, ra: MarsApparentRa, dec: MarsApparentDec},
+		"jupiter": {lo: JupiterApparentLo, bo: JupiterApparentBo, ra: JupiterApparentRa, dec: JupiterApparentDec},
+		"saturn":  {lo: SaturnApparentLo, bo: SaturnApparentBo, ra: SaturnApparentRa, dec: SaturnApparentDec},
+		"uranus":  {lo: UranusApparentLo, bo: UranusApparentBo, ra: UranusApparentRa, dec: UranusApparentDec},
+		"neptune": {lo: NeptuneApparentLo, bo: NeptuneApparentBo, ra: NeptuneApparentRa, dec: NeptuneApparentDec},
+	}
+
+	seen := make(map[string]bool, len(cases))
+	for _, sample := range samples {
+		tc, ok := cases[sample.Body]
+		if !ok {
+			t.Fatalf("unknown body %q", sample.Body)
+		}
+		if seen[sample.Body] {
+			t.Fatalf("duplicate body %q in apparent baseline", sample.Body)
+		}
+		seen[sample.Body] = true
+
+		date, err := time.Parse(time.RFC3339, sample.InputUTC)
+		if err != nil {
+			t.Fatalf("parse sample time %q: %v", sample.InputUTC, err)
+		}
+		jd := TD2UT(Date2JDE(date.UTC()), true)
+		prefix := sample.Body + "." + sample.InputUTC
+
+		assertPlanetApparentAngleClose(t, prefix+".RightAscension", tc.ra(jd), sample.RightAscension, 0.001)
+		assertPlanetPhaseClose(t, prefix+".Declination", tc.dec(jd), sample.Declination, 0.001)
+		assertPlanetApparentAngleClose(t, prefix+".EclipticLongitude", tc.lo(jd), sample.EclipticLongitude, 0.001)
+		assertPlanetPhaseClose(t, prefix+".EclipticLatitude", tc.bo(jd), sample.EclipticLatitude, 0.001)
+	}
+
+	for body := range cases {
+		if !seen[body] {
+			t.Fatalf("missing body %q in apparent baseline", body)
+		}
+	}
+}
+
+func assertPlanetApparentAngleClose(t *testing.T, name string, got, want, tolerance float64) {
+	t.Helper()
+	if diff := angleDiffAbs(got, want); diff > tolerance {
+		t.Fatalf("%s mismatch: got %.12f want %.12f diff %.12f tolerance %.12f", name, got, want, diff, tolerance)
+	}
+}

basic/planet_truncated.go

@@ -24,14 +24,8 @@ func planetXYZN(planetIndex int, jd float64, n int) (float64, float64, float64)
 }
 
 func planetApparentLoBoN(planetIndex int, jd float64, n int) (float64, float64) {
-	x, y, z := planetXYZN(planetIndex, jd, n)
-	to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
-	x, y, z = planetXYZN(planetIndex, jd-to, n)
-	lo := math.Atan2(y, x)
-	bo := math.Atan2(z, math.Sqrt(x*x+y*y))
-	lo = Limit360(lo*180/math.Pi) + Nutation2000Bi(jd)
-	bo = bo * 180 / math.Pi
-	return lo, bo
+	geo, _ := planetApparentGeocentricPositionN(planetIndex, jd, n)
+	return geo.lo, geo.bo
 }
 
 func planetApparentRaManualN(planetIndex int, jd float64, n int) float64 {
@@ -57,8 +51,7 @@ func planetApparentRaDecManualN(planetIndex int, jd float64, n int) (float64, fl
 }
 
 func planetEarthAwayN(planetIndex int, jd float64, n int) float64 {
-	x, y, z := planetXYZN(planetIndex, jd, n)
-	return math.Sqrt(x*x + y*y + z*z)
+	return planetEarthAwayExplicitN(planetIndex, jd, n)
 }
 
 func planetHeightN(jde, lon, lat, timezone float64, n int, apparentRaDec func(float64, int) (float64, float64)) float64 {

basic/saturn.go

@@ -87,53 +87,22 @@ func SaturnApparentRaDec(jd float64) (float64, float64) {
 }
 
 func EarthSaturnAway(jd float64) float64 {
-	x, y, z := ASaturnXYZ(jd)
-	to := math.Sqrt(x*x + y*y + z*z)
-	return to
+	return planetEarthAwayExplicitN(5, jd, -1)
 }
 
 func SaturnApparentLo(jd float64) float64 {
-	x, y, z := ASaturnXYZ(jd)
-	to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
-	x, y, z = ASaturnXYZ(jd - to)
-	lo := math.Atan2(y, x)
-	bo := math.Atan2(z, math.Sqrt(x*x+y*y))
-	lo = lo * 180 / math.Pi
-	bo = bo * 180 / math.Pi
-	lo = Limit360(lo)
-	//lo-=GXCLo(lo,bo,jd)/3600;
-	//bo+=GXCBo(lo,bo,jd);
-	lo += Nutation2000Bi(jd)
-	return lo
+	geo, _ := planetApparentGeocentricPositionN(5, jd, -1)
+	return geo.lo
 }
 
 func SaturnApparentBo(jd float64) float64 {
-	x, y, z := ASaturnXYZ(jd)
-	to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
-	x, y, z = ASaturnXYZ(jd - to)
-	//lo := math.Atan2(y, x)
-	bo := math.Atan2(z, math.Sqrt(x*x+y*y))
-	//lo = lo * 180 / math.Pi
-	bo = bo * 180 / math.Pi
-	//lo+=GXCLo(lo,bo,jd);
-	//bo+=GXCBo(lo,bo,jd)/3600;
-	//lo+=Nutation2000Bi(jd);
-	return bo
+	geo, _ := planetApparentGeocentricPositionN(5, jd, -1)
+	return geo.bo
 }
 
 func SaturnApparentLoBo(jd float64) (float64, float64) {
-	x, y, z := ASaturnXYZ(jd)
-	to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
-	x, y, z = ASaturnXYZ(jd - to)
-	lo := math.Atan2(y, x)
-	bo := math.Atan2(z, math.Sqrt(x*x+y*y))
-	lo = lo * 180 / math.Pi
-	bo = bo * 180 / math.Pi
-	lo = Limit360(lo)
-	//lo-=GXCLo(lo,bo,jd)/3600;
-	//bo+=GXCBo(lo,bo,jd);
-	lo += Nutation2000Bi(jd)
-	return lo, bo
+	geo, _ := planetApparentGeocentricPositionN(5, jd, -1)
+	return geo.lo, geo.bo
 }
 
 func SaturnMag(jd float64) float64 {

basic/saturn_events.go

@@ -40,6 +40,28 @@ func saturnSunLongitudeDeltaN(jde, degree float64, filter bool, n int) float64 {
 	return sub
 }
 
+func saturnRADerivative(jde, delta float64) float64 {
+	sub := SaturnApparentRa(jde+delta) - SaturnApparentRa(jde-delta)
+	if sub > 180 {
+		sub -= 360
+	}
+	if sub < -180 {
+		sub += 360
+	}
+	return sub / (2 * delta)
+}
+
+func saturnRADerivativeN(jde, delta float64, n int) float64 {
+	sub := SaturnApparentRaN(jde+delta, n) - SaturnApparentRaN(jde-delta, n)
+	if sub > 180 {
+		sub -= 360
+	}
+	if sub < -180 {
+		sub += 360
+	}
+	return sub / (2 * delta)
+}
+
 func saturnConjunctionFull(jde, degree float64, next uint8) float64 {
 	//0=last 1=next
 	daysPerDegree := SATURN_S_PERIOD / 360
@@ -94,113 +116,92 @@ func saturnConjunction(jde, degree float64, next uint8) float64 {
 }
 
 func LastSaturnConjunction(jde float64) float64 {
-	return saturnConjunction(jde, 0, 0)
+	return inclusiveLastPhaseEvent(jde, 0, saturnConjunction)
 }
 
 func NextSaturnConjunction(jde float64) float64 {
-	return saturnConjunction(jde, 0, 1)
+	return inclusiveNextPhaseEvent(jde, 0, saturnConjunction)
 }
 
 func LastSaturnOpposition(jde float64) float64 {
-	return saturnConjunction(jde, 180, 0)
+	return inclusiveLastPhaseEvent(jde, 180, saturnConjunction)
 }
 
 func NextSaturnOpposition(jde float64) float64 {
-	return saturnConjunction(jde, 180, 1)
+	return inclusiveNextPhaseEvent(jde, 180, saturnConjunction)
 }
 
 func NextSaturnEasternQuadrature(jde float64) float64 {
-	return saturnConjunction(jde, 90, 1)
+	return inclusiveNextPhaseEvent(jde, 90, saturnConjunction)
 }
 
 func LastSaturnEasternQuadrature(jde float64) float64 {
-	return saturnConjunction(jde, 90, 0)
+	return inclusiveLastPhaseEvent(jde, 90, saturnConjunction)
 }
 
 func NextSaturnWesternQuadrature(jde float64) float64 {
-	return saturnConjunction(jde, 270, 1)
+	return inclusiveNextPhaseEvent(jde, 270, saturnConjunction)
 }
 
 func LastSaturnWesternQuadrature(jde float64) float64 {
-	return saturnConjunction(jde, 270, 0)
+	return inclusiveLastPhaseEvent(jde, 270, saturnConjunction)
 }
 
-func saturnRetrograde(jde float64, searchBeforeOpposition bool) float64 {
-	//0=last 1=next
-	raRate := func(jde float64, delta float64) float64 {
-		sub := SaturnApparentRa(jde+delta) - SaturnApparentRa(jde-delta)
-		if sub > 180 {
-			sub -= 360
-		}
-		if sub < -180 {
-			sub += 360
-		}
-		return sub / (2 * delta)
-	}
-	jde = saturnConjunctionFull(jde, 180, 1)
+func saturnRetrogradeAroundOpposition(oppositionJD float64, searchBeforeOpposition bool) float64 {
+	oppositionTT := TD2UT(oppositionJD, true)
+	startTT := oppositionTT
+	endTT := oppositionTT
 	if searchBeforeOpposition {
-		jde -= 60
+		easternQuadratureUT := saturnConjunction(oppositionTT, 90, 0)
+		startTT = TD2UT(easternQuadratureUT, true)
 	} else {
-		jde += 60
+		westernQuadratureUT := saturnConjunction(oppositionTT, 270, 1)
+		endTT = TD2UT(westernQuadratureUT, true)
 	}
-	for {
-		currentRate := raRate(jde, 1.0/86400.0)
-		if math.Abs(currentRate) > 0.55 {
-			jde += 2
-			continue
-		}
-		break
-	}
-	estimateJD := jde
-	for {
-		prevJD := estimateJD
-		rateValue := raRate(prevJD, 2.0/86400.0)
-		rateSlope := (raRate(prevJD+15.0/86400.0, 2.0/86400.0) - raRate(prevJD-15.0/86400.0, 2.0/86400.0)) / (30.0 / 86400.0)
-		estimateJD = prevJD - rateValue/rateSlope
-		if math.Abs(estimateJD-prevJD) <= 30.0/86400.0 {
-			break
-		}
-	}
-	bestJD := eventZeroRefine(estimateJD, 15.0/86400.0, 0.5/86400.0, func(jd float64) float64 {
-		return raRate(jd, 0.5/86400.0)
+	bestJD := zeroEventInWindow(startTT, endTT, 2.0, 2.0, 30.0/86400.0, func(jd float64) float64 {
+		return saturnRADerivativeN(jd, 1.0/86400.0, saturnEventSearchN)
+	}, func(jd float64) float64 {
+		return saturnRADerivative(jd, 0.5/86400.0)
 	})
 	return TD2UT(bestJD, false)
 }
 
 func NextSaturnRetrogradeToPrograde(jde float64) float64 {
-	date := saturnRetrograde(jde, false)
-	if date < jde {
-		oppositionJD := saturnConjunctionFull(jde, 180, 1)
-		return saturnRetrograde(oppositionJD+10, false)
+	lastOppositionJD := saturnConjunctionFull(jde, 180, 0)
+	date := saturnRetrogradeAroundOpposition(lastOppositionJD, false)
+	if sameEventUTQueryTT(date, jde) || eventUTQueryAfterOrEqual(date, jde) {
+		return date
 	}
-	return date
+	nextOppositionJD := saturnConjunctionFull(jde, 180, 1)
+	return saturnRetrogradeAroundOpposition(nextOppositionJD, false)
 }
 
 func LastSaturnRetrogradeToPrograde(jde float64) float64 {
-	jde = saturnConjunctionFull(jde, 180, 0) - 10
-	date := saturnRetrograde(jde, false)
-	if date > jde {
-		oppositionJD := saturnConjunctionFull(jde, 180, 0)
-		return saturnRetrograde(oppositionJD-10, false)
+	lastOppositionJD := saturnConjunctionFull(jde, 180, 0)
+	date := saturnRetrogradeAroundOpposition(lastOppositionJD, false)
+	if sameEventUTQueryTT(date, jde) || eventUTQueryBeforeOrEqual(date, jde) {
+		return date
 	}
-	return date
+	previousOppositionJD := saturnConjunctionFull(eventUTLastQueryTT(lastOppositionJD), 180, 0)
+	return saturnRetrogradeAroundOpposition(previousOppositionJD, false)
 }
 
 func NextSaturnProgradeToRetrograde(jde float64) float64 {
-	date := saturnRetrograde(jde, true)
-	if date < jde {
-		oppositionJD := saturnConjunctionFull(jde, 180, 1)
-		return saturnRetrograde(oppositionJD+10, true)
+	nextOppositionJD := saturnConjunctionFull(jde, 180, 1)
+	date := saturnRetrogradeAroundOpposition(nextOppositionJD, true)
+	if sameEventUTQueryTT(date, jde) || eventUTQueryAfterOrEqual(date, jde) {
+		return date
 	}
-	return date
+	followingOppositionJD := saturnConjunctionFull(eventUTNextQueryTT(nextOppositionJD), 180, 1)
+	return saturnRetrogradeAroundOpposition(followingOppositionJD, true)
 }
 
 func LastSaturnProgradeToRetrograde(jde float64) float64 {
-	jde = saturnConjunctionFull(jde, 180, 0) - 10
-	date := saturnRetrograde(jde, true)
-	if date > jde {
-		oppositionJD := saturnConjunctionFull(jde, 180, 0)
-		return saturnRetrograde(oppositionJD-10, true)
+	nextOppositionJD := saturnConjunctionFull(jde, 180, 1)
+	date := saturnRetrogradeAroundOpposition(nextOppositionJD, true)
+	if sameEventUTQueryTT(date, jde) || eventUTQueryBeforeOrEqual(date, jde) {
+		return date
 	}
-	return date
+	lastOppositionJD := saturnConjunctionFull(jde, 180, 0)
+	return saturnRetrogradeAroundOpposition(lastOppositionJD, true)
 }

basic/station_truth_test.go

@@ -0,0 +1,207 @@
+package basic
+
+import (
+	"testing"
+	"time"
+)
+
+type stationEvent struct {
+	when time.Time
+	kind string
+}
+
+type stationTruthCase struct {
+	name    string
+	events  []stationEvent
+	lastR2P func(float64) float64
+	nextR2P func(float64) float64
+	lastP2R func(float64) float64
+	nextP2R func(float64) float64
+}
+
+func mustJST(value string) time.Time {
+	loc := time.FixedZone("JST", 9*3600)
+	t, err := time.ParseInLocation("2006-01-02 15:04", value, loc)
+	if err != nil {
+		panic(err)
+	}
+	return t
+}
+
+func toUTJD(t time.Time) float64 {
+	return TD2UT(Date2JDE(t.UTC()), true)
+}
+
+func TestStationTruthAgainstNAOJ(t *testing.T) {
+	cases := []stationTruthCase{
+		{
+			name: "Mars",
+			events: []stationEvent{
+				{when: mustJST("2024-12-08 05:59"), kind: "P2R"},
+				{when: mustJST("2025-01-16 11:39"), kind: "OPP"},
+				{when: mustJST("2025-02-24 18:35"), kind: "R2P"},
+				{when: mustJST("2027-01-12 01:10"), kind: "P2R"},
+				{when: mustJST("2027-02-20 00:51"), kind: "OPP"},
+				{when: mustJST("2027-04-03 02:33"), kind: "R2P"},
+			},
+			lastR2P: LastMarsRetrogradeToPrograde,
+			nextR2P: NextMarsRetrogradeToPrograde,
+			lastP2R: LastMarsProgradeToRetrograde,
+			nextP2R: NextMarsProgradeToRetrograde,
+		},
+		{
+			name: "Jupiter",
+			events: []stationEvent{
+				{when: mustJST("2024-10-09 16:13"), kind: "P2R"},
+				{when: mustJST("2024-12-08 05:58"), kind: "OPP"},
+				{when: mustJST("2025-02-04 22:07"), kind: "R2P"},
+				{when: mustJST("2025-11-12 04:54"), kind: "P2R"},
+				{when: mustJST("2026-01-10 17:42"), kind: "OPP"},
+				{when: mustJST("2026-03-11 11:44"), kind: "R2P"},
+				{when: mustJST("2026-12-13 21:03"), kind: "P2R"},
+				{when: mustJST("2027-02-11 09:29"), kind: "OPP"},
+				{when: mustJST("2027-04-13 15:17"), kind: "R2P"},
+			},
+			lastR2P: LastJupiterRetrogradeToPrograde,
+			nextR2P: NextJupiterRetrogradeToPrograde,
+			lastP2R: LastJupiterProgradeToRetrograde,
+			nextP2R: NextJupiterProgradeToRetrograde,
+		},
+		{
+			name: "Saturn",
+			events: []stationEvent{
+				{when: mustJST("2024-07-01 06:15"), kind: "P2R"},
+				{when: mustJST("2024-09-08 13:35"), kind: "OPP"},
+				{when: mustJST("2024-11-16 14:57"), kind: "R2P"},
+				{when: mustJST("2025-07-14 16:57"), kind: "P2R"},
+				{when: mustJST("2025-09-21 14:46"), kind: "OPP"},
+				{when: mustJST("2025-11-29 09:35"), kind: "R2P"},
+				{when: mustJST("2026-07-28 08:09"), kind: "P2R"},
+				{when: mustJST("2026-10-04 21:29"), kind: "OPP"},
+				{when: mustJST("2026-12-12 08:21"), kind: "R2P"},
+				{when: mustJST("2027-08-11 02:53"), kind: "P2R"},
+				{when: mustJST("2027-10-18 09:36"), kind: "OPP"},
+				{when: mustJST("2027-12-25 12:05"), kind: "R2P"},
+			},
+			lastR2P: LastSaturnRetrogradeToPrograde,
+			nextR2P: NextSaturnRetrogradeToPrograde,
+			lastP2R: LastSaturnProgradeToRetrograde,
+			nextP2R: NextSaturnProgradeToRetrograde,
+		},
+		{
+			name: "Uranus",
+			events: []stationEvent{
+				{when: mustJST("2024-01-27 19:50"), kind: "R2P"},
+				{when: mustJST("2024-09-02 00:44"), kind: "P2R"},
+				{when: mustJST("2024-11-17 11:45"), kind: "OPP"},
+				{when: mustJST("2025-01-31 04:04"), kind: "R2P"},
+				{when: mustJST("2025-09-06 13:55"), kind: "P2R"},
+				{when: mustJST("2025-11-21 21:25"), kind: "OPP"},
+				{when: mustJST("2026-02-04 13:37"), kind: "R2P"},
+				{when: mustJST("2026-09-11 03:19"), kind: "P2R"},
+				{when: mustJST("2026-11-26 07:41"), kind: "OPP"},
+				{when: mustJST("2027-02-08 23:03"), kind: "R2P"},
+				{when: mustJST("2027-09-15 17:50"), kind: "P2R"},
+				{when: mustJST("2027-11-30 18:22"), kind: "OPP"},
+			},
+			lastR2P: LastUranusRetrogradeToPrograde,
+			nextR2P: NextUranusRetrogradeToPrograde,
+			lastP2R: LastUranusProgradeToRetrograde,
+			nextP2R: NextUranusProgradeToRetrograde,
+		},
+		{
+			name: "Neptune",
+			events: []stationEvent{
+				{when: mustJST("2024-07-03 12:08"), kind: "P2R"},
+				{when: mustJST("2024-09-21 09:17"), kind: "OPP"},
+				{when: mustJST("2024-12-08 20:05"), kind: "R2P"},
+				{when: mustJST("2025-07-05 23:30"), kind: "P2R"},
+				{when: mustJST("2025-09-23 21:54"), kind: "OPP"},
+				{when: mustJST("2025-12-11 09:21"), kind: "R2P"},
+				{when: mustJST("2026-07-08 13:02"), kind: "P2R"},
+				{when: mustJST("2026-09-26 10:36"), kind: "OPP"},
+				{when: mustJST("2026-12-13 19:47"), kind: "R2P"},
+				{when: mustJST("2027-07-11 01:06"), kind: "P2R"},
+				{when: mustJST("2027-09-28 23:19"), kind: "OPP"},
+				{when: mustJST("2027-12-16 07:16"), kind: "R2P"},
+			},
+			lastR2P: LastNeptuneRetrogradeToPrograde,
+			nextR2P: NextNeptuneRetrogradeToPrograde,
+			lastP2R: LastNeptuneProgradeToRetrograde,
+			nextP2R: NextNeptuneProgradeToRetrograde,
+		},
+	}
+
+	for _, tc := range cases {
+		tc := tc
+		t.Run(tc.name, func(t *testing.T) {
+			for i, event := range tc.events {
+				switch event.kind {
+				case "P2R":
+					before := event.when.Add(-24 * time.Hour)
+					after := event.when.Add(24 * time.Hour)
+					nextP2R := JDE2DateByZone(tc.nextP2R(toUTJD(before)), event.when.Location(), false)
+					lastP2R := JDE2DateByZone(tc.lastP2R(toUTJD(after)), event.when.Location(), false)
+					if !sameMinute(nextP2R, event.when) {
+						t.Fatalf("%s next P2R mismatch: got %s want %s", tc.name, nextP2R, event.when)
+					}
+					if !sameMinute(lastP2R, event.when) {
+						t.Fatalf("%s last P2R mismatch: got %s want %s", tc.name, lastP2R, event.when)
+					}
+				case "R2P":
+					before := event.when.Add(-24 * time.Hour)
+					after := event.when.Add(24 * time.Hour)
+					nextR2P := JDE2DateByZone(tc.nextR2P(toUTJD(before)), event.when.Location(), false)
+					lastR2P := JDE2DateByZone(tc.lastR2P(toUTJD(after)), event.when.Location(), false)
+					if !sameMinute(nextR2P, event.when) {
+						t.Fatalf("%s next R2P mismatch: got %s want %s", tc.name, nextR2P, event.when)
+					}
+					if !sameMinute(lastR2P, event.when) {
+						t.Fatalf("%s last R2P mismatch: got %s want %s", tc.name, lastR2P, event.when)
+					}
+				case "OPP":
+					prev := nearestOfKindBefore(tc.events, i, "P2R")
+					next := nearestOfKindAfter(tc.events, i, "R2P")
+					if prev.IsZero() || next.IsZero() {
+						continue
+					}
+					query := event.when
+					lastP2R := JDE2DateByZone(tc.lastP2R(toUTJD(query)), query.Location(), false)
+					nextR2P := JDE2DateByZone(tc.nextR2P(toUTJD(query)), query.Location(), false)
+					if !sameMinute(lastP2R, prev) {
+						t.Fatalf("%s opposition last P2R mismatch: got %s want %s", tc.name, lastP2R, prev)
+					}
+					if !sameMinute(nextR2P, next) {
+						t.Fatalf("%s opposition next R2P mismatch: got %s want %s", tc.name, nextR2P, next)
+					}
+				}
+			}
+		})
+	}
+}
+
+func nearestOfKindBefore(events []stationEvent, idx int, kind string) time.Time {
+	for i := idx - 1; i >= 0; i-- {
+		if events[i].kind == kind {
+			return events[i].when
+		}
+	}
+	return time.Time{}
+}
+
+func nearestOfKindAfter(events []stationEvent, idx int, kind string) time.Time {
+	for i := idx + 1; i < len(events); i++ {
+		if events[i].kind == kind {
+			return events[i].when
+		}
+	}
+	return time.Time{}
+}
+
+func sameMinute(got, want time.Time) bool {
+	diff := got.Sub(want)
+	if diff < 0 {
+		diff = -diff
+	}
+	return diff <= 2*time.Minute
+}

basic/testdata/planet_apparent_baseline.json

@@ -0,0 +1,58 @@
+[
+  {
+    "body": "mercury",
+    "input_utc": "2026-01-01T00:00:00Z",
+    "right_ascension": 268.524035973,
+    "declination": -24.001291126,
+    "ecliptic_longitude": 268.6516112,
+    "ecliptic_latitude": -0.5700521
+  },
+  {
+    "body": "venus",
+    "input_utc": "2026-01-01T00:00:00Z",
+    "right_ascension": 280.056455581,
+    "declination": -23.622404581,
+    "ecliptic_longitude": 279.2064734,
+    "ecliptic_latitude": -0.5050645
+  },
+  {
+    "body": "mars",
+    "input_utc": "2026-01-01T00:00:00Z",
+    "right_ascension": 283.879610807,
+    "declination": -23.720007274,
+    "ecliptic_longitude": 282.6881475,
+    "ecliptic_latitude": -0.8911035
+  },
+  {
+    "body": "jupiter",
+    "input_utc": "2026-01-01T00:00:00Z",
+    "right_ascension": 113.124332352,
+    "declination": 21.979135798,
+    "ecliptic_longitude": 111.3575894,
+    "ecliptic_latitude": 0.2391257
+  },
+  {
+    "body": "saturn",
+    "input_utc": "2026-01-01T00:00:00Z",
+    "right_ascension": 357.380959207,
+    "declination": -3.596394732,
+    "ecliptic_longitude": 356.1672313,
+    "ecliptic_latitude": -2.2587419
+  },
+  {
+    "body": "uranus",
+    "input_utc": "2026-01-01T00:00:00Z",
+    "right_ascension": 55.737099009,
+    "declination": 19.509648526,
+    "ecliptic_longitude": 57.9492508,
+    "ecliptic_latitude": -0.1975992
+  },
+  {
+    "body": "neptune",
+    "input_utc": "2026-01-01T00:00:00Z",
+    "right_ascension": 0.077612938,
+    "declination": -1.418610222,
+    "ecliptic_longitude": 359.5068407,
+    "ecliptic_latitude": -1.3324096
+  }
+]

basic/uranus.go

@@ -87,53 +87,22 @@ func UranusApparentRaDec(jd float64) (float64, float64) {
 }
 
 func EarthUranusAway(jd float64) float64 {
-	x, y, z := AUranusXYZ(jd)
-	to := math.Sqrt(x*x + y*y + z*z)
-	return to
+	return planetEarthAwayExplicitN(6, jd, -1)
 }
 
 func UranusApparentLo(jd float64) float64 {
-	x, y, z := AUranusXYZ(jd)
-	to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
-	x, y, z = AUranusXYZ(jd - to)
-	lo := math.Atan2(y, x)
-	bo := math.Atan2(z, math.Sqrt(x*x+y*y))
-	lo = lo * 180 / math.Pi
-	bo = bo * 180 / math.Pi
-	lo = Limit360(lo)
-	//lo-=GXCLo(lo,bo,jd)/3600;
-	//bo+=GXCBo(lo,bo,jd);
-	lo += Nutation2000Bi(jd)
-	return lo
+	geo, _ := planetApparentGeocentricPositionN(6, jd, -1)
+	return geo.lo
 }
 
 func UranusApparentBo(jd float64) float64 {
-	x, y, z := AUranusXYZ(jd)
-	to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
-	x, y, z = AUranusXYZ(jd - to)
-	//lo := math.Atan2(y, x)
-	bo := math.Atan2(z, math.Sqrt(x*x+y*y))
-	//lo = lo * 180 / math.Pi
-	bo = bo * 180 / math.Pi
-	//lo+=GXCLo(lo,bo,jd);
-	//bo+=GXCBo(lo,bo,jd)/3600;
-	//lo+=Nutation2000Bi(jd);
-	return bo
+	geo, _ := planetApparentGeocentricPositionN(6, jd, -1)
+	return geo.bo
 }
 
 func UranusApparentLoBo(jd float64) (float64, float64) {
-	x, y, z := AUranusXYZ(jd)
-	to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
-	x, y, z = AUranusXYZ(jd - to)
-	lo := math.Atan2(y, x)
-	bo := math.Atan2(z, math.Sqrt(x*x+y*y))
-	lo = lo * 180 / math.Pi
-	bo = bo * 180 / math.Pi
-	lo = Limit360(lo)
-	//lo-=GXCLo(lo,bo,jd)/3600;
-	//bo+=GXCBo(lo,bo,jd);
-	lo += Nutation2000Bi(jd)
-	return lo, bo
+	geo, _ := planetApparentGeocentricPositionN(6, jd, -1)
+	return geo.lo, geo.bo
 }
 
 func UranusMag(jd float64) float64 {

basic/uranus_events.go

@@ -40,6 +40,28 @@ func uranusSunLongitudeDeltaN(jde, degree float64, filter bool, n int) float64 {
 	return sub
 }
 
+func uranusRADerivative(jde, delta float64) float64 {
+	sub := UranusApparentRa(jde+delta) - UranusApparentRa(jde-delta)
+	if sub > 180 {
+		sub -= 360
+	}
+	if sub < -180 {
+		sub += 360
+	}
+	return sub / (2 * delta)
+}
+
+func uranusRADerivativeN(jde, delta float64, n int) float64 {
+	sub := UranusApparentRaN(jde+delta, n) - UranusApparentRaN(jde-delta, n)
+	if sub > 180 {
+		sub -= 360
+	}
+	if sub < -180 {
+		sub += 360
+	}
+	return sub / (2 * delta)
+}
+
 func uranusConjunctionFull(jde, degree float64, next uint8) float64 {
 	//0=last 1=next
 	daysPerDegree := URANUS_S_PERIOD / 360
@@ -94,113 +116,92 @@ func uranusConjunction(jde, degree float64, next uint8) float64 {
 }
 
 func LastUranusConjunction(jde float64) float64 {
-	return uranusConjunction(jde, 0, 0)
+	return inclusiveLastPhaseEvent(jde, 0, uranusConjunction)
 }
 
 func NextUranusConjunction(jde float64) float64 {
-	return uranusConjunction(jde, 0, 1)
+	return inclusiveNextPhaseEvent(jde, 0, uranusConjunction)
 }
 
 func LastUranusOpposition(jde float64) float64 {
-	return uranusConjunction(jde, 180, 0)
+	return inclusiveLastPhaseEvent(jde, 180, uranusConjunction)
 }
 
 func NextUranusOpposition(jde float64) float64 {
-	return uranusConjunction(jde, 180, 1)
+	return inclusiveNextPhaseEvent(jde, 180, uranusConjunction)
 }
 
 func NextUranusEasternQuadrature(jde float64) float64 {
-	return uranusConjunction(jde, 90, 1)
+	return inclusiveNextPhaseEvent(jde, 90, uranusConjunction)
 }
 
 func LastUranusEasternQuadrature(jde float64) float64 {
-	return uranusConjunction(jde, 90, 0)
+	return inclusiveLastPhaseEvent(jde, 90, uranusConjunction)
 }
 
 func NextUranusWesternQuadrature(jde float64) float64 {
-	return uranusConjunction(jde, 270, 1)
+	return inclusiveNextPhaseEvent(jde, 270, uranusConjunction)
 }
 
 func LastUranusWesternQuadrature(jde float64) float64 {
-	return uranusConjunction(jde, 270, 0)
+	return inclusiveLastPhaseEvent(jde, 270, uranusConjunction)
 }
 
-func uranusRetrograde(jde float64, searchBeforeOpposition bool) float64 {
-	//0=last 1=next
-	raRate := func(jde float64, delta float64) float64 {
-		sub := UranusApparentRa(jde+delta) - UranusApparentRa(jde-delta)
-		if sub > 180 {
-			sub -= 360
-		}
-		if sub < -180 {
-			sub += 360
-		}
-		return sub / (2 * delta)
-	}
-	jde = uranusConjunctionFull(jde, 180, 1)
+func uranusRetrogradeAroundOpposition(oppositionJD float64, searchBeforeOpposition bool) float64 {
+	oppositionTT := TD2UT(oppositionJD, true)
+	startTT := oppositionTT
+	endTT := oppositionTT
 	if searchBeforeOpposition {
-		jde -= 60
+		easternQuadratureUT := uranusConjunction(oppositionTT, 90, 0)
+		startTT = TD2UT(easternQuadratureUT, true)
 	} else {
-		jde += 60
+		westernQuadratureUT := uranusConjunction(oppositionTT, 270, 1)
+		endTT = TD2UT(westernQuadratureUT, true)
 	}
-	for {
-		currentRate := raRate(jde, 1.0/86400.0)
-		if math.Abs(currentRate) > 0.55 {
-			jde += 2
-			continue
-		}
-		break
-	}
-	estimateJD := jde
-	for {
-		prevJD := estimateJD
-		rateValue := raRate(prevJD, 2.0/86400.0)
-		rateSlope := (raRate(prevJD+15.0/86400.0, 2.0/86400.0) - raRate(prevJD-15.0/86400.0, 2.0/86400.0)) / (30.0 / 86400.0)
-		estimateJD = prevJD - rateValue/rateSlope
-		if math.Abs(estimateJD-prevJD) <= 30.0/86400.0 {
-			break
-		}
-	}
-	bestJD := eventZeroRefine(estimateJD, 15.0/86400.0, 0.5/86400.0, func(jd float64) float64 {
-		return raRate(jd, 0.5/86400.0)
+	bestJD := zeroEventInWindow(startTT, endTT, 2.0, 2.0, 30.0/86400.0, func(jd float64) float64 {
+		return uranusRADerivativeN(jd, 1.0/86400.0, uranusEventSearchN)
+	}, func(jd float64) float64 {
+		return uranusRADerivative(jd, 0.5/86400.0)
 	})
 	return TD2UT(bestJD, false)
 }
 
 func NextUranusRetrogradeToPrograde(jde float64) float64 {
-	date := uranusRetrograde(jde, false)
-	if date < jde {
-		oppositionJD := uranusConjunctionFull(jde, 180, 1)
-		return uranusRetrograde(oppositionJD+10, false)
+	lastOppositionJD := uranusConjunctionFull(jde, 180, 0)
+	date := uranusRetrogradeAroundOpposition(lastOppositionJD, false)
+	if sameEventUTQueryTT(date, jde) || eventUTQueryAfterOrEqual(date, jde) {
+		return date
 	}
-	return date
+	nextOppositionJD := uranusConjunctionFull(jde, 180, 1)
+	return uranusRetrogradeAroundOpposition(nextOppositionJD, false)
 }
 
 func LastUranusRetrogradeToPrograde(jde float64) float64 {
-	jde = uranusConjunctionFull(jde, 180, 0) - 10
-	date := uranusRetrograde(jde, false)
-	if date > jde {
-		oppositionJD := uranusConjunctionFull(jde, 180, 0)
-		return uranusRetrograde(oppositionJD-10, false)
+	lastOppositionJD := uranusConjunctionFull(jde, 180, 0)
+	date := uranusRetrogradeAroundOpposition(lastOppositionJD, false)
+	if sameEventUTQueryTT(date, jde) || eventUTQueryBeforeOrEqual(date, jde) {
+		return date
 	}
-	return date
+	previousOppositionJD := uranusConjunctionFull(eventUTLastQueryTT(lastOppositionJD), 180, 0)
+	return uranusRetrogradeAroundOpposition(previousOppositionJD, false)
 }
 
 func NextUranusProgradeToRetrograde(jde float64) float64 {
-	date := uranusRetrograde(jde, true)
-	if date < jde {
-		oppositionJD := uranusConjunctionFull(jde, 180, 1)
-		return uranusRetrograde(oppositionJD+10, true)
+	nextOppositionJD := uranusConjunctionFull(jde, 180, 1)
+	date := uranusRetrogradeAroundOpposition(nextOppositionJD, true)
+	if sameEventUTQueryTT(date, jde) || eventUTQueryAfterOrEqual(date, jde) {
+		return date
 	}
-	return date
+	followingOppositionJD := uranusConjunctionFull(eventUTNextQueryTT(nextOppositionJD), 180, 1)
+	return uranusRetrogradeAroundOpposition(followingOppositionJD, true)
 }
 
 func LastUranusProgradeToRetrograde(jde float64) float64 {
-	jde = uranusConjunctionFull(jde, 180, 0) - 10
-	date := uranusRetrograde(jde, true)
-	if date > jde {
-		oppositionJD := uranusConjunctionFull(jde, 180, 0)
-		return uranusRetrograde(oppositionJD-10, true)
+	nextOppositionJD := uranusConjunctionFull(jde, 180, 1)
+	date := uranusRetrogradeAroundOpposition(nextOppositionJD, true)
+	if sameEventUTQueryTT(date, jde) || eventUTQueryBeforeOrEqual(date, jde) {
+		return date
 	}
-	return date
+	lastOppositionJD := uranusConjunctionFull(jde, 180, 0)
+	return uranusRetrogradeAroundOpposition(lastOppositionJD, true)
 }

basic/venus.go

@@ -87,53 +87,22 @@ func VenusApparentRaDec(jd float64) (float64, float64) {
 }
 
 func EarthVenusAway(jd float64) float64 {
-	x, y, z := AVenusXYZ(jd)
-	to := math.Sqrt(x*x + y*y + z*z)
-	return to
+	return planetEarthAwayExplicitN(2, jd, -1)
 }
 
 func VenusApparentLo(jd float64) float64 {
-	x, y, z := AVenusXYZ(jd)
-	to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
-	x, y, z = AVenusXYZ(jd - to)
-	lo := math.Atan2(y, x)
-	bo := math.Atan2(z, math.Sqrt(x*x+y*y))
-	lo = lo * 180 / math.Pi
-	bo = bo * 180 / math.Pi
-	lo = Limit360(lo)
-	//lo-=GXCLo(lo,bo,jd)/3600;
-	//bo+=GXCBo(lo,bo,jd);
-	lo += Nutation2000Bi(jd)
-	return lo
+	geo, _ := planetApparentGeocentricPositionN(2, jd, -1)
+	return geo.lo
 }
 
 func VenusApparentBo(jd float64) float64 {
-	x, y, z := AVenusXYZ(jd)
-	to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
-	x, y, z = AVenusXYZ(jd - to)
-	//lo := math.Atan2(y, x)
-	bo := math.Atan2(z, math.Sqrt(x*x+y*y))
-	//lo = lo * 180 / math.Pi
-	bo = bo * 180 / math.Pi
-	//lo+=GXCLo(lo,bo,jd);
-	//bo+=GXCBo(lo,bo,jd)/3600;
-	//lo+=Nutation2000Bi(jd);
-	return bo
+	geo, _ := planetApparentGeocentricPositionN(2, jd, -1)
+	return geo.bo
 }
 
 func VenusApparentLoBo(jd float64) (float64, float64) {
-	x, y, z := AVenusXYZ(jd)
-	to := 0.0057755183 * math.Sqrt(x*x+y*y+z*z)
-	x, y, z = AVenusXYZ(jd - to)
-	lo := math.Atan2(y, x)
-	bo := math.Atan2(z, math.Sqrt(x*x+y*y))
-	lo = lo * 180 / math.Pi
-	bo = bo * 180 / math.Pi
-	lo = Limit360(lo)
-	//lo-=GXCLo(lo,bo,jd)/3600;
-	//bo+=GXCBo(lo,bo,jd);
-	lo += Nutation2000Bi(jd)
-	return lo, bo
+	geo, _ := planetApparentGeocentricPositionN(2, jd, -1)
+	return geo.lo, geo.bo
 }
 
 func VenusMag(jd float64) float64 {

basic/venus_events.go

@@ -3,6 +3,7 @@ package basic
 import (
 	"math"
 
+	"b612.me/astro/planet"
 	. "b612.me/astro/tools"
 )
 
@@ -33,6 +34,23 @@ func venusSunLongitudeDeltaN(jde float64, n int) float64 {
 	return sub
 }
 
+func venusConjunctionAngleDelta(diff float64) float64 {
+	diff = Limit360(diff)
+	if diff > 180 {
+		diff -= 360
+	}
+	if diff < -180 {
+		diff += 360
+	}
+	return diff
+}
+
+func venusConjunctionHeliocentricDelta(jd, targetDeg float64, n int) float64 {
+	planetLo := planet.WherePlanetN(2, 0, jd, n)
+	earthLo := planet.WherePlanetN(-1, 0, jd, n)
+	return venusConjunctionAngleDelta(planetLo - earthLo - targetDeg)
+}
+
 func venusSunRADelta(jde float64) float64 {
 	sub := Limit360(VenusApparentRa(jde) - SunApparentRa(jde))
 	if sub > 180 {
@@ -66,13 +84,94 @@ func venusRADerivativeN(jde, val float64, n int) float64 {
 	return sub / (2 * val)
 }
 
+func venusRAContinuousForMax(jde float64) float64 {
+	ra := VenusApparentRa(jde)
+	if ra < 180 {
+		return ra + 360
+	}
+	return ra
+}
+
+func venusRAContinuousForMaxN(jde float64, n int) float64 {
+	ra := VenusApparentRaN(jde, n)
+	if ra < 180 {
+		return ra + 360
+	}
+	return ra
+}
+
+func venusRAContinuousForMin(jde float64) float64 {
+	ra := VenusApparentRa(jde)
+	if ra > 180 {
+		return ra - 360
+	}
+	return ra
+}
+
+func venusRAContinuousForMinN(jde float64, n int) float64 {
+	ra := VenusApparentRaN(jde, n)
+	if ra > 180 {
+		return ra - 360
+	}
+	return ra
+}
+
+func venusRAExtremumRefine(seed, start, end, step float64, fn func(float64) float64) float64 {
+	centerJD := clampFloat64(seed, start, end)
+	halfStep := step
+	bestJD := centerJD
+	bestVal := fn(centerJD)
+	for i := 0; i < 8; i++ {
+		leftJD := clampFloat64(centerJD-halfStep, start, end)
+		rightJD := clampFloat64(centerJD+halfStep, start, end)
+		leftVal := fn(leftJD)
+		centerVal := fn(centerJD)
+		rightVal := fn(rightJD)
+		if leftVal > bestVal {
+			bestVal = leftVal
+			bestJD = leftJD
+		}
+		if centerVal > bestVal {
+			bestVal = centerVal
+			bestJD = centerJD
+		}
+		if rightVal > bestVal {
+			bestVal = rightVal
+			bestJD = rightJD
+		}
+		denominator := leftVal - 2*centerVal + rightVal
+		if denominator == 0 {
+			centerJD = bestJD
+			halfStep /= 2
+			continue
+		}
+		vertexJD := centerJD + 0.5*halfStep*(leftVal-rightVal)/denominator
+		vertexJD = clampFloat64(vertexJD, leftJD, rightJD)
+		vertexVal := fn(vertexJD)
+		if vertexVal > bestVal {
+			bestVal = vertexVal
+			bestJD = vertexJD
+		}
+		centerJD = bestJD
+		halfStep /= 2
+	}
+	return bestJD
+}
+
 func venusSunElongationN(jde float64, n int) float64 {
 	lo1, bo1 := VenusApparentLoBoN(jde, n)
-	lo2 := SunApparentLo(jde)
+	lo2 := HSunApparentLoN(jde, n)
 	bo2 := HSunTrueBoN(jde, n)
 	return StarAngularSeparation(lo1, bo1, lo2, bo2)
 }
 
+func venusTrueElongationN(jde float64, n int) float64 {
+	earth := mercuryHelioN(-1, jde, n)
+	planetPos := mercuryHelioN(2, jde, n)
+	geo := mercuryGeocentric(planetPos, earth)
+	return StarAngularSeparation(geo.lo, geo.bo, HSunTrueLoN(jde, n), HSunTrueBoN(jde, n))
+}
+
 func venusElongationDerivative(jde, val float64) float64 {
 	sub := VenusSunElongation(jde+val) - VenusSunElongation(jde-val)
 	if sub > 180 {
@@ -96,91 +195,98 @@ func venusElongationDerivativeN(jde, val float64, n int) float64 {
 }
 
 func venusConjunction(jde float64, next uint8) float64 {
-	//0=last 1=next
-	nowSub := venusSunLongitudeDeltaN(jde, venusEventSearchN)
-	pos := math.Abs(venusSunLongitudeDeltaN(jde+1/86400.0, venusEventSearchN)) - math.Abs(nowSub)
-	if pos >= 0 && next == 1 && nowSub > 0 {
-		jde += VENUS_S_PERIOD/8.0 + 2
+	queryTT := jde
+	direction := -1.0
+	if next == 1 {
+		direction = 1
 	}
-	if pos >= 0 && next == 1 && nowSub < 0 {
-		jde += VENUS_S_PERIOD/6.0 + 2
-	}
-	if pos <= 0 && next == 0 && nowSub < 0 {
-		jde -= VENUS_S_PERIOD/8.0 + 2
-	}
-	if pos <= 0 && next == 0 && nowSub > 0 {
-		jde -= VENUS_S_PERIOD/6.0 + 2
-	}
-	for {
-		nowSub := venusSunLongitudeDeltaN(jde, venusEventSearchN)
-		pos := math.Abs(venusSunLongitudeDeltaN(jde+1/86400.0, venusEventSearchN)) - math.Abs(nowSub)
-		if math.Abs(nowSub) > 24 || (pos > 0 && next == 1) || (pos < 0 && next == 0) {
-			if next == 1 {
-				jde += 8
-			} else {
-				jde -= 8
-			}
-			continue
-		}
-		break
-	}
-	JD1 := jde
-	for {
-		JD0 := JD1
-		stDegree := venusSunLongitudeDelta(JD0)
-		stDegreep := (venusSunLongitudeDelta(JD0+0.000005) - venusSunLongitudeDelta(JD0-0.000005)) / 0.00001
-		JD1 = JD0 - stDegree/stDegreep
-		if math.Abs(JD1-JD0) <= 0.00001 {
-			break
+	left := queryTT
+	leftVal := venusSunLongitudeDeltaN(left, venusEventSearchN)
+	if math.Abs(leftVal) <= 30.0/86400.0 {
+		exact := eventZeroRefine(left, 1.0, 0.000005, venusSunLongitudeDelta)
+		if math.Abs(exact-queryTT) <= 1.0 {
+			return TD2UT(exact, false)
 		}
 	}
-	return TD2UT(JD1, false)
+	const step = 8.0
+	for i := 0; i < 80; i++ {
+		right := queryTT + direction*step*float64(i+1)
+		rightVal := venusSunLongitudeDeltaN(right, venusEventSearchN)
+		if leftVal == 0 || rightVal == 0 || leftVal*rightVal <= 0 {
+			center := (left + right) / 2.0
+			halfWindow := math.Abs(right-left) / 2.0
+			return TD2UT(eventZeroRefine(center, halfWindow, 0.000005, venusSunLongitudeDelta), false)
+		}
+		left = right
+		leftVal = rightVal
+	}
+	return TD2UT(eventZeroRefine(queryTT, VENUS_S_PERIOD, 0.000005, venusSunLongitudeDelta), false)
+}
+
+func venusConjunctionTypeAt(eventUT float64) bool {
+	return EarthVenusAway(eventUT) <= EarthAway(eventUT)
+}
+
+func nextVenusTypedConjunctionFromEvent(jde float64, inferior bool) float64 {
+	date := NextVenusConjunctionStrict(jde)
+	if venusConjunctionTypeAt(date) == inferior {
+		return date
+	}
+	return NextVenusConjunctionStrict(eventUTNextQueryTT(date))
+}
+
+func lastVenusTypedConjunctionFromEvent(jde float64, inferior bool) float64 {
+	date := LastVenusConjunctionStrict(jde)
+	if venusConjunctionTypeAt(date) == inferior {
+		return date
+	}
+	return LastVenusConjunctionStrict(eventUTLastQueryTT(date))
 }
 
 func LastVenusConjunction(jde float64) float64 {
-	return venusConjunction(jde, 0)
+	return inclusiveLastSimpleEvent(jde, LastVenusConjunctionStrict, NextVenusConjunctionStrict)
 }
 
 func NextVenusConjunction(jde float64) float64 {
+	return inclusiveNextSimpleEvent(jde, LastVenusConjunctionStrict, NextVenusConjunctionStrict)
+}
+
+func LastVenusConjunctionStrict(jde float64) float64 {
+	return venusConjunction(jde, 0)
+}
+
+func NextVenusConjunctionStrict(jde float64) float64 {
 	return venusConjunction(jde, 1)
 }
 
+func nextVenusTypedConjunction(jde float64, inferior bool) float64 {
+	return nextVenusTypedConjunctionFromEvent(jde, inferior)
+}
+
+func lastVenusTypedConjunction(jde float64, inferior bool) float64 {
+	return lastVenusTypedConjunctionFromEvent(jde, inferior)
+}
+
 func NextVenusInferiorConjunction(jde float64) float64 {
-	date := NextVenusConjunction(jde)
-	if EarthVenusAway(date) > EarthAway(date) {
-		return NextVenusConjunction(date + 2)
-	}
-	return date
+	return nextVenusTypedConjunction(jde, true)
 }
 
 func NextVenusSuperiorConjunction(jde float64) float64 {
-	date := NextVenusConjunction(jde)
-	if EarthVenusAway(date) < EarthAway(date) {
-		return NextVenusConjunction(date + 2)
-	}
-	return date
+	return nextVenusTypedConjunction(jde, false)
 }
 
 func LastVenusInferiorConjunction(jde float64) float64 {
-	date := LastVenusConjunction(jde)
-	if EarthVenusAway(date) > EarthAway(date) {
-		return LastVenusConjunction(date - 2)
-	}
-	return date
+	return lastVenusTypedConjunction(jde, true)
 }
 
 func LastVenusSuperiorConjunction(jde float64) float64 {
-	date := LastVenusConjunction(jde)
-	if EarthVenusAway(date) < EarthAway(date) {
-		return LastVenusConjunction(date - 2)
-	}
-	return date
+	return lastVenusTypedConjunction(jde, false)
 }
 
 func venusRetrograde(jde float64) float64 {
 	//0=last 1=next
-	lastHe := LastVenusConjunction(jde)
-	nextHe := NextVenusConjunction(jde)
+	lastHe := LastVenusConjunctionStrict(jde)
+	nextHe := NextVenusConjunctionStrict(jde)
 	nowSub := venusSunRADelta(jde)
 	if nowSub > 0 {
 		jde = lastHe + ((nextHe - lastHe) / 5.0 * 3.5)
@@ -213,152 +319,317 @@ func venusRetrograde(jde float64) float64 {
 }
 
 func NextVenusRetrograde(jde float64) float64 {
-	date := venusRetrograde(jde)
-	if date < jde {
-		nextHe := NextVenusConjunction(jde)
-		return venusRetrograde(nextHe + 2)
+	p2r := NextVenusProgradeToRetrograde(jde)
+	r2p := NextVenusRetrogradeToPrograde(jde)
+	if sameEventJD(p2r, r2p) {
+		return p2r
 	}
-	return date
+	if p2r < r2p {
+		return p2r
+	}
+	return r2p
 }
 
 func LastVenusRetrograde(jde float64) float64 {
-	lastHe := LastVenusConjunction(jde)
-	date := venusRetrograde(lastHe + 2)
-	if date > jde {
-		lastLastHe := LastVenusConjunction(lastHe - 2)
-		return venusRetrograde(lastLastHe + 2)
+	p2r := LastVenusProgradeToRetrograde(jde)
+	r2p := LastVenusRetrogradeToPrograde(jde)
+	if sameEventJD(p2r, r2p) {
+		return p2r
 	}
-	return date
+	if p2r > r2p {
+		return p2r
+	}
+	return r2p
+}
+
+func venusStationInWindow(start, end float64, progradeToRetrograde bool) float64 {
+	var best float64
+	if progradeToRetrograde {
+		guess := scanWindowForMax(start, end, 2.0, func(jd float64) float64 {
+			return venusRAContinuousForMaxN(jd, venusEventSearchN)
+		})
+		best = venusRAExtremumRefine(guess, start, end, 1.0, func(jd float64) float64 {
+			return venusRAContinuousForMax(jd)
+		})
+	} else {
+		guess := scanWindowForMax(start, end, 2.0, func(jd float64) float64 {
+			return -venusRAContinuousForMinN(jd, venusEventSearchN)
+		})
+		best = venusRAExtremumRefine(guess, start, end, 1.0, func(jd float64) float64 {
+			return -venusRAContinuousForMin(jd)
+		})
+	}
+	return TD2UT(best, false)
+}
+
+func venusProgradeToRetrogradeAroundInferior(inferior float64) float64 {
+	return venusStationInWindow(inferior-30.0, inferior-14.0, true)
+}
+
+func venusRetrogradeToProgradeAroundInferior(inferior float64) float64 {
+	return venusStationInWindow(inferior+14.0, inferior+24.0, false)
 }
 
 func NextVenusProgradeToRetrograde(jde float64) float64 {
-	date := NextVenusRetrograde(jde)
-	sub := Limit360(VenusApparentRa(date) - SunApparentRa(date))
-	if sub > 180 {
-		return NextVenusRetrograde(date + VENUS_S_PERIOD/2)
+	inferior := NextVenusInferiorConjunction(jde)
+	for {
+		date := venusProgradeToRetrogradeAroundInferior(inferior)
+		if eventUTQueryAfterOrEqual(date, jde) {
+			return date
+		}
+		inferior = NextVenusInferiorConjunction(eventUTNextQueryTT(inferior))
 	}
-	return date
 }
 
 func NextVenusRetrogradeToPrograde(jde float64) float64 {
-	date := NextVenusRetrograde(jde)
-	sub := Limit360(VenusApparentRa(date) - SunApparentRa(date))
-	if sub < 180 {
-		return NextVenusRetrograde(date + 12)
+	inferior := LastVenusInferiorConjunction(jde)
+	for {
+		date := venusRetrogradeToProgradeAroundInferior(inferior)
+		if eventUTQueryAfterOrEqual(date, jde) {
+			return date
+		}
+		inferior = NextVenusInferiorConjunction(eventUTNextQueryTT(inferior))
 	}
-	return date
 }
 
 func LastVenusProgradeToRetrograde(jde float64) float64 {
-	date := LastVenusRetrograde(jde)
-	sub := Limit360(VenusApparentRa(date) - SunApparentRa(date))
-	if sub > 180 {
-		return LastVenusRetrograde(date - 12)
+	inferior := NextVenusInferiorConjunction(jde)
+	for {
+		date := venusProgradeToRetrogradeAroundInferior(inferior)
+		if eventUTQueryBeforeOrEqual(date, jde) {
+			return date
+		}
+		inferior = LastVenusInferiorConjunction(eventUTLastQueryTT(inferior))
 	}
-	return date
 }
 
 func LastVenusRetrogradeToPrograde(jde float64) float64 {
-	date := LastVenusRetrograde(jde)
-	sub := Limit360(VenusApparentRa(date) - SunApparentRa(date))
-	if sub < 180 {
-		return LastVenusRetrograde(date - VENUS_S_PERIOD/2)
+	inferior := LastVenusInferiorConjunction(jde)
+	for {
+		date := venusRetrogradeToProgradeAroundInferior(inferior)
+		if eventUTQueryBeforeOrEqual(date, jde) {
+			return date
+		}
+		inferior = LastVenusInferiorConjunction(eventUTLastQueryTT(inferior))
 	}
-	return date
 }
 
 func VenusSunElongation(jde float64) float64 {
 	lo1, bo1 := VenusApparentLoBo(jde)
-	lo2 := SunApparentLo(jde)
+	lo2 := HSunApparentLo(jde)
 	bo2 := HSunTrueBo(jde)
 	return StarAngularSeparation(lo1, bo1, lo2, bo2)
 }
 
-func venusGreatestElongation(jde float64) float64 {
-	lastHe := LastVenusConjunction(jde)
-	nextHe := NextVenusConjunction(jde)
-	nowSub := venusSunRADelta(jde)
-	if nowSub > 0 {
-		jde = lastHe + ((nextHe - lastHe) / 5.0 * 2.5)
-	} else {
-		jde = lastHe + ((nextHe - lastHe) / 5.0)
-	}
-	for {
-		nowSub := venusElongationDerivativeN(jde, 1.0/86400.0, venusEventSearchN)
-		if math.Abs(nowSub) > 0.15 {
-			jde += 5
-			continue
-		}
-		break
-	}
-	JD1 := jde
-	for {
-		JD0 := JD1
-		stDegree := venusElongationDerivative(JD0, 2.0/86400.0)
-		stDegreep := (venusElongationDerivative(JD0+15.0/86400.0, 2.0/86400.0) - venusElongationDerivative(JD0-15.0/86400.0, 2.0/86400.0)) / (30.0 / 86400.0)
-		JD1 = JD0 - stDegree/stDegreep
-		if math.Abs(JD1-JD0) <= 30.0/86400.0 {
-			break
-		}
-	}
-	min := eventZeroRefine(JD1, 15.0/86400.0, 0.5/86400.0, func(jd float64) float64 {
-		return venusElongationDerivative(jd, 0.5/86400.0)
+func venusGreatestElongationInWindow(start, end float64) float64 {
+	best := maximizeInWindow(start, end, 5.0, func(jd float64) float64 {
+		return venusTrueElongationN(jd, venusEventSearchN)
+	}, func(jd float64) float64 {
+		return venusTrueElongationN(jd, -1)
 	})
-	//fmt.Println((min - lastHe) / (nextHe - lastHe))
-	return TD2UT(min, false)
+	return TD2UT(best, false)
+}
+
+func venusEastElongationWindowEndingAt(inferior float64) (float64, float64) {
+	lastSuperior := LastVenusSuperiorConjunction(eventUTLastQueryTT(inferior))
+	return lastSuperior + innerEventEpsilon, inferior - innerEventEpsilon
+}
+
+func venusWestElongationWindowEndingAt(superior float64) (float64, float64) {
+	lastInferior := LastVenusInferiorConjunction(eventUTLastQueryTT(superior))
+	return lastInferior + innerEventEpsilon, superior - innerEventEpsilon
+}
+
+func venusEastElongationWindowContaining(jde float64) (float64, float64) {
+	nextInferior := NextVenusInferiorConjunction(jde)
+	start, end := venusEastElongationWindowEndingAt(nextInferior)
+	if eventUTQueryBeforeOrEqual(start, jde) && eventUTQueryAfterOrEqual(end, jde) {
+		return start, end
+	}
+	currentInferior := LastVenusInferiorConjunction(jde)
+	return venusEastElongationWindowEndingAt(currentInferior)
+}
+
+func venusWestElongationWindowContaining(jde float64) (float64, float64) {
+	nextSuperior := NextVenusSuperiorConjunction(jde)
+	start, end := venusWestElongationWindowEndingAt(nextSuperior)
+	if eventUTQueryBeforeOrEqual(start, jde) && eventUTQueryAfterOrEqual(end, jde) {
+		return start, end
+	}
+	currentSuperior := LastVenusSuperiorConjunction(jde)
+	return venusWestElongationWindowEndingAt(currentSuperior)
+}
+
+func nextVenusGreatestElongationTyped(jde float64, east bool) float64 {
+	if east {
+		start, windowEnd := venusEastElongationWindowContaining(jde)
+		for {
+			date := venusGreatestElongationInWindow(start, windowEnd)
+			if eventUTQueryAfterOrEqual(date, jde) {
+				return date
+			}
+			nextInferior := NextVenusInferiorConjunction(eventUTNextQueryTT(windowEnd))
+			start, windowEnd = venusEastElongationWindowEndingAt(nextInferior)
+		}
+	}
+	start, windowEnd := venusWestElongationWindowContaining(jde)
+	for {
+		date := venusGreatestElongationInWindow(start, windowEnd)
+		if eventUTQueryAfterOrEqual(date, jde) {
+			return date
+		}
+		nextSuperior := NextVenusSuperiorConjunction(eventUTNextQueryTT(windowEnd))
+		start, windowEnd = venusWestElongationWindowEndingAt(nextSuperior)
+	}
+}
+
+func lastVenusGreatestElongationTyped(jde float64, east bool) float64 {
+	if east {
+		start, windowEnd := venusEastElongationWindowContaining(jde)
+		for {
+			date := venusGreatestElongationInWindow(start, windowEnd)
+			if eventUTQueryBeforeOrEqual(date, jde) {
+				return date
+			}
+			prevInferior := LastVenusInferiorConjunction(eventUTLastQueryTT(start))
+			start, windowEnd = venusEastElongationWindowEndingAt(prevInferior)
+		}
+	}
+	start, windowEnd := venusWestElongationWindowContaining(jde)
+	for {
+		date := venusGreatestElongationInWindow(start, windowEnd)
+		if eventUTQueryBeforeOrEqual(date, jde) {
+			return date
+		}
+		prevSuperior := LastVenusSuperiorConjunction(eventUTLastQueryTT(start))
+		start, windowEnd = venusWestElongationWindowEndingAt(prevSuperior)
+	}
+}
+
+func venusGreatestElongation(jde float64) float64 {
+	east := venusSunRADelta(jde) > 0
+	if east {
+		return nextVenusGreatestElongationTyped(jde, true)
+	}
+	return nextVenusGreatestElongationTyped(jde, false)
 }
 
 func NextVenusGreatestElongation(jde float64) float64 {
-	date := venusGreatestElongation(jde)
-	if date < jde {
-		nextHe := NextVenusConjunction(jde)
-		return venusGreatestElongation(nextHe + 2)
+	east := NextVenusGreatestElongationEast(jde)
+	west := NextVenusGreatestElongationWest(jde)
+	if sameEventJD(east, west) {
+		return east
 	}
-	return date
+	if east < west {
+		return east
+	}
+	return west
 }
 
 func LastVenusGreatestElongation(jde float64) float64 {
-	lastHe := LastVenusConjunction(jde)
-	date := venusGreatestElongation(lastHe + 2)
-	if date > jde {
-		lastLastHe := LastVenusConjunction(lastHe - 2)
-		return venusGreatestElongation(lastLastHe + 2)
+	east := LastVenusGreatestElongationEast(jde)
+	west := LastVenusGreatestElongationWest(jde)
+	if sameEventJD(east, west) {
+		return east
 	}
-	return date
+	if east > west {
+		return east
+	}
+	return west
+}
+
+func LastVenusInferiorConjunctionInclusive(jde float64) float64 {
+	date := LastVenusConjunction(jde)
+	if venusConjunctionTypeAt(date) {
+		return date
+	}
+	return LastVenusConjunction(eventUTLastQueryTT(date))
+}
+
+func NextVenusInferiorConjunctionInclusive(jde float64) float64 {
+	date := NextVenusConjunction(jde)
+	if venusConjunctionTypeAt(date) {
+		return date
+	}
+	return NextVenusConjunction(eventUTNextQueryTT(date))
+}
+
+func LastVenusSuperiorConjunctionInclusive(jde float64) float64 {
+	date := LastVenusConjunction(jde)
+	if !venusConjunctionTypeAt(date) {
+		return date
+	}
+	return LastVenusConjunction(eventUTLastQueryTT(date))
+}
+
+func NextVenusSuperiorConjunctionInclusive(jde float64) float64 {
+	date := NextVenusConjunction(jde)
+	if !venusConjunctionTypeAt(date) {
+		return date
+	}
+	return NextVenusConjunction(eventUTNextQueryTT(date))
+}
+
+func LastVenusRetrogradeInclusive(jde float64) float64 {
+	return inclusiveLastSimpleEvent(jde, LastVenusRetrograde, NextVenusRetrograde)
+}
+
+func NextVenusRetrogradeInclusive(jde float64) float64 {
+	return inclusiveNextSimpleEvent(jde, LastVenusRetrograde, NextVenusRetrograde)
+}
+
+func LastVenusProgradeToRetrogradeInclusive(jde float64) float64 {
+	return inclusiveLastSimpleEvent(jde, LastVenusProgradeToRetrograde, NextVenusProgradeToRetrograde)
+}
+
+func NextVenusProgradeToRetrogradeInclusive(jde float64) float64 {
+	return inclusiveNextSimpleEvent(jde, LastVenusProgradeToRetrograde, NextVenusProgradeToRetrograde)
+}
+
+func LastVenusRetrogradeToProgradeInclusive(jde float64) float64 {
+	return inclusiveLastSimpleEvent(jde, LastVenusRetrogradeToPrograde, NextVenusRetrogradeToPrograde)
+}
+
+func NextVenusRetrogradeToProgradeInclusive(jde float64) float64 {
+	return inclusiveNextSimpleEvent(jde, LastVenusRetrogradeToPrograde, NextVenusRetrogradeToPrograde)
+}
+
+func LastVenusGreatestElongationInclusive(jde float64) float64 {
+	return inclusiveLastSimpleEvent(jde, LastVenusGreatestElongation, NextVenusGreatestElongation)
+}
+
+func NextVenusGreatestElongationInclusive(jde float64) float64 {
+	return inclusiveNextSimpleEvent(jde, LastVenusGreatestElongation, NextVenusGreatestElongation)
+}
+
+func LastVenusGreatestElongationEastInclusive(jde float64) float64 {
+	return inclusiveLastSimpleEvent(jde, LastVenusGreatestElongationEast, NextVenusGreatestElongationEast)
+}
+
+func NextVenusGreatestElongationEastInclusive(jde float64) float64 {
+	return inclusiveNextSimpleEvent(jde, LastVenusGreatestElongationEast, NextVenusGreatestElongationEast)
+}
+
+func LastVenusGreatestElongationWestInclusive(jde float64) float64 {
+	return inclusiveLastSimpleEvent(jde, LastVenusGreatestElongationWest, NextVenusGreatestElongationWest)
+}
+
+func NextVenusGreatestElongationWestInclusive(jde float64) float64 {
+	return inclusiveNextSimpleEvent(jde, LastVenusGreatestElongationWest, NextVenusGreatestElongationWest)
 }
 
 func NextVenusGreatestElongationEast(jde float64) float64 {
-	date := NextVenusGreatestElongation(jde)
-	sub := Limit360(VenusApparentRa(date) - SunApparentRa(date))
-	if sub > 180 {
-		return NextVenusGreatestElongation(date + 1)
-	}
-	return date
+	return nextVenusGreatestElongationTyped(jde, true)
 }
 
 func NextVenusGreatestElongationWest(jde float64) float64 {
-	date := NextVenusGreatestElongation(jde)
-	sub := Limit360(VenusApparentRa(date) - SunApparentRa(date))
-	if sub < 180 {
-		return NextVenusGreatestElongation(date + 1)
-	}
-	return date
+	return nextVenusGreatestElongationTyped(jde, false)
 }
 
 func LastVenusGreatestElongationEast(jde float64) float64 {
-	date := LastVenusGreatestElongation(jde)
-	sub := Limit360(VenusApparentRa(date) - SunApparentRa(date))
-	if sub > 180 {
-		return LastVenusGreatestElongation(date - 1)
-	}
-	return date
+	return lastVenusGreatestElongationTyped(jde, true)
 }
 
 func LastVenusGreatestElongationWest(jde float64) float64 {
-	date := LastVenusGreatestElongation(jde)
-	sub := Limit360(VenusApparentRa(date) - SunApparentRa(date))
-	if sub < 180 {
-		return LastVenusGreatestElongation(date - 1)
-	}
-	return date
+	return lastVenusGreatestElongationTyped(jde, false)
 }