KSP open SSTO V1.0 空天飞机上升轨道计算器
为方便各位KSP爱好者在RSS环境下手搓SSTO，公开电脑里库存的自己练手写的matlab小程序，可以拟定各种参数，粗略估计一架SSTO能否入轨。如果在这里面都算出来不到第一宇宙速度，那在游戏里大概率是入不了轨了。
例如一架250t起飞重量、燃料占比85%的SSTO，用涡喷起飞，6°小角度爬升~10000m，随后推平到3.5°，并加速到~2Ma。假设涡喷比较屑，到~2Ma就无法继续加速，便启动火箭加速至~3Ma填补推力陷阱，直到RBCC启动速度，然后继续3.5°小角度爬升~30000m，加速到6~7Ma。3万米以上缺氧，所以拉起20°并切换到火箭模式入轨。
为保护小绿人不受摧残，限制z向过载3g，轴向过载3.5g。
本程序简化了升阻力系数（认为不随速度变化）、发动机比冲（也不随速度变化），所以仅供参考，不对小绿人的安全负责。
上升轨道如图：
代码如下
clear all
%%用涡喷起飞，小角度爬升~10000m，随后加速到~2Ma，启动火箭加速至~3Ma填补推力陷阱，
%随后RBCC启动，继续小角度爬升~30000m，加速到6~7Ma，拉起切换到火箭模式入轨。
%mode1是涡喷，mode2是RBCC，mode3是火箭
%基本参数
max_g_limit=3.5;%大于该过载，自动下降推力
thrust_turbo=1200000;%推力牛顿
isp_turbo=2500;%比冲s
thrust_rocket=3000000;
isp_rocket=455;
thrust_rbcc=2500000;
isp_rbcc=1000;
V_rocketstart=600;%t填补推力缺口，大于该速度后启动火箭加速至rbcc启动速度然后切换rbcc
V_rbcc=900;%rbcc启动的速度
m=250000;%起飞重量kg
m0=m*0.15;%空重
%aerodynamic
Cd=0.04;%爬升阶段阻力系数
Cl=0.2;%爬升最大升力系数
glide_Cl=1;%滑翔升力系数，如果不小心再入了，用得上这个值
glide_Cd=1;%滑翔阻力系数，如果不小心再入了，用得上这个值
S_wing=500;%翼面积
%launch parameter
pull_up_g=3;%拉起过载
push_down_g=-0.1;%推平过载
push_down_Cl=-0.05;%推平Cl
max_theta=6;%起飞最大爬升角
push_down_height=10000;%改平高度
push_down_angle=3.5;%改平后爬升角
pull_up_height=30000;%二次拉起高度
pull_up_angle=20;%二次拉起爬升角
H=0;X=0;%发射架初始位置，高度m，水平距离m，默认都是0
%air
rho=1.225*2.71828^(-1.389*10^(-4)*H);%空气密度
D=0;L=0;%升阻力初始化设置0，勿动
%Launch setting
dt=0.1;%时间步长
calculate_T=900;%计算时间总长度
Vx=0;Vy=0;V=0;%初始速度
theta=0;%初始化
k=1;
P=zeros(k,11);
for k=1:1:(calculate_T/dt)
rho=1.225*2.71828^(-1.389*10^(-4)*H);
k
t=k*dt;
g1=9.8-Vx^2/(6371000+H);
if m>=m0
if V<V_rocketstart%%用涡喷起飞爬升
Fbz=thrust_turbo;%恒定推力
m=m-Fbz/(isp_turbo*9.8)*dt;
mode=1;
end
if V>=V_rocketstart&&V<V_rbcc||H>30000%%大于涡喷最大速度，小于RBCC启动速度时，或高度大于30000，用火箭
Fbz=min(thrust_rocket,max_g_limit*9.8*m);%恒定推力
m=m-Fbz/(isp_rocket*9.8)*dt;
mode=3;
end
if V>=V_rbcc&&H<=30000
Fbz=thrust_rbcc;%恒定推力
m=m-Fbz/(isp_rbcc*9.8)*dt;
mode=2;
end
%%%%%%%%%%%%%%%%%%%%%%%%
else
Fbz=0;
end
%计算气动力
if Vy<0%%滑翔时增加Cl，Cd
Cl=glide_Cl;Cd=glide_Cd;
end
D=0.5*rho*V^2*S_wing*Cd;
%%%升力爬升曲线
if H<=push_down_height%%起飞阶段
if H<=0||L/m<=pull_up_g*9.8%地面滑跑
L=0.5*rho*V^2*S_wing*Cl;end
if L/m>pull_up_g*9.8&&theta<=max_theta/57.3 %%离地拉起升力不大于1.5g
L=min(9.8*pull_up_g*m,0.5*rho*V^2*S_wing*Cl);end
if theta>max_theta/57.3%%保持大角度火箭动力爬升
L=0;end
end
if H>push_down_height&&theta>push_down_angle/57.3&&H<pull_up_height%%推平
L=max(9.8*push_down_g*m,0.5*rho*V^2*S_wing*push_down_Cl);
end
if H>push_down_height&&theta<=push_down_angle/57.3&&H<pull_up_height%%保持小角度冲压爬升
L=min(9.8*pull_up_g*m,0.5*rho*V^2*S_wing*Cl);
end
if H>=pull_up_height&&theta<pull_up_angle/57.3%%重新大角度火箭动力爬升
L=min(9.8*pull_up_g*m,0.5*rho*V^2*S_wing*Cl);
end
if H>=pull_up_height&&theta>=pull_up_angle/57.3%%维持大角度火箭动力爬升
L=0;
end
%end
%%%
ax=((Fbz-D)*cos(theta)-L*sin(theta))/m;
ay=((Fbz-D)*sin(theta)+L*cos(theta))/m-g1;
acc_b=(Fbz-D)/m;
if H>=0&&Vx>=0
Vx=Vx+ax*dt;
Vy=Vy+ay*dt;
V=sqrt(Vx^2+Vy^2);
theta=asin(Vy/V);
X=X+Vx*dt;H=H+Vy*dt;
end
if H<0&&Vx>=0%%%%%%%%%%%%%%%%%%%%%地面滑跑
Vx=Vx+ax*dt;
Vy=0;
V=Vx;
theta=0;
X=X+Vx*dt;H=0;
end
if Vx<0
break
end
P(k,1)=X;
P(k,2)=H;
P(k,3)=L/(m*9.8);
P(k,4)=theta;
P(k,5)=V;
P(k,6)=Vx;
P(k,7)=Vy;
P(k,8)=D;
P(k,9)=t;
P(k,10)=ax;
P(k,11)=ay;
P(k,12)=mode;
end
range=P(:,1);height=P(:,2);az=P(:,3);theta=P(:,4);velocity=P(:,5);vx=P(:,6);vy=P(:,7);drag=P(:,8);time=P(:,9);ax=P(:,10);ay=P(:,11);
mode=P(:,12);
hold on
%plot(height,theta.*57.3)
%plot(range,height)
plot(subplot(2,1,1),height/1000,velocity/340,'r',height/1000,mode,'b')
title(subplot(2,1,1),'height vs Ma(Ma=V/340) ,mode')
xlabel(subplot(2,1,1),'height km ')
ylabel(subplot(2,1,1),'Mach or mode')
plot(subplot(2,1,2),time,velocity/1000,'r',time,height/10000,'b')
title(subplot(2,1,2),'time vs velocity and height')
xlabel(subplot(2,1,2),'time s ')
ylabel(subplot(2,1,2),'velocity km/s or height 10km')
%plot(time,velocity)
%plot(height,drag)
%plot(time,az);
%plot(time,sqrt(ax.^2+ay.^2))
% plot(height/1000,mode)
hold on
