右栏

中华农历论坛 → 知识讨论区 → 历法知识 → 天文算法讨论


  共有227799人关注过本帖树形打印

主题：天文算法讨论
帅哥哟，离线，有人找我吗？
xjw01
 小大  31楼 个性首页 | 信息 | 搜索 | 邮箱 | 主页 | UC


加好友 发短信
等级：蝙蝠侠 帖子：866 积分：3804 威望：3 精华：20 注册：2008/3/20 22:14:00 		  发帖心情 Post By：2008/5/29 22:09:00

第8贴算法测试程序——地球、太阳

这是精简的超高速的算法

<script language=javascript>
rad = 180*3600/Math.PI;
function rad2mrad(v){   //对超过0-2PI的角度转为0-2PI
  v=v % (2*Math.PI);
  if(v<0) return v+2*Math.PI;
  return v;
}
function rad2str(d,tim){ //将弧度转为字串
 //tim=0输出格式示例: -23°59' 48.23"
 //tim=1输出格式示例:  18h 29m 44.52s
 var s="+";
 var w1="°",w2="’",w3="”";
 if(d<0)  d=-d,s='-';
 if(tim){ d*=12/Math.PI; w1="h ",w2="m ",w3="s "; }
 else     d*=180/Math.PI;
 var a=Math.floor(d); d=(d-a)*60;
 var b=Math.floor(d); d=(d-b)*60;
 var c=Math.floor(d); d=(d-c)*100;
 d=Math.floor(d+0.5);
 if(d>=100) d-=100, c++;
 if(c>=60)  c-=60,  b++;
 if(b>=60)  b-=60,  a++;
 a="   "+a, b="0"+b, c="0"+c, d="0"+d;
 s+=a.substr(a.length-3,3)+w1;
 s+=b.substr(b.length-2,2)+w2;
 s+=c.substr(c.length-2,2)+".";
 s+=d.substr(d.length-2,2)+w3;
 return s;
}

//以下是地球黄经数据,最大误差0.25"
var E_L0=new Array(
33416565,4.6692568,6283.07584999,
348943,4.626102,12566.1517,34971,2.74412,5753.38488,34176,2.82887,3.52312,31359,3.62767,77713.77147,
26762,4.41808,7860.41939,23427,6.13516,3930.2097,13243,0.74246,11506.76977,12732,2.0371,529.69097,
11992,1.10963,1577.34354,9903,5.2327,5884.9268,9019,2.0451,26.2983,8572,3.5085,398.149,
7798,1.1788,5223.6939,7531,2.5334,5507.5532,5053,4.5829,18849.2275,4924,4.2051,775.5226,
3567,2.9195,0.0673,3171,5.849,11790.6291,2841,1.8987,796.298,2710,0.3149,10977.0788,
2428,0.3448,5486.7778,2062,4.8065,2544.3144,2054,1.8695,5573.1428,2023,2.4577,6069.7768,
1555,0.8331,213.2991,1322,3.4112,2942.4634,1262,1.083,20.7754,1151,0.6454,0.9803,
1029,0.636,4694.003,1019,0.9757,15720.8388,1017,4.2668,7.1135,992,6.21,2146.165,
976,0.681,155.42,858,5.983,161000.686,851,1.299,6275.962,847,3.671,71430.696,
796,1.808,17260.155,788,3.037,12036.461,747,1.755,5088.629,739,3.503,3154.687,
735,4.679,801.821,696,0.833,9437.763,624,3.978,8827.39,611,1.818,7084.897,
570,2.784,6286.599,561,4.387,14143.495,556,3.47,6279.553,520,0.189,12139.554,
516,1.333,1748.016,511,0.283,5856.478,490,0.487,1194.447,410,5.368,8429.241,
409,2.399,19651.048,392,6.168,10447.388,368,6.041,10213.286,366,2.57,1059.382,
360,1.709,2352.866,356,1.776,6812.767,333,0.593,17789.846,304,0.443,83996.847,
300,2.74,1349.867,254,3.165,4690.48,247,0.215,3.59,237,0.485,8031.092,
236,2.065,3340.612,228,5.222,4705.732,219,5.556,553.569,214,1.426,16730.464,
211,4.148,951.718,203,0.371,283.859,199,5.222,12168.003,199,5.775,6309.374,
191,3.822,23581.258,189,5.386,149854.4,179,2.215,13367.973,175,4.561,135.065,
162,5.988,11769.854,151,4.196,6256.778,144,4.193,242.729,143,3.724,38.028,
140,4.401,6681.225,136,1.889,7632.943,125,1.131,5.523,121,2.622,955.6,
120,1.004,632.784,113,0.177,4164.312,108,0.327,103.093,105,0.939,11926.254,
105,5.359,1592.596,103,6.2,6438.496,100,6.029,5746.271,98,1,11371.7,
98,5.24,27511.47,94,2.62,5760.5,92,0.48,522.58,92,4.57,4292.33,
90,5.34,6386.17,86,4.17,7058.6,84,3.3,7234.79,84,4.54,25132.3,
81,6.11,4732.03,81,6.27,426.6,80,5.82,28.45,79,1,5643.18,
78,2.96,23013.54,77,3.12,7238.68,76,3.97,11499.66,73,4.39,316.39,
73,0.61,11513.88,72,4,74.78,71,0.32,263.08,68,5.91,90955.55,
66,3.66,17298.18,65,5.79,18073.7,63,4.72,6836.65,62,1.46,233141.31,
61,1.07,19804.83,60,3.32,6283.01,60,2.88,6283.14,55,2.45,12352.85);

