节气朔望时刻计算和日食月食预测

节气朔望时刻计算和日食月食预测

中国农历的朔望月是农历历法的基础，而朔望月又是严格以日月合朔发生的那一天作为月首，因此日月合朔时间的计算是制定农历历法的关键。
二十四节气是中国人的独特创造，它根据太阳在黄道上的位置（黄经）将全年分为24个段落。由于二十四节气是根据太阳直射地球的某些纬度而确定的，所以它属于阳历的范畴，每个节气与阳历的日期基本对应。例如，春分、秋分，黄道和赤道平面相交，此时黄经分别为0度、180度，太阳直射赤道，昼夜相等。夏至，太阳直射北纬23.5度，黄经90度，北半球白昼最长。冬至，太阳直射南纬23.5度，黄经270度，北半球白昼最短。春分和秋分（二分）正处春秋两季中间，夏至和冬至（二至）正处夏冬两季中间。这样，一年就可用春分、夏至、秋分、冬至划为四段。如将每段再分成6小段，太阳在黄道上移动15度，每小段约15天左右，全年就可分为24小段，即24个节气，每15天为一个“节气”。这就是节气交节时刻的数学定义。
要计算日月合朔时间，首先要对日月合朔这一天文现象进行数学定义。朔望月是在地球上观察到的月相周期，平均长度约等于29.53059日，而恒星月（天文月）是月亮绕地球公转一周的时间，长度约27.32166日。月相周期长度比恒星月长大约两天，这是因为在月球绕地球旋转一周的同时，地球还带着它绕太阳旋转了一定的角度的缘故，所以月相周期不仅与月球运行有关，还和太阳运行有关。日月合朔的时候，太阳、月亮和地球三者接近一条直线，月亮未被照亮的一面对着地球，因此地球上看不到月亮，此时又被称为新月。

月亮绕太阳公转的白道面和地球绕太阳公转的黄道面存在一个最大约5°的夹角，因此大多数情况下，日月合朔时都不是严格在同一条直线上，不过也会发生在同一直线的情况，此时就会发生日食。
日月合朔的数学定义就是太阳和月亮的地心视黄经差为0的时刻。
另外一种情况是在太阳的黄道面上的另一边，这时太阳和月亮的地心视黄经差为180度，这就是望的时刻。俗称月亮最圆最亮的时候。如果此刻在黄道面上的夹角小于5度，就可能发生月食。
要计算日月合朔时刻，望时刻，节气交节时刻，其实就是计算太阳地心视黄经和月亮地心视黄经。计算方法有古代天文算法和近现代天文算法，还有各种民间的算法。

