计算月球位置的c语言程序

计算月球位置的c语言程序
////////////////////////////////////////////////////////////////////////////////////////////////////////
//                              名称：月球位置计算
//                              作者：胡铂(http://hubble.lamost.org)
//                              日期：2004－09－29
//                              说明：根据北京天文同好会提供的《Astronomy Algrithms》翻译版实现，并得到了
//                                    zjuglr的帮助。
//
///////////////////////////////////////////////////////////////////////////////////////////////////////
#include "math.h"
#include "stdio.h"
#define DE 3.141592654/180  
//////////////////////////////////计算儒略日历书时//////////////////////////////////////////////////////////////
float jde(int Y,int M,int D,int hour,int min,int sec)
{
int f,g;
double mid1,mid2,J,JDE,A;
 if(M>=3)
 {
  f=Y;
  g=M;
 }
 if(M==1||M==2)
 {
  f=Y-1;
  g=M+12;
 };
 mid1=floor(365.25*f);
 mid2=floor(30.6001*(g+1));
 A=2-floor(f/100)+floor(f/400);
 J=mid1+mid2+D+A+1720994.5;
 JDE=J+hour/24+min/1440+sec/86400;
 return JDE;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////
};
void main(void)
{
/////////////////////////////////////////变量定义///////////////////////////////////////////////////////////////////
int i,year,month,day,hour,min,sec;
double JDE,T,L1,D,M,M1,F,A1,A2,A3,E,SUML,lamda,SUMB,beta,SUMR,SIN1,SIN2,COS1,Dist;
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////数据/////////////////////////////////////////////////////////////////
static double La[60]={0,2,2,0,0,0,2,2,2,2,0,1,0,2,0,0,4,0,4,2,2,1,1,2,2,4,2,0,2,2,1,2,0,0,2,2,2,4,0,3,2,4,0,2,2,2,4,0,4,1,2,0,1,3,4,2,0,1,2,2};
static double Lb[60]={0,0,0,0,1,0,0,-1,0,-1,1,0,1,0,0,0,0,0,0,1,1,0,1,-1,0,0,0,1,0,-1,0,-2,1,2,-2,0,0,-1,0,0,1,-1,2,2,1,-1,0,0,-1,0,1,0,1,0,0,-1,2,1,0,0};
static double Lc[60]={1,-1,0,2,0,0,-2,-1,1,0,-1,0,1,0,1,1,-1,3,-2,-1,0,-1,0,1,2,0,-3,-2,-1,-2,1,0,2,0,-1,1,0,-1,2,-1,1,-2,-1,-1,-2,0,1,4,0,-2,0,2,1,-2,-3,2,1,-1,3,-1};
static double Ld[60]={0,0,0,0,0,2,0,0,0,0,0,0,0,-2,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,0,0,0,-2,2,0,2,0,0,0,0,0,0,-2,0,0,0,0,-2,-2,0,0,0,0,0,0,0,-2};
static double Sl[60]={6288774,1274027,658314,213618,-185116,-114332,58793,57066,53322,45758,-40923,-34720,-30383,15327,-12528,10980,10675,10034,8548,-7888,-6766,-5163,4987,4036,3994,3861,3665,-2689,-2602,2390,-2348,2236,-2120,-2069,2048,-1773,-1595,1215,-1110,-892,-810,759,-713,-700,691,596,549,537,520,-487,-399,-381,351,-340,330,327,-323,299,294,0};
static double Sr[60]={-20905355,-3699111,-2955968,-569925,48888,-3149,246158,-152138,-170733,-204586,-129620,108743,104755,10321,0,79661,-34782,-23210,-21636,24208,30824,-8379,-16675,-12831,-10445,-11650,14403,-7003,0,10056,6322,-9884,5751,0,-4950,4130,0,-3958,0,3258,2616,0,-2117,2354,0,0,0,0,0,0,0,-4421,0,0,0,0,1165,0,0,8752};
static double Sb[60]={5128122,280602,277693,173237,55413,46271,32573,17198,9266,8822,8216,4324,4200,-3359,2463,2211,2065,-1870,1828,-1794,-1749,-1565,-1491,-1475,-1410,-1344,-1335,1107,1021,833,777,671,607,596,491,-451,439,422,421,-366,-351,331,315,302,-283,-229,223,223,-220,-220,-185,181,-177,176,166,-164,132,-119,115,107};
static double Ba[60]={0,0,0,2,2,2,2,0,2,0,2,2,2,2,2,2,2,0,4,0,0,0,1,0,0,0,1,0,4,4,0,4,2,2,2,2,0,2,2,2,2,4,2,2,0,2,1,1,0,2,1,2,0,4,4,1,4,1,4,2};
static double Bb[60]={0,0,0,0,0,0,0,0,0,0,-1,0,0,1,-1,-1,-1,1,0,1,0,1,0,1,1,1,0,0,0,0,0,0,0,0,-1,0,0,0,0,1,1,0,-1,-2,0,1,1,1,1,1,0,-1,1,0,-1,0,0,0,-1,-2};
static double Bc[60]={0,1,1,0,-1,-1,0,2,1,2,0,-2,1,0,-1,0,-1,-1,-1,0,0,-1,0,1,1,0,0,3,0,-1,1,-2,0,2,1,-2,3,2,-3,-1,0,0,1,0,1,1,0,0,-2,-1,1,-2,2,-2,-1,1,1,-1,0,0};
static double Bd[60]={1,1,-1,-1,1,-1,1,1,-1,-1,-1,-1,1,-1,1,1,-1,-1,-1,1,3,1,1,1,-1,-1,-1,1,-1,1,-3,1,-3,-1,-1,1,-1,1,-1,1,1,1,1,-1,3,-1,-1,1,-1,-1,1,-1,1,-1,-1,-1,-1,-1,-1,1};
///////////////////////////////////////////////////////////////////////////////////////////////////////////////

/////////////////////////////////////计算日期和时间///////////////////////////////////////////////////////////
year=2004;
month=9;
day=4;
hour=0;
min=0;
sec=0;
//////////////////////////////////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////计算时间日期的儒略日历书时//////////////////////////////////////////////////
JDE=jde(year,month,day,hour,min,sec);
/////////////////////////////////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////计算自J2000.0开始的儒略世纪数///////////////////////////////////////
T=(JDE-2451545)/36525;
/////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////
L1=218.3164477+481267.88123421*T-0.0015786*T*T+T*T*T/538841-T*T*T*T/65194000;//月亮的平黄经
D=297.8501921+445267.1114034*T-0.0018819*T*T+T*T*T/545868-T*T*T*T/113065000;//月亮的平均太阳距角
M=357.5291092+35999.0502909*T-0.0001536*T*T+T*T*T/24490000;//太阳的平近点角
M1=134.9633964+477198.8675055*T+0.0087414*T*T+T*T*T/69699-T*T*T*T/14712000;//月亮的平近点角
F=93.2720950+483202.0175233*T-0.0036539*T*T-T*T*T/3526000+T*T*T*T/863310000;//月亮的黄纬参量(由升交点起算的月球平均距离)
A1=119.75+131.849*T;//金星的摄动
A2=53.09+479264.290*T;//木星的摄动
A3=313.45+481266.484*T;
E=1-0.002516*T-0.0000074*T*T;//计算反映地球轨道偏心率变化
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////计算月球地心黄经周期项；////////////////////////////////////////////////

SUML=0;
for(i=0;i<=59;i++)
{
SIN1=La*D+Lb*M+Lc*M1+Ld*F;
SUML=SUML+Sl*0.000001*sin(SIN1*DE)*pow(E,fabs(Lb));
};
//////////////////////////////////////////////////////////////////////////////////////////////////////////////

///////////////////////////////计算月球地心黄经///////////////////////////////////////////////////////////////
lamda=L1+SUML+(3958*sin(A1*DE)+1962*sin((L1-F)*DE)+318*sin(A2*DE))/1000000;
lamda=fmod(lamda,360);
//////////////////////////////////////////////////////////////////////////////////////////////////////////////

/////////////////////////////计算月球地心黄纬周期项///////////////////////////////////////////////////////////
SUMB=0;
for(i=0;i<=59;i++)
{
SIN2=Ba*D+Bb*M+Bc*M1+Bd*F;
SUMB=SUMB+Sb*0.000001*sin(SIN2*DE)*pow(E,fabs(Lb));
};
//////////////////////////////////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////计算月球地心黄纬////////////////////////////////////////////////////
beta=SUMB+(-2235*sin(L1*DE)   //(modified)
   +382*sin(A3*DE)+175*sin((A1-F)*DE)
   +175*sin((A1+F)*DE)+127*sin((L1-M1)*DE)
   -115*sin((L1+M1)*DE))/1000000;
//////////////////////////////////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////计算月球地心距离周期项////////////////////////////////////////////////
SUMR=0;
for(i=0;i<=59;i++)
{
COS1=La*D+Lb*M+Lc*M1+Ld*F;
SUMR=SUMR+Sr*0.001*cos(COS1*DE)*pow(E,fabs(Lb));
};
/////////////////////////////////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////计算月球地心距离///////////////////////////////////////////////////////
//Dist=385000.56+SUMR/1000;
Dist=385000.56+SUMR;
////////////////////////////////////////////////////////////////////////////////////////////////////////////
printf("\n%d-%d-%d %d:%d:%d\n",year,month,day,hour,min,sec);
printf("JDE=%f\n",JDE);
printf("T=%f\n",T);
printf("L1=%f\n",L1);
printf("D=%f\n",D);
printf("M=%f\n",M);
printf("M1=%f\n",M1);
printf("F=%f\n",F);
printf("A1=%f\n",A1);
printf("A2=%f\n",A2);
printf("A3=%f\n",A3);
printf("E=%f\n",E);
printf("SUML=%f\n",SUML);
printf("SUMB=%f\n",SUMB);
printf("SUMR=%f\n",SUMR);
printf("lamda=%f\n beta=%f\n Dist=%f\n",lamda,beta,Dist);

}

