Satellite.cpp

#include "Satellite.h"

int Satellite::satellite_number = 0;
const std::vector<double> Satellite::coef_time_step { 0.000000000000000000, 0.056262560526922147, //коэффиценты для разбития 
							     					  0.180240691736892365, 0.352624717113169637,  // шага на подинтервалы
                                                      0.547153626330555383, 0.734210177215410532,
										              0.885320946839095768, 0.977520613561287501, 1 };

double calculate_E(const double M, const double e) {
	double buffer_Ek{ M };
	double buffer_Ek1{ 0 };
	for (short i = 0; i < 4; i++) {
		//buffer_Ek1 = M + e * sin(buffer_Ek);
		buffer_Ek1 = buffer_Ek - ((buffer_Ek - e*sin(buffer_Ek) - M) / (1 - e * cos(buffer_Ek)));
		buffer_Ek = buffer_Ek1;
	}
	return buffer_Ek1;
}

int Satellite::calculate_r(const Matrix& DCS) {

	m_r = sqrt(pow(DCS(0, 0), 2) + pow(DCS(1, 0), 2) + pow(DCS(2, 0), 2));
	return 0;
}

int Satellite::calculate_density(const Matrix& DCS)
{
	calculate_r(DCS);
	m_h = m_r - Radius_Earth * (1 - compression_earth*pow((DCS(2,0)/m_r), 2));
	m_density = 2E-13 * exp((-m_h + 200) / 60);
	return 0;
}

int Satellite::calculate_V(const Matrix& DVS) {

	m_v = sqrt(pow(DVS(0, 0), 2) + pow(DVS(1, 0), 2) + pow(DVS(2, 0), 2));
	return 0;
}

Matrix& Satellite::calculate_F(const Matrix& i_DCS, const Matrix& i_DVS) {
	calculate_V(i_DVS);
	calculate_density(i_DCS);
	double Sb = m_param.cross_sectional_radius / m_param.mass;
	m_F = -Sb * m_density * m_v * i_DVS;
	return m_F;
}

void Satellite::step_calculate_F(const Matrix& i_DCS, const Matrix& i_DVS, std::vector<Matrix>& F, const double step_time, const short step_integration)
{
	F.at(step_integration) = m_moon->calculateFm(i_DCS, step_time) + m_sun->calculateFs(i_DCS, step_time) +
		m_sun->calculatePs(i_DCS, m_param, step_time) + m_earth->recalculateAccelerationCC(step_time, i_DCS) +
		calculate_F(i_DCS, i_DVS);
	//F.at(step_integration) = i_DCS + step_time ;
}



void calc_α(std::vector<Matrix>& α,const std::vector<Matrix>& F,  const std::vector<double>& step_time, int current_moment) {
	switch (current_moment)
	{
	case 1:
		α[0] = (F[1] - F[0]) / step_time[1];
		break;
	case 2:
		α[1] = (((F[2] - F[0]) / step_time[2]) - α[0]) / (step_time[2] - step_time[1]);
		break;
	case 3:
		α[2] = ((((F[3] - F[0]) / step_time[3]) - α[0]) / 
			(step_time[3] - step_time[1]) - α[1]) /
			(step_time[3] - step_time[2]);