返回

第八章 日月食原理

一、月食

月球被地影全部或部分遮掩的现象称为月食。当日月黄经相差180度(或附近)才可能发生月食。所以月食一般都发生在望日附近，即农历十五或十六日，这时地球运动至太阳和月球之间。 在月球轨道处，地球本影的直径比较大，是月球的2.5倍左右，月球常常被地球本影全部遮住，称为月全食。当月球始终只有部分被地球本影遮住时，这种天象称为月偏食。 如果月球进入半影区域，太阳的光也可以被遮掩掉一些，这种现象在天文上称为半影月食，不过月亮在半影区仍然很明亮，肉眼很难分辨。

月食过程基本描述

月食-月食的过程
半影食始：月球刚刚和半影区接触，这时肉眼觉察不到。
初亏：月食开始。月球由东缘慢慢进入地影，月球与地球本影第一次外切。
食既：月球的西边缘与地球本影的西边缘内切，月球刚好全部进入地球本影内。
食甚：月球的中心与地球本影的中心最近。
生光：月球东边缘与地球本影东边缘相内切，这时全食阶段结束。
复圆：月球的西边缘与地球本影东边缘相外切，这时月食全过程结束。
月球被食的程度叫“食分”，它等于食甚时月轮边缘深入地球本影最远距离与月球视经之比。
半影食终：月球离开半影，整个月食过程正式完结。

寿星月食算法

1、算出t时刻太阳、月亮的当日平春分点地心黄道视坐标(当然也可以旋转到赤道坐标)。只要日月使用统一的坐标就可以了，使用真分点或平分点的黄道或赤道坐标均可。

2、根据上图中的侧视图计算出t时刻本影及半影在距离地心d处的垂直于影轴的圆半径，d为此刻地月质心距。本影半径记为er，半影半径记为Er，月球视地心视半径记为mr。计算月食，实际上就是计算月亮被照的情况，不会因为观测者的位置改变而改变天象事实，所以与观测者的视差无关。计算方法如下：

设t时刻日地质心距为D，地球平均半径为r，太阳半径为R，那么由平面几何知识易得：er = r - (R-r)/D*d。天文计算时，习惯使用角度表达半径(视差法表示)，那么可以除以d，得到：

er = r/d - (R-r)/D，同理得到，Er = r/d + (R+r)/D

由于地球大气的影响，er和Er实际值比以上计算值要大百分2，所以结果取er=1.02*er，Er=1.02*Er；r取值为地球平均半径，即赤道半径的(1-1/298.25642/2)=0.99834，所以r=0.99834*6378.14千米；d的单位是千米；(R-r)/D实际上非常接近于：太阳视半径 - 太阳的地球赤道视差，所以R取值959.63角秒,r取8.794角秒,那么D使用天文单位。

t时刻月亮视半径的计算：mr = 358473400/d，单位是角秒。

3、根据上图中的侧视图计算出本影、半影中心坐标：t时刻太阳坐标相反的方向，即经度=太阳经度+180度，纬度=-太阳纬度，本(半)影中心到地心的距离 = t时刻地月质心距 = d千米。

4、t时刻影子中心与月球中心的距离计算：

t时刻月球坐标记为(J1,W1,d)，影子中心坐标为(J2,W2,d)，接下来计算出月球中心与地影中心的经度差和纬度差dJ和dW：dJ = J1-J2，dW = W1-W2

月食发生时，地影与月球相距很近，在地影与月影附近的区域可看作一个小平面，不必看做球面。该平面是天球在地影中心与月球中心连线的中心的天球切面。在这个小切面上创建一个临时坐标：x轴垂直于经线，y轴平行于纬线，原点建在地影中心。那么月心坐标为：x = dJ*cos( (W1+W2)/2 )，y = dW，单位统一使用弧度制。地影中心与月球中心角距离为pr = sqrt(x*x+y*y)，结果是弧度，sqrt是开方函数，

5、用图表即时表示出半影、本影、月球的位置关系（如上图中的正视图），并结合半径的数值显示，即可直观的找出初亏、食甚等关键数据。你先利用寿星万年历进行人工查找初亏、食甚，就可快速理解日月食的计算原理，那么只要你有一定的编程基础，就可迅速编制计算自动计算月食的计算机程序。寿星万年历中使用迭代法计算日月食甚。

二、日食

