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<=180
    E = atand(tanE)*2;
else
    E=  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>360
    int=floor(M_all1/360);
    M_all=M_all1-int*360;
else
    M_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>=0
    f=rad2deg(atan(tanf2)*2);
    f_all=f;
elseif E_rad>=pi
    f=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坐标系,其数据在作者的数据包里可以免费下载,最后得到四个图

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

如若内容造成侵权/违法违规/事实不符,请联系我们进行投诉反馈qq邮箱809451989@qq.com,一经查实,立即删除!

相关文章

创建__init__()方法

自学python如何成为大佬(目录):https://blog.csdn.net/weixin_67859959/article/details/139049996?spm1001.2014.3001.5501 在创建类后&#xff0c;可以手动创建一个__init__()方法。该方法是一个特殊的方法&#xff0c;类似Java语言中的构造方法。每当创建一个类的新实例时…

Java学习笔记 集合的使用

在实际的项目中开发过程中&#xff0c;会有很多的对象&#xff0c;如何高效、方便的管理这些对象&#xff0c;是影响程序性能与可维护性的重要环节。在Java语言中为这个问题提供了一套完美的解决方案&#xff0c;也就是接下来要介绍的集合框架。 1.1 集合框架的结构 从Collect…

【IC验证】一文速通多通道数据整型器(MCDF)

目录 01 README 02 MCDF设计结构 2.1 功能描述 2.2 设计结构 2.3 接口与时序 2.3.1 系统信号接口 2.3.2 通道从端接口 2.3.3 整形器接口 2.3.4 控制寄存器接口 2.3.4.1 接口时序图 2.3.4.2 各数据位信息 03 验证框图 3.1 reg_pkg 3.1.1 reg_trans 3.1.2 reg_driv…

[work] AI算法八股总结

一、深度学习面试宝典 amusi/Deep-Learning-Interview-Book: 深度学习面试宝典&#xff08;含数学、机器学习、深度学习、计算机视觉、自然语言处理和SLAM等方向&#xff09; (github.com)https://github.com/amusi/Deep-Learning-Interview-Book 深度学习八股https://github…

ChatTTS:开源最强文本转真人语音工具

目录 1.前言 2.详细介绍 2.1 什么是ChatTTS 2.2 项目地址: 2.3 应用特点: 3.如何安装和使用 3.1.谷歌colab 3.1.1.点击链接 3.1.2 进行保存 3.1.3 按照流程依次点击运行 3.1.4 填写自己需要转的文字 3.2 本地运行 3.2.1 下载或克隆项目源码到本地 3.2.2 …

JAVA多线程与IO流知识总结

文章目录 IO流体系字节流FileOutputStreamFileInputStreamtry-catch处理 字符流FileReaderFIleWriter原理分析 使用场景字节缓冲流BufferedInputStreamBufferedOutputStream 字符缓冲流BufferedReaderBufferedWriter 转换流InputStreamReaderOutputStreamWriter 序列化流Object…

SQL刷题笔记day8——SQL进阶——表与索引操作