var E_L1=new Array(
2060589,2.6782346,6283.07585,43034,2.63513,12566.1517,4253,1.5905,3.5231,1193,5.7956,26.2983,
1090,2.9662,1577.3435,935,2.592,18849.228,721,1.138,529.691,678,1.875,398.149,
673,4.409,5507.553,590,2.888,5223.694,560,2.175,155.42,454,0.398,796.298,
364,0.466,775.523,290,2.647,7.114,208,5.341,0.98,191,1.846,5486.778,
185,4.969,213.299,173,2.991,6275.962,162,0.032,2544.314,158,1.43,2146.165,
146,1.205,10977.079,125,2.834,1748.016,119,3.258,5088.629,118,5.274,1194.447,
115,2.075,4694.003,106,0.766,553.569,100,1.303,6286.599,97,4.24,1349.87,
95,2.7,242.73,86,5.64,951.72,76,5.3,2352.87,64,2.65,9437.76,
61,4.67,4690.48,58,1.77,1059.38,53,0.91,3154.69,52,5.66,71430.7,
52,1.85,801.82,50,1.42,6438.5,43,0.24,6812.77,43,0.77,10447.39,
41,5.24,7084.9,37,2,8031.09,36,2.43,14143.5,35,4.8,6279.55,
34,0.89,12036.46,34,3.86,1592.6,33,3.4,7632.94,32,0.62,8429.24,
32,3.19,4705.73,30,6.07,4292.33,30,1.43,5746.27,29,2.32,20.36,
27,0.93,5760.5,27,4.8,7234.79,25,6.22,6836.65,23,5,17789.85,
23,5.67,11499.66,21,5.2,11513.88,21,3.96,10213.29,21,2.27,522.58,
21,2.22,5856.48,21,2.55,25132.3,20,0.91,6256.78,19,0.53,3340.61,
19,4.74,83996.85,18,1.47,4164.31,18,3.02,5.52,18,3.03,5753.38,
16,4.64,3.29,16,6.12,5216.58,16,3.08,6681.22,15,4.2,13367.97,
14,1.19,3894.18,14,3.09,135.07,14,4.25,426.6,13,5.77,6040.35,
13,3.09,5643.18,13,2.09,6290.19,13,3.08,11926.25,12,3.45,536.8);

var E_L2=new Array(
87198,1.0721,6283.07585,3091,0.8673,12566.1517,273,0.053,3.523,163,5.188,26.298,
158,3.685,155.42,95,0.76,18849.23,89,2.06,77713.77,70,0.83,775.52,
51,4.66,1577.34,41,1.03,7.11,38,3.44,5573.14,35,5.14,796.3,
32,6.05,5507.55,30,1.19,242.73,29,6.12,529.69,27,0.31,398.15,
25,2.28,553.57,24,4.38,5223.69,21,3.75,0.98,17,0.9,951.72,
15,5.76,1349.87,14,4.36,1748.02,13,3.72,1194.45,13,2.95,6438.5,
12,2.97,2146.17,11,1.27,161000.69,10,0.6,3154.69,10,5.99,6286.6,
9,4.8,5088.63,9,5.23,7084.9,8,3.31,213.3,8,3.42,5486.78,
7,6.19,4690.48,7,3.43,4694,6,1.6,2544.31,6,1.98,801.82,
6,2.48,10977.08,5,1.44,6836.65,5,2.34,1592.6,5,1.31,4292.33,
5,3.81,149854.4,4,0.04,7234.79,4,4.94,7632.94,4,1.57,71430.7,
4,3.17,6309.37,3,0.99,6040.35,3,0.67,1059.38,3,3.18,2352.87,
3,3.55,8031.09,3,1.92,10447.39,3,2.52,6127.66,3,4.42,9437.76,
3,2.71,3894.18,3,0.67,25132.3,3,5.27,6812.77,3,0.55,6279.55);

var E_L3=new Array(
2892,5.8438,6283.0758,168,5.488,12566.152,30,5.2,155.42,
13,4.72,3.52,7,5.3,18849.23,6,5.97,242.73,4,3.79,553.57,1,4.3,6286.6,1,0.91,6127.66);

var E_L4=new Array(77,4.13,6283.08,8,3.84,12566.15,4,0.42,155.42);
var E_L5=new Array(2,2.77,6283.08,1,2.01,155.42);

//地球黄纬数据,误差0.2"
var E_B0=new Array(
2796,3.1987,84334.6616,1016,5.4225,5507.5532,804,3.88,5223.694,438,3.704,2352.866,
319,4,1577.344,227,3.985,1047.747);

var E_B1=new Array(90,3.9,5507.55,62,1.73,5223.69);
var E_B2=new Array(17,1.63,84334.66);

//地球向径数据,误差0.00001AU
var E_R0=new Array(
1000139888,0,0,16706996,3.09846351,6283.07584999,
139560,3.055246,12566.1517,30837,5.19847,77713.77147,
16285,1.17388,5753.38488,15756,2.84685,7860.41939,
9248,5.4529,11506.7698,5424,4.5641,3930.2097,
4721,3.661,5884.9268,3460,0.9637,5507.5532,
3288,5.8998,5223.6939,3068,0.2987,5573.1428,
2432,4.2735,11790.6291,2118,5.8471,1577.3435,
1858,5.0219,10977.0788,1748,3.0119,18849.2275);

var E_R1=new Array(1030186,1.1074897,6283.07585,17212,1.06442,12566.1517,7022,3.1416,0);
var E_R2=new Array(43594,5.78455,6283.07585,1236,5.5793,12566.1517,123,3.142,0,88,3.63,77713.77);
var E_R3=new Array(1446,4.2732,6283.0758,67,3.92,12566.15);
var E_R4=new Array(39,2.56,6283.08,3,2.27,12566.15);
var E_R5=new Array(1,1.22,6283.08);

 

