MATLAB script for TREL's drift model: Made minor changes

drift_sim.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Project           : PARACHUTE DRIFT MODEL
%
% Program name      : drift_sim.m
%
% Author            : John Steinman, Kadhir Umasankar (RE)
%
% Last revised      : 20191209
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc; clear; close all;

% DEFINE INTITIAL CONDITIONS ==================================
v_d = 314.286;  % (m/s) terminal velocity of main chute
v_m = 7.62;     % (m/s) terminal velocity of main chute
h_d = 10000;    % (m) height of drogue deployment
h_m = 5000;     % (m) height of main deployment
nsim = 2000;    % number of simulations
dt = 1;         % sec

% SIMULATE DRIFT ===================================================
drift = zeros(1, nsim);
x_lim = [0 0];  % storing the maximum extents of the graph in the x-axis

for i= 1:nsim
  hold on
  [x_d, y_d] = driftSim(h_d, h_m, 0, v_d, dt); % finding drift when just the drogue parachute is deployed
  [x_m, y_m] = driftSim(h_m, 0, x_d(end), v_m, dt); % finding drift when the main parachute is deployed
  drift(i) = x_m(end);
  x = [x_d x_m]; y = [y_d y_m]; % concatenating drogue drift and main drift
  x_lim = find_x_limits(x, x_lim); % using function to find xlim
  p1 = plot(x,y,'b-', 'Linewidth', .3, 'HandleVisibility','off'); % HandleVisibility keeps it from showing up in the legend
  p1.Color(4) = 0.1;
end
axis([x_lim(1) x_lim(2) 0 h_d+50])
xlabel('Drift (m)')
ylabel('Altitude (m)')
sigma = std(drift);
mean = mean(drift);
% plotting confidence intervals
line([-3*sigma+mean, 3*sigma+mean],[0,0],'Color','r','LineWidth',3)
line([-2*sigma+mean, 2*sigma+mean],[0,0],'Color','y','LineWidth',3)
line([-sigma+mean, sigma+mean],[0,0],'Color','g','LineWidth',3)
legend('99.7% confidence','95% confidence','68% confidence')

% FUNCTIONS =====================================================

function[x, y] = driftSim(h1, h2, x_0, v, dt)
    %% Function Name: driftSim
    %
    % Purpose:
    %   Find random paths that the rocket can take during downward
    %   trajectory
    %
    % Assumptions: 
    %   X-velocity of rocket goes to 0 when each parachute is
    %   deployed
    %   Y-velocity stays at terminal velocity throughout flight
    %
    % Inputs:
    %   h1: initial height
    %   h2: final height
    %   x_0: initial x-position
    %   v: vertical velocity
    %   dt: time differential
    %
    % Outputs:
    %   x: vector of x-values for this part of flight
    %   y: vector of y-values for this part of flight
    
    m   = 217.724; % (kg) weight of rocket
    A_n = 8.2296 * 0.381; % (m^2) area of rocket normal to wind
    vx = 0; % assuming the x-velocity ofrocket goes to 0 when each parachute is deployed
    h = h1;
    i = 1;
    y(1) = h1; x(1) = x_0;
    while(h > h2)
        i = i + 1;

       % Change in x pos
       wind_v = 2*(60/10000*h)*rand-(60/10000*h); % uniform

       % p = interp1(rho_alt, rho, h);
       % interpolating adds computational cost

       % accounting for changes in wind pressure with altitude
       if (h >= 9000 && h <= 10000), p = 0.4671;
       elseif (h >= 8000), p = 0.5258;
       elseif (h >= 7000), p = 0.59;
       elseif (h >= 6000), p = 0.6601;
       elseif (h >= 5000), p = 0.7364;
       elseif (h >= 4000), p = 0.8194;
       elseif (h >= 3000), p = 0.9093;
       elseif (h >= 1000), p = 1.112;
       elseif (h >= 2000), p = 1.007;
       elseif (h >= 0), p = 1.225;
       end
        
       % finding force due to wind, corresponding changes in
       % horizontal velocity and horizontal displacement, and new
       % x-velocity and x-position
       F_w = 0.5 * p * wind_v*abs(wind_v) * A_n;
       dvx = F_w/m * dt;
       dx = vx*dt + 0.5*F_w/m*dt^2;
       vx = vx + dvx; % find new velocity
       x(i) = x(i-1) + dx;

       % Change in altitude
       dh = abs(v)*dt;
       h = h - dh;
       y(i) = h;
    end
end

function x_lim = find_x_limits(x, x_lim)
    %% Function Name: find_x_limits
    %
    % Purpose: 
    %   To find the x limits of the graph. It was made into a function to declutter the code and keep it from detracting from the important part of the script
    %
    % Inputs:
    %   x: vector of x-values for that path
    %   x_lim: x-limits
    %
    % Outputs:
    %   x_lim: new x-limits

    if min(x) < x_lim(1)
        x_lim(1) = min(x);
    end
    if max(x) > x_lim(2)
        x_lim(2) = max(x);
    end
end