為方便各位KSP愛好者在RSS環(huán)境下手搓SSTO，公開電腦里庫存的自己練手寫的matlab小程序，可以擬定各種參數(shù)，粗略估計(jì)一架SSTO能否入軌。如果在這里面都算出來不到第一宇宙速度，那在游戲里大概率是入不了軌了。

例如一架250t起飛重量、燃料占比85%的SSTO，用渦噴起飛，6°小角度爬升~10000m，隨后推平到3.5°，并加速到~2Ma。假設(shè)渦噴比較屑，到~2Ma就無法繼續(xù)加速，便啟動(dòng)火箭加速至~3Ma填補(bǔ)推力陷阱，直到RBCC啟動(dòng)速度，然后繼續(xù)3.5°小角度爬升~30000m，加速到6~7Ma。3萬米以上缺氧，所以拉起20°并切換到火箭模式入軌。

為保護(hù)小綠人不受摧殘，限制z向過載3g，軸向過載3.5g。

本程序簡(jiǎn)化了升阻力系數(shù)（認(rèn)為不隨速度變化）、發(fā)動(dòng)機(jī)比沖（也不隨速度變化），所以僅供參考，不對(duì)小綠人的安全負(fù)責(zé)。

上升軌道如圖：

代碼如下

clear all

%%用渦噴起飛，小角度爬升~10000m，隨后加速到~2Ma，啟動(dòng)火箭加速至~3Ma填補(bǔ)推力陷阱，

%隨后RBCC啟動(dòng)，繼續(xù)小角度爬升~30000m，加速到6~7Ma，拉起切換到火箭模式入軌。

%mode1是渦噴，mode2是RBCC，mode3是火箭

%基本參數(shù)

max_g_limit=3.5;%大于該過載，自動(dòng)下降推力

thrust_turbo=1200000;%推力牛頓

isp_turbo=2500;%比沖s

thrust_rocket=3000000;

isp_rocket=455;

thrust_rbcc=2500000;

isp_rbcc=1000;

V_rocketstart=600;%t填補(bǔ)推力缺口，大于該速度后啟動(dòng)火箭加速至rbcc啟動(dòng)速度然后切換rbcc

V_rbcc=900;%rbcc啟動(dòng)的速度

m=250000;%起飛重量kg

m0=m*0.15;%空重

%aerodynamic

Cd=0.04;%爬升階段阻力系數(shù)

Cl=0.2;%爬升最大升力系數(shù)

glide_Cl=1;%滑翔升力系數(shù)，如果不小心再入了，用得上這個(gè)值

glide_Cd=1;%滑翔阻力系數(shù)，如果不小心再入了，用得上這個(gè)值

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;%發(fā)射架初始位置，高度m，水平距離m，默認(rèn)都是0

%air

rho=1.225*2.71828^(-1.389*10^(-4)*H);%空氣密度

D=0;L=0;%升阻力初始化設(shè)置0，勿動(dòng)

%Launch setting

dt=0.1;%時(shí)間步長

calculate_T=900;%計(jì)算時(shí)間總長度

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啟動(dòng)速度時(shí)，或高度大于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

?%計(jì)算氣動(dòng)力

?if Vy<0%%滑翔時(shí)增加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%%保持大角度火箭動(dòng)力爬升

? ? ? ? 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%%重新大角度火箭動(dòng)力爬升

? ? ? ? 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%%維持大角度火箭動(dòng)力爬升

? ? ? ? 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