本文章修改自YoungJune的博客

实践题目

题目一

天体星表位置：ICRS[α,δ]=14h34m16s.81183,−12°31′10′′.396ICRS[\alpha,\delta]=14^h34^m16^s.81183,-12\degree31^\prime10^{\prime\prime}.396
自行：μα=−354.45mas/y,μδ=+595.35mas/y\mu_\alpha=-354.45mas/y,\mu_\delta=+595.35mas/y
视差：0′′.164990^{\prime\prime}.16499，径向速度（退行速度）=0km/s=0km/s
观测时刻2023年5月1日，23h15m43s.55UTC2023年5月1日，23^h15^m43^s.55UTC
地球指向参数请在IERS网站自查
求：

  1. 该星地心视位置
  2. 测站坐标：W=70°44′11′′.560,S=30°14′26′′.731W=70\degree44^\prime11^{\prime\prime}.560,S=30\degree14^\prime26^{\prime\prime}.731，参考椭球面高度2378m2378m，求该星站心位置、地平坐标
    注：赤经以时分秒表示，赤纬以度分秒表示，精确到毫角秒

题目二

观测时刻2023年5月1日，23h15m43s.55UTC2023年5月1日，23^h15^m43^s.55UTC
地球指向参数请在IERS网站自查
求：

  1. 月球和八大行星的视位置
  2. 若观测站是慕士塔格站，求对应观测时刻月球和八大行星站心位置、地平坐标，并指出其中哪些天体可观测
    注：赤经以时分秒表示，赤纬以度分秒表示，精确到毫角秒

实践准备

查询地球指向参数

首先进入IERS官网
然后找到上方菜单栏中的"Data/Products/Tools"一栏，进入Earth Orientation Data页面
如果是查询近期数据，请找到 “Daily EOP Data Files” ；如果是查询稍早的数据，可以在 “Standard EOP data files” 中寻找
根据检索到的数据，我们可以确定2023年5月1日的地球指向参数为：DX=0.027600,DY=0.473427,UT1−UTC=−0.0342596DX=0.027600,DY=0.473427,UT1-UTC=-0.0342596

所需程序准备

前往NOVAS发布页即可下载三种版本的程序包（C,Fortran,Python）
本文主要阐述Fortran程序的使用
对于下载的文件中，最重要的为：

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

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

  • jplsubs.f

程序参数设置

观测时间设置

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

地球指向参数设置

1
2
 DATA LEAPS, UT1UTC, XP, YP /
.     37, -0.0342596D0, 0.027600D0,0.473427D0 /

观测地点设置

题目一：

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

题目二（慕士塔格测站）：

1
DATA GLON, GLAT, HT / 74.896667D0, 38.329722D0, 4520.D0 /

天体位置设置

1
2
3
4
5
6
7
8
9
10
!赤经赤纬(hours&degrees)
RA2000 = 14.571336619444445D0
DC2000 = -12.519554583333333D0
!自行数据(milliarcseconds per year)
PMRA   =  -354.45D0
PMDEC  = +595.35D0
!视差(milliarcseconds)
PARX   =   164.99D0
!径向速度(kilometers per second)
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

运行程序

CMD命令行中在程序目录下运行以下代码即可查看结果

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
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 /
     .     2023, 5, 1, 23.262097222222224D0 /
     
      DATA LEAPS, UT1UTC, XP, YP /
     .     37, -0.0342596D0, 0.027600D0,0.473427D0 /

      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
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
       2023           5           1   23.262097222222224

TT and UT1 Julian Dates and Delta-T:
 2460066.470055         2460066.469254          69.218259600000

The star geocentric positions:
  14.592467433622         -12.619216610229
  14.h    35.m   32.88276s
 -12.d    37.m    9.180s

The star topocentric positions:
  14.592468245527         -12.619200046096
  14.h    35.m   32.88568s
 -12.d    37.m    9.120s

The Star zenith distance and azimuth:
  76.011253730290          96.653658986587
  76.d     0.m   40.513s
  96.d    39.m   13.172s

题目二

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 /
     .     2023, 5, 1, 23.262097222222224D0 /
     
      DATA LEAPS, UT1UTC, XP, YP /
     .     37, -0.0342596D0, 0.027600D0,0.473427D0 /

      DATA GLON, GLAT, HT / 74.896667D0, 38.329722D0, 4520.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
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
176
177
178
179
180
181
182
183
184
185
NOVAS Calculations
-------------------------

Geodetic location:
  74.896667000000          38.329722000000        4520.000000000000

DATE and TIME
       2023           5           1   23.262097222222224

TT and UT1 Julian Dates and Delta-T:
 2460066.470055         2460066.469254          69.218259600000



