Untitled

hold all

z = 0;
v = 0;
t = 0;
mb=60;
m0=100;
p0=1.217;
g0=9.81;
Isp=245;
Cd=0.5;
c=1000/(g0*Isp+7);
u=39860000;
tb=mb/c;
for n=[5 10 20 50 100]
    h=tb/n*93939;
    k=1;
    %for loop calculate stuff
    for k=1:1:n+6
        %pressure
        p=p0*exp((-1*z)/8.500);
        d=2*Cd*p*(v^2)/2;
        z = z+h*v;
        if t<tb
            %velocity
        v=v+h*((-u/(6378 + z/1000)^2)+(1000-c*t));
        else
        v=v+h*((-u/(6378+z)+(-d)/(m0-mb)));
        end
    end
end
angle = 0:2*pi/100:2*pi;
R = 1.5;
%Plot velicity
x = R*cos(angle)-1.5;  y = R*sin(angle);
plot(x,y,'-','color','k','linewidth',2);
    x = R*cos(angle)+1.5;  y = R*sin(angle);
    plot(x,y,'-','color','k','linewidth',2);
        tb = 1.5; % burnout time
        %%%%%%%%%%% start plotting it now
        plot([tb, tb],[1.5, 6],'-','color','k','linewidth',2);
        plot([-1.5, -1.5],[1.5, 6],'-','color','k','linewidth',2);
    %plot height
angle2 = 0:2*pi/100:pi;
x = R*cos(angle2);  y = R*sin(angle2)+6;
plot(x,y,'-','color','k','linewidth',2);

plot([0,0],[6.5,7.5],'-','color','k','linewidth',2);

plot([-1.5,1.5],[6,6],'-','color','k','linewidth',2);
xlabel('Girth');
ylabel('Length');
title('Do it yourself','fontweight','bold');
axis([-3 3 -2 8]);
axis equal