寿星万年历---Lua实现

OStack程序员社区-中国程序员成长平台 › 门户 › 编程› LUA›LUA教程

原作者: [db:作者] 来自: [db:来源] 收藏邀请

由于最近刚好准备用lua脚本去写一个软件，需要用到非常精确的节气计算时间，寿星万年历是我在网上见到的一份极高高精度的万年历，其采用先进的算法实现，其精度堪比刘安国教授为中国科学院国家授时中心制作的日梭万年历。但网络上只有javascript版本。于是自己将其翻译为lua脚本程序，并公布于此，方便大家使用。

仅保留了节气计算部分，省略了一些塑望月计算、阴历计算的内容。

我之前的一篇文章：寿星万年历---java算法实现，是我用java实现的一个版本。

寿星万年历相关信息：http://www.fjptsz.com/xxjs/xjw/rj/113.htm 设计:许剑伟

日梭万年历相关信息：http://www.time.ac.cn/calendar/calendar.htm 作者：刘安国

关于寿星万年历相关简要描述如下：

寿星万年历是一款采用现代天文算法制作的农历历算程序，含有公历与回历信息，可以很方便的进行公、农、回三历之间的转换。提供公元-4712年到公元9999年的日期查询功能。其中1500年到1940农历数据已经与陈垣的《二十史朔闰表》核对；含有从公420元(南北朝/宋武帝元年)到今的基本年号。在过去几百年中，寿星万年历的误差是非常小的，节气时刻计算及日月合朔时刻的平均误差小于1秒，太阳坐标的最大可能误差为0.2角秒，月亮坐标的最大可能误差为3角秒，平均误差为误差的1/6。万年历中含有几百个国内城市的经纬度，并且用户可根据自已的需要扩展经纬度数据。

以下是本人用lua实现的代码：

--*******农历节气计算部分
--========角度变换===============
local rad = 180*3600/math.pi --每弧度的角秒数
local RAD = 180/math.pi      --每弧度的角度数
function int2(v) --取整数部分
    v=math.floor(v)
    if v<0 then return v+1
    else return v
    end
end

function rad2mrad(v)   --对超过0-2PI的角度转为0-2PI
    v=math.fmod(v ,2*math.pi)
    if v<0 then  return v+2*math.pi
    else return v
    end
end

function rad2str(d,tim) --将弧度转为字串
    ---tim=0输出格式示例: -23°59' 48.23"
    ---tim=1输出格式示例:  18h 29m 44.52s
    local s="+"
    local w1="°" w2="’"  w3="”"
    if d<0 then  d=-d  s='-'end
    if tim~= 0 then  d=d*12/math.pi w1="h " w2="m " w3="s "
    else     d=d*180/math.pi end
    local a=math.floor(d) d=(d-a)*60
    local b=math.floor(d) d=(d-b)*60
    local c=math.floor(d) d=(d-c)*100
    d=math.floor(d+0.5)
    if d>=100 then d=d-100 c=c+1 end
    if c>=60  then c=c-60  b=b+1 end
    if b>=60  then b=b-60  a=a+1 end
    a="   "+a b="0"+b c="0"+c d="0"+d
    local alen = string.len(a)
    local blen = string.len(b)
    local clen = string.len(c)
    local dlen = string.len(d)
    s = s .. string.sub(a, alen-3,alen)+w1
    s = s .. string.sub(b, blen-2,blen)+w2
    s = s .. string.sub(c, clen-2,clen)+"."
    s = s .. string.sub(d, dlen-2,dlen)+w3
    return s
end
--================日历计算===============
local J2000=2451545 --2000年前儒略日数(2000-1-1 12:00:00格林威治平时)