(一)月影轴、本影、伪本影、半影

(二)直角坐标与球面坐标的关系

(三)坐标变换
　球面转直角坐标llr2xyz()：
function llr2xyz(JW){
  var r=new Array(),J=JW[0],W=JW[1],R=JW[2];
  r[0]=R*cos(W)*cos(J);
  r[1]=R*cos(W)*sin(J);
  r[2]=R*sin(W);
  return r;
}
　直角转球面坐标xyz2llr()
function xyz2llr(xyz){
  var r=new Array(), x=xyz[0],y=xyz[1],z=xyz[2];
  r[2] = sqrt(x*x+y*y+z*z);
  r[1] = asin(z/r[2]);
  r[0] = rad2mrad( atan2(y,x) );
  return r;
}
　球面坐标旋转函数llrConv()：
function llrConv(A,E){
  var B=new Array(),J=A[0],W=A[1];
  B[0]=rad2mrad( atan2(sin(J)*cos(E) - tan(W)*sin(E), cos(J)) );
  B[1]=asin ( cos(E)*sin(W) + sin(E)*cos(W)*sin(J) );
  B[2]=A[2];
  return B;
}
　说明：rad2mrad()函数的作用是将角度转换到0到2PI
　用法：B = llrConv(A,E)
　功能：由A道坐标转到B道坐标，E为A道与B道的夹角。
　A与B的起算点为A与B的交点，交点有升交点和降交点之分，起算点由用户任选，但应注意：在A道上沿经度测量方向走，经过起算点时，如果纬度增加，则表示起算点为升交点，这时E取正，反之为降点E取负。
　例1：A为黄道坐标，转换后B为赤道坐标，如果经度起算点为春风点，那么沿A道(即黄道)，经过春风点时，黄纬是增加的，所以E取正，E是黄赤交角。
　例2：A为赤道坐标，转换后B为黄道坐标，如果经度起算点为春风点，那么沿A道(即赤道)，经过春风点时，赤纬是增加的，所以E取负，E是黄赤交角。

(四)贝塞尔坐标系

贝塞尔坐标系的定义：

z轴经过地心且平行于月影轴，往太阳那个方向为正。垂直于月影轴且经过地心的平面为基面。基面与赤道面的交线为x轴，xyz三轴构成右手坐标系。如果规定x轴正方向的赤经比z轴赤经多90度而不是-90度，这样y轴总是在北半球。

设z轴上某一点记坐B向量，写为直角坐标形式记作B[0,1,2]，jd时刻太阳赤道坐标为S[0,1,2]，月亮赤道坐标为M[0,1,2]，存在关系式B = S - M。应注意，这个关系式是一个矢量运算关系。当S、M为球面坐标时，请先将它们转为直角坐标，然后再运算。最后得到的B还可以转为球面坐标。因此，在球面坐标系中，z轴向量的赤道坐标表示：B = xyz2llr( llr2xyz(S) - llr2xyz(M) )。实际上，B的方向十分接近太阳，有中心线日食发生期间，B与S的角度差不到10角秒。

当z轴向量B求出后，贝赤交角立即得到：ζ= 90-B[1]，同时得到A点的赤经(起算点赤经)：υ = 90+B[0]。

已知天体赤道坐标P[]求其贝塞尔坐标：
　(1)先把P经度转为A点起算的经度P[0] = P[0]-υ
　(2)然后旋转P = llrConv(P,-ζ)。
　(3)因为在赤道上沿经度测量方向走，经过起算点A为降交点，所以这里ζ取负。通常，天体的贝塞尔坐标以直角坐标形式表达，所以还应计算P = llr2xyz(P)。
　例：已知月亮赤道坐标，求它的贝塞尔坐标。按以上三步骤即可。

已知天体贝塞尔坐标P[]求基赤道坐标：
　(1)P = xyz2llr(P)。
　(2)P = llrConv(P,ζ)
　(3)P[0] = P[0]+υ

用直角坐标进行贝、赤坐标变换

xyz是贝塞尔直角坐标系，XYZ是贝赤交线起算的赤道直角坐标系。x轴与X轴重合，y轴向量(0,P,Q),z轴向量(0,-Q,P)，其中P = cosζ，Q = sinζ

贝赤交线与贝黄交线的夹角