很不错.

能说说使用情况吗? 比如说精度, 适用的年限等.

我对《Astronomy Algrithms》感兴趣，在哪里能下载？

精度大约在秒级，具体的使用年限我不是很清楚，但是在2000年前后的2、3百年肯定可以用。
《Astronomy Algrithms》我也没有找到，不过在北京天文同好会有一系列他们翻译的的月球运动算法。

不错,稍后整理进精华

对了,Astronomy Algrithms这书我在网上没看到有下载的,倒是amazon.com和skypub.com有卖的,国内好像有翻译的发行本,不过不知道名字.

从国外买书很麻烦，国内不知道是否有翻译本，我没有听说过有，如果有的话，不知2Easy兄是否能查到名字，如果有名字就能找得到。

float jde(int Y,int M,int D,int hour,int min,int sec)/*计算儒略日历书时；*/
{
long int f,g,mid1,mid2;
float J,JDE,A;
 if(M>=3)
 {
  f=Y;
  g=M;
 }
 if(M==1||M==2)
 {
  f=Y-1;
  g=M=12;
 };
 mid1=floor(365.25*f);
   mid2=floor(30.6001*(g+1));
      A=2-floor(f/100)+floor(f/400);
        J=mid1+mid2+D+A+1720994.5;
           JDE=J+hour/24+min/1440+sec/86400;
 return JDE;

};