local JDate={ --日期元件
Y=2000, M=1, D=1, h=12, m=0, s=0,
dts = { --世界时与原子时之差计算表
-4000,108371.7,-13036.80,392.000, 0.0000, -500, 17201.0,  -627.82, 16.170,-0.3413,
-150, 12200.6,  -346.41,  5.403,-0.1593,  150,  9113.8,  -328.13, -1.647, 0.0377,
500,  5707.5,  -391.41,  0.915, 0.3145,  900,  2203.4,  -283.45, 13.034,-0.1778,
1300,   490.1,   -57.35,  2.085,-0.0072, 1600,   120.0,    -9.81, -1.532, 0.1403,
1700,    10.2,    -0.91,  0.510,-0.0370, 1800,    13.4,    -0.72,  0.202,-0.0193,
1830,     7.8,    -1.81,  0.416,-0.0247, 1860,     8.3,    -0.13, -0.406, 0.0292,
1880,    -5.4,     0.32, -0.183, 0.0173, 1900,    -2.3,     2.06,  0.169,-0.0135,
1920,    21.2,     1.69, -0.304, 0.0167, 1940,    24.2,     1.22, -0.064, 0.0031,
1960,    33.2,     0.51,  0.231,-0.0109, 1980,    51.0,     1.29, -0.026, 0.0032,
2000,    64.7,    -1.66,  5.224,-0.2905, 2150,   279.4,   732.95,429.579, 0.0158, 6000},
deltatT = function(JDate, y) --计算世界时与原子时之差,传入年
    local  i
    local d=JDate.dts
    for x=1,100, 5 do
        if y<d[x+5] or x==96 then  i=x break end
    end

    local t1=(y-d[i])/(d[i+5]-d[i])*10
    local t2=t1*t1
    local t3=t2*t1
    return d[i+1] +d[i+2]*t1 +d[i+3]*t2 +d[i+4]*t3
end,
deltatT2 = function(JDate,jd) --传入儒略日(J2000起算),计算UTC与原子时的差(单位:日)
    return JDate:deltatT(jd/365.2425+2000)/86400.0
end,
toJD = function(JDate, UTC) --公历转儒略日,UTC=1表示原日期是UTC
    local  y=JDate.Y m=JDate.M n=0 --取出年月
    if m<=2 then  m=m+12 y=y-1 end
    if JDate.Y*372+JDate.M*31+JDate.D>=588829 then --判断是否为格里高利历日1582*372+10*31+15
        n =int2(y/100) n =2-n+int2(n/4)--加百年闰
    end
    n = n + int2(365.2500001*(y+4716))    --加上年引起的偏移日数
    n = n + int2(30.6*(m+1))+JDate.D       --加上月引起的偏移日数及日偏移数
    n = n + ((JDate.s/60+JDate.m)/60+JDate.h)/24 - 1524.5
    if(UTC == 1) then return n+JDate.deltatT2(n-J2000) end
    return n
end,
setFromJD = function(JDate, jd,UTC) --儒略日数转公历,UTC=1表示目标公历是UTC
    if UTC==1 then  jd= jd - JDate:deltatT2(jd-J2000) end
    jd =jd+0.5
    local A=int2(jd) F=jd-A, D  --取得日数的整数部份A及小数部分F
    if A>2299161 then  D=int2((A-1867216.25)/36524.25) A=A+1+D-int2(D/4) end
    A = A + 1524 --向前移4年零2个月
    JDate.Y =int2((A-122.1)/365.25)--年
    D =A-int2(365.25*JDate.Y) --去除整年日数后余下日数
    JDate.M =int2(D/30.6001)       --月数
    JDate.D =D-int2(JDate.M*30.6001)--去除整月日数后余下日数
    JDate.Y=JDate.Y-4716 JDate.M=JDate.M-1
    if JDate.M>12 then  JDate.M=JDate.M - 12 end
    if JDate.M<=2 then  JDate.Y = JDate.Y+1 end
    --日的小数转为时分秒
    F=F*24 JDate.h=int2(F) F=F - JDate.h
    F=F*60 JDate.m=int2(F) F=F - JDate.m
    F=F*60 JDate.s=F
end,

setFromStr = function(JDate, s) --设置时间,参数例:"20000101 120000"或"20000101"
    JDate.Y=string.sub(s, 1,4)    JDate.M=string.sub(s, 5, 6)  JDate.D=string.sub(s,7, 8)
    JDate.h=string.sub(s, 10, 11) JDate.m=string.sub(s, 12,13) JDate.s=string.sub(s, 14,18)
end,
toStr = function(JDate) --日期转为串
    local Y="     " .. JDate.Y
    local M="0" .. JDate.M
    local D="0" .. JDate.D
    local h=JDate.h
    local m=JDate.m
    local s=math.floor(JDate.s+.5)
    if s>=60 then s=s-60 m=m+1 end
    if m>=60 then m=m-60 h=h+1 end
    h="0".. h m="0" .. m s="0" .. s
    local Ylen = string.len(Y)
    local Mlen = string.len(M)
    local Dlen = string.len(D)
    local hlen = string.len(h)
    local mlen = string.len(m)
    local slen = string.len(s)
    Y=string.sub(Y, Ylen -4,Ylen) M=string.sub(M, Mlen-1,Mlen) D=string.sub(D,Dlen-1, Dlen)
    h=string.sub(h, hlen-1, hlen) m=string.sub(m, mlen-1,mlen) s=string.sub(s, slen-1,slen)
    return Y .. "-" .. M .. "-" .. D .. " " .. h .. ":" .. m .. ":" .. s
end,

JQ = function(JDate) --输出节气日期的秒数
    local t = {}
    t.year=JDate.Y
    t.month=JDate.M
    t.day=JDate.D
    t.hour=JDate.h
    t.min=JDate.m
    t.sec=math.floor(JDate.s+.5)
    if t.sec>=60 then t.sec=t.sec-60 t.min=t.min+1 end
    if t.min>=60 then t.min=t.min-60 t.hour=t.hour+1 end
    return os.time(t)
end,

Dint_dec = function(JDate, jd,shiqu,int_dec) --算出:jd转到当地UTC后,UTC日数的整数部分或小数部分
    --基于J2000力学时jd的起算点是12:00:00时,所以跳日时刻发生在12:00:00,这与日历计算发生矛盾
    --把jd改正为00:00:00起算,这样儒略日的跳日动作就与日期的跳日同步
    --改正方法为jd=jd+0.5-deltatT+shiqu/24
    --把儒略日的起点移动-0.5(即前移12小时)
    --式中shiqu是时区,北京的起算点是-8小时,shiqu取8
    local u=jd+0.5-JDate.deltatT2(jd)+shiqu/24
    if int_dec~= 0 then  return math.floor(u) --返回整数部分
    else return u-math.floor(u)      --返回小数部分
    end
end,

d1_d2 = function(JDate, d1,d2) --计算两个日期的相差的天数,输入字串格式日期,如:"20080101"
    local Y=JDate.Y M=JDate.M D=JDate.D h=JDate.h m=JDate.m s=JDate.s --备份原来的数据
    JDate.setFromStr(string.sub(d1,1,8)+" 120000")    local jd1=JDate.toJD(0)
    JDate.setFromStr(string.sub(d2,1,8)+" 120000")    local jd2=JDate.toJD(0)

    JDate.Y=Y JDate.M=M JDate.D=D JDate.h=h JDate.m=m JDate.s=s --还原
    if jd1>jd2 then  return  math.floor(jd1-jd2+.0001)
    else        return -Math.floor(jd2-jd1+.0001)
    end
end,
}
--=========黄赤交角及黄赤坐标变换===========
local hcjjB = {84381.448, -46.8150, -0.00059, 0.001813}--黄赤交角系数表
local preceB= {0,50287.92262,111.24406,0.07699,-0.23479,-0.00178,0.00018,0.00001}--Date黄道上的岁差p

