TH方程学习(4)

一、背景介绍

在本节将会对TH方程打包成一个函数,通过输入目标星的轨道要素,追踪星在目标星VVLH坐标系下的相对位置和速度,以及预报的时间,得到预报时刻追踪星在VVLH坐标系下的位置和速度,以及整个状态转移矩阵。

合并(1)(2),得到归一化形式的XZ方向的状态转移矩阵\Phi_{XZ},(3),得到归一化形式的Y方向的状态转移矩阵\Phi_{Y}。所以合并后的形式可以写为\tilde{\Phi}=\begin{bmatrix} \Phi_{XZ}(1,1) & 0 &\Phi_{XZ}(1,2) & \Phi_{XZ}(1,3) & 0&\Phi_{XZ}(1,4) \\0 & \Phi_{Y}(1,1) &0& 0 & \Phi_{Y}(1,2) &0 \\ \Phi_{XZ}(2,1) & 0 & \Phi_{XZ}(2,2) & \Phi_{XZ}(2,3) & 0 &\Phi_{XZ}(2,4) \\ \Phi_{XZ}(3,1) & 0 & \Phi_{XZ}(3,2) & \Phi_{XZ}(3,3)& 0 & \Phi_{XZ}(3,4)\\ 0 & \Phi_{Y}(2,1) & 0 & 0 &\Phi_{Y}(2,2) & 0 \\ \Phi_{XZ}(4,1) & 0 & \Phi_{XZ}(4,2) & \Phi_{XZ}(4,3) &0 & \Phi_{XZ}(4,4) \end{bmatrix}

根据归一化,可以写成的矩阵形式如下所示

\tilde{\boldsymbol{X}}=\Phi_{\boldsymbol{X}}\boldsymbol{X}

\Phi_{X}=\begin{bmatrix} \rho & & & & & \\ & \rho & & & & \\ & & \rho & & & \\ -esin \theta & & & \frac{1}{k^2\rho} & & \\ & -e sin \theta & & & \frac{1}{k^2\rho} & \\ & & -e sin \theta & & & \frac{1}{k^2\rho} \end{bmatrix}

\Phi_{X}^{-1}=\begin{bmatrix} \frac{1}{\rho} & & & & & \\ & \frac{1}{\rho} & & & & \\ & & \frac{1}{\rho} & & & \\ k^2e\sin\theta & & & k^2\rho & & \\ & k^2e\sin\theta & & & k^2\rho & \\ & & k^2e\sin\theta & & & k^2\rho \end{bmatrix}

根据式子

\tilde{X}_{t}=\tilde{\Phi}\tilde{X}_0\\ \Phi_{X_{t}}X_t=\tilde{\Phi}\Phi_{X_0}X_0

所以状态转移矩阵为

\Phi=\Phi_{X_{t}}^{-1}\tilde{\Phi}_{X}\Phi_{X0}

实现代码如下所示

