DE405/406星历表算法

#pragma hdrstop
#pragma argsused //有入口参数的此行不能少
#include <stdio.h>
#include <stdlib.h>
#include <conio.h>
#include "sxwnl_eph.cpp"
/*===============================================================
                  [ DE星历表计算程序 ]
           实际测试过DE405及DE406,  xjw01@莆田十中 2009.2
·本程序在C++Builder6.0中编译通过
  如果VC6.0编译,请加上 #include "stdafx.h" 即可。但在在VC中，
  本程序的读写速度降低数倍。
·程序的功能：
  1.把DE星历表的文本星历库转为二进制格式
  2.进行星历计算测试并给出实现算法
  ·本程序根据JPL网站提供的DE星历表CD盘资料严格计算星历,与几万个标准
  数据比较，误差小于10^-13
·JPL网站上提供的DE星历表分为两部分：
  header.xxx头文件以及插值数据库，可以在以下地址下载
  ftp://ssd.jpl.nasa.gov/pub/eph/planets/ascii/
  该ftp中还提供了插值计算程序,有c、fortran、java等版,其中c程序无
  法在window下编译通过。fortran等版与其提供的头文件也不匹配,也无
  法调试成功。他们所提供的exe测试文件也无法在window平台下运行,所
  以编写些程序解决这些问题.
·插值系数数据库(以下简称数据库)由多个txt文件组成,每个文件覆盖一
  定的时间范围.如"ascm1900.406"数据库适用于BC1900到BC1800年,
  "ascp1900.406"适用范围是1900到2000年。不同版本的DE星历表,各数据
  库适用的范围不一定相同,还须根据数据库内部的数据进行分析。
·同一版本的数据库可以使用copy命令连起来,合并为一个数据库,如：
  copy ascm0100.406+ascp0000.406+ascp0100.406+ascp0200 data.406
  但应注意,应根据时间的前后关系按顺序连接。如果使把ascm0100.406与
  ascp0100.406直接连接起来,中间就出现了100年的断裂,将造本程序将识
  别成功。
·每个数据库中由若干个数据块构成,各数据块按时间顺序连接。不过使用
  copy命令连接后,在两个文的连接处有一个重复数据块,读取时应跳过连接
  处的重复块。各数块适用的时间范围一般只有几十天，如DE405每个数据
  块均为32天。同一版本的星历表，所有数据块的时间长度是相同的。一个
  数据块中含有太阳系各个天体坐标的插值系数。不同版本包含的天体个数
  不一定全部相同。
·各天体插值系数在数据块中的索引位置由头文件定义。文件结构详见"JPL
  星历表文件结构示意图"
·插值算法采用切比雪夫多项式插值
·本程序生成的二进制文件与JPL官方定义的不太相同。其数据结构是原先
  fortran的，用c语言生成相同格式的文件十分烦麻，所以本程序放弃官方
  定义的格式标准而自行定义数据结构，这使程序变得很短，实现相同的功
  能，代码量不到500行，比官方推荐的数千行c代码精简5倍以上。另一方
  面，官方程序也表明，其c程序不一定能够在window下运行(实际上根本没
  有调试过)，所以自行定义数据结构是一种较好的选择。
================================================================*/
void ascii2bin(char* headName,char* dataName,char* outName,double startjd,double stopjd){ //头文件名,数据文件名,输出文件,起止JD
if (stopjd<startjd){ printf("错误:终JD小于始JD"); exit(0); }
FILE *fp1, *fp2, *fp3;
if((fp3=fopen(outName,"r")) != NULL) {
   fclose(fp3);
   printf("目标文件 %s 已存在,是否覆盖?y/n\n", outName);
   if (getch()!='y') exit(0);
}
if((fp1=fopen(headName,"r"))==NULL) { printf("错误: 无法打开header文件 %s.\n", headName); exit(0); }
if((fp2=fopen(dataName,"r"))==NULL) { printf("错误: 无法打开 data 文件 %s.\n", dataName); exit(0); }
if((fp3=fopen(outName,"wb"))==NULL) { printf("错误: 无法创建输出文件： %s.\n", outName ); exit(0); }
//读取ascii文件头
DE_Header H;
DE_rAscHeader(fp1, H);
//读取ascii数据块.
int i, n = 0; //n写入的块数
double* block = new double[H.nn]; //系数数据块数组分配
double JDp; //前一块的终时间
for(i=0;1;i++){
   if(!DE_rAscBlock(fp2, H, block)) break; //读数据块
   if(i%40==0) printf("扫描第%d块,已写入%d块.\r",i,n);
   if(block[0] > stopjd) break; //星历始JD已超出输入JD
   if(block[1] < startjd) continue;
   if(block[1] - block[0] != H.Ta) { //每块天数检查(header.405为32天)
     printf("错误: 头文指定每块 %g 天, 与数据块1的 (%g to %g)不匹配.程序终止.\n", H.Ta, block[0], block[1]);
     exit(1);
   }
   if (n && JDp != block[0]) { //上下块之间时间不连续
     if(JDp-H.Ta==block[0]) continue; //JPL的ascii星历文件相连接后,存在一个重复块,应忽略
     printf("错误:数据文件中相邻块(%d与%d块)的时间不连续(前JD%f 后JD%f),程序终止.\n", n,n+1, JDp,block[0]);
     exit(1);
   }
   DE_wBinBlock(fp3, H.nn, n, block);
   if(!n) H.JD1 = block[0]; //保存第一块始时间
   H.JD2 = JDp = block[1];        //保存最后一块终时间
   n++;
}
DE_wBinHeader(fp3, &H); //写数据头
printf("写入%s文件%d字节. 其中文件头2块,数据%d块(每块%d个系数).\n\n", outName,(2+n)*H.nn*8, n, H.nn);
fclose(fp1); fclose(fp2); fclose(fp3);
delete[] block;
}

