From 645bf67427ca73e2a2e0ca2d5df33a4c15a1351a Mon Sep 17 00:00:00 2001 From: Don Cross Date: Mon, 6 Nov 2023 21:01:42 -0500 Subject: [PATCH] Go: moved Search function to a better location. --- generate/template/astronomy.go | 252 ++++++++++++++++----------------- source/golang/README.md | 2 +- source/golang/astronomy.go | 252 ++++++++++++++++----------------- 3 files changed, 253 insertions(+), 253 deletions(-) diff --git a/generate/template/astronomy.go b/generate/template/astronomy.go index 50f268c7..1ca00da9 100644 --- a/generate/template/astronomy.go +++ b/generate/template/astronomy.go @@ -2602,132 +2602,6 @@ 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. @@ -3002,6 +2876,132 @@ func quadInterp(tm, dt, fa, fm, fb float64) (bool, float64, float64) { return true, outT, outDfDt } +// 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 + } +} + //--- Search code ends here --------------------------------------------------------------------- // Calculates one of the 5 Lagrange points from body masses and state vectors. diff --git a/source/golang/README.md b/source/golang/README.md index 141c807a..3a025177 100644 --- a/source/golang/README.md +++ b/source/golang/README.md @@ -358,7 +358,7 @@ type AstroTime struct { ``` -### func [Search]() +### func [Search]() ```go func Search(context SearchContext, t1, t2 AstroTime, dtToleranceSeconds float64) (*AstroTime, error) diff --git a/source/golang/astronomy.go b/source/golang/astronomy.go index 5bea07bb..aa9fe535 100644 --- a/source/golang/astronomy.go +++ b/source/golang/astronomy.go @@ -2602,132 +2602,6 @@ 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. @@ -3002,6 +2876,132 @@ func quadInterp(tm, dt, fa, fm, fb float64) (bool, float64, float64) { return true, outT, outDfDt } +// 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 + } +} + //--- Search code ends here --------------------------------------------------------------------- // Calculates one of the 5 Lagrange points from body masses and state vectors.