fj2.htm

返回 第一章 关于寿星万年历 第二章 公历与儒略日 第三章 回历计算 第四章 天文坐标 第五章 近现代农历的天文学性质及其规则 第六章 离散序列的直线拟合算法

<script language=javascript>
function int2(v){ return Math.floor(v);}
var JD={ //日期元件
  Y:2000, M:1, D:1, h:12, m:0, s:0,
  toJD:function(){ //公历转儒略日
   var y=this.Y, m=this.M, n=0; //取出年月
   if(m<=2) m+=12,y--;
   if(this.Y*372+this.M*31+this.D>=588829)//判断是否为格里高利历日1582*372+10*31+15
     n =int2(y/100), n =2-n+int2(n/4);//加百年闰
   n +=int2(365.25*(y+4716)+0.01);    //加上年引起的偏移日数
   n +=int2(30.6*(m+1))+this.D;       //加上月引起的偏移日数及日偏移数
   n +=((this.s/60+this.m)/60+this.h)/24 - 1524.5;
   return n;
  },
  setFromJD:function(jd){ //儒略日数转公历
   jd+=0.5;
   var A=int2(jd), F=jd-A, D;  //取得日数的整数部份A及小数部分F
   if(A>=2299161) D=int2((A-1867216.25)/36524.25),A+=1+D-int2(D/4);
   A     +=1524; //向前移4年零2个月
   this.Y =int2((A-122.1)/365.25);//年
   D      =A-int2(365.25*this.Y); //去除整年日数后余下日数
   this.M =int2(D/30.6001);       //月数
   this.D =D-int2(this.M*30.6001);//去除整月日数后余下日数
   this.Y-=4716; this.M--;
   if(this.M>12) this.M-=12;
   if(this.M<=2) this.Y++;
   //日的小数转为时分秒
   F*=24; this.h=int2(F); F-=this.h;
   F*=60; this.m=int2(F); F-=this.m;
   F*=60; this.s=F;
  },
  toStr:function(){ //日期转为串
   var Y="     "+this.Y,M="0"+this.M, D="0"+this.D;
   var h=this.h,m=this.m,s=int2(this.s+.5);
   if(s>=60) s-=60,m++;
   if(m>=60) m-=60,h++;
   h="0"+h; m="0"+m; s="0"+s;
   Y=Y.substr(Y.length-5,5); M=M.substr(M.length-2,2); D=D.substr(D.length-2,2);
   h=h.substr(h.length-2,2); m=m.substr(m.length-2,2); s=s.substr(s.length-2,2);
   return Y+"-"+M+"-"+D+" "+h+":"+m+":"+s;
  }
};