function hcjj1 (t) --返回黄赤交角(常规精度),短期精度很高
    local t1=t/36525 t2=t1*t1  t3=t2*t1
    return (hcjjB[1] +hcjjB[2]*t1 +hcjjB[3]*t2 +hcjjB[4]*t3)/rad
end

function HCconv(JW,E) --黄赤转换(黄赤坐标旋转)
    --黄道赤道坐标变换,赤到黄E取负
    local HJ=rad2mrad(JW[1])  HW=JW[2]
    local sinE =math.sin(E) cosE =math.cos(E)
    local sinW=cosE*math.sin(HW)+sinE*math.cos(HW)*math.sin(HJ)
    local J=math.atan2( math.sin(HJ)*cosE-math.tan(HW)*sinE, math.cos(HJ) )
    JW[1]=rad2mrad(J)
    JW[2]=math.asin(sinW)
end

function addPrece(jd,zb) --补岁差
    local i t=1 v=0  t1=jd/365250
    for i=2,8 do t=t*t1 v=v+preceB[i]*t end
    zb[1]=rad2mrad(zb[1]+(v+2.9965*t1)/rad)
end

--===============光行差==================
local GXC_e={0.016708634, -0.000042037,-0.0000001267} --离心率
local GXC_p={102.93735/RAD,1.71946/RAD, 0.00046/RAD}  --近点
local GXC_l={280.4664567/RAD,36000.76982779/RAD,0.0003032028/RAD,1/49931000/RAD,-1/153000000/RAD} --太平黄经
local GXC_k=20.49552/rad --光行差常数
function addGxc(t,zb)--恒星周年光行差计算(黄道坐标中)
    local t1=t/36525
    local t2=t1*t1
    local t3=t2*t1
    local t4=t3*t1
    local L=GXC_l[1] +GXC_l[2]*t1 +GXC_l[3]*t2 +GXC_l[4]*t3 +GXC_l[5]*t4
    local p=GXC_p[1] +GXC_p[2]*t1 +GXC_p[3]*t2
    local e=GXC_e[1] +GXC_e[2]*t1 +GXC_e[3]*t2
    local dL=L-zb[1]
    local dP=p-zb[1]
    zb[1]=zb[1] - (GXC_k * (math.cos(dL)-e*math.cos(dP)) / math.cos(zb[2]))
    zb[2]=zb[2] - (GXC_k * math.sin(zb[2]) * (math.sin(dL)-e*math.sin(dP)))
    --print('aa', L,p,e,dL,dP,zb[1], zb[2])
    zb[1]=rad2mrad(zb[1])
