Undisturbed_movement/Matrix.cpp

 #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);
-	for (int i = 0; i < output.m_rows; i++) {
-		for (int j = 0; j < output.m_cols; j++)
+	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+(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 };
 	return (matrix_1+buff);
 }
 
-Matrix operator+(double in, Matrix& matrix_1){
+Matrix operator+(const double in, const Matrix& matrix_1){
 	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));
 	Matrix output{ matrix_1.m_rows, matrix_2.m_cols };
-	for (int i = 0; i < output.m_rows; i++) {
-		for (int j = 0; j < output.m_cols; j++) {
+	for (uint32_t i = 0; i < output.m_rows; i++) {
+		for (uint32_t j = 0; j < output.m_cols; j++) {
 			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);
 			output.at(i, j) = result;
 		}
@@ -33,20 +54,28 @@ Matrix operator*(Matrix& matrix_1, Matrix& matrix_2) {
 	return output;
 }
 
-Matrix operator*(Matrix& matrix_1, double in) {
-	Matrix output{ matrix_1.m_cols, matrix_1.m_rows };
-	for (int i = 0; i < output.m_rows; i++) {
-		for (int j = 0; j < output.m_cols; j++) {
+Matrix operator*(const Matrix& matrix_1, const double in) {
+	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) * in;
 		}
 	}
 	return output;
 }
 
-Matrix operator*(double in, Matrix& matrix_1) {
+Matrix operator*(const double in, const Matrix& matrix_1) {
 		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)
 {
 	if (this == &input_matrix) return *this;
@@ -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";
-	for (int i = 0; i < matrix_1.m_rows; i++) {
-		for (int j = 0; j < matrix_1.m_cols; j++) {
+	for (uint32_t i = 0; i < matrix_1.m_rows; i++) {
+		for (uint32_t j = 0; j < matrix_1.m_cols; j++) {
 			out << matrix_1(i, j) << " ";
 		}
 		out << "\n";
@@ -69,19 +98,24 @@ std::ostream& operator<<(std::ostream& out, Matrix& matrix_1) {
 	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";
-	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;
+	return (at(input_rows, input_cols));
 }
 
-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));
 }
@@ -90,8 +124,8 @@ double& Matrix::operator()(uint32_t input_rows, uint32_t input_cols)
 Matrix Matrix::Transpose()
 {
 	Matrix output{m_cols, m_rows };
-	for (int i = 0; i < output.m_rows; i++) {
-		for (int j = 0; j < output.m_cols; j++)
+	for (uint32_t i = 0; i < output.m_rows; i++) {
+		for (uint32_t j = 0; j < output.m_cols; j++)
 			output.at(i, j) = at(j, i);
 	}
 	return output;

Undisturbed_movement/Matrix.h

@@ -27,20 +27,35 @@ public:
 		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) {
 		m_matrix.push_back(input_number);
 		return 0;
 	}
 
 	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+(Matrix& matrix_1, double in);
-	friend Matrix operator+(double in, Matrix& matrix_1);
+	friend Matrix operator/(const Matrix& matrix_1, const double in);
+	friend Matrix operator/(const double in, const 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 Transpose();
@@ -48,8 +63,8 @@ public:
 	static Matrix RotationOY(double α);
 	static Matrix RotationOZ(double α);
 
-	friend std::ostream& operator<<(std::ostream& out, Matrix& i_matrix);
-	friend std::ofstream& operator<<(std::ofstream& out, Matrix& i_matrix);
+	friend std::ostream& operator<<(std::ostream& out, const Matrix& i_matrix);
+	friend std::ofstream& operator<<(std::ofstream& out, const Matrix& i_matrix);
 };
 
 

Undisturbed_movement/Planets.cpp

@@ -88,8 +88,9 @@ int Moon::calculate_β()
 	return 0;
 }
 
-int Moon::calculateGEC()
+int Moon::calculateGEC(double input_TDB)
 {
+	recalc(input_TDB);
 	calculateСorrecteFundamentalArg();
 	calculate_Δλ_c();
 	calculate_β();
@@ -97,6 +98,7 @@ int Moon::calculateGEC()
 	m_GEC(0, 0) = m_cfa.λ + RS * m_Δλ_c;
 	m_GEC(1, 0) = RS * m_β;
 	m_GEC(2, 0) = m_scalefactor * (1 / (RS * m_sin_πс));
+	calculateDC();
 	return 0;
 }
 
@@ -110,36 +112,61 @@ int Moon::calculateСorrecteFundamentalArg() {
 	return 0;
 }
 
-std::ostream& operator<<(std::ostream& out, Moon& i_moon)
+int Moon::calculateRm()
 {
-	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;
+	m_Rm = sqrt(pow(m_DC(0, 0), 2) + pow(m_DC(1, 0), 2) + pow(m_DC(2, 0), 2));
+	return 0;
+}
+
+int Moon::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;
 }
 
-std::ofstream& operator<<(std::ofstream& out, Moon& i_moon)
+Matrix Moon::calculateFm(const Matrix& DC_coord_object, double input_TDB)
 {
-	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;
+	calculateGEC(input_TDB);
+	calculateRm();
+	calculateRo(DC_coord_object);
+	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)));
+	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)));
+	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)));
+	return m_F;
 }
 
