已发布4.68

1、利用浪淘沙的算法，简化直线与椭球交点的计算

许兄过奖了，我的数学水平很一般。

 

我开始也在想着，把空间直线投影在赤道平面（即XY平面），及二个子午面（XZ平面，YZ平面）。后来，我用作图法试了一下，觉得很难解出交点坐标（X、Y、Z）。只好放弃了。

 

翻《数学手册》有现成的过空间二点的直线方程。但不是参数形式。

再参考手册上的参数形式直线方程，就得到我上面的直线方程了。

 

所以说我也只会翻书，不会自己想问题。

 

以后许兄要优化算法，可能我也帮不上大忙。当然，有什么想法时，提出来，大家互相探讨。

我想论坛里的许多朋友看到了，或许会提供一些好的思路。这样对许兄的程序优化，可能会有一点帮助吧。

我试着验证了几个日食数据。

新版与老版的计算结果是一致的。

还有一个问题：

在求解日出日没的初亏复圆线的时候，相当于求各时刻半影与贝塞尔基面上椭圆的交点。贝塞尔基面是经过地心且垂直于月影轴的平面，它与地球相交为椭圆。

这变成求解椭圆x^2/a^2+y^2/b^2=1与圆(x-x0)^2+(y-y0)^2的两个交点。《日月食计算概要》解半影问题，是把地球看成正球，球半径可考虑使用地球平均半径，这样，这个问题变成圆和圆的交点。不过，我还是想求椭圆与圆的交点，可是变成4次方程，有四解，而且解不出来。不知有没有简易的，且精度足够的近似方法。

 

再有一个问题：

 

如果考虑大气折射的影响，交点应在基面以下(34角分/(180*60角分))*3.14*地球半径=63公里，即z轴坐标为负63公里，不知如此考虑问题是否合适。现在还产生一个问题，日出日没的初亏复圆线是否有必要考虑这63公里的影响

在不考虑大气折射时，

第一个问题是否可以考虑用“晨昏圈（日夜分割线）”，即过地心且垂直于地日连线的平面与地球表面的交线。

考虑地球是旋转椭球体，这个晨昏圈方程不是特别简洁。

 

由于月影半径很小，只有晨昏圈上的一小段能够看到日食。

 

假若定义晨昏圈所在平面为晨昏面（即过地心垂直于地日连线的平面），则月亮的半影锥体在晨昏面上的投影为半影椭圆。

 

实际上需要求解的就是半影椭圆与晨昏圈的交点，至多只有二个交点。

交点处既是初亏，也正好是日出（或日落）。

交点之内一小段晨昏圈，实际初亏发生在日出后的一段时间。

 

解这二个椭圆的交点，似乎可以把其中一个椭圆（比如晨昏圈）表示为参数方程。

x=a * cost

y=b * sint

另一个椭圆用普通方程来表示，这个半影椭圆的长轴不一定与X轴平行或垂直了（方程式估计很复杂的）。

这个方程最终是关于 t 的二次方程，不会是四次方程。

 

但这个方程有正余弦的表达式，求解还是比较难的。

所以建议用“数值法”。而用“解析法”相对要难一些了。

用迭代法可能求得需要的解。

 

具体过程我没推导过。不知我这个想法能否管用。

 

注：晨昏圈不一定垂直于赤道面。所以参数方程式中的xy并非地球坐标的XY。还得进行坐标旋转的。

我想想就头痛了。

 

*******************

第二问题：

如果考虑大气折射，那么这个晨昏圈就非常复杂了，绝不会是简单的椭圆了。

 

所以解题时，似乎可先用理想曲线（无大气折射）求出一个解，然后加入大气折射进行修正，看看误差有多大。

 

1、我刚才认直考虑了一下，地球非正圆造成的影响比大气折射要小，以影响交点的纬度为主。半影刚切入地球时，交点的误比较大，但时间很短，因为半影很快进入进球，交点明确清析。

所以不要用椭圆，直接用正圆就可以了，反正这期间的误差还有其它因素，如大气折射。

2、楼上说到的晨昏圈，可能计算起来有难度。其实有简便的方法：