end

--===============章动计算==================
local nutB={--章动表
2.1824391966,   -33.757045954, 0.0000362262, 3.7340E-08,-2.8793E-10,-171996,-1742,92025, 89,
3.5069406862,  1256.663930738, 0.0000105845, 6.9813E-10,-2.2815E-10, -13187,  -16, 5736,-31,
1.3375032491, 16799.418221925,-0.0000511866, 6.4626E-08,-5.3543E-10,  -2274,   -2,  977, -5,
4.3648783932,   -67.514091907, 0.0000724525, 7.4681E-08,-5.7586E-10,   2062,    2, -895,  5,
0.0431251803,  -628.301955171, 0.0000026820, 6.5935E-10, 5.5705E-11,  -1426,   34,   54, -1,
2.3555557435,  8328.691425719, 0.0001545547, 2.5033E-07,-1.1863E-09,    712,    1,   -7,  0,
3.4638155059,  1884.965885909, 0.0000079025, 3.8785E-11,-2.8386E-10,   -517,   12,  224, -6,
5.4382493597, 16833.175267879,-0.0000874129, 2.7285E-08,-2.4750E-10,   -386,   -4,  200,  0,
3.6930589926, 25128.109647645, 0.0001033681, 3.1496E-07,-1.7218E-09,   -301,    0,  129, -1,
3.5500658664,   628.361975567, 0.0000132664, 1.3575E-09,-1.7245E-10,    217,   -5,  -95,  3}

function nutation(t) --计算黄经章动及交角章动
    local d={}
    d.Lon=0  d.Obl=0  t=t/36525
    local i,c
    local t1=t
    local t2=t1*t1
    local t3=t2*t1
    local t4=t3*t1
    local t5=t4*t1
    for i=1,#nutB,9 do
        c=nutB[i] +nutB[i+1]*t1 +nutB[i+2]*t2 +nutB[i+3]*t3 +nutB[i+4]*t4
        d.Lon=d.Lon + (nutB[i+5]+nutB[i+6]*t/10)*math.sin(c) --黄经章动
        d.Obl=d.Obl + (nutB[i+7]+nutB[i+8]*t/10)*math.cos(c) --交角章动
    end
    d.Lon=d.Lon/(rad*10000) --黄经章动
    d.Obl=d.Obl/(rad*10000) --交角章动
    return d
end

function nutationRaDec(t,zb) --本函数计算赤经章动及赤纬章动
    local Ra=zb[1]
    local Dec=zb[2]
    local E=hcjj1(t)
    local sinE=math.sin(E)
    local cosE=math.cos(E) --计算黄赤交角
    local d=nutation(t)                                  --计算黄经章动及交角章动
    local cosRa=math.cos(Ra)
    local sinRa=math.sin(Ra)
    local tanDec=math.tan(Dec)
    zb[1]=zb[1] + (cosE+sinE*sinRa*tanDec)*d.Lon-cosRa*tanDec*d.Obl --赤经章动
    zb[2]= zb[2] + sinE*cosRa*d.Lon+sinRa*d.Obl   --赤纬章动
    zb[1]=rad2mrad(zb[1])
end

--=================以下是月球及地球运动参数表===================
--[[***************************************
* 如果用记事本查看此代码,请在"格式"菜单中去除"自动换行"
* E10是关于地球的,格式如下:
*    它是一个数组,每3个数看作一条记录,每条记录的3个数记为A,B,C
*    rec=A*cos(B+C*t)  式中t是J2000起算的儒略千年数
*    每条记录的计算结果(即rec)取和即得地球的日心黄经的周期量L0
* E11格式如下: rec = A*cos*(B+C*t) *t,     取和后得泊松量L1
* E12格式如下: rec = A*cos*(B+C*t) *t*t,   取和后得泊松量L2
* E13格式如下: rec = A*cos*(B+C*t) *t*t*t, 取和后得泊松量L3
* 最后地球的地心黄经:L = L0+L1+L2+L3+...
* E20,E21,E22,E23...用于计算黄纬
* M10,M11等是关于月球的,参数的用法请阅读Mnn()函数
***************************************** --]]
--地球运动VSOP87参数
local E10={ --黄经周期项
1.75347045673, 0.00000000000,     0.0000000000,  0.03341656456, 4.66925680417,  6283.0758499914,  0.00034894275, 4.62610241759, 12566.1516999828,  0.00003417571, 2.82886579606,     3.5231183490,
0.00003497056, 2.74411800971,  5753.3848848968,  0.00003135896, 3.62767041758, 77713.7714681205,  0.00002676218, 4.41808351397,  7860.4193924392,  0.00002342687, 6.13516237631,  3930.2096962196