关于DE405星历表的调用

最新推荐文章于 2024-04-07 19:57:09 发布

Jingyu

最新推荐文章于 2024-04-07 19:57:09 发布

阅读量5.8k

点赞数 14

文章标签：星历表

本文链接：https://blog.csdn.net/weixin_43426685/article/details/102800349

版权

写在前面

关于星历表的调用网上的教程都太老了，大多是基于32位的操作系统做的，然后fortran的编写软件CVF6.5在win10上配置又极其麻烦，然后就是网上能把调用从头到尾说清楚的其实没有。
然后这个东西是最基础的部分，是为了算三体摄动用的，我在做本科毕业设计的时候就在看，但是一直没弄明白是怎么调用的（其实就是完全不知道网上这些个教程在说什么），后来为了偷懒选用matlab读表做也做掉了

唯一能讲的比较清楚的是这篇博文：
https://blog.csdn.net/ddd_1991/article/details/81738037
然后这两天又恰巧要用到这个东西又在看这个，发现这篇博文里很多东西写的不够详细并且存在一些错误，因此我在这篇文章的基础之上详细说一下星历表的调用

（其实是自己想做一个整理不然下次又忘了怎么弄）

偷懒的matlab方法（想看fortran调用的直接跳过这部分）

1.打开STK，按照以下截图操作

①

在这里插入图片描述
打开的时候选择create之后注意在这里选择你要做仿真的开始时间和结束时间

②

在这里插入图片描述

③

点击analysis中的 report & graph manager…
在这里插入图片描述

④

在这里插入图片描述
之后选择sun vector j2000，左侧选择use step size/time bound，步长根据实际仿真输出选择，time bound也根据实际要求进行选择

⑤

输出结果如下：
在这里插入图片描述
值得注意的是，这里的三轴信息是相对于该确定位置卫星的，因此需要在步骤4中选择卫星位置再输出一次，将两次得到的结合起来得到实际的太阳相对于地球的的位置信息
月球同理，选择lunar vector J2000即可

⑥

由于⑤中输出的信息可以直接复制，将其复制到excel中，分割之后确保一个数据占用一个单元格，由于matlab对数据操作的简便性，直接调用语句sun_poistion=xlsread(filename)即可
将sun_poistion看成矩阵，单次循环用一行即可，每次循环则矩阵行数加一知道循环结束即可

fortran调用星历表

1、准备工作1

vs和xe
这里就不讲如何安装这两个软件了，详细可以看这篇里面有讲：
添加链接描述

值得注意的是，xe版本一定要比vs高，不然不支持在vs中进行fortran操作的

2、准备工作2

在下面的网站中下载如图所示的所有文件
(JPL喷气推进实验室的ftp网站)
这里感谢我在美国的两位同学帮我下载全了这些东西
https://ssd.jpl.nasa.gov/?planet_eph_export
ftp://ssd.jpl.nasa.gov/pub/eph/planets/ascii/
ftp://ssd.jpl.nasa.gov/pub/eph/planets/fortran
在这里插入图片描述

3、开始

打开asc2eph.f

在这里插入图片描述
将红色画线上的的两条注释语句任选一条删掉前面的注释符号C，我在这里选择里第二条，因为有一个地方的注释写的4用在unix下，我不是很懂这部分所以所有关于nrecl的值都取了1
当然也有取4的，例如在李云飞师兄的文章中就取值为4，文章链接如下：
https://wenku.baidu.com/view/4cd0e24df7ec4afe04a1df13.html

打开testeph.f

将testeph.f中所有关于nrecl的值全部赋为1，并且将FSIZER3函数中的KSIZE赋值为2036
在这里插入图片描述
在testeph.f中找到state子函数，取消子函数中FSIZER3前的注释符号

在windows菜单中找到你所下XE对应版本的compiler
在这里插入图片描述
然后依次输入以下命令：

cd C:\xxxxxxxxxx

(就是你存放准备工作2中文件的文件夹路径)
在这里插入图片描述

copy header.405+ascp2000.405+ascp2020.405+ascp2040.405 infile.405

在这里插入图片描述

ifort asc2eph.f

在这里插入图片描述

asc2eph < infile.405

在这里插入图片描述

ifort testeph.f

在这里插入图片描述

testeph < testpo.405

在这里插入图片描述

exit

至此，第一部分操作完成,文件夹内容应如下图所示
在这里插入图片描述

4、IAU SOFA的下载

在下面的网站中：
http://www.iausofa.org/index.html
下载得到cal2jd.f和epv00.f这两个文件

5、项目创建

打开vs选择创建新项目选择如下画圈的的项目
在这里插入图片描述
然后项目所有文件如下图所示：

其中：
JPLEPH是步骤3生成的，将其复制到该vs项目下
cal2jd.for和epv00.for是步骤4下载的
testeph.f是下载得到的，这里要处理一下，删除KSIZER3前的所有代码，如图红线

