From 221ea1130a43ee1e30fa56aea04efd157bf4eccd Mon Sep 17 00:00:00 2001 From: Don Cross Date: Wed, 26 Jun 2019 17:36:29 -0400 Subject: [PATCH] Implemented Python Search function. Slight tweaks to C and JS versions. --- generate/template/astronomy.c | 5 ++- generate/template/astronomy.js | 2 +- generate/template/astronomy.py | 72 ++++++++++++++++++++++++++++++++++ source/c/astronomy.c | 5 ++- source/js/astronomy.js | 2 +- source/js/astronomy.min.js | 2 +- source/python/astronomy.py | 72 ++++++++++++++++++++++++++++++++++ 7 files changed, 153 insertions(+), 7 deletions(-) diff --git a/generate/template/astronomy.c b/generate/template/astronomy.c index 8d381c55..90545f61 100644 --- a/generate/template/astronomy.c +++ b/generate/template/astronomy.c @@ -2425,7 +2425,8 @@ astro_search_result_t Astronomy_Search( CALLFUNC(fq, tq); if (q_df_dt != 0.0) { - if (fabs(fq / q_df_dt) < dt_days) + dt_guess = fabs(fq / q_df_dt); + if (dt_guess < dt_days) { /* The estimated time error is small enough that we can quit now. */ result.time = tq; @@ -2434,7 +2435,7 @@ astro_search_result_t Astronomy_Search( } /* Try guessing a tighter boundary with the interpolated root at the center. */ - dt_guess = 1.2 * fabs(fq / q_df_dt); + dt_guess *= 1.2; if (dt_guess < dt/10.0) { astro_time_t tleft = Astronomy_AddDays(tq, -dt_guess); diff --git a/generate/template/astronomy.js b/generate/template/astronomy.js index 5b6eb95b..8d04963d 100644 --- a/generate/template/astronomy.js +++ b/generate/template/astronomy.js @@ -1942,7 +1942,7 @@ Astronomy.Search = function(func, t1, t2, options) { return func(t); } - const dt_days = dt_tolerance_seconds / SECONDS_PER_DAY; + const dt_days = Math.abs(dt_tolerance_seconds / SECONDS_PER_DAY); ++Perf.search; diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index de59bdf7..405060ff 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -1007,6 +1007,78 @@ def _QuadInterp(tm, dt, fa, fm, fb): df_dt = (2*Q*x + R) / dt return (x, t, df_dt) +def Search(func, t1, t2, dt_tolerance_seconds): + dt_days = abs(dt_tolerance_seconds / _SECONDS_PER_DAY) + f1 = func(t1) + f2 = func(t2) + iter = 0 + iter_limit = 20 + calc_fmid = True + while True: + iter += 1 + if iter > iter_limit: + raise Error('Excessive iteration in Search') + + dt = (t2.tt - t1.tt) / 2.0 + tmid = t1.AddDays(dt) + if abs(dt) < dt_days: + # We are close enough to the event to stop the search. + return tmid + + if calc_fmid: + fmid = func(tmid) + else: + # We already have the correct value of fmid from the previous loop. + calc_fmid = True + + # Quadratic interpolation: + # Try to find a parabola that passes through the 3 points we have sampled: + # (t1,f1), (tmid,fmid), (t2,f2). + q = _QuadInterp(tmid.ut, t2.ut - tmid.ut, f1, fmid, f2) + if q: + (q_x, q_ut, q_df_dt) = q + tq = Time(q_ut) + fq = func(tq) + if q_df_dt != 0.0: + dt_guess = abs(fq / q_df_dt) + if dt_guess < dt_days: + # The estimated time error is small enough that we can quit now. + return tq + + # Try guessing a tighter boundary with the interpolated root at the center. + dt_guess *= 1.2 + if dt_guess < dt / 10.0: + tleft = tq.AddDays(-dt_guess) + tright = tq.AddDays(+dt_guess) + if (tleft.ut - t1.ut)*(tleft.ut - t2.ut) < 0.0: + if (tright.ut - t1.ut)*(tright.ut - t2.ut) < 0.0: + fleft = func(tleft) + fright = func(tright) + if fleft < 0.0 and fright >= 0.0: + f1 = fleft + f2 = fright + t1 = tleft + t2 = tright + fmid = fq + calc_fmid = False + continue + + # Quadratic interpolation attempt did not work out. + # Just divide the region in two parts and pick whichever one appears to contain a root. + if f1 < 0.0 and fmid >= 0.0: + t2 = tmid + f2 = fmid + continue + + if fmid < 0.0 and 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 None + # END Search #---------------------------------------------------------------------------- diff --git a/source/c/astronomy.c b/source/c/astronomy.c index 1f9773b8..3e81ac81 100644 --- a/source/c/astronomy.c +++ b/source/c/astronomy.c @@ -3481,7 +3481,8 @@ astro_search_result_t Astronomy_Search( CALLFUNC(fq, tq); if (q_df_dt != 0.0) { - if (fabs(fq / q_df_dt) < dt_days) + dt_guess = fabs(fq / q_df_dt); + if (dt_guess < dt_days) { /* The estimated time error is small enough that we can quit now. */ result.time = tq; @@ -3490,7 +3491,7 @@ astro_search_result_t Astronomy_Search( } /* Try guessing a tighter boundary with the interpolated root at the center. */ - dt_guess = 1.2 * fabs(fq / q_df_dt); + dt_guess *= 1.2; if (dt_guess < dt/10.0) { astro_time_t tleft = Astronomy_AddDays(tq, -dt_guess); diff --git a/source/js/astronomy.js b/source/js/astronomy.js index caa17119..73a5c35b 100644 --- a/source/js/astronomy.js +++ b/source/js/astronomy.js @@ -2780,7 +2780,7 @@ Astronomy.Search = function(func, t1, t2, options) { return func(t); } - const dt_days = dt_tolerance_seconds / SECONDS_PER_DAY; + const dt_days = Math.abs(dt_tolerance_seconds / SECONDS_PER_DAY); ++Perf.search; diff --git a/source/js/astronomy.min.js b/source/js/astronomy.min.js index ff6258ec..15e40584 100644 --- a/source/js/astronomy.min.js +++ b/source/js/astronomy.min.js @@ -121,7 +121,7 @@ e.Equator=function(a,b,d,c,g){b=e.MakeTime(b);var h=V(b),k=.017453292519943295*d function(a,b,d){void 0===T&&(T=.017453292519943295*P(e.MakeTime(Y)).mobl,ba=Math.cos(T),ca=Math.sin(T));return aa(a,b,d,ba,ca)};e.GeoMoon=function(a){a=e.MakeTime(a);var b=E(a),d=b.distance_au*Math.cos(b.geo_eclip_lat);b=[d*Math.cos(b.geo_eclip_lon),d*Math.sin(b.geo_eclip_lon),b.distance_au*Math.sin(b.geo_eclip_lat)];var c=.017453292519943295*B(a);d=Math.cos(c);c=Math.sin(c);b=R(a.tt,[b[0],b[1]*d-b[2]*c,b[1]*c+b[2]*d],0);return new y(b[0],b[1],b[2],a)};e.HelioVector=function(a,b){b=e.MakeTime(b); if(a in z)return F(z[a],b);if(a in da){a:{var d=$jscomp.makeIterator(da[a]);for(a=d.next();!a.done;a=d.next()){a=a.value;var c=a.tt;var g=a.tt+a.ndays;c=(2*b.tt-(g+c))/(g-c);if(-1<=c&&1>=c){var h=[];for(g=0;3>g;++g){var k=1;var l=a.coeff[0][g];var m=c;l+=a.coeff[1][g]*m;for(d=2;dl;++l){g=e.HelioVector(a,k);d&&(c=F(z.Earth,k));h=new y(g.x-c.x,g.y-c.y,g.z-c.z,b);if("Sun"===a)return h;var m=b.AddDays(-h.Length()/173.1446326846693);g=Math.abs(m.tt-k.tt);if(1E-9>g)return h; -k=m}throw"Light-travel time solver did not converge: dt="+g;};e.Search=function(a,b,d,c){function g(b){++I.search_func;return a(b)}var h=(c&&c.dt_tolerance_seconds||1)/86400;++I.search;var k=c&&c.init_f1||g(b),l=c&&c.init_f2||g(d),m,f=0;c=c&&c.iter_limit||20;for(var n=!0;;){if(++f>c)throw"Excessive iteration in Search()";var p=new C(b.ut+.5*(d.ut-b.ut)),v=p.ut-b.ut;if(Math.abs(v)(v.ut-b.ut)*(v.ut-d.ut)&&0>(w.ut-b.ut)*(w.ut-d.ut))){q=g(v);var t=g(w);if(0>q&&0<=t){k=q;l=t;b=v;d=w;m=r;n=!1;continue}}}}if(0>k&&0<=m)d=p,l=m;else if(0>m&&0<=l)b=p,k=m;else return null}};e.SearchSunLongitude=function(a,b,d){b=e.MakeTime(b);d=b.AddDays(d);return e.Search(function(b){b=e.SunPosition(b);return J(b.elon-a)},b,d)};e.LongitudeFromSun=function(a,b){if("Earth"===a)throw"The Earth does not have a longitude as seen from itself."; b=e.MakeTime(b);a=e.GeoVector(a,b,!0);a=e.Ecliptic(a.x,a.y,a.z);b=e.GeoVector("Sun",b,!0);b=e.Ecliptic(b.x,b.y,b.z);for(b=a.elon-b.elon;0>b;)b+=360;for(;360<=b;)b-=360;return b};e.AngleFromSun=function(a,b){if("Earth"==a)throw"The Earth does not have an angle as seen from itself.";var d=e.GeoVector("Sun",b,!0);a=e.GeoVector(a,b,!0);return u(d,a)};e.EclipticLongitude=function(a,b){if("Sun"===a)throw"Cannot calculate heliocentric longitude of the Sun.";a=e.HelioVector(a,b);return e.Ecliptic(a.x,a.y, a.z).elon};var na=function(a,b,d,c,e,h,k,l){this.time=a;this.mag=b;this.phase_angle=d;this.phase_fraction=(1+Math.cos(.017453292519943295*d))/2;this.helio_dist=c;this.geo_dist=e;this.gc=h;this.hc=k;this.ring_tilt=l};e.Illumination=function(a,b){if("Earth"===a)throw"The illumination of the Earth is not defined.";var d=e.MakeTime(b),c=F(z.Earth,d);if("Sun"===a){var g=new y(-c.x,-c.y,-c.z,d);b=new y(0,0,0,d);c=0}else"Moon"===a?(g=e.GeoMoon(d),b=new y(c.x+g.x,c.y+g.y,c.z+g.z,d)):(b=e.HelioVector(a,b), diff --git a/source/python/astronomy.py b/source/python/astronomy.py index c13e830c..1cf7625d 100644 --- a/source/python/astronomy.py +++ b/source/python/astronomy.py @@ -1845,6 +1845,78 @@ def _QuadInterp(tm, dt, fa, fm, fb): df_dt = (2*Q*x + R) / dt return (x, t, df_dt) +def Search(func, t1, t2, dt_tolerance_seconds): + dt_days = abs(dt_tolerance_seconds / _SECONDS_PER_DAY) + f1 = func(t1) + f2 = func(t2) + iter = 0 + iter_limit = 20 + calc_fmid = True + while True: + iter += 1 + if iter > iter_limit: + raise Error('Excessive iteration in Search') + + dt = (t2.tt - t1.tt) / 2.0 + tmid = t1.AddDays(dt) + if abs(dt) < dt_days: + # We are close enough to the event to stop the search. + return tmid + + if calc_fmid: + fmid = func(tmid) + else: + # We already have the correct value of fmid from the previous loop. + calc_fmid = True + + # Quadratic interpolation: + # Try to find a parabola that passes through the 3 points we have sampled: + # (t1,f1), (tmid,fmid), (t2,f2). + q = _QuadInterp(tmid.ut, t2.ut - tmid.ut, f1, fmid, f2) + if q: + (q_x, q_ut, q_df_dt) = q + tq = Time(q_ut) + fq = func(tq) + if q_df_dt != 0.0: + dt_guess = abs(fq / q_df_dt) + if dt_guess < dt_days: + # The estimated time error is small enough that we can quit now. + return tq + + # Try guessing a tighter boundary with the interpolated root at the center. + dt_guess *= 1.2 + if dt_guess < dt / 10.0: + tleft = tq.AddDays(-dt_guess) + tright = tq.AddDays(+dt_guess) + if (tleft.ut - t1.ut)*(tleft.ut - t2.ut) < 0.0: + if (tright.ut - t1.ut)*(tright.ut - t2.ut) < 0.0: + fleft = func(tleft) + fright = func(tright) + if fleft < 0.0 and fright >= 0.0: + f1 = fleft + f2 = fright + t1 = tleft + t2 = tright + fmid = fq + calc_fmid = False + continue + + # Quadratic interpolation attempt did not work out. + # Just divide the region in two parts and pick whichever one appears to contain a root. + if f1 < 0.0 and fmid >= 0.0: + t2 = tmid + f2 = fmid + continue + + if fmid < 0.0 and 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 None + # END Search #----------------------------------------------------------------------------