rhosl= .002378;
rho15=.001496;
rho25=.001066;
rho35=.0007382;
W=520000;
S=4506;
b=200;
AR= b^2/S;
e=.86;
CDo=.0165;
Tsl=2*80000;
i=1;
for v=200:50:1400
CL= W*2/rho15/v^2/S;
k=1/(3.1416*e*AR);
Treq(i,1)= rho15/2*v^2*S*(CDo+ k*CL^2);
Tav(i,1)=rho15/rhosl*Tsl;
Preq(i,1)=Treq(i,1)*v;
Pav(i,1)=Tav(i,1);
V(i)=v;
EAS(i,1)=v/1.69*(rho15/rhosl)^.5;
i=i+1;
end
i=1;
for v=200:50:1400
CL= W*2/rho25/v^2/S;
k=1/(3.1416*e*AR);
Treq(i,2)= rho25/2*v^2*S*(CDo+ k*CL^2);
Tav(i,2)=rho25/rhosl*Tsl;
Preq(i,2)=Treq(i,2)*v;
Pav(i,2)=Tav(i,2);
V(i)=v;
EAS(i,2)=v/1.69*(rho25/rhosl)^.5;
i=i+1;
end
i=1;
for v=200:50:1400
CL= W*2/rho35/v^2/S;
k=1/(3.1416*e*AR);
Treq(i,3)= rho35/2*v^2*S*(CDo+ k*CL^2);
Tav(i,3)=rho35/rhosl*Tsl;
Preq(i,3)=Treq(i,2)*v;
Pav(i,3)=Tav(i,2);
V(i)=v;
EAS(i,3)=v/1.69*(rho35/rhosl)^.5;
i=i+1;
end
%h=figure;
%set(1,'PaperPosition' , [0.25 2.5 16 12])
%plot(V,Treq15,'o',V,Tav15,'o',V,Treq25,'+',V,Tav25,'+',V,Treq35,'s',V,Tav35,'s');
subplot(2,3,1);
plot(V(1:i-1),Preq(:,1),'ko-',V(1:i-1),Preq(:,2),'ks-',V(1:i-1),Preq(:,3),'kd-',V(1:i-1),pav(:,1),'ko:',V(1:i-1),Tav(:,2),'ks:',V(1:i-1),Tav(:,3),'kd:') ;
subplot(2,3,5);
plot(V(1:i-1),Treq(:,1),'ko-',V(1:i-1),Treq(:,2),'ks-',V(1:i-1),Treq(:,3),'kd-',V(1:i-1),Tav(:,1),'ko:',V(1:i-1),Tav(:,2),'ks:',V(1:i-1),Tav(:,3),'kd:') ;
legend('15000', '25000', '35000','15000', '25000', '35000');
Preq(:,1)=Treq(:,1).*V(1:i-1)';
Preq(:,1)=Treq(:,1).*V(1:i-1)';
Preq(:,1)=Treq(:,1).*V(1:i-1)';
Preq(:,2)=Treq(:,2).*V(1:i-1)';
Preq(:,2)=Treq(:,2).*V(1:i-1)';
Preq(:,2)=Treq(:,2).*V(1:i-1)';
Preq(:,3)=Treq(:,3).*V(1:i-1)';
Preq(:,3)=Treq(:,3).*V(1:i-1)';
Preq(:,3)=Treq(:,3).*V(1:i-1)';
Pav(:,1)=Tav(:,1).*V(1:i-1)';
Pav(:,1)=Tav(:,1).*V(1:i-1)';
Pav(:,1)=Tav(:,1).*V(1:i-1)';
Pav(:,2)=Tav(:,2).*V(1:i-1)';
Pav(:,2)=Tav(:,2).*V(1:i-1)';
Pav(:,2)=Tav(:,2).*V(1:i-1)';
Pav(:,3)=Tav(:,3).*V(1:i-1)';
Pav(:,3)=Tav(:,3).*V(1:i-1)';
Pav(:,3)=Tav(:,3).*V(1:i-1)';
xlabel('flight speed, fps');
ylabel('Thrust');
RC(:,1)=([Pav(:,1)]-[Preq(:,1)])*60/W;
RC(:,2)=([Pav(:,2)]-[Preq(:,2)])*60/W;
RC(:,3)=([Pav(:,3)]-[Preq(:,3)])*60/W;


%plot(V(1:i-1),RC(:,1:3)/60,'o-');
%plot(V(1:i-1),RC(:,1)/60,'ro-',RC(:,2)/60,'ro-',RC(:,3)/60,'ro-');
%figure
subplot(2,3,4)
plot(V(1:i-1),RC(:,1),'ko-',V(1:i-1),RC(:,2),'ks-',V(1:i-1),RC(:,3),'kd-');
xlabel('flight speed, fps');
ylabel('R/C, fpm');
legend('15000', '25000', '35000');
grid;
%figure;
%plot(EAS(1:i-1),RC(:,1:3)/60,'o-');
%figure;
subplot(2,3,2);
plot(EAS(1:i-1,1),RC(:,1),'ko-',EAS(1:i-1,2),RC(:,2),'ks-',EAS(1:i-1,3),RC(:,3),'kd-');
xlabel('EAS, knots');
ylabel('R/C, fpm');
legend('15000', '25000', '35000');

grid;
alt=[15000 25000 35000];
subplot(2,3,3);
plot(RC(5,1:3),alt(1:3),'ko-');
xlabel('R/C, fpm');
ylabel('Altitude, ft');
subplot(2,3,6);
T=[1/RC(5,1) 1/RC(5,2) 1/RC(5,3)];
plot(alt(1:3),T(1:3),'ks');

ylabel('1/RC');
xlabel('Altitude, ft');

RC(5,1:3)
grid