目录 1 创建一张新表 2 修改表 3 删除表 4 创建索引 5 删除索引 1 创建一张新表 我的答案 create table if not exists user_info_vip (id int(11) primary key auto_increment Comment自增ID, # 有了主键就不用写not nul了 uid int(11) unique not null Comment用户ID, …

MMPose-RTMO推理详解及部署实现(下)

目录 前言一、RTMO推理(Python)1. RTMO预测2. RTMO预处理3. RTMO后处理4. RTMO推理 二、RTMO推理(C)1. ONNX导出2. RTMO预处理3. RTMO后处理4. RTMO推理 三、RTMO部署1. 源码下载2. 环境配置2.1 配置CMakeLists.txt2.2 配置Makefile 3. ONNX导出4. engine生成5. 源码修改6. 运行…

Git 恢复已删除的branch

六一节晚上改了点code, 做完之后commit, 然后误删了这个branch, 并且新建了branch. 那么怎样恢复已删除的branch呢&#xff1f; 网上查询一番&#xff0c;找到了答案&#xff1a; 1. git reflog 找到被删的branch中最后一笔commit, 记录它的SHA1。 怎么看SHA1是被删除的bra…

Mac修改Mysql8.0密码

转载请标明出处&#xff1a;http://blog.csdn.net/donkor_/article/details/139392605 文章目录 前言修改密码Step1:修改my.conf文件Step2:添加配置skip-grant-tablesStep3:重启mysql服务Step4:进入mysqlStep5:刷新权限Step6:修改密码Step7:再次刷新权限Step8:删除/注释 skip-…

本地安装AI大模型

使用ollmam安装llmama3等模型 1.打开ollmam下载对应系统的软件&#xff0c;安装即可 官网&#xff1a;Ollama&#xff0c; 安装直接点就就行了&#xff0c;没有其他操作 2.安装模型 在官网找到对于的模型下载命令 记录命令:ollama run llama3 打开一个cmd窗口&#xff0c;输…

opencv-python(三)

马赛克 face img[162:428,297:527] # 人脸坐标区域face face[::10,::10] # 每10个中取出一个像素&#xff0c;马赛克face np.repeat(face, 10, axis0) # 行方向重复10次face np.repeat(face, 10, axis1) # 列方向重复10次img[162:428,297:527] face[:266,:230] # 填充&a…

54. 螺旋矩阵【rust题解】

题目 给你一个 m 行 n 列的矩阵 matrix &#xff0c;请按照 顺时针螺旋顺序 &#xff0c;返回矩阵中的所有元素。 示例 示例 1 输入&#xff1a;matrix [[1,2,3],[4,5,6],[7,8,9]] 输出&#xff1a;[1,2,3,6,9,8,7,4,5] 示例 2 输入&#xff1a;matrix [[1,2,3,4],[5,6,…

贴片和直插型IRM红外遥控接收头引脚定义和规格参数及使用注意事项

红外遥控接收头使用注意事项 引脚定义存在不同 红外遥控接收头大量使用在家用电器的遥控中&#xff0c;属于价廉物美的一种光电接收器件&#xff0c;批量价格约0.3元左右。 多数遥控接收头的引脚定义是OUT,GND,VCC&#xff0c;另有引脚定义不同为OUT,VCC,GND&#xff0c;记住…

基于STM32的水库预警系统的Proteus仿真

文章目录 一、水库预警系统1.题目要求2.思路2.1 OLED显示汉字2.2 水质传感器等等2.3 步进电机2.4 驱动水泵 3.仿真图3.1 未仿真时3.2 开始仿真&#xff0c;OLED开始显示3.3 提高水位&#xff0c;开启阀门和预警3.4 通过按键增大水位阈值&#xff0c;取消报警 4.仿真程序4.1 程序…

车联网安全入门——ICSim模拟器使用

文章目录 车联网安全入门——ISCim模拟器使用介绍主要特点&#xff1a;使用场景&#xff1a; 安装使用捕获can流量candumpcansnifferwiresharkSavvyCAN主要特点&#xff1a;使用场景&#xff1a; 重放can报文cansendSavvyCAN 总结 车联网安全入门——ISCim模拟器使用 &#x1…

前端表单校验完成之后,点击确认功能无反应FormInstance, FormRules

**产生原因&#xff1a;可能是在el-form 中添加的ref 前面加了“&#xff1a;”&#xff0c;也可能是ref中的值写错了** FormInstance, FormRules

Unity + 雷达 粒子互动(待更新)

效果预览: 花海(带移动方向) VFX 实例 脚本示例 使用TouchScript,计算玩家是否移动,且计算移动方向 using System.Collections; using System.Collections.Generic; using TouchScript; using TouchScript.Pointers; using UnityEngine; using UnityEngine.VFX;public …

装饰器,状态管理和if判断(HarmonyOS学习第六课)

Builder装饰器-自定义构建函数 前面介绍了如何创建一个自定义组件。该自定义组件内部UI结构固定&#xff0c;仅与使方法进行数据传递。ArkUI还提供了一种更轻量的UI 元素复用机制Builder&#xff0c;Builder 所装饰的函数遵循build( )函数语法规则&#xff0c;开发者可以将重…

网页安全登陆的设计思路

对于Web网站来讲,不管是企业内容信息化系统,还是公共站点(博客、音视频站等),都有需要用户注册和登录的功能。用以识别用户、信息交互、信息隔离以及商业行为等场景。用户数据已成为网站的重要资产。保护用户信息(数据)是网站安全运行的关键任务。本文以用户安全登录的场…