0%

NOVAS实践

关于使用NOVAS进行操作的部分过程。全称为Naval Observatory Vector Astrometry Software 。可用来计算天体位置等等。下面对其运用进行例子介绍。

题目内容

恒星观测[题目一]

观测时间:2020年5月1日,$23^h15^m43^s.55\mathrm{UTC}$

测站坐标:$\mathbf{W}70^\circ44’11’’.560 \quad \mathbf{S}30^\circ14’26’’.731$ 参考椭球面高度$2378 \mathrm{m} $

目标星体:$\mathbf{ICRS}[\alpha,\delta]=14^\circ34’16{‘’}.81183 \quad -12^\circ31’10{‘’}.3965 $

恒星自行:

地球指向参数:? 从IERS网站查询

求该星的:

  • 地心视位置
  • 站心坐标
  • 地平坐标

系内观测[题目二]

观测时间:2020年5月1日,$23^h15^m43^s.55\mathrm{UTC}$

测站坐标:兴隆观测站

目标星体:月球和八大行星

地球指向参数:? 从IERS网站查询

求这些目标的:

  • 地心视位置
  • 站心坐标
  • 地平坐标,且指出哪些天体可以观测。

实验准备

查询地球指向参数

首先访问IERS官网https://www.iers.org/IERS/EN/Home/home_node.html

在其Earth orientation data项目中,由于我们所查询的日期是较近的,由此选择为Daily EOP data files.

其中共有三项,

  • finals.daily(IAU1980)
  • finals.daily(IAU2000)
  • gpsrapid.daily

第三项为GPS所用数据,而我们需要使用的是第二项。

最终可得,在2020年5月1日这一天

$\mathrm{DX=0.075815[arcsec]\quad DY=0.435887[arcsec] \ UT1-UTC=-0.242517s}$

此外我们还可以查询到Leap second announcements in UTC。

即为从2017.1.1至今,$\mathrm{UTC - TAI = -37\, s}$

NOVAS程序包下载

首先需要找到NOVAS程序包。

我们可以发现好像官网坏掉了,其中,

http://aa.usno.navy.mil/software/novas/novas_info.php

http://aa.usno.navy.mil

只要是aa开头的网址无论如何也连不上。

通过朋友GasinAn指导,来到官网使用站内搜索,

https://www.usno.navy.mil/search?SearchableText=NOVAS

即可在最下方找到NOVAS_FORTRAN软件包。即为

https://www.usno.navy.mil/USNO/astronomical-applications/software-products/novas/novas-fortran/novasf3.1.tar.Z

NOVAS测试使用

对于下载的文件中,最重要的为:

  • example.f
  • NOVAS_F3.1.f
  • NOVAS_F3.1_solsys2.f

其中solsys2是利用jpl历表进行计算,因此需要使用jpl星表中的子程序。

  • jplsubs.f

此外最常用的为指导手册

  • NOVAS_F3.1_Guide.pdf

以上就是所需素材。

1
2
$ gfortran example.f NOVAS_F3.1.f NOVAS_F3.1_solsys2.f jplsubs.f -o example.out
$ ./example.out

可得显示结果:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
JPL ephemeris DE405 open. Start JD = 2451536.50  End JD = 2525008.50

NOVAS Sample Calculations
-------------------------

Geodetic location:
-70.0000000000 42.0000000000 0.0000000000

TT and UT1 Julian Dates and Delta-T:
2454580.942629 2454580.941871 65.57184500000

FK6 1307 geocentric and topocentric positions:
11.8915509892 37.6586357955
11.8915479153 37.6586695456

Moon geocentric and topocentric positions:
17.1390774264 -27.5374448869 0.002710296515
17.1031967646 -28.2902502967 0.002703785126
17.1031967646 -28.2902502967 0.002703785126

Moon zenith distance and azimuth:
81.6891016502 219.2708903405

Greenwich and local sidereal time and Earth Rotation Angle:
0.79362134148 20.12695467481 11.7956158462

Mars heliocentric ecliptic longitude and latitude and radius vector:
148.0032235906 1.8288284075 1.664218258879

Direction of zenith vector (RA & Dec) in GCRS:
20.1221838608 41.9769823554

如上所述,即为测试成功。

示例程序

实际上,通过阅读示例程序,可以发现已经写得很好,我们所做的就是改改参数就行了。由于我做这个实践时是用的3.0版本,由此以3.0版本为例的sample.f进行说明,并将其完整代码放在文末。

建立时间参数