var gan=Array("甲","乙","丙","丁","戊","己","庚","辛","壬","癸");
gan['甲']=0;gan['乙']=1;gan['丙']=2;gan['丁']=3;gan['戊']=4;
gan['己']=5;gan['庚']=6;gan['辛']=7;gan['壬']=8;gan['癸']=9;
function getsjb(lun,nd,ra,rb){//取数据表
 var i,j,k,p,c,v,jn;
 var n1=-1,n2=-1;
 nd=nd.split(" "); for(i=0;i<6;i++) nd[i]-=0; if(!nd[6]) nd[6]="";
 for(i=0;i<lun.length;i++){
  if(lun[i].substr(0,2)!="年份") continue;
  v=lun[i].substr(2,5)-0;
  if(v==nd[0]) n1=i;
  if(v==nd[3]){n2=i; break;}
 }
 if(n1==-1||n2==-1) return "年代中的年份参数有误";
 n1 += nd[1],  n2 += nd[4];
 ra.length=rb.length=0;
 //提取朔日和节气
 for(i=0,k=0,p=0;i<=n2-n1;i++){
   v=lun[i+n1]; if(v.substr(0,2)=="年份") continue;
   v=v.split(" ");
   ra[p++]=v[1]; //月大小
   j=3,jn=v.length; //节气行内始终
   if(i==0) j=nd[2]+3-1;
   if(i==n2-n1) { jn=nd[5]+3; if(jn>v.length) jn=v.length; }
   for(;j<jn;j++) rb[k++]=v[j];
   if(!i){
     //朔起始日期
     ra.lifa=nd[6]; //历法名称
     JD.Y=nd[0]-0; JD.M=v[2].substr(0,2)-0; JD.D=v[2].substr(2,2)-0; JD.h=12; JD.m=JD.s=0;
     if(nd[1]<5&&JD.M>8) JD.Y--;
     if(nd[1]>7&&JD.M<4) JD.Y++;
     ra.jd=JD.toJD();
     ra.rq0 = JD.toStr().substr(0,11);
     ra.n0 =Math.floor((ra.jd+10-2451545)/29.5306);

     //节气起始日期
     var x=rb[0].substr(0,2),xx=0;
     if(x=="一日") xx=1; if(x=="二日") xx=2; if(x=="三日") xx=3; if(x=="四日") xx=4; if(x=="五日") xx=5;
     if(x=="六日") xx=6; if(x=="七日") xx=7; if(x=="八日") xx=8; if(x=="九日") xx=9; if(x=="十日") xx=10;
     if(x=="十一") xx=11;if(x=="十二") xx=12;if(x=="十三") xx=13;if(x=="十四") xx=14;if(x=="十五") xx=15;
     if(x=="十六") xx=16;if(x=="十七") xx=17;if(x=="十八") xx=18;if(x=="十九") xx=19;if(x=="二十") xx=20;
     if(x=="廿一") xx=21;if(x=="廿二") xx=22;if(x=="廿三") xx=23;if(x=="廿四") xx=24;if(x=="廿五") xx=25;
     if(x=="廿六") xx=26;if(x=="廿七") xx=27;if(x=="廿八") xx=28;if(x=="廿九") xx=29;if(x=="三十") xx=30;
     rb.jd = ra.jd+xx-1;
     JD.setFromJD(rb.jd);
     rb.rq0 = JD.toStr().substr(0,11);
     rb.n0 = Math.floor((rb.jd+6-2451545+292)/365.2422*24);
   }
 }
 //历月处理
 for(i=ra.length-1;i>0;i--)
  if(ra[i-1]=="大") ra[i]=30; else ra[i]=29;
 ra[0]=0;
 for(i=1;i<ra.length;i++) ra[i]+=ra[i-1];
 //节气处理
 for(i=0;i<rb.length;i++) rb[i]=gan[rb[i].substr(2,1)];
 for(i=rb.length-1;i>0;i--){
   v=rb[i]-rb[i-1];
   if(v<0)  v+=10;
   if(v<10) v+=10;
   rb[i]=v;
 }
 rb[0]=0;
 for(i=1;i<rb.length;i++) rb[i]+=rb[i-1];
}


//平气朔处理
function nihe(r,k,b){//拟合计算,k、b为初始估值,r为序列数组
 var i,kk,v;
 var A=1, a=A-2e-4;  //调试时可把高度放宽，以便发现问题，取 1+1e-7，最后输出时高度减小一些，取1-2e-4。

 var k1=k-1, k2=k+1, f1, f2, x1;
 for(kk=0;kk<10;kk++){ //查找次数
  f1=-1e9, f2=1e9, x1=-1;
  for(i=0;i<r.length;i++){
   v=r[i]-(k*i+b);     //直线至少要平移的量
   if(v>f1) f1=v,x1=i;
   if(v<f2) f2=v;
  }
  if(f2+a-f1>0){ r.gh = (f2+A-f1)/2; b+=f1+r.gh; r.gk = k; r.gb = b; return "ok"; }
  for(i=0;i<r.length;i++){
   if(i==x1) continue;
   v=( r[i]+a - r[x1] ) / (i-x1);
   if(i<x1 && v>k1) k1=v;
   if(i>x1 && v<k2) k2=v;
   if(k1>k2) return "无解";
  }
  k=(k1+k2)/2;
 }
 return "查找"+kk+"次未找到解";
}


