You are an expert in writing code in many different languages. Your goal is to perform code completion for the following code, keeping in mind the rest of the code and the file meta data.
    Metadata details:
    filename: elastic_3D_FDTD_O22.m
    path of file: ./repos/WaveProp_in_MATLAB/elastic_3D_FDTD_wave_propagation
    
    the code of the file until where you have to start completion: % Finite-differences in time domain(FDTD) elastic wave propagation in 3D
% isotropic medium surrounded by simple sponge boundaries with exponential
% decay (Cerjan, 1985).
%
% We solve second order wave equation in time domain and displacement
% formulation getting wavefield in terms of displacement vector [ux, uy, uz].
%
% Elastic medium is parametrized by density and two velocities rho, vp
% and vs . We show CFL condition and number of points per wavelength
% prior running loop over time steps.
%
% We initiate source on a sphere with exponential decay.
%
% Conventional FD star-stencils deliver accuracy O(2,2)
% Stencils: [1 -2 1]/dx2 and [1 -1 -1 1]/4dxdz
% --------------------------------------------------------------
% The code is intentionally writen in a single file
% to simplify start up.
%
% The program does not save any files, add such option manually if needed.
% Drawing the wavefield is computationally demanding. Increase
% IT_DISPLAY value to reduce output and accelerate computation.
%
% The goal is to provide a simple example of elastic wave propagation
% in 3D isotropic medium.
%
% --------------------------------------------------------------
% Oleg Ovcharenko and Vladimir Kazei, 2018
%
% oleg.ovcharenko@kaust.edu.sa
% vladimir.kazei@kaust.edu.sa
%
% King Abdullah University of Science and Technology
% Thuwal, Saudi Arabia
% --------------------------------------------------------------

close all;
% Output ever
    
function (second derivative of a Gaussian):

force_x = sin(angle_force * deg2rad) * source_signal * dt2rho_src / (dx * dy * dz);
force_y = cos(angle_force * deg2rad) * source_signal * dt2rho_src / (dx * dy * dz);
force_z = sin(angle_force * deg2rad) * source_signal * dt2rho_src / (dx * dy * dz);
% Comment the line below if need 3 component force source
force_z = zeros(size(force_z));

% Moment tensor source signature is the same as of force 
% use moment tensor source = 1, use force source = 0
mt.flag = 0;

% Moment tensor components
mt.xx = 1;
mt.yy = 1;
mt.zz = 1;
mt.xy = 0;
mt.xz = 0;
mt.yz = 0;
mt.factor = factor;
mt.source_signal = source_signal;

%% OTHER
CFL = max(vp(:)) * dt / min(dx,dz);     % Courant number, should be < 1
min_wavelengh = min(vs(vs>0.1))/f0;     % shortest wavelength bounded by velocity in the air

%% DISTRIBUTED SOURCE OVER A SPHERE
% Radius of the spherical source [grid nodes]
szb = ceil(nx/(10 * 2));
% Create a 3d array with the ones shaping a sphere
SPb = strel('sphere',szb); sphere_b = double(SPb.Neighborhood);
% Uncomment the following 3 lines to make an empty sphere
% szs = szb - 3; dsz = szb - szs;
% SPs = strel('sphere',szs); sphere_s = double(SPs.Neighborhood);
% sphere_s  = padarray(sphere_s,[dsz dsz dsz]);
sphere_s = zeros(size(sphere_b));
sphere_e = sphere_b - sphere_s;

% Add zero padding so the source is of model size
sphere_e = padarray(sphere_e,[nz-ksrc-szb ny-jsrc-szb nx-isrc-szb]);

% Distance from given source location to each point of the distributed
% source sphere_e
dist = zeros(size(sphere_e));
% True point source location
src0 = [ksrc, jsrc, isrc];
for k = 1:size(sphere_e,1)
    for j = 1:size(sphere_e,2)
        for i = 1:size(sphere_e,3)
            point = [k,j,i];
            dist(k,j,i) = sqrt(sum((src0 - point).^2));
        end