function [rvt,Phi,rvt_phi] = TH_solver(ecc,Perigee,TA,delta_r,delta_v,t)
%TH_SOLVER 本函数旨在通过输入目标轨道的轨道要素,追踪星的初始相对位置和速度
%  输入
%  ecc    :  目标星的偏心率           0<ecc<1
%  Perigee:  目标星的近地点高度       (km)
%  TA     :  目标星的初始真近地点角   (deg) (0<TA<=360)
%  delta_r:  追踪星在目标星的VVLH位置 (km)
%  delta_v:  追踪星在目标星的VVLH速度 (km/s)
%  t      :  预报时间长度             (s)
%  输出
%  rvt    :  t时刻后,追踪星在目标星VVLH坐标系的位置和速度 (km,km/s)
%  Phi    :  状态转移矩阵               
%  rvt_phi:  通过状态转移矩阵计算出来的追踪星在目标星VVLH坐标系的位置和速度(用来与rvt做一个对比)
%% 参考文献<New State Transition Matrix for Relative Motion on an Aribitrary Elliptical Orbit>
%% 作者 Yamanaka Koji
global miu Re
miu = 3.986e5;
Re  = 6378.137;
% 计算半长轴
sma =    (Re+Perigee)/(1-ecc);
% 计算半通径
p  =     sma*(1-ecc^2);
% 计算角动量
h  =     sqrt(p*miu);
% 计算归一化的相对位置速度
rho     = 1+ecc*cosd(TA);
r_norm  = rho*delta_r;                                   %式子87
k       = sqrt(h/p^2);
v_norm  = -ecc*sind(TA)*delta_r+(1/(k^2*rho))*delta_v;   %式子87 % 计算不同时刻的真近地点角
tanf=tand(TA/2);
ee=sqrt((1+ecc)/(1-ecc));
tanE=tanf/ee;
% 求偏近地点角
if TA<=180E = atand(tanE)*2;
elseE=  atand(tanE)*2+360;
end
% 求平近地点角
M = rad2deg(deg2rad(E)-ecc*sind(E));
% 求平均角速度
n = sqrt(miu/sma^3);
% 平近地点角
M_all1=M+rad2deg(n*t);
if M_all1>360int=floor(M_all1/360);M_all=M_all1-int*360;
elseM_all=M_all1;
end% 利用开普勒方程求解每个时刻的偏近地点角,真近地点角
MM=deg2rad(M_all);
E_rad=Kepler_function(ecc,MM);
tanf2=sqrt((1+ecc)/(1-ecc))*tan(E_rad/2);
if E_rad<pi && E_rad>=0f=rad2deg(atan(tanf2)*2);f_all=f;
elseif E_rad>=pif=rad2deg(atan(tanf2)*2)+360;f_all=f;
end% 计算X-Z矩阵
s0            = rho*sind(TA); c0= rho*cosd(TA);
Phi0_inv      =  zeros(4,4);
Phi0_inv(1,1) =  1-ecc^2;
Phi0_inv(1,2) =  3*ecc*(s0/rho)*(1+1/rho);
Phi0_inv(1,3) =  -ecc*s0*(1+1/rho);
Phi0_inv(1,4) =  -ecc*c0+2;
Phi0_inv(2,2) =  -3*(s0/rho)*(1+ecc^2/rho);
Phi0_inv(2,3) =  s0*(1+1/rho);
Phi0_inv(2,4) =  c0-2*ecc;
Phi0_inv(3,2) =  -3*(c0/rho+ecc);
Phi0_inv(3,3) =  c0*(1+1/rho)+ecc;
Phi0_inv(3,4) =  -s0;
Phi0_inv(4,2) =  3*rho+ecc^2-1;
Phi0_inv(4,3) =  -rho^2;
Phi0_inv(4,4) =  ecc*s0;
Phi0_inv1     =  Phi0_inv.*(1/(1-ecc^2));% 根据f_all 计算每个时刻对应的转移矩阵 X-Z
xz0     =   [r_norm(1);r_norm(3);v_norm(1);v_norm(3)];
y0      =   [r_norm(2);v_norm(2)];
% 计算转移初始值
hatxz0  =   Phi0_inv1*xz0;
% 转移矩阵
% 已知 t 时刻的真近地点角
theta   =    f_all;
% 计算矩阵的参数
rho1    =    1+ecc*cosd(theta);
s1      =    rho1*sind(theta);
c1      =    rho1*cosd(theta);
ds1     =    cosd(theta)+ecc*cosd(2*theta);
dc1     =    -(sind(theta)+ecc*sind(2*theta));
J1      =    k^2*t;
%  转移矩阵的参数
Phi_theta       =  zeros(4);
Phi_theta(1,1)  =  1;
Phi_theta(1,2)  =  -c1*(1+1/rho1);
Phi_theta(1,3)  =  s1*(1+1/rho1);
Phi_theta(1,4)  =  3*rho1^2*J1;
Phi_theta(2,2)  =  s1;
Phi_theta(2,3)  =  c1;
Phi_theta(2,4)  =  2-3*ecc*s1*J1;
Phi_theta(3,2)  =  2*s1;
Phi_theta(3,3)  =  2*c1-ecc;
Phi_theta(3,4)  =  3*(1-2*ecc*s1*J1);
Phi_theta(4,2)  =  ds1;
Phi_theta(4,3)  =  dc1;
Phi_theta(4,4)  =  -3*ecc*(ds1*J1+s1/rho1^2);
%  利用XZ转移矩阵求出该时刻的位置和速度
hatxz           =   Phi_theta*hatxz0;
xt              =   hatxz(1,1)/rho1;
zt              =   hatxz(2,1)/rho1;
vxt             =   k^2*(hatxz(1,1)*ecc*sind(theta)+hatxz(3,1)*rho1);
vzt             =   k^2*(hatxz(2,1)*ecc*sind(theta)+hatxz(4,1)*rho1);
%  Y转移矩阵的参数
Phi_thetay      =   zeros(2);
Phi_thetay(1,1) =   cosd(theta-TA);
Phi_thetay(1,2) =   sind(theta-TA);
Phi_thetay(2,1) =   -sind(theta-TA);
Phi_thetay(2,2) =    cosd(theta-TA);
haty(1,:)       =    Phi_thetay*y0;
yt              =    haty(1,1)/rho1;
vyt             =    k^2*(haty(1,1)*ecc*sind(theta)+haty(1,2)*rho1);
rvt             =    [xt,yt,zt,vxt,vyt,vzt];
% 求出归一化形式的状态转移矩阵
% 首先X-Z形式的矩阵
Phi_XZ  = Phi_theta*Phi0_inv1;
Phi_Y   = Phi_thetay;
% 归一化的状态转移矩阵
Phi_hat = zeros(6,6);
Phi_hat(1,1) = Phi_XZ(1,1);
Phi_hat(1,3) = Phi_XZ(1,2);
Phi_hat(1,4) = Phi_XZ(1,3);
Phi_hat(1,6) = Phi_XZ(1,4);
Phi_hat(3,1) = Phi_XZ(2,1);
Phi_hat(3,3) = Phi_XZ(2,2);
Phi_hat(3,4) = Phi_XZ(2,3);
Phi_hat(3,6) = Phi_XZ(2,4);
Phi_hat(4,1) = Phi_XZ(3,1);
Phi_hat(4,3) = Phi_XZ(3,2);
Phi_hat(4,4) = Phi_XZ(3,3);
Phi_hat(4,6) = Phi_XZ(3,4);
Phi_hat(6,1) = Phi_XZ(4,1);
Phi_hat(6,3) = Phi_XZ(4,2);
Phi_hat(6,4) = Phi_XZ(4,3);
Phi_hat(6,6) = Phi_XZ(4,4);
% Y
Phi_hat(2,2) = Phi_Y(1,1);
Phi_hat(2,5) = Phi_Y(1,2);
Phi_hat(5,2) = Phi_Y(2,1);
Phi_hat(5,5) = Phi_Y(2,2);
% 求出初始时间的状态转移矩阵
PhiX0=zeros(6,6);
PhiX0(1,1) = rho;
PhiX0(2,2) = rho;
PhiX0(3,3) = rho;
PhiX0(4,1) = -ecc*sind(TA);
PhiX0(5,2) = -ecc*sind(TA);
PhiX0(6,3) = -ecc*sind(TA);
PhiX0(4,4) = 1/(k^2*rho);
PhiX0(5,5) = 1/(k^2*rho);
PhiX0(6,6) = 1/(k^2*rho);
% 求出末端时间的状态转移矩阵
PhiXt = zeros(6,6);
PhiXt(1,1) = 1/rho1;
PhiXt(2,2) = 1/rho1;
PhiXt(3,3) = 1/rho1;
PhiXt(4,1) = k^2*ecc*sind(theta);
PhiXt(5,2) = k^2*ecc*sind(theta);
PhiXt(6,3) = k^2*ecc*sind(theta);
PhiXt(4,4) = k^2*rho1;
PhiXt(5,5) = k^2*rho1;
PhiXt(6,6) = k^2*rho1;
Phi        = PhiXt*Phi_hat*PhiX0;
rvt_phi=Phi*[delta_r;delta_v];
end

