function tether_balls ( edrag, ecol, eg )

%*****************************************************************************80
%
%% tether_balls() simulates the motion of tethered balls in 2D.
%
%  Discussion:
%
%    Verlet integration is used to simulate the change in position 
%    and velocity.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    09 September 2022
%
%  Author:
%
%    Algorithm by Nick Berry.
%    This version by John Burkardt.
%
%  Reference:
%
%    Nick Berry,
%    Verlet integration,
%    July 2018,
%    https://datagenetics.com/blog/july22018/index.html
%
%  Input:
%
%    real edrag, should be slightly less than 1 to simulate drag.
%
%    real ecol, should be slightly less than 1 to simulate energy
%    loss in collisions.
%
%    real eg, the gravitational acceleration.
%
  if ( nargin < 1 )
    edrag = 0.999;
  end

  if ( nargin < 2 )
    ecol = 0.90;
  end

  if ( nargin < 3 )
    eg = 0.0;
  end

  blob_r = 5.0;
  ct = 0;
  gravity = false;
  height = 490.0;
  mode = 0;

  n = ?

  while ( true )

    [ px, py, px_old, py_old ] = update_points ( px, py, px_old, ...
      py_old, ecol, edrag, eg, width, height, blob_r );

    render_points ( );

  end

  return
end
function [ px, py, px_old, py_old ] = update_points ( px, py, px_old, ...
  py_old, ecol, edrag, eg, width, height, blob_r )

%*****************************************************************************80
%
%% update_points()...
%
  n = size ( px, 1 );

  for i = 1 : n
    vx = ( px(i) - px_old(i) ) * edrag;
    vy = ( py(i) - py_old(i) ) * edrag;
    px_old(i) = px(i);
    py_old(i) = py(i);
    px(i) = px(i) + vx;
    py(i) = py(i) + vy + eg;

    if ( ( width - blob_r) < px(i) ) 
      px(i) = width - blob_r;
      px_old(i) = px(i) + vx * ecol;
    elseif ( px(i) < blob_r )
      px(i) = blob_r;
      px_old(i) = px(i) + vx * ecol;
    end
  
    if ( height - blob_r < py(i) ) 
      py(i) = height - blob_r;
      py_old(i) = py(i) + vy * ecol;
    elseif ( py(i) < blob_r )
      py(i) = blob_r;
      py_old(i) = py(i) + vy * ecol;
    end

  end

  return
end