function Enn(F,t,n){ //计算E_L0或E_L1或E_L2等
  n=Math.floor(n*F.length+0.5); //按百分比取项数
  var i,v=0;
  for(i=0;i<n;i+=3) v+=F[i]*Math.cos(F[i+1] +t*F[i+2]);
  return v;
}
function earthLon(t,n){ //地球经度计算,返回Date分点黄经,传入千年数
  if(n<0) n=1; else n = 3*n/E_L0.length;
  var t2=t*t, t3=t2*t, t4=t3*t, t5=t4*t;
  var v  = 1753469512 + 6283319653318*t + 529674*t2 + 432*t3 - 1157*t4 - 9*t5; //地球平黄经(已拟合DE406)
      v += 630 * Math.cos(6+3*t); //拟合DE406
      v += Enn( E_L0, t, n );
      v += Enn( E_L1, t, n )*t;
      v += Enn( E_L2, t, n )*t2;
      v += Enn( E_L3, t, n )*t3;
      v += Enn( E_L4, t, n )*t4;
      v += Enn( E_L5, t, n )*t5;
  return v/1000000000;
}
function earthCoor(t,re,n){ //返回地球坐标
  var t2=t*t, t3=t2*t, t4=t3*t, t5=t4*t;
  re[0] = earthLon(t,n);
  re[1] = Enn( E_B0, t, 1 )
        + Enn( E_B1, t, 1 )*t
        + Enn( E_B2, t, 1 )*t2;
  re[2] = Enn( E_R0, t, 1 )
        + Enn( E_R1, t, 1 )*t
        + Enn( E_R2, t, 1 )*t2
        + Enn( E_R3, t, 1 )*t3
        + Enn( E_R4, t, 1 )*t4
        + Enn( E_R5, t, 1 )*t5;
  re[1]/=1000000000;
  re[2]/=1000000000;
}
function earthV(t){ //地球速度,t是千年数,误差小于万分3
 var f=6283.07585*t;
 return v = 6283.32
  +210 *Math.sin(1.527+f)
  +4.4 *Math.sin(1.48+f*2)
  +12.9*Math.sin(5.82+f)*t
  +0.55*Math.sin(4.21+f)*t*t;
}
function earth_t(W){ //已知黄经求时间
  var t,dW,v= 6283.319653318;
  t  = ( W - 1.75347       )/v; v=earthV(t,1);   //v的精度0.03%，详见原文
  t += ( W - earthLon(t,10))/v;
  t += ( W - earthLon(t,-1))/v;
  return t;
}

//迭代算法测试
t=1;
L=earthLon(t,-1);
t2=earth_t(L);
document.write("迭代算法测试<br>");
document.write("输入时间(千年):"+t+"<br>");
document.write("地球黄经(弧度):"+L+"<br>");
document.write("返算时间(千年):"+t2+"<br>");
document.write("迭代误差(秒):"+(t2-t)*365250*86400+"<br><br>");

 

 

//=======================
// 以下是视位置计算
//=======================
var nutB=new Array(
  2.1824,  -33.75705, 36e-6,-1720,920,
  3.5069, 1256.66393, 11e-6,-132, 57,
  1.3375,16799.4182, -51e-6, -23, 10,
  4.3649,  -67.5141,  72e-6,  21, -9,
  0.04,   -628.302,   0,     -14,  0,
  2.36,   8328.691,   0,       7,  0,
  3.46,   1884.966,   0,      -5,  2,
  5.44,  16833.175,   0,      -4,  2,
  3.69,  25128.110,   0,      -3,  0,
  3.55,    628.362,   0,       2,  0);
function nutation(t){ //章动计算
 var i,c,a, t2=t*t,dL=0,dE=0;
 for(i=0;i<nutB.length;i+=5){
  c = nutB[i]+nutB[i+1]*t+nutB[i+2]*t2;
  if(i==0) a=-1.742*t; else a=0;
  dL+=(nutB[i+3]+a)*Math.sin(c);
  dE+= nutB[i+4]   *Math.cos(c);
 }
 this.dL=dL/100;  //黄经章动
 this.dE=dE/100;  //交角章动
}
function gxc_sun(t,L){ //太阳光行差(黄经),L为太阳黄经,t是千年数
 var v = L-(282.93735+17.1946*t+ 0.046*t*t)/180*Math.PI; //真近点角
 var e = 0.016708634-0.00042037*t-0.00001267*t*t;
 return -20.49552 * (1+e*Math.cos(v));
}
function HCconv(JW,E){ //黄赤转换(黄赤坐标旋转)
  //黄道赤道坐标变换,赤到黄E取负
  var HJ=rad2mrad(JW[0]), HW=JW[1];
  var sinE =Math.sin(E),cosE =Math.cos(E);
  var sinW=cosE*Math.sin(HW)+sinE*Math.cos(HW)*Math.sin(HJ);
  var J=Math.atan2( Math.sin(HJ)*cosE-Math.tan(HW)*sinE, Math.cos(HJ) );
  JW[0]=rad2mrad(J);
  JW[1]=Math.asin(sinW);
}

function hcjj (t){ //返回黄赤交角,t是千年数
 var t2=t*t, t3=t2*t,t4=t3*t;
 return (84381.4088 -468.36051*t -0.01667*t2 -1.99911*t3-0.00523*t4)/rad;
}

 

z=new Array();
da=1;
testT=0.008+(da-1.5)/365250;    //世纪数
earthCoor(testT,z,-1);         //地球黄经
z[0]+=Math.PI,z[1]=-z[1];      //太阳黄经黄纬

//L=z[0]-50287.9*(testT-0.008)/rad;
//L=rad2mrad(L);
//document.write(rad2str(L,0)+"<br>");

gx= gxc_sun(testT,z[0]);       //太阳光行差
nu= new nutation(testT*10);    //章动计算
z[0] += (nu.dL+gx)/rad;        //补章动与光行差
E = hcjj(testT)+nu.dE/rad;     //得到真黄赤交角
HCconv(z,E); //转到赤道坐标