1
2
3
4
5
6
7
8
9
*     ESTABLISH TIME ARGUMENTS
CALL JULDAT ( IYEAR, MONTH, IDAY, HOUR, UTCJD )
TTJD = UTCJD + ( LEAPS + 32.184D0 ) / 86400.D0
UT1JD = UTCJD + UT1UTC / 86400.D0
DELTAT = 32.184D0 + LEAPS - UT1UTC
CALL SETDT ( DELTAT )
WRITE ( *, * )
WRITE ( *, * ) 'TT and UT1 Julian Dates and Delta-T: '
WRITE ( *, 9020 ) TTJD, UT1JD, DELTAT

这一部分对时间进行修正,最开始输入的时间为正常的公历。

通过使用JULDAT获得此刻UTC。

且有:

由于刚查询到$\mathrm{UTC - TAI = -37\, s}$ 。对应此处LEAPS即为$\mathrm{37\, s}$

再有先前查询的$\mathrm{UT1-UTC=-0.242517s}$

此处便可得到UT1,TT二者时间。这也是后面所主要用的。

恒星测量

1
2
3
4
5
6
7
8
9
10
11
12
13
14
*     APPARENT AND TOPOCENTRIC PLACE OF STAR FK6 1307 = GRB 1830
RA2000 = 11.88299133D0
DC2000 = 37.71867646D0
PMRA = 4003.27D0
PMDEC = -5815.07D0
PARX = 109.21D0
RV = -98.8D0
CALL APSTAR ( TTJD, IEARTH, RA2000, DC2000, PMRA, PMDEC, PARX, RV,
. RA, DEC )
CALL TPSTAR ( UT1JD, GLON, GLAT, HT, RAT, DECT )
WRITE ( *, * )
WRITE ( *, * ) 'FK6 1307 geocentric and topocentric positions:'
WRITE ( *, 9010 ) RA, DEC
WRITE ( *, 9010 ) RAT, DECT

APSTAR为测量恒星的视位置,TPSTAR为测量恒星的站心位置。

视位置为在地心观测,没有大气的情况时的天体位置。对应的是真春分点、真赤道。

站心位置于站心观测,没有大气的情况时的天体位置。对应的是真春分点、真赤道。

两者差别仅为观测点所在地。

对于视位置测量,所需为J2000时赤经赤纬、自行、视差、视向速度。最终返回视位置赤经赤纬。

TPSTAR在紧接着APSTAR后面调用时,会利用先前的的相同数据,且需要提供观测站地理赤经赤纬与高度。最终返回站心位置赤经赤纬。

因此,在处理自己的问题时,只需要将对应的观测时间、赤经赤纬、自行、视差、视向速度、观测站地即可理位置修改即可。

月球测量

1
2
3
4
5
6
7
8
*     APPARENT AND TOPOCENTRIC PLACE OF THE MOON
MOON = IDSS ( 'MOON' )
CALL APPLAN ( TTJD, MOON, IEARTH, RA, DEC, DIS )
CALL TPPLAN ( UT1JD, GLON, GLAT, HT, RAT, DECT, DIST )
WRITE ( *, * )
WRITE ( *, * ) 'Moon geocentric and topocentric positions:'
WRITE ( *, 9030 ) RA, DEC, DIS
WRITE ( *, 9030 ) RAT, DECT, DIST

与先前恒星测量略有不同。

IDSS函数为通过输入观测目标的字母串,然后返回对应的编号值,本质即为字典。

根据手册,IEARTH作用不大,几乎只起到了占位符的作用,在处理自己的问题时对其保持不变即可。

通过APPLAN、TPPLAN函数对太阳系内主要天体进行计算。

1
2
3
4
5
6
7
8
9
10
*     TOPOCENTRIC PLACE OF THE MOON -- SHOULD BE SAME AS ABOVE
OBSERV(1) = GLON
OBSERV(2) = GLAT
OBSERV(3) = HT
C OBSERV(4) = ( TTJD - UT1JD ) * 86400.D0
CALL PLACE ( TTJD, 'MOON', 1, 1, STAR, OBSERV, SKYPOS )
RAT = SKYPOS(4)
DECT = SKYPOS(5)
DIST = SKYPOS(6)
WRITE ( *, 9030 ) RAT, DECT, DIST

此为另一种方法计算月球的站心位置。实际上类似于先前的APSTAR、APPLAN等函数,都是调用了最基本的PLACE函数。

地平坐标计算

1
2
3
4
5
6
7
*     POSITION OF THE MOON IN LOCAL HORIZON COORDINATES
* (POLAR MOTION IGNORED HERE)
CALL ZDAZ ( UT1JD, 0.D0, 0.D0, GLON, GLAT, HT, RAT, DECT, 1,
. ZD, AZ, RAR, DECR )
WRITE ( *, * )
WRITE ( *, * ) 'Moon zenith distance and azimuth:'
WRITE ( *, 9010 ) ZD, AZ