根据Paul Schlyter, Stockholm, Sweden中等精度的行星位置计算方法，实现了其中太阳位置的计算方法。

/////////////////////////////////////////////////////////////////////////////////////////////////////
// 名称：太阳位置计算程序
// 作者：胡铂（http://hubble.lamost.org）
// 日期 2004－09－30
// 说明：根据Paul Schlyter, Stockholm, Sweden的中等精度计算方法实现
////////////////////////////////////////////////////////////////////////////////////////////////////
#include "stdio.h"
#include "math.h"
#define Pi 3.14159265
#define DE Pi/180
/////////////////////////////////////////////////////////////////
float jde(int Y,int M,int D,int hour,int min,int sec)
{
int f,g;
double mid1,mid2,J,JDE,A;
if(M>=3)
{
f=Y;
g=M;
}
if(M==1||M==2)
{
f=Y-1;
g=M+12;
};
mid1=floor(365.25*f);
mid2=floor(30.6001*(g+1));
A=2-floor(f/100)+floor(f/400);
J=mid1+mid2+D+A+1720994.5;
JDE=J+hour/24+min/1440+sec/86400;
return JDE;
}

/////////////////////////////////////变量声明///////////////////////////////
void main()
{ int i,year,month,day,hour,min,sec;
double d,w,a,e,M,oblecl,L,E,xe,ye,r,v,lon,x,y,z,xequt,yequt,zequt,dist,RA,Decl;
// scanf("%d,%d,%d,%d,%d,%d",&year,&month,&day,&hour,&min,&sec);
// printf("\n");
///////////////////////////////测试数据////////////////////////////////////////////////////
year=1990;
month=4;
day=19;
hour=0;
min=0;
sec=0;
////////////////////////////////////////////轨道根数//////////////////////////////////////
d=jde(year,month,day,hour,min,sec)-2451543.5;//相对儒略日；
w=282.9404+4.70935*0.00001*d;//升交点经度
a=1;
e=0.016709-1.151*0.000000001*d;//偏心率
M = 356.0470 + 0.9856002585 * d;//平近点角
oblecl = 23.4393-3.563*0.0000001 * d;//黄赤交角
////////////////////////////////////////////////////////////////////////////////////////////
L=w+M;//太阳的平均精度
L=fmod(L,360);
E=M + (180/Pi) * e * sin(M*DE) * (1 + e * cos(M*DE));//（开普勒方程的近似解）
E=fmod(E,360);
xe=cos(E*DE) - e;//椭圆轨道上的直角坐标x
ye=sin(E*DE) * sqrt(1 - e*e);//椭圆轨道上的直角坐标y
r=sqrt(xe*xe+ye*ye);//距离
v=atan2(ye,xe)*180/Pi;//真近点角
//lon=fmod((v+w),360);//太阳的精度
//太阳的黄道直角坐标
lon=v+w;
x=r*cos(lon*Pi/180);
y=r*sin(lon*Pi/180);
z=0;
//太阳的赤道直角坐标
xequt=x;
yequt=y*cos(oblecl*Pi/180);
zequt = y* sin(oblecl*Pi/180);
//日地距离、赤经赤纬
dist=sqrt(xequt*xequt+yequt*yequt);
RA=atan2(yequt,xequt)*180/Pi;
RA=fmod(RA,360);
Decl=asin(zequt/r)*180/Pi;
Decl=fmod(Decl,360);
////////////////////////////////////////////////////////////////////
printf("d=%f\n",d);
printf("w=%f\n",w);
printf("e=%f\n",e);
printf("M=%f\n",M);
printf("L=%f\n",L);
printf("E=%f\n",E);
printf("oblecl=%f",oblecl);
printf("xe=%f,ye=%f\n",xe,ye);
printf("r=%f,v=%f\n",r,v);
printf("lon=%f\n",lon);
printf("x=%f,y=%f\n",x,y);
printf("xequt=%f,yequt=%f,zequt=%f\n",xequt,yequt,zequt);
printf("dist=%f\n",dist);
printf("RA=%f,Decl=%f\n",RA,Decl);
}

其实参数可以放在一个专门的文件里呀。这样列在程序里很冗长的，引用的时候打开那个文件就可以了，更重要的是可以更改参数。

我只是演示一下如何计算，和如果写在软件里还是不一样的。编程格式什么的我都没有讲究：）

我想得到他们分别在WGS84的坐标怎么进行转换  没找到材料呀！知道的麻烦告诉下！

同问，如何通过楼主给出的计算方法得到月球的地心黄经、黄纬后转换成地心系直角坐标系下的坐标。

//太阳的黄道直角坐标
lon=v+w;
x=r*cos(lon*Pi/180);
y=r*sin(lon*Pi/180);
z=0;
//太阳的赤道直角坐标
xequt=x;
yequt=y*cos(oblecl*Pi/180);
zequt = y* sin(oblecl*Pi/180);
这段程序就是将太阳位置由地心黄道直角坐标转换到地心直角坐标，月球的同理。就是一个不同坐标系的旋转问题。若想知道原理，建议翻一下高中的解析几何，里面有两个直角坐标系旋转的计算。