+//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()
 {
 	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() {
 	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()
 {
 	m_M.at(0) = 0.87167007 + 1021.32292307 * m_ts;
@@ -344,13 +409,15 @@ int Sun::calculatePerturbation()
 	return 0;
 }
 
-int Sun::calculateGEC()
+int Sun::calculateGEC(double input_TDB)
 {
+	recalc(input_TDB);
 	calculateAverageAnomalies();
 	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(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));
+	calculateDC();
 	return 0;
 }
 
@@ -372,11 +439,15 @@ Matrix Earth::fill_matrix(const char fileName[])
 	
 	str >> buff1;
 	rows = stoi(buff1);
+	
 	str >> 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; 
-	cols = 0;
+	cols = 2;
 
 	while (!fin.eof()) {
 		std::getline(fin, buff);
@@ -386,37 +457,202 @@ Matrix Earth::fill_matrix(const char fileName[])
 			output(rows, cols) = stod(buff1);
 			cols++;
 		}
-		cols = 0;
+		cols = 2;
 		rows++;
-		
 	}
 
 	fin.close();
 	return output;
 }
 
-Earth::Earth(double time_MJD, Matrix input_GCC) : ConvertSC{ time_MJD }, m_GDC{ CC_to_TC(m_precession, m_nutation, input_GCC) } {
-	r = sqrt(pow(m_GDC(0,0), 2) + pow(m_GDC(1,0), 2) + pow(m_GDC(2,0), 2));
+int Earth::setF()
+{
+	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
-	startPartialDeriv(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);
-	startPartialDeriv(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(0, 0) = -(m_GDC(0, 0)) / pow(m_r, 3);
+	PartialDeriv(1, 0) = (1 / m_r) - pow(m_GDC(0, 0), 2) / pow(m_r, 3);
+	PartialDeriv(2, 0) = -(m_GDC(0, 0) * m_GDC(1, 0)) / pow(m_r, 3);
+	PartialDeriv(3, 0) = -(m_GDC(0, 0) * m_GDC(2, 0)) / pow(m_r, 3);
 
 	//вычисление частныхм производных по y
-	startPartialDeriv(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);
-	startPartialDeriv(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(0, 1) = -(m_GDC(1, 0)) / pow(m_r, 3);
+	PartialDeriv(1, 1) = -(m_GDC(0, 0) * m_GDC(1, 0)) / pow(m_r, 3);
+	PartialDeriv(2, 1) = (1 / m_r) - pow(m_GDC(1, 0), 2) / pow(m_r, 3);
+	PartialDeriv(3, 1) = -(m_GDC(1, 0) * m_GDC(2, 0)) / pow(m_r, 3);
 
 	//вычисление частных производных по z
-	startPartialDeriv(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);
-	startPartialDeriv(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(0, 1) = -(m_GDC(2, 0)) / pow(m_r, 3);
+	PartialDeriv(1, 1) = -(m_GDC(0, 0) * m_GDC(2, 0)) / pow(m_r, 3);
+	PartialDeriv(2, 1) = -(m_GDC(1, 0) * m_GDC(2, 0)) / 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;
 }
+
+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

@@ -3,20 +3,21 @@
 
 class Moon : private ConvertSC {
 private:
-	double m_sin_πс;
-	double m_Δλ_c;
-	double m_β;
-	const double m_scalefactor = 6378.43817285;
+	const double gravitational_parameter_Moon{ 4902.799 }; // гравитационный парметр луны
+	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_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 };
-	
 
-public:
-	Moon(double input_TDB = 0) :ConvertSC{ input_TDB } {
-		calculateDC();
-	};
-	~Moon() {};
+	Matrix m_F{ 3, 1 }; // вектор ускорения вызванного гравитацией луны
 
 	int calculate_s();
 	int calculate_r();
@@ -24,7 +25,20 @@ public:
 	int calculate_sin_πс();
 	int calculate_Δλ_c();
 	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::ofstream& operator<<(std::ofstream& out, Moon& i_moon);
@@ -32,44 +46,138 @@ public:
 
 class Sun:private ConvertSC {
 private:
-	const double AU{ 149597870.691 };
-	std::vector<double> m_M{ 0.0,0.0,0.0,0.0,0.0 };
+	const double AU{ 149597870.691 }; // астрономическая единица
+	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_F{ 3, 1 }; //вектор ускорения вызванного притяжением солнца
+	Matrix m_P{ 3, 1 }; // вектор ускрения возванного давлением солнечного свет 
+
 	//возмущения в долготе обусловленные Венерой
 	double DILong_V();
+	//возмущение в долготе обусловленные Марсом
 	double DILong_M();
+	//Возмущения в долготе обусловленные Юпитером и Сатурном
 	double DILong_JS();
+	//Возмущения в широте
 	double DILat();
+	//Возмущения в расстоянии обусловленные Венрой 
 	double DIR_V();
+	//Возмушения в расстоянии обусловленный Марсом
 	double DIR_M();
+	//Возмущения в расстоянии обусловленные Юпитером и Свтурном
 	double DIR_JS();
 
+	int calculateRs();
+	int calculateRo(const Matrix& DC_coord_object);
+	int calculateAverageAnomalies();
+	int calculatePerturbation();
+
 public:
 	Sun(double input_TDB = 0) : ConvertSC{ input_TDB } {};
 	~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();
-	int calculatePerturbation();
-	int calculateGEC();
+struct ratio {
+	double Rr;
+	double Xr;
+	double Yr;
+	double Zr;
+};
+struct coord {
+	double x;
+	double y;
+	double z;
 };
 
 class Earth : private ConvertSC {
 private:
-	const double  gravitational_parameter_Earth = 398600.4415;
-	uint32_t m_cols;
-	uint32_t m_rows;
-	Matrix m_GDC;
-	double m_r;
-	Matrix startPartialDeriv{4,3}; 
-
+	const double  gravitational_parameter_Earth = 398600.4415; // гравитационный парматер Земли
+	const double Radius_Earth = 6378.13630; // средний радиус земли
+	double m_r; // радиус до объекта в текущий момент времени 
+	double m_h; // высота полета
+	ratio m_ratio; 
+	coord m_u;
+	uint32_t m_cols; //n
+	uint32_t m_rows; //k
 	
 	Matrix coefficient_С = fill_matrix("coef_c.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 []);
+	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:
 	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() {};
-	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

+#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

+#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

@@ -74,6 +74,13 @@ int ConvertSC::rewrite() {
 	return 0;
 }
 
+int ConvertSC::recalc(double input_TDB)
+{
+	settime(input_TDB);
+	calculateDC();
+	return 0;
+}
+
 int ConvertSC::calculateDC()
 {
 	Matrix buff1 = AEC_to_CC(m_precession);
@@ -87,50 +94,50 @@ int ConvertSC::calculateDC()
 Matrix Earth_rot_matrix(double GTST) {
 	return Matrix::RotationOZ(GTST);
 }
-
-std::ostream& operator<<(std::ostream& out, precession_parameters& input_pp)
-{
-	out << "-------------------------------------------------------------"
-		<< "\nPrecession parameters\n" << "ζα = " << input_pp.m_ζα
-		<< "\nθα = " << input_pp.m_θα << "\nm_zα = " << input_pp.m_zα
-		<< "\n-------------------------------------------------------------";
-	return out;
-}
-
-std::ofstream& operator<<(std::ofstream& out, precession_parameters& input_pp)
-{
-	out << input_pp;
-	return out;
-}
-
-std::ostream& operator<<(std::ostream& out, fundamental_arguments& input_fa)
-{
-	out << "\n-------------------------------------------------------------"
-		<< "\nfundamental arguments"
-		<< "\nλ = " << input_fa.λ << "\nl = " << input_fa.l
-		<< "\nL = " << input_fa.L << "\nF = " << input_fa.F
-		<< "\nD = " << input_fa.D
-		<< "\n-------------------------------------------------------------";
-	return out;
-}
-
-std::ofstream& operator<<(std::ofstream& out, fundamental_arguments& input_fa)
-{
-	out << input_fa;
-	return out;
-}
-
-std::ostream& operator<<(std::ostream& out, ConvertSC& input_CSC)
-{
-	out << "\n*************************************************************"
-		<< "\nFields of the ConvertCS class"
-		<< "\nm_ε = " << input_CSC.m_ε << "\nΔψ = " << input_CSC.Δψ
-		<< "\nm_Δε = " << input_CSC.Δε << "\nm_t = " << input_CSC.m_t
-		<< "\nm_ts = " << input_CSC.m_ts << "\nprecession parameters = " << input_CSC.m_pp
-		<< "\nprecession matrix = " << input_CSC.m_precession << "\nnutation matrix = " << input_CSC.m_nutation
-		<< "\nfundamental arguments = " << input_CSC.m_fa;
-	return out;
-}
+//
+//std::ostream& operator<<(std::ostream& out, precession_parameters& input_pp)
+//{
+//	out << "-------------------------------------------------------------"
+//		<< "\nPrecession parameters\n" << "ζα = " << input_pp.m_ζα
+//		<< "\nθα = " << input_pp.m_θα << "\nm_zα = " << input_pp.m_zα
+//		<< "\n-------------------------------------------------------------";
+//	return out;
+//}
+//
+//std::ofstream& operator<<(std::ofstream& out, precession_parameters& input_pp)
+//{
+//	out << input_pp;
+//	return out;
+//}
+//
+//std::ostream& operator<<(std::ostream& out, fundamental_arguments& input_fa)
+//{
+//	out << "\n-------------------------------------------------------------"
+//		<< "\nfundamental arguments"
+//		<< "\nλ = " << input_fa.λ << "\nl = " << input_fa.l
+//		<< "\nL = " << input_fa.L << "\nF = " << input_fa.F
+//		<< "\nD = " << input_fa.D
+//		<< "\n-------------------------------------------------------------";
+//	return out;
+//}
+//
+//std::ofstream& operator<<(std::ofstream& out, fundamental_arguments& input_fa)
+//{
+//	out << input_fa;
+//	return out;
+//}
+//
+//std::ostream& operator<<(std::ostream& out, ConvertSC& input_CSC)
+//{
+//	out << "\n*************************************************************"
+//		<< "\nFields of the ConvertCS class"
+//		<< "\nm_ε = " << input_CSC.m_ε << "\nΔψ = " << input_CSC.Δψ
+//		<< "\nm_Δε = " << input_CSC.Δε << "\nm_t = " << input_CSC.m_t
+//		<< "\nm_ts = " << input_CSC.m_ts << "\nprecession parameters = " << input_CSC.m_pp
+//		<< "\nprecession matrix = " << input_CSC.m_precession << "\nnutation matrix = " << input_CSC.m_nutation
+//		<< "\nfundamental arguments = " << input_CSC.m_fa;
+//	return out;
+//}
 
 std::ofstream& operator<<(std::ofstream& out, ConvertSC& input_CSC)
 {
@@ -155,7 +162,7 @@ Matrix TEC_to_CC(Matrix& TEC)
 Matrix CC_to_TC(Matrix& Precession, Matrix& Nutation, Matrix& ERM)
 {
 	Matrix buf = CC_to_TEC(Precession, Nutation);
-	return (ERM * buf);
+	return (buf * ERM);
 }
 
 Matrix TC_to_CC(Matrix& CC_to_TC)

Undisturbed_movement/SystemCoords.h

@@ -31,8 +31,8 @@ protected:
 	double m_t{ 0 }, m_ts{ 0 };
 	Matrix m_precession{ 3,3 };
 	Matrix m_nutation{ 3,3 };
-	Matrix m_GEC{ 3, 1 };
-	Matrix m_DC{ 3, 1 };
+	Matrix m_GEC{ 3, 1 }; //геоцентрические эклептические координаты
+	Matrix m_DC{ 3, 1 }; //прямоугольные координаты в небесной СК
 	precession_parameters m_pp{ 0,0,0 };
 	fundamental_arguments m_fa{ 0,0,0,0,0 };
 
@@ -59,6 +59,7 @@ public:
 	int calculateNutationParameters();
 	int calculateMatrixNutation();
 	int rewrite();
+	int recalc(double input_TDB = 0);
 
 	double getΔψ() {
 		return Δψ;
@@ -66,15 +67,10 @@ public:
 	double getε() {
 		return m_ε;
 	}
-	int calculateGEC() {};
+	
 	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::ofstream& operator<<(std::ofstream& out, ConvertSC& input_CSC);

Undisturbed_movement/Undisturbed_movement.cpp

-#include "Planets.h"
+#include "Satellite.h"
 int main()
 {
-    Matrix a{ 3, 1, 5.6 };
-    Earth m_earth{18657, a};
+    Matrix a{ 3, 1, {1000, 1000, 1000} };
+    Matrix b{ 3, 1, {100, 100, 100} };
+    Satellite m_s{ 59777.67019, a, b };
+    m_s.integration_step(59777.77019);
     return 0;
 }
 

Undisturbed_movement/Undisturbed_movement.vcxproj

@@ -145,6 +145,7 @@
     <ClCompile Include="orbit_time.cpp" />
     <ClCompile Include="orbit.cpp" />
     <ClCompile Include="Planets.cpp" />
+    <ClCompile Include="Satellite.cpp" />
     <ClCompile Include="SystemCoords.cpp" />
     <ClCompile Include="Undisturbed_movement.cpp" />
   </ItemGroup>
@@ -153,6 +154,7 @@
     <ClInclude Include="orbit.h" />
     <ClInclude Include="orbit_time.h" />
     <ClInclude Include="Planets.h" />
+    <ClInclude Include="Satellite.h" />
     <ClInclude Include="SystemCoords.h" />
   </ItemGroup>
   <Import Project="$(VCTargetsPath)\Microsoft.Cpp.targets" />

Undisturbed_movement/Undisturbed_movement.vcxproj.filters

@@ -33,6 +33,9 @@
     <ClCompile Include="Planets.cpp">
       <Filter>Исходные файлы</Filter>
     </ClCompile>
+    <ClCompile Include="Satellite.cpp">
+      <Filter>Исходные файлы</Filter>
+    </ClCompile>
   </ItemGroup>
   <ItemGroup>
     <ClInclude Include="orbit.h">
@@ -50,5 +53,8 @@
     <ClInclude Include="Planets.h">
       <Filter>Файлы заголовков</Filter>
     </ClInclude>
+    <ClInclude Include="Satellite.h">
+      <Filter>Файлы заголовков</Filter>
+    </ClInclude>
   </ItemGroup>
 </Project>
\ No newline at end of file

Undisturbed_movement/coef_s.txt

-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
 -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