节气交节时刻计算

节气交节时刻计算

我是编程爱好者，我的第一个设计的程序是万年历。程序的公历部分很容易做，只要计算月首1日是星期几就能打印出月历。问题来了，光有公历的万年历怎么拿得出手。于是就寻找添加农历的算法。等打印出有农历信息的万年历，发现没有节气信息。不管是挂历还是台历都有节气时间信息。我也要有。于是就开始了研究节气时间计算。节气在哪天比较容易计算。而节气交节时刻计算就不是一般的难了。

期间我研究过许多方法，有曾次亮先生的节气算法改成公式计算，有用数据表取数据计算，有用C++ VSOP87的天文算法，有许剑伟的天文算法等等。
天文算法精度很高但太复杂。数据表方法数据量太大，做个200年万年历还行。简单点的办法是采用天文算法的太阳黄经度数的简化算法。

二十四节气是中国人的独特创造，它根据太阳在黄道上的位置（黄经）将全年分为24个段落。由于二十四节气是根据太阳直射地球的某些纬度而确定的，所以它属于阳历的范畴，每个节气与阳历的日期基本对应。例如，春分、秋分，黄道和赤道平面相交，此时黄经分别为0度、180度，太阳直射赤道，昼夜相等。夏至，太阳直射北纬23.5度，黄经90度，北半球白昼最长。冬至，太阳直射南纬23.5度，黄经270度，北半球白昼最短。春分和秋分（二分）正处春秋两季中间，夏至和冬至（二至）正处夏冬两季中间。这样，一年就可用春分、夏至、秋分、冬至划为四段。如将每段再分成6小段，太阳在黄道上移动15度，每小段约15天左右，全年就可分为24小段，即24个节气，每15天为一个“节气”。

下面的源码是用简单C语言写的。myspringc编译器v2.7编译调试。可用该编译器做成安卓手机桌面应用程序。

附言：另有节气朔望计算，日食月食预测计算源码。需要的可联系我。

//*******************************************
//**** 节气交节时刻计算
//**** 简单C语言 myspringc v2.7 编译
//**** 张纯叔 (micelu@126.com)
//********************************************
string sBarDes[10];
int nBarId[10];
float pi;
string s; //return string
int n,i;
int dy,dm,dd; //date
double dy1,dm1,dd1;
string dy2,dm2,dd2; //output date yy-mm-dd
int hh,mm,ss; //time
double hh1,mm1,ss1;
string hh2,mm2,ss2; //output times hh:mm:ss
double sw; //sun 黄经
double jd; //儒略日
string jqnames; //节气名称
int sv; //设置求节气黄经度数
string dat$;//retunesolardat$
string jqs1,jqs2; //节气，中气
string jnum; //节气编码

main(){
pi=3.1415926535;
jqnames=“小寒大寒立春雨水惊蛰春分清明谷雨立夏小满芒种夏至小暑大暑立秋处暑白露秋分寒露霜降立冬小雪大雪冬至”;
sBarDes[0]="求? 节气 ";
nBarId[0]=100;
sBarDes[3]=“全年节气表”;
nBarId[3]=103;
sBarDes[4]=“退出程序”;
nBarId[4]=104;
setToolBarHeight(10);
setButtonTextSize(15);
setToolBarBackgroundColor(255,192,192,192);
setButtonColor(255,0,0,240);
setButtonTextColor(255,255,245,0);
setToolBar(100,myToolBarProc,sBarDes,nBarId,5);
setTitle(“节气交节时刻计算”);
dy=2000; askjq ();
while (){}
}//main ()

getdate (){ //输入日期
int ds[3];
string sDat[101];
getDate(ds);
pickDate(“输入日期:”,ds);
sDat[n]=intToString(ds[0])+"-"+intToString(ds[1])+"-"+intToString(ds[2]);
clearOutput ();
dy=ds[0];
dm=ds[1];
dd=ds[2];
print “InputDate = “,dy,”-”,dm,"-",dd;
//*** 日期转儒略日 GDtoJD
jd=(int)(365.25*(dy+4716)+0.01)+(int)(30.6001*(dm+1))+dd-10-1524.5-0.5;
print “计算结果如下：”;
print "jd = ",jd;
}//getdate ()

myToolBarProc(int nBtn,int nContext){
if(nBtn100){//输入日期，计算儒略日
getdate();
print " ";
ask ();
}
if(nBtn103){//计算全年节气
getdate ();
askjq(); // this year all solarterm
}
if(nBtn==104){//退出程序
clearOutput();
setDisplay (0);
exit (0);
}
}//my toolbar()

askjq (){//calculate and show year 24 solarterm
//打印全年节气表 ( 1900–2100 )
if (dy<1900)dy=2000;
clearOutput ();
print " ";
print " “,dy,” 年 节气交节时刻表 ";
print " ";
for (dm=1;dm<13;dm++){
ask (); }
}//askjq ()

ask (){ //input dm calculate solarterm ********
//输入年份月份 dy-dm 计算节气时间
n=dm*2-1;
sv=(n-1)15+285;
if (sv>360)sv=sv-360;
jnum=intToString(n);
cal_jq ();
print s; //打印节气
n=dm2;
sv=(n-1)*15+285;
if (sv>359.9)sv=sv-360;
jnum=intToString (n);
cal_jq ();
print s; //打印中气
}//ask()