ZDAZ函数,其中第二第三项为极移修正,此处示例忽略了极移因此为0.D0,在稍后的操作中我们可以将极移的影响计算进去,修改为所查询到的DX、DY,代入即可。

此外需要提供测站的地理位置,观测目标的站心位置,第9项的1意为考虑大气折射,忽略则设置为0即可。最后返回天顶距与方位角。


示例程序后面还计算了格林尼治时间、火星黄道坐标、地球极移等等,但与这次试验目标无关,也就不介绍了。

试验过程

题目一

完整内容附于最后。

首先修改参数,将示例程序中的数据替换为题目给出与查询的数据:

1
2
3
4
5
6
7
 DATA IYEAR, MONTH, IDAY, HOUR /
. 2020, 5, 1, 23.262097222222224D0 /

DATA LEAPS, UT1UTC, XP, YP /
. 37, -0.2425171D0, 0.075815D0,0.435887D0 /

DATA GLON, GLAT, HT / -70.736544D0, -30.24075861D0, 2378.D0 /
1
2
3
4
5
6
RA2000 = 14.571336619444445D0
DC2000 = -12.519554583333333D0
PMRA = -354.45D0
PMDEC = +595.35D0
PARX = 164.99D0
RV = 0D0

由于题目涉及修改结果的单位,因此添加一个用于修改单位的子程序,将以小时、度为单位的赤经赤纬修改为以小时、分、秒与度、角分、角秒表示。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
SUBROUTINE thetaA ( RA, rah, ram, ras )
DOUBLE PRECISION RA
DOUBLE PRECISION rah,ram,ras

IF (RA.gt.0) THEN
rah = int(RA)
ram = ( RA - rah ) * 60
ras = ( ram - int(ram) )*60
ram = int(ram)
ENDIF

IF (RA.lt.0) THEN
RA = -RA
rah = int(RA)
ram = ( RA - rah ) * 60
ras = ( ram - int(ram) )*60
ram = int(ram)
rah = -rah
RA = -RA
ENDIF
RETURN
END

由此也可以发现,Fortran在调用子程序时,与python略有不同,即使是最开始的值,也有可能被修改,在将输入值取负号后,还要再取一次负号防止更改了初始值。

题目二

完整内容附于最后。

与题一的初始数据不同为观测站的坐标。

通过查询兴隆观测站官网可得数据。http://www.xinglong-naoc.org/html/jggk/detail-2.html

1
DATA GLON, GLAT, HT / 117.577222D0, 40.395833D0, 900.D0 /

为了程序的美观,用循环进行处理。首先建立字符串数组,然后按顺序抽取即可。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
      CHARACTER*3 PLANETS_name(9),PLANET_name

PLANETS_name= (/'MOO','MER','VEN','MAR',
.'JUP','SAT','URA','NEP','PLU'/)

PLANET_num=1

DO 30 WHILE ( PLANET_num.LT.10 )

PLANET_name=PLANETS_name(PLANET_num)
NAME = IDSS ( PLANET_name )
CALL APPLAN ( TTJD, NAME, IEARTH, RA, DEC, DIS )
PLANET_num = PLANET_num + 1
30 CONTINUE

对是否能观测到,此处理解为是否处于地平线以上,由此有

1
2
3
4
IF (zdd.GT.90.0) THEN
WRITE ( *, * )
WRITE ( *, * ) '!!!!!!',PLANET_name,' INVISIBLE'
ENDIF

实验结果

题目一

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
NOVAS Calculations
-------------------------

Geodetic location:
-70.736544000000 -30.240758610000 2378.000000000000

DATE and TIME
2020 5 1 23.262097222222224

TT and UT1 Julian Dates and Delta-T:
2458971.470055 2458971.469251 69.426517100000

The star geocentric positions:
14.589657618243 -12.604673088883
14.h 35.m 22.76743s
-12.d 36.m 16.823s

The star topocentric positions:
14.589658498483 -12.604656578923
14.h 35.m 22.77059s
-12.d 36.m 16.764s

The Star zenith distance and azimuth:
75.371843281056 96.278676337388
75.d 22.m 18.636s
96.d 16.m 43.235s

题目二

由结果可知位于兴隆观测站,此刻月球在地平线以下,不可见。

但由于观测时间为2020年5月1日,$23^h15^m43^s.55\mathrm{UTC}$

北京时间为2020年5月2日,$7^h15^m43^s.55\mathrm{UTC+8}$

