Error with if loop within a while loop.

h(1)=150000;            %initial height
a(1)=(40*10^7)/(6371+h(1))^2;        %initial acceleration dependant on height
dt=0.005;        %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
m=850;   %Mass
c=0.7
p=(h/71)+1.4;        %Air Density
Fd=0.5*p*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)=150000-(g(i)*(t(i+1))^2);
    if h>100000
        g(i+1)=(40*10^7)/(6371+h(i+1))^2;
        a(i+1)=g(i+1);
    else
        p(i+1)=(h(i+1)/71)+1.4
        Fd(i+1)=0.5*(p(i+1))*c*As*(v(i+1))^2;
        a(i+1)=(g(i+1))-((Fd(i+1))/m);
    end
    v=(a(i+1))*(t(i+1));
    i=i+1;
end

Befor I added the if loop I could plot graphs for (t,v)and (t,a) whilst h>100000. I have added the if loop so i can see data while h>0, but an error appears when I try to run it. Any help would be appricieted.

[Geoff Hayes]

Sam - what is the error message? Please copy and paste the full message to this question. Although, when I run your code I don't see any errors...it just continues for 30000+ iterations (which is perhaps the error?). Also

if h>100000

may not make sense since h is an array. Do you mean

if h(i+1) > 100000

instead?

[Sam Potter]

Thank you for the reply. The error that pops up after running is this:Index exceeds the number of array elements (1). However, this error only pops up mid way through the loop. When i run it, the loop counter is the same as it was for when I was only working out h(end)>=100000 before I added the if loop, which I think means the erro is somewhere in the if loop? I also changed it to if h(i+1) but the erro still persisted. Do you think it is becuase i dont have an initial value for p?.

Answer 1

[Ameer Hamza]

As pointed out by Geoff, that you should correct the line like this

if h(i+1) > 100000

However, this was not causing the error. The error was caused by the line

Fd(i+1)=0.5*(p(i+1))*c*As*(v(i+1))^2;

Here you are using v(i+1), but in your loop, you are not saving values of v.

Following code will fix this issue, also it appears that you missed the calculation of g inside else, I also added it

clear
h(1)=150000;            %initial height
a(1)=(40*10^7)/(6371+h(1))^2;        %initial acceleration dependant on height
dt=0.005;        %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
m=850;   %Mass
c=0.7;
p=(h/71)+1.4;        %Air Density
Fd=0.5*p*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)=150000-(g(i)*(t(i+1))^2);
    if h(end)>100000
        g(i+1)=(40*10^7)/(6371+h(i+1))^2;
        a(i+1)=g(i+1);
    else
        p(i+1)=(h(i+1)/71)+1.4;
        Fd(i+1)=0.5*(p(i+1))*c*As*(v(i))^2;
        g(i+1)=(40*10^7)/(6371+h(i+1))^2;
        a(i+1)=(g(i+1))-((Fd(i+1))/m);
    end
    v(i+1)=(a(i+1))*(t(i+1));
    i=i+1;
end