function calc(){
 var i,j,c="",C="",c2,c3,c4, se=0,qe=0;
 var ra=new Array(), rb=new Array();
 var nd=fen.innerText, lun=nlb.innerText;
 nd =nd.split("\r\n"); //年代分段表
 lun=lun.split("\r\n"); //农历表
 for(i=0;i<nd.length;i++){
  if(!nd[i]) continue;
  getsjb(lun,nd[i],ra,rb);
  if(rb.jd<1947158){// 619-01-11,平朔结束
   c2=nihe(ra,29.5306,0);
   if(c2=="ok") c += (ra.jd-0.5+ra.gb).toFixed(6)+","+ra.gk.toFixed(8)+", //"+ra.rq0+" h="+ra.gh.toFixed(5)+" "+ ra.lifa+"\r\n";     else c += "//"+ra.rq0+c2+"\r\n" ;
   if(c2=="ok") se = (ra.jd-0.5+ra.gb+ra.gk*ra.length).toFixed(2);
  }
  c3=nihe(rb,365.2422/24,0); if(c3=="ok") c4 = " 回归年="+(rb.gk*24).toFixed(5);
  if(c3=="ok") C += (rb.jd-0.5+rb.gb).toFixed(6)+","+rb.gk.toFixed(9)+", //"+rb.rq0+" h="+rb.gh.toFixed(5)+" "+ ra.lifa+c4+"\r\n";  else C += "//"+ra.rq0+c2+"\r\n";
  if(c3=="ok") qe = (rb.jd-0.5+rb.gb+rb.gk*rb.length).toFixed(2);
 }
 JD.setFromJD(se-0); c+=se+"//"+JD.toStr().substr(0,11)+"\r\n";
 JD.setFromJD(qe-0); C+=qe+"//"+JD.toStr().substr(0,11)+"\r\n";
 var c0="//气朔日期计算公式：D = k*n + b  , 式中n=0,1,2,3,...,N-1\r\n";
 c0+="//h表示k不变b允许的误差,如果b不变则k许可误差为h/N\r\n";
 c0+="//每行第1个参数为b,第2参数为b\r\n\r\n";
 out.innerText=c0+"//朔直线拟合参数\r\n"+c+"//气直线拟合参数\r\n"+C;
}



//定气朔处理
function suo_low(W){ //低精度定朔计算,在2000年至600，误差在2小时以内
  var v = 7771.37714500204;
  var t  = ( W + 1.08472 )/v, L;
  t -= ( -0.0000331*t*t
+ 0.10976 *Math.cos( 0.785 + 8328.6914*t)
+ 0.02224 *Math.cos( 0.187 + 7214.0629*t)
- 0.03342 *Math.cos( 4.669 +  628.3076*t ) )/v
+ (32*(t+1.8)*(t+1.8)-20)/86400/36525;
  return t*36525 + 2451545;
}
function qi_low(W){ //最大误差小于30分钟，平均5分
  var t,L,v= 628.3319653318;
  t =  ( W - 4.895062166 )/v; //第一次估算,误差2天以内
  t -= ( 53*t*t + 334116*Math.cos( 4.67+628.307585*t) + 2061*Math.cos( 2.678+628.3076*t)*t )/v/10000000; //第二次估算,误差2小时以内

  L = 48950621.66 + 6283319653.318*t + 53*t*t //平黄经
 +334166 * Math.cos( 4.669257+  628.307585*t)  //地球椭圆轨道级数展开
   +3489 * Math.cos( 4.6261  + 1256.61517*t )  //地球椭圆轨道级数展开
 +2060.6 * Math.cos( 2.67823 +  628.307585*t ) * t  //一次泊松项
   - 994 - 834*Math.sin(2.1824-33.75705*t); //光行差与章动修正

   t -= (L/10000000 -W )/628.332 + (32*(t+1.8)*(t+1.8)-20)/86400/36525;
   return t*36525 + 2451545;
}

