function dsq = brownian_displacement_simulation ( k, n, m, d, t )

%*****************************************************************************80
%
%% brownian_displacement_simulation() simulates Brownian displacement.
%
%  Discussion:
%
%    Thanks to Feifei Xu for pointing out a missing factor of 2 in the
%    stepsize calculation, 08 March 2016.
%
%    Thanks to Wang Zhe for pointing out a missing factor of m in the
%    stepsize calculation, 23 December 2018.
%
%    This function computes the square of the distance of the Brownian
%    particle from the starting point, repeating this calculation 
%    several times.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    26 December 2018
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer K, the number of repetitions.
%    The default is 20;
%
%    integer N, the number of time steps to take, plus 1.
%    This might be 1001.
%
%    integer M, the spatial dimension.  Typically, this is 2.
%
%    real D, the diffusion coefficient.  This might be 10.0.
%    Computationally, this is simply a scale factor between time and space.
%
%    real T, the total time, which defaults to 1.0.
%
%  Output:
%
%    real DSQ(K,N), the displacements over time for each repetition.
%    DSQ(:,1) is 0.0, because we include the displacement at the initial time. 
%

%
%  Supply default values if necessary.
%
  if ( nargin < 5 )
    t = 1.0;
  end

  if ( nargin < 4 )
    d = 10.0;
  end

  if ( nargin < 3 )
    m = 2;
  end

  if ( nargin < 2 )
    n = 1001;
  end

  if ( nargin < 1 )
    k = 20;
  end
%
%  Set the time step.
%
  dt = t / ( n - 1 );

  dsq = zeros ( k, n );

  for i = 1 : k
%
%  Compute the individual steps.
%
    x = zeros ( m, n + 1 );
%
%  Stepsize is normal.
%  Thanks to Wang Zhe for pointing out the missing "m" factor, 
%  25 December 2018.
%
    s = sqrt ( 2.0 * m * d * dt ) * randn ( 1, n - 1 );
%
%  Direction is random.
%
    if ( m == 1 )
      dx(1:m,1:n - 1) = s;
    else
      a = randn ( m, n - 1 );
      v = s ./ sqrt ( sum ( a.^2 ) );
      b = spdiags ( v', 0, n - 1, n - 1 );
      dx(1:m,1:n-1) = a * b;
    end
%
%  Each position is the sum of the previous steps.
%
    x(1:m,2:n) = cumsum ( dx(1:m,1:n-1), 2 );

    if ( m == 1 )
      dsq(i,1:n) = x(1,1:n).^2;
    else
      dsq(i,1:n) = sum ( x(1:m,1:n).^2 );
    end

  end

  return
end
