农历及日月食函数
//------农历及日月食------////角度函数
function ang(x,t,c1,t0,t2,t3){
return tail(c1*x)*2*PI()+t0-t2*t*t-t3*t*t*t;
}
//返回农历日数及日月食信息的函数，如-324.57923415，负号表示闰月，百位3表示月偏食(2为月全食,1为日食0为无食),百位及十位表示日数,小数部分是朔望时刻(单位为天,若该天不朔或望则小数部分为零)
function lunDate(y,m,d,calType){
var t=(y-1899.5)/100;
var ms=floor((y-1900)*12.3685);
var rpi=180/PI();
var zone=8; //时区
var f0=ang(ms,t,0,0.75933,2.172e-4,1.55e-7)+0.53058868*ms-8.37e-4*t+zone/24+0.5;
var fc=0.1734-3.93e-4*t;
var j0=693595+29*ms;
var aa0=ang(ms,t,0.08084821133,359.2242/rpi,0.0000333/rpi,0.00000347/rpi);
var ab0=ang(ms,t,7.171366127999999e-2,306.0253/rpi,-0.0107306/rpi,-0.00001236/rpi);
var ac0=ang(ms,t,0.08519585128,21.2964/rpi,0.0016528/rpi,0.00000239/rpi);
var leap=0; //闰月数,0则不闰
var ecli=0; //日月食
var lunD=-1; //农历日数
var shuoD=0; //本阴历月的阴历朔日数
var shuoT=0; //本阴历月的朔时刻
var wangD=0; //本阴历月的阴历望日数
var wangT=0; //本阴历月的望时刻
for(var k=-1;k<=13;k+=0.5){ //k=整数为朔,k=半整数为望
var aa=aa0+0.507984293*k;
var ab=ab0+6.73377553*k;
var ac=ac0+6.818486628*k;
var f1=f0+1.53058868*k+fc*sin(aa)-0.4068*sin(ab)+0.0021*sin(2*aa)+0.0161*sin(2*ab)+0.0104*sin(2*ac)-0.0074*sin(aa-ab)-0.0051*sin(aa+ab);
var j=j0+28*k+f1; //朔或望的等效标准天数及时刻
//记录当前日期的j值
var lunD0=ESD(y,m,d,calType)-floor(j); //当前日距朔日的差值
if(k==floor(k)&&lunD0>=0&&lunD0<=29){
var k1=k; //记录当前时间对应的k值
shuoT=tail(j);
lunD=lunD0+1;
}
if(k==(k1+0.5)){
wangT=tail(j);
wangD=floor(j)-(ESD(y,m,d,calType)-lunD+1)+1;
}
//判断日月食
if((lunD==1&&k==k1)||(lunD==wangD&&k==(k1+0.5))){
if(abs(sin(ac))<=0.36){
var s=5.19595-0.0048*cos(aa)+0.002*cos(2*aa)-0.3283*cos(ab)-0.006*cos(aa+ab)+0.0041*cos(aa-ab);
var r=0.207*sin(aa)+0.0024*sin(2*aa)-0.039*sin(ab)+0.0115*sin(2*ab)-0.0073*sin(aa+ab)-0.0067*sin(aa-ab)+0.0117*sin(2*ac);
var p=abs(s*sin(ac)+r*cos(ac));
var q=0.0059+0.0046*cos(aa)-0.0182*cos(ab)+0.0004*cos(2*ab)-0.0005*cos(aa+ab);
if(p-q<=1.5572){
ecli=1; //日食
if(k!=floor(k)){
if(p+q>=1.0129)
ecli=3; //月偏食
else
ecli=2; //月全食
}
}
}
}
}
//k循环结束
// var v=lunD; //返回值
// if(v==1)
// v+=shuoT //朔日则返回朔的时刻
// if(v==wangD)
// v+=wangT; //望日则返回望的时刻
//return(v+ecli*100) 目录
第7章 儒略日
一般说明
儒略日的计算
由儒略日推算历日
给定日期间的相隔天数
星期几的问题
第8章 复活节日期
格里高利复活节
儒略历复活节
第9章 力学时和世界时
第7章 儒略日
本章中我们将给出一种方法，将儒略历或格里高利历1中的日期转换成相应的儒略日数（JD），以及反向转换。
一般说明
儒略日数(Julian Day number)，或简单地说，儒略日∗（Julian Day, JD），是指从公元前4712年开始连续计算日数得出的天数及不满一日的尾数。传统上儒略日的计数是从格林尼治平午，即世界时12点开始的；若以力学时（DT）（或历书时）为标准，我们则通常以历书儒略日(Julian Ephemiris Day, JDE**)来表示。例如，
1977年4月26.4日 UT = JD 2443 259.9
1977年4月26.4日 TD = JDE 2443 259.9
在下面将要描述的方法中，我们得考虑将其它历日换算为格里高利历日期。因而，1582年10月4日（儒略历）的下一日为1582年10月15日（格里高利历）。
格里高利历并非在所有国家迅速得到一致采用，这一点在历史研究中要引起注意。例如在大不列颠，晚至公元1752年才变更历法，而土耳其则要等到1927年。
儒略历由尤里乌斯•恺撒于公元前45年在罗马帝国创立，而其最终形式确立于公元8年前后，尽管如此，我们仍可以借助天文学家的演算无止境地向前推算儒略历。比如在这一系统中，我们可以说某次日食发生在公元前1203年8月28日，虽然在那个遥远的年代，罗马帝国根本还未被创建，而8月这个月份更有待设立。
对于公元1年之前的年份如何计数，天文学家同历史学家并不一致。在本书中，“公元前”的年份以天文方法计数。这样，+1年的前一年为0年，再之前才是-1年。所以历史学家所说的公元前585年实际上是-584年。
天文上以负数计数年份只是为算术目的起见。比如，在历史学计数中的，可被4整除的年份为儒略历闰年这个规则不再有效了。虽然像公元前1，5，9，13，…这些年份的确也在天文学上的闰年序列中，然而它们却被记为0，-4，-8，-12 …，它们都能被4整除。
我们以INT(x)来表示数x的整数部分，即小数点前的整数。例如：
INT(7/4) = 1 INT(5.02) = 5
INT(8/4) = 2 INT(5.9999) = 5
1 即通常所用阳历 --译者注
∗ 某些书中的称法是“儒略日期”（Julian Date）而非“儒略日”（Julian Day）。就本书而言，儒略日期指的是儒略历中的日期，正如格里高利日期所指为格里高利历中的日期一般。而儒略日与儒略历并不相关。
** 而并非像有时所写作的JED。E只是对JD的一种表征。
天文算法
这里负数可能会造成些问题。在某些计算机或程序语言中，INT(x)为小于等于x的最大整数。那样的话就有INT(-7.83) = -8，因为-7的确比-7.83大。
但在其它一些语言中，INT为所给数的整数部分，也即小数点前的部分。这样就有INT(-7.83) = -7，这被称为舍位(Truncation)。某些程序包含了这两项功能：表示前一种涵义的INT(x)以及TRUNC(x)。
因此，对负数进行取整时要特别注意（对正数，二者所得结果一致）。在本书给出的算式中，INT的使用总是针对正数的。
儒略日的计算
下面的方法对正数年和负数年都是有效的，负的儒略日数除外。
设Y为给定年份，M为月份，D为该月日期（可以带小数）。
􀁺 若M > 2，Y和M不变.
若 M =1或2，以Y–1代Y，以M+12代M.
换句话说，如果日期在1月或2月，则被看作是在前一年的13月或14月。
􀁺 对格里高利历，有
对儒略历，取B = 0.
􀁺 要求的儒略日即为
实际上用30.6（取代30.6001）即可得出正确结果，但这里用30.6001以确保总能取得适当的整数。[事实上可用30.601甚至30.61来取代30.6001.] 打个比方，5乘30.6刚好是153，然而大多数计算机不能精确表示出30.6 — 参见第二章关于BCD的说明 — 这导致得出一个152.999 9998的结果，它的整数部分为152，如此算出的JD就不正确了。
例7.a — 计算1957年10月4.81日，即苏联发射第一颗人造卫星的日期相应的儒略日。
这里我们取 Y =1957，M = 10，D = 4.81.
因 M > 2，则Y与M不变.
日期为格里高利历，于是算得
例 7.b — 计算333年1月27日12时相应的儒略日。
因 M = 1，有 Y = 333 -1 =132 且 M = 1 +12 =13.
因日期为儒略历日期，故 B = 0.
下表列出了一些历日所对应的儒略日，可作测试程序之用。
2000年1月1.5日 2451 545.0
1987年1月27.0日 2446 822.5
1987年6月19.5日 2446 966.0
1988年1月17.0日 2447 187.5
1988年6月19.5日 2447 332.0
1900年1月1.0日 2415 020.5
1600年1月1.0日 2305 447.5
1600年12月31.0日 2305 812.5
837年4月10.3日 2026 871.8
-1000年7月12.5日 1356 001.0
-1000年2月29.0日 1355 866.5
-1001年8月17.9日 1355 671.4
-4712年1月1.5日 0.0
如果你只对1900年3月1日至2100年二月28日之间的日期感兴趣，那么在公式（7.1）中可取B = -13。
在一些应用中需要知道某一给定年1月0.0日的儒略日JD0，它与上一年的12月31.0日是同一日。对格里高利历中的年份可以如下计算。
特别地，对1901年至2099年间的年份，上式简化为
闰年是如何设定的？
在儒略历中，能被4整除的年份为闰年，这一年有366天，其它年份为平年（365天）。
如900年和1236年为闰年，而750年和1429年为平年。
格里高利历法也采用这一规则，但下列年份除外：不能被400整除的百年为平年，如1700年，1800年，1900年和2100年。其余能被400整除的百年则为闰年，如1600年，2000年和2400年。
现代工作中有时会用到修正儒略日（Modified Julian Day, MJD），比如在提到人造卫星轨道根数时。对比儒略日，修正儒略日从格林威治时间平午夜开始计算，相当于
MJD = JD – 2400000.5
因而MJD = 0.01指向的是1858年11月17日世界时0时。
由儒略日推算历日
下面的方法对正数年和负数年都是有效的，负儒略日数除外。
将JD加上0.5，令 Z 为其整数部分，F 为尾数（小数）部分。
若 Z < 2299 161，取A = Z.
若 Z 大于等于2299 161，计算
然后计算
该月日期（带小数部分）则为
月份m为
E – 1 若 E < 14
E – 13 若 E = 14或15
年份为
C – 4716 若 m > 2
C – 4715 若 m =1 或 2
对比早先公式(7.1)的情况，这个公式里求E时用的数30.6001不能代之以30.6，哪怕计算机没有先前所说的问题。否则，你得到的结果会是2月0日而不是1月31日，或者4月0日而不是3月31日。
例 7.c — 计算JD 2436 116.31所对应的历日。
2436 116.31 + 0.5 = 2436 116.81
Z = 2436 116 且 F = 0.81
因 Z > 2299 161，故有
于是我们求得
B = 2437 653 C = 6637 D = 2437 313 E = 11
该月日期 = 4.81
月份 m = E – 1 = 10 (因 E < 14 )
年份 = C – 4716 = 1957 (因 m > 2 )
因此，求得的日期为1957年11月4.81日。
练习：计算JD = 1842 713.0 和 JD = 1507 900.13 所对应的历日。
答案：333年1月27.5日 和 -584年5月28.63日。
给定日期间的相隔天数
两个历日间的相隔天数可以通过计算它们相应的儒略日之差求得。
例 7.d — 周期彗星哈雷于1910年4月20日及1986年2月9日两次过近日点。这两次经过的时间间隔是多少？
1910年4月20.0日对应 JD 2418 781.5
1986年2月9.0日对应 JD 2446 470.5
其差值为27 689天。
练习： 1991年7月11日之后恰相隔10 000天的日期。
答案： 2018年11月26日。
星期几的问题
给定日期是星期几可由以下方法获得 — 计算该日0时的儒略日，加上1.5，再除以7 ，所得余数将指示出星期几：若余数为0，则为星期日，1为星期一，2为星期二，3为星期三，4为星期四，5为星期五，6为星期六。
儒略历到格里高利历的换算并不影响星期。因而，在1582年，11月4日星期四接下来的一天便是11月15日星期五。
例 7.e — 找出1954年6月30日是星期几。
1954年6月30.0日对应 JD 2434 923.5
2434 923.5 + 1.5 = 2434 925
2434 925除以7所得的余数为3.
所以这一天是星期三。
年内的序数日
年内的序数日 N 可由以下公式得出[1]。
此处 M 为月份，D为该月日期，另有
K = 1 若为闰年
K = 2 若为平年
N取整数，自1月1日开始取值1，直至12月31日取值365（或闰年取值366）。
例 7.f — 1978年11月14日。
平年，M = 11, D = 14, K = 2.
得出 N = 318.
例 7.g — 1988年4月22日。
闰年， M = 4， D = 22，K = 1.
得出 N = 113.
现在让我们反过来考虑： 年内序数日N已知，要求相应的日期，M为月份，D为该月日期。下面的算法是由比利时民间天文协会的A. Pouplier发现的[2]。
如上所述，令
K = 1 若为闰年
K = 2 若为平年
若 N < 32 ，则令 M = 1
参考书目
[1] 《Almanac for Computers for the Year 1978》第B2页，特区华盛顿，美国海军天文台，航海天文年历办公室。
[2] A. Pouplier致Jean Meeus的信件， 1987年4月10日。
第8章 复活节日期
本章我们将给出计算指定年基督教复活节星期日（Easter Sunday）日期的方法 — 注意并非犹太教的逾越节1。
格里高利复活节
以下方法由Spencer Jones在他的《普通天文学》（General Astronomy）(1922年版73-74页)中给出。这一方法后来重新印刷在《英国天文协会期刊》（Journal of the British Astronomical Association）第88卷91页（1977年12月）中，并说明这一方法发明于1876年，曾出现在《Butcher’s Ecclesiastical Calendar》上。
和高斯给出的公式不同，这种方法适用于格里高利历中始自1583年的所有年份，并无例外。找出复活节日期的程序如下：
则有 n = 月份 （3 = 3月，4 = 4月），
p + 1 = 本月内复活节星期日所在的日期。
如果所用的计算机语言没有模（modulo）或求余(remainder) 功能，则计算余数时必须特别注意。假设要求34除以30的余数，在HP-67及HP-41C便携式计算器上我们得到
34/30 = 1.133 333 333
尾数部分为0. 133 333 333，将它乘上30时，结果变成了3.999 999 990，这个结果不等于原本的正确答案4。最终这可能导致得到一个错误的复活节日期。
1逾越节 (出埃及记 12; 利未记 23:5) 开始於尼散月十四日（宗教历的正月，即阳历3月至4月期间），连接七天，每天要献许多祭，头一日和第七日休息，给百姓守严肃会。
在正月十四日逾越节开始那天，耶稣被钉十架 (约翰福音 19:14)，三天后的星期日身体复活，始为复活节。 — 译者注
可用以下年份测试你的程序：
1991年 → 3月31日 1954年 → 4月18日
1992年 → 4月19日 2000年 → 4月23日
1993年 → 4月11日 1818年 → 3月22日
复活节日期总在3月22日（如1818年和2285年）到4月25日（如1886年，1943年和2038年）之间。
复活节星期日的推算规则众所周知：它为春分后（含春分）的第一次满月后的第一个星期日。事实上，推算规则在很久以前就由基督教神职人员确定了。为应用这一规则所使用的满月日期是通过一种教会推算得出的结果而不是真实的天文学上的满月；同样，春分日也总是被定在3月21日，但实际上则可能会早一到两天。
例如，1967年的春分在3月21日，满月在3月26日（UT），3月26日之后的第一个星期日是4月2日。然而，复活节星期日却在3月26日。
在1900至2100年间，使用纯粹天文规则推算出了一些不同于教会规则的复活节星期日，这些年份是：1900，1903，1923，1924，1927，1943，1954，1962，1967，1974，1981，2038，2049，2069，2076，2089，2095，2096。
格里高利历的复活节日期要5 700 000年才循环一次。有人发现格复活节最常见的日期是4月19日。
儒略历复活节
儒略历中的复活节日期可按以下方法推算。
则有 f = 月份 （3 = 3月，4 = 4月），
g + 1 = 本月内复活节星期日所在的日期。
儒略历的复活节日期有为期532年的循环。比如，179年，711年和1243年的复活节都在4月12日。
第9章 力学时和世界时
世界时（UT），即格林威治地方时，基于地球的自转。它对日常生活及与地方时角（local hour angle）相关的天文计算都是必需的。
然而，地球自转一直在变缓，再加之速度变化难以预测，这使得世界时成了一种非均匀的时间系统。
但天文学家们需要一个严格均匀的时间尺度来进行精确计算（如天体力学，轨道，历表）。1960年至1983年间，诸如《Astronomical Ephemeris》等一些大型的天文年历中使用了一种定义于运动定律的称为历书时（ET）的均匀时间尺度：它基于行星的运动。1984年，历书时又被秒长定义于原子钟的力学时取代。实质上力学时是历书时的一个拓展。
这里还有一个质心力学时（TDB）和地球力学时（TDT）的区分。这两个系统最多相差0.0017秒，此种差异与地球以椭圆轨道绕日运动有关(相对论效应)。因这一差异小到可以被大多数实际应用忽略，故此处我们对质心力学时和地球力学时不加区分，统称力学时(TD)。
力学时和世界时之间的精确差值 ΔT = TD - UT 只能由天文台推算。表9.A列出了一些年份年初时的 ΔT值，除最后两个值外均取自1988年的《Astronomical Ephemeris》[1]。
对于不远的将来1，我们可以继续推算表9.A中的值。例如，可使用预备值
ΔT = + 60 秒 1993年
ΔT = + 67 秒 2000年
ΔT = + 80 秒 2010年
对表9.A以外的世纪， Morrison及Stephenson则给出以下关系式来推导ΔT的近似值[2]：
此处的year可以带小数。这个公式还可写成以下形式
此处T为从公元2000年算起的以世纪为单位的时间差。或者使用儒略日表示
根据这些表达式，可以推算到公元前4000年时，世界时将出现多达2小时的差异。将来对公式的改进将有益于TD和UT间的换算，但它不会影响算法、程序、历法，以及目前以TD为均匀时间尺度的表格。
1984年Stephenson和Morrison[3]又发表了另两个抛物线表达式来推算以往年代的ΔT。公元前390至公元后1600年适用以下两种抛物线拟合：
1 指对本书出版时的1991年而言。 — 译者注
表 9.A
- 390年至 +948年：
+948年至 +1600年：
T为从公元1800年算起相差的世纪数，得到的ΔT以秒为单位。
两年后，Stephenson和Houlden[4]给出了另两个推算以往年代ΔT的表达式：
(i) 公元948年以前的时段： ΔT=1830-405E+6.5E2
(ii) 公元948年至1600年： ΔT=22.5t2
E为从公元948年算起的世纪差数，t为从公元1850年算起的世纪差数。
公式(i)和(ii)等价于以下表达式，该处T为从公元2000.0年算起的世纪差数( T < 0 )。
公元948年以前的时段：
公元948年至1600年：
公元1871年至1901年之间的ΔT值为负。值得注意的是，无论对遥远的过去还是遥远的将赤，ΔT值都是正的
。
1871年至1901年间的时段是例外，那时对世界时和力学时显示的同一个字面时刻，世界时总比力学时要晚。比如，1990年21月27日0时（UT）比1990年1月27日0时（TD）要晚上57秒。我们使用的是 UT = TD - ΔT。
例9.a — 某次新月开始于1977年2月18日3h37m40s（力学时）（见例题47.a）。
该时刻的ΔT值等于+48秒。于是该月相相应的世界时为
例9.b — 假设需要计算+333年2月6日日6时（UT）时的水星位置。
我们有
代入公式（9.1）得出ΔT = +7074秒，即118分钟。因此，TD = 6h + 118分钟 = 7h58m，故应使用力学时333年2月6日7h58m进行上述计算。
2 1990年并非1871年至1901年间的年份，此处疑有误。— 译者注
Schmadel与Zech[5]构造了以下渐近法来模拟ΔT值，它对1800-1988年间所有时段有效。该式得出的结果与表9.A最多有1.9秒的误差。
上式中ΔT以天为单位，θ为从1900.0年算起经过的儒略世纪数。
Schmadel与Zech也提供了计算更短时间段的表达式。对于1800-1899年间，以下表达式给出的ΔT值（以天为单位）至多误差1.0秒：
对1900-1987年间，以下表达式给出的ΔT值（以天为单位）至多误差1.0秒：
这里θ的指代与上述第一个表达式相同。
请注意以上3式均为经验表达式，它们在各自的定义域外是无效的。
参考书目
[1]《Astronomical Almanac for 1988》（特区华盛顿） K8-K9页。
[2] L. V. Morrison与F. R. Stephenson著《Sun and Planetary System》（Reidel，Dordrecht，1982出版），96卷第73页。— 由P. Bretagnon及J. –L. Simon引用于《Planetary Programs and Tables from -4000 to +2800》(Willmann-Bell，Richmond，1986出版)，第5页。
[3] F. R. Stephenson与L. V. Morrison著《Long-term Changes in the Rotation of the Earth》（Phil. Trans. Royal Soc.，1984年出版），A，313卷第47-70页。
[4] F. R. Stephenson及M. A. Houlden著《Atlas of Historical Eclipse Maps》，剑桥大学出版社，英格兰，1986年出版，第x页。
[5] L. D. Schmadel及G. Zech著《Empirical Transformations from U.T. to E.T. for the period 1800-1988》，Astronomische Nachrichten，1988年出版，309卷第219-331页。
__
34# 大 中 小 发表于 2008-3-29 09:22 只看该作者
坛主今夏江南四城市之行纪实
牧夫官方淘宝店拥有各种优质望远镜和天文器材！
月球位置(第45章)
［ 许剑伟,2008-02-27日,译于莆田十中］
译者注1:因本人不很了解球面天文学的相关术语,所以下文用到一个自创名词"地心Date平黄道分点",意思包含 1)是黄道坐标系;(2)是瞬时黄道;(3)是平黄道,不含黄经章动;(4)黄经从黄道赤道升交点起算;(5)黄纬不会受章动引影;(6)右手坐标系,即逆旋为正;(7)是坐标原点建立在地心
译者注2:T^2表示T的2次方,同理T^3表示T的3次方
译者注3:表格计算并没有直接翻译,而是自已写了一段更详细的半数学化的文字进行表述,原文讲述得过于简单。
为了准确计算出某时刻月球的准确位置,须计算月球黄经黄纬及距离的数百个周期项。这已超出本书的范围,这里仅考虑主要的周期项,得到的黄经精度是10",纬度精度是4"。
利用本章描述的算法,可得到地心Date平黄道分点(译者注:平黄道与平赤道的升交点,近似春风点)的月心位置坐标:黄纬(λ)、黄纬(β)及地心到月心距离(Δ千米)。
赤道地平视差π由下式获得:
sinπ=6378.14/Δ
一、计算方法:
本章的周期项是基于ELP-2000/82月球理论。但L',D,M,M',F平参数使用Chapront的改进表达式。
T使用21.1式计算,T表达为J2000起算的世纪数,并取足够的小数位数(至少9位,每0.000 000 001世纪月球移动1.7角秒)。
使用以下表达式计算角度L',D,M,M',F,角度单位是度。为避免出现大角度,最后结果还应转为0—360度。
月球平黄经:
L'=218.3164591+481267.88134236T-0.0013268T^2+T^3/538841-T^4/65194000
月球距角(从地心看月日在天球上的角距离):
D =297.8502042+445267.1115168T-0.0016300T^2+T^3/545868-T^4/113065000
太阳平近点角:
M=357.5291092+35999.0502909T-0.0001536T^2+T^3/24490000
月亮平近点角:
M'=134.9634114+477198.8676313T+0.0089970T^2+T^3/69699-T^4/14712000
月球纬度参数(升交点起算的平角距):
F =93.2720993+483202.0175273T-0.0034029T^2-T^3/3526000+T^4/863310000
三个必要的参数:
A1=119.75+131.849T
A2= 53.09+479264.290T
A3=313.45+481266.484T
取和计算45.A表中各项(ΣI及Σr),取和计算45.B表中各项(Σb)。ΣI与Σb是正弦项取和,Σr是余弦项取和。正余弦项表达为A*sin(θ)或A*cos(θ)，式中的θ是表中D、M、M'、F的线性组合,组合系数在表45.A及45.B相应的列中,A是振幅。
以表45.A第8行为例:
I8 = A*sin(θ) = +57066 * sin( 2D-M-M'+0 )
r8 = A*cos(θ) = -152138 * cos( 2D-M-M'+0 )
同理可计算第1、2、3、4....各行,得到I1、I2、I3...及r1、r2、r3...
最后ΣI=I1+I2+I3+...；Σr=r1+r2+r3+...
然而,表中的这些项包含了了M(太阳平近点角),它与地球公转轨道的离心率有关,就目前而言离心率随时间不断减小。由于这个原因,振幅A实际上是个变量(并不是表中的常数),角度中含M或-M时,还须乘上E,含2M或-2M时须乘以E的平方进行修正。E的表达式如下:
E=1 - 0.002516T - 0.0000074T^2
此外,还要处理主要的行星摄动问题。A1与金星摄动相关,A2与木星摄动相关,L'与地球扁率摄动相关。
ΣI += +3958 * sin( A1 )
+1962 * sin( L' - F )
+318 * sin( A2 )
Σb += -2235 * sin( L' )
+382 * sin( A3)
+175 * sin( A1 - F )
+175 * sin( A1 + F )
+127 * sin( L' - M')
-115 * sin( L' + M')
最后得到月球的坐标如下:
λ = L'+ ΣI/1000000 (黄经单位:度)
β = Σb/1000000 (黄纬单位:度)
Δ = 385000.56 + Σr/1000 (距离单位:千米)
因45.A及45.B表中的振幅系数的单位是10^-6度及10^-3千米,所以上式计算时除以1000000和1000。
二、两个计算用的表:
［表45.A］
月球黄经周期项(ΣI)及距离(Σr).
黄经单位:0.000001度,距离单位:0.001千米.
--------------------------------------------------
角度的组合系数 ΣI的各项振幅A Σr的各项振幅A
D M M' F (正弦振幅) (余弦振幅)
--------------------------------------------------
0 0 1 0 6288744 -20905355
2 0 -1 0 1274027 -3699111
2 0 0 0 658314 -2955968
0 0 2 0 213618 -569925
0 1 0 0 -185116 48888
0 0 0 2 -114332 -3149
2 0 -2 0 58793 246158
2 -1 -1 0 57066 -152138
2 0 1 0 53322 -170733
2 -1 0 0 45758 -204586
0 1 -1 0 -40923 -129620
1 0 0 0 -34720 108743
0 1 1 0 -30383 104755
2 0 0 -2 15327 10321
0 0 1 2 -12528 0
0 0 1 -2 10980 79661
4 0 -1 0 10675 -34782
0 0 3 0 10034 -23210
4 0 -2 0 8548 -21636
2 1 -1 0 -7888 24208
2 1 0 0 -6766 30824
1 0 -1 0 -5163 -8379
1 1 0 0 4987 -16675
2 -1 1 0 4036 -12831
2 0 2 0 3994 -10445
4 0 0 0 3861 -11650
2 0 -3 0 3665 14403
0 1 -2 0 -2689 -7003
2 0 -1 2 -2602 0
2 -1 -2 0 2390 10056
1 0 1 0 -2348 6322
2 -2 0 0 2236 -9884
//继续
0 1 2 0 -2120 5751
0 2 0 0 -2069 0
2 -2 -1 0 2048 -4950
2 0 1 -2 -1773 4130
2 0 0 2 -1595 0
4 -1 -1 0 1215 -3958
0 0 2 2 -1110 0
3 0 -1 0 -892 3258
2 1 1 0 -810 2616
4 -1 -2 0 759 -1897
0 2 -1 0 -713 -2117
2 2 -1 0 -700 2354
2 1 -2 0 691 0
2 -1 0 -2 596 0
4 0 1 0 549 -1423
0 0 4 0 537 -1117
4 -1 0 0 520 -1571
1 0 -2 0 -487 -1739
2 1 0 -2 -399 0
0 0 2 -2 -381 -4421
1 1 1 0 351 0
3 0 -2 0 -340 0
4 0 -3 0 330 0
2 -1 2 0 327 0
0 2 1 0 -323 1165
1 1 -1 0 299 0
2 0 3 0 294 0
2 0 -1 -2 0 8752
--------------------------------------------------
［表45.B］
月球黄纬周期项(ΣI).单位:0.000001度.
-------------------------------------
角度的组合系数 ΣI的各项振幅A
D M M' F (正弦振幅)
-------------------------------------
0 0 0 1 5128122
0 0 1 1 280602
0 0 1 -1 277693
2 0 0 -1 173237
2 0 -1 1 55413
2 0 -1 -1 46271
2 0 0 1 32573
0 0 2 1 17198
2 0 1 -1 9266
0 0 2 -1 8822
2 -1 0 -1 8216
2 0 -2 -1 4324
2 0 1 1 4200
2 1 0 -1 -3359
2 -1 -1 1 2463
2 -1 0 1 2211
2 -1 -1 -1 2065
0 1 -1 -1 -1870
4 0 -1 -1 1828
0 1 0 1 -1794
0 0 0 3 -1749
0 1 -1 1 -1565
1 0 0 1 -1491
0 1 1 1 -1475
0 1 1 -1 -1410
0 1 0 -1 -1344
1 0 0 -1 -1335
0 0 3 1 1107
4 0 0 -1 1021
4 0 -1 1 833
0 0 1 -3 777
4 0 -2 1 671
2 0 0 -3 607
2 0 2 -1 596
2 -1 1 -1 491
2 0 -2 1 -451
0 0 3 -1 439
2 0 2 1 422
2 0 -3 -1 421
2 1 -1 1 -366
2 1 0 1 -351
4 0 0 1 331
2 -1 1 1 315
2 -2 0 -1 302
0 0 1 3 -283
2 1 1 -1 -229
1 1 0 -1 223
1 1 0 1 223
0 1 -2 -1 -220
2 1 -1 -1 -220
1 0 1 1 -185
2 -1 -2 -1 181
0 1 2 1 -177
4 0 -2 -1 176
4 -1 -1 -1 166
1 0 1 -1 -164
4 0 1 -1 132
1 0 -1 -1 -119
4 -1 0 -1 115
2 -2 0 1 107
-------------------------------------
三、计算举例:
例45.a, 计算月球的地心黄经、黄纬、距离及赤道视差,时间1992年4月0时(力学时), 结果如下:
JDE = 2448724.5(儒略日) A1 = 109°.57
T = -0.077221081451 A2 = 123°.78
L'= 134°.290186 A3 = 229°.53
D = 113°.842309 E = 1.000194
M = 97°.643514 ΣI =-1127527 (含A1,A2等项)
M'= 5°.150839 Σb =-3229127 (含A1,A2等项)
F = 219°.889726 Σr =-16590875
从以上算出:
λ = 134°.290186 - 1°.127527 = 133°.162659
β = -3°.229127 = -3°13'45"
Δ = 385000.56 - 16590.875 = 368409.7 km
π = arcsine(6378.14/368409.7)=0°.991990=0°59'31".2
要获得地心视黄经,还应加上黄经章动(Δψ),Δψ=+16".595=+0°.004610。
λ视=133°.162659 + 0°.004610
=133°.167269
=133°10'02"
瞬时黄赤交角=平黄赤交角(εo)+交角章动(Δε)
ε=εo + Δε=23°26'26".29 = 23°.440636
(注:章动计算详见21章)
这样就可得到月球的地心视赤经和视赤纬:
α = 134°.688473 = 8h 58m 45s.2
δ = +13°.768366 =+13°46' 06"
利用完整的ELP-2000/82月球理论获得的准确值是(注:不妨同以上计算结果比较):
λ = 133°10'00" α = 8h 58m 45s.1
β = -3°13'45" δ = +13°46' 06"
Δ = 368405.6 km π = 0°59' 31".2
四:月球的升交点和近地点
根据Chapront[2],月球升交点(平)黄经Ω 及(平)近点角π,可由以下二式计算(单位是度)
Ω = 125.0445550 - 1934.1361849T + 0.0020762T^2 + T^3/467410 - T4/60616000
π = 83.3532430 + 4069.0137111T - 0.0103238T^2 - T^3/80053 + T4/18999000
式中T的单位与上文的相同(即:J2000起算的世纪数).这些经度是指黄经(Date平黄道分点起算的经度)。
从Ω的公式中,我们可以找到升(或降)交点等于春风点的瞬时,即Ω=0°或180°。在1910至2110期间,这种情况发生在如下日期:
Ω=0° Ω=180°
----------------------------
1913年05月27 1922年09月16
1932年01月06 1941年04月27
1950年08月17 1959年12月07
1969年03月29 1978年07月19
1987年11月08 1997年02月27
2006年06月19 2015年10月10
2025年01月29 2034年05月21
2043年09月10 2052年12月30
2062年04月22 2071年08月12
2080年12月01 2090年03月23
2099年07月13 2108年11月03
36# 大 中 小 发表于 2008-3-29 09:57 只看该作者
坛主今夏江南四城市之行纪实
牧夫官方淘宝店拥有各种优质望远镜和天文器材！
大气折射(第15章)
[许剑伟 于莆田十中 2008年3月]
[译者注]以下所述的纬度均指"地平纬度"或"高度角",我不会译,就这么将就吧。
大气折射是光线通过地球大气层时光线发生弯曲。光线经过密度不断增加的大气时发生连续弯曲。这造成观测到的星体位置比真实位置高。在天顶,大气折射是零,越接近地平,折射越大。地平纬度是45度时,折射为约为1'，在地平上,大约是35'。因此,太阳和月亮升起的时候,它们实际在地平线之下。由于折射率的变化,在低纬时我们看到的是椭圆形的太阳。当确定位置是,必须做大气折射修正，有以下两种情形：
•观测到星体地平纬度是ho,我们应找到一个适当的R值修正ho,真纬度(真高度角)是h=ho-R
•从天体坐标中得到没有空气情况下地平真纬度h,找一个适当的R修正h,,得到视纬度ho=h+R
我们遇到的大部分公式都是针对第一种情况的(已知观测值求真值)。但是这里,我们将考虑两种情况。
通常,我们可以使用'平均'的方法。然而,接近地平是的反常折射则不行,变形的夕阳造诉我们,在低纬度时,无法得到很高的精度。当天体的纬度大于15度,以下两个公式可供你选择,你可根据实际情况选择其一:
R=58".294tan(90-ho)-0".0668tan3(90°-ho)
R=58".276tan(90-h) -0".0824tan3(90°-h)
第1个公式是Smart提供的,第2个公式是由第一个公式导出的。当纬度低于15度,这两个表达式将变得不准确,甚至毫无意义。 从以上公式看出,在高纬度区,折射与90-h的正切值成正比。
New South Wales大学的G.GBennett结出了一个出人意料的简单的折射公式,在0到90度范围内有很好的精度。
如果R表达为以分为单位,Bennett's公式是:
R=1/tan(ho+7.31/(ho+4.4))...式1
式中ho是视纬度,单位是度。在0到90度范围内,精度是0.07'=4.2"。应当注意的是:当ho=0时,R=-0".08(即0.0013515分),而不是0。可用第二项公式修正,先算出R,接下来利用下式修正R,
dR=-0.06sin(14.7R+13),
结果的单位是分。括号中的表达式单位是度。修正后,在ho=0到90°范围内,最大误差0.015'=0.9"。注意,在ho=90°时,计算的结果是R=-0.89",不作第二项修正反而更好。
逆问题,已知真纬度,求折射的影响。有以下公式:
R=1.02/tan(h+10.3/(h+5.11))...式2
该式与Bennett的公式勿合到4"。同样,h=90°时,该式算得R不等于零。差值是0.0019279。
以上公式假设观测式是在海平面,大气压是1010毫巴,温度10度摄氏度。
当气压增加温度下降,折射增加。设地表气压为P毫巴,气温是T摄氏度，那么以上各式R的值应乘以下式:
(P/1010)*(283/(273+T))
然而,这只是大约修正。因为折射率还与光的波长有关。这些表达式适用于黄光,它对人眼的灵敏度最高。
例15.a:表面光滑的太阳圆盘下边沿视纬度是30'。设太阳的真直径是32',气温及大气压为常规条件。求真位置。
(1)ho1=30'=0.5°,由式1算得R=28'.754,则
(2)下边沿真地平纬度h1=30'-28'.754=1'.246,
(3)上边沿真地平纬度h2=1'.246+32'=33'.245
(4)利用式2算得上边沿的视地平纬度（高度）ho2=57'.864
日圆盘的垂直方向直径与与水平方向的直径比 ho2-ho1)/32=(57.864-30)/32=0.871 [分享]彗星轨道根数算法文件（β版）
Attribute VB_Name = "TrackEquation"
'彗星轨道根数算法文件
Option Explicit
Private Const pi = 3.14159265358979 '给定π的值
Public Function ComputePosition(year As Long, month As Long, day As Double, q As Double, e As Double, w As Double, n As Double, i As Double, cType As Long) As Single
'定义数据库中轨道参数相应的变量
'定义名称，年，月，日，q，e，ω，Ω，i， H，G
'Dim name As String
'定义变量
Dim epsilon As Double '定义黄赤交角变量(角度)
Dim longradii As Double '定义长半径变量
epsilon = 23.439375 '给定2000.0的黄赤交角ε常量值
epsilon = epsilon * pi / 180
i = i * pi / 180 '将角度转换为弧度 (弧度＝角度×pi/180)
'(角度＝弧度×180/pi)
n = n * pi / 180
w = w * pi / 180
'在sin之类的函数计算中，sin(变量)一式中的变量应该为弧度值
'所有角度值应先转换为弧度值后进行计算
'定义6个高斯常数A,B,C,a,b,c
Dim aa As Double, bb As Double, cc As Double '定义A,B,C变量
Dim a As Double, b As Double, c As Double '定义a,b,c变量
longradii = q / (1 - e) '计算长半径数值
'---A--------------------------------------------------
aa = Atn(-Cos(n) / Cos(i) / Sin(n)) '计算A的值(弧度)
aa = aa * 180 / pi '弧度转换为角度
'---B--------------------------------------------------
bb = Atn(Sin(n) * Cos(epsilon) / (Cos(i) * Cos(epsilon) * Cos(n) - Sin(i) * Sin(epsilon))) '计算B的值(弧度)
bb = bb * 180 / pi '弧度转换为角度
'---C--------------------------------------------------
cc = Atn(Sin(n) * Sin(epsilon) / (Cos(i) * Cos(n) * Sin(epsilon) + Sin(i) * Cos(epsilon)))
cc = cc * 180 / pi '弧度转换为角度
'_________________ ____________________
' 验算 开始
Dim sina As Double, sinb As Double, sinc As Double '定义sina,sinb,sinc
sina = Cos(n) / Sin(aa * pi / 180)
sinb = Sin(n) * Cos(epsilon) / Sin(bb * pi / 180)
sinc = Sin(n) * Sin(epsilon) / Sin(cc * pi / 180)
'如果sa*sa+sb*sb+sc*sc-2的绝对值大于0.00001,说明数据计算有问题
If Abs(sina * sina + sinb * sinb + sinc * sinc - 2) > 0.0000001 Then ComputePosition = Abs(sina * sina + sinb * sinb + sinc * sinc - 2)
'_________________验算 完成____________________'验算公式二和三可以不必使用，通常用第一式即可
'定义6个辅助量Px,Py,Pz,Qx,Qy,Qz
Dim px As Double, py As Double, pz As Double
Dim qx As Double, qy As Double, qz As Double
px = sina * Sin(aa * pi / 180 + w)
py = sinb * Sin(bb * pi / 180 + w)
pz = sinc * Sin(cc * pi / 180 + w)
qx = sina * Cos(aa * pi / 180 + w)
qy = sinb * Cos(bb * pi / 180 + w)
qz = sinc * Cos(cc * pi / 180 + w)
'----验算-----
If Abs(px * px + py * py + pz * pz - 1) > 0.0000001 Then ComputePosition = Abs(qx * qx + qy * qy + qz * qz - 1)
'___________________________________________________
Dim cosfi As Double '定义cosφ变量
cosfi = Sqr(1 - e * e) '计算出cosφ的数值
'定义6个大写坐标Ax,Ay,Az,Bx,By,Bz
Dim ax As Double, ay As Double, az As Double
Dim bx As Double, by As Double, bz As Double
longradii = q / (1 - e) '计算出轨道半长径
'下面分别计算Ax,Ay,Az,Bx,By,Bz的数值
ax = longradii * px
ay = longradii * py
az = longradii * pz
bx = longradii * cosfi * qx
by = longradii * cosfi * qy
bz = longradii * cosfi * qz
Dim GS As Double '定义天体平均运动速度 n 变量 其值为sqr(GS/longradii的3次方)
'计算出的GS的值是 x.xxxxxx度/日
GS = Sqr(0.000295912208 / (longradii * longradii * longradii)) * 180 / pi
'_____________________________________________________________________________
Dim M As Double '定义平近点角M变量 该变量为开普勒方程的一个常量
'需要判断已知日期同计算日期之间的时间差
M = days360("1998-11-21", "1999-1-12") * GS
M = 7.71
'此处应计算两个日期之间的时间差，精确到小数点后5位
'_____________________________________________________________________________
'通过渐进法(牛顿迭代法)计算开普勒方程并找出其解来
Dim ee As Double '定义偏近点角E变量
Dim e0 As Double '定义中间变量E0
'将M转换为弧度值，如果M仍旧是角度值时，需要根据上面的计算而知
M = M * pi / 180
e0 = M '将M赋值给E0
count: ee = M + e * Sin(e0)
If ee <> e0 Then e0 = ee: GoTo count
If ee = e0 Then ee = e0
'变量E通过计算并已经给出数据
'ee(数据是弧度数)
'定义cosE，sinE，cosE-e
Dim cose As Double, sine As Double, cosee As Double
cose = Cos(ee)
sine = Sin(ee)
cosee = Cos(ee) - e
Dim R As Double '定义彗日距离，单位为A.U.
R = longradii * (1 - e * cose)
'定义日心赤道坐标(x',y',z')
Dim lx As Double, ly As Double, lz As Double
lx = ax * cosee + bx * sine
ly = ay * cosee + by * sine
lz = az * cosee + bz * sine
'定义X,Y,Z太阳直角坐标变量 需要从天文年历中提取数据
'计算公式
Dim xx As Double, yy As Double, zz As Double
xx = 0.35782
yy = -0.84049
zz = -0.3644
'定义赤经,赤纬(α,δ)变量
Dim alpha As Double, delta As Double
alpha = Atn((ly + yy) / (lx + xx))
delta = Atn((lz + zz) * Sin(alpha) / (ly + yy))
Select Case cType
Case 1
ComputePosition = alpha * 180 / pi '角度值
Case 2
ComputePosition = delta * 180 / pi '角度值
End Select
End Function
Public Function days360(dt1 As Date, dt2 As Date) As Long '计算两个日期间的天数
Dim z1 As Long, z2 As Long
Dim d1 As Long, d2 As Long
Dim m1 As Long, m2 As Long
Dim y1 As Long, y2 As Long
d1 = day(dt1)
m1 = month(dt1)
y1 = year(dt1)
d2 = day(dt2)
m2 = month(dt2)
y2 = year(dt2)
If d1 = 31 Then
z1 = 30
Else
z1 = d1
End If
If d2 = 31 And d1 >= 30 Then
z2 = 30
Else
z2 = d2
End If
days360 = (y2 - y1) * 360 + (m2 - m1) * 30 + (z2 - z1)
End Function
页:
[1]