Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • GASpiridonov/satellite-motion-simulation
1 result
Show changes
Commits on Source (2)
Hide whitespace changes
Inline Side-by-side
Showing
with 759 additions and 165 deletions
Undisturbed_movement/Matrix.cpp
View file @ aa4949cc
#include "Matrix.h" #include "Matrix.h"
Matrix operator+(Matrix& matrix_1, Matrix& matrix_2) Matrix operator+(const Matrix& matrix_1, const Matrix& matrix_2)
{ {
Matrix output(matrix_1.m_rows, matrix_1.m_cols); Matrix output(matrix_1.m_rows, matrix_1.m_cols);
for (int i = 0; i < output.m_rows; i++) { for (uint32_t i = 0; i < output.m_rows; i++) {
for (int j = 0; j < output.m_cols; j++) for (uint32_t j = 0; j < output.m_cols; j++)
output.at(i, j) = matrix_1.at(i, j) + matrix_2.at(i,j); output.at(i, j) = matrix_1.at(i, j) + matrix_2.at(i,j);
} }
return output; return output;
} }
Matrix operator+(Matrix& matrix_1, double in){ Matrix operator+(const Matrix& matrix_1, const double in){
Matrix buff{ matrix_1.m_rows, matrix_1.m_cols, in }; Matrix buff{ matrix_1.m_rows, matrix_1.m_cols, in };
return (matrix_1+buff); return (matrix_1+buff);
} }
Matrix operator+(double in, Matrix& matrix_1){ Matrix operator+(const double in, const Matrix& matrix_1){
return (matrix_1 + in); return (matrix_1 + in);
} }
Matrix operator*(Matrix& matrix_1, Matrix& matrix_2) { Matrix operator-(const Matrix& matrix_1, const Matrix& matrix_2)
{
Matrix output(matrix_1.m_rows, matrix_1.m_cols);
for (uint32_t i = 0; i < output.m_rows; i++) {
for (uint32_t j = 0; j < output.m_cols; j++)
output.at(i, j) = matrix_1.at(i, j) - matrix_2.at(i, j);
}
return output;
}
Matrix operator-(const Matrix& matrix_1, const double in) {
Matrix buff{ matrix_1.m_rows, matrix_1.m_cols, in };
return (matrix_1 - buff);
}
Matrix operator-(const double in, const Matrix& matrix_1) {
Matrix buff{ matrix_1.m_rows, matrix_1.m_cols, in };
return (buff - matrix_1);
}
Matrix operator*(const Matrix& matrix_1, const Matrix& matrix_2) {
assert((matrix_1.m_cols == matrix_2.m_rows)); assert((matrix_1.m_cols == matrix_2.m_rows));
Matrix output{ matrix_1.m_rows, matrix_2.m_cols }; Matrix output{ matrix_1.m_rows, matrix_2.m_cols };
for (int i = 0; i < output.m_rows; i++) { for (uint32_t i = 0; i < output.m_rows; i++) {
for (int j = 0; j < output.m_cols; j++) { for (uint32_t j = 0; j < output.m_cols; j++) {
double result = 0; double result = 0;
for (int k = 0; k < matrix_1.m_cols; k++) for (uint32_t k = 0; k < matrix_1.m_cols; k++)
result += matrix_1.at(i, k) * matrix_2.at(k, j); result += matrix_1.at(i, k) * matrix_2.at(k, j);
output.at(i, j) = result; output.at(i, j) = result;
} }
...@@ -33,20 +54,28 @@ Matrix operator*(Matrix& matrix_1, Matrix& matrix_2) { ...@@ -33,20 +54,28 @@ Matrix operator*(Matrix& matrix_1, Matrix& matrix_2) {
return output; return output;
} }
Matrix operator*(Matrix& matrix_1, double in) { Matrix operator*(const Matrix& matrix_1, const double in) {
Matrix output{ matrix_1.m_cols, matrix_1.m_rows }; Matrix output{ matrix_1.m_rows, matrix_1.m_cols };
for (int i = 0; i < output.m_rows; i++) { for (uint32_t i = 0; i < output.m_rows; i++) {
for (int j = 0; j < output.m_cols; j++) { for (uint32_t j = 0; j < output.m_cols; j++) {
output.at(i, j) = matrix_1.at(i, j) * in; output.at(i, j) = matrix_1.at(i, j) * in;
} }
} }
return output; return output;
} }
Matrix operator*(double in, Matrix& matrix_1) { Matrix operator*(const double in, const Matrix& matrix_1) {
return (matrix_1*in); return (matrix_1*in);
} }
Matrix operator/(const Matrix& matrix_1, const double in) {
return (matrix_1 * (1 / in));
}
Matrix operator/(const double in, const Matrix& matrix_1) {
return ((1 / in) * matrix_1);
}
Matrix Matrix::operator=(const Matrix& input_matrix) Matrix Matrix::operator=(const Matrix& input_matrix)
{ {
if (this == &input_matrix) return *this; if (this == &input_matrix) return *this;
...@@ -58,10 +87,10 @@ Matrix Matrix::operator=(const Matrix& input_matrix) ...@@ -58,10 +87,10 @@ Matrix Matrix::operator=(const Matrix& input_matrix)
} }
std::ostream& operator<<(std::ostream& out, Matrix& matrix_1) { std::ostream& operator<<(std::ostream& out, const Matrix& matrix_1) {
out << "\n"; out << "\n";
for (int i = 0; i < matrix_1.m_rows; i++) { for (uint32_t i = 0; i < matrix_1.m_rows; i++) {
for (int j = 0; j < matrix_1.m_cols; j++) { for (uint32_t j = 0; j < matrix_1.m_cols; j++) {
out << matrix_1(i, j) << " "; out << matrix_1(i, j) << " ";
} }
out << "\n"; out << "\n";
...@@ -69,19 +98,24 @@ std::ostream& operator<<(std::ostream& out, Matrix& matrix_1) { ...@@ -69,19 +98,24 @@ std::ostream& operator<<(std::ostream& out, Matrix& matrix_1) {
return out; return out;
} }
std::ofstream& operator<<(std::ofstream& out, Matrix& i_matrix) //std::ofstream& operator<<(std::ofstream& out, const Matrix& i_matrix)
//{
// out << "\n";
// for (uint32_t i = 0; i < i_matrix.m_rows; i++) {
// for (uint32_t j = 0; j <i_matrix.m_cols; j++) {
// out << i_matrix.at(i, j) << " ";
// }
// out << "\n";
// }
// return out;
//}
double& Matrix::operator()(uint32_t input_rows, uint32_t input_cols)
{ {
out << "\n"; return (at(input_rows, input_cols));
for (int i = 0; i < i_matrix.m_rows; i++) {
for (int j = 0; j <i_matrix.m_cols; j++) {
out << i_matrix(i, j) << " ";
}
out << "\n";
}
return out;
} }
double& Matrix::operator()(uint32_t input_rows, uint32_t input_cols) double Matrix::operator()(uint32_t input_rows, uint32_t input_cols) const
{ {
return (at(input_rows, input_cols)); return (at(input_rows, input_cols));
} }
...@@ -90,8 +124,8 @@ double& Matrix::operator()(uint32_t input_rows, uint32_t input_cols) ...@@ -90,8 +124,8 @@ double& Matrix::operator()(uint32_t input_rows, uint32_t input_cols)
Matrix Matrix::Transpose() Matrix Matrix::Transpose()
{ {
Matrix output{m_cols, m_rows }; Matrix output{m_cols, m_rows };
for (int i = 0; i < output.m_rows; i++) { for (uint32_t i = 0; i < output.m_rows; i++) {
for (int j = 0; j < output.m_cols; j++) for (uint32_t j = 0; j < output.m_cols; j++)
output.at(i, j) = at(j, i); output.at(i, j) = at(j, i);
} }
return output; return output;
......
Undisturbed_movement/Matrix.h
View file @ aa4949cc
...@@ -27,20 +27,35 @@ public: ...@@ -27,20 +27,35 @@ public:
return m_matrix.at((input_rows * m_cols) + input_cols); return m_matrix.at((input_rows * m_cols) + input_cols);
} }
double at(uint32_t input_rows, uint32_t input_cols) const {
assert(input_rows < m_rows&& input_cols < m_cols);
return m_matrix.at((input_rows * m_cols) + input_cols);
}
int add(double input_number) { int add(double input_number) {
m_matrix.push_back(input_number); m_matrix.push_back(input_number);
return 0; return 0;
} }
double& operator() (uint32_t input_rows, uint32_t input_cols); double& operator() (uint32_t input_rows, uint32_t input_cols);
double operator() (uint32_t input_rows, uint32_t input_cols) const;
friend Matrix operator+(const Matrix& matrix_1, const Matrix& matrix_2);
friend Matrix operator+(const Matrix& matrix_1, const double in);
friend Matrix operator+(const double in, const Matrix& matrix_1);
friend Matrix operator-(const Matrix& matrix_1, const Matrix& matrix_2);
friend Matrix operator-(const Matrix& matrix_1, const double in);
friend Matrix operator-(const double in, const Matrix& matrix_1);
friend Matrix operator*(const Matrix& matrix_1, const Matrix& matrix_2);
friend Matrix operator*(const Matrix& matrix_1, double in);
friend Matrix operator*(double in, const Matrix& matrix_1);
friend Matrix operator+(Matrix& matrix_1, Matrix& matrix_2); friend Matrix operator/(const Matrix& matrix_1, const double in);
friend Matrix operator+(Matrix& matrix_1, double in); friend Matrix operator/(const double in, const Matrix& matrix_1);
friend Matrix operator+(double in, Matrix& matrix_1);
friend Matrix operator*(Matrix& matrix_1, Matrix& matrix_2);
friend Matrix operator*(Matrix& matrix_1, double in);
friend Matrix operator*(double in, Matrix& matrix_1);
Matrix operator= (const Matrix& input_matrix); Matrix operator= (const Matrix& input_matrix);
Matrix Transpose(); Matrix Transpose();
...@@ -48,8 +63,8 @@ public: ...@@ -48,8 +63,8 @@ public:
static Matrix RotationOY(double α); static Matrix RotationOY(double α);
static Matrix RotationOZ(double α); static Matrix RotationOZ(double α);
friend std::ostream& operator<<(std::ostream& out, Matrix& i_matrix); friend std::ostream& operator<<(std::ostream& out, const Matrix& i_matrix);
friend std::ofstream& operator<<(std::ofstream& out, Matrix& i_matrix); friend std::ofstream& operator<<(std::ofstream& out, const Matrix& i_matrix);
}; };
......
Undisturbed_movement/Planets.cpp
View file @ aa4949cc
...@@ -88,8 +88,9 @@ int Moon::calculate_β() ...@@ -88,8 +88,9 @@ int Moon::calculate_β()
return 0; return 0;
} }
int Moon::calculateGEC() int Moon::calculateGEC(double input_TDB)
{ {
recalc(input_TDB);
calculateСorrecteFundamentalArg(); calculateСorrecteFundamentalArg();
calculate_Δλ_c(); calculate_Δλ_c();
calculate_β(); calculate_β();
...@@ -97,6 +98,7 @@ int Moon::calculateGEC() ...@@ -97,6 +98,7 @@ int Moon::calculateGEC()
m_GEC(0, 0) = m_cfa.λ + RS * m_Δλ_c; m_GEC(0, 0) = m_cfa.λ + RS * m_Δλ_c;
m_GEC(1, 0) = RS * m_β; m_GEC(1, 0) = RS * m_β;
m_GEC(2, 0) = m_scalefactor * (1 / (RS * m_sin_πс)); m_GEC(2, 0) = m_scalefactor * (1 / (RS * m_sin_πс));
calculateDC();
return 0; return 0;
} }
...@@ -110,36 +112,61 @@ int Moon::calculateСorrecteFundamentalArg() { ...@@ -110,36 +112,61 @@ int Moon::calculateСorrecteFundamentalArg() {
return 0; return 0;
} }
std::ostream& operator<<(std::ostream& out, Moon& i_moon) int Moon::calculateRm()
{ {
out << "\n******************************************" m_Rm = sqrt(pow(m_DC(0, 0), 2) + pow(m_DC(1, 0), 2) + pow(m_DC(2, 0), 2));
<< "\n class Moon" return 0;
<< "\n******************************************" }
<< "\nm_s = " << i_moon.m_s.at(0) << " " << i_moon.m_s.at(1) << " " << i_moon.m_s.at(2) << " "
<< i_moon.m_s.at(3) << " " << i_moon.m_s.at(4) << " " << i_moon.m_s.at(5) << " " int Moon::calculateRo(const Matrix& DC_coord_object)
<< i_moon.m_s.at(6) << " " {
<< "\nr_s = " << i_moon.m_r.at(0) << " " << i_moon.m_r.at(1) << " " m_Ro = sqrt(pow(m_DC(0, 0) - DC_coord_object(0, 0), 2) +
<< i_moon.m_r.at(2) << " " << i_moon.m_r.at(3) << " " << i_moon.m_r.at(4) pow(m_DC(1, 0) - DC_coord_object(1, 0), 2) +
<< i_moon.m_cfa << "\nm_Δλ_c = " << i_moon.m_Δλ_c << "\n\nm_β = " << i_moon.m_β << "\n\nm_sin_πс = " << i_moon.m_sin_πс pow(m_DC(2, 0), 2) - DC_coord_object(2, 0));
<< "\n\nm_GEC: " << i_moon.m_GEC << "\nm_dc:" << i_moon.m_DC; return 0;
return out;
} }
std::ofstream& operator<<(std::ofstream& out, Moon& i_moon) Matrix Moon::calculateFm(const Matrix& DC_coord_object, double input_TDB)
{ {
out << "\n******************************************" calculateGEC(input_TDB);
<< "\n class Moon" calculateRm();
<< "\n******************************************" calculateRo(DC_coord_object);
<< "\nm_s = " << i_moon.m_s.at(0) << " " << i_moon.m_s.at(1) << " " << i_moon.m_s.at(2) << " " m_F(0, 0) = gravitational_parameter_Moon * (((m_DC(0, 0) - DC_coord_object(0, 0)) / pow(m_Ro, 3)) - (m_DC(0, 0) / pow(m_Rm, 3)));
<< i_moon.m_s.at(3) << " " << i_moon.m_s.at(4) << " " << i_moon.m_s.at(5) << " " m_F(1, 0) = gravitational_parameter_Moon * (((m_DC(1, 0) - DC_coord_object(1, 0)) / pow(m_Ro, 3)) - (m_DC(1, 0) / pow(m_Rm, 3)));
<< i_moon.m_s.at(6) << " " m_F(2, 0) = gravitational_parameter_Moon * (((m_DC(2, 0) - DC_coord_object(2, 0)) / pow(m_Ro, 3)) - (m_DC(2, 0) / pow(m_Rm, 3)));
<< "\nr_s = " << i_moon.m_r.at(0) << " " << i_moon.m_r.at(1) << " " return m_F;
<< i_moon.m_r.at(2) << " " << i_moon.m_r.at(3) << " " << i_moon.m_r.at(4)
<< i_moon.m_cfa << "\nm_Δλ_c = " << i_moon.m_Δλ_c << "\n\nm_β = " << i_moon.m_β << "\n\nm_sin_πс = " << i_moon.m_sin_πс
<< "\n\nm_GEC: " << i_moon.m_GEC << "\nm_dc:" << i_moon.m_DC;
return out;
} }
//std::ostream& operator<<(std::ostream& out, Moon& i_moon)
//{
// out << "\n******************************************"
// << "\n class Moon"
// << "\n******************************************"
// << "\nm_s = " << i_moon.m_s.at(0) << " " << i_moon.m_s.at(1) << " " << i_moon.m_s.at(2) << " "
// << i_moon.m_s.at(3) << " " << i_moon.m_s.at(4) << " " << i_moon.m_s.at(5) << " "
// << i_moon.m_s.at(6) << " "
// << "\nr_s = " << i_moon.m_r.at(0) << " " << i_moon.m_r.at(1) << " "
// << i_moon.m_r.at(2) << " " << i_moon.m_r.at(3) << " " << i_moon.m_r.at(4)
// << i_moon.m_cfa << "\nm_Δλ_c = " << i_moon.m_Δλ_c << "\n\nm_β = " << i_moon.m_β << "\n\nm_sin_πс = " << i_moon.m_sin_πс
// << "\n\nm_GEC: " << i_moon.m_GEC << "\nm_dc:" << i_moon.m_DC;
// return out;
//}
//
//std::ofstream& operator<<(std::ofstream& out, Moon& i_moon)
//{
// out << "\n******************************************"
// << "\n class Moon"
// << "\n******************************************"
// << "\nm_s = " << i_moon.m_s.at(0) << " " << i_moon.m_s.at(1) << " " << i_moon.m_s.at(2) << " "
// << i_moon.m_s.at(3) << " " << i_moon.m_s.at(4) << " " << i_moon.m_s.at(5) << " "
// << i_moon.m_s.at(6) << " "
// << "\nr_s = " << i_moon.m_r.at(0) << " " << i_moon.m_r.at(1) << " "
// << i_moon.m_r.at(2) << " " << i_moon.m_r.at(3) << " " << i_moon.m_r.at(4)
// << i_moon.m_cfa << "\nm_Δλ_c = " << i_moon.m_Δλ_c << "\n\nm_β = " << i_moon.m_β << "\n\nm_sin_πс = " << i_moon.m_sin_πс
// << "\n\nm_GEC: " << i_moon.m_GEC << "\nm_dc:" << i_moon.m_DC;
// return out;
//}
double Sun::DILong_V() double Sun::DILong_V()
{ {
std::vector<double> earth_c{ 2 * m_M.at(1), 3 * m_M.at(1), 4 * m_M.at(1), 5 * m_M.at(1), 6 * m_M.at(1), 7 * m_M.at(1) }; std::vector<double> earth_c{ 2 * m_M.at(1), 3 * m_M.at(1), 4 * m_M.at(1), 5 * m_M.at(1), 6 * m_M.at(1), 7 * m_M.at(1) };
...@@ -317,6 +344,44 @@ double Sun::DIR_JS() { ...@@ -317,6 +344,44 @@ double Sun::DIR_JS() {
return output; return output;
} }
int Sun::calculateRs()
{
m_Rs = sqrt(pow(m_DC(0, 0), 2) + pow(m_DC(1, 0), 2) + pow(m_DC(2, 0), 2));
return 0;
}
int Sun::calculateRo(const Matrix& DC_coord_object)
{
m_Ro = sqrt(pow(m_DC(0, 0) - DC_coord_object(0, 0), 2) +
pow(m_DC(1, 0) - DC_coord_object(1, 0), 2) +
pow(m_DC(2, 0), 2) - DC_coord_object(2, 0));
return 0;
}
Matrix& Sun::calculateFs(const Matrix& DC_coord_object, double input_TDB)
{
calculateGEC(input_TDB);
calculateRs();
calculateRo(DC_coord_object);
m_F(0, 0) = gravitational_parameter_Sun * (((m_DC(0, 0) - DC_coord_object(0, 0)) / pow(m_Ro, 3)) - (m_DC(0, 0) / pow(m_Rs, 3)));
m_F(1, 0) = gravitational_parameter_Sun * (((m_DC(1, 0) - DC_coord_object(1, 0)) / pow(m_Ro, 3)) - (m_DC(1, 0) / pow(m_Rs, 3)));
m_F(2, 0) = gravitational_parameter_Sun * (((m_DC(2, 0) - DC_coord_object(2, 0)) / pow(m_Ro, 3)) - (m_DC(2, 0) / pow(m_Rs, 3)));
return m_F;
}
Matrix& Sun::calculatePs(const Matrix& DC_coord_object, double kr, double As, double ms, double input_TDB)
{
calculateGEC(input_TDB);
calculateRs();
calculateRo(DC_coord_object);
double C = (0.001 * light_pressure * kr * (As / ms)) * pow(((1.4959787 * pow(10, 8)) / m_Ro), 2);
m_P(0, 0) = C * ((DC_coord_object(0,0) - m_DC(0, 0))/ m_Ro);
m_P(1, 0) = C * ((DC_coord_object(1, 0) - m_DC(0, 0)) / m_Ro);
m_P(2, 0) = C * ((DC_coord_object(2, 0) - m_DC(2, 0)) / m_Ro);
return m_P;
}
int Sun::calculateAverageAnomalies() int Sun::calculateAverageAnomalies()
{ {
m_M.at(0) = 0.87167007 + 1021.32292307 * m_ts; m_M.at(0) = 0.87167007 + 1021.32292307 * m_ts;
...@@ -344,13 +409,15 @@ int Sun::calculatePerturbation() ...@@ -344,13 +409,15 @@ int Sun::calculatePerturbation()
return 0; return 0;
} }
int Sun::calculateGEC() int Sun::calculateGEC(double input_TDB)
{ {
recalc(input_TDB);
calculateAverageAnomalies(); calculateAverageAnomalies();
calculatePerturbation(); calculatePerturbation();
m_GEC(0, 0) = 4.93823996 + m_M.at(1) + RS * ((6191.2 + 1.1 * m_ts) * m_ts + m_perturbation(0, 0)); m_GEC(0, 0) = 4.93823996 + m_M.at(1) + RS * ((6191.2 + 1.1 * m_ts) * m_ts + m_perturbation(0, 0));
m_GEC(1, 0) = RS * m_perturbation(1, 0); m_GEC(1, 0) = RS * m_perturbation(1, 0);
m_GEC(2, 0) = AU * (1.0001398 - 0.0000007 * m_ts + pow(10, -6) * m_perturbation(2, 0)); m_GEC(2, 0) = AU * (1.0001398 - 0.0000007 * m_ts + pow(10, -6) * m_perturbation(2, 0));
calculateDC();
return 0; return 0;
} }
...@@ -372,11 +439,15 @@ Matrix Earth::fill_matrix(const char fileName[]) ...@@ -372,11 +439,15 @@ Matrix Earth::fill_matrix(const char fileName[])
str >> buff1; str >> buff1;
rows = stoi(buff1); rows = stoi(buff1);
str >> buff1; str >> buff1;
cols = stoi(buff1); cols = stoi(buff1);
Matrix output{ rows, cols};
m_rows = (rows >= cols) ? rows : cols;
m_cols = (rows >= cols) ? rows : cols;
Matrix output{ m_rows, m_cols};
rows = 0; rows = 0;
cols = 0; cols = 2;
while (!fin.eof()) { while (!fin.eof()) {
std::getline(fin, buff); std::getline(fin, buff);
...@@ -386,37 +457,202 @@ Matrix Earth::fill_matrix(const char fileName[]) ...@@ -386,37 +457,202 @@ Matrix Earth::fill_matrix(const char fileName[])
output(rows, cols) = stod(buff1); output(rows, cols) = stod(buff1);
cols++; cols++;
} }
cols = 0; cols = 2;
rows++; rows++;
} }
fin.close(); fin.close();
return output; return output;
} }
Earth::Earth(double time_MJD, Matrix input_GCC) : ConvertSC{ time_MJD }, m_GDC{ CC_to_TC(m_precession, m_nutation, input_GCC) } { int Earth::setF()
r = sqrt(pow(m_GDC(0,0), 2) + pow(m_GDC(1,0), 2) + pow(m_GDC(2,0), 2)); {
for(short i = 0; i < 3; i ++)
m_F(i, 0) = -gravitational_parameter_Earth * (m_GDC(i, 0) / (pow(m_r, 3)));
return 0;
}
int Earth::setRatio()
{
m_ratio.Rr = Radius_Earth / m_r;
m_ratio.Xr = m_GEC(0, 0) / m_r;
m_ratio.Yr = m_GEC(1, 0) / m_r;
m_ratio.Zr = m_GEC(2, 0) / m_r;
return 0;
} }
int Earth::calculateStartPartialDeriv() int Earth::calculatePartialDeriv()
{ {
//вычисление частных производных по x //вычисление частных производных по x
startPartialDeriv(0, 0) = -(m_GDC(0, 0)) / pow(m_r, 3); PartialDeriv(0, 0) = -(m_GDC(0, 0)) / pow(m_r, 3);
startPartialDeriv(1, 0) = (1 / m_r) - pow(m_GDC(0, 0), 2) / pow(m_r, 3); PartialDeriv(1, 0) = (1 / m_r) - pow(m_GDC(0, 0), 2) / pow(m_r, 3);
startPartialDeriv(2, 0) = -(m_GDC(0, 0) * m_GDC(1, 0)) / pow(m_r, 3); PartialDeriv(2, 0) = -(m_GDC(0, 0) * m_GDC(1, 0)) / pow(m_r, 3);
startPartialDeriv(3, 0) = -(m_GDC(0, 0) * m_GDC(2, 0)) / pow(m_r, 3); PartialDeriv(3, 0) = -(m_GDC(0, 0) * m_GDC(2, 0)) / pow(m_r, 3);
//вычисление частныхм производных по y //вычисление частныхм производных по y
startPartialDeriv(0, 1) = -(m_GDC(1, 0)) / pow(m_r, 3); PartialDeriv(0, 1) = -(m_GDC(1, 0)) / pow(m_r, 3);
startPartialDeriv(1, 1) = -(m_GDC(0, 0) * m_GDC(1, 0)) / pow(m_r, 3); PartialDeriv(1, 1) = -(m_GDC(0, 0) * m_GDC(1, 0)) / pow(m_r, 3);
startPartialDeriv(2, 1) = (1 / m_r) - pow(m_GDC(1, 0), 2) / pow(m_r, 3); PartialDeriv(2, 1) = (1 / m_r) - pow(m_GDC(1, 0), 2) / pow(m_r, 3);
startPartialDeriv(3, 1) = -(m_GDC(1, 0) * m_GDC(2, 0)) / pow(m_r, 3); PartialDeriv(3, 1) = -(m_GDC(1, 0) * m_GDC(2, 0)) / pow(m_r, 3);
//вычисление частных производных по z //вычисление частных производных по z
startPartialDeriv(0, 1) = -(m_GDC(2, 0)) / pow(m_r, 3); PartialDeriv(0, 1) = -(m_GDC(2, 0)) / pow(m_r, 3);
startPartialDeriv(1, 1) = -(m_GDC(0, 0) * m_GDC(2, 0)) / pow(m_r, 3); PartialDeriv(1, 1) = -(m_GDC(0, 0) * m_GDC(2, 0)) / pow(m_r, 3);
startPartialDeriv(2, 1) = -(m_GDC(1, 0) * m_GDC(2, 0)) / pow(m_r, 3); PartialDeriv(2, 1) = -(m_GDC(1, 0) * m_GDC(2, 0)) / pow(m_r, 3);
startPartialDeriv(3, 1) = (1 / m_r) - pow(m_GDC(2, 0), 2) / pow(m_r, 3); PartialDeriv(3, 1) = (1 / m_r) - pow(m_GDC(2, 0), 2) / pow(m_r, 3);
return 0;
}
int Earth::calculateStartParameters()
{
m_R(0, 2) = (gravitational_parameter_Earth / m_r) * pow(m_ratio.Rr, 2);
m_Partial_R(0, 2) = 3 * gravitational_parameter_Earth * pow(m_ratio.Rr, 2);
m_Z(0, 1) = m_ratio.Zr;
m_Z(0, 2) = (3 / 2) * pow(m_ratio.Zr, 2) - (1 / 2);
m_Partial_Z(0, 1) = 1;
m_Partial_Z(0, 2) = 3 * m_ratio.Zr;
m_X(0, 0) = 1;
m_Y(0, 0) = 0;
return 0;
}
int Earth::calculateAccelerationByEarth()
{
calculatePartialDeriv();
calculateStartParameters();
setRatio();
setF();
calculateRn();
calculatePartialRn();
calculateZn0();
calculatePartialZn0();
for (int k = 1; k < m_rows; k++) {
calculateZn_k(k);
calculatePartialZn_k(k);
calculateXY_k(k);
claculatePartialXxy_k(k);
claculatePartialYxy_k(k);
calculateQn_k(k);
calculatePartialQn_k(k);
calculateF(k);
}
return 0;
}
int Earth::calculateRn()
{
for (int i = 3; i < m_cols; i++)
m_R(0, i) = m_ratio.Rr * m_R(0, i - 1);
return 0;
}
int Earth::calculatePartialRn()
{
for (int i = 3; i < m_cols; i++)
m_Partial_R(0, i) = (i + 1) * gravitational_parameter_Earth * pow(m_ratio.Rr, i);
return 0; return 0;
} }
int Earth::calculateZn0()
{
for (int i = 3; i < m_cols; i++)
m_Z(0, i) = (double(2 * i - 1) / i) * m_ratio.Zr * m_Z(0, i - 1) - (double(i - 1) / i) * m_Z(0, i - 2);
return 0;
}
int Earth::calculatePartialZn0()
{
for (int i = 3; i < m_cols; i++)
m_Partial_Z(0, i) = i * m_Z(0, i - 1) + m_ratio.Zr * m_Partial_Z(0, i - 1);
return 0;
}
int Earth::calculateZn_k(int k)
{
for (int i = 2; i < m_cols; i++)
m_Z(k, i) = m_Partial_Z(k - 1, i);
return 0;
}
int Earth::calculatePartialZn_k(int k)
{
for (int i = 0; i <= k; i++)
m_Partial_Z(k, i) = 0;
for (int i = k + 1; i < m_cols; i++)
m_Partial_Z(k, i) = (2 * i - 1) * m_Z(k, i - 1) + m_Partial_Z(k, i - 2);
return 0;
}
int Earth::calculateXY_k(int k)
{
m_X(k, 0) = m_X(k - 1, 0) * m_ratio.Xr - m_Y(k - 1, 0) * m_ratio.Yr;
m_Y(k, 0) = m_Y(k - 1, 0) * m_ratio.Xr + m_X(k - 1, 0) * m_ratio.Yr;
return 0;
}
int Earth::claculatePartialXxy_k(int k)
{
m_Partial_Xx(k, 0) = m_Partial_Xx(k - 1, 0) * m_ratio.Xr + m_X(k - 1, 0) - m_Partial_Yx(k - 1, 0) * m_ratio.Yr;
m_Partial_Xy(k, 0) = m_Partial_Xy(k - 1, 0) * m_ratio.Xr - m_Partial_Yy(k - 1, 0) * m_ratio.Yr - m_Y(k-1, 0);
return 0;
}
int Earth::claculatePartialYxy_k(int k) {
m_Partial_Yx(k, 0) = m_Partial_Yx(k - 1, 0) * m_ratio.Xr + m_Y(k - 1, 0) - m_Partial_Xx(k - 1, 0) * m_ratio.Yr;
m_Partial_Yy(k, 0) = m_Partial_Yy(k - 1, 0) * m_ratio.Xr + m_Partial_Xy(k - 1, 0) * m_ratio.Yr + m_X(k - 1, 0);
return 0;
}
int Earth::calculateQn_k(int i_k)
{
int i = (i_k > 2) ? i_k : 2;
for (; i < m_cols; i++)
m_Q(i_k, i) = coefficient_С(i_k, i) * m_X(i_k, 0) + coefficient_S(i_k, i) * m_Y(i_k, 0);
return 0;
}
int Earth::calculate_r()
{
m_r = sqrt(pow(m_GDC(0, 0), 2) + pow(m_GDC(1, 0), 2) + pow(m_GDC(2, 0), 2));
return 0;
}
int Earth::calculatePartialQn_k(int i_k)
{
int i = (i_k > 2) ? i_k : 2;
for (; i < m_cols; i++) {
m_Partial_Qx(i_k, i) = coefficient_С(i_k, i) * m_Partial_Xx(i_k, 0) + coefficient_S(i_k, i) * m_Partial_Yx(i_k, 0);
m_Partial_Qy(i_k, i) = coefficient_С(i_k, i) * m_Partial_Xy(i_k, 0) + coefficient_S(i_k, i) * m_Partial_Yy(i_k, 0);
}
return 0;
}
int Earth::calcuclateUkn(int k, int n)
{
m_u.x = m_Partial_R(0, n) * PartialDeriv(0, 0) * m_Z(k, n) * m_Q(k, n)
+ m_R(0, 1) * m_Partial_Z(k, n) * PartialDeriv(3, 0) * m_Q(k, n)
+ m_R(0, 1) * m_Z(k, n) * (m_Partial_Qx(k, n) * PartialDeriv(1, 0) + m_Partial_Qy(k, n) * PartialDeriv(2, 0));
m_u.y = m_Partial_R(0, n) * PartialDeriv(0, 1) * m_Z(k, n) * m_Q(k, n)
+ m_R(0, n) * m_Partial_Z(k, n) * PartialDeriv(3, 1) * m_Q(k, n)
+ m_R(0, n) * m_Z(k, n) * (m_Partial_Qx(k, n) * PartialDeriv(1, 1) + m_Partial_Qy(k, n) * PartialDeriv(2, 1));
m_u.y = m_Partial_R(0, n) * PartialDeriv(0, 2) * m_Z(k, n) * m_Q(k, n)
+ m_R(0, n) * m_Partial_Z(k, n) * PartialDeriv(3, 2) * m_Q(k, n)
+ m_R(0, n) * m_Z(k, n) * (m_Partial_Qx(k, n) * PartialDeriv(1, 2) + m_Partial_Qy(k, n) * PartialDeriv(2, 2));
return 0;
}
int Earth::calculateF(int k) {
double SumX = 0, SumY = 0, SumZ = 0;
for (int i = k; i < m_cols; i++) {
calcuclateUkn(k, i);
SumX += m_u.x;
SumY += m_u.x;
SumZ += m_u.x;
}
m_F(0, 0) += SumX;
m_F(1, 0) += SumY;
m_F(2, 0) += SumZ;
return 0;
}
\ No newline at end of file
Undisturbed_movement/Planets.h
View file @ aa4949cc
...@@ -3,20 +3,21 @@ ...@@ -3,20 +3,21 @@
class Moon : private ConvertSC { class Moon : private ConvertSC {
private: private:
double m_sin_πс; const double gravitational_parameter_Moon{ 4902.799 }; // гравитационный парметр луны
double m_Δλ_c; const double m_scalefactor{ 6378.43817285 };
double m_β;
const double m_scalefactor = 6378.43817285; double m_sin_πс{0};
double m_Δλ_c{0};
double m_β{0};
double m_Rm{ 0 }; // расстояние до луны в текущий момент времени
double m_Ro{ 0 }; // расстояние от луны до объекта в ткущий момент времени
std::vector<double> m_s{ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }; std::vector<double> m_s{ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
std::vector<double> m_r{ 0.0, 0.0, 0.0, 0.0, 0.0 }; std::vector<double> m_r{ 0.0, 0.0, 0.0, 0.0, 0.0 };
fundamental_arguments m_cfa{ 0.0, 0.0, 0.0, 0.0, 0.0 }; fundamental_arguments m_cfa{ 0.0, 0.0, 0.0, 0.0, 0.0 };
public: Matrix m_F{ 3, 1 }; // вектор ускорения вызванного гравитацией луны
Moon(double input_TDB = 0) :ConvertSC{ input_TDB } {
calculateDC();
};
~Moon() {};
int calculate_s(); int calculate_s();
int calculate_r(); int calculate_r();
...@@ -24,7 +25,20 @@ public: ...@@ -24,7 +25,20 @@ public:
int calculate_sin_πс(); int calculate_sin_πс();
int calculate_Δλ_c(); int calculate_Δλ_c();
int calculate_β(); int calculate_β();
int calculateGEC();
int calculateRm();
int calculateRo(const Matrix& DC_coord_object);
public:
Moon(double input_TDB = 0) :ConvertSC{ input_TDB } {
calculateDC();
};
~Moon() {};
int calculateGEC(double input_TDD);
Matrix calculateFm(const Matrix& DC_coord_object, double input_TDB);
friend std::ostream& operator<<(std::ostream& out, Moon& i_moon); friend std::ostream& operator<<(std::ostream& out, Moon& i_moon);
friend std::ofstream& operator<<(std::ofstream& out, Moon& i_moon); friend std::ofstream& operator<<(std::ofstream& out, Moon& i_moon);
...@@ -32,44 +46,138 @@ public: ...@@ -32,44 +46,138 @@ public:
class Sun:private ConvertSC { class Sun:private ConvertSC {
private: private:
const double AU{ 149597870.691 }; const double AU{ 149597870.691 }; // астрономическая единица
std::vector<double> m_M{ 0.0,0.0,0.0,0.0,0.0 }; const double gravitational_parameter_Sun{ 132712438000 }; // гравитационный параметр солнца
const double light_pressure{ 4.5606e-6 };
double m_Rs{0}; // расстояние до солнца в текущий момент времени
double m_Ro{0}; // расстояний от солнца до обхекта в текщий момент времени
std::vector<double> m_M{ 0.0,0.0,0.0,0.0,0.0 };
Matrix m_perturbation{ 3,1 }; // возмущение в эклептической долготе, возмущение в эклептической широте, возмущение в расстоянии до Солнца Matrix m_perturbation{ 3,1 }; // возмущение в эклептической долготе, возмущение в эклептической широте, возмущение в расстоянии до Солнца
Matrix m_F{ 3, 1 }; //вектор ускорения вызванного притяжением солнца
Matrix m_P{ 3, 1 }; // вектор ускрения возванного давлением солнечного свет
//возмущения в долготе обусловленные Венерой //возмущения в долготе обусловленные Венерой
double DILong_V(); double DILong_V();
//возмущение в долготе обусловленные Марсом
double DILong_M(); double DILong_M();
//Возмущения в долготе обусловленные Юпитером и Сатурном
double DILong_JS(); double DILong_JS();
//Возмущения в широте
double DILat(); double DILat();
//Возмущения в расстоянии обусловленные Венрой
double DIR_V(); double DIR_V();
//Возмушения в расстоянии обусловленный Марсом
double DIR_M(); double DIR_M();
//Возмущения в расстоянии обусловленные Юпитером и Свтурном
double DIR_JS(); double DIR_JS();
int calculateRs();
int calculateRo(const Matrix& DC_coord_object);
int calculateAverageAnomalies();
int calculatePerturbation();
public: public:
Sun(double input_TDB = 0) : ConvertSC{ input_TDB } {}; Sun(double input_TDB = 0) : ConvertSC{ input_TDB } {};
~Sun() {}; ~Sun() {};
int calculateGEC(double input_TDB);
Matrix& calculateFs(const Matrix& DC_coord_object, double input_TDB);
Matrix& calculatePs(const Matrix& DC_coord_object, double kr, double As, double ms, double input_TDB);
};
int calculateAverageAnomalies(); struct ratio {
int calculatePerturbation(); double Rr;
int calculateGEC(); double Xr;
double Yr;
double Zr;
};
struct coord {
double x;
double y;
double z;
}; };
class Earth : private ConvertSC { class Earth : private ConvertSC {
private: private:
const double gravitational_parameter_Earth = 398600.4415; const double gravitational_parameter_Earth = 398600.4415; // гравитационный парматер Земли
uint32_t m_cols; const double Radius_Earth = 6378.13630; // средний радиус земли
uint32_t m_rows; double m_r; // радиус до объекта в текущий момент времени
Matrix m_GDC; double m_h; // высота полета
double m_r; ratio m_ratio;
Matrix startPartialDeriv{4,3}; coord m_u;
uint32_t m_cols; //n
uint32_t m_rows; //k
Matrix coefficient_С = fill_matrix("coef_c.txt"); Matrix coefficient_С = fill_matrix("coef_c.txt");
Matrix coefficient_S = fill_matrix("coef_s.txt"); Matrix coefficient_S = fill_matrix("coef_s.txt");
Matrix m_GDC; // прямуоугольные координаты объекта в земной СК
Matrix PartialDeriv{4,3};
Matrix m_R{ 1, m_cols };
Matrix m_Partial_R{ 1, m_cols };
Matrix m_Z{ m_rows, m_cols };
Matrix m_Partial_Z{ m_rows, m_cols };
Matrix m_X{ m_rows, 1 };
Matrix m_Y{ m_rows, 1 };
Matrix m_Partial_Xx{ m_rows, 1 };
Matrix m_Partial_Xy{ m_rows, 1 };
Matrix m_Partial_Yx{ m_rows, 1 };
Matrix m_Partial_Yy{ m_rows, 1 };
Matrix m_Q{ m_rows, m_cols };
Matrix m_Partial_Qx{ m_rows, m_cols };
Matrix m_Partial_Qy{ m_rows, m_cols };
Matrix m_F{3, 1}; // компонеты ускорения вызванного притяжением земли в земных координатах
Matrix m_FCS{ 3, 1 }; // компонеты ускорения вызванного притяжением земли в небесных координатах
int calculate_r();
int calculate_v(const Matrix& GEV);
Matrix fill_matrix(const char fileName []); Matrix fill_matrix(const char fileName []);
int setRatio();
int setF();
int calculateRn();
int calculatePartialRn();
int calculateZn0();
int calculatePartialZn0();
int calculateZn_k(int k);
int calculatePartialZn_k(int k);
int calculateXY_k(int k);
int claculatePartialXxy_k(int k);
int claculatePartialYxy_k(int k);
int calculateQn_k(int k);
int calculatePartialQn_k(int k);
int calcuclateUkn(int k, int n);
int calculateF(int k);
int calculatePartialDeriv();
int calculateStartParameters();
int calculateAccelerationByEarth();
public: public:
Earth() {}; Earth() {};
Earth(double , Matrix); Earth(double input_TDB, Matrix input_GCC) : ConvertSC{ input_TDB}, m_GDC{ CC_to_TC(m_precession, m_nutation, input_GCC) }
{
calculate_r();
calculateAccelerationByEarth();
m_FCS = TC_to_CC(m_F);
}
~Earth() {}; ~Earth() {};
int calculateStartPartialDeriv();
Matrix& recalculateAccelerationCC(double input_TDB, Matrix i_GDC) {
m_GDC = i_GDC;
recalc(input_TDB);
calculateAccelerationByEarth();
Matrix m_CC_to_TC{ CC_to_TC(m_precession, m_nutation, m_GDC) };
m_FCS = m_F * TC_to_CC(m_CC_to_TC);
return m_FCS;
}
}; };
\ No newline at end of file
Undisturbed_movement/Satellite.cpp 0 → 100644
View file @ aa4949cc
#include "Satellite.h"
int Satellite::calculate_r() {
m_r = sqrt(pow(m_DCS(0, 0), 2) + pow(m_DCS(1, 0), 2) + pow(m_DCS(2, 0), 2));
return 0;
}
int Satellite::calculate_density()
{
m_h = m_r - Radius_Earth * (1 - compression_earth*pow((m_DCS(2,0)/m_r), 2));
m_density = 2E-13 * exp(-(m_h - 200) / 60);
return int();
}
int Satellite::calculate_V() {
m_v = sqrt(pow(m_DVS(0, 0), 2) + pow(m_DVS(1, 0), 2) + pow(m_DVS(2, 0), 2));
return 0;
}
Matrix& Satellite::calculate_F(){
calculate_V();
double Sb = m_As / m_ms;
m_F = -Sb * m_density * m_v * m_DVS;
return m_F;
}
Matrix& Satellite::recalculate_F(const Matrix& i_DVS) {
double v = sqrt(pow(i_DVS(0, 0), 2) + pow(i_DVS(1, 0), 2) + pow(i_DVS(2, 0), 2));
double Sb = m_As / m_ms;
m_F = -Sb * m_density * v * i_DVS;
return m_F;
}
Matrix Satellite::step_calculate_F(const Matrix& i_DCS, const Matrix& i_DVS, const double step_time)
{
return (m_moon.calculateFm(i_DCS, step_time) + m_sun.calculateFs(i_DCS, step_time) + m_sun.calculatePs(i_DCS, m_ks, m_As, m_ms, step_time)
+ m_earth.recalculateAccelerationCC(step_time, i_DCS)
+ recalculate_F(i_DVS));
}
int calc_α(std::vector<Matrix>& α,const std::vector<Matrix>& F, const std::vector<double>& step_time, int current_moment) {
switch (current_moment)
{
case 1:
α.at(0) = ((F.at(1) - F.at(0)) / step_time.at(1));
break;
case 2:
α.at(1) = (((F.at(2) - F.at(0)) / step_time.at(2) - α.at(0)) / (step_time.at(2) - step_time.at(1)));
break;
case 3:
α.at(2) = ((((F.at(3) - F.at(0)) / step_time.at(3) - α.at(0)) / (step_time.at(3) - step_time.at(1)) - α.at(1)) /
(step_time.at(3) - step_time.at(2)));
break;
case 4:
α.at(3) = (((((F.at(4) - F.at(0)) / step_time.at(4) - α.at(0)) / (step_time.at(4) - step_time.at(1)) - α.at(1)) / (step_time.at(4) - step_time.at(2)) - α.at(2))
/ (step_time.at(4) - step_time.at(3)));
break;
case 5:
α.at(4) = ((((((F.at(5) - F.at(0)) / (step_time.at(5)) - α.at(0)) / (step_time.at(5) - step_time.at(1)) - α.at(1))
/ (step_time.at(5) - step_time.at(2)) - α.at(2)) / (step_time.at(5) - step_time.at(3)) - α.at(3))
/ (step_time.at(5) - step_time.at(4)));
break;
case 6:
α.at(5) = (((((((F.at(6) - F.at(0)) / step_time.at(6) - α.at(0)) / (step_time.at(6) - step_time.at(1)) - α.at(1))
/ (step_time.at(6) - step_time.at(2)) - α.at(2)) / (step_time.at(6) - step_time.at(3)) - α.at(3))
/ (step_time.at(6) - step_time.at(4)) - α.at(4)) / (step_time.at(6) - step_time.at(5)));
break;
case 7:
α.at(6) = ((((((((F.at(7) - F.at(0)) / step_time.at(7) - α.at(0)) / (step_time.at(7) - step_time.at(1)) - α.at(1))
/ (step_time.at(7) - step_time.at(2)) - α.at(2)) / (step_time.at(7) - step_time.at(3)) - α.at(3))
/ (step_time.at(7) - step_time.at(4)) - α.at(4)) / (step_time.at(7) - step_time.at(5)) - α.at(5))
/ (step_time.at(7) - step_time.at(6)));
break;
}
return 1;
}
int Satellite::integration_step(double i_time) {
std::vector<double> step_time{ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
std::vector<Matrix> F{}; //массив хранящий ускорения по подшагам
std::vector<Matrix> A{}; //массив хранящий коэффициенты полинома
std::vector<Matrix> α{}; //массив хранящий параметры коэффициентов полинома
Matrix buff_DCS{ m_DCS };
Matrix buff_DVS { m_DVS };
Matrix coef_c{8,8};
for (short i = 0; i < 8; i++)
step_time.at(i) = (i_time - m_start_time) * coef_time_step.at(i);
coef_c(0, 0) = 1;
for (short i = 0; i < 8; i++) {
for (short j = 1; j < 8; j++) {
if (i == j)
coef_c(i, j) = 1;
else if (i == 0)
coef_c(j, i) = -step_time.at(j) * coef_c(j - 1, 0);
else
coef_c(i, j) = coef_c(i-1,j-1) - step_time.at(i) * coef_c(i - 1, 0);
}
}
for (int i = 0; i < 8; i++) {
F.push_back(Matrix{ 3, 1 });
A.push_back(Matrix{ 3, 1 });
α.push_back(Matrix{ 3, 1 });
}
for (short i = 0; i < 8; i++) {
buff_DCS = (m_DCS + step_time.at(i) * m_DVS + (1 / 2) * F.at(0) * pow(step_time.at(i), 2)
+ (1 / 6) * A.at(0) * pow(step_time.at(i), 3) + (1 / 12) * A.at(1) * pow(step_time.at(i), 4)
+ (1 / 20) * A.at(2) * pow(step_time.at(i), 5) + (1 / 30) * A.at(3) * pow(step_time.at(i), 6)
+ (1 / 42) * A.at(4) * pow(step_time.at(i), 7));
buff_DVS = (m_DVS + F.at(0) * step_time.at(i) + (1 / 2) * A.at(0) * pow(step_time.at(i), 2)
+ (1 / 3) * A.at(1) * pow(step_time.at(i), 3) + (1 / 4) * A.at(2) * pow(step_time.at(i), 4)
+ (1 / 5) * A.at(3) * pow(step_time.at(i), 5) + (1 / 6) * A.at(4) * pow(step_time.at(i), 5)
+ (1 / 7) * A.at(5) + pow(step_time.at(i), 7));
F.push_back(step_calculate_F(buff_DCS, buff_DVS, step_time.at(i)));
calc_α(α, F, step_time, i);
A.at(0) = α.at(0) + coef_c(1, 0) * α.at(1) + coef_c(2, 0) * α.at(2) + coef_c(3, 0) * α.at(3) + coef_c(4, 0) * α.at(4)
+ coef_c(5, 0) * α.at(5) + coef_c(6, 0) * α.at(6);
A.at(1) = α.at(1) + coef_c(2, 1) * α.at(2) + coef_c(3, 1) * α.at(3) + coef_c(4, 1) * α.at(4) + coef_c(5, 1) * α.at(5)
+ coef_c(6, 1) * α.at(6);
A.at(2) = α.at(2) + coef_c(3, 2) * α.at(3) + coef_c(4, 2) * α.at(4) + coef_c(5, 2) * α.at(5) + coef_c(6, 2) * α.at(6);
A.at(3) = α.at(3) + coef_c(4, 3) * α.at(4) + coef_c(5, 3) * α.at(5) + coef_c(6, 3) * α.at(6);
A.at(4) = α.at(4) + coef_c(5, 4) * α.at(4) + coef_c(6, 4) * α.at(6);
A.at(5) = α.at(5) + coef_c(6, 5) * α.at(6);
}
buff_DCS = (m_DCS + i_time * m_DVS + (1 / 2) * F.at(0) * pow(i_time, 2)
+ (1 / 6) * A.at(0) * pow(i_time, 3) + (1 / 12) * A.at(1) * pow(i_time, 4)
+ (1 / 20) * A.at(2) * pow(i_time, 5) + (1 / 30) * A.at(3) * pow(i_time, 6)
+ (1 / 42) * A.at(4) * pow(i_time, 7));
buff_DVS = (m_DVS + F.at(0) * i_time + (1 / 2) * A.at(0) * pow(i_time, 2)
+ (1 / 3) * A.at(1) * pow(i_time, 3) + (1 / 4) * A.at(2) * pow(i_time, 4)
+ (1 / 5) * A.at(3) * pow(i_time, 5) + (1 / 6) * A.at(4) * pow(i_time, 5)
+ (1 / 7) * A.at(5) + pow(i_time, 7));
m_DCS = buff_DCS;
m_DVS = buff_DVS;
std::cout << "buff_DCS = " << buff_DCS << " buff_DVS = " << buff_DVS;
return 0;
}
Undisturbed_movement/Satellite.h 0 → 100644
View file @ aa4949cc
#pragma once
#include "Planets.h"
class Satellite {
private:
const double Radius_Earth = 6378.13630; //
const double compression_earth = 1 / 298.257; //
const std::vector<double> coef_time_step{ 0.000000000000000000, 0.056262560526922147, //
0.180240691736892365, 0.352624717113169637, //
0.547153626330555383, 0.734210177215410532,
0.885320946839095768, 0.977520613561287501 };
double m_start_time;
double m_As = 0; //
double m_ms = 1; //
double m_ks = 1.5; //
double m_r; // -
double m_v; //
double m_h; //
double m_density; //
Matrix m_DCS{ 3, 1 }; //
Matrix m_DVS{ 3, 1 }; //
Matrix m_F{ 3, 1 }; //
Moon m_moon{m_start_time};
Sun m_sun{m_start_time};
Earth m_earth{ m_start_time, m_DCS };
int calculate_r();
int calculate_density();
int calculate_V();
Matrix& calculate_F();
Matrix& recalculate_F(const Matrix& i_DVS);
Matrix step_calculate_F(const Matrix& i_DCS, const Matrix& i_DVS, const double step_time);
public:
Satellite(double time_MJD, Matrix& i_DCS, Matrix& i_DVS) : m_start_time{time_MJD}, m_DCS { i_DCS }, m_DVS{ i_DVS } {};
int integration_step(double i_time);
~Satellite() {};
};
\ No newline at end of file
Undisturbed_movement/SystemCoords.cpp
View file @ aa4949cc
...@@ -74,6 +74,13 @@ int ConvertSC::rewrite() { ...@@ -74,6 +74,13 @@ int ConvertSC::rewrite() {
return 0; return 0;
} }
int ConvertSC::recalc(double input_TDB)
{
settime(input_TDB);
calculateDC();
return 0;
}
int ConvertSC::calculateDC() int ConvertSC::calculateDC()
{ {
Matrix buff1 = AEC_to_CC(m_precession); Matrix buff1 = AEC_to_CC(m_precession);
...@@ -87,50 +94,50 @@ int ConvertSC::calculateDC() ...@@ -87,50 +94,50 @@ int ConvertSC::calculateDC()
Matrix Earth_rot_matrix(double GTST) { Matrix Earth_rot_matrix(double GTST) {
return Matrix::RotationOZ(GTST); return Matrix::RotationOZ(GTST);
} }
//
std::ostream& operator<<(std::ostream& out, precession_parameters& input_pp) //std::ostream& operator<<(std::ostream& out, precession_parameters& input_pp)
{ //{
out << "-------------------------------------------------------------" // out << "-------------------------------------------------------------"
<< "\nPrecession parameters\n" << "ζα = " << input_pp.m_ζα // << "\nPrecession parameters\n" << "ζα = " << input_pp.m_ζα
<< "\nθα = " << input_pp.m_θα << "\nm_zα = " << input_pp.m_zα // << "\nθα = " << input_pp.m_θα << "\nm_zα = " << input_pp.m_zα
<< "\n-------------------------------------------------------------"; // << "\n-------------------------------------------------------------";
return out; // return out;
} //}
//
std::ofstream& operator<<(std::ofstream& out, precession_parameters& input_pp) //std::ofstream& operator<<(std::ofstream& out, precession_parameters& input_pp)
{ //{
out << input_pp; // out << input_pp;
return out; // return out;
} //}
//
std::ostream& operator<<(std::ostream& out, fundamental_arguments& input_fa) //std::ostream& operator<<(std::ostream& out, fundamental_arguments& input_fa)
{ //{
out << "\n-------------------------------------------------------------" // out << "\n-------------------------------------------------------------"
<< "\nfundamental arguments" // << "\nfundamental arguments"
<< "\nλ = " << input_fa.λ << "\nl = " << input_fa.l // << "\nλ = " << input_fa.λ << "\nl = " << input_fa.l
<< "\nL = " << input_fa.L << "\nF = " << input_fa.F // << "\nL = " << input_fa.L << "\nF = " << input_fa.F
<< "\nD = " << input_fa.D // << "\nD = " << input_fa.D
<< "\n-------------------------------------------------------------"; // << "\n-------------------------------------------------------------";
return out; // return out;
} //}
//
std::ofstream& operator<<(std::ofstream& out, fundamental_arguments& input_fa) //std::ofstream& operator<<(std::ofstream& out, fundamental_arguments& input_fa)
{ //{
out << input_fa; // out << input_fa;
return out; // return out;
} //}
//
std::ostream& operator<<(std::ostream& out, ConvertSC& input_CSC) //std::ostream& operator<<(std::ostream& out, ConvertSC& input_CSC)
{ //{
out << "\n*************************************************************" // out << "\n*************************************************************"
<< "\nFields of the ConvertCS class" // << "\nFields of the ConvertCS class"
<< "\nm_ε = " << input_CSC.m_ε << "\nΔψ = " << input_CSC.Δψ // << "\nm_ε = " << input_CSC.m_ε << "\nΔψ = " << input_CSC.Δψ
<< "\nm_Δε = " << input_CSC.Δε << "\nm_t = " << input_CSC.m_t // << "\nm_Δε = " << input_CSC.Δε << "\nm_t = " << input_CSC.m_t
<< "\nm_ts = " << input_CSC.m_ts << "\nprecession parameters = " << input_CSC.m_pp // << "\nm_ts = " << input_CSC.m_ts << "\nprecession parameters = " << input_CSC.m_pp
<< "\nprecession matrix = " << input_CSC.m_precession << "\nnutation matrix = " << input_CSC.m_nutation // << "\nprecession matrix = " << input_CSC.m_precession << "\nnutation matrix = " << input_CSC.m_nutation
<< "\nfundamental arguments = " << input_CSC.m_fa; // << "\nfundamental arguments = " << input_CSC.m_fa;
return out; // return out;
} //}
std::ofstream& operator<<(std::ofstream& out, ConvertSC& input_CSC) std::ofstream& operator<<(std::ofstream& out, ConvertSC& input_CSC)
{ {
...@@ -155,7 +162,7 @@ Matrix TEC_to_CC(Matrix& TEC) ...@@ -155,7 +162,7 @@ Matrix TEC_to_CC(Matrix& TEC)
Matrix CC_to_TC(Matrix& Precession, Matrix& Nutation, Matrix& ERM) Matrix CC_to_TC(Matrix& Precession, Matrix& Nutation, Matrix& ERM)
{ {
Matrix buf = CC_to_TEC(Precession, Nutation); Matrix buf = CC_to_TEC(Precession, Nutation);
return (ERM * buf); return (buf * ERM);
} }
Matrix TC_to_CC(Matrix& CC_to_TC) Matrix TC_to_CC(Matrix& CC_to_TC)
......
Undisturbed_movement/SystemCoords.h
View file @ aa4949cc
...@@ -31,8 +31,8 @@ protected: ...@@ -31,8 +31,8 @@ protected:
double m_t{ 0 }, m_ts{ 0 }; double m_t{ 0 }, m_ts{ 0 };
Matrix m_precession{ 3,3 }; Matrix m_precession{ 3,3 };
Matrix m_nutation{ 3,3 }; Matrix m_nutation{ 3,3 };
Matrix m_GEC{ 3, 1 }; Matrix m_GEC{ 3, 1 }; //геоцентрические эклептические координаты
Matrix m_DC{ 3, 1 }; Matrix m_DC{ 3, 1 }; //прямоугольные координаты в небесной СК
precession_parameters m_pp{ 0,0,0 }; precession_parameters m_pp{ 0,0,0 };
fundamental_arguments m_fa{ 0,0,0,0,0 }; fundamental_arguments m_fa{ 0,0,0,0,0 };
...@@ -59,6 +59,7 @@ public: ...@@ -59,6 +59,7 @@ public:
int calculateNutationParameters(); int calculateNutationParameters();
int calculateMatrixNutation(); int calculateMatrixNutation();
int rewrite(); int rewrite();
int recalc(double input_TDB = 0);
double getΔψ() { double getΔψ() {
return Δψ; return Δψ;
...@@ -66,15 +67,10 @@ public: ...@@ -66,15 +67,10 @@ public:
double getε() { double getε() {
return m_ε; return m_ε;
} }
int calculateGEC() {};
int calculateDC(); int calculateDC();
int recalc(double input_time = 0) {
settime(input_time);
calculateGEC();
calculateDC();
return 0;
}
friend std::ostream& operator<<(std::ostream& out, ConvertSC& input_CSC); friend std::ostream& operator<<(std::ostream& out, ConvertSC& input_CSC);
friend std::ofstream& operator<<(std::ofstream& out, ConvertSC& input_CSC); friend std::ofstream& operator<<(std::ofstream& out, ConvertSC& input_CSC);
......
Undisturbed_movement/Undisturbed_movement.cpp
View file @ aa4949cc
#include "Planets.h" #include "Satellite.h"
int main() int main()
{ {
Matrix a{ 3, 1, 5.6 }; Matrix a{ 3, 1, {1000, 1000, 1000} };
Earth m_earth{18657, a}; Matrix b{ 3, 1, {100, 100, 100} };
Satellite m_s{ 59777.67019, a, b };
m_s.integration_step(59777.77019);
return 0; return 0;
} }
......
Undisturbed_movement/Undisturbed_movement.vcxproj
View file @ aa4949cc
...@@ -145,6 +145,7 @@ ...@@ -145,6 +145,7 @@
<ClCompile Include="orbit_time.cpp" /> <ClCompile Include="orbit_time.cpp" />
<ClCompile Include="orbit.cpp" /> <ClCompile Include="orbit.cpp" />
<ClCompile Include="Planets.cpp" /> <ClCompile Include="Planets.cpp" />
<ClCompile Include="Satellite.cpp" />
<ClCompile Include="SystemCoords.cpp" /> <ClCompile Include="SystemCoords.cpp" />
<ClCompile Include="Undisturbed_movement.cpp" /> <ClCompile Include="Undisturbed_movement.cpp" />
</ItemGroup> </ItemGroup>
...@@ -153,6 +154,7 @@ ...@@ -153,6 +154,7 @@
<ClInclude Include="orbit.h" /> <ClInclude Include="orbit.h" />
<ClInclude Include="orbit_time.h" /> <ClInclude Include="orbit_time.h" />
<ClInclude Include="Planets.h" /> <ClInclude Include="Planets.h" />
<ClInclude Include="Satellite.h" />
<ClInclude Include="SystemCoords.h" /> <ClInclude Include="SystemCoords.h" />
</ItemGroup> </ItemGroup>
<Import Project="$(VCTargetsPath)\Microsoft.Cpp.targets" /> <Import Project="$(VCTargetsPath)\Microsoft.Cpp.targets" />
......
Undisturbed_movement/Undisturbed_movement.vcxproj.filters
View file @ aa4949cc
...@@ -33,6 +33,9 @@ ...@@ -33,6 +33,9 @@
<ClCompile Include="Planets.cpp"> <ClCompile Include="Planets.cpp">
<Filter>Исходные файлы</Filter> <Filter>Исходные файлы</Filter>
</ClCompile> </ClCompile>
<ClCompile Include="Satellite.cpp">
<Filter>Исходные файлы</Filter>
</ClCompile>
</ItemGroup> </ItemGroup>
<ItemGroup> <ItemGroup>
<ClInclude Include="orbit.h"> <ClInclude Include="orbit.h">
...@@ -50,5 +53,8 @@ ...@@ -50,5 +53,8 @@
<ClInclude Include="Planets.h"> <ClInclude Include="Planets.h">
<Filter>Файлы заголовков</Filter> <Filter>Файлы заголовков</Filter>
</ClInclude> </ClInclude>
<ClInclude Include="Satellite.h">
<Filter>Файлы заголовков</Filter>
</ClInclude>
</ItemGroup> </ItemGroup>
</Project> </Project>
\ No newline at end of file
Undisturbed_movement/coef_s.txt
View file @ aa4949cc
12 11 13 11
0 0 0 0 0 0 0 0 0 0 0
1.54309997E-9 2.68011894E-7 -4.49459935E-7 -8.06634638E-8 2.11646643E-8 6.93698935E-8 4.01997816E-8 1.42365696E-8 -8.12891488E-8 -1.64654645E-8 -2.37844845E-8 1.54309997E-9 2.68011894E-7 -4.49459935E-7 -8.06634638E-8 2.11646643E-8 6.93698935E-8 4.01997816E-8 1.42365696E-8 -8.12891488E-8 -1.64654645E-8 -2.37844845E-8
-9.03868073E-7 -2.11402398E-7 1.48155457E-7 -5.23267240E-8 -4.65039481E-8 9.28231439E-9 5.38131641E-9 -2.22867920E-9 -3.05712916E-9 -5.09736032E-9 1.41642228E-9 -9.03868073E-7 -2.11402398E-7 1.48155457E-7 -5.23267240E-8 -4.65039481E-8 9.28231439E-9 5.38131641E-9 -2.22867920E-9 -3.05712916E-9 -5.09736032E-9 1.41642228E-9
0 1.97201324E-7 -1.20112918E-8 -7.10087714E-9 1.84313369E-10 -3.06115024E-9 -8.72351950E-10 -5.63392145E-10 -8.98933286E-10 -6.86352078E-10 9.15457482E-11 0 1.97201324E-7 -1.20112918E-8 -7.10087714E-9 1.84313369E-10 -3.06115024E-9 -8.72351950E-10 -5.63392145E-10 -8.98933286E-10 -6.86352078E-10 9.15457482E-11
......