cal_jq (){
//求？节气时刻：二分法求直线方程根
jd=dy* (365.2423112-0.000000000000064*(dy-100)(dy-100)- 0.00000003047(dy-100))+15.218427*n+1721050.71301;
for (i=0;i<10;i++){ //试算10次内得解
cal_sun ();
if (sw<0)sw=sw+360;
if (sw>sv)jd=jd-(sw-sv)*0.96;
if (sw<sv)jd=jd+(sv-sw)*0.96;
if (absd(sw-sv)<0.00001)break;
}
jd=jd+0.333333; //UTC+8小时 >> 北京时间
jdtoGD ();
jqs1=subString (jqnames,(n-1)*2,2);
if (sw<sv)sw=sw+0.00001;
string sws; //节气黄经度数
sws=doubleToString (sw);
if (sv<100)sws=" “+sws;
if (n<10)jnum=” “+jnum;
s=jnum+” “+jqs1+” " +dat$+" "+sws;
//return solar date and times
}//cal_jq ()

cal_sun(){//input jd, calculate solarterm
//输入儒略日 jd 计算太阳黄经
double t,t2,t3,t4;
double Lo;
double m,e,c,th;
double ang;
double v,R; //太阳黄经计算 …
ang=pi/180;
//太阳几何平黄经：
t=(jd-2451545.0)/36525;
t2=tt;
t3=ttt;
t4=tttt;
Lo=280.46645+36000.76983t+0.0003032t2;
//（Date平分点起算）
Lo=Lo-(int)(Lo/360)*360;
//print "太阳几何平黄经 Lo = ",Lo;

//太阳平近点角:
m=357.5291+35999.0503t-0.0001559t2-0.00000048*t3;
m=m-(int)(m/360)*360;

//地球轨道离心率:
e=0.016708617-0.000042037t-0.0000001236t2;
//太阳中间方程:
c=(1.9146-0.004817* t-0.000014t2)sin(angm)+ (0.019993-0.000101t)sin(2angm)+0.00029 sin(3angm);
//那么，太阳的真黄经是：
th=Lo+c;
//真近点角是:
v = m + c;
//日地距离的单位是"天文单位",距离表达为：
R=1.000001018*(1-ee)/(1+ecos(ang*v));

//print "太阳平近点角 m = ",m;
//print "地球轨道离心率 e = ",e;
//print "太阳中间方程 c = ",c;
//print "太阳真黄经 th = ",th;
//print "真近点角 v = ", v;
//print "日地距离 R = ",R;

//章动修正
double zd;
zd=125.04-1934.136t; // '章动修正
sw=th-0.00569-0.00478sin(ang*zd); //光行差修正
if (sw>360)sw=sw-360;
if (sw<0)sw=sw+360;
//print "太阳视黄经 sw = ",sw;
}//cal_sun ()

jdtoGD(){ //return dat$ ********
// jd to GD 儒略日转公历日期
int a,b,c,d,e;
double F;
double allss;
F=jd-(int)(jd);
//print " "; //调试
//print " JD = ",jd;
//print " 时分秒 日小数 = ",F;

a=(int)(jd+0.5);
b=a+1537;
c=(int)((b-122.1)/365.25);
d=(int)(365.25c);
e=(int)((b-d)/30.6001);
dd=b-d-(int)(30.6001e);
dm=e-1-(int)((e/14)*12);
dy1=c-4715-(int)((7+dm)/10);
//print a," “,b,” “,c,” “,d,” “,e; //调试
dy2=intToString (dy);
dm2=intToString (dm);
dd2=intToString (dd);
if(dm<10)dm2=“0”+dm2;
if(dd<10)dd2=“0”+dd2;
dm2=subString (dm2,0,2);
dd2=subString (dd2,0,2);
//print dy2+” 年 “+dm2+” 月 “+dd2+” 日 ";

//日allss 的小数转为时分秒 ************
allss=(int)((jd-a)86400+43200.5);
//print "allss = ", allss; //调试
hh1=(int)(allss/3600);
mm1=(int)((allss-hh13600)/60);
ss1=(int)(allss-hh13600-mm160);
if(ss1>=60){
ss1=ss1-60;
mm1=mm1+1;}
if(mm1>=60){
mm1=mm1-60;
hh1=hh1+1;}
//print "JD 转为 GD，计算结果：”;
hh2=doubleToString(hh1);
mm2=doubleToString(mm1);
ss2=doubleToString(ss1);
//时间格式输出 hh:mm:ss 00:00:00
if(hh1<10){
hh2=“0”+doubleToString(hh1);}
if(mm1<10){
mm2=“0”+doubleToString(mm1);}
if(ss1<10){
ss2=“0”+doubleToString(ss1);}
hh2=subString (hh2,0,2);
mm2=subString (mm2,0,2);
ss2=subString (ss2,0,2);

dat$=dm2+"-"+dd2+""+hh2+":"+mm2+":"+ss2;//returndat$
}//jdtoGD() **************
//**** End *************

//**** 另附曾次亮先生节气算法 ********
cal_jqz (){ //input n to calculate solarterm
double juD,tht,yrD,shuoD,vs,dalT;
juD=dy*(365.2423112-0.000000000000064*(dy-100)(dy-100)- 0.00000003047(dy-100))+15.218427*n+1721050.71301;

tht=0.0003dy-0.372781384-0.2617913325n;

yrD=(1.945sin(tht)-0.01206sin(2tht))(1.048994-0.00002583 *dy);

shuoD=-0.0018sin(2.313908653dy-0.439822951-3.0443n);
//‘vs = juD //’** 平气
vs = (juD + yrD + shuoD); // '** 定气
dalT=-15+(jd-2382148)(jd-2382148)/41048480;
dalT = dalT / 86400;
jd=vs-0.5-dalT;
jdtoGD ();
jqs1=subString (jqnames,(n-1)*2,2);
if (sw<sv)sw=sw+0.00001;
string sws; //节气黄经度数
sws=doubleToString (sw);
if (sv<100)sws=" “+sws;
if (n<10)jnum=” “+jnum;
s=jnum+” “+jqs1+” " +dat$+" "+sws;
｝