此刻已经天亮,想观察其他行星也较为困难。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
NOVAS Calculations
-------------------------

Geodetic location:
117.577222000000 40.395833000000 900.000000000000

DATE and TIME
2020 5 1 23.262097222222224

TT and UT1 Julian Dates and Delta-T:
2458971.470055 2458971.469251 69.426517100000



MOO--------------------------
geocentric positions:
10.002164584844 16.760759164598 0.248441208653D-02
10.h 0.m 7.79251s
16.d 45.m 38.733s

geocentric positions:
9.999294039988 15.945468700467 0.250783067720D-02
9.h 59.m 57.45854s
15.d 56.m 43.687s

zenith distance and azimuth:
123.582884688801 356.368501682819
123.d 34.m 58.385s
356.d 22.m 6.606s

!!!!!!MOO INVISIBLE



MER--------------------------
geocentric positions:
2.427435986723 13.785296374688 0.132758063451D+01
2.h 25.m 38.76955s
13.d 47.m 7.067s

geocentric positions:
2.427527905700 13.784245926636 0.132756306462D+01
2.h 25.m 39.10046s
13.d 47.m 3.285s

zenith distance and azimuth:
65.598629569233 92.406335639752
65.d 35.m 55.066s
92.d 24.m 22.808s



VEN--------------------------
geocentric positions:
5.211194189733 27.798656381788 0.422630830955D+00
5.h 12.m 40.29908s
27.d 47.m 55.163s

geocentric positions:
5.211501824013 27.794582885383 0.422628475382D+00
5.h 12.m 41.40657s
27.d 47.m 40.498s

zenith distance and azimuth:
86.527930865671 55.627412476625
86.d 31.m 40.551s
55.d 37.m 38.685s



MAR--------------------------
geocentric positions:
21.680947417071 -15.630191999995 0.122176552998D+01
21.h 40.m 51.41070s
-15.d 37.m 48.691s

geocentric positions:
21.680949113246 -15.631844547812 0.122174163011D+01
21.h 40.m 51.41681s
-15.d 37.m 54.640s

zenith distance and azimuth:
56.025813476944 181.891602158404
56.d 1.m 32.929s
181.d 53.m 29.768s



JUP--------------------------
geocentric positions:
19.937474157665 -20.888675351049 0.482906570787D+01
19.h 56.m 14.90697s
-20.d 53.m 19.231s

geocentric positions:
19.937465607544 -20.889112954678 0.482904862752D+01
19.h 56.m 14.87619s
-20.d 53.m 20.807s

zenith distance and azimuth:
66.485135849914 208.345030217478
66.d 29.m 6.489s
208.d 20.m 42.109s



SAT--------------------------
geocentric positions:
20.277781114579 -19.840446580966 0.979417307885D+01
20.h 16.m 40.01201s
-19.d 50.m 25.608s

geocentric positions:
20.277780351685 -19.840666223154 0.979415418266D+01
20.h 16.m 40.00927s
-19.d 50.m 26.398s

zenith distance and azimuth:
63.799228019693 203.832889256775
63.d 47.m 57.221s
203.d 49.m 58.401s



URA--------------------------
geocentric positions:
2.315087711060 13.398555892598 0.208070404851D+02
2.h 18.m 54.31576s
13.d 23.m 54.801s

geocentric positions:
2.315095143519 13.398475560271 0.208070222072D+02
2.h 18.m 54.34252s
13.d 23.m 54.512s

zenith distance and azimuth:
64.557956607156 93.864638027835
64.d 33.m 28.644s
93.d 51.m 52.697s



NEP--------------------------
geocentric positions:
23.429610814108 -4.817465182521 0.305438653615D+02
23.h 25.m 46.59893s
-4.d 49.m 2.875s

geocentric positions:
23.429616636412 -4.817518811878 0.305438382057D+02
23.h 25.m 46.61989s
-4.d 49.m 3.068s

zenith distance and azimuth:
50.518358929758 147.499103457954
50.d 31.m 6.092s
147.d 29.m 56.772s



PLU--------------------------
geocentric positions:
19.806107014962 -21.983642709766 0.337215325831D+02
19.h 48.m 21.98525s
-21.d 59.m 1.114s

geocentric positions:
19.806109277109 -21.983716516423 0.337215166924D+02
19.h 48.m 21.99340s
-21.d 59.m 1.379s

zenith distance and azimuth:
68.214885202733 209.697854624252
68.d 12.m 53.587s
209.d 41.m 52.277s