void testEph(char*testpo, char* bin){
int i; char buf[DE_maxL+1];
FILE* fp;  DE_coord cd;
if( !cd.open(bin) ) { printf("打开 %s 失败。",bin); exit(0); }
if((fp=fopen(testpo,"r")) == NULL) { printf("%s 无法打开.",testpo); exit(0); }
// testpo中的行变量定义
char Adate[30]; // 日期
double Ajd;     // 测试时间(JD)
double Acd;     // 参考坐标
int Ade,At,Acenter,Ax; //版本号,天体号,坐标中心,坐标号
//临时变量
double r[6];      //天体的 x, xdot, y, ydot, z, zdot
double wd = 1e-13;//容许误差值,如果程序计算的结果与testpo相差大于wd,则视为有错
int    wn = 0;    //误差较大数据个数的计数
//测试开始
printf("\nde  ----date--- ---jed--- t# c# x ---testpo坐标---- ---本程序计算---- --出差--\n");
for(i=0;i<8;i++) fgets(buf,DE_maxL,fp); //读取前面的几个无效行
for(i=0;i<300000;i++){
  if(!fgets(buf,DE_maxL,fp)) break;
  if(i%300==0) printf("正在测试第%d个\r",i+1);
  sscanf(buf, "%d %s %lf %d %d %d %lf", &Ade, Adate, &Ajd, &At, &Acenter, &Ax, &Acd);
  if(cd.calc(Ajd)) {
    cd.getCoordOne(At-1,Acenter-1,r);
    double dd=Acd-r[Ax-1];//误差
    if(fabs(dd)>wd){ //打印出错数据
      printf("%3d %11s %9.1f %2d %2d %1d %17.13f %17.13f %8.2e\n", Ade, Adate, Ajd, At, Acenter, Ax, Acd, r[Ax-1], dd);
      wn++;
    }
  }else printf("JD%9.1f不在给定的星历表范围内\n",Ajd);
}
printf("\n测试报告：共测试 %d 个坐标, 误差大于 %7.1e 的数据有 %d 个.\n",i-1,wd,wn);
fclose(fp);
}
void debugPro(){ //asc转bin以及星历计算调试函数
printf("\n===== DE星历表调试程序 最后修改：2009.2 xjw01. =====\n");
printf("\n调试一：header.405 + ascp2000.405 转换为二进制文件 debug.405 ...\n\n");
ascii2bin("header.405","ascp2000.405","debug.405",2400000,2500000); //ascii转bin测试
//getch(); return;
printf("\n调试二：读取二进制debug.405并计算一个坐标 ...\n\n");
DE_coord cd;
double testjd=2453164.5, r1[6],r2[6]; //r1,r2两个天天体的坐标(x,y,z,vx,vy,vz)
if( !cd.open("debug.405") ) { printf("打开debug.405失败。"); exit(0); }
if( !cd.checkJD(testjd)   ) { printf("错误: JD%15.2f不在星历表范围内.\n", testjd); exit(0); }
//开始计算
if(cd.calc(testjd)) {
   cd.getCoordOne(DE_sun,   DE_earth, r1); //取太阳地心坐标
   cd.getCoordOne(DE_venus, DE_earth, r2); //取金星地心坐标
   printf("  时间 2004.06.08 00:00:00，JD=%15.7f\n\n", testjd);
   printf("         2004天文年历(AU)     计算结果(AU)          速度(AU/天)\n");
   printf("  太阳x     0.2199082    %20.17f  %20.17f\n", r1[0],r1[3]);
   printf("  太阳y     0.9091691    %20.17f  %20.17f\n", r1[1],r1[4]);
   printf("  太阳z     0.3941587    %20.17f  %20.17f\n", r1[2],r1[5]);
   printf("  太阳d     1.0150415    %20.17f\n", sqrt(r1[0]*r1[0]+r1[1]*r1[1]+r1[2]*r1[2]) );
   printf("  金星d     0.2888953    %20.17f\n", sqrt(r2[0]*r2[0]+r2[1]*r2[1]+r2[2]*r2[2]) );
   printf("  d = sqrt(x*x + y*y + z*z)\n\n");
}
}
void main(int argc, char *argv[]) {
//调试
//ascii2bin("header.406","data.406","bin.406",-9999999,9999999); return;
//debugPro(); getch();
//testEph("testpo.406","bin.406"); getch();
//return;
if(argc<2){
   printf("  ========================================================\n");
   printf("  ----------- 《寿星万年历》JPL星历表工具1.0  ------------\n");
   printf("  --------------------------------------------------------\n\n");
   printf("  1、命令范例--文本转二进制：\n\n");
   printf("    myJPL a header.405 data.405 bin.405 startJD stopJD\n");
   printf("    myJPL a header.406 data.406 bin.406 startJD stopJD\n");
   printf("    myJPL a header.406 data.406 bin.406\n\n");
   printf("    不指定JD范围则全文件转换\n\n");
   printf("  2、命令范例--星历计算测试(与标准testpo.xxx比对)：\n\n");
   printf("    myJPL b testpo.405 bin.405\n");
   printf("    myJPL b testpo.406 bin.406\n\n");
   printf("  3、命令范例--算法调试：jpl c\n\n");
   printf("    算法调试时,请先确认header.405与ascp2000.405存在\n\n");
   printf("  ===== 许剑伟 [e-mail] 莆田十中 2008.2 =====\n");
   return;
}
if(argv[1][0]=='c') debugPro();
if(argv[1][0]=='b'){ //星历计算测式
   if(argc!=4){ printf("错误：参数个数不正确\n"); return; }
   testEph(argv[2],argv[3]);
}
if(argv[1][0]=='a'){ //ascii转bin
   if(argc!=5&&argc!=7) {printf("错误：参数个数不正确\n"); return; }
   double startJD,stopJD;
   if(argc==5) startJD = -99999999,      stopJD = 99999999;
   if(argc==7) startJD = atof(argv[5]), stopJD = atof(argv[6]);
   ascii2bin(argv[2],argv[3],argv[4],startJD,stopJD);
}
}


