-- 作者：xjw01
-- 发布时间：2008/8/18 7:52:00
--
以下是我调用swiss星历函数验证算法正确性的程序，你不妨参考一下 #ifndef __sweH
#define __sweH
#include "swephexp.h" //瑞士星历表 /****************************************************************
IMPLIB是C++Builder自带的工具,在命令行中执行它可得到相关dll的LIB连接库
对再瑞士星历表请执行：
(1)命令: IMPLIB -f -c swe32bor.lib 路径\\swedll32.dll
(2)Options->Projects->Directories/Conditionals:
在 conditional define 中加入 USE_DLL
(3)#include "swiss\\\\swephexp.h"是必要的
*****************************************************************/ /*********************
星历计算函数 swe_calc_ut()与swe_calc()
1.swe_calc (double tjd_ut, int p, int iflag, double *x2, char *serr)
2.swe_calc_ut(double tjd_ut, int p, int iflag, double *x2, char *serr)
这个函数的用法相同,同是的是tjd_ut的含意不同
double tjd_ut 是儒略日数,力学时(函数1)或UT时(函数2)
int p 是星体编号,int
double x2[6] 是返回值,含坐标或速度
int iflag 决定计算的方式
char serr[256] 错误信息
iflgret = swe_calc_ut(tjd_ut, p, iflag, x2, serr);
·函数返回置(记为iflgret),整型,正常计算时 iflgret与iflag相同
·如果些计算不能实现，则iflgret与iflag不相同
·如果iflgret不于零，说明计算失败,x2全部值零,错误信息置于serr中
·iflag的值：以在以下值中选取你所需要的，并相加起来
如果没有比特被设置, 即 if iflag == 0,
swe_calc() 计算普通星历(如果书店的)，即天体的地球黄道坐标中的
视位置(黄经,黄纬,距离)，这是相对于当日真分点的。
·如果需要天体的速度，应设置 iflag = SEFLG_SPEED
·对于一些数学意义的点，如月球轨道交点或近地点，就没有视位置，则返回真位置。
·要进行其它类型的计算，则要用到以下参数(取自swephexp.h)，将它们组合("相加"或"|"):
例如 iflag = SEFLG_SPEED | SEFLG_TRUEPOS;
或者 iflag = SEFLG_SPEED + SEFLG_TRUEPOS;
使用这个iflag值, swe_calc() 将计算真位置(即不含光行时)和速度
·iflag的可选值
#define SEFLG_JPLEPH 1L // 使用 JPL 星历
#define SEFLG_SWIEPH 2L // 使用 SWISSEPH 星历, 默认
#define SEFLG_MOSEPH 4L // 使用 Moshier 星历(半解析法)
#define SEFLG_HELCTR 8L // 返回日心坐标
#define SEFLG_TRUEPOS 16L // 返回真位置(不是视位置)
#define SEFLG_J2000 32L // 不含岁差,即返回 J2000分点坐标
#define SEFLG_NONUT 64L // 不含章动，即date平分点
#define SEFLG_SPEED3 128L // 从3点得到速度 (别使用它, SEFLG_SPEED 更快且更精确)
#define SEFLG_SPEED 256L // 高精度速度 (analyt. comp.)
#define SEFLG_NOGDEFL 512L // 关闭引力偏转
#define SEFLG_NOABERR 1024L // 关闭恒星周年光行差
#define SEFLG_EQUATORIAL 2048L // 赤道坐标
#define SEFLG_XYZ 4096L // 直角坐标，而不是球面坐标
#define SEFLG_RADIANS 8192L // 坐标单位是弧度而不是度
#define SEFLG_BARYCTR 16384L // 太阳系质心坐标
#define SEFLG_TOPOCTR (32*1024L) // 地平坐标位置
#define SEFLG_SIDEREAL (64*1024L) // 恒星位置 详见swephexp.h
------------------
计算儒略日数
double swe_julday(int year, int month, int day, double hour, int gregflag);
*******************/ class SWISS{ //瑞士星历表
public:
int p; //星体序号
int UT;
int iflag,iflgret;
char serr[256]; //星体名、错误信息
SWISS(){
swe_set_ephe_path("src\\\\eph"); //路径为空值使用内置的半解析法
p=1;
iflag=0;
UT=0;
}
char* getname(){ //取星体名称
static char snam[10][40],nameK=-1;
nameK++; if(nameK>9) nameK=0;
swe_get_planet_name(p, snam[nameK]);//取出行星p的名称
return snam[nameK];
}
double julday(int year,int mon,int day,double hour){ //返回儒略日数
swe_julday(year,mon,day,hour,1); }
void cal(double jd,double z[6]){ //星历计算
jd+=2451545.0;
if(UT) iflgret = swe_calc_ut(jd, p, iflag, z, serr);
else iflgret = swe_calc (jd, p, iflag, z, serr);
}
double deltat(double jd) { return swe_deltat(jd); } //取deltatT
void close(){swe_close();};
~SWISS(){swe_close();}
}; #endif
|
-- 作者：xjw01
-- 发布时间：2008/8/18 7:53:00
--
#include "src\\\\ptsz.h"
#include "src\\\\swe.h"
#include <conio.h>
void swe_ex1(){ //瑞士星历表使用范例1
char *sp, sdate[AS_MAXCH], snam[40], serr[AS_MAXCH];
int jday = 1, jmon = 1, jyear = 2000;
double jut = 12.0;
double tjd_ut, te, x2[6];
long iflag, iflgret;
int p;
iflag = SEFLG_SPEED;
swe_set_ephe_path("swiss\\\\eph");
while (TRUE) {
printf("\\n日期 (d.m.y) ?");
gets(sdate); /* stop if a period . is entered */
if (*sdate == \'.\') return;
if (sscanf (sdate, "%d%*c%d%*c%d", &jday,&jmon,&jyear) < 1) exit(1);
//把年月日转为儒略日数
tjd_ut = swe_julday(jyear,jmon,jday,jut,SE_GREG_CAL); //注意,这里仍是力学时
//tjd_ut+=swe_deltat(tjd_ut);
printf("日期: %02d.%02d.%d at 0:00 UT时\\n", jday, jmon, jyear);
printf("行星 \\t 黄经 \\t 黄纬 \\t 距离 \\t 黄经速度\\n");
for (p = SE_SUN; p <= SE_CHIRON; p++) { //遍历所有行星
if (p == SE_EARTH) continue;
iflgret = swe_calc_ut(tjd_ut, p, iflag, x2, serr); //计算行星p的坐标
if (iflgret<0) printf("错误:%s\\n", serr); //如果返回负值表明有错误，错误信息在serr
swe_get_planet_name(p, snam);//取出行星p的名称
printf("%10s\\t%11.7f\\t%10.7f\\t%10.7f\\t%10.7f\\n",//列印坐标
snam,x2[0],x2[1],x2[2],x2[3]);
}
}
swe_close();
}
void swe_ex2(){ //ELP与swiss星历比较
double z[6],z2[6],r[6],r2[6],jd,c;
SWISS K; K.p=1;
semiA E; E.readfileABC("G_DE");
int i;
printf("2008年1月份,月球视位置\\r\\n");
printf("0时 视黄经/视赤经 视黄纬/视赤纬 方法\\r\\n");
for(i=0;i<30;i++){
jd=365.25*8+i-0.5; //时间 E.evaluate(jd,z2); //ELP计算
MPP405_FitSwiss(jd,z2); //转到J2000坐标
addB03Prece(jd,z2,0); //加岁差
YQgxc_CD(jd,z2,0); //月球光行时修正
r2[0]=z2[0]; r2[1]=z2[1]; r2[2]=z2[2];
addNutationLon (jd,z2); //补黄经章动得到黄道视坐标
hcConv(jd,r2,r2,1); //得到赤道视坐标 K.iflag = SEFLG_RADIANS; K.cal(jd,z); //swiss黄道坐标
K.iflag += SEFLG_EQUATORIAL; K.cal(jd,r); //swiss赤道坐标 printf("%3d日 %s %s SWISS\\r\\n",i+1, rad2str(z[0],0) ,rad2str(z[1],0));
printf(" %s %s ELP/MPP02_405\\r\\n",rad2str(z2[0],0),rad2str(z2[1],0));
printf(" %s %s SWISS\\r\\n", rad2str(r[0],0), rad2str(z[1],0));
printf(" %s %s ELP/MPP02_405\\r\\n",rad2str(r2[0],0),rad2str(z2[1],0));
printf("\\r\\n");
}
getch();
}
void swe_ex3(){
SWISS K;
double c,c2,nian;
int i;
for(i=0;i<30;i++){
nian=i*1;
c=K.deltat(J2000+365.25*nian)*86400;
c2=deltatT(nian+2000);
printf("%6.1f %f %f %f\\r\\n",nian+2000,c2,c,c2-c);
}
getch();
}
//以下是测试程序
void ex_swe_vsop(){//vsop87与swe的比较(地球的日心坐标比较)
/*
#define SEFLG_JPLEPH 1L // 使用 JPL 星历
#define SEFLG_SWIEPH 2L // 使用 SWISSEPH 星历, 默认
#define SEFLG_MOSEPH 4L // 使用 Moshier 星历(半解析法)
#define SEFLG_HELCTR 8L // 返回日心坐标
#define SEFLG_TRUEPOS 16L // 返回真位置(不是视位置,不含光行时)
#define SEFLG_J2000 32L // 不含岁差,即返回 J2000分点坐标
#define SEFLG_NONUT 64L // 不含章动，即date平分点
#define SEFLG_SPEED 256L // 高精度速度 (analyt. comp.)
#define SEFLG_NOGDEFL 512L // 关闭引力偏转
#define SEFLG_NOABERR 1024L // 关闭恒星周年光行差
#define SEFLG_EQUATORIAL 2048L // 赤道坐标
#define SEFLG_BARYCTR 16384L // 太阳系质心坐标
#define SEFLG_RADIANS (8*1024) //弧度(不是角度制) 对于行星和月球,光行时与周年光行差应同时计算,即TRUEPOS与NOABERR应同时有或无
*/
double dw=M_PI/180;
double z1[6],z2[6];
double jd;
semiA vsop; SWISS K;
vsop.readfileABC("Dear"); K.p=14; //地球
//vsop.readfileABC("Bven"); K.p=3; //金星
//vsop.readfileABC("Bjup"); K.p=5; //木星
K.iflag = SEFLG_HELCTR + SEFLG_TRUEPOS //日心,无光行时
+ SEFLG_J2000 + SEFLG_NONUT //无岁差,无章动
+ SEFLG_NOABERR+ SEFLG_NOGDEFL;//关闭周年光行差和引力偏转(其实日心坐标会自动关闭)
int i,j; const NC=830*2;
double av=0,max=0,c,dy=0,T=6;
double y[NC],x[NC];
for(i=0;i<NC;i++){
//坐标计算
jd=T*i/NC-T/2;
vsop.evaluate(jd*365250.0,z1);
K.cal(jd*365250.0,z2);
double t=jd,t2=t*t,t3=t2*t,t4=t3*t,t5=t4*t;
//以下对地球拟合到DE406,黄经误差:
//修正前,J2000+-3000年1.0",+-2000年0.67",+-1000年0.39",+-500年0.28",+-250年0.17"
//修正后,J2000+-3000年0.6",+-2000年0.25",+-1000年0.12",+-500年0.02",+250年0.012"
//距离与黄纬不再修正,其误差为:
//距离：+-250年3E-8AU,+-1000年1E-7AU,+-3000年1E-7AU
//黄纬：+-250年0.01", +-1000年0.024",+-3000年0.08"
//对地球Date球面坐标的黄经拟合DE406+瑞士星历表
//J2000+-3000年0.48",+-2000年0.22",+-1000年0.09",+-500年0.016",+-250年0.012"
dy=0.195+ 2.92*t -0.1*t2 -0.017*t3 -0.0034*t4-0.13*cos(6+3*t);
z1[0]-=dy/rad;
z1[0] = rad2mrad(z1[0]-precess(t*365250,Pr_p,Pmo_B03));
llr2xyz(z1,z1); //转为直角
Pxyz2Jxyz(jd*365250,z1,z1); //旋转到J2000黄道
xyz2llr(z1,z1); //直角转球面
// dy = 0.2 -0.22*t -0.054*t2 + 0.1*cos(t)*t - 0.11*cos(1.5*t)*t -0.13*cos(3*t); //地球黄经
// dy = 0.07-0.18*t+0.023*t2 -0.026*t3; //金星黄经
/* 土星黄经,误差+-1000年0.5",+-2000年0.9",+-3000年7"*/
/* dy = 0.2862-0.515110*t -0.384146*t2 -0.091959*t3; //土星黄经
if(t<-2.4) dy+=-1224.7 -919*t -171*t2;
if(t>=2.0&&t<2.5) dy+=-211+188*t -41.3*t2;
if(t>=2.5&&t<2.8) dy+=-947.4+729*t -139.718270*t2;
if(t>=2.8) dy+=57-21*t;*/
// z1[0]-=dy/rad;
//printf("VSOP %11.7f\\t%10.7f\\t%10.7f\\r\\n",z1[0]/dw,z1[1]/dw,z1[2]);
//printf("SWISS %11.7f\\t%10.7f\\t%10.7f\\r\\n",z2[0],z2[1],z2[2]);
c=(z1[0]/dw-z2[0])*3600;
y[i]=c,x[i]=t;
// if(i<200)
printf("%3d %11.5f\\r\\n",i-NC/2*0,c*100);
c=fabs(c); if(c>1000) c=fabs(c-3600*360);
av+=c; if(c>max) max=c;
}
printf("\\r\\n%f %f\\r\\n",av/NC,max);
const NB=6; double abc[NB];
// DXSzxec(y+161,x+161,5,NB,abc);
DXSzxec(y,x,NC,NB,abc);
for(i=0;i<NB;i++) printf("%f ",abc[i]);
getch();
} void ex_swe_mpp02(){//ELP/MPP02与swe的比较(月的地心坐标比较)
/*
#define SEFLG_JPLEPH 1L // 使用 JPL 星历
#define SEFLG_SWIEPH 2L // 使用 SWISSEPH 星历, 默认
#define SEFLG_MOSEPH 4L // 使用 Moshier 星历(半解析法)
#define SEFLG_HELCTR 8L // 返回日心坐标
#define SEFLG_TRUEPOS 16L // 返回真位置(不是视位置)
#define SEFLG_J2000 32L // 不含岁差,即返回 J2000分点坐标
#define SEFLG_NONUT 64L // 不含章动，即date平分点
#define SEFLG_SPEED 256L // 高精度速度 (analyt. comp.)
#define SEFLG_NOGDEFL 512L // 关闭引力偏转
#define SEFLG_NOABERR 1024L // 关闭恒星周年光行差
#define SEFLG_EQUATORIAL 2048L // 赤道坐标
#define SEFLG_BARYCTR 16384L // 太阳系质心坐标
*/
double dw=M_PI/180;
double z1[6],z2[6];
double jd;
ELP mpp; SWISS K;
ELP82 elp;
mpp.init(0); K.p=1; //月球
//K.iflag = SEFLG_TRUEPOS //无光行时
// + SEFLG_J2000 + SEFLG_NONUT //无岁差,无章动
// + SEFLG_NOABERR+ SEFLG_NOGDEFL;//关闭周年光行差和引力偏转
K.iflag=0; //天文年历中的视位置
K.iflag+=SEFLG_EQUATORIAL;
int i,j; const NC=50*2;
double av=0,max=0,c,dy,T=1;
double y[NC],x[NC];
for(i=0;i<NC;i++){
//坐标计算
jd=T*i/NC-T/2+0.0;
x[i]=jd; mpp.evaluate(jd*36525,z1,z1+3,1); MPP405_FitSwiss(jd*36525,z1);
addB03Prece(jd*36525,z1,0);
// addNutationLon (jd*36525,z1); //补黄经章动
YQgxc_CD(jd*36525,z1,0);
hcConv(jd*36525,z1,z1,1);
//Dllr2Jllr(jd*36525,z1,z1,0);
K.cal(jd*36525,z2); //printf("ELP %11.7f\\t%10.7f\\t%10.7f\\r\\n",z1[0]/dw,z1[1]/dw,z1[2]);
//printf("SWISS %11.7f\\t%10.7f\\t%10.7f\\r\\n",z2[0],z2[1],z2[2]);
c=(z1[0]/dw-z2[0])*3600;
y[i]=c;
printf("%3d %11.5f %11.5f\\r\\n",i,c,z2[0]);
c=fabs(c); if(c>1000) c-=2*M_PI;
av+=c; if(c>max) max=c;
}
const NB=3; double abc[NB];
DXSzxec(y,x,NC,NB,abc);
for(i=0;i<NB;i++) printf("%f ",abc[i]); printf("\\r\\n");
double max2=0;
for(i=0;i<NC;i++){
double t=1; int j;
for(j=0,c=0;j<NB;j++) c+=abc[j]*t,t*=x[i];
c=y[i]-c; printf("%f\\r\\n",c);
c=fabs(c); if(c>1000) c=fabs(c-360*3600);
if(c>max2) max2=c;
}
printf("\\r\\n%f %f %f",av/NC,max,max2);
getch();
} void ex_swe_ips(){//ips与swe的比较(地球的日心坐标比较)
/*
#define SEFLG_JPLEPH 1L // 使用 JPL 星历
#define SEFLG_SWIEPH 2L // 使用 SWISSEPH 星历, 默认
#define SEFLG_MOSEPH 4L // 使用 Moshier 星历(半解析法)
#define SEFLG_HELCTR 8L // 返回日心坐标
#define SEFLG_TRUEPOS 16L // 返回真位置(不是视位置)
#define SEFLG_J2000 32L // 不含岁差,即返回 J2000分点坐标
#define SEFLG_NONUT 64L // 不含章动，即date平分点
#define SEFLG_SPEED 256L // 高精度速度 (analyt. comp.)
#define SEFLG_NOGDEFL 512L // 关闭引力偏转
#define SEFLG_NOABERR 1024L // 关闭恒星周年光行差
#define SEFLG_EQUATORIAL 2048L // 赤道坐标
#define SEFLG_BARYCTR 16384L // 太阳系质心坐标
*/
double dw=M_PI/180;
double z1[6],z2[6];
double jd;
semiA vsop; SWISS K;
vsop.readfileABC("Jven"); K.p=3; //金星
K.iflag = SEFLG_HELCTR + SEFLG_TRUEPOS //日心,无光行时
+ SEFLG_J2000 + SEFLG_NONUT //无岁差,无章动
+ SEFLG_NOABERR+ SEFLG_NOGDEFL;//关闭周年光行差和引力偏转(其实日心坐标会自动关闭)
int i,j,NC=83*2;
double av=0,max=0,c,dy,T=5;
for(i=0;i<NC;i++){
//坐标计算
jd=T*i/NC-T/2+0.0;
vsop.evaluate(jd*365250.0,z1); xyz2llr(z1,z1);
K.cal(jd*365250.0,z2);
//printf("IPS %11.7f\\t%10.7f\\t%10.7f\\r\\n",z1[0]/dw,z1[1]/dw,z1[2]);
//printf("SWISS %11.7f\\t%10.7f\\t%10.7f\\r\\n",z2[0],z2[1],z2[2]);
c=(z1[0]/dw-z2[0])*3600;
printf("%3d %11.5f\\r\\n",i,c*100);
c=fabs(c); if(c>1000) c=fabs(c-3600*360);
av+=c; if(c>max) max=c;
}
printf("\\r\\n%f %f",av/NC,max);
getch();
}
|
-- 作者：ymy111
-- 发布时间：2008/8/21 14:50:00
--
多谢 下面几个函数有什么分别呢---英文我实在太差（对照金山快译也没有明白）
Private Declare Function swe_houses Lib "swedll32.dll" _
Alias "_swe_houses@36" ( _
ByVal tjd_ut As Double, _
ByVal geolat As Double, _
ByVal geolon As Double, _
ByVal ihsy As Long, _
ByRef hcusps As Double, _
ByRef ascmc As Double _
) As Long \' hcusps must be first of 13 array elements
\' ascmc must be first of 10 array elements
Private Declare Function swe_houses_d Lib "swedll32.dll" _
Alias "_swe_houses_d@24" ( _
ByRef tjd_ut As Double, _
ByRef geolat As Double, _
ByRef geolon As Double, _
ByVal ihsy As Long, _
ByRef hcusps As Double, _
ByRef ascmc As Double _
) As Long \' hcusps must be first of 13 array elements
\' ascmc must be first of 10 array elements Private Declare Function swe_houses_ex Lib "swedll32.dll" _
Alias "_swe_houses_ex@40" ( _
ByVal tjd_ut As Double, _
ByVal iflag As Long, _
ByVal geolat As Double, _
ByVal geolon As Double, _
ByVal ihsy As Long, _
ByRef hcusps As Double, _
ByRef ascmc As Double _
) As Long \' hcusps must be first of 13 array elements
\' ascmc must be first of 10 array elements Private Declare Function swe_houses_ex_d Lib "swedll32.dll" _
Alias "_swe_houses_ex_d@28" ( _
ByRef tjd_ut As Double, _
ByVal iflag As Long, _
ByRef geolat As Double, _
ByRef geolon As Double, _
ByVal ihsy As Long, _
ByRef hcusps As Double, _
ByRef ascmc As Double _
) As Long \' hcusps must be first of 13 array elements
\' ascmc must be first of 10 array elements Private Declare Function swe_houses_armc Lib "swedll32.dll" _
Alias "_swe_houses_armc@36" ( _
ByVal armc As Double, _
ByVal geolat As Double, _
ByVal eps As Double, _
ByVal ihsy As Long, _
ByRef hcusps As Double, _
ByRef ascmc As Double _
) As Long \' hcusps must be first of 13 array elements
\' ascmc must be first of 10 array elements Private Declare Function swe_houses_armc_d Lib "swedll32.dll" _
Alias "_swe_houses_armc_d@24" ( _
ByRef armc As Double, _
ByRef geolat As Double, _
ByRef eps As Double, _
ByVal ihsy As Long, _
ByRef hcusps As Double, _
ByRef ascmc As Double _
) As Long \' hcusps must be first of 13 array elements
\' ascmc must be first of 10 array elements Private Declare Function swe_house_pos Lib "swedll32.dll" _
Alias "_swe_house_pos@36" ( _
ByVal armc As Double, _
ByVal geolat As Double, _
ByVal eps As Double, _
ByVal ihsy As Long, _
ByRef xpin As Double, _
ByVal serr As String _
) As Double
\' xpin must be first of 2 array elements Private Declare Function swe_house_pos_d Lib "swedll32.dll" _
Alias "_swe_house_pos_d@28" ( _
ByRef armc As Double, _
ByRef geolat As Double, _
ByRef eps As Double, _
ByVal ihsy As Long, _
ByRef xpin As Double, _
ByRef hpos As Double, _
ByVal serr As String _
) As Long
\' xpin must be first of 2 array elements
[此贴子已经被作者于2008-8-21 16:26:04编辑过]
|