document.write("坐标精度测试<br>");
document.write("testT：" + testT+"<br>");
document.write("视赤经：" + rad2str(z[0],1)+"<br>");
document.write("视赤纬：" + rad2str(z[1],0)+"<br>");
document.write("距离：" + z[2]+"<br>");
document.write("请与《2008中国天文年历》太阳视赤经赤纬比对。该年的第 "+da+" 天0h TD。注意：几乎分秒不差！<br>");
document.write("本程序中VSOP87地球数据序列以做了长期项拟合到DE405/406<br>");

</script>

[此贴子已经被作者于2008-5-30 8:05:57编辑过]
支持(0) 中立(0) 反对(0)  单帖管理 | 引用 | 回复 回到顶部
帅哥哟，离线，有人找我吗？
xjw01
 小大  32楼 个性首页 | 信息 | 搜索 | 邮箱 | 主页 | UC


加好友 发短信
等级：蝙蝠侠 帖子：866 积分：3804 威望：3 精华：20 注册：2008/3/20 22:14:00 		  发帖心情 Post By：2008/6/1 21:17:00

完整日月测试程序

适用于公元前3000年到公元5000年,超出此范围误差增加很多

 下载信息  [文件大小：   下载次数： ]
图片点击可在新窗口打开查看点击浏览该文件:
支持(0) 中立(0) 反对(0)  单帖管理 | 引用 | 回复 回到顶部
帅哥哟，离线，有人找我吗？
xjw01
 小大  33楼 个性首页 | 信息 | 搜索 | 邮箱 | 主页 | UC


加好友 发短信
等级：蝙蝠侠 帖子：866 积分：3804 威望：3 精华：20 注册：2008/3/20 22:14:00 		  发帖心情 Post By：2008/6/1 21:19:00

定气测试结果

视黄经0: 2004-03-20 14:48:39
视黄经15: 2004-04-04 18:43:20
视黄经30: 2004-04-20 01:50:25
视黄经45: 2004-05-05 12:02:29
视黄经60: 2004-05-21 00:59:13
视黄经75: 2004-06-05 16:13:47
视黄经90: 2004-06-21 08:56:52
视黄经105: 2004-07-07 02:31:17
视黄经120: 2004-07-22 19:50:10
视黄经135: 2004-08-07 12:19:36
视黄经150: 2004-08-23 02:53:15
视黄经165: 2004-09-07 15:12:54
视黄经180: 2004-09-23 00:29:50
视黄经195: 2004-10-08 06:49:17
视黄经210: 2004-10-23 09:48:48
视黄经225: 2004-11-07 09:58:32
视黄经240: 2004-11-22 07:21:40
视黄经255: 2004-12-07 02:48:57
视黄经270: 2004-12-21 20:41:37
视黄经285: 2005-01-05 14:03:00
视黄经300: 2005-01-20 07:21:35
视黄经315: 2005-02-04 01:43:02
视黄经330: 2005-02-18 21:31:58
视黄经345: 2005-03-05 19:45:10
视黄经360: 2005-03-20 20:33:26
视黄经375: 2005-04-05 00:34:16
视黄经390: 2005-04-20 07:37:14
视黄经405: 2005-05-05 17:52:50
视黄经420: 2005-05-21 06:47:24
视黄经435: 2005-06-05 22:01:53
视黄经450: 2005-06-21 14:46:07
视黄经465: 2005-07-07 08:16:34
视黄经480: 2005-07-23 01:40:42
视黄经495: 2005-08-07 18:03:22
视黄经510: 2005-08-23 08:45:27
视黄经525: 2005-09-07 20:56:39
视黄经540: 2005-09-23 06:23:11
视黄经555: 2005-10-08 12:33:17
视黄经570: 2005-10-23 15:42:20
视黄经585: 2005-11-07 15:42:25
视黄经600: 2005-11-22 13:14:56
视黄经615: 2005-12-07 08:32:40
视黄经630: 2005-12-22 02:34:55
视黄经645: 2006-01-05 19:46:57
视黄经660: 2006-01-20 13:15:16
视黄经675: 2006-02-04 07:27:14
视黄经690: 2006-02-19 03:25:33
视黄经705: 2006-03-06 01:28:39
视黄经720: 2006-03-21 02:25:34
视黄经735: 2006-04-05 06:15:30
视黄经750: 2006-04-20 13:26:03
视黄经765: 2006-05-05 23:30:40
视黄经780: 2006-05-21 12:31:33
视黄经795: 2006-06-06 03:36:59
视黄经810: 2006-06-21 20:25:52
视黄经825: 2006-07-07 13:51:28
视黄经840: 2006-07-23 07:17:43
视黄经855: 2006-08-07 23:40:48
视黄经870: 2006-08-23 14:22:34
视黄经885: 2006-09-08 02:39:01
视黄经900: 2006-09-23 12:03:22
视黄经915: 2006-10-08 18:21:23
视黄经930: 2006-10-23 21:26:28
视黄经945: 2006-11-07 21:34:51
视黄经960: 2006-11-22 19:01:45
视黄经975: 2006-12-07 14:26:48
视黄经990: 2006-12-22 08:22:04
视黄经1005: 2007-01-06 01:40:08
视黄经1020: 2007-01-20 19:00:48
视黄经1035: 2007-02-04 13:18:11
视黄经1050: 2007-02-19 09:08:56
视黄经1065: 2007-03-06 07:17:59
视黄经1080: 2007-03-21 08:07:26
视黄经1095: 2007-04-05 12:04:39
视黄经1110: 2007-04-20 19:07:04
视黄经1125: 2007-05-06 05:20:23
视黄经1140: 2007-05-21 18:11:56
视黄经1155: 2007-06-06 09:27:04
视黄经1170: 2007-06-22 02:06:25
视黄经1185: 2007-07-07 19:41:42
视黄经1200: 2007-07-23 13:00:10
视黄经1215: 2007-08-08 05:31:14
视黄经1230: 2007-08-23 20:07:57
视黄经1245: 2007-09-08 08:29:29
视黄经1260: 2007-09-23 17:51:14
视黄经1275: 2007-10-09 00:11:30
视黄经1290: 2007-10-24 03:15:24
视黄经1305: 2007-11-08 03:24:00
视黄经1320: 2007-11-23 00:49:52
视黄经1335: 2007-12-07 20:14:04
视黄经1350: 2007-12-22 14:07:47
视黄经1365: 2008-01-06 07:24:49
视黄经1380: 2008-01-21 00:43:30
视黄经1395: 2008-02-04 19:00:22
视黄经1410: 2008-02-19 14:49:32
视黄经1425: 2008-03-05 12:58:48
视黄经1440: 2008-03-20 13:48:17
视黄经1455: 2008-04-04 17:45:51
视黄经1470: 2008-04-20 00:51:08
视黄经1485: 2008-05-05 11:03:26
支持(0) 中立(0) 反对(0)  单帖管理 | 引用 | 回复 回到顶部
帅哥哟，离线，有人找我吗？
xjw01
 小大  34楼 个性首页 | 信息 | 搜索 | 邮箱 | 主页 | UC