设两圆的交点为(x1,y1,0)，(x2,y2,0)，这就日出(或日没)时刻看见初亏的地方，当然，半影进入地球过半后，一般为日出(或日没)看见复圆的地方。

由于大气折射，在交点下方63公里出就已看到初亏或复圆。所以交点改正为(x1,y1,-63)和(x2,y2,-63)即可

最后将这两个坐标旋转到赤道坐标，再转换为地理坐标。

把每一时刻的两圆交点在地标描点，并连线，就得到日出没的初亏复圆线。

求两正圆的交点就容易一些。

第一个圆的圆心为(0,0)，半径为R

第二个圆的圆心为(x0,y0),半径为r

 

求交点。

 

两圆心距离d=sqrt(x0^2+y0^2)

第一圆心为O，第二圆心为A，交点为B，A与x轴的夹角为 θ，则cosθ=x0/d，sinθ=y0/d

角BOA为α，则利用余弦定理得cosα=(R^2+d^2-r^2)/(2Rd)，sinα=+-sqrt(1-cosα^2)

那么在以OA为x轴的直角坐标系中，B的坐标为(x',y')=(Rcosα，Rsinα)

将坐标轴旋轴-θ角度，则回到原来的坐标系中，则B坐标为

x = x'cosθ - y'sinθ = Rcosα cosθ-Rsinα sinθ

y = x'sinθ + y'cosθ = Rcosα sinθ-Rsinα cosθ

 

x = x'cosθ - y'sinθ = Rcosα cosθ-Rsinα sinθ

这里,“-”是不是笔误？

我计算的结果是“+”，即：

x = x'cosθ + y'sinθ = Rcosα cosθ+Rsinα sinθ

***************************

y = x'sinθ - y'cosθ = Rcosα sinθ-Rsinα cosθ

 

 

还有，实际上二圆相交于二个点，

即应当有二个坐标。

另一个坐标是不是：

x = x'cosθ - y'sinθ = Rcosα cosθ-Rsinα sinθ

***************************

y = x'sinθ + y'cosθ = Rcosα sinθ+Rsinα cosθ

 

如果采用旋转椭球体进行计算，得到的精度只能提高一二秒钟时间，那么还真不如使用正球体来解日食的初亏时间了。

 

实际上，人们在观测日食时，并不太关心初亏是12：05：05发生，还是在12：05：08发生。

至少在观测前，都已经做好各种准备了。就等着日食的发生，多等几秒钟，没什么关系的。

这个不像发射火箭，误差一秒钟，可能轨道会偏几百甚至几千米。

看来，日出日没的初亏复圆线用两圆交点最好求了。

第一交点

x = x'cosθ - y'sinθ = Rcosα cosθ-Rsinα sinθ

y = x'sinθ + y'cosθ = Rcosα sinθ+Rsinα cosθ

第二交点

x = x'cosθ - y'sinθ = Rcosα cosθ+Rsinα sinθ

y = x'sinθ + y'cosθ = Rcosα sinθ-Rsinα cosθ

 

关于日出日没食甚线解法

 

贝塞尔基面是经过地心且垂直于月影轴的平面，设他与地球相交为圆（下文称为基圆），基圆半径为地球的平均半径。

在求解日出日没的食甚线的时候，相当于求各时刻在基圆上看到日食甚对应的点。

某时刻，月影轴在基面上的坐标为(x0,y0)，也就是半影的圆心坐标。

月影轴在基面上扫过的轨道几乎是一条直线（下文称为影足线），设此刻直线斜率为k。

在基圆上，地球自轴引起的站点速度是基本垂直于基面的，也就是说，观测点在基面上的分速度接近于0。因此，影轴相对于观测点的速度可以看做影轴速度，方向为与影足线相同。

当基面上的观测点与半影圆心的连线垂直于影足线，这个观测点看到日食甚。

 

基于以上原理，经过半影圆心做垂直于影足线的直线，得到与基圆的交点，就是日出(或日没)且食甚的点，当然，这个交点与到影足(圆心)的距离不能超过本影半径。所有各时刻的日出(日没)食甚点转到地理坐标，并连线构成日出日没食甚线。

 

此算法也有点麻烦，不知有没有更简单的方法。