MOO--------------------------
 geocentric positions:
  11.867927667533           4.374562203087        0.264943349375D-02
  11.h    52.m    4.53960s
   4.d    22.m   28.424s

 geocentric positions:
  11.821376351285           3.793633520953        0.265664960559D-02
  11.h    49.m   16.95486s
   3.d    47.m   37.081s

 zenith distance and azimuth:
 100.160856819185         283.142744082747
 100.d     9.m   39.085s
 283.d     8.m   33.879s

!!!!!!MOO   INVISIBLE



MER--------------------------
 geocentric positions:
   2.578536927061          15.900518772120        0.563556414096D+00
   2.h    34.m   42.73294s
  15.d    54.m    1.868s

 geocentric positions:
   2.578748500124          15.897530481371        0.563562942754D+00
   2.h    34.m   43.49460s
  15.d    53.m   51.110s

 zenith distance and azimuth:
  98.727010767764          61.664636731864
  98.d    43.m   37.239s
  61.d    39.m   52.692s

!!!!!!MER   INVISIBLE



VEN--------------------------
 geocentric positions:
   5.538268782949          25.740064814090        0.974181104028D+00
   5.h    32.m   17.76762s
  25.d    44.m   24.233s

 geocentric positions:
   5.538314712842          25.737857697822        0.974197974572D+00
   5.h    32.m   17.93297s
  25.d    44.m   16.288s

 zenith distance and azimuth:
 113.147641256223          19.942632549962
 113.d     8.m   51.509s
  19.d    56.m   33.477s

!!!!!!VEN   INVISIBLE



MAR--------------------------
 geocentric positions:
   7.427469111192          23.753297018925        0.174725237665D+01
   7.h    25.m   38.88880s
  23.d    45.m   11.869s

 geocentric positions:
   7.427452978442          23.752072604868        0.174727215423D+01
   7.h    25.m   38.83072s
  23.d    45.m    7.461s

 zenith distance and azimuth:
 117.470433469730         351.784421625528
 117.d    28.m   13.560s
 351.d    47.m    3.918s

!!!!!!MAR   INVISIBLE



JUP--------------------------
 geocentric positions:
   1.668595019943           9.254589341921        0.592146579491D+01
   1.h    40.m    6.94207s
   9.d    15.m   16.522s

 geocentric positions:
   1.668615524196           9.254316594161        0.592146820476D+01
   1.h    40.m    7.01589s
   9.d    15.m   15.540s

 zenith distance and azimuth:
  93.191702488666          75.557408466648
  93.d    11.m   30.129s
  75.d    33.m   26.670s

!!!!!!JUP   INVISIBLE



SAT--------------------------
 geocentric positions:
  22.520919121747         -10.820654369663        0.101733733204D+02
  22.h    31.m   15.30884s
 -10.d    49.m   14.356s

 geocentric positions:
  22.520932292854         -10.820809878686        0.101733590842D+02
  22.h    31.m   15.35625s
 -10.d    49.m   14.916s

 zenith distance and azimuth:
  70.551901948775         122.173063036610
  70.d    33.m    6.847s
 122.d    10.m   23.027s



URA--------------------------
 geocentric positions:
   3.074245732695          17.025077571673        0.206508473512D+02
   3.h     4.m   27.28464s
  17.d     1.m   30.279s

 geocentric positions:
   3.074248556030          17.024975915445        0.206508569424D+02
   3.h     4.m   27.29480s
  17.d     1.m   29.913s

 zenith distance and azimuth:
 102.899636066746          55.669469421038
 102.d    53.m   58.690s
  55.d    40.m   10.090s

!!!!!!URA   INVISIBLE



NEP--------------------------
 geocentric positions:
  23.833155397077          -2.377021552298        0.306191000461D+02
  23.h    49.m   59.35943s
  -2.d    22.m   37.278s

 geocentric positions:
  23.833160697695          -2.377068639681        0.306190919369D+02
  23.h    49.m   59.37851s
  -2.d    22.m   37.447s

 zenith distance and azimuth:
  79.015367069923         101.922500235627
  79.d     0.m   55.321s
 101.d    55.m   21.001s



PLU--------------------------
 geocentric positions:
  20.209334272430         -22.516544810705        0.345529677071D+02
  20.h    12.m   33.60338s
 -22.d    30.m   59.561s

 geocentric positions:
  20.209340369760         -22.516596092512        0.345529486383D+02
  20.h    12.m   33.62533s
 -22.d    30.m   59.746s

 zenith distance and azimuth:
  63.570373943039         159.664217809630
  63.d    34.m   13.346s
 159.d    39.m   51.184s