【C++】利用四阶龙格库塔轨道外推

最新推荐文章于 2024-05-16 11:10:23 发布

叶洋不太冷

最新推荐文章于 2024-05-16 11:10:23 发布

阅读量466

点赞数 1

分类专栏： C++卫星轨道学习文章标签： c++

本文链接：https://blog.csdn.net/m0_46253997/article/details/131577438

版权

C++卫星轨道学习专栏收录该内容

2 篇文章 0 订阅

订阅专栏

卫星的运动微分方程式如下所示：

对卫星进行轨道外推即利用龙格库塔法解上图所示的二阶与一阶微分方程

对于初始时刻的位置position[x,y,z]和速度velocity[vx,vy,vz]，都是常数，故以上微分方程是常微分方程，而位置的一阶导是速度，故接下来就是利用四阶龙格库塔法解上面的六个常微分方程。

代码如下所示：

#include<iostream>
#include <vector>
#include <cmath>

using namespace std;

vector<double> yijieweifen(vector<double> position,vector<double> velocity){
	vector<double> dy;
	double miu=398600.44;
	dy.push_back(-(miu*position[0])/pow(pow(position[0],2)+pow(position[1],2)+pow(position[2],2),1.5));
	dy.push_back(-(miu*position[1])/pow(pow(position[0],2)+pow(position[1],2)+pow(position[2],2),1.5));
	dy.push_back(-(miu*position[2])/pow(pow(position[0],2)+pow(position[1],2)+pow(position[2],2),1.5));
	dy.push_back(velocity[0]);
	dy.push_back(velocity[1]);
	dy.push_back(velocity[2]);

	return dy;



}

vector<double> qiuK(int h, vector<double>position, vector<double>velocity)//求K
{
	vector<double>k1;
	k1.push_back(yijieweifen( position,velocity)[0]);
	k1.push_back(yijieweifen( position,velocity)[1]);
	k1.push_back(yijieweifen( position,velocity)[2]);
	k1.push_back(yijieweifen( position,velocity)[3]);
	k1.push_back(yijieweifen( position,velocity)[4]);
	k1.push_back(yijieweifen( position,velocity)[5]);

	vector<double> position1;
	vector<double> velocity1;
	velocity1.push_back(velocity[0]+h*k1[0]/2);
	velocity1.push_back(velocity[1]+h*k1[1]/2);
	velocity1.push_back(velocity[2]+h*k1[2]/2);
	position1.push_back(position[0]+h*k1[3]/2);
	position1.push_back(position[1]+h*k1[4]/2);
	position1.push_back(position[2]+h*k1[5]/2);
	vector<double>k2;
	k2.push_back(yijieweifen( position1,velocity1)[0]);
	k2.push_back(yijieweifen( position1,velocity1)[1]);
	k2.push_back(yijieweifen( position1,velocity1)[2]);
	k2.push_back(yijieweifen( position1,velocity1)[3]);
	k2.push_back(yijieweifen( position1,velocity1)[4]);
	k2.push_back(yijieweifen( position1,velocity1)[5]);


	vector<double> position2;
	vector<double> velocity2;
	velocity2.push_back(velocity[0]+h*k2[0]/2);
	velocity2.push_back(velocity[1]+h*k2[1]/2);
	velocity2.push_back(velocity[2]+h*k2[2]/2);
	position2.push_back(position[0]+h*k2[3]/2);
	position2.push_back(position[1]+h*k2[4]/2);
	position2.push_back(position[2]+h*k2[5]/2);
	vector<double>k3;
	k3.push_back(yijieweifen( position2,velocity2)[0]);
	k3.push_back(yijieweifen( position2,velocity2)[1]);
	k3.push_back(yijieweifen( position2,velocity2)[2]);
	k3.push_back(yijieweifen( position2,velocity2)[3]);
	k3.push_back(yijieweifen( position2,velocity2)[4]);
	k3.push_back(yijieweifen( position2,velocity2)[5]);

	vector<double> position3;
	vector<double> velocity3;
	velocity3.push_back(velocity[0]+h*k3[0]);
	velocity3.push_back(velocity[1]+h*k3[1]);
	velocity3.push_back(velocity[2]+h*k3[2]);
	position3.push_back(position[0]+h*k3[3]);
	position3.push_back(position[1]+h*k3[4]);
	position3.push_back(position[2]+h*k3[5]);
	vector<double>k4;
	k4.push_back(yijieweifen( position3,velocity3)[0]);
	k4.push_back(yijieweifen( position3,velocity3)[1]);
	k4.push_back(yijieweifen( position3,velocity3)[2]);
	k4.push_back(yijieweifen( position3,velocity3)[3]);
	k4.push_back(yijieweifen( position3,velocity3)[4]);
	k4.push_back(yijieweifen( position3,velocity3)[5]);

	//以上都是求K的值
	//对速度和位置进行更新
	velocity[0]=velocity[0]+h*(k1[0]+2*k2[0]+2*k3[0]+k4[0])/6;
	velocity[1]=velocity[1]+h*(k1[1]+2*k2[1]+2*k3[1]+k4[1])/6;
	velocity[2]=velocity[2]+h*(k1[2]+2*k2[2]+2*k3[2]+k4[2])/6;
	position[0]=position[0]+h*(k1[3]+2*k2[3]+2*k3[3]+k4[3])/6;
	position[1]=position[1]+h*(k1[4]+2*k2[4]+2*k3[4]+k4[4])/6;
	position[2]=position[2]+h*(k1[5]+2*k2[5]+2*k3[5]+k4[5])/6;

	vector<double> data;
	data.push_back(velocity[0]);
	data.push_back(velocity[1]);
	data.push_back(velocity[2]);
	data.push_back(position[0]);
	data.push_back(position[1]);
	data.push_back(position[2]);

	return data;

}
void main()
{
	int h=60;int t=60;
	vector<double> velocity;
	vector<double> position;
	velocity.push_back(-0.535832);
	velocity.push_back(6.773181);
	velocity.push_back(3.677537);
	position.push_back(6662.055595);
	position.push_back(407.044782);
	position.push_back(221.005843);

	for(int i=0;i<10;i++)
	{
		vector<double> data;
		data=qiuK(h,position,velocity);
		t+=60;
		cout<< "time" << " vx" << " vy" << " vz" << " x" << " y" << " z"<< endl; 
		cout<< t << " " << data[0] <<" " << data[1] <<" " << data[2] <<" " << data[3] <<" " << data[4] <<" " << data[5]<< endl;
		//更新初始值
		velocity[0]=data[0];
		velocity[1]=data[1];
		velocity[2]=data[2];
		position[0]=data[3];
		position[1]=data[4];
		position[2]=data[5];

	}

	system("pause");
}