加好友 发短信
等级：蝙蝠侠 帖子：866 积分：3804 威望：3 精华：20 注册：2008/3/20 22:14:00 		  发帖心情 Post By：2008/6/1 21:24:00

定朔测试结果

月-日黄经差0
2004-12-12 09:29:02
2005-01-10 20:02:49
2005-02-09 06:28:02
2005-03-10 17:10:21
2005-04-09 04:31:59
2005-05-08 16:45:25
2005-06-07 05:55:06
2005-07-06 20:02:30
2005-08-05 11:04:46
2005-09-04 02:45:25
2005-10-03 18:27:54
2005-11-02 09:24:38
2005-12-01 23:00:56
2005-12-31 11:11:44
2006-01-29 22:14:37
2006-02-28 08:30:45
2006-03-29 18:15:16
2006-04-28 03:43:53
2006-05-27 13:25:35
2006-06-26 00:05:16
2006-07-25 12:30:57
2006-08-24 03:09:47
2006-09-22 19:45:02
2006-10-22 13:14:03
2006-11-21 06:18:01
2006-12-20 22:00:47
2007-01-19 12:00:41
2007-02-18 00:14:17
2007-03-19 10:42:34
2007-04-17 19:36:01
2007-05-17 03:27:13
2007-06-15 11:13:07
2007-07-14 20:03:46
2007-08-13 07:02:31
2007-09-11 20:44:14
2007-10-11 13:00:38
2007-11-10 07:03:04
2007-12-10 01:40:25
2008-01-08 19:37:10
2008-02-07 11:44:30
2008-03-08 01:14:14
2008-04-06 11:55:17
2008-05-05 20:18:12
2008-06-04 03:22:33
2008-07-03 10:18:34
2008-08-01 18:12:31
2008-08-31 03:58:01
2008-09-29 16:12:14
2008-10-29 07:13:50
2008-11-28 00:54:34
2008-12-27 20:22:24
2009-01-26 15:55:18
2009-02-25 09:35:03
2009-03-27 00:05:54
2009-04-25 11:22:33
2009-05-24 20:10:59
2009-06-23 03:34:56
2009-07-22 10:34:35
2009-08-20 18:01:34
2009-09-19 02:44:16
2009-10-18 13:33:03
2009-11-17 03:13:42
2009-12-16 20:02:04
2010-01-15 15:11:20
2010-02-14 10:51:17
2010-03-16 05:01:05
2010-04-14 20:28:53
2010-05-14 09:04:21
2010-06-12 19:14:33
2010-07-12 03:40:26
2010-08-10 11:08:08
2010-09-08 18:29:48
2010-10-08 02:44:26
2010-11-06 12:51:43
2010-12-06 01:35:39
2011-01-04 17:02:33
2011-02-03 10:30:36
2011-03-05 04:45:50
2011-04-03 22:32:18
2011-05-03 14:50:40
2011-06-02 05:02:34
2011-07-01 16:53:53
2011-07-31 02:39:47
2011-08-29 11:04:04
2011-09-27 19:08:40
2011-10-27 03:55:46
2011-11-25 14:09:39
2011-12-25 02:06:21
2012-01-23 15:39:14
2012-02-22 06:34:33
2012-03-22 22:37:05
2012-04-21 15:18:24
2012-05-21 07:47:00
2012-06-19 23:02:03
2012-07-19 12:23:59
2012-08-17 23:54:26
2012-09-16 10:10:39
2012-10-15 20:02:30
2012-11-14 06:08:00
2012-12-13 16:41:37
支持(0) 中立(0) 反对(0)  单帖管理 | 引用 | 回复 回到顶部
帅哥哟，离线，有人找我吗？
xjw01
 小大  35楼 个性首页 | 信息 | 搜索 | 邮箱 | 主页 | UC


加好友 发短信
等级：蝙蝠侠 帖子：866 积分：3804 威望：3 精华：20 注册：2008/3/20 22:14:00 		  发帖心情 Post By：2008/8/30 19:33:00

《天文算法》第2章 关于精度

  [许剑伟 于莆田十中 2008年4月27日]

  本章所要讨论的主题是：具体问题所需要的精度；程序设计语言的精度；最后结果发表的精度。

