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); //asc