function yasuo(s){
 s=s.replace(/0{60}/g,'A'); s=s.replace(/0{50}/g,'B');
 s=s.replace(/0{40}/g,'C'); s=s.replace(/0{30}/g,'D');
 s=s.replace(/0{20}/g,'E'); s=s.replace(/0{10}/g,'F');

 s=s.replace(/0{9}1/g,'a');
 s=s.replace(/0{8}1/g,'b');
 s=s.replace(/0{7}1/g,'c');
 s=s.replace(/0{6}1/g,'d');
 s=s.replace(/0{5}1/g,'e');
 s=s.replace(/00001/g,'f');
 s=s.replace(/0001/g,'g');
 s=s.replace(/001001/g,'h');
 s=s.replace(/001/g,'i');
 s=s.replace(/0101/g,'j');
 s=s.replace(/01/g,'k');

 s=s.replace(/0{9}2/g,'l');
 s=s.replace(/0{8}2/g,'m');
 s=s.replace(/0{7}2/g,'n');
 s=s.replace(/0{6}2/g,'o');
 s=s.replace(/0{5}2/g,'p');
 s=s.replace(/0{4}2/g,'q');
 s=s.replace(/0{3}2/g,'r');
 s=s.replace(/0{2}2/g,'s');
 s=s.replace(/02/g,'t');

 s=s.replace(/0{5}/g,'G');
 s=s.replace(/0{4}/g,'H');
 s=s.replace(/0{3}/g,'I');
 s=s.replace(/0{2}/g,'J');
 return s;
}
function fenghang(s,cn,na){ //分行
 var i,r="";
 for(i=0;i<s.length;i+=cn){
   if(!i) r += na+' ="';
   else   r += na+'+="';
   r += s.substr(i,cn)+'";'+"\r\n";
 }
 return r;
}

function suobj(r){
 var i,j,s="",S="",c="",C="",c2,d,d2,D;
 var ra=new Array(), rb=new Array();
 var nd="619 1 1 1960 12 1"; //年代
 var lun=nlb.innerText; lun=lun.split("\r\n"); //农历表
 getsjb(lun,nd,ra,rb);

 for(i=0,d=ra.jd; i<ra.length; i++){
   d = ra.jd+ra[i];
   d2 = suo_low((ra.n0+i)*Math.PI*2)+8/24;
   D = int2(d2+0.5);
   c2=d-D;
   if(c2>1||c2<-1) { alert(ra.n0);c+="err"; break; }
   if(c2<0) c2=2;
   c+=c2; if(c.length>200) C+=c,c="";
   if( D!=d ){
    JD.setFromJD(d2);
    s += JD.toStr() + " "+(d-D)+"\r\n";
    if(s.lengty>200) S+=s,s="";
   }
 }
 S+=s; C+=c; C=fenghang(yasuo(C),100,"suoS");
 out.innerText="//"+ra.rq0+"开始"+ra.length+"个朔日修正表 d0="+ra.jd+"\r\n" + C+"\r\n"+s;
}

function qibj(r){
 var i,j,s="",S="",c="",C="",c2,d,d2,D;
 var ra=new Array(), rb=new Array();
 var nd="1645 1 1 1960 12 1"; //年代
 var lun=nlb.innerText; lun=lun.split("\r\n"); //农历表
 getsjb(lun,nd,ra,rb);

 for(i=0,d=rb.jd; i<rb.length; i++){
   d = rb.jd+rb[i];
   d2 = qi_low((rb.n0+i)*Math.PI/12)+8/24;

   D = int2(d2+0.5);
   c2=d-D;
   if(c2>1||c2<-1) { alert(rb.n0);c+="err"; break; }
   if(c2<0) c2=2;
   c+=c2; if(c.length>200) C+=c,c="";
   if( D!=d ){
    JD.setFromJD(d2);
    s += JD.toStr() + " "+(d-D)+"\r\n";
    if(s.lengty>200) S+=s,s="";
   }
 }
 S+=s; C+=c; C=fenghang(yasuo(C),100,"qiS");
 out.innerText="//"+rb.rq0+"开始"+rb.length+"个节气修正表 d0="+rb.jd+"\r\n" + C + "\r\n" + s;
}


