Commit 8ca60907 authored by Alexander Lapshin's avatar Alexander Lapshin

bom fix and interpolator

parent cfdb31bd
...@@ -319,12 +319,16 @@ void XInterpolator::create(const SInitOrbit& orbit, const TimeJD& beg, const Tim ...@@ -319,12 +319,16 @@ void XInterpolator::create(const SInitOrbit& orbit, const TimeJD& beg, const Tim
finalize(*m_ephi); finalize(*m_ephi);
} }
bool XInterpolator::get_pos(const TimeJD& date, Vect6& res) const Vect6 XInterpolator::get_pos(const TimeJD& date) const
{ {
assert(m_ephi != 0); if(!m_ephi){
assert(date >= m_beg); throw Exp() << "interpolation not initialised\n";
assert(date <= m_end); }
if(!(date >= m_beg && date <= m_end)){
throw Exp() << "interpolation date out of range\n";
}
double outvec[6], backvec[6], aberout[3], aberback[3]; double outvec[6], backvec[6], aberout[3], aberback[3];
double relcorrout, relcorrback; double relcorrout, relcorrback;
int fvel = 1; int fvel = 1;
...@@ -333,7 +337,11 @@ bool XInterpolator::get_pos(const TimeJD& date, Vect6& res) const ...@@ -333,7 +337,11 @@ bool XInterpolator::get_pos(const TimeJD& date, Vect6& res) const
double sec = date.GetFraction() * 86400; double sec = date.GetFraction() * 86400;
cinterp_pvac(*m_ephi, jdi, sec, outvec, backvec, aberout, aberback, &relcorrout, &relcorrback, fvel, &ircode); cinterp_pvac(*m_ephi, jdi, sec, outvec, backvec, aberout, aberback, &relcorrout, &relcorrback, fvel, &ircode);
res = Vect6( if(ircode != 0){
throw Exp() << "interpolation fail\n";
}
return Vect6(
outvec[0], outvec[0],
outvec[1], outvec[1],
outvec[2], outvec[2],
...@@ -341,6 +349,4 @@ bool XInterpolator::get_pos(const TimeJD& date, Vect6& res) const ...@@ -341,6 +349,4 @@ bool XInterpolator::get_pos(const TimeJD& date, Vect6& res) const
outvec[4], outvec[4],
outvec[5] outvec[5]
); );
return true;
} }
...@@ -13,7 +13,7 @@ public: ...@@ -13,7 +13,7 @@ public:
XInterpolator(const SInitOrbit& orbit, const TimeJD& beg, const TimeJD& end); XInterpolator(const SInitOrbit& orbit, const TimeJD& beg, const TimeJD& end);
~XInterpolator(); ~XInterpolator();
void create(const SInitOrbit& orbit, const TimeJD& beg, const TimeJD& end); void create(const SInitOrbit& orbit, const TimeJD& beg, const TimeJD& end);
bool get_pos(const TimeJD& date, Vect6& res) const; Vect6 get_pos(const TimeJD& date) const;
int get_size() const { return m_positions.size(); } int get_size() const { return m_positions.size(); }
TimeJD get_beg() const { return m_beg; } TimeJD get_beg() const { return m_beg; }
TimeJD get_end() const { return m_end; } TimeJD get_end() const { return m_end; }
......
#include "orbBlocks.h" #include "orbBlocks.h"
#include "TimeFunctions.h" #include "TimeFunctions.h"
#include "Propagator.h" #include "Propagator.h"
#include "xalgorithm.h" #include "xalgorithm.h"
...@@ -169,18 +169,17 @@ void OrbBlocksStore::init_min_dist(const SInitOrbits& orbits, const TimeJD& beg, ...@@ -169,18 +169,17 @@ void OrbBlocksStore::init_min_dist(const SInitOrbits& orbits, const TimeJD& beg,
class MinDistEvaluator class MinDistEvaluator
{ {
public: public:
MinDistEvaluator(const XInterpolator& ip1, const XInterpolator& ip2, const TimeJD& base) MinDistEvaluator(const XInterpolator& ip1, const XInterpolator& ip2, const TimeJD& base):
:m_ip1(ip1), m_ip2(ip2), m_base(base) m_ip1(ip1),
m_ip2(ip2),
m_base(base)
{ {
} }
double calc(double time) const { double calc(double time) const {
TimeJD date = ShiftDate(m_base, time * 86400); TimeJD date = ShiftDate(m_base, time * 86400);
Vect6 p1; Vect6 p1 = m_ip1.get_pos(date);
m_ip1.get_pos(date, p1); Vect6 p2 = m_ip2.get_pos(date);
Vect6 p2;
m_ip2.get_pos(date, p2);
return sqrt( return sqrt(
sqr(p1.x - p2.x) + sqr(p1.x - p2.x) +
...@@ -232,11 +231,11 @@ void OrbBlocksStore::calc_min_dist(const SInitOrbit& orb1, const SInitOrbit& orb ...@@ -232,11 +231,11 @@ void OrbBlocksStore::calc_min_dist(const SInitOrbit& orb1, const SInitOrbit& orb
result_dist = eval.calc(shift); result_dist = eval.calc(shift);
} }
bool OrbBlocksStore::get_pos(const TimeJD& date, Vect6& result) const Vect6 OrbBlocksStore::get_pos(const TimeJD& date) const
{ {
for (std::vector<XInterpolator*>::const_iterator it = m_interpolators.begin(); it != m_interpolators.end(); ++it) { for (std::vector<XInterpolator*>::const_iterator it = m_interpolators.begin(); it != m_interpolators.end(); ++it) {
if (date >= (*it)->get_beg() && date <= (*it)->get_end()) { if (date >= (*it)->get_beg() && date <= (*it)->get_end()) {
return (*it)->get_pos(date, result); return (*it)->get_pos(date);
} }
} }
......
#pragma once #pragma once
#include "TM.h" #include "TM.h"
#include <vector> #include <vector>
#include "orbBlock.h" #include "orbBlock.h"
...@@ -14,7 +14,7 @@ public: ...@@ -14,7 +14,7 @@ public:
OrbBlocksStore() {}; OrbBlocksStore() {};
~OrbBlocksStore(); ~OrbBlocksStore();
void init(const SInitOrbits& orbits, const TimeJD& beg, const TimeJD& end); void init(const SInitOrbits& orbits, const TimeJD& beg, const TimeJD& end);
bool get_pos(const TimeJD& date, Vect6& result) const; Vect6 get_pos(const TimeJD& date) const;
const SInitOrbit& get_orbit(const TimeJD& date) const; const SInitOrbit& get_orbit(const TimeJD& date) const;
const OrbBlocks& get_blocks() const { return m_blocks; } const OrbBlocks& get_blocks() const { return m_blocks; }
private: private:
......
...@@ -8,8 +8,7 @@ TrjInter::TrjInter(const OrbBlocksStore& orbit_processor) ...@@ -8,8 +8,7 @@ TrjInter::TrjInter(const OrbBlocksStore& orbit_processor)
bool TrjInter::Predict(const TimeJD &t, const IFrame &fr, Vect6 &f) bool TrjInter::Predict(const TimeJD &t, const IFrame &fr, Vect6 &f)
{ {
Vect6 p; Vect6 p = m_orbit_processor.get_pos(t);
m_orbit_processor.get_pos(t, p);
p.x /= 1E6; p.x /= 1E6;
p.y /= 1E6; p.y /= 1E6;
p.z /= 1E6; p.z /= 1E6;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment