Issue with if loop within kinematics while loop

Hi, apologies i you have already seen this code, I am havng issues with this step of it. For context it is a spacesip falling from 150000 metres. At 3000 metres a parachute opens meaning that the area changes from As (spaceship=5) to Ap(parachute = 150)

h(1)=150000;            %initial height
a(1)=(40*10^7)/(6371+h(1))^2;        %initial acceleration dependant on height
dt=5;        %time step
t(1)=0;         %initial time
v(1)=a(1)*t(1);         %velocity
g(1)=((40*10^7)/(6371+h(1))^2);     %Downward force h>100000, upward force = 0 above 100000
As= 5;      %Area of spaceship
Ap = 150;  %area of parachute
m=850;   %Mass
c=0.7;
p(1)=-(100/71)+1.4;        % Initial Air Density (Air density occurs at h=100000, from then p=(h/71)+1.4)
Fd(1)=0.5*p(1)*c*As*v^2;     %Downward force h<=100000
i=1; %loop counter
while h(end)>=0
    t(i+1)=t(i)+dt;
    h(i+1)=h(i)-(v(i)*dt); % Find the height of previous time increment
    g(i+1)=(40*10^7)/((6371+h(i+1))^2);
    if h(i+1)>100000
        a(i+1)= g(i+1)  %Acceleration=Gravity-(Fd/m)
        v(i+1)=v(i)+(a(i+1)*dt);
    elseif h(i+1)>=3000
        p(i+1)=-((h(i+1)/1000)/71)+1.4; 
        Fd(i+1)=0.5*c*(p(i+1))*As*(v(i))^2;
        a(i+1)=g(i+1)-(Fd(i+1)/m);   %Acceleration=Gravity-(Fd/m)
        v(i+1)=v(i)+(a(i+1)*dt);
    else
        p(i+1)=-((h(i+1)/1000)/71)+1.4; 
        Fd(i+1)=0.5*c*(p(i+1))*Ap*(v(i))^2;  %Air resistanace open parachute
        a(i+1)=g(i+1)-(Fd(i+1)/m);   %Acceleration=Gravity-(Fd/m)
        v(i+1)=v(i)+(a(i+1)*dt);
    end
    i=i+1;
end

It runs fine until the parachute opens. The acceration becomes massively negative, causinng the spaceship to go miles into the atmosphere. What should happen is a drop in acceration to a reasonable negative acceleration, a small but positive velocity, and the height to lowly drop to 0. The difference between the elseif and the else loop if that in the Fd statement I have changed As(=5) to Ap(=150). Any help would be apprecieted.

[darova]

[Sam Potter]

Hi, Is this better?

[Cris LaPierre]

I've reformatted the code in your post above as code for you using the button darova pointed out.

It might be more helpful if you stick to one post rather than creating new posts for each step of your problem. It looks like you are struggling with if-statements and loops. Have you had a look at Chapter 11 of MATLAB Onramp, which explains how they work?

[Cris LaPierre]

A quick look suggest units are still an issue! Your equation for a (and g) expects h to be in km, not meters.

Answer 1

[darova]

Everything is OK. The problem is in the parachute:

when h < 3000, Fd is large (because of parachute of course). Fd → a → v (large velocity)

Next step: Fd=Fd(v) is function of v^2. So Fd is incredibly large

Solution: change time step after parachute opens

cla
hold on
while h(end)>3000
    % your code
    i=i+1;
    plot(h(end),a(end),'.r')
    pause(0.01)
end
dt = 0.1;
i = i - 1;
while h(end) > 0
    % code if h < 3000
    i = i + 1;
    plot(h(end),a(end),'.b')
    pause(0.01)
end
hold off