sample.f内容

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
*
* NOVAS SAMPLE CALCULATIONS
*
* SEE CHAPTER 3 OF USER'S GUIDE FOR EXPLANATION
* REQUIRES SOLSYS VERSION 1 OR 2
*
IMPLICIT DOUBLE PRECISION ( A-H, O-Z )
IMPLICIT INTEGER ( I-N )

DIMENSION STAR(6), OBSERV(6), SKYPOS(7),
. POS(3), VEL(3), POSE(3), VTER(3), VCEL(3)

PARAMETER ( PI = 3.14159265358979324D0 )
PARAMETER ( DEGRAD = PI / 180.D0 )

DATA STAR, OBSERV / 12 * 0.D0 /

DATA IYEAR, MONTH, IDAY, HOUR, LEAPS, UT1UTC, XP, YP /
. 2008, 4, 24, 10.605D0, 33, -0.387845D0, -0.002D0, +0.529D0 /

DATA GLON, GLAT, HT / -70.D0, 42.D0, 0.D0 /


* FORMAT STATEMENTS FOR OUTPUT

9010 FORMAT(1X, 3(F17.12, 8X))
9020 FORMAT(1X, 2(F15.6, 8X), 1F17.12)
9030 FORMAT(1X, 2(F17.12, 8X), 1D18.12)

WRITE ( *, * )
WRITE ( *, * ) 'NOVAS Sample Calculations'
WRITE ( *, * ) '-------------------------'

* SETUP CALLS
* HIGH ACCURACY AND EQUINOX MODE ARE NOVAS DEFAULTS.
CALL HIACC
CALL EQINOX
IEARTH = IDSS ( 'EARTH' )

* WRITE OUT ASSUMED LONGITUDE, LATITUDE, HEIGHT (ITRS = WGS-84)
WRITE ( *, * )
WRITE ( *, * ) 'Geodetic location: '
WRITE ( *, 9010 ) GLON, GLAT, HT

* ESTABLISH TIME ARGUMENTS
CALL JULDAT ( IYEAR, MONTH, IDAY, HOUR, UTCJD )
TTJD = UTCJD + ( LEAPS + 32.184D0 ) / 86400.D0
UT1JD = UTCJD + UT1UTC / 86400.D0
DELTAT = 32.184D0 + LEAPS - UT1UTC
CALL SETDT ( DELTAT )
WRITE ( *, * )
WRITE ( *, * ) 'TT and UT1 Julian Dates and Delta-T: '
WRITE ( *, 9020 ) TTJD, UT1JD, DELTAT

* APPARENT AND TOPOCENTRIC PLACE OF STAR FK6 1307 = GRB 1830
RA2000 = 11.88299133D0
DC2000 = 37.71867646D0
PMRA = 4003.27D0
PMDEC = -5815.07D0
PARX = 109.21D0
RV = -98.8D0
CALL APSTAR ( TTJD, IEARTH, RA2000, DC2000, PMRA, PMDEC, PARX, RV,
. RA, DEC )
CALL TPSTAR ( UT1JD, GLON, GLAT, HT, RAT, DECT )
WRITE ( *, * )
WRITE ( *, * ) 'FK6 1307 geocentric and topocentric positions:'
WRITE ( *, 9010 ) RA, DEC
WRITE ( *, 9010 ) RAT, DECT

* APPARENT AND TOPOCENTRIC PLACE OF THE MOON
MOON = IDSS ( 'MOON' )
CALL APPLAN ( TTJD, MOON, IEARTH, RA, DEC, DIS )
CALL TPPLAN ( UT1JD, GLON, GLAT, HT, RAT, DECT, DIST )
WRITE ( *, * )
WRITE ( *, * ) 'Moon geocentric and topocentric positions:'
WRITE ( *, 9030 ) RA, DEC, DIS
WRITE ( *, 9030 ) RAT, DECT, DIST

* TOPOCENTRIC PLACE OF THE MOON -- SHOULD BE SAME AS ABOVE
OBSERV(1) = GLON
OBSERV(2) = GLAT
OBSERV(3) = HT
C OBSERV(4) = ( TTJD - UT1JD ) * 86400.D0
CALL PLACE ( TTJD, 'MOON', 1, 1, STAR, OBSERV, SKYPOS )
RAT = SKYPOS(4)
DECT = SKYPOS(5)
DIST = SKYPOS(6)
WRITE ( *, 9030 ) RAT, DECT, DIST

* POSITION OF THE MOON IN LOCAL HORIZON COORDINATES
* (POLAR MOTION IGNORED HERE)
CALL ZDAZ ( UT1JD, 0.D0, 0.D0, GLON, GLAT, HT, RAT, DECT, 1,
. ZD, AZ, RAR, DECR )
WRITE ( *, * )
WRITE ( *, * ) 'Moon zenith distance and azimuth:'
WRITE ( *, 9010 ) ZD, AZ

