Skip to content

C
common
Commit cbaaa495 authored by Alexander Lapshin's avatar Alexander Lapshin
Browse files
Options

.

parent 76023efb
Hide whitespace changes
Inline Side-by-side
Showing with 392 additions and 313 deletions
bnolib/BNOAstroBase.cpp
View file @ cbaaa495
......@@ -660,43 +660,43 @@ void BNOAstroBase::GetPosGelio(TPlanetId pl, const TimeJD &t, Vect3 *crd, Vect3
}
}
//-----------------------------------------------------------------------------
//
void BNOAstroBase::GetPosGeo(TPlanetId pl, const TimeJD &t, Vect3 *crd, Vect3 *vel,
Vect3 *acc) const
{
if (pl == plIdEarth || pl == plIdMoon)
//
void BNOAstroBase::GetPosGeo(TPlanetId pl, const TimeJD &t, Vect3 *crd, Vect3 *vel,
Vect3 *acc) const
{
if (pl == plIdEarth || pl == plIdMoon)
{
const IEphemeride *iph = GetPlanet(pl)->GetEphemeride();
TAssert(iph->GetCS() == csIdJ2000);
iph->GetPosition(t, crd, vel, acc);
const IEphemeride *iph = GetPlanet(pl)->GetEphemeride();
TAssert(iph->GetCS() == csIdJ2000);
iph->GetPosition(t, crd, vel, acc);
}
else
else
{
Vect3 cem, vem, aem;
Vect3 cm, vm, am;
Vect3 cs, vs, as;
GetPosGelio(pl, t, crd, vel, acc);
GetPosGelio(plIdEarthMoon, t, crd == NULL ? NULL : &cem, vel == NULL ? NULL : &vem,
acc == NULL ? NULL : &aem);
GetPosGeo(plIdMoon, t, crd == NULL ? NULL : &cm, vel == NULL ? NULL : &vm,
acc == NULL ? NULL : &am);
if (crd != NULL)
{
GetGelio2GeoShift(cem, cm, cs);
(*crd) += cs;
}
if (vel != NULL)
{
GetGelio2GeoShift(vem, vm, vs);
(*vel) += vs;
}
if (acc != NULL)
{
GetGelio2GeoShift(aem, am, as);
(*acc) += as;
}
Vect3 cem, vem, aem;
Vect3 cm, vm, am;
Vect3 cs, vs, as;
GetPosGelio(pl, t, crd, vel, acc);
GetPosGelio(plIdEarthMoon, t, crd == NULL ? NULL : &cem, vel == NULL ? NULL : &vem,
acc == NULL ? NULL : &aem);
GetPosGeo(plIdMoon, t, crd == NULL ? NULL : &cm, vel == NULL ? NULL : &vm,
acc == NULL ? NULL : &am);
if (crd != NULL)
{
GetGelio2GeoShift(cem, cm, cs);
(*crd) += cs;
}
if (vel != NULL)
{
GetGelio2GeoShift(vem, vm, vs);
(*vel) += vs;
}
if (acc != NULL)
{
GetGelio2GeoShift(aem, am, as);
(*acc) += as;
}
}
}
}
//-----------------------------------------------------------------------------
//
bool BNOAstroBase::Init(bool noAtm, bool loadPlanets)
......
bnolib/IModel.h
View file @ cbaaa495
......@@ -11,130 +11,132 @@
#include "DVectors.h"
#include "Ids.h"
///////////////////////////////////////////////////////////////////////////////
// ,
//
// ,
//
class IEphemeride
{
public:
//
virtual ~IEphemeride()
{} ;
//
{
public:
//
virtual ~IEphemeride()
{};
//
virtual bool LoadData(const TimeJD &s, const TimeJD &e) = 0;
//
//
virtual void UnLoadData() = 0;
// .
virtual void GetPosition(const TimeJD &t, Vect3 *pos = NULL,
Vect3 *vel = NULL, Vect3 *acc = NULL) const = 0;
//
// .
virtual void GetPosition(const TimeJD &t, Vect3 *pos = NULL,
Vect3 *vel = NULL, Vect3 *acc = NULL) const = 0;
//
virtual TCoordSysId GetCS() const = 0;
};
};
///////////////////////////////////////////////////////////////////////////////
//
class TGeometryModel
{
public:
//
class TGeometryModel
{
public:
double Re; //
double alpha; //
};
};
///////////////////////////////////////////////////////////////////////////////
//
class IAtmosphereModel
{
public:
// .
virtual ~IAtmosphereModel()
{} ;
// .
//
virtual double Density(const TimeJD &t, const Vect3 &v, const Vect3 &vSun)
const = 0;
//
//
class IAtmosphereModel
{
public:
// .
virtual ~IAtmosphereModel()
{};
// .
//
virtual double Density(const TimeJD &t, const Vect3 &v, const Vect3 &vSun)
const = 0;
//
virtual TAtmosphereId GetID() const = 0;
// ?
// ?
virtual bool IsSunPosNeed() const = 0;
// ?
// ?
virtual bool IsTableNeed() const = 0;
};
};
///////////////////////////////////////////////////////////////////////////////
// ,
// ,
class TGeoSolarFlux
{
public:
{
public:
double f135; // 135 - 77
double f90; // 90 - ???
double f81; // 81 - ???
double f135w; // 135
double f81w; // 81 - 90
double f01; // 1 -
double kp ; // ( )
int ap ; //
double kp; // ( )
int ap; //
double three_hours_kp[8]; //
bool iskp3;
};
};
///////////////////////////////////////////////////////////////////////////////
//
//
class ITableGeoSolarFlux
{
public:
//
virtual ~ITableGeoSolarFlux()
{} ;
//
{
public:
//
virtual ~ITableGeoSolarFlux()
{};
//
virtual const TGeoSolarFlux *Get_F107(const TimeJD &tm) const = 0;
//
//
virtual const double *Get_Kp3hour(const TimeJD &tm) const = 0;
//
//
virtual double Get_Kp3hourMod(const TimeJD &tm) const = 0;
};
};
///////////////////////////////////////////////////////////////////////////////
//
class IGravityModel
{
public:
//
virtual ~IGravityModel()
{} ;
//
//
class IGravityModel
{
public:
//
virtual ~IGravityModel()
{};
//
virtual double c(unsigned int n, unsigned int m) const = 0;
virtual double d(unsigned int n, unsigned int m) const = 0;
//
//
virtual unsigned int N() const = 0;
virtual unsigned int M() const = 0;
//
//
virtual TGravityModelId GetID() const = 0;
virtual const double* const* get_cc() const { return 0; }
virtual const double* const* get_dd() const { return 0; }
double mu; //
double Re; //
double alpha; //
};
};
///////////////////////////////////////////////////////////////////////////////
//
//
class IPlanet
{
public:
//
{
public:
//
IPlanet()
{};
//
{};
//
virtual ~IPlanet()
{};
{};
//
virtual const char *Name() const = 0;
//
//
virtual const char *Name() const = 0;
//
virtual TPlanetId Id() const = 0;
//
//
virtual const IGravityModel *GetGravityModel() const = 0;
//
//
virtual const IAtmosphereModel *GetAtmosphere() const = 0;
//
//
virtual IEphemeride *GetEphemeride() = 0;
virtual const IEphemeride *GetEphemeride() const = 0;
double omega; //
double RShad; //
TGeometryModel Geometry;
};
};
///////////////////////////////////////////////////////////////////////////////
#endif
bnolib/ODESystem.cpp
View file @ cbaaa495
......@@ -7,6 +7,8 @@
#include "ODESystem.h"
#include "GeoCoords.h"
#include "IAstroBase.h"
#include "EarthGravity.h"
///////////////////////////////////////////////////////////////////////////////
//
//-----------------------------------------------------------------------------
......@@ -197,7 +199,7 @@ bool ODEPhaseSystem::GetRP(double t, const double *ph, double *f)
//-----------------------------------------------------------------------------
void ODEPhaseSystem::EarthNonCentralAcc(Vect3 &acc, const Vect3 &pos)
{
double r = pos.Len();
double r = pos.Len();
double ro = EarthGravity->Re / r;
double xr = pos.x / r;
double yr = pos.y / r;
......@@ -215,19 +217,23 @@ void ODEPhaseSystem::EarthNonCentralAcc(Vect3 &acc, const Vect3 &pos)
double x2 = 0;
double p0 = p0_1;
double cn0, t1;
const double* const* CC = EarthGravity->get_cc();
const double* const* DD = EarthGravity->get_dd();
int n;
for (n = 2; n <= Ng; n++)
{
ro1 *= ro;
Ron[n] = ro1;
cn0 = EarthGravity->c(n, 0) * ro1;
x0 += cn0 * (n+1) * p0;
//cn0 = EarthGravity->c(n, 0) * ro1;
cn0 = CC[0][n] * ro1;
x0 += cn0 * (n + 1) * p0;
x2 += cn0 * Pm[n];
if (n == Ng)
break;
t1 = (2 * n + 1) * zr;
Pm[n+1] = (t1 * Pm[n] - (n + 1) * Pm[n-1]) / n;
Pm[n + 1] = (t1 * Pm[n] - (n + 1) * Pm[n - 1]) / n;
p0 = (t1 * p0_1 - n * p0_2) / (n + 1);
p0_2 = p0_1;
p0_1 = p0;
......@@ -244,7 +250,7 @@ void ODEPhaseSystem::EarthNonCentralAcc(Vect3 &acc, const Vect3 &pos)
double z12 = 0;
int m;
double cm=0.0, sm=0.0, cc;
double cm = 0.0, sm = 0.0, cc;
for (m = 1; m <= Mg; m++)
{
double x01 = 0;
......@@ -260,14 +266,19 @@ void ODEPhaseSystem::EarthNonCentralAcc(Vect3 &acc, const Vect3 &pos)
n = m < 2 ? 2 : m;
const double *pn = Nm1;
for (; n <= Mg; n++)
{
t1 = Ron[n];
double cnm = EarthGravity->c(n, m) * t1;
double dnm = EarthGravity->d(n , m) * t1;
t1 = Pm[n];
//double cnm = EarthGravity->c(n, m) * t1;
double cnm = CC[m][n] * t1;
//double dnm = EarthGravity->d(n, m) * t1;
double dnm = DD[m][n] * t1;
t1 = Pm[n];
double t2 = cnm * t1;
x11 += t2;
t1 *= dnm;
......@@ -310,7 +321,7 @@ void ODEPhaseSystem::EarthNonCentralAcc(Vect3 &acc, const Vect3 &pos)
// ???
double *Pmt = Pm;
Pm = Pm1;
Pm = Pm1;
Pm1 = Pmt;
};
......@@ -319,7 +330,7 @@ void ODEPhaseSystem::EarthNonCentralAcc(Vect3 &acc, const Vect3 &pos)
acc.x = x11 * (z11 - xr * x01);
acc.y = x11 * (z12 - yr * x01);
acc.z = x11 * (z2 - zr * x01);
acc.z = x11 * (z2 - zr * x01);
// . J2000
......
bnolib/ODESystem.h
View file @ cbaaa495
......@@ -110,4 +110,3 @@ class ODEVarSystem : public ODEPhaseSystem
///////////////////////////////////////////////////////////////////////////////
#endif
bnolib/SGravity.h
View file @ cbaaa495
......@@ -11,157 +11,227 @@
///////////////////////////////////////////////////////////////////////////////
// ( )
///////////////////////////////////////////////////////////////////////////////
class SGravityModel : public IGravityModel
{
public:
//
class SGravityModel : public IGravityModel
{
public:
//
SGravityModel()
{};
//
{}
//
SGravityModel(double m, double r, double a)
{ mu = m; Re = r; alpha = a;};
{
mu = m; Re = r; alpha = a;
}
//
//
virtual double c(unsigned int /*n*/, unsigned int /*m*/) const
{ return 0;};
virtual double d(unsigned int /*n*/, unsigned int /*m*/) const
{ return 0;};
{
return 0;
}
//
virtual double d(unsigned int /*n*/, unsigned int /*m*/) const
{
return 0;
}
//
virtual unsigned int N() const
{ return 0;};
{
return 0;
}
virtual unsigned int M() const
{ return 0;};
{
return 0;
};
//
//
virtual TGravityModelId GetID() const
{ return grIdMu;};
};
{
return grIdMu;
};
};
///////////////////////////////////////////////////////////////////////////////
// PZ90
///////////////////////////////////////////////////////////////////////////////
class EarthPZ90GravityModel : public IGravityModel
{
public:
//
class EarthPZ90GravityModel : public IGravityModel
{
public:
//
EarthPZ90GravityModel();
//
virtual double c(unsigned int n, unsigned int m) const
{return CC[m][n];};
virtual double d(unsigned int n, unsigned int m) const
{ return DD[m][n]; };
//
virtual double c(unsigned int n, unsigned int m) const
{
return CC[m][n];
};
virtual double d(unsigned int n, unsigned int m) const
{
return DD[m][n];
};
virtual const double* const* get_cc() const { return CC; }
//
virtual unsigned int N() const
{ return 36; };
virtual unsigned int M() const
{ return 36; };
virtual const double* const* get_dd() const { return DD; }
//
//
virtual unsigned int N() const
{
return 36;
};
virtual unsigned int M() const
{
return 36;
};
//
virtual TGravityModelId GetID() const
{ return grIdPZ90;};
{
return grIdPZ90;
};
protected:
//
protected:
//
double **CC;
double **DD;
};
};
///////////////////////////////////////////////////////////////////////////////
// ( )
///////////////////////////////////////////////////////////////////////////////
class EarthPZ902GravityModel : public IGravityModel
class EarthPZ902GravityModel : public IGravityModel
{
public:
//
EarthPZ902GravityModel();
//
virtual double c(unsigned int n, unsigned int m) const
{return CC[m][n];};
virtual double d(unsigned int n, unsigned int m) const
{ return DD[m][n]; };
//
virtual unsigned int N() const
{ return 36; };
virtual unsigned int M() const
{ return 36; };
//
virtual TGravityModelId GetID() const
{ return grIdPZ902s;};
//
EarthPZ902GravityModel();
double** getCC() const {
return CC;
}
double** getDD() const {
return DD;
}
//
virtual double c(unsigned int n, unsigned int m) const
{
return CC[m][n];
};
virtual double d(unsigned int n, unsigned int m) const
{
return DD[m][n];
};
virtual const double* const* get_cc() const { return CC; }
virtual const double* const* get_dd() const { return DD; }
//
virtual unsigned int N() const
{
return 36;
};
virtual unsigned int M() const
{
return 36;
};
//
virtual TGravityModelId GetID() const
{
return grIdPZ902s;
};
protected:
//
double **CC;
double **DD;
//
double **CC;
double **DD;
};
///////////////////////////////////////////////////////////////////////////////
// ( )
///////////////////////////////////////////////////////////////////////////////
class GravityModel : public IGravityModel
{
public:
// ,
//
GravityModel(TGravityModelId Id, unsigned int N, unsigned int M,
const IGravityModel *cd);
//
class GravityModel : public IGravityModel
{
public:
// ,
//
GravityModel(TGravityModelId Id, unsigned int N, unsigned int M,
const IGravityModel *cd);
//
GravityModel();
//
//
virtual ~GravityModel()
{ freeAll();};
//
virtual double c(unsigned int n, unsigned int m) const
{return CC[m][n];};
virtual double d(unsigned int n, unsigned int m) const
{ return DD[m][n]; };
//
double c(unsigned int n, unsigned int m, double value)
{ return (CC[m][n] = value); };
double d(unsigned int n, unsigned int m, double value)
{ return (DD[m][n] = value); };
//
virtual unsigned int N() const
{ return N_; };
virtual unsigned int M() const
{ return M_; };
//
{
freeAll();
};
//
virtual double c(unsigned int n, unsigned int m) const
{
return CC[m][n];
};
virtual double d(unsigned int n, unsigned int m) const
{
return DD[m][n];
};
//
double c(unsigned int n, unsigned int m, double value)
{
return (CC[m][n] = value);
};
double d(unsigned int n, unsigned int m, double value)
{
return (DD[m][n] = value);
};
virtual const double* const* get_cc() const { return CC; }
virtual const double* const* get_dd() const { return DD; }
//
virtual unsigned int N() const
{
return N_;
};
virtual unsigned int M() const
{
return M_;
};
//
virtual TGravityModelId GetID() const
{ return Id;};
{
return Id;
};
//
//
//
//
void Init(unsigned int N, unsigned int M, const IGravityModel *c_d);
//
//
void initEmpty(TGravityModelId Id);
//
//
void null();
//
//
void print(char *comment, FILE *out = stdout) const;
//
//
TGravityModelId Id;
protected:
//
protected:
//
double **CC;
double **DD;
//
//
unsigned int N_;
unsigned int M_;
//
//
bool isInit;
//
//
void freeAll();
//
//
void SetSize(unsigned int N, unsigned int M);
};
};
///////////////////////////////////////////////////////////////////////////////
#endif
bnolib/ephemeris.cpp
View file @ cbaaa495
......@@ -151,103 +151,104 @@ bool Ephemeride::ReadData(const TimeJD &start, const TimeJD &end)
//---------------------------------------------------------------------------
void Ephemeride::CalculatePolynom(double t_cheb, double *coeff, double *coord,
double *coord_dt, double *coord_dt2) const
{
double b2 = 0.0;
double b1 = 0.0;
double b;
double b2_s = 0.0;
double b1_s = 0.0;
double b_s = 0.0;
double b2_ss = 0.0;
double b1_ss = 0.0;
double b_ss = 0.0;
for (int i = Order[BodyID]; i > 0; --i)
{
double b2 = 0.0;
double b1 = 0.0;
double b;
double b2_s = 0.0;
double b1_s = 0.0;
double b_s = 0.0;
double b2_ss = 0.0;
double b1_ss = 0.0;
double b_ss = 0.0;
for (int i = Order[BodyID]; i > 0; --i)
{
b = coeff[i] - b2 + 2.0 * b1 * t_cheb;
b2 = b1;
b1 = b;
if (coord_dt || coord_dt2)
{
b_s = 2.0 * (b1_s * t_cheb + b2) - b2_s;
b2_s = b1_s;
b1_s = b_s;
}
if (coord_dt2)
{
b_ss = 2.0 * (b1_ss * t_cheb + 2.0 * b2_s) - b2_ss;
b2_ss = b1_ss;
b1_ss = b_ss;
}
b = coeff[i] - b2 + 2.0 * b1 * t_cheb;
b2 = b1;
b1 = b;
if (coord_dt || coord_dt2)
{
b_s = 2.0 * (b1_s * t_cheb + b2) - b2_s;
b2_s = b1_s;
b1_s = b_s;
}
if (coord_dt2)
{
b_ss = 2.0 * (b1_ss * t_cheb + 2.0 * b2_s) - b2_ss;
b2_ss = b1_ss;
b1_ss = b_ss;
}
}
if (coord)
*coord = coeff[0] - b2 + b1 * t_cheb;
if (coord_dt)
*coord_dt = b1 + b_s * t_cheb - b2_s;
if (coord_dt2)
*coord_dt2= 2.0 * b1_s + b1_ss * t_cheb - b2_ss;
}
if (coord)
*coord = coeff[0] - b2 + b1 * t_cheb;
if (coord_dt)
*coord_dt = b1 + b_s * t_cheb - b2_s;
if (coord_dt2)
*coord_dt2 = 2.0 * b1_s + b1_ss * t_cheb - b2_ss;
}
//---------------------------------------------------------------------------
void Ephemeride::GetPosition(const TimeJD &t, Vect3 *pos,
Vect3 *vel, Vect3 *acc) const
{
TimeJD tdb(t);
TAssert(XCoefs != NULL);
TAssert(tdb >= StartTime && tdb <= EndTime);
if (GetFromCash(t, pos, vel, acc))
return;
double diff = (tdb - StartTime).GetDouble();
double tch = 2.0 * fmod(diff / DInterval[BodyID], 1.0) - 1.0;
int offset = (int)floor(diff / DInterval[BodyID]) * (Order[BodyID] + 1);
CalculatePolynom(tch, XCoefs + offset,
pos == NULL ? NULL : &pos->x,
vel == NULL ? NULL : &vel->x,
acc == NULL ? NULL : &acc->x);
CalculatePolynom(tch, YCoefs + offset,
pos == NULL ? NULL : &pos->y,
vel == NULL ? NULL : &vel->y,
acc == NULL ? NULL : &acc->y);
CalculatePolynom(tch, ZCoefs + offset,
pos == NULL ? NULL : &pos->z,
vel == NULL ? NULL : &vel->z,
acc == NULL ? NULL : &acc->z);
double dt_scale = 2 / (86400.0 * DInterval[BodyID]);
double dt2_scale = dt_scale * dt_scale * 1.e3;
if (pos)
void Ephemeride::GetPosition(const TimeJD &t, Vect3 *pos,
Vect3 *vel, Vect3 *acc) const
{
// TAssert(XCoefs != NULL);
// TAssert(tdb >= StartTime && tdb <= EndTime);
if (GetFromCash(t, pos, vel, acc))
return;
TimeJD tdb(t);
double diff = (tdb - StartTime).GetDouble();
double tch = 2.0 * fmod(diff / DInterval[BodyID], 1.0) - 1.0;
int offset = (int)floor(diff / DInterval[BodyID]) * (Order[BodyID] + 1);
CalculatePolynom(tch, XCoefs + offset,
pos == NULL ? NULL : &pos->x,
vel == NULL ? NULL : &vel->x,
acc == NULL ? NULL : &acc->x);
CalculatePolynom(tch, YCoefs + offset,
pos == NULL ? NULL : &pos->y,
vel == NULL ? NULL : &vel->y,
acc == NULL ? NULL : &acc->y);
CalculatePolynom(tch, ZCoefs + offset,
pos == NULL ? NULL : &pos->z,
vel == NULL ? NULL : &vel->z,
acc == NULL ? NULL : &acc->z);
double dt_scale = 2 / (86400.0 * DInterval[BodyID]);
double dt2_scale = dt_scale * dt_scale * 1.e3;
if (pos)
{
pos->x *= 1.e-3;
pos->y *= 1.e-3;
pos->z *= 1.e-3;
pos->x *= 1.e-3;
pos->y *= 1.e-3;
pos->z *= 1.e-3;
}
if (vel)
if (vel)
{
vel->x *= dt_scale;
vel->y *= dt_scale;
vel->z *= dt_scale;
vel->x *= dt_scale;
vel->y *= dt_scale;
vel->z *= dt_scale;
}
if (acc)
if (acc)
{
acc->x *= dt2_scale;
acc->y *= dt2_scale;
acc->z *= dt2_scale;
acc->x *= dt2_scale;
acc->y *= dt2_scale;
acc->z *= dt2_scale;
}
PutToCash(t, pos, vel, acc);
PutToCash(t, pos, vel, acc);
return;
}
return;
}
//---------------------------------------------------------------------------
//
bool Ephemeride::GetFromCash(const TimeJD &t, Vect3 *pos, Vect3 *vel, Vect3 *acc) const
......
utils/TextFunctions.h
View file @ cbaaa495
......@@ -185,7 +185,7 @@ static bool ReadFile(std::string path, std::string &str)
size_t res = fread(buf,len,1,fp); //read into buffer
fclose(fp);
str = std::string(buf);
delete buf;
delete[] buf;
if (!res){
return false;
......
utils/interpolate.cpp
View file @ cbaaa495
......@@ -104,10 +104,6 @@ void XInterpolator::calc_positions(const SInitOrbit& orbit, const TimeJD& beg, c
TimeJD beg_ext = ShiftDate(beg, -n*step);
TimeJD end_ext = ShiftDate(end, +n*step);
// std::string s_beg = DateToStr(beg);
// std::string s_date = DateToStr(orbDate);
// std::string s_end = DateToStr(end);
if (orbit.GetDate() <= beg_ext){
propagate(orbit, beg_ext, end_ext, 1, result);
}else if (orbit.GetDate() >= end_ext){
......
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