// sxwnl_eph.cpp 文件

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
const DE_maxL = 128;  //ASCII文本每行最字符数设置
struct DE_Header { //文件头结构
int nn;           // 每个数据块含有切比雪夫(Chebyshev)系数的个数,它是ascii头文件中的NCOEFF
char ver[211];    // 星历表版本信息串
int nCon;         // 常数个数
char cCon[300][7];// 常名称
double vCon[300]; // 常数值
double au;        // 天文单位大小 km/1AU
double emrat;     // 地月质量比
double clight;    // 光速, km/s
int DEver;        // DE星历表版本号
int LEver;        // LE星历表版本号
double JD1,JD2;   // 始历元,终历元
double Ta;        // 每数据块对应的时间长度
int p1[13];       // 切比雪夫系数数组在文件中的位置索引辅助量(快速定位用的),详见结构图
int p2[13];       // 作用同上
int p3[13];       // 作用同上
};
//把形如 33D+3 的科学记数转为标准的 33e+3
void d2e(char *p,int n=1000){
int i;
for(i=0;i<n;i++){
  if(!p[i]) break;
  if(p[i]=='D'||p[i]=='d') p[i]='e';
}
}
//文件的读写
int  DE_rBinHeader(FILE *fp, DE_Header *h); // 读二进制文件头
int  DE_wBinHeader(FILE *fp, DE_Header *h); // 写二进制文件头
void DE_rAscHeader(FILE *fp, DE_Header &h); // 读ascii头文件
int  DE_rBinBlock (FILE *fp, int nn, int n, double *r); // 读二进制数据块
int  DE_wBinBlock (FILE *fp, int nn, int n, double* r); // 写二进制数据块
int  DE_rAscBlock (FILE *fp, DE_Header &h, double *r);  // 按顺序读ascii数据块
/*=============================
JPL文本星历文件读取函数
===============================*/
//在header文件(ascii)中跳过一个组号标识
void gotoNextGroup(char *group, FILE *fp) {
char g[13]; //"GROUP"信息缓存
char buf[DE_maxL + 1];
fgets(buf, DE_maxL, fp); //读空行
fgets(buf, DE_maxL, fp); //读 "GROUP   组号\n"
strncpy(g, buf, 12); g[12] = '\0';
if(strncmp(g,group,12) != 0) { printf("头文件有误;%s未找到.\n", group); exit(1); }
fgets(buf, DE_maxL, fp); //读空行
}
//读JPL的ASSCII星历文件头
void DE_rAscHeader(FILE *fp, DE_Header &h) {
int i, j, k, n;
char buf[DE_maxL + 1];
for(i=0; i<211; i++) h.ver[i]=0;
for(i=0; i<300; i++){
  for(j=0; j<7; j++) h.cCon[i][j]=0;
  h.vCon[i]=0;
}
rewind(fp);
//header文件第1行: KSIZE= 数字1 NCOEFF= 数字2
fgets(buf, DE_maxL, fp);
sscanf(buf, "%*6s%6d%*11s%6d", &i, &h.nn); //格式串中的*是丢弃输入的作用,ncoeff为切比雪夫系数的个数
if (i != 2*h.nn) {
   fprintf(stderr, "header文件格式错误：KSIZE <> 2*NCOEFF\n");
   exit(1);
}
//GROUP 1010: 星历标题, DE/LE版本号, 始JD,终JD
gotoNextGroup("GROUP   1010", fp);
fgets(buf, DE_maxL, fp); strcpy(h.ver,buf); // JPL Ephemeris 标题行
fgets(buf, DE_maxL, fp); strcat(h.ver,buf); // 始时间,该版本适的适用JD范围
fgets(buf, DE_maxL, fp); strcat(h.ver,buf); // 终时间
if (strncmp(h.ver, "JPL", 3)) { printf("错误: 该文件不是JPL ASCII头文件.\n"); exit(1); }
//GROUP 1030: 始JD, 终JD, 每块时间范围.
gotoNextGroup("GROUP   1030", fp);
fgets(buf, DE_maxL, fp);
sscanf(buf, " %le %le %le", &h.JD1, &h.JD2, &h.Ta);
//GROUP 1040: 常数名称
  gotoNextGroup("GROUP   1040", fp);
  fgets(buf, DE_maxL, fp);
  h.nCon = atoi(buf); //读到header.405的156
  // 读常数名称, 每行10个, 每个6字节
  for(i=0; i<h.nCon; i++){
   if (i%10==0) fgets(buf, DE_maxL, fp);
   sscanf(buf+(i%10)*8,"%s",h.cCon+i);
  }
//GROUP 1041: 常数值.
gotoNextGroup("GROUP   1041", fp);
fgets(buf, DE_maxL, fp);
if (atoi(buf) != h.nCon) { printf("错误: 常数个数与常数值的个数不相等.\n"); exit(1); }
for (i=0; i<h.nCon; i+=3) { // 读常数值, 每行3个, 每个26字节.
   fgets(buf, DE_maxL, fp);
   d2e(buf);
   sscanf(buf, "%le %le %le", h.vCon+i, h.vCon+i+1, h.vCon+i+2);
}
//GROUP 1050: 常数值.
gotoNextGroup("GROUP   1050", fp);
for(i=0;i<3;i++){
  fgets(buf, DE_maxL, fp);
  for (j=0; j<13; j++){
   if(i==0) h.p1[j] = atoi(buf + 6*j); //读p1
   if(i==1) h.p3[j] = atoi(buf + 6*j); //读p2
   if(i==2) h.p2[j] = atoi(buf + 6*j); //读p3
  }
}
h.au = h.emrat = 0;
h.DEver = h.LEver = 0;
for (i = 0; i < h.nCon; i++) {
   if      (strcmp(h.cCon[i], "AU")    == 0) h.au    = h.vCon[i]; //1AU距离的千米数
   else if (strcmp(h.cCon[i], "EMRAT") == 0) h.emrat = h.vCon[i]; //地月地量比
   else if (strcmp(h.cCon[i], "DENUM") == 0) h.DEver = h.vCon[i]; //DE星历表版本号(如200,403,405,406)
   else if (strcmp(h.cCon[i], "CLIGHT")== 0) h.clight= h.vCon[i]; //光速
   else if (strcmp(h.cCon[i], "LENUM") == 0) h.LEver = h.vCon[i]; //LE星历表版本号
}
if (!h.LEver) h.LEver = h.DEver;
}
//-----------------------------------------------------
//从ASCII星历文件中读一个系数数据块
int DE_rAscBlock( FILE *fp, DE_Header &h, double *r) {//按顺序读,数据块返在r中
int i,k, n=0;       //n为本次读入系数的个数
char buf[DE_maxL+1];//读缓存
double v[3];        //用于读取每行的3个double
//ASCII数据库中,第1行是块号及系数个数
if (!fgets(buf, DE_maxL, fp)) return n;
sscanf(buf, "%*d %d", &k);
if (k != h.nn) {
   printf("错误: ascii数据文件每块%d个系数,而头文件定义每块%d个, 二者不相同\n", k,h.nn);
   exit(1);
}
for(;;){
   if(!fgets(buf, DE_maxL, fp)) return n;
   d2e(buf); //d2e把readbuf转为标准的科学记数
   sscanf( buf, "%le %le %le", v, v+1, v+2);
   for(i=0; i<3; i++){
     r[n++] = v[i];
     if(n>=k) return n;
   }
}
}
/*==============================================
JPL二进制星历文件写入
==============================================*/
//写JPL星历二进制文件头
int DE_wBinHeader(FILE *outfp, DE_Header *h) {
rewind(outfp); //文件指针复位
int i, hsize=2*8*h->nn, hsize2=sizeof(DE_Header); //文件头大小
fwrite(h,hsize2,1,outfp);
for (i = hsize2; i < hsize; i++)  fputc('\0', outfp); //补成两块
return 0;
}
//-----------------------------------------------------
//写二进制系数数据块
int DE_wBinBlock(FILE *fp, int nn,int n, double* r) {
fseek(fp, 0, SEEK_END);
int i, fn=ftell(fp); //fn当前文件大小
int pos = (n + 2) * nn * 8; //计算写入点
for (i=0; i < pos-fn; i++) fputc('\0', fp); //如果文件大小不够,则补空字节
fseek(fp, pos, 0);//文件指针定位在写入点
fwrite(r, nn*8, 1, fp);
return 1;
}
/*==============================================
JPL二进制星历文件读取
==============================================*/
//读JPL二进制星历文件头
int DE_rBinHeader(FILE *fp, DE_Header *h) {
rewind(fp); //文件指针复位
fread(h,sizeof(DE_Header),1,fp);
if(strncmp(h->ver,"JPL",3)){ printf("当前读取的不是JPL文件二进制文件."); exit(0); }
return 0;
}
//-----------------------------------------------------
//读JPL星历二进制数据块
int DE_rBinBlock(FILE *fp, int nn, int n, double* r){//每块数据个数nn,块号n,返回r
fseek(fp, (n+2)*nn*8, 0);
if( !fread(r, nn*8, 1, fp) ) return -1;
return 1;
}
/*==============================================
JPL坐标计算
==============================================*/
//星体索引号
const DE_mercury = 0;
const DE_venus   = 1;
const DE_earth   = 2;
const DE_mars    = 3;
const DE_jupiter = 4;
const DE_saturn  = 5;
const DE_uranus  = 6;
const DE_neptune = 7;
const DE_pluto   = 8;
const DE_moon    = 9;  //月亮的太阳系质心坐标编号
const DE_sun     = 10;
const DE_SSBARY  = 11; //太阳系质心
const DE_EMBARY  = 12; //地月质心
const DE_nutation= 13;
const DE_LIBRATION=14;
const DE_GEOMOON = 15; //月亮地心坐标编号(插值计算的结果)
const DE_numBody = 16; //以上天体坐标个数
//-----------------------------------------------------
class DE_coord{
private: FILE *fp; //星历文件指针
public:
DE_Header h; double* block;  //星历文件头及数据块
double pv[DE_numBody][6];    //calc()返回值
//------初始化(文件操作、内存管理等)-----
DE_coord(){ fp=NULL, block=NULL; } //构造函数
int open(char* fname){ //打开星历表并读取文件头
   close();
   if( (fp=fopen(fname,"rb")) == NULL ) return 0;
   DE_rBinHeader(fp,&h); //读取文件头
   block = new double[h.nn];
   return 1;
}
void close(){ //关闭星历文件
   if(fp)    fclose(fp);
   if(block) delete[] block;
   fp=NULL, block=NULL;
}
~DE_coord(){ close(); } //析构
//--------------星历计算-----------------
int checkJD(double jd){ // 检查jd是否在星历表范围之内
   if(jd<h.JD1 || jd>h.JD2) return 0;
   return 1;
}
void getCoordOne(int a, int bCenter, double* r) { //在pv中提取一个星体的坐标
   int i,j;
   double *pa=pv[a], *pb=pv[bCenter];
   for (i=0; i<6; i++) {
    if(a==13||a==14) r[i] = pa[i]; // 提取章动(13)和天平动(14)
    else r[i] = pa[i]-pb[i];
   }
   if(a==14) r[4]=r[5]=0;
}
void cheby(double x, double* y, int nC, int n, double* r); //切比雪夫插值
int calc(double jd, int center=DE_SSBARY); // 坐标计算
};
//--------------------------------------------------------------------------
void DE_coord::cheby( double x,double* y, // Chebyshev 插值的x值,介于[-1,1], y值即Chebyshev多项式系数序列
   int nC,     // 坐标个数
   int n,      // 每坐标系数组的个数
   double *r){ // 返回在r数组中, 前一半是位置,后一半是速度
int i,j,k;
static double s[20],v[20]; //位置和速度的Tn(x)多项式,最大cheby系数不超过20个(de406中最多只有14个)
static double x0=2;        //保存最后一次计算时x的值,初值定为一个不可能使用的值
if (x0 != x) { //对于一个x值,其Tn(x)多项式值是固定的,即同一个JD,Tn(x)只需计算一次
   x0 = x;
   //初始化位置Tn多项式的值
   s[0] = 1.0; // T0(x) = 1
   s[1] = x;   // T1(x) = x
   for (i=2; i<20; i++)  s[i] = 2*x*s[i-1] - s[i-2];  //Tn(x)=2*x*T[n-1]-T[n-2],参见实用数学手册第2版
   //初始化Tn'多项式的值(Tn对x求导数),用于速计算
   v[0] = 0.0;  // T0' = 0
   v[1] = 1.0;  // T1' = 1
   for (i=2; i<20; i++) v[i] = 2*x*v[i-1] + 2*s[i-1] - v[i-2];
}
for (i=0; i<nC; i++) { // 遍历x,y,z
   k = i+nC, r[i]=r[k] = 0;;
   for(j=0; j<n; j++){
     r[i] += s[j] * y[i*n+j]; //位置累和计算(得位置插值结果)
     r[k] += v[j] * y[i*n+j]; //速度累和计算(得速度插值结果)
   }
}
}
//--------------------------------------------------------------------------
int DE_coord::calc(double jd, int center){ //坐标计算
if (jd<h.JD1 || jd>h.JD2) return 0;
int i, j;
int blockN = (jd-h.JD1)/h.Ta; //JD所在块号
DE_rBinBlock(fp, h.nn, blockN, block); //读取数据块
//接下来通过插值计算出各天体的位置和速度
//p1是各天体0段0组在块中的起始位置,p2为段数, p3为每组系数个数
double t = jd-block[0]; //子块内起算的时间
for (i=0; i<13; i++) {
    if(!h.p2[i]||!h.p3[i]) { for(j=0;j<6;j++) pv[i][j]=0; continue; }
    int    nC = (i==11 ? 2:3); // 坐标个数,章动2个,其它3个
    double Tb = h.Ta/h.p2[i]; // 每段时间
    int    nD = t/Tb;      // 段号(取整)
    double t2 = t - nD*Tb; // t2段内时间
    int    pD = h.p1[i]-1 + nD*nC*h.p3[i]; // 得到段的起点位置
    double tx = t2*2/Tb - 1;                         //切比雪夫插用的时间介于[-1,1]
    cheby(tx, block+pD, nC, h.p3[i], pv[i]); //插值计算
    for(j=0;j<nC;j++) pv[i][j+nC]*=2/Tb; //还须把速度单位变换回来(求导数起的)
}
// 插值计算结束,到此已算出太阳系质心坐标系中的行星及太阳坐标,其中地球以日月质心坐标给出,月亮以地心坐标给出
// 接下来补算3个坐标(月亮、地球、太阳系质心三者的太阳系质心坐标),保存位置也做一些调整
for (j=0; j<6; j++) {
    // 把月亮地心坐标、天平动、章动、地月质心移到后面去
    pv[15][j] = pv[ 9][j]; // 月亮地心坐标
    pv[14][j] = pv[12][j]; // 此刻天平动
    pv[13][j] = pv[11][j]; // 此时章动
    pv[12][j] = pv[2][j];  // 把地月质心移到pv[12],而准备用pv[2]放地球坐标
    // 计算地球和月亮的太阳系质心坐标
    pv[11][j] =0; //太阳系质心坐标
    pv[2][j] -= pv[9][j]/(1.0+h.emrat); // 把pv[2]地月质心坐标归算为地球坐标(借助月亮位置来计算)
    pv[9][j] += pv[2][j]; // 把月亮的坐标由地心改为太阳系质心
}
for (i=0; i<13; i++ ) // 把km改为AU
   for (j=0; j<6; j++) pv[i][j] /= h.au;
// 如果必要,以下将太阳系质心坐标改为其它中心坐标(如日心或地心)
if (center!=DE_SSBARY) {
   for (i=0; i<13; i++) {
     if (i == center) i++; // 中心天体跳过
     if (i >= 13) continue;
     for (j=0; j<6; j++)
       pv[i][j] -= pv[center][j];
   }
   for(i=0;i<6;i++) pv[center][i] = 0; // 中心天体位置及速度置0
}
return 1;
}
//--------------------------------------------------------------------------
许老师为什么不用瑞士星历表？
Swiss Ephemeris是由Astrodienst基于美国宇航局JPL实验室发布的
DE406星表拓展而来的一个高精度的星历表.它不是一个面对终端用户的产
品,而是面对天文软件程序开发者的一个工具集.它给我们提供了如下一
些计算功能:行星和恒星的计算;食现象和行星现象的计算;日期和时间转
换功能;初始化,安装和关闭的功能;黄道十二宫的计算以及一些辅助功
能.
Swiss Ephemeris与其他星历表相比,具有如下一些优势:
1.Swiss Ephemeris是在NASA的JPL实验室最新的行星和月亮表DE405和
DE406基础上拓展而来的.DE405星表的时间跨度为公元前3000年到公元
3000年共6000年,它的数据总共需要550M的磁盘空间.DE406是它的压缩
版本并且提高了精度.这些数据被Astrodient利用复杂的压缩技术进一步
压缩,对于6000年的行星数据只需5M磁盘空间,月亮数据需要13M空间,
并且达到了0.001角秒的精度.此外他们还将时间跨度拓展为从公元前5400
年到公元5400年.
2.完成了从惯性参照系到当以当日真春分点为原点的天文坐标系的转
换.所有需要考虑的修正如相对论性偏差,光在太阳重力场中的弯曲等都
以极高的精度加以考虑以使整体转换精度达到0.001角秒.
3.可以利用三种不同的星表数据(DE406,Swiss Ephemeris,Moshier
model)进行计算,该软件包将根据所安装的数据文件自动切换到精度最高
的星历表进行计算.
4.它还包含了所有已知的小行星甚至是假想天体的计算.它包含了
55000颗小行星的数据.
20
5.它不仅精度高,并且运算速度很快.在一台1000MHz的奔三处理器
的测试机器上,计算10000颗星星的位置只需9秒钟.