* GREENWICH APPARENT SIDEREAL TIME AND EARTH ROTATION ANGLE
CALL SIDTIM ( UT1JD, 0.D0, 1, GAST )
ST = GAST + GLON / 15.D0
IF ( ST .GE. 24.D0 ) ST = ST - 24.D0
IF ( ST .LT. 0.D0 ) ST = ST + 24.D0
CALL EROT ( UT1JD, 0.D0, ERA )
WRITE ( *, * )
WRITE ( *, * ) 'Greenwich and local sidereal time and ',
. 'Earth Rotation Angle:'
WRITE ( *, 9010 ) GAST, ST, ERA

* HELIOCENTRIC POSITION OF MARS IN BCRS
TDBJD = TTJD
* ABOVE IS AN APPROXIMATION GOOD TO 0.0017 SECONDS
* APPROXIMATION COULD LEAD TO ERROR OF ~50 M IN POSITION OF MARS
MARS = IDSS ( 'MARS' )
CALL SOLSYS ( TDBJD, MARS, 1, POS, VEL, IERR )
IF ( IERR .EQ. 0 ) THEN
CALL EQEC ( 0.D0, 0, POS, POSE )
CALL ANGLES ( POSE, ELON, ELAT )
WRITE ( *, * )
WRITE ( *, * ) 'Mars heliocentric ecliptic longitude and ',
. 'latitude and radius vector: '
R = DSQRT (POSE(1)**2 + POSE(2)**2 +POSE(3)**2)
WRITE ( *, 9010 ) ELON * 15.D0, ELAT, R
ELSE
WRITE ( *, * ) 'For Mars, SOLSYS returns error ', IERR
END IF

* TERRESTRIAL TO CELESTIAL TRANSFORMATION
SINLON = DSIN ( GLON * DEGRAD )
COSLON = DCOS ( GLON * DEGRAD )
SINLAT = DSIN ( GLAT * DEGRAD )
COSLAT = DCOS ( GLAT * DEGRAD )
* FORM VECTOR TOWARD LOCAL ZENITH (ORTHOGONAL TO ELLIPSOID) IN ITRS
VTER(1) = COSLAT * COSLON
VTER(2) = COSLAT * SINLON
VTER(3) = SINLAT
* TRANSFORM VECTOR TO GCRS
CALL TERCEL ( UT1JD, 0.D0, XP, YP, VTER, VCEL )
CALL ANGLES ( VCEL, RA, DEC )
WRITE ( *, * )
WRITE ( *, * ) 'Direction of zenith vector (RA & Dec) in GCRS:'
WRITE ( *, 9010 ) RA, DEC

WRITE ( *, * )
STOP
END

star.f内容[题目一]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
      IMPLICIT DOUBLE PRECISION ( A-H, O-Z )
IMPLICIT INTEGER ( I-N )

DIMENSION STAR(6), OBSERV(6), SKYPOS(7),
. POS(3), VEL(3), POSE(3), VTER(3), VCEL(3)

PARAMETER ( PI = 3.14159265358979324D0 )
PARAMETER ( DEGRAD = PI / 180.D0 )


DATA STAR, OBSERV / 12 * 0.D0 /

DATA IYEAR, MONTH, IDAY, HOUR /
. 2020, 5, 1, 23.262097222222224D0 /

DATA LEAPS, UT1UTC, XP, YP /
. 37, -0.2425171D0, 0.075815D0,0.435887D0 /

DATA GLON, GLAT, HT / -70.736544D0, -30.24075861D0, 2378.D0 /


* FORMAT STATEMENTS FOR OUTPUT

9010 FORMAT(1X, 3(F17.12, 8X))
9020 FORMAT(1X, 2(F15.6, 8X), 1F17.12)
9030 FORMAT(1X, 2(F17.12, 8X), 1D18.12)
9040 FORMAT(1X, F5.0,'d', 2X,F5.0,'m' 2X, 1F7.3,'s')
9050 FORMAT(1X, F5.0,'h', 2X,F5.0,'m' 2X, 1F9.5,'s')

WRITE ( *, * )
WRITE ( *, * ) 'NOVAS Calculations'
WRITE ( *, * ) '-------------------------'

* SETUP CALLS
* HIGH ACCURACY AND EQUINOX MODE ARE NOVAS DEFAULTS.
CALL HIACC
CALL EQINOX
IEARTH = IDSS ( 'EARTH' )

* WRITE OUT ASSUMED LONGITUDE, LATITUDE, HEIGHT (ITRS = WGS-84)
WRITE ( *, * )
WRITE ( *, * ) 'Geodetic location: '
WRITE ( *, 9010 ) GLON, GLAT, HT

