#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();
}