总觉得，瑞士星历表没有您想像的那么好用，原因有3个：(1)代码太长不易加
入自己的工程中；(2)不易搞清楚它的坐标的具体函义，及其调用方法，
我在使用瑞士星历表时，仍然须要结合《2008年天文年历》、elm/mpp02、
VSOP87、DE405/406星历表、岁差章动计算才能彻底知道瑞士星历表调用
参数的关系。注意：如果不考虑文件存取的问题，DE405/DE406比瑞士星历表
的计算要简单得多!!!!!说实在，文件存取及相关的数据结构，是计算机专业
的事。(3)瑞士星历表把星历扩展到+-5400年没有太大的意义，如果有必要，
DE406当时也会提供这么长的数据，因为数据积分方法不适合长时间的外推。
我曾用科威尔方法积分，用80位的长双精度表示数字(不是64位的)，发现
几千年后的截断积累误差非常大，不是这毫角秒数量级，可能达到角秒数
量级，在这种情况下计算太阳位置，千后以后的误差可能达到1分钟，如果
是这样，远期的的星历推算还不如使用VSOP87。实际上，从数值积分的原
理看来，如果外推10年产生1ms的误差，那么外推100年的误差就不是10ms，
而是更大，比如30ms，千年后可能达到角秒，误差积累与时间不是线性关系
的！所以说，瑞士星历表扩展了范围，只是为了提高适用范围，并不是说
它的精度很高。