* ESTABLISH TIME ARGUMENTS
CALL JULDAT ( IYEAR, MONTH, IDAY, HOUR, UTCJD )
TTJD = UTCJD + ( LEAPS + 32.184D0 ) / 86400.D0
UT1JD = UTCJD + UT1UTC / 86400.D0
DELTAT = 32.184D0 + LEAPS - UT1UTC
WRITE ( *, * )
WRITE ( *, * ) 'DATE and TIME'
WRITE ( *, * ) IYEAR,MONTH,IDAY,HOUR
WRITE ( *, * )
WRITE ( *, * ) 'TT and UT1 Julian Dates and Delta-T: '
WRITE ( *, 9020 ) TTJD, UT1JD, DELTAT

* APPARENT AND TOPOCENTRIC PLACE OF STAR FK6 1307 = GRB 1830
RA2000 = 14.571336619444445D0
DC2000 = -12.519554583333333D0
PMRA = -354.45D0
PMDEC = +595.35D0
PARX = 164.99D0
RV = 0D0
CALL APSTAR ( TTJD, IEARTH, RA2000, DC2000, PMRA, PMDEC, PARX, RV,
. RA, DEC )
WRITE ( *, * )
WRITE ( *, * ) 'The star geocentric positions:'
WRITE ( *, 9010 ) RA, DEC
CALL thetaA(RA,rah,ram,ras)
WRITE(*,9050) rah,ram,ras
CALL thetaA(DEC,decd,decm,decs)
WRITE(*,9040) decd,decm,decs

CALL TPSTAR ( UT1JD, GLON, GLAT, HT, RAT, DECT )
WRITE ( *, * )
WRITE ( *, * ) 'The star topocentric positions:'
WRITE ( *, 9010 ) RAT, DECT
CALL thetaA(RAT,rath,ratm,rats)
WRITE ( *, 9050 ) rath,ratm,rats
CALL thetaA(DECT,dectd,dectm,dects)
WRITE ( *, 9040 ) dectd,dectm,dects


* POSITION OF THE STAR IN LOCAL HORIZON COORDINATES

CALL ZDAZ ( UT1JD, XP, YP, GLON, GLAT, HT, RAT, DECT, 1,
. ZD, AZ, RAR, DECR )
WRITE ( *, * )
WRITE ( *, * ) 'The Star zenith distance and azimuth:'
WRITE ( *, 9010 ) ZD, AZ
CALL thetaA(ZD,zdd,zdm,zds)
WRITE ( *, 9040 ) zdd,zdm,zds
CALL thetaA(AZ,azd,azm,azs)
WRITE ( *, 9040 ) azd,azm,azs

END


SUBROUTINE thetaA ( RA, rah, ram, ras )
DOUBLE PRECISION RA
DOUBLE PRECISION rah,ram,ras

IF (RA.gt.0) THEN
rah = int(RA)
ram = ( RA - rah ) * 60
ras = ( ram - int(ram) )*60
ram = int(ram)
ENDIF

IF (RA.lt.0) THEN
RA = -RA
rah = int(RA)
ram = ( RA - rah ) * 60
ras = ( ram - int(ram) )*60
ram = int(ram)
rah = -rah
RA = -RA
ENDIF
RETURN
END

planet.f内容[题目二]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
      IMPLICIT DOUBLE PRECISION ( A-H, O-Z )
IMPLICIT INTEGER ( I-N )

DIMENSION STAR(6), OBSERV(6), SKYPOS(7),
. POS(3), VEL(3), POSE(3), VTER(3), VCEL(3)

PARAMETER ( PI = 3.14159265358979324D0 )
PARAMETER ( DEGRAD = PI / 180.D0 )

INTEGER PLANET_num
CHARACTER*3 PLANETS_name(9),PLANET_name

DATA STAR, OBSERV / 12 * 0.D0 /

DATA IYEAR, MONTH, IDAY, HOUR /
. 2020, 5, 1, 23.262097222222224D0 /

DATA LEAPS, UT1UTC, XP, YP /
. 37, -0.2425171D0, 0.075815D0,0.435887D0 /

DATA GLON, GLAT, HT / 117.577222D0, 40.395833D0, 900.D0 /

PLANETS_name= (/'MOO','MER','VEN','MAR',
.'JUP','SAT','URA','NEP','PLU'/)



* FORMAT STATEMENTS FOR OUTPUT

