.

utils/interpolate.cpp

@@ -302,9 +302,6 @@ void XInterpolator::create(const SInitOrbit& orbit, const TimeJD& beg, const Tim
 	m_ephi = new ei;
 	m_beg = m_positions.front().get_date();
 	m_end = m_positions.back().get_date();
-
-//	std::string sdate = DateToStr(m_beg);
-//	std::string edate = DateToStr(m_end);
 	
 	init(*m_ephi, m_positions.size());
 

utils/orbBlocks.cpp

 #include "orbBlocks.h"
 #include "TimeFunctions.h"
 #include "Propagator.h"
+#include "xalgorithm.h"
 
 template<bool dir, typename vartype> struct cmporb : public std::binary_function<vartype,vartype,bool> 
 {
@@ -24,7 +25,7 @@ OrbBlocksStore::~OrbBlocksStore()
 void OrbBlocksStore::init(const SInitOrbits& orbits, const TimeJD& beg, const TimeJD& end)
 {
 	init_min_dist(filter(orbits), beg, end);
-
+	
 	for (OrbBlocks::const_iterator it = m_blocks.begin(); it != m_blocks.end(); it++){
 		m_interpolators.push_back(new XInterpolator(it->get_orbit(), it->get_beg(), it->get_end()));
 	}
@@ -163,47 +164,70 @@ void OrbBlocksStore::init_min_dist(const SInitOrbits& orbits, const TimeJD& beg,
 	}
 }
 
-void OrbBlocksStore::calc_min_dist(const SInitOrbit& orb1, const SInitOrbit& orb2, const TimeJD& beg, const TimeJD& end, TimeJD& result_jd, double& result_dist)
+class MinDistEvaluator
 {
-//	Interpolator ip(orb1, beg, end);
-//	Interpolator ip(orb2, beg, end);
+public:
+	MinDistEvaluator(const XInterpolator& ip1, const XInterpolator& ip2, const TimeJD& base)
+		:m_ip1(ip1), m_ip2(ip2), m_base(base)
+	{
+	}
+	double calc(double time) const {
+		TimeJD date = ShiftDate(m_base, time*86400);
 
-	Propagator prop1(orb1, false);
-	Propagator prop2(orb2, false);
+		Vect6 p1;
+		m_ip1.get_pos(date, p1);
+
+		Vect6 p2;
+		m_ip2.get_pos(date, p2);
+
+		return sqrt(
+			sqr(p1.x - p2.x) + 
+			sqr(p1.y - p2.y) + 
+			sqr(p1.z - p2.z)) / 1000.0;
+	}
+private:
+	TimeJD m_base;
+	const XInterpolator& m_ip1;
+	const XInterpolator& m_ip2;
+};
 
-	TimeJD date = beg;
+void OrbBlocksStore::calc_min_dist(const SInitOrbit& orb1, const SInitOrbit& orb2, const TimeJD& beg, const TimeJD& end, TimeJD& result_jd, double& result_dist)
+{
+	XInterpolator ip1(orb1, beg, end);
+	XInterpolator ip2(orb2, beg, end);
+	MinDistEvaluator eval(ip1, ip2, beg);
+
+	double step = (int)(std::min(orb1.GetKeplerData().GetPeriod(), orb2.GetKeplerData().GetPeriod()) * 1000 / 20);
+	double time = 0;
+	double time_end = DaysInterval(beg, end);
 	double min = -1;
-	TimeJD minjd;
+	double time_min;
 
 	bool exit = false;
 	while(!exit){
 
-		//std::string date_s = DateToStr(date); 
-
-		if(date >= end){
-			date = end;
+		if(time >= time_end){
+			time = time_end;
 			exit = true;
 		}
-
-		PhasePoint6D vec1, vec2;
-		prop1.Propagate(date, vec1);
-		prop2.Propagate(date, vec2);
-
-		double dist = sqrt(
-			sqr(vec1.CoordsVel.x - vec2.CoordsVel.x) + 
-			sqr(vec1.CoordsVel.y - vec2.CoordsVel.y) + 
-			sqr(vec1.CoordsVel.z - vec2.CoordsVel.z));
-
+		
+		double dist = eval.calc(time);
+		
 		if(min == -1 || min > dist){
 			min = dist;
-			minjd = date;
+			time_min = time;
 		}
 
-		date = ShiftDate(date, 600);
+		time += step / 86400.0;
 	}
 
-	result_jd = minjd;
-	result_dist = min;
+	double tol = 60 / 86400.0;
+	double left = std::max(0.0, time_min - step / 86400.0 /2.0);
+	double right = std::min(time_end, time_min + step / 86400.0 /2.0);
+	double shift = golden_section_search(eval, left, right, tol);
+	
+	result_jd = ShiftDate(beg, shift * 86400);
+	result_dist = eval.calc(shift);
 }
 
 bool OrbBlocksStore::get_pos(const TimeJD& date, Vect6& result) const

utils/xalgorithm.cpp

+#include "xalgorithm.h"
+

utils/xalgorithm.h

+#pragma once
+
+template<typename T> double golden_section_search(const T& f, double a, double b, double tol)
+{
+	double invphi = (sqrt(5.0) - 1) / 2.0; 
+	double invphi2 = (3 - sqrt(5.0)) / 2.0; 
+
+	a = std::min(a, b);
+	b = std::max(a, b);
+
+	double h = b - a;
+
+	if (h <= tol){
+		return (a + b) / 2.0;
+	}
+
+	int n = ceil( log(tol/h) / log(invphi) );
+
+	double c = a + invphi2 * h;
+	double d = a + invphi * h;
+	double yc = f.calc(c);
+	double yd = f.calc(d);
+
+	for(int k=0; k < n - 1; k++){
+		if (yc < yd){
+			b = d;
+			d = c;
+			yd = yc;
+			h = invphi * h;
+			c = a + invphi2 * h;
+			yc = f.calc(c);
+		}else{
+			a = c;
+			c = d;
+			yc = yd;
+			h = invphi * h;
+			d = a + invphi * h;
+			yd = f.calc(d);
+		}
+	}
+
+	if (yc < yd){
+		return (a + d) / 2.0;
+	}else{
+		return (c + b) / 2.0;
+	}
+}
+