我的想法：如果不打算做太专业的计算，使用vsop87算法就可以了，改正一下坐标的起算点(只是一个常数)，计算地内行星，近几十年的精度高达10ms(这是指最大可能的误差)，平均误差大约为几毫角秒，瑞士星历表的精度仅仅高了3至5倍。当然vsop87计算地外行星，误差大一些。如果想做专业一些的计算，很有必要了解DE星历表，瑞士星历表倒不是很有必要(它本身也在使用DE星历表)，因为《中国天文年历》已经开始使用DE星历表。顺便说明一下，DE星历表使用的是它的专用坐标系(ICRF)，swiss把它旋转到了J2000惯性平黄道及J2000平赤道坐标系，二者有几十毫角秒的差异，中国天文年历也使用这个J2000坐标系。
不要过份介意远期的精度！几千年后的星历精度一般认为只能精钟到分钟，不可能精确到毫秒。
DE星历表被很多国家采用，我觉得不单单是因为它的精度高，还有一个原因是它的算法简单，维护起来容易，DE102、DE200、DE403、DE405、DE406等等，算法几乎相同。DE星历表的算法难度与vsop87差不多，只是DE星历表体积很大，需要一些文件操作，造成了一些麻烦


以上程序的javascript调用测试结果
寿星星历activeX测试
测试1：Javascript仿c语言printf()格式化输出函数
printf("%s %4.1d %6e", "test", -0.058, 0.00008233458); 结果如下
test -0.1 8.233458e-5