9010 FORMAT(1X, 3(F17.12, 8X))
9020 FORMAT(1X, 2(F15.6, 8X), 1F17.12)
9030 FORMAT(1X, 2(F17.12, 8X), 1D18.12)
9040 FORMAT(1X, F5.0,'d', 2X,F5.0,'m' 2X, 1F7.3,'s')
9050 FORMAT(1X, F5.0,'h', 2X,F5.0,'m' 2X, 1F9.5,'s')

WRITE ( *, * )
WRITE ( *, * ) 'NOVAS Calculations'
WRITE ( *, * ) '-------------------------'

* SETUP CALLS
* HIGH ACCURACY AND EQUINOX MODE ARE NOVAS DEFAULTS.
CALL HIACC
CALL EQINOX
IEARTH = IDSS ( 'EARTH' )

* WRITE OUT ASSUMED LONGITUDE, LATITUDE, HEIGHT (ITRS = WGS-84)
WRITE ( *, * )
WRITE ( *, * ) 'Geodetic location: '
WRITE ( *, 9010 ) GLON, GLAT, HT

* ESTABLISH TIME ARGUMENTS
CALL JULDAT ( IYEAR, MONTH, IDAY, HOUR, UTCJD )
TTJD = UTCJD + ( LEAPS + 32.184D0 ) / 86400.D0
UT1JD = UTCJD + UT1UTC / 86400.D0
DELTAT = 32.184D0 + LEAPS - UT1UTC
WRITE ( *, * )
WRITE ( *, * ) 'DATE and TIME'
WRITE ( *, * ) IYEAR,MONTH,IDAY,HOUR
WRITE ( *, * )
WRITE ( *, * ) 'TT and UT1 Julian Dates and Delta-T: '
WRITE ( *, 9020 ) TTJD, UT1JD, DELTAT

PLANET_num=1

DO 30 WHILE ( PLANET_num.LT.10 )

PLANET_name=PLANETS_name(PLANET_num)

NAME = IDSS ( PLANET_name )

WRITE ( *, * )
WRITE ( *, * )
WRITE ( *, * )
CALL APPLAN ( TTJD, NAME, IEARTH, RA, DEC, DIS )
WRITE ( *, * ) PLANET_name,'--------------------------'
WRITE ( *, * ) ' geocentric positions:'
WRITE ( *, 9030 ) RA, DEC, DIS
CALL thetaA(RA,rah,ram,ras)
WRITE(*,9050) rah,ram,ras
CALL thetaA(DEC,decd,decm,decs)
WRITE(*,9040) decd,decm,decs

CALL TPPLAN ( UT1JD, GLON, GLAT, HT, RAT, DECT, DIST )
WRITE ( *, * )
WRITE ( *, * ) ' geocentric positions:'
WRITE ( *, 9030 ) RAT, DECT, DIST
CALL thetaA(RAT,rath,ratm,rats)
WRITE ( *, 9050 ) rath,ratm,rats
CALL thetaA(DECT,dectd,dectm,dects)
WRITE ( *, 9040 ) dectd,dectm,dects




* POSITION OF THE PLANET IN LOCAL HORIZON COORDINATES

CALL ZDAZ ( UT1JD, XP, YP, GLON, GLAT, HT, RAT, DECT, 1,
. ZD, AZ, RAR, DECR )
WRITE ( *, * )
WRITE ( *, * ) ' zenith distance and azimuth:'
WRITE ( *, 9010 ) ZD, AZ
CALL thetaA(ZD,zdd,zdm,zds)
WRITE ( *, 9040 ) zdd,zdm,zds
CALL thetaA(AZ,azd,azm,azs)
WRITE ( *, 9040 ) azd,azm,azs

IF (zdd.GT.90.0) THEN
WRITE ( *, * )
WRITE ( *, * ) '!!!!!!',PLANET_name,' INVISIBLE'
ENDIF

PLANET_num = PLANET_num + 1
30 CONTINUE

END


SUBROUTINE thetaA ( RA, rah, ram, ras )
DOUBLE PRECISION RA
DOUBLE PRECISION rah,ram,ras

IF (RA.gt.0) THEN
rah = int(RA)
ram = ( RA - rah ) * 60
ras = ( ram - int(ram) )*60
ram = int(ram)
ENDIF

IF (RA.lt.0) THEN
RA = -RA
rah = int(RA)
ram = ( RA - rah ) * 60
ras = ( ram - int(ram) )*60
ram = int(ram)
rah = -rah
RA = -RA
ENDIF
RETURN
END

ps:

hexo公式总是容易出毛病。

1
$\mathrm{\mu_ \alpha=-354.45 mas/y  \quad \mu_ \alpha=+595.35 mas/y}$

这句话有问题吗? 就是识别不出来。