Python in Astronomy (5.数据的处理与可视化（二））

2.3日弧

        从地球上观察者的角度来看，天体在天球表面的运动是沿着地平线上方的弧线，这条在地平线上方的弧线也被称为日弧，对于观测者而言，天体只有在这条日弧上才能被观测到，如果天体进入到地平线以下便不能被该观测者观测到。天体在日弧上的位置信息可以由一个被称为时角的物理量来描述。一天24个小时的时角对应于一个平行于天球赤道的360度圆弧，正如下图中完整的红色圆圈。基于此，时角可以等价的用度数或者弧度表示。

        celestial meridian 也就是所谓的天球子午圈，在时角坐标系中，也被称为观测者子午圈。celestial equator 就是所谓的天赤道，天体达到地平线上最高的高度时，对应图中的就是周日弧与天球子午圈相交的一点。根据定义，当天体到达这点时，时角为零。

        图中的对应的是北天极与地平面之间的距角，由于地球的自转轴垂直于天赤道，这里所谓的自转轴就是北天极与南天极的连线。赤道面与地平面之间的倾角可以表示为，天体在子午线上处于最高点时，最大距角可以表示为，这里的即天体的赤纬。在这个时角坐标系中，天体在时刻开始上升，并在时达到最高值，然后天体开始落下。当天体从地平线上落下时，正好有。

        根据球面三角的知识，我们可以求出：

其中，是天体的赤纬，是观测者在地球上的位置的纬度。因此，就是天体，原则上在天空中可见的时间，它也被称为日弧的长度。

        例如，让我们考虑猎户座中的参宿四恒星，它是一颗红巨星，是天空中最亮的恒星之一。在astropy.coordinates库的帮助下很容易找到它的位置。

from astropy.coordinates import SkyCoord, EarthLocation
betelgeuse = SkyCoord.from_name('Betelgeuse')
print(betelgeuse)

然而，当第一次遇到这个输出时，可能需要一些解读：

这告诉我们，这个被称为参宿四的天体，它的赤经和赤纬分别为88.793度和7.407度。在第2行中定义的变量Betelgeuse，参宿四不仅仅代表一个天体，它是一个Python对象，更具体地说，它是一个属于SkyCoord类的对象。对象的dec属性允许我们直接引用赤纬。

delta = betelgeuse.dec
print(delta)

赤纬的输出格式分别为(d)度、弧分（m）和弧秒（s），这是天文学中表示角坐标的首选格式：

即，，其中，。假设，我们想确定在汉堡Hamburg天文台可以看到的参宿四的日弧长度（），除了要知道恒星的赤纬外，我们还需要观测者的位置。与SkyCoord对应的天体坐标的对应物是EarthLocation（也在第1行中导入了），它允许我们设置地球上一个位置的地理纬度和经度：

import astropy.units as u

# 观测者的地理位置
obs = EarthLocation(lat=53*u.deg+28*u.arcmin+49*u.arcsec, lon=10*u.deg+14*u.arcmin+23*u.arcsec)

# 获取纬度
phi = obs.lat

在上面的表达式中，u.deg相当于1度，u.arcmin 相当于1分， u.arcsec 相当于1秒。astropy中的units模块制定执行物理与天文单位系统间的转化的规则。接下来，我们利用math模块中的三角函数来计算h：

import math
h = math.acos(-math.tan(delta.radian) * math.tan(phi.radian))

在应用三角函数进行数学运算之前，有必要将角度δ和φ转换为弧度。对于天体坐标，我们所需要做的就是使用角度的radian属性。为了获得以小时为单位的T，我们需要记住的是，一个360度的角对应一个恒星日，就是地球完整的一个自转，这比太阳日短四分钟。这个是地球轨道运动的结果。units模块使转换变得容易。

T =(math.degrees(2*h)/360)*u.sday
print("T = {:.2f}".format(T.to(u.h)))

第一个T是以恒星日（u. sday）定义的，它相当于24小时或360度。然后，我们通过将该方法应用于第2行中的to（）来转换为太阳时（u.h）。结果是

如果没有太阳，参宿四可以在天文台看到13小时。当然，这颗星星只有在这段时间和夜晚的重叠期间才能看到，这取决于观测的日期。我们将在下一节回到这个问题。

        太阳的日弧在我们的日常生活中起着核心作用，因为它决定了白天的时间。我们在上一节中引入了太阳的的近似值，我们可以计算出一年中白天的长度是如何变化的。首先，我们需要计算N从0到364天太阳的。使用NumPy模块，这很容易做到。.在下面的示例中，表达式np.arange（365）表示用从0开始到小于365的最大数（364）的整数序列填充数组。除此之外，代码的工作方式类似于上节中基于NumPy的分点和至点的计算。

import numpy as np

N = np.arange(365) # 包含元素0、1、2、…、364的数组
omega =2*math.pi/365.24 # 地球的角速度，单位是rad/day
ecl = math.radians(23.44) # 黄道倾角

# 计算一年中所有日子的太阳赤纬
delta =-np.arcsin(math.sin(ecl)* np.cos(omega*(N+10)))

现在我们可以使用numpy模块中的函数计算数组delta中所有值的日长T（与参宿四的代码示例比较并检查进行了哪些更改）：

# 以太阳时计算日长
h = np.arccos(-np.tan(delta)* math.tan(phi.radian))
T =(np.degrees(2*h)/360) * u.sday.to(u.h)

在这里，仍然是上面定义的天文台的纬度。当然，现在T是一个具有365个元素的数组。由于我们想要以太阳时为单位的日长，乘以u.sday.to（u.h），便可得到以太阳时为单位的恒星日的长度。在处理大型数据集（即多个值）时，大多数情况下最好提取一些统计数据或以图形形式显示数据。为了显示日长的年变化，绘制T与N的关系图。我们可以使用Python库中的matplotlib（参见matplotlib.org），这个库提供了一个名为pyplot的模块，用于在数组中绘制数据：

import matplotlib.pyplot as plt

plt.plot(N, T)
plt.xlabel("Day")
plt.ylabel("Day length [hr]")
plt.savefig("day_length.pdf")

函数plt.plot（）（其中plt是matplotlib.pyplot的常用别名）生成一个显示数组N（x轴）和T（y轴）给出的数据点（x，y）的图像。默认情况下，点通过线连接以生成连续图形。上面的代码还给图像
添加了轴标签，也将绘图保存到PDF格式的day_length.pdf文件中。文件的位置取决于启动Python的目录，如果想将绘图存储在其他地方，可以指定完整路径。如

Path = 'D:\pythonfiles\GRID\G02_tte_2101211840_2101211901_v03_00.fits'

正如预期的那样，1月的日长很短，在第一个至日（N = 171）达到最大值，然后减少，直到第二个至日达到N = 355。最小和最大日长可以在NumPy数组的min（）和max（）方法的帮助下推断。

print("Minimum day length = {:5.2f} h".format(T.min()))
print("Maximum day length = {:5.2f} h".format(T.max()))

最小日长和最大日长之间的差异随纬度的增加而增加。除了极圈以外，日长在0 ~ 24 h之间变化。另一个例子是朗伊尔城Longyearbyen ，位于挪威北部的斯匹次卑尔根岛上，。由于我们只需要知道一个精度在角分内的（精度受到近似的十进制数的限制），使用EarthLocation来完整地描述地理位置会做得太过了。相反，我们简单地使用1 =（1/60）π，转换为弧度：

phi = math.radians(78+13/60) # latitude of Longyearbyen

h = np.arccos(-np.tan(delta)*math.tan(phi))
T =(np.degrees(2*h)/360) * u.sday.to(u.h)

当你执行这段代码时，会得到

事实证明，Python能够处理这个问题，但我们仍然应该尝试理解这里的错误。 

请记住，太阳赤纬的范围是。又这个公式。

当时，也就是位置处于极区，则右侧小于- 1或大于1。如果◦，其参数的绝对值大于1，则反余弦未定义。这对应于极夜或极日，在此期间太阳从不升起或落下。当上式的右侧落在区间[- 1,1]之外时，可以通过设置等于±1来解决这个问题。使用Numpy模块，这很容易通过np.clip()函数实现：

tmp = np.clip(-np.tan(delta)*math.tan(phi), -1.0, 1.0)
h = np.arccos(tmp)
T =(np.degrees(2*h)/360) * u.sday.to(u.h)

正如预期的那样，你会发现朗伊尔城的人们在冬天必须面对极夜(整天黑暗)，而在夏天则必须面对持续24小时的极日。

        在计算了汉堡Hamburg和朗伊尔城Longyearbyen的日长之后，用一个图来比较这两个图将是有益的。显然，我们同时需要两个位置的数据，在上面的示例中，汉堡Hamburg的数据被朗伊尔城Longyearbyen的计算覆盖了。解决此问题的一种直接方法是使用不同命名的数组，但在Python中有一种更优雅、更方便的替代方法。数据除了储存在数组中之外，还可以储存在Python的字典中。与传统意义上的字典类似，字典中的数据项是通过关键字而不是数字索引引用的，字典的数据项可以是任何东西，包括数组。让我们建立一个将位置与其纬度相关联的字典：

phi ={ 'Hamburg': obs.lat.radian,
       'Longyearbyen' : math.radians(78 + 13/60) }

在Python中，字典是由用花括号括起来的关键字和项对定义的。每个关键字必须是字符串(例如' Hamburg ')，用冒号与相应的项(在Hamburg的情况下为表达式obs.lat.radian)分隔，然后后面跟着逗号。我们可以通过下面的代码访问单个项目，如Hamburg

print(phi['Hamburg'])

以弧度打印汉堡的纬度。语法类似于访问数组元素，不同之处在于使用关键字而不是索引。向字典中添加新条目是很容易的。例如，下面的代码将纽约和曼谷添加到我们的字典中：

phi['New York']= math.radians(40 + 43/60)
phi['Bangkok']= math.radians(13 + 45/60)

实际上，字典现在有四个条目，可以通过打印len(phi)来检查。下面的代码将打印所有项目，计算每个地点的日长，并将它们合并到一个图中:

for key in phi:
    print(key + ": {:.2f} deg".format(math.degrees(phi[key])))
    h = np.arccos(np.clip(-np.tan(delta)*math.tan(phi[key]),-1.0, 1.0))
    T =(np.degrees(2*h)/360) * u.sday.to(u.h)    
    plt.plot(N, T, label=key)

plt.xlabel("Day")
plt.xlim(0,364)
plt.ylabel("Day length [hr]")
plt.ylim(0,24)
plt.legend(loc='upper right')
plt.savefig("daylength.pdf")

如下图所示。正如预期的那样，朗伊尔城的一天长度在0到24小时之间变化，而位于热带地区的曼谷全年的一天长度约为12小时。从冬天到夏天，温带地区的白昼长度增加了几个小时。

我们还可以在字典中添加一个南极附近的纬度，如

phi['Antarctic']= math.radians(-78 + 43/60)

由上面的事例，我们可以容易得到世界上任何纬度的日长。

参考文献

[1]. W. Schmidt and M. Völschow, Numerical Python in Astronomy and Astrophysics , Undergraduate Lecture Notes in Physics,