测试2：activeX路径：E:\天文包\剑\万年历\actX\

测试3：解方程组
　x+y+z=6
　x+2y+z=8
　x+y+3z=12
　结果：(x,y,z)=1,2,3

=============================行星测试==================================
测试4：JD=36525(J2000起算),地球J2000坐标(VSOP87日心球面)：
　黄经: 99°43'13.746"黄纬:- 0°00'45.554"距离:0.9833513122
　黄经: 99°43'13.692"黄纬:- 0°00'45.555"距离:0.9833513125修正了平黄经

测试5：JD=36525(J2000起算),地球J2000坐标(swiss地心球面)：
　黄经: 99°43'13.692"黄纬:- 0°00'45.559"距离:0.9833513062

测试6：JD=36525(J2000起算),地球J2000坐标(DE406地心球面)：
　黄经: 99°43'13.694"黄纬:- 0°00'45.556"距离:0.9833513065

=============================月历测试==================================
测试7：JD=2926.5(J2000起算) 月球J2000坐标(ELP/mpp02 fit DE405 地心球面)
　黄经: 256°47'45.133"黄纬:- 4°52'10.387"距离:401817.6342
测试8：JD=2926.5(J2000起算) 月球J2000坐标(swiss, 地心球面)
　黄经: 256°47'45.133"黄纬:- 4°52'10.342"距离:401817.6341
测试9：JD=2926.5(J2000起算) 月球J2000坐标(DE406, 地心球面)
　黄经: 256°47'45.134"黄纬:- 4°52'10.386"距离:401817.6337
　黄经: 256°47'45.134"黄纬:- 4°52'10.341"距离:401817.6337 补0.047ms交角
视黄经: 256°54'36.317"黄纬:- 4°52'14.134"距离:401817.6711