具体给定问题所需的精度

  所需的计算精度，取决于具体问题的目标要求。例如，如果想计算行星的升或降的时刻，计算到0.01度已足够。原因是明显的，天球周日运动，每转过1度对应4分钟，所以0.01度的位置误差，引起的升或降时刻的误差大约为0.04分钟。为了解决这个问题，如果我们考虑数百个周期项，把位置精度控制到了0".01，那是在浪费精力和计算机的时间。

  但是，如果要计算行星引起的恒星食，那么就要求行星位置的计算精度高于1"，因为行星圆面的尺寸很小。

  为某一具体问题编写的程序，不一定适用于另于相似的问题。假设，为了计算恒星的位置，一个程序使用了第20章的“低精方法”得到的岁差。观测者利用望远镜寻找星体，使用这种精度的位置是足够的。但是，在计算“星食”或“会合”时，象这样低精度的位置是没有价值的。

  如果所需的精度已确定，我们必须使用一种适当的算法，满足其精度要求。John Mosley 提到了一个计算行星位置的商用价值的程序，但由于没有考虑摄动，即使在最近的几个角秒，土星、天王星、海王星的位置误差达1度。

  为了取得更高的精度，不是简单地对大约值保留更多的有效数字，常常需要使用其它方法。例如，如果想获得精确到0.01度的火星位置，使用无摄的椭圆轨道(开普勒运动)就已足够。然而，如果对火星位置精度要求达10"或更高，就必须考虑其它行星的摄动，程序将变得比较长。

  对于程序员，必须清楚他的公式及具体问题的精度要求，充分考虑各种可省略的项让程序变得更简洁。例如，涉及Date平分点的太阳几何平黄经：

  L = 280°27'59".244+129602771".380T+1".10915T2

  式中T是从历元2000年1月1.5日起算的儒略世纪数(36525个历书日)。在这一表达式中，如果|T|<0.95，即在1905年到2095年，最后一项(太阳的世纪加速度)小于1"。如果1"精度已足够，那么，在这一时期内，T2项可以忽略。但时，在公元100年，T=-19，最后一项的值是394"，大于0.1度。

计算机的精度

  这是一个很复杂的问题。程序设计语言可工作在足够位数的有效数字。注意：与小数位是不同的！例如，0.0000183有7位小数，但只有3个有效数字。有效数字指，除去左边前导的所有的"0"后，剩下数字的位数。

  一台机器，四舍五入到6位有效数字，结果1000000+2就是1000000。

  计算两个十分接近数的差值时，是很危险的。假设执行了以下减法：6.92736 - 6.92735 = 0.00001.每个数字的有效数字是6位，但是，相减后的结果却只有一位有效数字。另外，两个给定的数可能已做了四舍五入，那么情况将更糟糕。假设两个数的精确值是6.9273649和6.9273451，那么正确的结果是0.0000198，它几乎是刚才的两倍。

  6位或8位有效数字（早期的微机或如今的“单精度”），均不能满足天文计算的要求。

  在很多应用中，要求计算机的有效数字位数比最后结果的有效数字大很多位。例如，让我们考虑任意时刻的月亮平黄经公式，单位是度(参见第45章)：

L' = 218.3164591 + 481267.88134236T - 0.0013268T2 + 0.0000019T3

  式中T是历元2000年1月1.5TD起算的儒略世纪数(36525日)。假设，我们要得到0.001度的月亮平黄经精度。因为经度限制在0——360度，你可能会想到，计算机内部精度只需6位就已足够(小数点之前3位，小数点之后3位)。然而，这个问题不是这样的，因为，在转到0——360度之前，L'的值可能很大。

  例如，让我们计算2040年1月1日12h TD，即T=0.4时，L'的值。我们得到L'=192725.469度，转换后变为125.469度。但是，如果电脑工作在6位有效数字，得到的不再是L'=192725.469度，而192725度(6位！)，转换后变为125度，所以最后的结果精度只有1度，误差0.569度或28'，而且这仅发生在历元后的40年。在这种情况下，不可能用于计算日月食或星食。

  以下BASIC小程序可以找出电脑的内部精度：

10 X=1
20 J=0
30 X=X*2
40 IF X+1<>X THEN 60
50 GOTO 80
60 J=J+1
70 GOTO 30
80 PRINT J,J*0.30103
90 END

  这里，J是浮点数的尾数部分的比特数，而0.30103J是10进制数的有效数字。常数0.30103是log10(2)。例如，HP-85计算机的J=39，因此是11.7位。在HP-UX Technical Basic 5.0，工作在HP-Integral微机，我们得到J=52，因此是15.6位。QUICK-BASIC 4.5中得到的是J=63，因此是19.0位。

  不过，这个精度指简单的算术运算的精度，不是指三角函数。虽然日本东芝公司带GWBASIC电脑，J=55，有16.6位有效数字，但计算正弦时，只有7位是正确的，至少9位完全错误。

  一种迅速检测三角函数计算精度的方法是PRINT 4*ATN(1)。如果计算机工作在弧度制，应得到一个著名的数PI=3.14159265358979...或者，可以计算一个某个角的正弦值，该正弦值精确已知的，例如：

  sin(0.61弧度) = 0.57286746010048...

  计算机的舍入误差不是可避免的。例如，考虑这个值，1/3=0.33333333...因为计算机不能处理无限小数，一个数字必须在某处截断。

  从某次计算到下一次计算，舍入误差可能被积累。多数情况下，这是不要紧的，因为这种误差几乎可以忽略。但是，在某些算术应用中，这种误差积累可能无限制的增加，虽然这已超出本书的范围，不过，我们还要提到两种情况。

  考虑以下程序：

10 X=1/3
20 FOR I=1 TO 30
30 X=(9*X+1)*X-1
40 PRINT I,X
50 NEXT I
60 END

  第30行，实际上是将x赋值给它本身(1/3赋值给x)。可是，大部分计算机将得到无穷大的结果。用刚才提到的HP-UX Technical Basic得到：