Get_ephemeris.for为主程序，代码如下：（这部分我只是在儒略日的使用上进行了修改，其他代码还是用的别人的）
(之前那个人的代码里关于儒略日的使用是有问题的，无论是IAU还是de405都用的是儒略日而不是简约儒略日)


      program Get_ephemeris

      implicit none
C++++++++++++++++++++++++++++
C     For routine iau_cal2jd(year,month,day,djm0,djm,djm_status)
C     year,month,day    integer     year, month, day in Gregorian calendar    
C     djm0        d     MJD zero-point: always 2400000.5
C     djm         d     Julian Date for 0 hrs
C     djm_status  i     status:
      INTEGER(4) year,month,day,djm_status
      DOUBLE PRECISION djm0,djm

C++++++++++++++++++++++++++++
C     For  routine iau_EPV00(djm1,djm2,pvh,pvb,jstat)  
C     djm1     d        TDB date part A (Note 1)
C     djm2     d        TDB date part B (Note 1)
C     pvh      d(3,2)   heliocentric Earth position/velocity (au,au/day)
C     pvb      d(3,2)   barycentric Earth position/velocity (au,au/day)
C     jstat    i        status: 0 = OK
C                              +1 = warning: date outside 1900-2100 AD   
      INTEGER jstat
      DOUBLE PRECISION djm1,djm2
      DIMENSION pvh(3,2),pvb(3,2)
      DOUBLE PRECISION pvh,pvb

C++++++++++++++++++++++++++++
C     For routine PLEPH(djm3,ntarg,ncent,rrd)
C     ntarg = INTEGER NUMBER OF 'TARGET' POINT.
C
C     ncent = INTEGER NUMBER OF CENTER POINT.
C
C            THE NUMBERING CONVENTION FOR 'NTARG' AND 'NCENT' IS:
C
C                1 = MERCURY           8 = NEPTUNE
C                2 = VENUS             9 = PLUTO
C                3 = EARTH            10 = MOON
C                4 = MARS             11 = SUN
C                5 = JUPITER          12 = SOLAR-SYSTEM BARYCENTER
C                6 = SATURN           13 = EARTH-MOON BARYCENTER
C                7 = URANUS           14 = NUTATIONS (LONGITUDE AND OBLIQ)
C                            15 = LIBRATIONS, IF ON EPH FILE
C     rrd =  OUTPUT 6-WORD D.P. ARRAY CONTAINING POSITION AND VELOCITY
C            OF POINT 'NTARG' RELATIVE TO 'NCENT'. 
      DOUBLE PRECISION djm3
      DIMENSION rrd(6)
      DOUBLE PRECISION rrd 

      year = 2011
      month = 5
      day = 1  

      call iau_cal2jd(year,month,day,djm0,djm,djm_status)

      djm1 = djm
      djm2 = djm0
      
      
      call iau_EPV00(djm1,djm2,pvh,pvb,jstat)

      djm3 = djm0+djm
      
      call PLEPH(djm3,3,11,rrd)

      write(*,*) "查询地球相对于太阳位置速度的时间为（年/月/日）："
      write(*,*) year,month,day
      write(*,100)
      write(*,*) "地球相对于太阳的位置和速度为(au,au/y) (by IAU EPV00):"
      write(*,*) pvh(1,1),pvh(2,1),pvh(3,1),pvh(1,2),pvh(2,2),pvh(3,2)
      write(*,100)
      write(*,*) "地球相对于太阳的位置和速度为(au,au/y) (by JPL PLEPH):"
      write(*,*) rrd(1),rrd(2),rrd(3),rrd(4),rrd(5),rrd(6)
      
100   format(/)

      stop
      end

测试结果如下：
在这里插入图片描述
可以看出和之前那位博主的结果一样，测试成功

当然你在做实际仿真的时候是要精确到秒的，那么修改儒略日的计算方式即可，加入hour和second就行了，关于儒略日的计算我就不详细讲了，应该都学过，没学过的百度学一下也很快

在C++中调用

这部分我还没弄完，原理是将fortran编译成动态链接库再在c++中进行调用，详细的可以看下面这篇文章，我做完了也会发的
https://www.cnblogs.com/xuelixue/articles/10170037.html

有任何使用上的问题欢迎和我联系，邮箱在这里：ljy624996209@qq.com

Jingyu

关注

14
点赞
踩
46

收藏

觉得还不错? 一键收藏
4
评论
关于DE405星历表的调用

写在前面关于星历表的调用网上的教程都太老了，大多是基于32位的操作系统做的，然后fortran的编写软件CVF6.5在win10上配置又极其麻烦，然后就是网上能把调用从头到尾说清楚的其实没有。然后这个东西是最基础的部分，是为了算三体摄动用的，我在做本科毕业设计的时候就在看，但是一直没弄明白是怎么调用的（其实就是完全不知道网上这些个教程在说什么），后来为了偷懒选用matlab读表做也做掉了唯一...
复制链接

扫一扫