2.利用主函数调动子函数的方法来实现功能

clc;clear
%% 本代码旨在验证TH方程的正确性
%% 验证ecc=0.7的案例
delta_r=[0.1;0.01;0.01];
delta_v=[0.0001;0.0001;0.0001];
t=0:60:1200*60;
vr=zeros(length(t),6);
vrt_phi=zeros(length(t),6);
ecc=0.7;Perigee=500;
for j=1:length(t)[vrt,Phi,vrt_phi1]=TH_solver(0.7,500,45,delta_r,delta_v,t(j));vr(j,:)=vrt;vrt_phi(j,:)=vrt_phi1';
end
figure(1)
load STKVVLh2.mat
plot(Tar_x(1:1200),Tar_z(1:1200),'b--');
axis([-20 5 -5 4])
set(gca,'XDir','reverse');
set(gca,'YDir','reverse');
set(gca,'FontName','Times New Roman')
xlabel('X axis(km)','FontName','Times New Roman')
ylabel('Z axis(km)','FontName','Times New Roman')
title(['e=',num2str(ecc),' Perigee=',num2str(Perigee),'km'],'FontName','Times New Roman')
grid on
hold on
plot(vr(:,1),vr(:,3),'r-')
plot(vrt_phi(:,1),vrt_phi(:,3),'g--')
legend(['Numercial'],['TH-Numberical'],['TH-State Transition Matrix'],'Location','southeast')
hold off
figure(2)
plot(Tar_x(1:1200),Tar_y(1:1200));
axis([-20 5 -0.8 0.6])
set(gca,'XDir','reverse');
set(gca,'FontName','Times New Roman')
xlabel('X axis(km)','FontName','Times New Roman')
ylabel('Y axis(km)','FontName','Times New Roman')
title(['e=',num2str(ecc),' Perigee=',num2str(Perigee),'km'],'FontName','Times New Roman')
grid on
hold on
plot(vr(:,1),vr(:,2),'r-')
plot(vrt_phi(:,1),vrt_phi(:,2),'g--')
legend(['Numercial'],['TH-Numberical'],['TH-State Transition Matrix'],'Location','southeast')
%% 验证ecc=0.1的案例
t=0:60:220*60;
vr=zeros(length(t),6);
vrt_phi=zeros(length(t),6);
ecc=0.1;Perigee=500;
for j=1:length(t)[vrt,Phi]=TH_solver(0.1,500,45,delta_r,delta_v,t(j));vr(j,:)=vrt;vrt_phi(j,:)=vrt_phi1';
end
figure(3)
load STKVVLh.mat
plot(Tar_x(1:220),Tar_z(1:220),'b--');
axis([-3.5 0.5 -0.6 0.3])
set(gca,'XDir','reverse');
set(gca,'YDir','reverse');
set(gca,'FontName','Times New Roman')
xlabel('X axis(km)','FontName','Times New Roman')
ylabel('Z axis(km)','FontName','Times New Roman')
title(['e=',num2str(ecc),' Perigee=',num2str(Perigee),'km'],'FontName','Times New Roman')
grid on
hold on
plot(vr(:,1),vr(:,3),'r-')
plot(vrt_phi(:,1),vrt_phi(:,3),'g--')
legend(['Numercial'],['TH-Numberical'],['TH-State Transition Matrix'],'Location','southeast')
hold off
figure(4)
plot(Tar_x(1:220),Tar_y(1:220),'b--');
axis([-3.5 0.5 -0.2 0.15])
set(gca,'XDir','reverse');
set(gca,'FontName','Times New Roman')
xlabel('X axis(km)','FontName','Times New Roman')
ylabel('Y axis(km)','FontName','Times New Roman')
title(['e=',num2str(ecc),' Perigee=',num2str(Perigee),'km'],'FontName','Times New Roman')
grid on
hold on
plot(vr(:,1),vr(:,2),'r-')
plot(vrt_phi(:,1),vrt_phi(:,2),'g--')
legend(['Numercial'],['TH-Numberical'],['TH-State Transition Matrix'],'Location','southeast')

