-- 作者:春光
-- 发布时间:2006/6/19 13:28:00
--
月球位置计算很复杂,就摄动项就有1650项之多。要想准确计算其位置并不容易。现在介绍这方面的书是很少的。就连和它有关的天体力学也是很不容易找到。我在网上找到一些关于这些资料,请参考。 计算月球位置的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;//计算反映地球轨道偏心率变化 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////计算月球地心黄经周期项;////////////////////////////////////////////////
|
-- 作者:春光
-- 发布时间:2006/6/19 13:29:00
--
SUML=0; for(i=0;i<=59;i++) { SIN1=La[i]*D+Lb[i]*M+Lc[i]*M1+Ld[i]*F; SUML=SUML+Sl[i]*0.000001*sin(SIN1*DE)*pow(E,fabs(Lb[i])); }; ////////////////////////////////////////////////////////////////////////////////////////////////////////////// ///////////////////////////////计算月球地心黄经/////////////////////////////////////////////////////////////// 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[i]*D+Bb[i]*M+Bc[i]*M1+Bd[i]*F; SUMB=SUMB+Sb[i]*0.000001*sin(SIN2*DE)*pow(E,fabs(Lb[i])); }; ////////////////////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////计算月球地心黄纬//////////////////////////////////////////////////// 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[i]*D+Lb[i]*M+Lc[i]*M1+Ld[i]*F; SUMR=SUMR+Sr[i]*0.001*cos(COS1*DE)*pow(E,fabs(Lb[i])); }; ///////////////////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////计算月球地心距离/////////////////////////////////////////////////////// //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); } 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; }; ///////////////////////////////////////////////////////////////////////////////////////////////////// // 名称:太阳位置计算程序 // 作者:胡铂(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); }
|