Play for Fun

0%

NOVAS实践

发表于 Valine：

关于使用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}$

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