我研究过曾次亮先生的《四千年气朔交食速算法》。他是把数据放表格里，查表找数计算，类似于以前没计算器时用的6位数学用表。能算出结果，但是时间误差太大。精确的算法是现代理论天文算法，是外国人的，网上介绍有许剑伟老师的java算法，还有c++ 语言编写的。精度很高，但太复杂，使初学者望而却步。
本文下面介绍的是简化的天文算法。比较通俗易懂，精度差点，误差5分钟。能计算预测出日食月食时间。对于精度要求不高的可参考一下。
计算的要点是：已知黄经度数如何反求时间。
我的方法很简单直观，在已知度数的正负两边加减调整，逼近正解。
for (i=0;i<10;i++){ //试算5次内求得解
cal_sun (); //求太阳黄经sw , 已知节气黄经度 sv
if (sw<0)sw=sw+360;
if (sw>sv)jd=jd-(sw-sv);
if (sw<sv)jd=jd+(sv-sw);
if (absd(sw-sv)<0.00001)break; //计算精度
//print sw; }
朔望时刻求解方法类此法。

关于日月食无需太多的计算，就可精确获得日月食的主要特点。日食情况是比较复杂，在地面上不同的观测位置，事件发生的时间是不同的。而对于月食，所有的观测者将在同一时刻看到相同的月相。
由于这个原因，我们将不考虑各地的日食情况。有兴趣的读者可以使用贝塞尔根计算各地的日食情况。
计算的一般数据是这样：首先，利用公式(49.1)和(47.3)计算平新月和平满月时刻(JDE)。记住，对于新月(日食)，k必须是一个整数，对应满月(月食)，k为整数加0.5。然后，利用公式(47.4)到(47.7)计算此时的角度M、M′、F和Ω，并利用公式(45.6)式计算E。
F的值将给出发生日食或月食的首要信息。如果F与整倍数的180度的差值小于13°.9，那么肯定要发生日月食，如果大于21°.0，那就没有日月食。如果在这两个值之间，则不能确定是否有日月食，还须进一步检查。可以使用以下规则：如果|sin F|>0.36，则没有日月食。
注意，在一个阴历月以后，F增加了30°.6705。
如果F接近0°或360°，那么日月食出现在升交点。如果F接近180度，那么食点出现在月亮轨道的降交点。

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

//*******************************************
//**** 计算： 节气时刻 朔望时刻 日食月食 预测
//**** 简单的C语言 myspringc v2.7 编译调试
//**** 张纯叔 (micelu@126.com)
//********************************************
//**** myspringc v2.7 源码 solar.c
Canvas cs; //绘图画布，图略
string sBarDes[10]; //menu & bar
int nBarId[10];
float src[4]; //ClearDraw (cls) clear screen
float pi; //draw mode picture
int i,k;
int n;
string s; //临时string
string s0,s1,s2;
int dy,dm,dd; //date yy-mm-dd
double dy1,dm1,dd1;
string dy2,dm2,dd2;
int hh,mm,ss; //time hh:mm:ss
double hh1,mm1,ss1; //time hh:mm:ss
string hh2,mm2,ss2; //output times $
double t,t2,t3,t4; //calculate jd, t …
double ang,de;
double ml,mw,tw,sw; //moon,sun 黄经
double D,m,m1,f;
double a1,a2,a3;
double e,sum1; //参数 …
double Lo; // 黄经
double c,th;
double v,R; //黄经计算 …
double jd; //儒略日
int swd; //select show朔日望日
string jqnames; //节气名称
int sv; //节气黄经度数
string dat$;//retunesolardat$ = jname2;
string jqs1,jqs2; //2 solar in a month
string jnum; // solar number

main(){
pi=3.1415926535;
jqnames=“小寒大寒立春雨水惊蛰春分清明谷雨立夏小满芒种夏至小暑大暑立秋处暑白露秋分寒露霜降立冬小雪大雪冬至”;
de=pi/180;
sBarDes[0]=“求? 节气”;
nBarId[0]=100;
sBarDes[1]=“求? 朔日”;
nBarId[1]=101;
sBarDes[2]=“求? 望日”;
nBarId[2]=102;
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=2021, dm=5 ; //调试 有月全食
while (){}
}//main ()

getdate (){ //input date
int ds[3];
string sDat[101];
getDate(ds);
pickDate(“输入日期:”,ds); //安卓日期输入控件
sDat[n]=intToString(ds[0])+"-"+intToString(ds[1])+"-"+intToString(ds[2]);
clearOutput ();
//print ds[0],ds[1],ds[2];
dy=ds[0];
dm=ds[1];
dd=ds[2];
print “Input Date = “,dy,”-”,dm,"-",dd;
//*** 日期转儒略日 GDtoJD
//*** 此公式( 1900 - 2100 ) 有效，调整系数 -10 , -13 , 谨慎用
jd=(int)(365.25*(dy+4716)+0.01)+(int)(30.6001*(dm+1))+dd-10-1524.5-0.5;
print “计算结果如下：”;
print "jd = ",jd;
}//getdate ()

cal_shuo (){ //求？朔时刻 截弦逼近求直线方程根
swd=0;
// clearOutput ();
getdate (); //输入日期，转儒略日 jd
jdtoGD ();
print dat$;//打印输入日期ca{l}_{s}un();//求日黄经ca{l}_{m}oon();//求月黄经for(k=0;k<50;k++)//试算>>>jd=jd+(sw-mw)\ast 0.0722;if(sw>330)jd=jd+(sw-mw)\ast 0.0422;ca{l}_{s}un();ca{l}_{m}oon();//printsw,"",mw;if(absd(sw-mw)<0.00001)break;print"times=",k+1;print"求？朔日";print"日黄经=",sw;print"月黄经=",mw;print"黄经差=",sw-mw;jd=jd+0.333333;//转北京时间8/24print"jd=",jd;jdtoGD();printdat$; //print date$
cal_swd (); //判断日月食
}//cal_shuo ()

cal_wang (){ //求？望时刻
double sv0; //设定黄经度
swd=1;
clearOutput ();
getdate ();
jdtoGD ();
print dat$; //date $
cal_sun ();
cal_moon ();
for (k=0;k<50;k++){
sv0=absd(sw-mw);
sv0=absd(sv0-180);
jd=jd+sv0*0.0549;
cal_sun ();
cal_moon ();
//print sw," ",mw;
if(sw>mw)sv0=absd(sw-mw-180);
if (sw<mw)sv0=absd (mw-sw-180);
// print sv0;
if (sv0<0.00001)break;
}
print "times = ",k+1;
print "求？ 望日 ";
print "日黄经 = ",sw;
print "月黄经 = ",mw;
print "黄经差 = ",sw-mw;
jd=jd+0.333333; //转北京时间 8/24
print "jd = ",jd;
jdtoGD();
print dat$;//printdate$
cal_swd (); //判断日月食
}//cal_wang()

cal_swd(){ //判断日月食
double k;
double jde;
double t,t2,t3,t4;
double f,f1,m,m1,om;
double o,n1,h;
double a1,p,q,w,y,u;
double dt,d,bf,all,by;
print " ";
print "日月食计算判断: >>>>>> ";
jd=jd-0.333333;
print "jd = ",jd;
dy1=stringToDouble (dy2);
dm1=stringToDouble (dm2);
dd1=stringToDouble (dd2);
//print " “;
//print dy1,” - “,dm1,” - ",dd1;

k=(dy1+(dm1-1)/12+dd1/365-2000)*12.3685;
//print " ";
print "k = ",k;

if (swd0)k=(int)(k);
if (swd1)k=(int)(k)+0.5;
if (dy<2000)k=k-1;
print "k = ",k;

t=k/1236.85;
t2=tt;
t3=t2t;
t4=t3*t;
print "t = k/1236.85 = ",t;

jde=2451550.09765+29.530588953 * k+0.0001337t2-0.00000015t3
+0.00000000073t4; //(47.1)
if (jde-jd>12)jde=jde-354;
e=1-0.002516t-0.0000074*t2;
//print "jde = ",jde;

m=2.5534+29.10535669k-0.0000218t2
-0.00000011t3; //(47.3)
m=m-(int)(m/360)360;
if(m<0)m=m+360;
m1=201.5643+385.81693528k
+0.0107438t2+0.00001239t3
-0.000000058t4; //(47.4)
m1=m1-(int)(m1/360)*360;
if(m1<0)m1=m1+360;

f=160.7108+390.67050274k-0.0016341t2
-0.00000227t3+0.00000011t4; //(47.5)
f=f-(int)(f/360)*360;
if(f<0)f=f+360;
print "f = ",f;
print " ";

double df,df1;
if (f>200){df=360-f;}else{df=180-f;}
if (swd0)s=" (日食) ";
if (swd1)s=" (月食)";
print "df = “,df;
//print " “;
print " 判断： 如果 df < 13.9 则有食, df > 21 则无食”;
if (df<13.9){print " 判断结果： 有食”+s;}else{print " 判断结果：无食 ";}
print " ";
if (df>21)return;

if (swd0)s1="（日食）";
if (swd1)s1="（月食）";

if(df>13.9&&df<21){//check
print " DF>13.9 & DF<21 “+” 需核查 ";
print “核查：如果 | sin(f) | > 0.36，则没有日月食。”;
df1=absd (sin (f));
print “| sin(f) | = “,df1;
if (df1>0.36){print " 无食”;}else{print " 有食”+s;}
}//check

jdtoGD ();
if (swd0)s="朔日： ";
if (swd1)s="望日： ";
//print " ";
print s+dat$+"TD";//date$
//return;
print “····················································”;
//check *************
om=124.7746+(-1.56375580k ) +0.0020691t2+0.00000215*t3; //(47.7)
om=om-(int)(om/360)*360;
if (om<0)om=om+360;

f1=f-0.02665sin(om);
a1=290.77+0.107408k-0.009173t2;
p=0.2070esin(m)
+0.0024esin(2m)-0.0392sin(m1)
+0.0116sin(2m1)-0.0073esin(m1+m)
+0.0067esin(m1-m)+0.0118sin(2f1);
q=5.2207-0.0048ecos(m)
+0.0020ecos(2m)
-0.3299cos(m1)-0.0060ecos(m1+m)
+0.0041e*cos(m1-m);

w=absd(cos(f1));
y=(pcos(f1)+qsin(f1))(1-0.0048w);
u=0.0059+0.0046ecos(m1)
-0.0182cos(m1)+0.0004cos(2m1)
-0.0005cos(m+m1);

p=1.2985+u; //半影 p=1.2848+u;1.2985
o=0.7403-u; //本影
p=1.0157-u ; //p=1.0128-u; 1.0157
dt=0.4678-u;

n1=0.5458+0.0400cos(m1);
h=1.5573+u;
bf=(60/n1)sqrt(pp-yy); //部分被食
all=(60/n1)sqrt (dtdt-yy); //全部被食 sfa
by=(60/n1)sqrt(hh-yy); //半影食相 sfb
bf=60/n1;

double sfb,sfa;
if (swd==1){//望日 月食
if (bf>0)s=" 月食：部分食";
if (all>0)s=" 月食：全食";
sfb=(1.5573+u-absd (y))/0.5450;
sfa=(1.0128-u-absd(y))/0.5450;
print "半影食分 ",sfb;
print "本影食分 ",sfa;
}

//日食 : ? 如果 u<0，是全食
//如果 u>+0.0047，是环食
//如果 u介于0到+0.0047，是环食或全环食(混合食)
//w1=0.00464*(sqrt (1-yy))
//u<w1 全环食 ， u>w1 环食
//(1.5433+u - |y|) / (0.5461 + 2u ) (52.5)食分
double w1;
double sf;
if(swd==0){//日食
sf=(1.5433+u-absd (y))/(0.5461+2u);
w1=0.00464(sqrt(absd (1-y*y)));
if (w1>0){
if (u<0)s=" 日食：全食 “;
if (u>0.0047)s=” 日食：环食 “;
if (u<w1)s=” 日食：全环食 “;
if (u>w1)s=” 日食：环食 ";}
print "w1 = ",w1;
print s;
print "食分 = ",sf;
}
print " ";
print "e = ",e;
print "m = ",m;
print "m1 = ",m1;
print "f = ",f;
print "om = ",om;
print "f1 = ",f1;
print "a1 = ",a1;
print "p = ",p;
print "q = ",q;
print "w = ",w;
print "y = ",y;
print "u = ",u;
print "dt = ",dt;
print "n1 = ",n1;
print "h = ",h;
print “bf = “,(int)(bf),” 分钟”;
print “all = “,(int)(all),” 分钟”;
print “by = “,(int)(by),” 分钟”;
}//cal_swd ()

myToolBarProc(int nBtn,int nContext){
if(nBtn100){//输入日期，计算儒略日
// clearOutput ();
getdate();
print " ";
ask ();
}
if(nBtn101){//求？朔时刻
cal_shuo();
}
if(nBtn102){//求？望时刻
cal_wang ();
}
if(nBtn103){//计算全年节气
// clearOutput ();
getdate ();
askjq(); // this year all solarterm
}
if(nBtn==104){//退出程序
clearOutput();
cs.ClearDraw (0,src);
setDisplay (0);
exit (0);
}
}//my toolbar()

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

ask (){ //input dm calculate solarterm ********
n=dm*2-1;
sv=(n-1)15+285;
if (sv>360)sv=sv-360;
jnum=intToString(n);
cal_jq ();
//print sv;
print s;
//print jd;
n=dm2;
sv=(n-1)*15+285;
if (sv>359.9)sv=sv-360;
jnum=intToString (n);
cal_jq ();
//print sv;
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++){ //5次内求得解
cal_sun ();
if (sw<0)sw=sw+360;
if (sw>sv)jd=jd-(sw-sv);
if (sw<sv)jd=jd+(sv-sw);
if (absd(sw-sv)<0.00001)break;
//print sw;
}
jd=jd+0.333333; //转北京时间 8/24
jdtoGD ();
jqs1=subString (jqnames,(n-1)*2,2);
if (sw<sv)sw=sw+0.00001;
string sws;
//print “jd = “,jd;
sws=doubleToString (sw);
if (sv<100)sws=” “+sws;
if (n<10)jnum=” “+jnum;
s=jnum+” “+jqs1+” " +dat$+” "+sws;
//print s;
//return solar date and times
}//cal_jq ()

cal_sun(){//input jd, calculate solarterm
ang=pi/180;
//太阳几何平黄经：
t=(jd-2451545.0)/36525;
t2=tt;
t3=ttt;
t4=tttt;
//print " ";
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;
//print "月球黄经 mw = ",mw;
}//cal_sun ()

//calculate jd to GD ********************
jdtoGD(){ //return date $ ********
int a,b,c,d,e;
double F;
double allss;
F=jd-(int)(jd);
//print " ";
//print " JD = ",jd;
//print " 时分秒 日小数 = ",F;
//jd=jd+0.333;

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 " ";
//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);
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;//if(len(jnum)==1)jnum=""+jnum;s=dat$;
// return dat$
}//jdtoGD() **************

cal_moon (){//calculate moon return mw$
//月球平黄经:
ml = 218.3164591 + 481267.88134236 * t - 0.0013268 * t2 + t3 / 538841 - t4 / 65194000;
ml=ml-(int)(ml/360)*360;

//月日距角:
D = 297.8502042 + 445267.1115168 * t - 0.00163 * t2 + t3 / 545868 - t4 / 113065000;
D=D-(int)(D/360)*360;

//太阳平近点角:
m = 357.5291092 + 35999.0502909 * t - 0.0001536 * t2 + t3 / 24490000;
m=m-(int)(m/360)*360;

//月亮平近点角:
m1 = 134.9634114 + 477198.8676313 * t + 0.008997 * t2 + t3 / 69699 - t4 / 14712000;
m1=m1-(int)(m1/360)*360;

//月球经度参数 (到升交点的平角距离):
f = 93.2720993 + 483202.0175273 * t - 0.0034029 * t2 - t3 / 3526000 + t4 / 863310000;

//三个必要的参数:
a1 = 119.75 + 131.849 * t;
a2 = 53.09 + 479264.29 * t;
a3 = 313.45 + 481266.484 * t;
e = 1 - 0.002516 * t - 0.0000074 * t2;

cal_sum1 ();
sum1 = sum1+3958 * sin(a1) + 1962 * sin(ml - f) + 318 * sin(a2);
mw = ml + sum1 / 1000000-0.016;
mw=mw-(int)(mw/360)*360;
if (mw<0)mw=mw+360;
}//cal_moon () return mw

cal_sum1 (){ //calculate moon data
double La[60];
double Lb[60];
double Lc[60];
double Ld[60];
double sL[60];
double siu[60];
double sin1;
//double de;
string bas;
string bbs,ccs,dds;
string si[60];
// de=pi/180;
bas=“0,2,2,0,0,0,2,2,2,2,0,1,0,2,0,0,4,0,4,2,2,1,1,2,2,4,2,0,2,2,1,2,0,0,2,2,2,4,0,3,2,4,0,2,2,2,4,0,4,1,2,0,1,3,4,2,0,1,2,2,”;
bbs=“00000000010000-100-10100010000000000000001010001-10000000100-100-20102-20000-1000001-1020201-10000-1000100010000-102010000”;
ccs=“01-100020000-2-10100-10001000101-103-2-100-100010200-3-2-1-201000200-10100-102-101-2-1-1-200010400-2000201-2-30201-103-1”;
dds=“00000000000200000000000000-202-200000000000000000000000002000000000000-2020002000000000000-200000000-2-200000000000000-2”;

si[0]=“6288744”; //6288774
si[1]=“1274027”;
si[2]=“658314”;
si[3]=“213618”;
si[4]="-185116";
si[5]="-114332";
si[6]=“58793”;
si[7]=“57066”;
si[8]=“53322”;
si[9]=“45758”;
si[10]="-40923";
si[11]="-34720";
si[12]="-30383";
si[13]=“15327”;
si[14]="-12528";
si[15]=“10980”;
si[16]=“10675”;
si[17]=“10034”;
si[18]=“8548”;
si[19]="-7888";
si[20]="-6766";
si[21]="-5163";
si[22]=“4987”;
si[23]=“4036”;
si[24]=“3994”;
si[25]=“3861”;
si[26]=“3665”;
si[27]="-2689";
si[28]="-2602";
si[29]=“2390”;
si[30]="-2348";
si[31]=“2236”;
si[32]="-2120";
si[33]="-2069";
si[34]=“2048”;
si[35]="-1773";
si[36]="-1595";
si[37]=“1215”;
si[38]="-1110";
si[39]="-892";
si[40]="-810";
si[41]=“759”;
si[42]="-713";
si[43]="-700";
si[44]=“691”;
si[45]=“596”;
si[46]=“549”;
si[47]=“537”;
si[48]=“520”;
si[49]="-487";
si[50]="-399";
si[51]="-381";
si[52]=“351”;
si[53]="-340";
si[54]=“330”;
si[55]=“327”;
si[56]="-323";
si[57]=“299”;
si[58]=“294”;
si[59]=“0”;

double as,bs,cs,ds;
double sisum;
double susum;
sum1=0;
for (i=0;i<60;i++){
s=subString (bas,i2,1);
La[i]=stringToDouble (s);
as=as+La[i];
s=subString (bbs,i2,2);
Lb[i]=stringToDouble (s);
bs=bs+Lb[i];
s=subString (ccs,i2,2);
Lc[i]=stringToDouble (s);
cs=cs+Lc[i];
s=subString (dds,i2,2);
Ld[i]=stringToDouble (s);
ds=ds+Ld[i];
sin1=La[i]*D+Lb[i]*m+Lc[i]*m1+Ld[i]*f;
s=si[i]; //sumL60$
sL[i]=stringToDouble (s);
sisum=sisum+sL[i];
susum=susum+sin1;
sum1=sum1+sL[i]sin(sin1de);
}

//print "as = ",as; //查看 调试
//print "bs = ",bs;
//print "cs = ",cs;
//print "ds = ",ds;
//print "sumLs (0-59) = ",susum;
//print "si (60) sum= ",sisum;

}//cal_sum1

//***** end *****