var1=-.5+rand;  %standard proportional variation
for ii=1:length(z)
    tp(ii)=t(ii)+(var1)*td(ii);
    rp(ii)=r(ii)+(var1)*rd(ii);
    pp(ii)=p(ii)+(var1)*pd(ii);
    if pp(ii)<0 %remove negative pressure values 
        pp(ii)=0;
    else 
    end 
    if rp(ii)<0 %remove negative density values
        rp(ii)=0;
    else 
    end 
end

jitr=.25;    %ratio base variation 
base=1-jitr; %ratio variation every step 

for ii=1:length(z)
    tpv(ii)=t(ii)+(base*var1)*td(ii)+(jitr*(-.5+rand)*td(ii));
    rpv(ii)=r(ii)+(base*var1)*rd(ii)+(jitr*(-.5+rand)*rd(ii));
    ppv(ii)=p(ii)+(base*var1)*pd(ii)+(jitr*(-.5+rand)*pd(ii));
    if ppv(ii)<0
        ppv(ii)=0;
    else 
    end 
    if rpv(ii)<0
        rpv(ii)=0;
    else 
    end 
end

            if tim(step) >= totburnt && tim(step) < totburnt + dt
                if (flag < length(tburn) + 1)
                m_shot = m_shot - mstage(flag)*0.85;
                end
                flag = flag + 1;
                if (flag < length(tburn) + 1)
                totburnt = totburnt + tburn(flag);
                Fthrx = Fthr(flag)*cos(theta)*cos(phi);
                Fthry = Fthr(flag)*cos(theta)*sin(phi);
                Fthrz = Fthr(flag)*sin(theta);
                end
            end
            

            if tim(step) < totburnt
                m_shot = m_shot - ((0.15)/tburn(flag))*dt*mstage(flag);
            else
                Fthrx = 0;
                Fthry = 0;
                Fthrz = 0;
            end