/*----------------------------------------------------------------
卫星绕地球运动的模拟程序(地心坐标系并且设卫星的质量极小)
<改进29楼程序的时分秒处理>
----------------------------------------------------------------*/
#include <stdio.h>
#include <math.h>

const double dt = 10e-3; //步长10毫秒(可调整)
const double TPI=2*3.1415926535897932;
const double GM=6.672e-11 * 5.976e+24;
const double Re=6371000; //地球平均半径(米);

double g(double r) //距地心r处的重力加速度
{
return GM/(r*r);
}

void hms(double tim,int*hh,int*mm,double*ss)
// 将tim中的秒数化成 hh小时 mm分 ss秒
{
*hh = (int)(tim/3600);
tim = tim - *hh*3600 ;
*mm = (int)(tim / 60);
*ss = tim - *mm * 60 ;
}

int main( )
{
double H; //卫星距地球表面的高度
double V; //卫星相对于地心的速率
double alpha; //卫星运动的仰角
double tim=0; //计时器(秒)
double the=0; //卫星的真近点角
double S=0.0; //卫星的当前里程
double r; //卫星的地心距离
double dth; //真近点角的增量
int second=0; //存放前一秒钟
int hh,mm; //小时、分
double ss; //秒(含有小数)
alpha=0; //仰角的缺省值

{ //这些红色代码可验证本模拟程序,删除掉也可以
double Hmax,Hmin,Hav,aa,bb,cc,ee,Vmax,Vmin,T;
printf("近地点高度(千米)、远地点高度(千米)＝");
scanf("%lf%*c%lf",&Hmin,&Hmax);
Hmin*=1e3;Hmax*=1e3;
Hav=(Hmax+Hmin)/2;
aa=Re+Hav;
cc=(Hmax-Hmin)/2;
bb=sqrt(aa*aa-cc*cc);
ee=cc/aa;
printf("\nａ＝%.0f千米，ｂ＝%.0f千米，ｃ＝%.0f千米，ｅ＝%.3f\n\n",
aa/1e3,bb/1e3,cc/1e3,ee);
printf(" %.3f \n",bb*bb/aa/1e3);
printf("极坐标方程为ｒ＝────────\n");
printf(" １＋%0.3fcosθ\n\n",ee);

Vmax=sqrt(GM*(2/(aa-cc)-1/aa)); //用活力公式计算近地点速度
Vmin=sqrt(GM*(2/(aa+cc)-1/aa)); //用活力公式计算远地点速度
printf("Ｖmax＝%.3f(米/秒)，Ｖmin＝%.3f(米/秒)\n\n",Vmax,Vmin);

T=TPI/sqrt(GM)*pow(aa,1.5); //理论周期(秒)
hms(T,&hh,&mm,&ss);
printf("Ｔ＝%02d:%02d:%.3f\n\n",hh,mm,ss);
printf("press Enter key to continue....");
getchar();getchar();
}

printf("\n高度(千米) 速度(米/秒) 仰角(度) ＝");
scanf("%lf%*c%lf%*c%lf",&H,&V,&alpha);
alpha*=TPI/360;
H*=1000;
r=Re+H;
printf("\n高度(千米) 速度(米/秒) 仰角(度) θ角(度) hh:mm:ss.cc 里程(千米)\n");
do
{
double ds=V*dt;
double sa=sin(alpha),ca=cos(alpha),gt=g(r)*dt;
double r0=r;
r=sqrt(r*r+ds*ds+2*r*ds*sa)-0.5*gt*dt;
dth=V*ca*dt/r;//坐标系随卫星旋转dθ角
S+=sqrt(pow(r-r0,2)+r*r0*dth*dth);
V-=gt*sa; //t＋dt时刻的速度
alpha+=dth-gt*ca/V; //刷新仰角
the+=dth; //刷新卫星的近点角
tim+=dt; //刷新计时器的秒数
hms(tim,&hh,&mm,&ss);
if(fabs(ss-second)>=1){ second=(int)ss;
printf("%-11.2f%-12.3f%-9.3f%-9.3f%02d:%02d:%5.2f %-9.2f\n",
(r-Re)/1e3,V,alpha/TPI*360,the/TPI*360,hh,mm,ss, S/1e3);}
}
while(the<TPI);
S*=TPI/the;
tim*=TPI/the;
printf("轨道周长＝%.3f千米\n",S/1e3);
printf("轨道周期＝%.3f秒\n",tim);
printf("平均速度＝%.3f米/秒\n",S/tim);
return 0;
}

