一、前提条件：

假定所有的几何条件都是理想的：

1 地球是理想球体

2 地球轨道是理想圆

3 自转轴与黄道面夹角始终不变

 

二、参数定义

1 参照坐标：

O点：地心

x轴：自转轴在黄道面的投影方向的垂直方向，太阳方向

y轴：自转轴在黄道面的投影方向，北极方向

z轴：与黄道面的垂直方向，北半球方向

xy平面即为黄道面

2 位置角度

α：地心到太阳中心连线与x轴的夹角

β：自转轴与黄道面夹角

γ：某地纬度值

θ：某地经度值

3 切平面

过地心且与地心太阳中心连线相垂直的平面

4 纬度面

过同一纬度的平面与地球相交的圆面

纬度面在xy平面的投影为椭圆，圆心为O’

5 半径

R：地球半径

r：纬度面的半径。r=R*cosγ

 

三、原理

能知道某地的日照时间，即，切平面所分割的纬度面的弧度，可推知某地太阳升起的时刻。

 

四、推导

有两种推导方法，已验证其结果一致。示意图如下

http://s5/middle/6db19107tab4e9db29824&690

1 平面解析几何方法

方法：

纬度面投影椭圆与切平面投影线在xy平面的交点为A，B。

通过求解纬度面的投影椭圆方程与切平面投影线方程的方程组，得到A与B的横坐标值，x1、x2。从图知，日照时间=x1对应的弧度 + x2对应的弧度。

 

求解：

纬度面投影椭圆方程：  x²+(y-R*cosβ*sinγ)²/sin²β=R²*cos²γ

切平面交线方程：      y=x*ctgα

将y代入椭圆方程，可得一元二次方程：

x²*(ctg²α+sin²β)-2*x*R*ctgα*cosβ*sinγ +R²*(cos²β*sin²γ-sin²β*cos²γ)=0

求得：

x1=R*(ctgα*cosβ*sinγ+√(ctg²α*sin²β*cos²γ+sin²β*(sin²β-sin²γ)))/(ctg²α+sin²β)

x2=R*(ctgα*cosβ*sinγ-√(ctg²α*sin²β*cos²γ+sin²β*(sin²β-sin²γ)))/(ctg²α+sin²β)

cosΔ1=x1/R*cosγ

cosΔ2=x2/R*cosγ

 

2 立体解析几何方法

方法：

纬度面与切平面的交点为A，B。

通过求解纬度面与切平面交线方程的方程组，得到A与B的距离，即纬度面内AB的弦长，再根据三角形内角和公式求得AB对应弧度。

假定A，B坐标为：A(x1,y1,z1)，B(x2,y2,z2)，AB的距离为D，则：

D=√((x1-x2)²+(y1-y2)²+(z1-z2)²)

 

求解：

地球面方程：  x²+y²+z²=R²

纬度面方程：  y*cosβ+z*sinβ=R*sinγ

切平面方程：  y=x*ctgα

中间推导过程比较复杂，省略。

求得：

D = 2*R*√((sin²β+cos²γ*ctg²α-sin²γ)/(ctg²α+sin²β))

根据弦长的弧度公式：cosΔ =1-a²/(2*r²)，求得：

cosΔ=1-2*(sin²β+cos²γ*ctg²α-sin²γ)/((ctg²α+sin²β)*cos²γ)

 

五、最终结果

假定太阳升起的时刻为T，某地所在时区值为t，得到：

T=12*Δ/(2*π)+24*θ/360-t

T=6*Δ/π+θ/12-t

 

六、验证

下面贴出VBA的验证程序，两种方法的输出结果一致，但与授时中心的结果有差异。春分、夏至差异小些，秋分、冬至差异比较大。应该是跟近日点，远日点与春分点的位置有关。

 

平面解析几何方法的VBA的验证程序：

http://s2/middle/6db19107tab4e742d29f1&690

 

立体解析几何方法的VBA的验证程序：

http://s10/middle/6db19107t7787d8d1ccc9&690

七、后记

这几天没事浏览博客，发现这里 太阳几点升起？这个有点难 的计算公式有问题，且文字说明比较乱,自己看着都费劲。花了点时间重新推导了一下，再整理文字，配上示意图。应该可以看明白了。