</script>

<a href=index.htm>返回</a>
<input type=button onclick="calc()" value=平气朔>
<input type=button onclick="qibj()" value=定气>
<input type=button onclick="suobj()" value=定朔><br>
<textarea id=fen rows=15 cols=60>
-103 1 1 84 12 2 汉书·律历志(太初历)平气平朔
85 1 1 236 12 2 后汉书·律历志(四分历)
237 1 1 444 12 2 晋书·律历志(景初历)
445 1 1 509 12 2 宋书·律历志(何承天元嘉历)
510 1 1 589 12 2 宋书·律历志(祖冲之大明历)
590 1 1 596 12 2 随书·律历志(开皇历)
597 1 1 618 12 2 随书·律历志(大业历)
619 1 1 665 12 2 新唐书·历志(戊寅元历)平气定朔
666 1 1 728 12 2 新唐书·历志(麟德历)
729 1 1 762 10 2 新唐书·历志(大衍历,至德历)
762 11 1 783 12 2 新唐书·历志(五纪历)
784 1 1 821 12 2 新唐书·历志(正元历,观象历)
822 1 1 892 12 2 新唐书·历志(宣明历)
893 1 1 955 12 2 新唐书·历志(崇玄历)
956 1 1 963 12 2 旧五代·历志(钦天历)
964 1 1 983 12 2 宋史·律历志(应天历)
983 1 1 1000 12 2 宋史·律历志(乾元历)
1001 1 1 1064 11 2 宋史·律历志(仪天历,崇天历)
1064 12 1 1067 12 2 宋史·律历志(明天历)
1068 1 1 1074 12 2 宋史·律历志(崇天历)
1075 1 1 1093 12 2 李锐补修(奉元历)
1094 1 1 1102 12 2 宋史·律历志
1103 1 1 1105 12 2 李锐补修(占天历)
1106 1 1 1135 12 2 宋史·律历志(纪元历)
1136 1 1 1190 12 2 宋史·律历志(统元历,乾道历,淳熙历)
1191 1 1 1198 12 2 宋史·律历志(会元历)
1199 1 1 1207 12 2 宋史·律历志(统天历)
1208 1 1 1251 3 2 宋史·律历志(开禧历)
1252 4 1 1252 12 2 淳祐历
1253 1 1 1270 12 2 会天历
1271 1 1 1276 12 2 宋史·律历志(成天历)
1277 1 1 1279 12 2 本天历
1280 1 1 1644 12 2 元史·历志(郭守敬授时历)
</textarea><br>
<textarea id=out rows=15 cols=120></textarea><br>





<textarea rows=15 cols=120>
对《三千五百年历日天象》校准如下

一、朔的错误
-26年正月和二月的大小有误
18年的十、十一月的大小有误
747 年十一月日期错误，应是12.7
1095年十二月日期错误
1732年八、十月的大小有误

二、以下是节气订正
-26年 二月十四 改 十五
-26年 二月廿九 改 三十
18年 十一月十三 改 十二
18年 十一月廿九 改 廿八
28年 九月十七 改 十八
44年 十二月十六 改 十七
81年 九月三日 改 四日
89年 五月十二 改 十三
91年 十二月十一 改 十二
134年 二月廿八 改 廿九
365年 十二月九日 改 十日
466年 二月十八 改 十九
616年 四月十四 改 十三
627年 四月二七 改 十七
876年 八月十四 改 十日
1115年 十月十三 改 十二
1149年 八月十三 改 十二
1427年 十月九日 改 十日
1719年 十一月十一 改 十二
1728年 十二月六日 癸改壬
1732年 九月五日 戊改己, 二十 癸改甲
1732年 十月五日 戊改己, 二十 癸改甲
1754年 二月廿八 壬改戊

698 十一正及十月分别第一个节气减2日和1日，变为壬日
724 十二月十三改十二