你初中学过积分了吗?

基于固定的直角坐标系的小质量卫星轨道模拟程序设计

[题]地球半径R=6378千米，引力常数与地球质量乘积GM = 398603×10⁶(千米³/千秒²)。用固定的
直角坐标系模拟小质量卫星的运行轨道，仅考虑地球对卫星的吸引力。卫星轨道的近地点高Hmin和
远地点高Hmax作为输入参数。
[解]就本题而言，轨道形状是已知的：椭圆。像解析几何那样，把二维直角坐标系的原点定义在
椭圆的对称中心，从而让轨道的方程为
(x/a)²＋(y/b)²＝1
式中，半长径a＝R＋(Hmax＋Hmin)/2
半短径b＝(a²－c²)^1/2
焦距c＝(Hmax－Hmin)/2
我们将地球放到右焦点上，即令地心的直角坐标为(c ,0)。这样一来，近地点直角坐标为(a ,0)，
远地点则为(-a ,0)。下面制订ｔ→ｔ+△ｔ时段内，卫星的飞行轨迹以便建立迭代方案。
设ｔ时刻卫星的位置矢量是r、速度矢量是v，ｔ+△t时刻到达r＋△r。
最终，确立下列循环求解方案

⑴输入(或由理论确定)卫星的初始位置r与初速度v

-GM·(r－re)
⑵计算加速度a₁(r)＝─────── 式中re表示地心的位置矢量
|r－re|³

⑶粗算ｔ+△ｔ时刻卫星位置 r＋＝v△t (暂不考虑含△t平方的项)

⑷再算此位置的加速度a₂(r)，公式不变，但卫星的位置r有变化

⑸求出速度增量△v＝(a₁＋a₂)·(△t/2)

⑹ｔ+△t时刻的位置 r＋＝△v·(△t/2)

⑺ｔ+△t时刻的速度 v＋＝△v

⑻如果还没有绕行完一圈，则回到第⑵步

/*----------------------------------------------------------------
　　 小卫星绕地球运动的模拟程序(基于二维直角坐标系)

程序流程：输入近地点和远地点高度，由解析几何确定椭圆轨道的
　　方程，再用“活力公式”确定卫星过这两个点时的理论速度。模拟
　　程序将地心设定在右焦点(c,0)把近地点(a,0)作为模拟运动的初始
　　位置，把近地点的理论速度Vmax作为初速度(仰角等于零)。卫星在
　　XOY平面内沿椭圆轨道逆时针运动。

　　　　　 以下用的是“千米•吨•千秒”制
-----------------------------------------------------------------*/
#include <stdio.h>
#include <math.h>

const double dt = 10e-6; //步长10毫秒
const double TPI=2*3.1415926535897932;
const double GM=398603e6;//6.672e-11*5.976e+21吨
const double Re=6378; //地球平均半径6371(千米)

struct Vector
{ double x,y; } re; //地心的直角坐标

double distance(struct Vector rp,struct Vector rs)
// 从源点rs到场点rp的距离
{
return sqrt(pow(rp.x-rs.x,2)+pow(rp.y-rs.y,2));
}

struct Vector acc( struct Vector r )
// 由地心引力产生的加速度(矢量)
{ static
struct Vector acc;
double temp = -GM/pow(distance(r,re),3);
acc.x = temp*(r.x-re.x);
acc.y = temp*(r.y-re.y);
return acc;
}

void hms(double sec,int*hh,int*mm,double*ss)
// 将 sec 中的秒数化成 hh 小时 mm 分 ss 秒
{
*hh = (int)(sec/3600);
sec = sec - *hh*3600 ;
*mm = (int)(sec / 60);
*ss = sec - *mm * 60 ;
}

double Atan2( double y,double x )
// 带累计功能的反正切(可大于2π)
{
static double old=0,turns=0;
double now=atan2(y,x);
if(now<old)turns+=TPI;
old=now;
return now+turns;
}

