%
% $Id: roller.tex,v 0.0 2002/09/06 19:54:54 nicb Exp $
%
\documentclass[a4paper]{article}
%\includepackage{}
\usepackage{amssymb,amsmath,epsfig}
\setcounter{page}{1}
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\beqa}{\begin{eqnarray}}
\newcommand{\eeqa}{\end{eqnarray}}
\newcommand{\bit}{\begin{itemize}}
\newcommand{\eit}{\end{itemize}}
\newcommand{\bea}{\begin{array}}
\newcommand{\eea}{\end{array}}
\newcommand{\vsp}{\vspace{.5 cm}}
\newcommand{\hsp}{\hspace{.5 cm}}
\newcommand{\rowsp}{\rule[-.3cm]{0cm}{.6cm}}
\newcommand{\ds}{\displaystyle}
\newcommand{\neht}{\Leftarrow}
\newcommand{\then}{\Rightarrow}
\newcommand{\ioif}{\Leftrightarrow}
\newcommand{\mbf}[1]{\mbox{\boldmath ${#1}$}}
\newcommand{\R}{{\sf R\hspace*{-0.9ex}\rule{0.15ex}{1.5ex}\hspace*{0.9ex}}}
\newcommand{\N}{{\sf N\hspace*{-1.0ex}\rule{0.15ex}{1.3ex}\hspace*{1.0ex}}}
\newcommand{\Q}{{\sf Q\hspace*{-1.1ex}\rule{0.15ex}{1.5ex}\hspace*{1.1ex}}}
\newcommand{\C}{{\sf C\hspace*{-0.9ex}\rule{0.15ex}{1.3ex}\hspace*{0.9ex}}}
\newcommand{\Comment}[1]{}
\newcommand{\define}{\stackrel{\triangle}{=}}


\title{Rolling a ball along an arbitrarily shaped slope}
\author{Davide Rocchesso - Zueca, Venezia
}



\begin{document}
\maketitle
\begin{abstract}
I simulate a ball rolling down a slope with arbitrary shape
\end{abstract}

\section{Equations}

Call $p$ the position of the ball along the curve, and $v$ its velocity. Let $\alpha(x)$ be the angle formed by the slope at point $x$ with the horizontal line. We have 
\beqa
\dot{p}  & = & v; \nonumber \\
\dot{v}  & = & g \sin{\alpha(x)} \define \sigma(p). \label{dinamica}
\eeqa
Assuming that $p, v, \sigma$, as functions of time, are Laplace transformable, we can easily express the equations ~(\ref{dinamica}) as algebraic equalities in the variable $s$. Then, we can apply the bilinear transformation and come up with the following discrete-time equations:
\beqa
p(n)  & = & p(n-1) + \frac{1}{h} \left( v(n) + v(n-1)\right); \nonumber \\
v(n)  & = & v(n-1) + \frac{1}{h} \left( \sigma(n) + \sigma(n-1)\right). \label{dinamicadiscreta}
\eeqa
The equations~(\ref{dinamicadiscreta}) can be equivalently obtained by 1-step Adams-Moulton discretization of~(\ref{dinamica}).

The problem here is that there is an instantaneous dependency between $\sigma$ and $p$.

Defined ${\bf w} = \left[p \; v \right]'$, we can rewrite~(\ref{dinamicadiscreta}) as
\be
{\bf w}(n) = \left({\bf A}^{-1}\right)^2 {\bf w}(n-1) + {{\bf A}^{-1}} {\bf B} \sigma(n-1) + {{\bf A}^{-1}} {\bf C} \sigma(n) \; ,
\ee
where ${\bf B} = {\bf C} = \left[0 \dots 1/h \right]'$, and 
\be
{\bf A} = \left[ \bea{ll} 
1 & -1/h \\
0 & 1 
\eea \right] \; .
\ee

The $n$-th value of the variable $\sigma$ can be computed from
\beqa
\lefteqn{\sigma(n) =  f \left(p(n) \right) = \nonumber} \\ 
& = & f \left( p(n-1) + 2/h v(n-1) + 1/h^2 \sigma(n-1) \right. \nonumber \\ & & \hsp \left. + 1/h^2 \sigma(n) \right) \\ 
& \define & f \left(\tilde{p}(n) + 1/h^2 \sigma(n) \right) \nonumber \; ,
\eeqa
where the computable and known terms have been separated from the instantaneous dependency.

Instead of using the nonlinear function $f$, we use the sheared version 
\be
\sigma(n) = f_1\left( \tilde{p}(n) \right) \; .
\ee

\section{Equations with viscous damping}
\label{eqsdamp}

Just to be pedantic, I repeat the equations when the rolling ball of mass $m$ is affected by viscous damping with coefficient $g$.
Call $p$ the position of the ball along the curve, and $v$ its velocity. Let $\alpha(x)$ be the angle formed by the slope at point $x$ with the horizontal line. We have 
\beqa
\dot{p}  & = & v; \nonumber \\
\dot{v}  & = & g \sin{\alpha(x)} - \gamma v \define \sigma(p) - \gamma v , \label{dinamicag}
\eeqa
where $\gamma \define g / m$ and it has dimensions of $1/\mbox{s}$.
Assuming that $p, v, \sigma$, as functions of time, are Laplace transformable, we can easily express the equations ~(\ref{dinamicag}) as algebraic equalities in the complex variable $s$. Then, we can apply the bilinear transformation and come up with the following discrete-time equations:
\beqa
p(n)  & = & p(n-1) + \frac{1}{h} \left( v(n) + v(n-1)\right); \nonumber \\
v(n)  & = & \frac{h}{h + \gamma}v(n-1) + \frac{1}{h + \gamma} \left( \sigma(n) + \sigma(n-1)\right). \label{dinamicadiscretag}
\eeqa
The equations~(\ref{dinamicadiscretag}) can be equivalently obtained by 1-step Adams-Moulton discretization of~(\ref{dinamicag}).

The problem here is that there is an instantaneous dependency between $\sigma$ and $p$.

Defined ${\bf w} = \left[p \; v \right]'$, we can rewrite~(\ref{dinamicadiscretag}) as
\be
{\bf w}(n) = {\bf A}^{-1}{\bf A}_h {\bf w}(n-1) + {{\bf A}^{-1}} {\bf B} \sigma(n-1) + {{\bf A}^{-1}} {\bf C} \sigma(n) \; ,
\ee
where ${\bf B} = {\bf C} = \left[0 \; 1/(h + \gamma) \right]'$,  
\be
{\bf A} = \left[ \bea{ll} 
1 & -1/h \\
0 & 1 
\eea \right] \; ,
\ee
and
\be
{\bf A}_h = \left[ \bea{ll} 
1 & -1/h \\
0 & h/(h + \gamma) 
\eea \right] \; ,
\ee

The $n$-th value of the variable $\sigma$ can be computed from
\beqa
\lefteqn{\sigma(n) =  f \left(p(n) \right) = \nonumber} \\ 
& = & f \left( p(n-1) + \frac{2 h + \gamma}{h \left( h + \gamma \right)} v(n-1) + 1/(h (h + \gamma)) \sigma(n-1) \right. \nonumber \\ & & \hsp \left. + 1/(h (h + \gamma)) \sigma(n) \right) \\ 
& \define & f \left(\tilde{p}(n) + 1/(h (h + \gamma)) \sigma(n) \right) \nonumber \; ,
\eeqa
where the computable and known terms have been separated from the instantaneous dependency.

Instead of using the nonlinear function $f$, we use the sheared version 
\be
\sigma(n) = f_1\left( \tilde{p}(n) \right) \; .
\ee


\section{Octave script}
The following octave script implements the simulation of the system presented in section~\ref{eqsdamp}. Here the curve is computed and tabulated in parametric form ({\tt t} is the parameter), and $\sigma$ is also computed for each value of the parameter. The sheared nonlinear function is obtained by using $\tilde{p}$ as an  argument, for each value of the parameter {\tt t}.

{\tt
\begin{verbatim}
Fs = 100; % in samples per second: try a very small Fs to see the difference between p and ptilde
          % A very small Fs may produce ringing near equilibrium (non-passivity)
h = 2 * Fs; % constant for bilinear transformation (i.e., Adams-Moulton 1 step integration)
simulationlength = 100; % in seconds
inct = 0.01; % linear increment to express the curve in parametric form
          % too few points may produce ringing near equilibrium (non-passivity)
tmax = 0.5; % maximum value of the parameter (for parametric curve), in meters
tdipmax = 0.05; % half width of the bottom tip
gravity = 9.8; % vertical acceleration in m/s
gamma = 0.6; % damping over mass [1/seconds]
t = [-tmax:inct:tmax];

% Parametric curve: (x(t), y(t))
x = t;
y = log(abs(t.*(abs(t)>=tdipmax))+1);
y = y + (y==0).*log(tdipmax+1); % curve with a flat tip
y=y.^2;
figure(1); xlabel('m'); ylabel('m'); plot(x,y,';y(x);'); 

% Array of spatial increments between adjacent point of the parametric curve
pinc = sqrt((x(2:length(x)) - x(1:length(x)-1)).^2 + (y(2:length(y)) - y(1:length(y)-1)).^2);

% pslope = dy/dx (local slope at each point computed by finite difference)
pslope = zeros(1, length(pinc)+1);
pslope(2:length(pslope)) = diff(y) ./ diff(x);

% p = linear coordinate along the curve (for each parameter value)
p = zeros(1, length(pinc)+1);
for i=1:length(pinc)
  p(i+1) = p(i) + pinc(i);
end;
pslope(1) = pslope(2);

% Precalculation of sigma = f(p)
sigma = -gravity*sin(atan(pslope));

% Sheared variable ptilde, to access the sheared function sigma = f1(ptilde)
ptilde = p - 1/(h * (h + gamma)) * sigma;

figure(2);
xlabel('m'); ylabel('m/s^2'); plot(p,sigma,'r;sigma=f(p);'); hold on;
plot(ptilde, sigma,'g;sigma=f1(ptilde);'); hold off;

figure(3); 
xlabel('m'); ylabel('m');
plot(t,p,';p(t);'); hold on; plot(t,ptilde,';ptilde(t);'); hold off;

instants = simulationlength * Fs;
pn = zeros(1,instants); vn = zeros(1,instants); sigman = zeros(1,instants);
ptilden = zeros(1,instants); 

% Simulation cycles
for n = 2:instants
  ptilden(n) = pn(n-1) + (2*h + gamma)/(h * (h + gamma)) * vn(n-1) + 1/(h * (h + gamma)) * sigman(n-1);

  % find the index of vector ptilde corresponding to value ptilden(n)
  [minimo, indminimo] = min(abs(ptilde - ptilden(n)));

  % sigman(n) = f1(ptilden(n))
  sigman(n) = sigma(indminimo);

  % update tangential velocity
  vn(n) = h / (h + gamma) * vn(n-1) + 1/(h + gamma) * (sigman(n) + sigman(n-1));

  % update position along the curve
  pn(n) = pn(n-1) + 1/h * (vn(n) + vn(n-1));

%  printf("ptilden=%f, indminimo=%d, sigman=%f, vn=%f, pn=%f\n", ptilden(n), indminimo, sigman(n), vn(n), pn(n));
%  pause;

end;

figure(4);
ylabel('m'); xlabel('s'); 
plot([1:instants]/Fs, pn, ';curvilinear position in time;'); 
hold on; 
ylabel('m/s'); xlabel('s');
plot([1:instants]/Fs, vn, ';curvilinear velocity in time;'); 
hold off;
\end{verbatim}
}

\section{Application}
The first release of the invisiball will be one-dimensional. Consider a user pushing the membrane at a certain place. The depth of the push can be detected and, from this, a deflection profile can be deduced (available in closed form from Morse-Ingard). Given the initial position of the ball, we can compute its distance from the pushing point, so that we can locate the ball along the curve and start the simulation that will move this phantom ball along the straight line connecting the initial ball position to the pushing point.

The octave script seems to give reasonable results. In the lossless case ($\gamma = 0$), even at ridicolous sample rates ($1 Hz$ or so), hundreds of oscillations can be produced without blowing up or down the system. An artifact appears in the octave prototype at very low sample rates or for very small tables for {\tt ptilde}: the ball keeps oscillating around the bottom of the bowl. I thought this could be due to the line
{\tt
\begin{verbatim}
  % find the index of vector ptilde corresponding to value ptilden(n)
  [minimo, indminimo] = min(abs(ptilde - ptilden(n)));
\end{verbatim}
}
used to look for the element having value nearest to {\tt ptilden(n)} in {\tt ptilde}. 
Indeed, it would be safer to look for the nearest value in a ``passive'' sense, i.e. approximating always with the element closest to the middle point of the table.
In the C translation of the code, we implemented this search function as a ``passive'' binary search, using the assumption that the table will be monotonically increasing. However, the artifact seems to persist for very small sample rates (e.g., 10 Hz).

In the final application, a model such as this will control the velocity of Matthias' rolling ball algorithm. This velocity control module can be run at a very slow rate, if needed.

 \end{document}