4 步之后 0.333 333 333 333 308
14步之后 0.333 326 162 117 054
19步之后 0.215 899 338 763 055
24步之后 286.423...
30步之后,其值为10217

  不同的计算机的计算结果是不相同的，你甚至要以用手持式计算器做个测试：对1.0000001连续平方，27次之后，10位有效数字的结果应是674530.4707。以下是一些电脑或程序语言得到的结果：

674 494.06 在HP-67计算器
674 514.87 在HP-85微机
674 520.61 在TI-58计算器
674 530.4755 在HP-Integral(HP-UX Techn. Basic)
674 530.4755 在QUICKBASIC 4.5

  但是，故事还没结束。计算机用的数值信息有两种不同的基本表示方法。有些计算机，如老式的HP-85，使用BCD(二进制编码的10进制数)方案表示内部数值，而在其它的多数电脑中，使用二进制表示一个数。

  BCD编码方案是：对十进制数的每一位单独编码。这样，对于给定的计算机或程序设计言，所有的数字都可以被准确的表示了。另一方面，二进制编码，是用2的各次方的组合表示一个数。在二进制中，分数也是用2的各次方表示的，所以，不能用2的负次方的组合的的分数不是能用二进制表示的。例如1/10是不能用2的负次方组合得到，因为1/10 = 1/16+1/32+1/128...

  二进制算术运算比相应的BCD运算快得多，但是，有个不便之处就是，有些数(甚至是只有很少小数位的数)不能精确表示。

  其后果是，算术运算的结果可能不正确，哪怕是只有一点点的误差。假设x=4.34，H=INT(100*(x-INT(x)))的正确结果是34，然而，很多计算机语言得到的结果是33。原因是，此时计算机内部表示x的值为4.3399999998，或其它类似的值。

  另一个吃惊的例子：2+0.2+0.2+0.2+0.2+0.2-3

  在很多计算机中，结果不是0！在HP-Integral中，使用HP-UX Technical Basic 5.0，得到的结果是8.88*10-16。而在同样的电脑中，计算0.2+0.2+0.2+0.2+0.2+2-3的结果却是0，可见计算的顺序也是很重要的。让人惊呀的是，在HP-Integral中，使用以下程序计算2+(5*0.2)-3,我们们得到的结果零:

A=0.2+0.2+0.2+0.2+0.2
B=2+A
C=B-3
PRINT C

  考虑以下程序：

10 FOR I=0 TO 100 STEP 0.1
20 U=I
30 NEXT I
40 PRINT U
50 END

  这里，I和U取值依次是0到100，步长为0.1，最后得到的U应是100。在HP-85中的确是100，但在HP-Integral却得到99.9999999999986，这在某些应用中可能造成灾难性的后果。这个误差是由于0.1转为二进制表示，相当于0.0999999...，它与0.1相差很小，但是这里执行了1000步，最后的误差将放大了1000倍。在这种情况下，补救方法是，使用整数步长：

10 FOR J=0 TO 1000
20 I=J/10
30 U=I
40 NEXT J
50 PRINT U
60 END

  我们发现另一个吃惊的结果：

A=3*(1/3)
PRINT INT(A)

  在某些计算机中可能得到正确的结果1，有些却得到0。还可以试试这个例子：A = 0.1，PRINT INT(1000*A0)。

  另一个有趣的测试：

INPUT A
B=A/10
C=10*B
PRINT A-C

  结果应当是0，但不同的A值，结果可能有点差别。

  有一种简单的方法，判断计算机语言是否工作在BCD：通过查找整型变量所表示的最大整数值来判断。如果这个最大整数值是个“漂亮的整数”，那么表明机器工作在BCD编码。例如，在HP-85的最大整数是99999(或10^5-1)。但是，如果这个最大整数值是个陌生值(实际上它的2的n方，再减1)，那么机器不是工作在BCD。在早期的TRS-80，最大整数是32767(或2^15-1)，HP-Integral电脑上的HP-UX Technical Basic 5.0其最大整数是2147483647(或2^31-1)。

  算术运算中的四舍五入，可能导至其它吃惊的结果。在多数微机中，sqrt(25)-5的结果不是0！如果用这个结果作为关键测试值将是有问题的。根号25是个整数值吗？答案是否定的，因为计算机告诉我们sqrt(25)-INT(sqr(25))不等零。

最后结果的舍入

  结果应做舍入，保留正确的有意义的位数。

  “舍入”是取“最近”的近似值。例如，15.88取值为15.9或16，不是15。然而，日期或年的除外。例如，3月15.88日，指“一个时刻”附加在3月15日上：指3月15日0h之后的0.88日。所以，一个事件发生在3月15.88日，是指发生在3月15日，而不是3月16日。类似的，1977.69指1977年的某一时刻，而不是1978年。

  只有“有意义”的数字应当保留。例如，计算木星的可视星等的Muller公式：