int main( )
{
struct Vector rr,vv,dv,ac1,ac2;
double H; //卫星距地球表面的高度(千米)
double V; //卫星相对于地心的速率(米/秒)
double alpha; //卫星运动的仰角
double tim=0; //计时器(秒)
double the=0; //卫星的真近点角
double S=0.0; //卫星的当前里程(千米)
double r; //卫星的地心距离(千米)
int second=0; //存放前一秒钟
int hh,mm; //小时、分
double ss; //秒(含小数部)
double dt2=dt/2; //步长之半
double Hmax,Hmin,Hav,aa,bb,cc,ee,Vmax,Vmin,T;

printf("近地点高度(千米)、远地点高度(千米)＝");
scanf("%lf%*c%lf",&Hmin,&Hmax);
Hav=(Hmax+Hmin)/2;
aa=Re+Hav;
cc=(Hmax-Hmin)/2;
re.x=cc;re.y=0;
bb=sqrt(aa*aa-cc*cc);
ee=cc/aa;
printf("\nａ＝%.0f千米，ｂ＝%.0f千米，ｃ＝%.0f千米，ｅ＝%.4f\n\n",aa,bb,cc,ee);
printf(" %.3f千米 \n",bb*bb/aa);
printf("极坐标方程为ｒ＝────────\n");
printf(" 1＋%0.4fcosθ\n\n",ee);
Vmax=sqrt(GM*(2/(aa-cc)-1/aa)); //用活力公式计算近地点速度
Vmin=sqrt(GM*(2/(aa+cc)-1/aa)); //用活力公式计算远地点速度
printf("Ｖmax＝%.3f(米/秒)，Ｖmin＝%.3f(米/秒)\n\n",Vmax,Vmin);
T=TPI/sqrt(GM)*pow(aa,1.5); //理论周期(千秒)
hms(T*1e3,&hh,&mm,&ss);
printf("Ｔ＝%02d:%02d:%.3f\n\n",hh,mm,ss);
printf("press Enter key to continue....");
getchar();getchar();

alpha=0; //近地点处仰角等于零
r=Re+Hmin; //近地点到地心(千米)
rr.x=cc+r; rr.y=0; //得到位置矢量的初值
vv.x=0; vv.y=Vmax; //得到速度矢量的初值
printf("\n高度(千米) 速度(米/秒) 仰角(度) θ角(度) hh:mm:ss.cc 里程(千米)\n");
do
{
double product;
struct Vector r0=rr;
ac1 = acc(rr); //起点加速度矢量
rr.x+=vv.x*dt;
rr.y+=vv.y*dt;
ac2 = acc(rr); //终点加速度矢量(有误差)
dv.x=(ac1.x+ac2.x)*dt2;
dv.y=(ac1.y+ac2.y)*dt2;
rr.x+=dv.x*dt2;
rr.y+=dv.y*dt2;
vv.x+=dv.x;
vv.y+=dv.y;
product=vv.x*(rr.x-re.x)+vv.y*(rr.y);
r=distance(rr,re); //t＋dt时刻地心距
H=r-Re; //t＋dt时刻的高度
V=sqrt(vv.x*vv.x+vv.y*vv.y);//t＋dt时刻的速度
alpha=asin(product/(r*V)); //t＋dt时刻的仰角
the=Atan2(rr.y,rr.x-re.x); //t＋dt时刻近点角
S+=distance(rr,r0); //t＋dt时刻的里程
tim+=dt; //t＋dt时刻千秒数
hms(tim*1e3,&hh,&mm,&ss);
if(fabs(ss-second)>=1){ second=(int)ss;
printf("%-11.2f%-12.3f%-9.3f%-9.3f%02d:%02d:%5.2f %-9.2f\n",
H,V,alpha/TPI*360,the/TPI*360,hh,mm,ss,S); }
}
while(the<TPI);
S*=TPI/the;
tim*=TPI/the;
printf("轨道周长＝%.3f千米\n",S);
printf("轨道周期＝%.3f秒\n",tim*1e3);
printf("平均速度＝%.3f米/秒\n",S/tim);
return 0;
}

/*
本程序的准确度高于31楼的，在轨道形状、轨道周期
和能量守恒诸方面都与理论值吻合得非常好。这就为
多天体相互绕行模拟程序的设计奠定了一个良好基础
*/

[此贴子已经被作者于2007-11-21 20:16:05编辑过]