其中 STKVVLh,STKVVLh2的数据通过STK计算得到的VVLH坐标系,其数据在作者的数据包里可以免费下载,最后得到四个图

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/web/20888.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

替换所有的问号 ---- 模拟

题目链接 题目: 分析: 我们只需要遍历字符串, 将所有?进行修改即可但是需要判断, 修改的字符不能和前面后面重复同时, 有一个细节需要处理, 就是当?在最前面时, 没有前面的符号需要判断 在最后面的时候, 没有后面的字符需要判断 代码: class Solution {public String mod…

Django ORM深度游:探索多对一、一对一与多对多数据关系的奥秘与实践

系列文章目录 Django入门全攻略&#xff1a;从零搭建你的第一个Web项目Django ORM入门指南&#xff1a;从概念到实践&#xff0c;掌握模型创建、迁移与视图操作Django ORM实战&#xff1a;模型字段与元选项配置&#xff0c;以及链式过滤与QF查询详解Django ORM深度游&#xff…

打造智能化未来:智能运维系统架构解析与应用实践

在数字化转型的大背景下&#xff0c;智能运维系统成为了企业提升效率、降低成本、增强安全性的关键利器。本文将深入探讨智能运维系统的技术架构&#xff0c;介绍其核心要素和应用实践&#xff0c;帮助读者全面了解智能运维系统的概念、优势和应用价值。 ### 1. 智能运维系统的…

如何进入 -MGMTDB目录

想想这个问题&#xff0c;大家可能觉得很简单吧&#xff0c;先不看答案&#xff0c;你试一试如何进去 1.问题现象&#xff1a; 我直接进入&#xff1a; cd -MGMTDB 直接就报错了&#xff1a; [gridhost03 _mgmtdb]$ cd -MGMTDB -bash: cd: -M: invalid option cd: usage: c…

Mujoco仿真【xml文件的学习 4】

在学习Mujoco仿真的过程中&#xff0c;mujoco的版本要选择合适。先前我将mujoco的版本升级到了mujoco-3.1.4&#xff0c;在运行act的仿真代码时遇到了问题&#xff0c;撰写了博客&#xff1a; Aloha机械臂的mujoco仿真问题记录-CSDN博客 下面在进行mujoco仿真时&#xff0c;统…

十_信号3-可重入函数

如上图所示链表&#xff0c;在插入节点的时候捕获到了信号&#xff0c;并且该信号的自定义处理方式中也调用了插入节点的函数。 在main函数中&#xff0c;使用insert向链表中插入一个节点node1&#xff0c;在执行insert的时&#xff0c;刚让头节点指向node1以后(如上图序号1)&…

linux系统下,mysql增加用户

首先&#xff0c;在linux进入mysql mysql -u root -p 然后查看当前用户&#xff1a; select user,host from user; 增加用户语句&#xff1a; CREATE USER 用户名host范围 IDENTIFIED BY 密码;

【C++】C++入门1.0

鼠鼠最近在学C&#xff0c;那么好&#xff0c;俺来做笔记了&#xff01; 目录 1.C关键字(C98) 2.命名空间 2.1.命名空间定义 2.2.命名空间的使用 3.C的输入&&输出 4.缺省参数 4.1缺省参数概念 4.2.缺省参数的分类 5.函数重载 5.1.函数重载概念 5.2.C支持函数…

深度学习-04-数值的微分

深度学习-04-数值的微分 本文是《深度学习入门2-自製框架》 的学习笔记&#xff0c;记录自己学习心得&#xff0c;以及对重点知识的理解。如果内容对你有帮助&#xff0c;请支持正版&#xff0c;去购买正版书籍&#xff0c;支持正版书籍不仅是尊重作者的辛勤劳动&#xff0c;也…

所有人都可以做的副业兼职,短剧推广,1天挣几百,附详细方法!

自从上次向大家介绍了短剧掘金项目以来&#xff0c;便陆续收到了众多朋友的询问&#xff1a;现在是否还能加入短剧掘金的大军&#xff1f;答案是肯定的。目前&#xff0c;无论是各大视频平台还是其他渠道&#xff0c;短剧掘金项目都呈现出蓬勃发展的态势。而且&#xff0c;相关…

SSM高校社团管理系统-计算机毕业设计源码86128

目 录 摘要 1 绪论 1.1研究背景与意义 1.2开发现状 1.3研究方法 1.4 ssm框架介绍 1.5论文结构与章节安排 2 高校社团管理系统系统分析 2.1 可行性分析 2.2 系统流程分析 2.2.1数据增加流程 2.2.2数据修改流程 2.2.3数据删除流程 2.3 系统功能分析 2.3.1 功能性分…

QT实现动态翻译切换

1、实现QT动态中英文切换效果 效果如下: 2、原理 因为软件本身就是中文版,所以只需准备一个英文版的翻译即可,,那就是将所有需要翻译的地方用tr包裹,然后首先执行lupdate更新一下,接着用qt的翻译软件 Qt Linguist打开ts文件进行翻译,然后保存,最后使用 lrelease发布一…

位运算专题

常见位运算总结: 1. 基础位运算 << 左移 >> 右移 ~ 按位取反 & 按位与 口诀: 有0则0 | 按位或 口诀: 有1则1 ^ 异或 口诀:相同为0,相异为1 / 无进位相加 2. 位运算的优先级 做题时, 能加括号就加括号, 无需管优先级如何 3. 给一个二进制数n, 确定…

为什么总是卡在验证真人这里无法通过验证?

最近总是在浏览某些网站的时候卡在这个“确认你是真人”的验证页面&#xff0c;无法通过真人验证&#xff0c;这是怎么回事儿&#xff1f;如何解决呢&#xff1f; 首先&#xff0c;出现这个“确认您是真人”的验证一般都是这个网站使用了 CloudFlare 的安全防护 WAF 规则才会出…

Ray Tracing in one Weekend But on CUDA

Ray Tracing in one Weekend But on CUDA 环境说明项目代码项目内容思路实现方法效果 环境说明 代码运行在Visual Studio 2019环境&#xff0c;显卡为NVIDIA GeForce GTX 1650&#xff0c;CUDA版本为11.6&#xff0c;cuDNN版本为8.4.0。具体配置方式见CUDA C/C 从入门到入土 第…

洛谷P2370yyy2015c01 的 U 盘

传送门——P2370 yyy2015c01 的 U 盘 题解&#xff1a;题目意思很好理解&#xff0c;就是说&#xff0c;当能够达到预期的U盘的最小接口&#xff08;接口越大&#xff0c;能传递的文件越大&#xff09;&#xff0c;然后我们就需要先看题目了&#xff0c;有n个文件&#xff0c;每…

【数据结构(邓俊辉)学习笔记】图02——搜索

文章目录 0. 概述1. 广度优先搜索1.1 策略1.2 实现1.3 可能情况1.4 实例1.5 多联通1.6 复杂度1.7 最短路径 2. 深度优先搜索2.1 算法2.2 框架2.3 细节2.4 无向边2.5 有向边2.6 多可达域2.7 嵌套引理 3 遍历算法的应用 0. 概述 此前已经介绍过图的基本概念以及它在计算机中的表…

vector实现后半部分

一.迭代器失效 1.定义 指原迭代器在扩容/缩容/修改后指向无效元素或无效地址处 erase的迭代器失效 2.原因&#xff1a; 1.有的编译器实现erase会缩容拷贝 2.删除最后一个后&#xff0c;其指向无效元素 VS中不允许再次使用erase完的迭代器&#xff0c;为了让编写的代码移植…

Spring系列-SpringMvc父子容器启动原理解析

1、Spring整合SpringMVC 特性&#xff1a; 说到Spring整合SpringMVC唯一的体现就是父子容器&#xff1a; 通常我们会设置父容器&#xff08;Spring&#xff09;管理Service、Dao层的Bean, 子容器(SpringMVC)管理Controller的Bean .子容器可以访问父容器的Bean, 父容器无法访…

LLM——深入探索 ChatGPT在代码解释方面的应用研究

1.概述 OpenAI在自然语言处理&#xff08;NLP&#xff09;的征途上取得了令人瞩目的进展&#xff0c;这一切得益于大型语言模型&#xff08;LLM&#xff09;的诞生与成长。这些先进的模型不仅是技术创新的典范&#xff0c;更是驱动着如GitHub Copilot编程助手和Bing搜索引擎等广…