Go: added function Search.

cosinekitty · Nov 6, 2023 · 7c34c1a · 7c34c1a
1 parent 0fa1343
commit 7c34c1a
Show file tree

Hide file tree

Showing 6 changed files with 274 additions and 4 deletions.
diff --git a/generate/template/astronomy.cs b/generate/template/astronomy.cs
@@ -5802,7 +5802,7 @@ public static AstroTime SearchSunLongitude(double targetLon, AstroTime startTime
         /// `Search` uses a combination of bisection and quadratic interpolation
         /// to minimize the number of function calls. However, it is critical that the
         /// supplied time window be small enough that there cannot be more than one root
-        /// (ascedning or descending) within it; otherwise the search can fail.
+        /// (ascending or descending) within it; otherwise the search can fail.
         /// Beyond that, it helps to make the time window as small as possible, ideally
         /// such that the function itself resembles a smooth parabolic curve within that window.
         ///

diff --git a/generate/template/astronomy.go b/generate/template/astronomy.go
@@ -2602,6 +2602,132 @@ type SearchContext interface {
 	Eval(time AstroTime) (float64, error)
 }
 
+// Searches for a time at which a function's value increases through zero.
+// Certain astronomy calculations involve finding a time when an event occurs.
+// Often such events can be defined as the root of a function:
+// the time at which the function's value becomes zero.
+//
+// Search finds the *ascending root* of a function: the time at which
+// the function's value becomes zero while having a positive slope. That is, as time increases,
+// the function transitions from a negative value, through zero at a specific moment,
+// to a positive value later. The goal of the search is to find that specific moment.
+//
+// Search uses a combination of bisection and quadratic interpolation
+// to minimize the number of function calls. However, it is critical that the
+// supplied time window be small enough that there cannot be more than one root
+// (ascending or descending) within it; otherwise the search can fail.
+// Beyond that, it helps to make the time window as small as possible, ideally
+// such that the function itself resembles a smooth parabolic curve within that window.
+//
+// If an ascending root is not found, or more than one root
+// (ascending and/or descending) exists within the window `t1`..`t2`,
+// the search will return `null`.
+//
+// If the search does not converge within 20 iterations, it will return an error.
+func Search(context SearchContext, t1, t2 AstroTime, dtToleranceSeconds float64) (*AstroTime, error) {
+	const iterLimit = 20
+	dtDays := math.Abs(dtToleranceSeconds / SecondsPerDay)
+	f1, e1 := context.Eval(t1)
+	if e1 != nil {
+		return nil, e1
+	}
+	f2, e2 := context.Eval(t2)
+	if e2 != nil {
+		return nil, e2
+	}
+	iter := 0
+	calcFmid := true
+	fmid := 0.0
+	for {
+		iter++
+		if iter > iterLimit {
+			return nil, errors.New("Search did not converge within 20 iterations")
+		}
+
+		dt := (t2.Tt - t1.Tt) / 2.0
+		tmid := t1.AddDays(dt)
+		if math.Abs(dt) < dtDays {
+			// We are close enough to the event to stop the search.
+			return &tmid, nil
+		}
+
+		if calcFmid {
+			if f3, e3 := context.Eval(tmid); e3 != nil {
+				return nil, e3
+			} else {
+				fmid = f3
+			}
+		} else {
+			calcFmid = true // we already have the correct value of fmid from the previous loop
+		}
+
+		// Quadratic interpolation:
+		// Try to find a parabola that passes through the 3 points we have sampled:
+		// (t1,f1), (tmid,fmid), (t2,f2)
+
+		if success, qUt, qDfDt := quadInterp(tmid.Ut, t2.Ut-tmid.Ut, f1, fmid, f2); success {
+			tq := TimeFromUniversalDays(qUt)
+			fq, eq := context.Eval(tq)
+			if eq != nil {
+				return nil, eq
+			}
+			if qDfDt != 0.0 {
+				dtGuess := math.Abs(fq / qDfDt)
+				if dtGuess < dtDays {
+					// The estimated time error is small enough that we can quit now.
+					return &tq, nil
+				}
+
+				// Try guessing a tighter boundary with the interpolated root at the center.
+				dtGuess *= 1.2
+				if dtGuess < dt/10.0 {
+					tleft := tq.AddDays(-dtGuess)
+					tright := tq.AddDays(+dtGuess)
+					if (tleft.Ut-t1.Ut)*(tleft.Ut-t2.Ut) < 0.0 {
+						if (tright.Ut-t1.Ut)*(tright.Ut-t2.Ut) < 0.0 {
+							fleft, eleft := context.Eval(tleft)
+							if eleft != nil {
+								return nil, eleft
+							}
+							fright, eright := context.Eval(tright)
+							if eright != nil {
+								return nil, eright
+							}
+							if fleft < 0.0 && fright >= 0.0 {
+								f1 = fleft
+								f2 = fright
+								t1 = tleft
+								t2 = tright
+								fmid = fq
+								calcFmid = false // save a little work -- no need to re-calculate fmid next time around the loop
+								continue
+							}
+						}
+					}
+				}
+			}
+		}
+
+		// After quadratic interpolation attempt.
+		// Now just divide the region in two parts and pick whichever one appears to contain a root.
+		if f1 < 0.0 && fmid >= 0.0 {
+			t2 = tmid
+			f2 = fmid
+			continue
+		}
+
+		if fmid < 0.0 && f2 >= 0.0 {
+			t1 = tmid
+			f1 = fmid
+			continue
+		}
+
+		// Either there is no ascending zero-crossing in this range
+		// or the search window is too wide (more than one zero-crossing).
+		return nil, nil
+	}
+}
+
 func findAscent(depth int, context SearchContext, maxDerivAlt float64, t1, t2 AstroTime, a1, a2 float64) ascentInfo {
 	// See if we can find any time interval where the altitude-diff function
 	// rises from non-positive to positive.

diff --git a/source/csharp/README.md b/source/csharp/README.md
@@ -1842,7 +1842,7 @@ reports a solution outside this time window.
 `Search` uses a combination of bisection and quadratic interpolation
 to minimize the number of function calls. However, it is critical that the
 supplied time window be small enough that there cannot be more than one root
-(ascedning or descending) within it; otherwise the search can fail.
+(ascending or descending) within it; otherwise the search can fail.
 Beyond that, it helps to make the time window as small as possible, ideally
 such that the function itself resembles a smooth parabolic curve within that window.
 

diff --git a/source/csharp/astronomy.cs b/source/csharp/astronomy.cs
@@ -6936,7 +6936,7 @@ public static AstroTime SearchSunLongitude(double targetLon, AstroTime startTime
         /// `Search` uses a combination of bisection and quadratic interpolation
         /// to minimize the number of function calls. However, it is critical that the
         /// supplied time window be small enough that there cannot be more than one root
-        /// (ascedning or descending) within it; otherwise the search can fail.
+        /// (ascending or descending) within it; otherwise the search can fail.
         /// Beyond that, it helps to make the time window as small as possible, ideally
         /// such that the function itself resembles a smooth parabolic curve within that window.
         ///

diff --git a/source/golang/README.md b/source/golang/README.md
@@ -28,6 +28,7 @@ It provides a suite of well\-tested functions for calculating positions of the S
 - [type AstroMoonQuarter](<#AstroMoonQuarter>)
 - [type AstroSearchFunc](<#AstroSearchFunc>)
 - [type AstroTime](<#AstroTime>)
+  - [func Search\(context SearchContext, t1, t2 AstroTime, dtToleranceSeconds float64\) \(\*AstroTime, error\)](<#Search>)
   - [func TimeFromCalendar\(year, month, day, hour, minute int, second float64\) AstroTime](<#TimeFromCalendar>)
   - [func TimeFromTerrestrialDays\(tt float64\) AstroTime](<#TimeFromTerrestrialDays>)
   - [func TimeFromUniversalDays\(ut float64\) AstroTime](<#TimeFromUniversalDays>)
@@ -356,6 +357,23 @@ type AstroTime struct {
 }
 ```
 
+<a name="Search"></a>
+### func [Search](<https://github.com/cosinekitty/astronomy/blob/golang/source/golang/astronomy.go#L2627>)
+
+```go
+func Search(context SearchContext, t1, t2 AstroTime, dtToleranceSeconds float64) (*AstroTime, error)
+```
+
+Searches for a time at which a function's value increases through zero. Certain astronomy calculations involve finding a time when an event occurs. Often such events can be defined as the root of a function: the time at which the function's value becomes zero.
+
+Search finds the \*ascending root\* of a function: the time at which the function's value becomes zero while having a positive slope. That is, as time increases, the function transitions from a negative value, through zero at a specific moment, to a positive value later. The goal of the search is to find that specific moment.
+
+Search uses a combination of bisection and quadratic interpolation to minimize the number of function calls. However, it is critical that the supplied time window be small enough that there cannot be more than one root \(ascending or descending\) within it; otherwise the search can fail. Beyond that, it helps to make the time window as small as possible, ideally such that the function itself resembles a smooth parabolic curve within that window.
+
+If an ascending root is not found, or more than one root \(ascending and/or descending\) exists within the window \`t1\`..\`t2\`, the search will return \`null\`.
+
+If the search does not converge within 20 iterations, it will return an error.
+
 <a name="TimeFromCalendar"></a>
 ### func [TimeFromCalendar](<https://github.com/cosinekitty/astronomy/blob/golang/source/golang/astronomy.go#L384>)
 
@@ -896,7 +914,7 @@ type StateVector struct {
 ```
 
 <a name="LagrangePointFast"></a>
-### func [LagrangePointFast](<https://github.com/cosinekitty/astronomy/blob/golang/source/golang/astronomy.go#L2905>)
+### func [LagrangePointFast](<https://github.com/cosinekitty/astronomy/blob/golang/source/golang/astronomy.go#L3031>)
 
 ```go
 func LagrangePointFast(point int, majorState StateVector, majorMass float64, minorState StateVector, minorMass float64) (*StateVector, error)

diff --git a/source/golang/astronomy.go b/source/golang/astronomy.go
@@ -2602,6 +2602,132 @@ type SearchContext interface {
 	Eval(time AstroTime) (float64, error)
 }
 
+// Searches for a time at which a function's value increases through zero.
+// Certain astronomy calculations involve finding a time when an event occurs.
+// Often such events can be defined as the root of a function:
+// the time at which the function's value becomes zero.
+//
+// Search finds the *ascending root* of a function: the time at which
+// the function's value becomes zero while having a positive slope. That is, as time increases,
+// the function transitions from a negative value, through zero at a specific moment,
+// to a positive value later. The goal of the search is to find that specific moment.
+//
+// Search uses a combination of bisection and quadratic interpolation
+// to minimize the number of function calls. However, it is critical that the
+// supplied time window be small enough that there cannot be more than one root
+// (ascending or descending) within it; otherwise the search can fail.
+// Beyond that, it helps to make the time window as small as possible, ideally
+// such that the function itself resembles a smooth parabolic curve within that window.
+//
+// If an ascending root is not found, or more than one root
+// (ascending and/or descending) exists within the window `t1`..`t2`,
+// the search will return `null`.
+//
+// If the search does not converge within 20 iterations, it will return an error.
+func Search(context SearchContext, t1, t2 AstroTime, dtToleranceSeconds float64) (*AstroTime, error) {
+	const iterLimit = 20
+	dtDays := math.Abs(dtToleranceSeconds / SecondsPerDay)
+	f1, e1 := context.Eval(t1)
+	if e1 != nil {
+		return nil, e1
+	}
+	f2, e2 := context.Eval(t2)
+	if e2 != nil {
+		return nil, e2
+	}
+	iter := 0
+	calcFmid := true
+	fmid := 0.0
+	for {
+		iter++
+		if iter > iterLimit {
+			return nil, errors.New("Search did not converge within 20 iterations")
+		}
+
+		dt := (t2.Tt - t1.Tt) / 2.0
+		tmid := t1.AddDays(dt)
+		if math.Abs(dt) < dtDays {
+			// We are close enough to the event to stop the search.
+			return &tmid, nil
+		}
+
+		if calcFmid {
+			if f3, e3 := context.Eval(tmid); e3 != nil {
+				return nil, e3
+			} else {
+				fmid = f3
+			}
+		} else {
+			calcFmid = true // we already have the correct value of fmid from the previous loop
+		}
+
+		// Quadratic interpolation:
+		// Try to find a parabola that passes through the 3 points we have sampled:
+		// (t1,f1), (tmid,fmid), (t2,f2)
+
+		if success, qUt, qDfDt := quadInterp(tmid.Ut, t2.Ut-tmid.Ut, f1, fmid, f2); success {
+			tq := TimeFromUniversalDays(qUt)
+			fq, eq := context.Eval(tq)
+			if eq != nil {
+				return nil, eq
+			}
+			if qDfDt != 0.0 {
+				dtGuess := math.Abs(fq / qDfDt)
+				if dtGuess < dtDays {
+					// The estimated time error is small enough that we can quit now.
+					return &tq, nil
+				}
+
+				// Try guessing a tighter boundary with the interpolated root at the center.
+				dtGuess *= 1.2
+				if dtGuess < dt/10.0 {
+					tleft := tq.AddDays(-dtGuess)
+					tright := tq.AddDays(+dtGuess)
+					if (tleft.Ut-t1.Ut)*(tleft.Ut-t2.Ut) < 0.0 {
+						if (tright.Ut-t1.Ut)*(tright.Ut-t2.Ut) < 0.0 {
+							fleft, eleft := context.Eval(tleft)
+							if eleft != nil {
+								return nil, eleft
+							}
+							fright, eright := context.Eval(tright)
+							if eright != nil {
+								return nil, eright
+							}
+							if fleft < 0.0 && fright >= 0.0 {
+								f1 = fleft
+								f2 = fright
+								t1 = tleft
+								t2 = tright
+								fmid = fq
+								calcFmid = false // save a little work -- no need to re-calculate fmid next time around the loop
+								continue
+							}
+						}
+					}
+				}
+			}
+		}
+
+		// After quadratic interpolation attempt.
+		// Now just divide the region in two parts and pick whichever one appears to contain a root.
+		if f1 < 0.0 && fmid >= 0.0 {
+			t2 = tmid
+			f2 = fmid
+			continue
+		}
+
+		if fmid < 0.0 && f2 >= 0.0 {
+			t1 = tmid
+			f1 = fmid
+			continue
+		}
+
+		// Either there is no ascending zero-crossing in this range
+		// or the search window is too wide (more than one zero-crossing).
+		return nil, nil
+	}
+}
+
 func findAscent(depth int, context SearchContext, maxDerivAlt float64, t1, t2 AstroTime, a1, a2 float64) ascentInfo {
 	// See if we can find any time interval where the altitude-diff function
 	// rises from non-positive to positive.