m = -8.93 + 5 log(rΔ)

  式中r是木星到太阳的距离，Δ它到地球的距离，二者的单位是天文单位，对数计算是以10为底的。在1992年5月14日0h TD，我们有：
    r = 5.417149,Δ= 5.125382
  因此，m=-1.712514898。但是，如果以电脑算出这样的结果为借口，保留了所有的数字，这是十分可笑的，而且读者可能误解为精度非常的高。因为Muller公式中的常数-8.93的精度在0.01量级，你不能期望结果的精度比它更高。另外，由于木星大气的气象变化也将造成行星的星等的精度不超过0.01(甚至是0.1)。

  另一个例子，John Mosley 提到的用于计算天体升或降时刻的商用程序，精度是0.1秒钟，这也是不可能达到的精度。

  一些“感觉”及相关天文知识是很需要的。例如，计算月亮圆面被照亮部分的比例时，完全没有必要精确到0.000000001。

  执行舍入处理，应在所有的计算完成之后，而不是在程序开始或数据输入时。

  例如，计算1.4+1.4，结果取整。如果我们一开始就执行舍入，我们得到 1+1 = 2。实际上 1.4+1.4 = 2.8，四舍五入后得3。

  这里有另外一个例子。海王星“冲”的日期1996年7月18日，赤纬δ=-20°24'，那么行该星经过南子午圈时，地平纬度hm是多少？设观测站在德国Someberg天文台，纬度是φ=+50°23'，结果精确到度。使用以下公式：

  hm = 90°- φ + δ

  答案是: hm = 90°- 50°23' - 20°24' = 19°13' 因此是19°

  如果在计算之前对φ和δ做舍入计算，将得到不正确的结果：90°-50°-20°=20°.

  再看一个类似的错误：先把距离舍入到英里，再转换为千米。例如，在这种情况下，17km永远没有了,因为：
  10英里 = 16.09km，舍入后得16km
  11英里 = 17.70km，舍入后得18km

  赤经和赤纬。——因为24小时对应360度，1小时对应15度，时间1分钟对应角度15分(注意，是15m而不是15')，1秒钟对应角度1秒(即1s而不是1")，在1秒钟的时间内，地球自转了15"。

  正是由于这个原因，如果给出的天体赤纬精度为1"，赤经的精度则应舍入到十分之一时秒，否则赤纬的精度将高于赤经的精度。下表是赤经(α)和赤纬(δ)精度的大约对应关系。例如，如果δ的精度是1'，那么α应给出0.1m的精度。做为示例，表中还给出了Nova Cygni 1975的不同精度表示。

αδ例(Nova Cygni 1975)
1m
0m.1
1s
0s.1		0.1°
1'
0.1'
1"		α=21h10m
   21h09m.9
   21h09m53s
   21h09m52s.89		δ=+47°.9
   +47°57'
   +47°56'.7
   +47°56'41"

  最后注意，让我们讲述一下尾部的“0”，这很重要。例如，18.0与18是不同的。前者意味着真值介于17.95到18.05之间，而后者被舍入为整数，真值介于17.5到18.5之间。由于这个原因，为了表示精度，尾部的0必须给出：恒星的星等为7，它与7.00是不同的。

参考资料

1、John Mosley，《天空和望远镜》，卷78，第300页(1989年9月)
2、F.Gruenberger，《科学美国人》，“电脑娱乐”，卷250，第10页(1984年4月)
3、John Mosley，《天空和望远镜》，卷81，第201页(1991年2月)
支持(0) 中立(0) 反对(0)  单帖管理 | 引用 | 回复 回到顶部
帅哥哟，离线，有人找我吗？
xjw01
 小大  36楼 个性首页 | 信息 | 搜索 | 邮箱 | 主页 | UC


加好友 发短信
等级：蝙蝠侠 帖子：866 积分：3804 威望：3 精华：20 注册：2008/3/20 22:14:00 		  发帖心情 Post By：2008/8/30 19:37:00

《天文算法》第3章插值

 下载信息  [文件大小：   下载次数： ]
图片点击可在新窗口打开查看点击浏览该文件:
支持(0) 中立(0) 反对(0)  单帖管理 | 引用 | 回复 回到顶部
帅哥哟，离线，有人找我吗？
xjw01
 小大  37楼 个性首页 | 信息 | 搜索 | 邮箱 | 主页 | UC


加好友 发短信
等级：蝙蝠侠 帖子：866 积分：3804 威望：3 精华：20 注册：2008/3/20 22:14:00 		  发帖心情 Post By：2008/8/30 19:44:00

第4章 曲线拟合

 下载信息  [文件大小：   下载次数： ]
图片点击可在新窗口打开查看点击浏览该文件:

 下载信息  [文件大小：   下载次数： ]
图片点击可在新窗口打开查看点击浏览该文件:
支持(0) 中立(0) 反对(0)  单帖管理 | 引用 | 回复 回到顶部
帅哥哟，离线，有人找我吗？
kmlch
 小大  38楼 个性首页 | 信息 | 搜索 | 邮箱 | 主页 | UC


加好友 发短信
等级：论坛游民 帖子：58 积分：3673 威望：0 精华：0 注册：2008/8/6 9:17:00 		  发帖心情 Post By：2008/9/6 16:43:00

先收藏，好贴！

支持(0) 中立(0) 反对(0)  单帖管理 | 引用 | 回复 回到顶部
帅哥哟，离线，有人找我吗？
hanxiao129
 小大  39楼 个性首页 | 信息 | 搜索 | 邮箱 | 主页 | UC


加好友 发短信
等级：论坛游民 帖子：59 积分：583 威望：0 精华：0 注册：2007/12/1 11:11:00 		  发帖心情 Post By：2008/9/10 9:13:00

学习一下，不错的

支持(0) 中立(0) 反对(0)  单帖管理 | 引用 | 回复 回到顶部
帅哥哟，离线，有人找我吗？
peigong
 小大  40楼 个性首页 | 信息 | 搜索 | 邮箱 | 主页 | UC


加好友 发短信
等级：新手上路 帖子：2 积分：192 威望：0 精华：0 注册：2008/10/17 15:55:00 		  发帖心情 Post By：2008/10/17 16:00:00

准备做一个计算机易学方面的网站。

许老师的帖子，不能不注册了过来顶一下。

有没有人把许老师的这组帖子整理成一个文档啊，没有的话我来干一下，有的话发我一份，呵呵。
支持(0) 中立(0) 反对(0)  单帖管理 | 引用 | 回复 回到顶部
总数 115 上一页 1 2 3 4 5 6 7 8 9 10 下一页 ..12

返回版面帖子列表

天文算法讨论
签名：
RSS2.0　　Xhtml无图版　　Xslt无图版 Copyright ©2005 - 2018 NongLi.Com Email:chinanongli@tom.com   京ICP备05002438号 Powered By Dvbbs8.2.0