测试10：取swiss星历表中编号为1的星体名称：Moon
测试11：取swiss星历表中 2008年 deltatT： 65.48秒
测试12：DE406测试 JED 2440400.5 全部星体太阳系质心(SSB)坐标
水星      0.3572602064473    -0.0915490424305   -0.0859810399869 AU
          0.0033678456622     0.0248893428422    0.0129440715868 AU/日
金星      0.6082494331856    -0.3491324431959   -0.1955443457854 AU
          0.0109524201099     0.0156125067399    0.0063288764517 AU/日
地球      0.1160247289670    -0.9265813177178   -0.4017930667289 AU
          0.0168043165197     0.0017451662484    0.0007570133973 AU/日
火星     -0.1146885824391    -1.3283665308335   -0.6061551894194 AU
          0.0144820048079     0.0002372854924   -0.0002837498361 AU/日
木星     -5.3842094069921    -0.8312476561611   -0.2250947570335 AU
          0.0010923632912    -0.0065232941912   -0.0028230122672 AU/日
土星      7.8898899338228     4.5957107269260    1.5584315167251 AU
         -0.0032172034911     0.0043306322336    0.0019264174638 AU/日
天王星  -18.2699008149783    -1.1627115802190   -0.2503695407426 AU
          0.0002215401656    -0.0037676535582   -0.0016532438049 AU/日
