s.m

北京到广州铁路路线仿真

clear all;
bj=[116.46,39.92,0];
zhzh=[113.7,34.5,0];
wch=[114.32,30.35,0];
chsh=[112.92,28.33,0];
gzh=[113.25,23.13,0];
%选取北京、郑州、武昌、长沙、广州五个城市为站点
%设火车出站从零开始加速，速度达到180km/h即50m/s时，开始匀速运动
%假设加速时间为10分钟
v=50;%m/s
Re=6378137;%米
f=1/298.257;%扁平率
trans=2*pi/360;
%------计算北京到郑州的距离---------
[Xu_bj,Yu_bj,Zu_bj]=user_ecef(bj(1),bj(2),bj(3));
[Xu_zhzh,Yu_zhzh,Zu_zhzh]=user_ecef(zhzh(1),zhzh(2),zhzh(3));
Disbj_zhzh=sqrt((Xu_bj-Xu_zhzh)^2+(Yu_bj-Yu_zhzh)^2+(Zu_bj-Zu_zhzh)^2);
[Xug_zhzh,Yug_zhzh,Zug_zhzh]=ecef_g(bj(1),bj(2),Xu_zhzh,Yu_zhzh,Zu_zhzh,Xu_bj,Yu_bj,Zu_bj);
%计算在北京的用户水平坐标系下的郑州的位置
theta=pi+atan(Yug_zhzh/Xug_zhzh);
% vn=v*sin(theta);
% ve=v*cos(theta);%将速度分解为东向和北向
v0=0;
t0=0;
dt=10;%s
Tacc=600;%s
a=v/Tacc;
while v0<v
    long0=bj(1);
    lat0=bj(2); 
    Rm=Re*(1+2*f-3*f*sin(trans*lat0)^2);
    Rn=Re/((1-f)^2*sin(trans*lat0)^2+cos(trans*lat0)^2);   
    deltalong=(v0*dt+0.5*a*dt)*cos(theta)/(Rm*cos(trans*lat0)*trans);%十秒间隔内经度的增加量
    deltalat=(v0*dt+0.5*a*dt)*sin(theta)/(Rn*trans);%10秒间隔内纬度的增加量
    long0=long0+deltalong;
    lat0=lat0+deltalat;
    t0=t0+dt; 
    v0=v0+a*dt;%
end
   long1=long0;
   lat1=lat0;
    %[Xujs_zhzh,Yujs_zhzh,Zujs_zhzh]=user_ecef(long0,lat0,0);
 % 北京到郑州匀速
 while v0==v0;%m/s
     vn=v0*sin(theta);
     ve=v0*cos(theta);
     Rm=Re*(1+2*f-3*f*sin(trans*lat1)^2);
     Rn=Re/((1-f)^2*sin(trans*lat1)^2+cos(trans*lat1)^2);   
     deltalong=ve*10/(Rm*cos(trans*lat1)*trans);
     deltalat=vn*10/(Rn*trans);
     long1=long1+deltalong;
     lat1=lat1+deltalat;
     [Xu,Yu,Zu]=user_ecef(long1,lat1,0);
     Dis=sqrt((Xu_bj-Xu)^2+(Yu_bj-Yu)^2+(Zu_bj-Zu)^2);
     if Dis<=Disbj_zhzh-0.5*a*Tacc^2
          t0=t0+dt;
          v0=v;
     else
         t2=t0;
         break;
     end     
   end
long2=long1;
lat2=lat1;
while v0>0%减速进郑州站
    Rm=Re*(1+2*f-3*f*sin(trans*lat2)^2);
    Rn=Re/((1-f)^2*sin(trans*lat2)^2+cos(trans*lat2)^2);   
    deltalong=(v0*dt-a*dt*(t0-t2)-0.5*a*dt^2)*cos(theta)/(Rm*cos(trans*lat2)*trans);
    deltalat=(v0*dt-a*dt*(t0-t2)-0.5*a*dt^2)*sin(theta)/(Rn*trans);
    long2=long2+deltalong;
    lat2=lat2+deltalat;
    t0=t0+dt;
    v0=v-a*(t0-t2);
%     [Xu,Yu,Zu]=user_ecef(long2,lat2,0);
%     Dis=sqrt((Xu_bj-Xu)^2+(Yu_bj-Yu)^2+(Zu_bj-Zu)^2);
end