《天文算法》第1章 注释与提示

第一章 注释与提示 解释如何计算和编程的内容并不属于本书范围，读者可以自行查阅读有关书籍。然而能写出好的程序并非一日之功，这是一门需要循序渐进的艺术，只有通过实践锻炼才能写出既短又好的程序。超出0360范围的角度的三角函数 超出0360范围的角度经常出现在天文计算中，在例24.a中我们可以看到1992年10月13日太阳的平经度角是2318.19281。快速运动的天体，比如月球、木星的伽利略卫星或者行星的自转中甚至还会出现更大的角度（可以看例41.a中第9步中的角w的计算实例）。 把角度转化到0360范围之内有时是很必要的，因为一些计算器或者程序设计语言对于大的角度的三角函数计算结果是不准确的。例如，你可以试试计算3600030的正弦。正确结果应该为0.5。 译者注：事实上，对于现代计算机而言，上一行所述的问题基本不存在。不过，很多情况下，我们仍需把一个角度转换到0到360度。角度的表示方法 算机不能直接计算出以度分秒方式表示的角度的三角函数。在使用三角函数之前，应该把角度转换成以度为单位的十进制小数形式。因此，在计算232649的余弦时要先把角度转化成为23.44694444,然后再用余弦函数计算。 遗憾的是，几乎所有的电脑都是用弧度而不是度来计算，所以还应把度转为弧度单位，通常是一件麻烦的事情。赤经 赤经通常用时、分、秒方式来表示。如果需要计算赤经的三角函数，需要把赤经转换成用度为单位来表示（然后再转换成以弧度为单位），请注意1h对应于15。 例1.a 计算 = 9h 14m 55s.8的tan 。我们首先把转化成以时为单位的十进制小数：9h 14m 55s.8 914/60+55.8/3600=9.248833333时然后再乘以15。138.73250然后在除以180/得到以弧度为单位的57.295779513。然后得出 tan =-0.877517修正角所在的象限 当已知一个角的正弦、余弦或正切值时，可以通过三角函数对应的反函数如正弦对应的反正弦（arcsin），余弦对应的反余弦(arccos)，正切对应的反正切(arctan)来得出角的大小。但是请注意在一些计算机上和一些程序设计语言中，尤其是大部分早期的微型计算机中都没有提供反正弦和反余弦函数。反三角函数并不是单值函数，例如，如果sin=0.5，那么可以是30、150、390等等。基于这个原因，计算机中的反三角函数的取值范围只有0360的一半：反正弦和反正切的取值范围是在9090，而反余弦的取值范围则是0180。 例如，计算cos147，结果是0.8387，用反余弦函数计算0.8387的结果正是147但是，cos213的结果也是0.8387，而我们用反余弦函数计算的结果则是147。因此当使用反正弦、反余弦和反正切的时候，必要时需要通过一个或多个值来弄清楚它代表的角度，消除结果的不确定性。另外，每个问题都要单独检查。 例如，公式（12.4）和（24.7）给出了天体赤纬的正弦。因为赤纬的取值范围在9090，所以反正弦函数可以在正确的象限算出赤纬，因此这里就不必进行检验。公式（16.1）给出了角度差的余弦也是同样情况，实际上角度差的取值范围在0180，这正与反余弦函数的取值范围一致。 但是请看从赤经（）赤纬（）转换到黄经（）黄纬（）的公式：cossin=sinsin+coscossincoscos=coscos 令第一个方程为A，第二个方程为B，用A式除以B式，我们可以得到tan=A/B，则对A/B使用反正切函数可以求出角，该角的取值在9090范围内，角度结果可能会相差180（由于正切函数的周期为180）。确定角所在的正确象限可以通过如下测试：如果B30，这个公式就不再有效，比如T= 3305.8，公式的结果将是e=1。一个认为“计算机从不犯错的”人可能会认为T= 3305.8时，天王星轨道是抛物线，进而认为天王星起源于太阳系之外这显然是伪科学。 实际上尽管行星轨道的偏心率e在超过了定义的时间上限后变化并不是有规律的，但是时间在很少的几个千年纪之内，偏心率是可以用像（1.1）那样的多项式精确表示的。 进一步的观察我们可以发现公式中有周期项（公式中的正弦和余弦项，在几个世纪内变化很小）和长期项的不同（如公式中含有T，T2，T3，的项，它随着时间的增加快速增大）。当T很小的时候，T2项会变得很小，但是当|T|值很大的时候这一项会变得非常重要。因此当|T|值比较大的时候考虑含有T2等项的周期项是没有意义的，在计算中也不用考虑。避免幂计算 假设我们计算这样一个多项式：Y=A+BX+CX2+DX3+EX4 其中A，B，C，D，E是常数，X是变量。现在可以在计算机中一项项直接相加来求出每一个给定X的多项式的值。然而可以采用一些聪明的方法来避免计算X的幂，比如： Y=A+X(B+X(C+X(D+EX) 在这个式子中幂计算都消失了，采用了乘积来替代幂的计算。这种多项式的表示方法被称为Horner方法，这种方法因为避免了幂计算，所以特别适合自动计算。 不用幂计算，而采用计算A*A的方法来计算A的平方也是一个聪明的办法，我们使用这样一段程序在HP-85计算机上计算前200个正整数的平方：FOR I=1 TO 200 K=I2NEXT I 完成计算需要费时10.75秒。但是，当我们把第二行换成K=I*I的话，完成整个计算只需费时0.96秒！缩短一个程序 把程序写尽可能短小通常不仅是代表着艺术，而且在计算机内存受限情况下也是必须的。即使对于简单的计算，也存在一些把程序缩短的技巧。假设我们要计算下面多项相加的和S：S=0.0003233sin(2.6782+15.54204T)+0.0000984sin(2.6351+79.62980T)+0.0000721sin(1.5905+77.55226T)+0.0000198sin(3.2588+21.32993T)+ 首先因为正弦的系数都是很小的数，可以通过采用以一个常数作为计数单位（在这个例子中是107）来避免输入那么多的数字，比如我们用3233来代替0.0003233。因此在计算了所有项之后，我们再把和除以107。 其次，在程序中声明全部数值项也是不明智的。相反，我们应该采用所谓的循环来完成计算。上面A*sin(BCT)的每一项的A，B，C值应作为程序的数据部分。假如有50项，程序可以这样写：doubles=0, data50*3=3233,2.6782,15.54204,; /每3个一组int i,p;for(i=0;i50;i+) p = i*3;s += datap*sin( datap+1+datap+2*T );s/=1e-7; 译者注：如果数据量较大，应考虑将数据与程序合理的分离；数据可以放在程序中的“数组”之中，如果数据量多达几万行或更多，多数情况下我们会考虑放在外部文件之中。安全性测试 在可能出现“不可能”出现的情况下，需要进行安全性测试。例如，在迭代到特定数量之后计算停止却没有达到要求的精度。 或者考虑月掩星的情况，在程序中根据当地的环境来计算被掩恒星消失和再次出现的时间。然而，可能在给定的地点根本看不到这颗恒星被掩。在这种情况下，初切时间和终切时间根本就不存在，试着计算这两个时间会碰到计算一个负数的平方根的情况。为了避免出现这个问题，程序应该首先计算这颗恒星到月面中心的最短距离（从给定的地点观看），而且当且仅当这个距离小于月面半径的情况下才计算初切和终切时间。调试 在程序写完之后必须要检查被称为Bug的错误。定位和修改程序中Bug的过程被称为调试。在无论使用什么程序设计语言编程，可能会有如下几种类型的Bug： a. 语法错误：不符合程序设计语言的规则，比如拼写错误、遗漏了括号，或者编程语言的保留字。例如，在BASIC中：A=SIM(B) 应该是 A=SIN(B)P=SQR(A+B) 应该是 P=SQR(A+B) b. 语义错误：比如遗漏程序行，例如在程序中不存在800行的情况下有GOTO 800语句。 c. 运行错误：在程序运行中出现的错误。例如：A=SQR(B)，在程序中计算B的平方根，但是B的值为负值。 d. 其他编程错误：下面的几种错误经常会发生： 错把字母“o”当作数字“0”输入，或者把字母“I”当作数字“1”输入。 同一变量名在程序中使用了两次（变量代表不同意思）。 输入数值常数的的错误（比如错把127.3当作127.03，或者错把15当作0.15），把*错当+。 使用单位错误。例如一个角没有用弧度而是用度来作为单位，或者赤经用时做单位而没有转换成度或者弧度。 角定位在错误象限。请见“修正角所在的象限”一节。 数值舍入错误。例如一个人在计算d的余弦值，在d很小的情况下使用cos()效果不好。实际上如果d非常小，它的余弦值就几乎等于1，而且除弦值随d变化非常缓慢。在这种情况下计算出来余弦值是不准确的。 比如，cos15=0.99999997，而cos0的准确值为1。如果想计算角度d的值非常小，可以通过其他方法来计算d的余弦值。比如，可以参考第16章。 一个不能保证覆盖各种情况的迭代过程。请见第5章（迭代）和第29章（Kepler方程的解法）。 使用了错误的计算方法。例如，交换两个数值X，Y，就需要使用一个中间变量A。不正确的方法：X=Y；Y=X；正确的方法：A=X；X=Y；Y=Z； 在QuickBASIC，GWBASIC或者其他BASIC版本中，有一个SWAP函数：SWAP（X,