海王星  -16.0595450919245   -23.9429482908795   -9.4004227803540 AU
          0.0026431227916    -0.0015034920808   -0.0006812710049 AU/日
冥王星  -30.4878221121555    -0.8732454301967    8.9112969841848 AU
          0.0003225621959    -0.0031487479276   -0.0010801779316 AU/日
日心月    0.1152165516391    -0.9285759477194   -0.4028803293898 AU
          0.0174054013363     0.0015777207878    0.0006714512523 AU/日
太阳      0.0000000000000     0.0000000000000    0.0000000000000 AU
          0.0000000000000     0.0000000000000    0.0000000000000 AU/日
SSB      -0.0045025081562    -0.0007670747009   -0.0002660568052 AU
          0.0000003517482    -0.0000051776254   -0.0000022291019 AU/日
EMB       0.1160149091392    -0.9266055536404   -0.4018062776070 AU
          0.0168116200522     0.0017431316880    0.0007559737671 AU/日
地心月   -0.0008081773279    -0.0019946300016   -0.0010872626608 AU
          0.0006010848167    -0.0001674454606   -0.0000855621450 AU/日
JED 2440400.5 的DE405初始标准坐标,ICRF参考架(不是J2000赤道&分点),JPL说明书中的:
水星     0.35726020644727541518 -0.09154904243051842990 -0.08598103998694037053 位置(AU)
         0.00336784566219378527  0.02488934284224928990  0.01294407158679596663 速度(AU/日)
金星     0.60824943318560406033 -0.34913244319590053792 -0.19554434578540693592
         0.01095242010990883868  0.01561250673986247038  0.00632887645174666542
EMB      0.11601490913916648627 -0.92660555364038517604 -0.40180627760698804496
         0.01681162005220228976  0.00174313168798203152  0.00075597376713614610
火心    -0.11468858243909270380 -1.32836653083348816476 -0.60615518941938081574
         0.01448200480794474793  0.00023728549236071136 -0.00028374983610239698
木星    -5.38420940699214597622 -0.83124765616108382433 -0.22509475703354987777
         0.00109236329121849772 -0.00652329419119226767 -0.00282301226721943903
土星     7.88988993382281817537  4.59571072692601301962  1.55843151672508969735
        -0.00321720349109366378  0.00433063223355569175  0.00192641746379945286
天王星 -18.26990081497826660524 -1.16271158021904696130 -0.25036954074255487461
         0.00022154016562741063 -0.00376765355824616179 -0.00165324380492239354
海王星 -16.05954509192446441763 23.94294829087950141524 -9.40042278035400838599
         0.00264312279157656145 -0.00150349208075879462 -0.00068127100487234772
冥王星 -30.48782211215555045830 -0.87324543019672926542  8.91129698418475509659
         0.00032256219593323324 -0.00314874792755160542 -0.00108017793159369583
太阳     0.00450250884530125842  0.00076707348146464055  0.00026605632781203556
        -0.00000035174423541454  0.00000517762777222281  0.00000222910220557907
月亮    -0.00080817732791148419 -0.00199463000162039941 -0.00108726266083810178
         0.00060108481665912983 -0.00016744546061515148 -0.00008556214497398616