应用三垂线定理易证得自转轴总是在yz平面。由于太阳距我们很远，月球则要近得多，所以月球的移动对月影轴方向的影响是轻微的，月影轴的方向主要有太阳决定，基本与太阳相同，因此，如果太阳赤纬为负时，上图中y应画在北极左边。这就是说随着太阳运转，赤纬变化，北极在y轴上左右移动。

由于太阳黄纬小于2角秒且太阳视差小于10角秒，所以基面基本垂直于黄道面，当月影轴经过地球，基面与黄轴的夹角不超过10角秒。虽然，实际计算不推荐看成垂直，但建立这样一种概念有利于坐标变换。贝黄交角：约为90度。

下图则给出贝塞尔坐标、黄道、赤道、基面等的示意图：

图中∠BOA可按以下方法得到：由黄赤变换公式sinδ = sinβ cosε + cosβ sinε sinλ，A点赤纬是δ=0，B点黄经λ=υ，代入得tanβ = -tanε sinλ = -tanε sin(υ)。贝圆几乎垂直于黄道，所以A的黄纬β近似为∠BOA，则有tan∠BOA = -tanε sinυ

求贝圆的短半径：

从贝塞尔基面的定义知，贝塞尔基面与椭球形地球的交线为椭圆，下文称贝圆。贝圆的长半径在x轴，长度为地球赤道半径a。短半径在y轴，如何求其大小呢？我们引入下图，这是地球的一个剖面图，贝塞尔坐标的y轴在该落在剖面上，OY正是贝圆的半短径b₂，今设赤道半径为a，南北极半径为b。

依上图建立椭圆方程 (x/a)²+(y/b)²=1，把x = -r*cosζ，y = r*sinζ代入得：
　b₂ = r = b/sqrt(1-c*cos²ζ)，式中c = 1-b²/a²
　前面已讲到ζ的计算，约等于“90-太阳赤纬”，因此c*cos²ζ是个小量，因此上式可以简化为更简洁的公式：
　因为cosζ≈sinδ≈sinε*sinλ，(ε为黄赤交角，λ为太阳黄经，δ为太阳赤纬，因z轴与太阳几乎同方向，所以约等号成立)，所以
　b₂ = r ≈ b(1+c*sin²λ) ，式中c=(1-b²/a²)*sin²ε/2
　a取值单位1，地球偏率f=1/298.25642，b = 1-f，ε=84381.448角秒，则有
　b₂ = 0.99664718(1+0.0005296sin²λ)
　当影轴刚好扫过地球边缘，为了判断有无中心线，b₂通常它的最大可能值0.9972即可，不必取均值，其实南极大陆平均海拨也有2千米左右。

在贝塞尔坐标系中求解地球表面某点的速度：

地理坐标与贝塞尔坐标的转换：

我们在一个实际问题中求解。贝塞尔坐标中的点P(x,y,z)落在地球表面，z轴赤经赤纬是(υ,90-ζ)，格林尼治恒星时为g，地球赤道半径为a与极半径为b，若z未知，试求z的值。

空间直线与地球的交点

以下示例程序计算直线与地球的交点，直线经过(x1,y1,z1)和(x2,y2,z2)两点，e为极赤径比，r为赤半径：
function lineEll(x1,y1,z1, x2,y2,z2, e,r){ //求空间两点连线与地球的交点(靠近点x1的交点)
  var dx=x2-x1, dy=y2-y1, dz=z2-z1, e2=e*e, A,B,C,D,R,t, p=new Object();
  A = dx*dx + dy*dy + dz*dz/e2;
  B = x1*dx + y1*dy + z1*dz/e2;
  C = x1*x1 + y1*y1 + z1*z1/e2 - r*r;
  p.D = B*B-A*C; if(p.D<0) return p; //判别式小于0无解
  D = sqrt(p.D); if(B<0) D = -D; //只求靠近x1的交点
  t = (-B+D)/A;
  p.x=x1+dx*t, p.y=y1+dy*t, p.z=z1+dz*t;
  R = sqrt(dx*dx + dy*dy + dz*dz);
  p.R1 = R*abs(t), p.R2 = R*(t-1); //单位是千米
  return p;
}

求某一地标(含海拔)的贝塞尔坐标

求影轴上某点对应的本影和半影半径

影子中心相对于观测点的速度

求某时刻的中心点

半影与贝圆的交点(求两圆交点,贝圆用地球平均半径)

食界图的数值表达