PDA

View Full Version : Equations of Motion (Physics of Spin/Roll)

tailuge
04-25-2006, 06:48 PM
Hi,

I would really like to have an equation that can
predict the position of a ball at time t based
on an initial velocity vector and rotation vector.

pos(t)= func(t,V,W);
rot(t)= func(t,V,W);

I have read an article by Shepard which comes close
to telling me what I need, but I could not quite
get it working when writing a program. What I need
is an explicit equation for the force acting on a
spining sliding ball.

I have also had experience of solving the time at which
two balls travelling on parabolas will impact by using
the quartic equation - and discovered that it reveals
limitations of floating point representations in
C++ code , has anyone else encountered such issues?

There are many great resources online for billiards
but I would also like to place the equations of motion
for others to use.

T,

Three cushion billiards fan.

pigbrain
04-25-2006, 08:29 PM
i am so interested in such a program.
i wonder if i can have a copy. : )

Jal
04-26-2006, 02:14 AM
For a ball with velocity V and spin W (radians/sec), and where R is the vector from the center of the ball to the contact point with the surface of the table, the surface speed of the ball at the contact point is:

Vs = V + W X R

where W X R denotes the vector cross product.

In terms of components:

Vsx = -RWy + Vx
Vsy = +Rwx + Vy

From this we get its linear acceleration components:

Ax = -(Vsx/Vs)ug
Ay = -(Vsy/Vs)ug

where

Vs = +Sqrt(Vsx^2 + Vsy^2),
u = coefficient of friction, approx. = .2 or .22 or whatever,
g = gravitational acceleration.

Its angular acceleration components are:

Agx = +(5/2)(1/R)Ay
Agy = -(5/2)(1/R)Ax

So if it starts out with speed V0 (V0x and V0y), and spin W0 (W0x and W0y), at any later time t up until natural roll sets in, its velocity components will be:

Vx = V0x + (Ax)t
Vy = V0y + (Ay)t

Its spin components will be:

Wx = W0x + (Agx)t
Wy = W0y + (Agy)t

And its position will be:

X = X0 + (V0x)t + (1/2)(Ax)t^2
Y = Y0 + (V0y)t + (1/2)(Ay)t^2

Note: We're assuming a right-handed coordinate system. R is alway positive in the above even though it may be viewed as pointing in the negative z-direction. The signs on the terms take this into account.

An implicit assumption is that the direction of the friction force remains the same throughout. This may or may not be obvious (it isn't to me), but it can be shown that it does. We're also assuming that the coefficient of friction remains the same at all surface speeds, and that the ball is not bouncing. Also, we're ignoring the z-component of its spin which will cause a very slight masse' action due to its interaction with the cloth (which it has to climb over, sort of), but this is very tiny. If you don't ignore all of these things, you would have to do a timewise numeric integration I guess.

If the ball reaches natural roll before colliding with something, its speed components just as it reaches roll will be:

Vx = (5/7)[V0x + (2/5)(R)W0y]
Vy = (5/7)[V0y - (2/5)(R)W0x]

Dr. Dave derives similar equations in one of his technical proofs "TPA.4 Post-impact cue ball trajectory..." if you want to double check. (At least I hope they're the same but I haven't checked them against these in every detail yet.) See:

http://www.engr.colostate.edu/~dga/pool/

Jim

DickLeonard
04-26-2006, 06:53 AM
Jal now I know I was playing the wrong game.####

Jal
04-26-2006, 10:01 AM
<blockquote><font class="small">Quote DickLeonard:</font><hr> Jal now I know I was playing the wrong game.#### <hr /></blockquote>But then we wouldn't have your great stories, amongst other things. Besides, after a judicious application of the above, I still havent the foggiest where the cueball will decide to go a good part of the time.

Jim

wolfdancer
04-26-2006, 01:56 PM
Dick, did you know that willie Mays had a Doctorate in Ballistics?....that's how he was able to make that overhead catch against Vic Wertz

My dog has a masters in aerodynamics...that's how he can catch a Frisbee.
I always though Mays throw was at least as amazing as the catch.

dr_dave
04-26-2006, 02:22 PM
tailuge,

You can find all of the details and derivations in TP A.4 (http://www.engr.colostate.edu/~dga/pool/technical_proofs/new/TP_A-4.pdf).

Regards,
Dave

<blockquote><font class="small">Quote tailuge:</font><hr> Hi,

I would really like to have an equation that can
predict the position of a ball at time t based
on an initial velocity vector and rotation vector.

pos(t)= func(t,V,W);
rot(t)= func(t,V,W);

I have read an article by Shepard which comes close
to telling me what I need, but I could not quite
get it working when writing a program. What I need
is an explicit equation for the force acting on a
spining sliding ball.

I have also had experience of solving the time at which
two balls travelling on parabolas will impact by using
the quartic equation - and discovered that it reveals
limitations of floating point representations in
C++ code , has anyone else encountered such issues?

There are many great resources online for billiards
but I would also like to place the equations of motion
for others to use.

T,

Three cushion billiards fan.

<hr /></blockquote>

dr_dave
04-26-2006, 02:24 PM
tailuge,

FYI, if you want to check out other pool physics resources available online and in print, check out the Physics Resources section of my website (http://www.engr.colostate.edu/~dga/pool/physics/index.html).

Regards,
Dave

<blockquote><font class="small">Quote tailuge:</font><hr> Hi,

I would really like to have an equation that can
predict the position of a ball at time t based
on an initial velocity vector and rotation vector.

pos(t)= func(t,V,W);
rot(t)= func(t,V,W);

I have read an article by Shepard which comes close
to telling me what I need, but I could not quite
get it working when writing a program. What I need
is an explicit equation for the force acting on a
spining sliding ball.

I have also had experience of solving the time at which
two balls travelling on parabolas will impact by using
the quartic equation - and discovered that it reveals
limitations of floating point representations in
C++ code , has anyone else encountered such issues?

There are many great resources online for billiards
but I would also like to place the equations of motion
for others to use.

T,

Three cushion billiards fan.

<hr /></blockquote>

tailuge
04-27-2006, 03:02 AM
Jal(Jim) Dr_dave,

Thanks, I think the key phrase that made it
easy for me was 'surface speed'. That allows
me to combine the effect of rotation and linear
velocity in my mind to get one force.

I coded up a few functions (refactored back into
3d vectors) for anyone to take a look. I've added
'time to natural roll' and also not sure if the
sign of 'angularAcceleration' is correct.

To double check, I added some tests taken from
Sheperd, e.g. plain ball -&gt; v=v0*5/7 at Tnr.
And these tests all pass.

Many thanks, I even think the dr_dave web site
should have a section for code since this is
often the modern day way to 'feel' physiscs.

T.

Code + Tests
(to compile use openscenegraph vector class and quicktests)

#include &lt;osg/Geometry&gt;
#include "quicktest/quicktest.h"
#include "ballMath.h"

/*
* Equations of Motion for Sliding/Rolling ball
*/

/*
* V linear velocity
* W angular velocity
* P position
* Pa angular positoin
* ^ is crossproduct
*/

const vector surfaceVelocity(const vector&amp; V,const vector&amp; W)
{
return V + (W ^ vector(0,0,-1.0)) * bconst::radius;
}

const vector linearAcceleration(const vector&amp; V,const vector W)
{
vector Vs = surfaceVelocity(V,W);
return -(Vs/Vs.length()) * bconst::g * bconst::muslide;
}

const vector angularAcceleration(const vector&amp; V,const vector W)
{
vector La = linearAcceleration(V,W);
}

const vector velocityAtTime(const vector&amp; V, const vector&amp; W, double t)
{
return V + linearAcceleration(V,W)*t;
}

const vector angularVelocityAtTime(const vector&amp; V, const vector&amp; W, double t)
{
return W + angularAcceleration(V,W)*t;
}

const vector positionAtTime(const vector&amp; P, const vector&amp; V, const vector&amp; W, double t)
{
return P + V*t + linearAcceleration(V,W)*t*t*0.5;
}

const vector angluarPositionAtTime(const vector&amp; Pa, const vector&amp; V, const vector&amp; W, double t)
{
return Pa + W*t + angularAcceleration(V,W)*t*t*0.5;
}

const double timeToNaturalRoll(const vector&amp; V, const vector&amp; W)
{
vector Vs = surfaceVelocity(V,W);
return (2.0 * Vs.length()) / (7.0 * bconst::g * bconst::muslide);
}

/*
* unit tests
*/

QT_TEST(plainBallToFiveSeventhsV0atTnr)
{
vector v0(-1.0,2.0,0.0);
vector w0(0.0,0.0,0.0);
double tnr = timeToNaturalRoll(v0,w0);
vector v = velocityAtTime(v0,w0,tnr);
QT_CHECK_CLOSEF(v0.length()*5.0/7.0,v.length(),0.001);
}

QT_TEST(dragBallToThreeSeventhsV0atTnr)
{
vector v0(1.0,0.0,0.0);
double tnr = timeToNaturalRoll(v0,w0);
vector v = velocityAtTime(v0,w0,tnr);
QT_CHECK_CLOSEF(v0.length()*3.0/7.0,v.length(),0.001);
}

QT_TEST(masseBallToTwoSeventhsW0atTnr)
{
vector v0(0.0,0.0,0.0);
vector w0(1.0,0.0,0.0);
double tnr = timeToNaturalRoll(v0,w0);
vector v = velocityAtTime(v0,w0,tnr);
}

QT_TEST(surfaceVelAtTnrZero)
{
vector v0(1.0,-3.0,0.0);
vector w0(0.2,4.0,0.0);
double tnr = timeToNaturalRoll(v0,w0);
vector v = velocityAtTime(v0,w0,tnr);
vector w = angularVelocityAtTime(v0,w0,tnr);
vector Vs = surfaceVelocity(v,w);
QT_CHECK_CLOSEF(Vs.length(),0.0,0.001);
}

DickLeonard
04-27-2006, 04:56 AM
Wolfdancer I had read that Joe DiMaggio had the same gift, he knew where the ball was going as soon as it hit the bat and he would lope over and catch the ball. He never made spectacular catches, he was always two steps ahead of the ball and made them look routine.####

wolfdancer
04-27-2006, 06:01 AM
I remember reading something like that...seems like he knew from the sound where the ball was heading...
Folks have speculated about what if... Joe had been traded for Ted Williams....Joe hitting to the wall in fenway, Ted to the short "porch" at Yankee stadium.
Might have worked....but some crazed fan a shot one or both of the owners.......

Jal
04-27-2006, 03:11 PM
<blockquote><font class="small">Quote tailuge:</font><hr>...I coded up a few functions (refactored back into
3d vectors) for anyone to take a look. I've added
'time to natural roll' and also not sure if the
sign of 'angularAcceleration' is correct...<hr /></blockquote>You're most welcome and thanks for taking the time to present the code (although C++ is greek to me).

You appear to be using the factor (- A X R) in the angular acceleration calculation, which looks good to me (correct sign). So does the equation for the time to natural roll, which I neglected to include.

To be honest, although I had to think things through, I believe I cribbed the elegant expression for the surface speed from Dr. Dave.

Welcome to the forum, by the way. I just noticed that these were your first posts.

Jim

tailuge
04-28-2006, 12:51 AM
Jim,

Thanks, I like the idea of using the equations
of motion that ignore 'english' in a straight
line because collisions can be calculated apriori
by using the quartic equation.

I would of course like to understand and implement
the more accurate model but I'm sure the maths is beyond me.

Infact I'm thinking of having an independent vector
for the spin about the z axis and expressing this only
at ball collisions with cushions and other balls.

To this end I've seen some equations about spin induced
throw in some papers online but not sure if I've seen
much about ball-cushion interaction. My experience of
playing with 'check side'(backup english) off a cushion
is that there is a bite point (angle of incidence) beyond which the ball slides
rather than grips.

T.

(only just discovered 3C billiards and

Jal
04-28-2006, 03:40 PM
<blockquote><font class="small">Quote tailuge:</font><hr>
Thanks, I like the idea of using the equations
of motion that ignore 'english' in a straight
line because collisions can be calculated apriori
by using the quartic equation.<hr /></blockquote>You're welcome again Tailuge. /ccboard/images/graemlins/smile.gif It is a very small effect so perhaps best ignored, but if you want my treatment of it I'll be happy to pass it along. It's almost rude to foist math on someone who hasn't requested it specifically, especially when the foister, such as myself, doesn't have the credentials so you can at least have a high degree of confidence in the results (I get things wrong on occasion). So I won't offer up anything unless explicitly asked.

<blockquote><font class="small">Quote tailuge:</font><hr>
I would of course like to understand and implement
the more accurate model but I'm sure the maths is beyond me.<hr /></blockquote>Getting the equations of motion while treating the coefficient of friction u as variable might be possible and maybe even relatively easy - I don't know since haven't tried it. I have done it for throw, where it is highly sensitive to surface speed. However, it requires a method of succesive approximations (Newton's method) to find it. It converges very fast (around six iterations to yield succesive estimates to within 1E-12), but the results are pretty close to Dr. Dave's numbers (see below) and hardly worth the effort given the variable surface conditions on balls anyway. Besides, I don't know if anyone has performed experiments to find the surface speed relationship for the ball/cloth coefficient.

<blockquote><font class="small">Quote tailuge:</font><hr>
Infact I'm thinking of having an independent vector
for the spin about the z axis and expressing this only
at ball collisions with cushions and other balls.

To this end I've seen some equations about spin induced
throw in some papers online but not sure if I've seen
much about ball-cushion interaction. <hr /></blockquote>Dr. Dave's TP A.14 is probably the best treatment you'll see on throw. He incorporates the idea of the cueball reaching roll across the object ball under certain conditions, which is necessary to explain the measured results. He computes a initial value for u (mu) based on the surface speed just before the collision. My version treats it as variable during the collision, but as I said, there's not much difference in the final results. If you want a quick and dirty summary of his equations, let me know.

I started working on cushion collision a while ago but was sidetracked. Of course the big complication is that there is more than one point of contact. Except under very special circimstances, and this isn't one of them, the directions of the forces evolve during impact. I think the simplest approach is to do a numeric integration over time. But since the interaction with the cushion is dominant, you can treat the problem very much like throw, and either ignore or estimate the effects of friction with the bed. This seems to give fairly realistic numbers, but I've only checked it against crudely done tests. Bob Jewett has some carefully done bank results at his web site which I hope to use as a check, but the theory needs more work. I don't have Wayland Marlow's book which might have a thorough treatment.

<blockquote><font class="small">Quote tailuge:</font><hr>My experience of
playing with 'check side'(backup english) off a cushion
is that there is a bite point (angle of incidence) beyond which the ball slides
rather than grips.<hr /></blockquote>When the ball is stunned into the cushion (no draw or follow), I've estimated this angle to be about 55 degrees with respect to the perpendicular. It does of course slide at lesser angles, but reaches a rolling state before it leaves the cushion. Just as with throw, this should mean that the amount of slowing in the parallel direction is quasi-independent of the amount of friction with the cushion, but not with the bed. Or more accurately, while the range of incident angles for which this is true is dependent on the coefficient, when you're well within this range, differences in it have no effect (ignoring bed friction). As I said, I figured the range to be 0-55 degrees for u = .2, and no sidespin (less with check spin), but my undertanding of the whole thing is minimal and ignores various interactions.

Jim

tailuge
04-30-2006, 06:15 AM
Jim,

When I have the mental stamina to go
through a better model for spin I know who

It would be nice to put a library of code
sompleace on the net for others to use and then
we might all get the benefits of better physics
on the games that pop up now and then.

Just to round off my last post I've added the equations
of motion for the rolling without slipping state.

I might post again with the quartic for collision
of two balls moving under the equations described
in this post. If you have any equations on that I'd
like to see it.

It seems to boil down to:

(balls at positions A,B collide at time t when)

resulting in a polynomial with t^4 in it,
the only way I know how to solve that is
with the quartic equation.

T.

/*
* Rolling without slip equations
*/

const vector rollingLinearAcceleration(const vector&amp; V)
{
return -(V/V.length()) * bconst::g * bconst::muroll;
}

const vector rollingAngularAcceleration(const vector&amp; V)
{
vector La = rollingLinearAcceleration(V);
}

const vector rollingVelocityAtTime(const vector&amp; V, double t)
{
return V + rollingLinearAcceleration(V)*t;
}

const vector rollingAngularVelocityAtTime(const vector&amp; V, const vector&amp; W, double t)
{
return W + rollingAngularAcceleration(V)*t;
}

const vector rollingPositionAtTime(const vector&amp; P, const vector&amp; V, double t)
{
return P + V*t + rollingLinearAcceleration(V)*t*t*0.5;
}

const vector rollingAngularPositionAtTime(const vector&amp; Pa, const vector&amp; V, const vector&amp; W,double t)
{
return Pa + W*t + rollingAngularAcceleration(V)*t*t*0.5;
}

const double rollingTimeToMotionless(const vector&amp; V)
{
return V.length()/(bconst::g * bconst::muroll);
}

Jal
05-01-2006, 03:24 PM
Hi Tailuge,

My apology for the delay.

<blockquote><font class="small">Quote tailuge:</font><hr>It would be nice to put a library of code
sompleace on the net for others to use and then
we might all get the benefits of better physics
on the games that pop up now and then.<hr /></blockquote>Not a bad idea. I don't know how much interest there would be, but it could save people from having to re-invent the wheel over and over.

<blockquote><font class="small">Quote tailuge:</font><hr>I might post again with the quartic for collision
of two balls moving under the equations described
in this post. If you have any equations on that I'd
like to see it.

It seems to boil down to:

(balls at positions A,B collide at time t when)

resulting in a polynomial with t^4 in it,
the only way I know how to solve that is
with the quartic equation.<hr /></blockquote>I don't have any equations since it looks like a massive algebraic problem. But what you've indicated above as the starting point looks right to me, but if you could give a few more details on how you arrive at the quartic, I would be interested. I have tried to see if the problem could be simplified, but it would take more sophisticated tools than I possess to pull it off, if it can be done.

Just a few thoughts.

You should be able to use the quadratics in t to see if any two balls, i and j, will get close enough to each other to possibly collide. If they will, although it doesn't guarantee it, there should be a real root for one of the following:

Xi - Xj = + or - 2R
Yi - Yj = + or - 2R

unless the balls are initially within 2R of each other along both axis to begin with, and will never separate by more than that. You then have an approximate time at which impact will occur. You can then use Newton's method to find the root of the quartic, as an alternative if you like, and if it exists.

Frankly, I would calculate their positions at small time increments and check to see how it was going at each step. Have you considered this? I'm not saying it's the best way, just how I would do it.

<blockquote><font class="small">Quote tailuge:</font><hr>
vector La = rollingLinearAcceleration(V);
return (La ^ vector(0,0,-1.0)) *(1.0/bconst::radius);}<hr /></blockquote>Everything else looks spot on but this doesn't seem right. Shouldn't this be the same formula as you used above (with the minus sign and (5/2)factor)?

Edit: Forget the 5/2 factor, a brain blooper on my part.

Jim

tailuge
05-02-2006, 06:12 AM
Hi Jim,

Calculating positions in small steps ends up missing
glancing impacts. To get around that you can use
straight line approximations of the motion between
increments, this approach is probably best if you
want to have some realy fancy physics model for spin.

I've written that kind of system before but this time
I'm creating an 'eventTree' that has nodes at transitions
of state e.g.
t=0.0 ball1 sliding,
t=1.4 ball1 rolling,
t=1.7 ball1 hits ball2
etc etc..

There is a kind of simplicity in this and I'd like a system
where I can make calculations of outcomes quickly without
many iterations (my goal is to measure variation in outcomes
of some 3C shots based on small variation on inputs)

The Quartic (for balls P and Q):
standard equations of motion for balls:
P(t) = pa*t^2 + pv*t + ps
Q(t) = qa*t^2 + qv*t + qs
let D(t)= P(t)-Q(t)
|D(t)|=sqrt( x component squared + y component squared)
distance^2 between balls is Dx(t)^2+Dy(t)^2
Dx(t) is a quadratic so Dx(t)^2 is quartic.
collision are roots of Dx(t)^2+Dy(t)^2 - 4r^2 = 0

tools to see if there was a better approach than
using the quartic equation but my skills were insufficient
to get anywhere.

I've pasted in the code that I use to calculate the
coefficients of this quartic and another function
to evaluate the quartic at time t. I have not included
the quartic solver, but there are many on the net.

T.

/*
In the following balls equations of motion are given by:
P(t) = pa*t^2 + pv*t + ps
Q(t) = qa*t^2 + qv*t + qs
Note: std::vector is an array, vector is a 3d vector
appologies for the name clash.
*/

std::vector&lt;double&gt; quarticCoefficients(
const vector&amp; pa, const vector&amp; pv, const vector&amp; ps,
const vector&amp; qa, const vector&amp; qv, const vector&amp; qs,
const double r)
{
std::vector&lt;double&gt; coeff(5);

vector da = pa-qa;
vector dv = pv-qv;
vector ds = ps-qs;

coeff[0]= ds.x()*ds.x() + ds.y()*ds.y() - 4*r*r;
coeff[1]= 2.0*(ds.x()*dv.x() + ds.y()*dv.y());
coeff[2]= 2.0*(da.x()*ds.x() + da.y()*ds.y()) + dv.x()*dv.x() + dv.y()*dv.y();
coeff[3]= 2.0*(da.x()*dv.x() + da.y()*dv.y());
coeff[4]= da.x()*da.x() + da.y()*da.y();

return coeff;
}

const double evaluateQuartic(const std::vector&lt;double&gt; coeff, double t)
{
return coeff[4]*t*t*t*t + coeff[3]*t*t*t + coeff[2]*t*t + coeff[1]*t + coeff[0];
}

/*
example test, also shows how balls initial conditions
become transformed in to quartic coefficients
*/

QT_TEST(ballsMovingTogetherHaveCollideTime)
{
vector qs(0,0,0);
vector qv(1,0,0);
vector qw(0,0,0);

vector pv(-1,0,0);
vector pw(0,0,0);

std::vector&lt;double&gt; coeffs = quarticCoefficients(qs,qv,linearAcceleration(qv,qw )*0.5,

double t = minPositiveRoot(coeffs);
QT_CHECK_GREATER(t,0.0);
}

tailuge
05-02-2006, 03:10 PM
Appologies,

A little testing revealed I got the argument
list in the wrong order, should be:

std::vector&lt;double&gt; quarticCoefficients(
const vector&amp; ps, const vector&amp; pv, const vector&amp; pa,
const vector&amp; qs, const vector&amp; qv, const vector&amp; qa,
const double r)
{
std::vector&lt;double&gt; coeff(5);

vector da = pa-qa;
vector dv = pv-qv;
vector ds = ps-qs;

coeff[0]= ds.x()*ds.x() + ds.y()*ds.y() - 4*r*r;
coeff[1]= 2.0*(ds.x()*dv.x() + ds.y()*dv.y());
coeff[2]= 2.0*(da.x()*ds.x() + da.y()*ds.y()) + dv.x()*dv.x() + dv.y()*dv.y();
coeff[3]= 2.0*(da.x()*dv.x() + da.y()*dv.y());
coeff[4]= da.x()*da.x() + da.y()*da.y();

return coeff;
}

Jal
05-04-2006, 03:42 PM
<blockquote><font class="small">Quote tailuge:</font><hr>
The Quartic (for balls P and Q):
standard equations of motion for balls:
P(t) = pa*t^2 + pv*t + ps
Q(t) = qa*t^2 + qv*t + qs
let D(t)= P(t)-Q(t)
|D(t)|=sqrt( x component squared + y component squared)
distance^2 between balls is Dx(t)^2+Dy(t)^2
Dx(t) is a quadratic so Dx(t)^2 is quartic.
collision are roots of Dx(t)^2+Dy(t)^2 - 4r^2 = 0

tools to see if there was a better approach than
using the quartic equation but my skills were insufficient
to get anywhere.<hr /></blockquote>Hi Tailuge,

Again, sorry for the delay and thanks for posting this and the code. It's very neat and succint...bravo. By now you may have gone to annother approach to solving it other than the quartic formula(s), or maybe have resolved the issues with it. But I've been doing some further reading on Newton's method and thought I'd pass on a couple of things. In my opinion, this is much easier than the quartic formula(s).

A quick review:

In general, a root of the equation f(x) = 0 can be obtained with successive approximations, x0, x1, x2, ... using:

xn+1 = xn - f(xn)/f'(xn)

where xn+1 indicates the next better approximation, xn is the current approximation, f(xn) is the function evaluated at xn, and f'(xn) is the function's first derivative at xn.

So you start with some initial guess x0. If x0 is sufficiently accurate, you get two good digits with every iteration (quadratic convergence).

For the quartic, letting t stand for the current estimate:

f(t) = c4(t^4) + c3(t^3) + c2(t^2) + c1(t) + c0 and

f'(t) = 4c4(t^3) + 3c3(t^2) + 2c2(t) + c1

There are extensions to Newton's method which produce faster convergence, such as 3, 4, 5 or more digits per iteration. But the simple method above seems to beat them all (for the quartic) in terms of the number of arithmetic operations it takes to reach a certain precision. I counted the number of them as follows (it includes the initial overhead to compute the modified coefficients for the derivative):

6 digits: 33X,24+-,3/
8: 43X,32+-,4/
10: 53X,40+-,5/
12: 63X,48+-,6/

You might want to try to compare this to the number of operations for the quartic formula you're using. This will also involve recursion in evaluating the radicals. Without knowing the algorithms, we can only say that it will add to it.

In the vast majority of cases, an initial estimate of t = 0 will probably be good enough, although I haven't tested it. In general, when two balls collide, there should be two real roots, one for the initial collision and the other as they "pass through" and exit each other. (If they're moving very slowly so that they'll come to a stop while overlapping, then only only one real positive root will exist.) So the smallest positive root will be the desired one, if one exists. Newton's method should converge to it with an initial guess of t = 0, unless the following is possible.

They may collide in negative time under some circumstances, so you could get convergence to one of these roots. If other considerations eliminate the possibility of a positive root, fine, no collision will take place in real time. Otherwise you'll need a better intial guess at t, perhaps using the quadratic condition indicated earlier.

A watch has to be kept on f/f' - it tells you if a collision will happen. When one ball passes another, the derivative f' will approach zero while f will remain relatively large. In the case of a grazing collision, f/f' should behave itself, but I'm not sure if this will be true if it passes at exactly 2R (within computer precision)). Some testing or more research is in order.

There is a test which is supposed to guarantee that a root exists and that Newton's method will converge on it. I don't know much about it but you can start at the following if you're interested:

http://numbers.computation.free.fr/Constants/Algorithms/newton.html

I expect the problems you're having with the quartic involves distinguishing the real from the complex roots (but you probably know more about it than me). I think Newton's method is much more straight forward and possibly more computationaly efficient.

Just some thoughts.

Jim

tailuge
05-05-2006, 03:56 AM
Jim,

Firstly I can imagine a situation with 4 real positive
roots. Imagine you shoot a masse shot 1 foot up the table
and it spins back to your starting position. Place an
object ball half way along this path and you will get
an entry/exit collision on the way up, and another
entry/exit collision on the way back.

I would like to use Newtons method to find the
roots but it seems like it could easily wander off
to -ve roots or 2nd positive root etc, I guess this

As for the quartic equation, its much like the quadratic
equation, a few lines of code (see wikipedia for a good
example). I think the problem I'm hitting now is
numerical computational because even though I can
get the 1st positive root, when I place the balls
at the time where that root predicts they are not
exactly in contact.

I did a test where one ball is shot into another
and the object ball is put progressively further away
so the root(collision) gets larger each time.

As you can see from the data, the quartic evaluates to
zero at the root, but the balls do not have 0 seperation
at that time (ball positions are quadratics in t)

Thats because the numerical accuracy of 64-bit floats
shows up a diffence between the quadratic and quartic
evaluation of the same underlying calculation.

(The last two lines have no root because
the object ball was too far away)

So what I conclude I need is a 'BlackBox' root refiner
that takes the function of seperation between the balls
and uses that to get to a value for t where the balls
are close to within some tolerence.

Once again its fiddly because the 'refiner' has to avoid
wandering off into another root.

Any thoughts would be most welcomed, I'll try to reply with
by black box code - we have decended from the purity of
spheres and maths into ...bits.

T.

// in the following each line has the
// object ball further away, until no
// collisions occur. i is just the line count
// t = minPositiveQuarticRoot(coeffs);
// quartic = evaluateQuartic(coeffs,t);
// seperation = seperationOfBallsAtTime(P,Q,t);

[---------------- RUNNING TESTS ----------------]

i= 1 t=0.3115508386 quartic=-0.0000000000 seperation=-0.0611997085
i= 2 t=0.5922976179 quartic=0.0000000000 seperation=-0.1976227453
i= 3 t=0.8976896748 quartic=0.0000000000 seperation=-0.3058225917
i= 4 t=1.2358063598 quartic=-0.0000000000 seperation=-0.0392693843
i= 5 t=1.6205537415 quartic=0.0000000000 seperation=0.8031743960
i= 6 t=2.0801449652 quartic=0.0000000000 seperation=2.2924946255
i= 7 t=2.6950482176 quartic=0.0000000000 seperation=4.9631832987
i= 8 t=-1.0000000000 quartic=1070.3850904217 seperation=31.0874944457
i= 9 t=-1.0000000000 quartic=1240.2948804608 seperation=33.5870448561

tailuge
05-05-2006, 11:47 PM
Update,

I wrote a refiner that succesfuly gets actual
ball seperation down to a given tolerence and
seems quite stable.

After some perplexing results I found an error
in one of my seperation calculations and once
fixed the quartic roots and seperation results
were all much closer.

From here what I would like is some hard test cases
where I could see if the machine agrees with some
outside calculation. Perhaps I could try one of my
algebra packages and see if I can generate a test case.

T.

//
// Look for solution to sepObj(t)=0
// such that sepObj(t) &gt;= 0 and sepObj(t) &lt; tolerence
// sepObg must be continuous
//

const double refineMinPositiveValue(sepFunctor* sepObj,const double t0,const double tolerence)
{
// test tolerence

double f0 = (*sepObj)(t0);
if ((f0&gt;=0.0) &amp;&amp; (f0 &lt; tolerence)) return t0;

double t1 = t0 * 0.999;
double f1 = (*sepObj)(t1);

if (f1==f0) return -1; // cannot refine flat line

// use gradient and sign of f0 to see if solution &lt; or &gt; t0

double searchside = 0;

if ((f0 &lt; 0) xor (gradient &lt; 0)) searchside = +1;
else searchside = -1;

// now must find a point p where sign(fPtr(p)) != sign(f0)
// then use bisection method

// gradient gives good search step choice

double guess = t0;

for(int i=0; i&lt;4; i++)
{
guess += tstep*searchside;
double v = (*sepObj)(guess);

if ((v &gt;= 0) xor (f0 &gt; 0))
return bisectToTolerence(sepObj,guess,t0,tolerence);
}

// error did not find candidate bisector point

return -1.0;
}

Jal
05-07-2006, 01:24 PM
Hi Tailuge,

<blockquote><font class="small">Quote tailuge:</font><hr>As for the quartic equation, its much like the quadratic
equation, a few lines of code (see wikipedia for a good
example). I think the problem I'm hitting now is
numerical computational because even though I can
get the 1st positive root, when I place the balls
at the time where that root predicts they are not
exactly in contact.

I did a test where one ball is shot into another
and the object ball is put progressively further away
so the root(collision) gets larger each time.

As you can see from the data, the quartic evaluates to
zero at the root, but the balls do not have 0 seperation
at that time (ball positions are quadratics in t)

Thats because the numerical accuracy of 64-bit floats
shows up a diffence between the quadratic and quartic
evaluation of the same underlying calculation.

(The last two lines have no root because
the object ball was too far away)

So what I conclude I need is a 'BlackBox' root refiner
that takes the function of seperation between the balls
and uses that to get to a value for t where the balls
are close to within some tolerence.

Once again its fiddly because the 'refiner' has to avoid
wandering off into another root.

Any thoughts would be most welcomed, I'll try to reply with
by black box code - we have decended from the purity of
spheres and maths into ...bits.

T.

// in the following each line has the
// object ball further away, until no
// collisions occur. i is just the line count
// t = minPositiveQuarticRoot(coeffs);
// quartic = evaluateQuartic(coeffs,t);
// seperation = seperationOfBallsAtTime(P,Q,t);

[---------------- RUNNING TESTS ----------------]

i= 1 t=0.3115508386 quartic=-0.0000000000 seperation=-0.0611997085
i= 2 t=0.5922976179 quartic=0.0000000000 seperation=-0.1976227453
i= 3 t=0.8976896748 quartic=0.0000000000 seperation=-0.3058225917
i= 4 t=1.2358063598 quartic=-0.0000000000 seperation=-0.0392693843
i= 5 t=1.6205537415 quartic=0.0000000000 seperation=0.8031743960
i= 6 t=2.0801449652 quartic=0.0000000000 seperation=2.2924946255
i= 7 t=2.6950482176 quartic=0.0000000000 seperation=4.9631832987
i= 8 t=-1.0000000000 quartic=1070.3850904217 seperation=31.0874944457
i= 9 t=-1.0000000000 quartic=1240.2948804608 seperation=33.5870448561 <hr /></blockquote>This is kind of perpexing although I see below that you may have fixed it. The separation is the square root of the quartic, so how is it possible that you're getting such differing values? I really don't think it's a floating point precision problem - I would think the separation values would be much closer to zero even with single precision (32-bit) arithmetic. When you take the square root of the quartic, you only lose about half the leading zeros on the right side of the decimal point. There must be some inconsistency in the way you're figuring the separation.

Jim

Jal
05-07-2006, 01:43 PM
<blockquote><font class="small">Quote tailuge:</font><hr> ...I wrote a refiner that succesfuly gets actual
ball seperation down to a given tolerence and
seems quite stable.

After some perplexing results I found an error
in one of my seperation calculations and once
fixed the quartic roots and seperation results
were all much closer.

From here what I would like is some hard test cases
where I could see if the machine agrees with some
outside calculation. Perhaps I could try one of my
algebra packages and see if I can generate a test case.
<hr /></blockquote>It's not clear if you still need this "refiner" to get acceptable results after fixing the separation equation. I'm pretty sure it shouldn't be necessary but maybe you can supply more details. The way you arrived at the quartic in the first place seemed right to me, which as mentioned above, is the square of the separation.

Also, could you tell me how you get the roots of the quartic. Are you using a math package and does it return them in complex form, or does it only give the real roots if they exist?

As far as testing, how about setting the balls up so they're just touching at various points along their circumferences, giving them various speeds and spins, then back computing their states at some negative time, and then starting with that for their initial states? They should collide every time. /ccboard/images/graemlins/smirk.gif

Jim

tailuge
05-07-2006, 06:24 PM
Jim,

I'm currently using the c++ quartic solver

http://www.pearsoned-asia.com/rao/pl2.htm

It returns 4 roots, with real and imaginary
components and I just choose the least
positive real root.

The refiner code is not required now
that my bug is fixed but I still use it
to ensure that the balls do not have any
overlap before begining the subsequent calculation
of collision velocities.

I also use the refiner to ensure the ball-cush
collision time has no overlap (after solving the

I have a bunch of unit tests that give me a basic
level of correctness but I've yet to draw anything
on screen to get that final reassurance nothing
odd happens.

I now have an 'event exapander' that takes some initial
events (ball A sliding towards ball B) and constructs a
tree of subsequent events (with each event being a
record of initial physical conditions at a time)

A:slide B:atrest
=&gt;
A:slide-&gt;A:roll-&gt;A:impact B
A:slide B:slide
A:slide-&gt;A:roll-&gt;A:atrest
B:slide-&gt;B:roll-&gt;B:atrest

It's nice because there is very little iteration
required to find out the outcome of a particular
shot.

How would you write the cushion reflection equations?
Currently for a ball with velocity V, angular velocity W
impacting a cushion aligned with the Y axis I do

V.x=-V.x * cushioncoeff (e.g. 0.9)
W.y= W.y * cushioncoeffspin (e.g. 0.3)

(.x means x component - computer speak not math sorry)

This allows for 'banana shots' off cusions and feels
right but I would much prefere a model grounded in some physics.

T.

Jal
05-07-2006, 09:49 PM
Hi tailuge,

A quick response for a change, but of course no answer to your cushion problem yet. Let me get back to you, within a day (or two) I should hope.

Thanks for answering my questions. I take it then that you look for the roots where the imaginary part is below some threshold?

As far as cushion reflection, at this point I can only offer something very similar to throw, which I'm pretty familiar with, and which ignores bed friction and its feedback with the cushion. But I would like to look it over a little first.

Typically, .7 is cited as the coefficient of restitution in the perpendicular direction (Vx ---&gt; -.7Vx). But after doing measurements on some of Dr. Dave's high speed videos involving straight on shots into the cushion, although the measurements were scattered quite a bit, the averages seem to point to a coefficient of .8, with bed friction accounting for another drop of .1. So effectively it's .7 for straight on shots, but then if you accept say, .8, it should be between .7 and .8 for angled approaches because of the varying directions of the bed friction (which also changes during impact). These figures are for pool tables, and specifically for the one that Dr. Dave demonstrations were performed on. I'll have to look up billiard table specs.

As far as friction with the cushion in the parallel direction, it's going to depend on spin, as with throw, and a ball will often slide for a while and then end up rolling across it. If it retains some topspin it will also curve after leaving the cushion. The "vertical" spin that it develops from cushion friction will also contribute to the curve since it's not perfectly vertical.

I have sort of worked it out for the case of a pure stun shot into a cushion, but the treatment of the bed friction is not too simple and still completely disregards feedback with the cushion.

So give me a moment and I should be able to give you the simple version shortly. Of course, I'm sure you wouldn't mind someone else chiming in here for a change.

Jim

tailuge
05-11-2006, 08:51 AM
Jim,

I've been busy writing the code to view the
balls moving and immediatly uncovered some
problems.

roots of the quartic (imag=0). The big issue that
took me this long to track down is how the algorithm
behaves when the t^4 coefficient is zero (I believe
the correct term is 'depressed quartic' certainly made
me feel down for a day or two).

This can happen when one ball is rolling quickly into the back of another rolling more slowly. The quartic reduces to
a quadratic and the code must deal with this case.

One other issue was that refining the seperation had to
produce a +ve seperation. seperation=0 caused my setup to
stick balls together.

I have a system that appears to work now, I would realy
like to plug in any accurate models you have for throw
on ball/ball ball/cush. If they are
presented in terms of (pos,vel,angualrv) at point of
contact then I might just be able to understand them.

Thanks for all your help so far, I can share more of the
code with you at your request.

T.

dr_dave
05-12-2006, 09:37 AM
T,

Concerning solution of the quartic equations, I think the Canadian research group developing the pool-playing robot has the equations and their solution worked out in detail. Greenspan (the principal investigator) might have some papers and/or code that he is willing to share. See "Canadian pool-playing robot" under "videos" in the threads summary area of my website (http://www.engr.colostate.edu/~dga/pool/threads.html) for more info. At the same webpage, you can also find links to postings concerning the coefficient of restitution (under "bank and kick shots") that might answer some of your questions there.

Regards,
Dave

Jal
05-13-2006, 01:16 AM
Hi Tailuge,

This time a double apology for the delay, since I said it would only be a day or two. At this point I'll just say I'm sorry in advance and you can assume it will take me forever to respond.

I've been working on extending the treatment of a stun shot into a cushion to the general case of any incoming spin. When I started I hoped it would be a quick job, but not so. It's not that the math is that tough but sometimes (often) my cogitating is slow. As of today I've gotten the equations down, but would like to test them against some of Bob Jewett's data. If I don't have anything to offer soon (say tomorrow evening?), I'll at least give you the simple treatment.

In the meantime, as far as ball/ball throw, perhaps you might peruse Dr. Dave's TP A.14. What I will offer is his treatment, but generalized for two balls colliding while in motion. I think it's in the form you're looking for.

Since you combined the earlier component equations into vectors, I'll stick with them here. They are in uppercase bold. Unit vectors are in lower case bold, and vector magnitudes will either be in normal text or surrounded by |'s, as in |V|, as seems appropriate. The vertical unit vector is z.

Let's arbitrarily give one of the balls "cueball" status, and the other, "object" ball status. At the moment of collision they have position Pc, Po, velocity Vc, Vo, and spin Wc, Wo, respectively.

The radius vector from the center of the cueball to the contact point is then:

R = (1/2)(Po - Pc), and its unit vector is
r = R/r

It's convenient to switch to a frame of reference moving along with the object ball at velocity Vo, so that the object ball is at rest (but likely spinning) in it. The cueball's velocity in this frame is:

V = Vc - Vo and its unit vector is
v = V/|V|

The sine and cosine of the cut angle are (where * indicates the vector dot product):

sin(theta) = (v X r)*z
cos(theta) = v*r

If sin(theta) is positive, it's a cut to the left, otherwise to the right. The cosine is always positive.

The relative surface speed at the contact point is:

Vs = V + Wc X R - Wo X (- R) or eliminating the minus signs

Vs = V + Wc X R + Wo X R

Their relative surface speed in the tangential (friction) plane is:

Vsf = r X Vs X r, with unit vector
vsf = Vsf/|Vsf|

The double cross product can be executed in any order, it doesn't matter in this case.

The horizontal tangent line unit vector in the friction plane and the cosine of the angle the relative surface speed makes with it are:

t = r X z
cos(phi) = vsf*t (again, dot product)

Okay, with that we can figure the post impact states of the two balls. The tangent of the throw angle gamma is:

tan(gamma) = uqcos(phi)

where u is the coefficient of friction (calculated as below), q is the fraction of the compression impulse involved in the actual throw (see below), and cos(phi) was just calculated. The sign of cos(phi) will produce the correct sign for tan(gamma).

The balls may slide during the entire impact period or end up rolling across each other (after first sliding). They'll slide for the entire period if:

u|V|cos(theta) &lt; (1/7)|Vsf|

in which case q = 1. Otherwise,

q = (1/7)|Vsf|/[u|V|cos(theta)]

and the amount of throw is independent of u. The first case where q = 1 is the more common occurence. This is essentially Dr. Dave's treatment in TP A.14. In it, he calculates u according to:

u = a + b(exp(-c|Vsf|))

where exp() is the function taking the base of the natural logarithm to the power of the argument. The constants a,b,c are:

a = 9.951E-3, b = .108, c = 1.088

based on Wayland Marlow's data.

Now for the post impact states of the balls. Let e stand for the coefficient of restitution between them. It's approximately .95, but according to tests done by another poster here (Mac, aka Cushioncrawler), it varies a bit with impact speed, |V|cos(theta) in our case.

The object ball's post impact speed along the line of centers (r) in our moving frame is:

Vor = (1/2)(1+e)|V|cos(theta)

The cueball will retain a small amount of speed in this direction as well, which is:

Vcr = (1/2)(1-e)|V|cos(theta)

Their post impact states in the stationary (table) frame of reference are then:

Vo' = Vo + Vor(r + tan(gamma)t)

Vc' = Vo + Vcrr + |V|sin(theta) - (Vor)tan(gamma)t

Wo' = Wo - (5/2)(1/r)(uq)(Vor)r X vsf

Wc' = Wc - (5/2)(1/r)(uq)(Vor)r X vsf

The change in spin of both balls is the same. The uq factor is missing the cos(phi) in the spin equations because we want the change in surface speed in the direction of the surface speed, not in the horizontal direction as per throw. The 5/2 factor comes from the fact that the change in surface speed of each ball is 5/2 the change in its translational speed. The total change in surface speed due to changes in both balls' speed and spin is 7 times the change in one ball's translational speed (5/2+5/2+1+1). That's where the number 1/7|Vsr| arises when determining whether or not the balls will roll across each other.

While we're at it, there is something else which affects their mutual departure. Because of the finite contact time, the direction between the ball's centers as well as the tangent line shift during impact. As a result the cut angle increases a little. By my figuring this can be more than the throw angle at fairly severe cut angles. This was brought up on another forum (RSB) by Bob Jewett and the numbers he cited we're much smaller than mine. However, I still feel the following is correct, but we shall see. Without derivation (but will supply on request), the added cut angle alpha is given by:

alpha = arctan[sin(theta)/(2R/(|V|T) - cos(theta))]

where T is the contact time, approx. .0002 sec.

This is the amount by which the line of centers and the tangent line are rotated. We can incorporate this into the above equations by rotating the unit vectors t and r by alpha. We do this by figuring their components in a coordinate system rotated by - alpha. So let

al = - alpha

And let the components of t and r be tx,ty,rx and ry, while i and j are the standard unit vectors. After rotation they are:

t = (txcos(al) + tysin(al))i + (tycos(al) - txsin(al))j

r = (rxcos(al) + rysin(al))i + (rycos(al) - rxsin(al))j

These are to be applied before using them in any of the above formulas of course.

I'm glad you're making progress with your code. The reason I was curious about the quartic is that in my CRC handbook, the real roots of the cubic and quartic can contain complex numbers. The imaginary parts cancel but it takes some extra effort to do so. It's clear now that this wasn't your problem.

I haven't really looked at your refiner yet in detail, so I can't offer applause, as I'm sure it's due. This supposes that I could understand it. The C++ code, as I indicated earlier, is pretty foreign to me. All I can do is to try to get the gist.

If you have questions on any of the above, fire away. Please be critical of it. If anything doesn't seem to make sense, let me know.

Jim

Jal
05-13-2006, 11:32 AM
<blockquote><font class="small">Quote jal:</font><hr>Vc' = Vo + Vcrr + |V|sin(theta) - (Vor)tan(gamma)t<hr /></blockquote>The third term has no unit vector which should be t. I forgot the enclosing parenthesis. So it should read:

Vc' = Vo + Vcrr + (|V|sin(theta) - (Vor)tan(gamma))t

Jim

tailuge
05-14-2006, 12:56 AM
Jim,

An amazing post, thanks. I wanted to try it out
straight away and mostly just had to cut and paste
to get it running.

My main test case is:

Ball C fired full ball into ball O with no spin
on either ball. Balls should remain on straight line
after collision (no throw).

This causes Vsf to be zero and the computer has a bit
of trouble normalising that so I did the following:

vector Vsf = r ^ Vs ^ r;
vector vsf = Vsf/Vsf.length();
if (Vsf.length()==0) vsf = vector(0,0,0);

Which I think is ok.

Other than that looking at what happens to the balls
when knocking them around seems to be fine. I can
switch back to the simpler model for collisions and I
will attempt to draw the outcome for a half ball

1) simple model
2) plain ball throw model
3) topspin throw model
4) backspin throw model

(might be the best I can offer in terms of testing)

half of your code but will try that as soon as I can get
some simple way of viewing the results. I'll post
the code when I've done some more tests.

As for refining the roots, I think given that we are
looking at 'first' impact it should be the case that
the gradient has to be negative at that root and so the
code could be simplified further.

When it comes to drawing the balls on screen I'm struggling
with getting the 'angular position' converted into a
matrix rotation. My current code looks like this:

vector axis = angularpos/angularpos.length();
osg::Matrix rt = osg::Matrix::rotate(rev,axis);

The rotate matrix rotates by 'rev' radians around 'axis'.
It doesnt look right on screen - I'm sure its wrong.

plenty to do and its been fascinating so far.

T.

SPetty
05-14-2006, 08:37 AM
<blockquote><font class="small">Quote tailuge:</font><hr>An amazing post, thanks.<hr /></blockquote>I just gotta say that this entire thread is amazing. /ccboard/images/graemlins/smirk.gif

Fran Crimi
05-14-2006, 10:11 AM
<blockquote><font class="small">Quote SPetty:</font><hr> <blockquote><font class="small">Quote tailuge:</font><hr>An amazing post, thanks.<hr /></blockquote>I just gotta say that this entire thread is amazing. /ccboard/images/graemlins/smirk.gif <hr /></blockquote>

I agree. /ccboard/images/graemlins/smirk.gif Just so happens I'm presently giving lessons to the Chairman of the Physics dept. of a famous university (confidential). He said his area of expertise is 'chaos'. I told him to come to the CCB. /ccboard/images/graemlins/grin.gif

I figure there may be some stuff here that might interest him. Plus, maybe he can translate some of it.

Fran

Jal
05-15-2006, 12:41 AM
<blockquote><font class="small">Quote tailuge:</font><hr> Jim,

An amazing post, thanks. I wanted to try it out
straight away and mostly just had to cut and paste
to get it running.

My main test case is:

Ball C fired full ball into ball O with no spin
on either ball. Balls should remain on straight line
after collision (no throw).

This causes Vsf to be zero and the computer has a bit
of trouble normalising that so I did the following:

vector Vsf = r ^ Vs ^ r;
vector vsf = Vsf/Vsf.length();
if (Vsf.length()==0) vsf = vector(0,0,0);

Which I think is ok.

Other than that looking at what happens to the balls
when knocking them around seems to be fine. I can
switch back to the simpler model for collisions and I
will attempt to draw the outcome for a half ball

1) simple model
2) plain ball throw model
3) topspin throw model
4) backspin throw model

(might be the best I can offer in terms of testing)

half of your code but will try that as soon as I can get
some simple way of viewing the results. I'll post
the code when I've done some more tests.

As for refining the roots, I think given that we are
looking at 'first' impact it should be the case that
the gradient has to be negative at that root and so the
code could be simplified further.

When it comes to drawing the balls on screen I'm struggling
with getting the 'angular position' converted into a
matrix rotation. My current code looks like this:

vector axis = angularpos/angularpos.length();
osg::Matrix rt = osg::Matrix::rotate(rev,axis);

The rotate matrix rotates by 'rev' radians around 'axis'.
It doesnt look right on screen - I'm sure its wrong.

plenty to do and its been fascinating so far.

T.

<hr /></blockquote>Thanks for your kind remark Tailuge.

I also think that setting vsf to zero is the way to go.

As far as the gradient in the refiner code, I still haven't gone over it yet so I can't remark at this time. I hope to do so soon.

Just one thing concerning the position rotation. It looks like you're dividing the angular position by two times the radius to get the angle. While I'm not sure what's going on here, shouldn't this be just one times the radius?

This is beginning to get very interesting. I had no idea you were already to the point where you have the balls displayed on screen. I would love to see the results, if there's any way to do this, whenever you feel it's in a form you're satisfied with. (Crude is okay with me.) Not necessary, just an idle wish.

Jim

tailuge
05-16-2006, 05:57 AM
Progress so far is that I can visualise the outcome
of a shot based on a small error in aim/spin/power.

http://www.geocities.com/luktay/t2.gif

Not happy with the 'feel' currently, I have:

const double tablewidth = tablelength/2.0;
const double muroll = (1.0/130.0);
const double muslide = (1.0/5.0);
const double maxvel = 120.0 * radius;

T.

(oh and I just discovered Page2
in this forum, sorry for not replying sooner)

Jal
05-16-2006, 10:49 AM
Thanks for that image. It's looking pretty good - you even have everything in perspective!

If you know you should be dividing by 2*R, fine. I just wasn't sure what you were doing there.

I sent you a private message concerning your refiner code. Check "My Home" at the top of the page.

The program to test the cushion reflection is just about done, so it should be soon, wihtin a day I hope, but you know how reliable that is by now.

Jim

Bob_Jewett
05-16-2006, 03:38 PM
<blockquote><font class="small">Quote tailuge:</font><hr> ... I would really like to have an equation that can
predict the position of a ball at time t based
on an initial velocity vector and rotation vector.

pos(t)= func(t,V,W);
rot(t)= func(t,V,W);

... <hr /></blockquote>
Do you want to include the fact that ball-table friction is probably not constant vs. surface velocity? I think this means that the path before smooth rolling is not a parabola.

You could just integrate the differential equations. That way your solution can be more correct, even including the fact that during a throw shot, the relative velocities of the two surfaces can go to zero during contact (as Dr. Dave has pointed out), and that the friction between the balls depends on surface velocity, as Marlow measured.

You may also want to include the fact that on nearly all shots the cue ball leaves the surface of the table, and on many shots, the object ball leaves the table. The first fact is well illustrated on the DVD by Efler and Leitner in which you can see a bouncing ball with draw (or follow?) move in straight line segments as viewed from above rather than a parabola. This effect is modelled in Virtual Pool.

What Virtual Pool does not model well is what happens when the cue ball hits two balls that are touching. I think integration would let you handle this case as well.

wolfdancer
05-16-2006, 03:51 PM
What's really amazing to me...is I haven't a clue about what they are talking about.....and I got my GED on only the third try....or was it the fourth?
It sounds very interesting though, and I can see the appeal for eng/math/physics/ persons.

Bob_Jewett
05-16-2006, 04:15 PM
<blockquote><font class="small">Quote tailuge:</font><hr> ... Three cushion billiards fan.

<hr /></blockquote>
I just noticed that part. If you are interested in simulating carom billiards, you may want to contact Regis Petit who has written a book on billiard physics called "Billard -- Theorie du Jeu" ISBN 2-7027-0573-1, Editions CHIRON/CASTEILLA.

That work and Coriolis have been implemented in a simulator called "Coriolis 3D" by Guy Grasland and Jean-Luc Frantz. See http://vrbilliard.tripod.com/practice.html for some more info and example screen shots.

tailuge
05-16-2006, 05:53 PM
Jim,

Drawing everything is the easiest part, just using
the open source project 'openscenegraph'. When you
write your programs what language are you using?

I'd like to place the source for what I've written for
others to improve or use, but its still a bit of a mess.

As for the ball rotation on screen, it's not correct.
The angular position describes how the ball is oriented
about the x,y axis. But you cant combine the rotations
because they are order dependent i.e.

Matrix(rotx)*Matrix(roty) != Matrix(roty)*Matrix(rotx)

The problem I need to solve is:

If we say that the position of the red dot on
the white ball can be described completely
by the angular position then there must be a way
to work out the (x,y,z) of the dot relative to the
balls center from (angularx,angulary).

T.

Jal
05-17-2006, 05:39 PM
Hi Tailuge,

My brain turned to mush while working on something else. Then I read your dot problem....oh boy, things just got more complicated. I think you'll need to do a timewise numeric integration for the general case where the spin axis is changing direction. I have some equations but I'm wondering if there's not a better way to do it than what I currently have down. Let me think about it some more. Maybe you've resolved the issue by now.

I just wanted to note a couple of things regarding the earlier formulas.

First, when calculating u (friction coefficient), the units for the surface speed are meters/sec.

Second, if you use the added cut angle alpha, it's sign should the same as sin(theta), and then al has the opposite sign. And I think you could apply the rotation to R when first calculated, saving having to use it on t, but you might want to give that more thought.

Jim

tailuge
05-18-2006, 06:45 AM
Jim,

&gt; timewise numeric integration

I should really learn about that, seems
to crop up over and over in pool modelling.

I think I've got a poor mans solution that
amounts to the same thing. I keep a second
angular position which starts as a unit vector
which I rotate about an axis defined by angular
velocity by amount proportional to time elapsed
since last frame. (looks correct)

The problem is interesting because the angular position
we maintain by the equations in this post should be
enough if only we knew the mapping. Its not spherical
coordinates, not sure what to search for in google.

Or is it that the angular position as I posted in my
original code is nonsense? An angular position with
z always zero is a bit suspect.

T.

Jal
05-18-2006, 11:46 PM
Hello Tailuge,

The problem with using:

angle position=a0+vt+(1/2)at^2

when the axis is changing orientation is that the reference frame (axis) has moved. When you find the magnitude of the position you get some angle, but with respect to what? So I think you need to do something like this numeric integration in order to keep track of a point on the ball's surface.

Let the position of the point be Ps at some time t. The ball's position is P. The vector from the ball's center to this point is:

Rs = Ps - P

The component of this perpendicular to the spin vector W is:

Rz = w X Rs X w with unit vector
rz = Rz/|Rz|

where w is the unit vector of W.

The speed of the point is:

Vs = V + W X Rz

Its acceleration As is:

As = dVs/dt = A + dW/dt X Rz + W X dRz/dt

The first term is the acceleration of the ball. The second is the point's acceleration due to torque from friction with the bed. Evaluating it:

dW/dt X Rz = (5/2)(1/r)(rb X A) X Rz

where rb is the unit vector to the tables surface and r is the ball's radius.

The third term is the centripetal acceleration and we can use one of these expressions:

W X dRz/dt = W X Vs = -(|Vs|^2/|Rz|)rz

Combining the above we have for the acceleration of the point:

As = A + (5/2)(1/r)(rb X A)X Rz - (|Vs|^2/|Rz|)rz

To do the integration let d_t be some small interval of time, say 1/10000 sec, or whatever works best. Then the position of the point after d_t will be:

Ps' = Ps + Vs(d_t) + (1/2)As(d_t^2)

We're treating the acceleration as constant over this small time interval. It really isn't but the smaller the time step the more true this is. We could differentiate the expression for As to get, say, As' and then a better approximation to Ps' would be:

Ps' = Ps + Vs(d_t) + (1/2)As(d_t^2) + (1/6)As'(d_t^3)

but you get a very long expression for As', unless there is some simplification which I overlooked.

Once the new position is obtained, we don't want the point to leave the surface of the ball as small errors accumulate. We can make sure it doesn't by constraining it to a distance r from the ball's center in the direction of its unit vector.

Ps' ---&gt; (r)(Ps' - P')/|Ps' - P'| + P'

Now we have to go back to the beginning and recompute Rs, Rz,As, etc, to obtain the next position using the updated value of W . This continues until something of significance happens.

If the axis has a fixed orientation but W is accelerating, you can keep track of the angular position with the simpler method. But this is a minority situation since some sidespin will usually be present, which decays with time.

Jim

tailuge
05-20-2006, 11:48 PM
Jim,

Currently I have a very simple bit of code
to maintain the rotation. It calculates how much
the ball rotates per video frame (from angularv)
and accumulates that rotation across frames.

I didnt implement your suggestion yet, because I
also think that there must be a better way. If
linear position can be calculated (pos(t)) simply from
the equations of motion then surely angular position
must have an analagous form. (But there my brain
quits)

My code:

void ball::updateRendered(const double deltat)
{
// translation matrix
osg::Matrix mt = osg::Matrix::translate(pos);

// rotation required for this frame
vector raxis = angularv;
raxis.normalize();
osg::Matrix rt = osg::Matrix::rotate(deltat * angularv.length(),raxis);

// accumulate rotation
totalrot = totalrot * rt;
osg::Matrix combined = totalrot * mt;

// rotate and translate ball but only translate shadow
xformBall-&gt;setMatrix(combined);
}

In one of the many worlds I made a three cushion shot:
http://www.geocities.com/luktay/hit3c.gif

It really needs a 'spin about z' to make it vaguely
realistic.

T.

tailuge
05-21-2006, 11:36 PM
Jim,

Now that I have a good way to visualise the balls
rolling I've noticed that using the derivation of
throw given in your post (without second part) the
angular velocity of the balls after impact has
a component about the z axis. It all 'looks' very
convincing i.e. the balls pick up english as they
slide past and grip each other.

I've added the code but I think its correct.

I'm not sure how the equations of motion we have
used so far will behave with z component in angularv
non zero. Somehow this rotation does not decay to
zero at natural roll with the code I'm using, I think
this is the most obvious error.

T.

//
// given 2 ball events that collide at t+deltat seconds,
// calculate the velocities/spins after that impact
//
//
// based directly on post by Jim in billiards digest
//

eventListType ballEvent::newBallVelocityAtCollisionThrowModel(co nst ballEvent&amp; be1,const ballEvent&amp; be2,double deltat)
{
// only applicable to events with same starting time

traceAssert(be1.t == be2.t);

// advance both balls by same amount to point of impact

ballEvent r1= be1.progressEventToTime(be1.t+deltat);
ballEvent r2= be2.progressEventToTime(be2.t+deltat);

// to be consistent with Jim's posting extract vectors as follows

vector Pc = r1.pos;
vector Po = r2.pos;
vector Vc = r1.vel;
vector Vo = r2.vel;
vector Wc = r1.angularv;
vector Wo = r2.angularv;
vector z = vector(0,0,-1);

// The radius vector from the center of the cueball to the contact point is then

vector R = (Po - Pc)*(1.0/2.0);
R.normalize();

// switch to a frame of reference moving along with the object ball
// The cueball's velocity in this frame is

vector V = Vc - Vo;
vector v = V/V.length();

double sinTheta = (v ^ r)*z;
double cosTheta = (v * r);

tracePrint(1,"sinTheta=%f cosTheta=%f\n",sinTheta,cosTheta);

// The relative surface speed at the contact point is:

vector Vs = V + (Wc ^ R) + (Wo ^ R);

// Their relative surface speed in the tangential (friction) plane is

vector Vsf = r ^ Vs ^ r;
vector vsf = Vsf/Vsf.length();
if (Vsf.length()==0) vsf = vector(0,0,0);

tracePrint(1,"Vs =%s\n",prettyPrint(Vs).c_str());
tracePrint(1,"Vsf=%s\n",prettyPrint(Vsf).c_str());
tracePrint(1,"vsf=%s\n",prettyPrint(vsf).c_str());

// The horizontal tangent line unit vector in the friction plane
// and the cosine of the angle the relative surface speed makes with it are

vector t = r ^ z;
double cosPhi = vsf * t;

tracePrint(1,"t=%s\n",prettyPrint(t).c_str());
tracePrint(1,"cosPhi=%f\n",cosPhi);

// u is the coefficient of friction
// This is essentially Dr. Dave's treatment in TP A.14. In it, he calculates u according to:
// where exp() is the function taking the base of the natural logarithm to the power of the argument.
// The constants a,b,c are based on Wayland Marlow's data.

double a = 9.951E-3;
double b = 0.108;
double c = 1.088;
double u = a + b*(exp(-c*Vsf.length()));

// The balls may slide during the entire impact period or end up rolling across each other
// (after first sliding).

double q = 1;
if (u*V.length()*cosTheta &gt;= (1.0/7.0)*Vsf.length())
q = (1.0/7.0)*Vsf.length()/(u*V.length()*cosTheta);

// Tangent of throw angle

double tanGamma = u*q*cosPhi;

tracePrint(1,"tanGamma=%f u=%f q=%f\n",tanGamma,u,q);

// Let e stand for the coefficient of restitution between them

double e = 0.95;

// The object ball's post impact speed along the line of centers (r) in our moving frame is

double Vor = (1.0/2.0)*(1.0+e)*V.length()*cosTheta;

// The cueball will retain a small amount of speed in this direction as well, which is

double Vcr = (1.0/2.0)*(1.0-e)*V.length()*cosTheta;

tracePrint(1,"Vor=%f Vcr=%f\n",Vor,Vcr);

// Their post impact states in the stationary (table) frame of reference are then

vector newVo = Vo + (r + t * tanGamma)*Vor;
vector newVc = Vo + r*Vcr + t*(V.length()*sinTheta - Vor*tanGamma);
vector newWo = Wo - (r ^ vsf)*((5.0/2.0)*(1.0/bconst::radius)*(u*q)*(Vor));
vector newWc = Wc - (r ^ vsf)*((5.0/2.0)*(1.0/bconst::radius)*(u*q)*(Vor));

r1.vel = newVc;
r2.vel = newVo;
r1.angularv = newWc;
r2.angularv = newWo;

r1.desc = std::string("sliding");
r2.desc = std::string("sliding");

eventListType res;
res.push_back(r1);
res.push_back(r2);

return res;
}

Jal
05-22-2006, 02:21 PM
<blockquote><font class="small">Quote Jal:</font><hr>...
The third term is the centripetal acceleration and we can use one of these expressions:

W X dRz/dt = W X Vs = -(|Vs|^2/|Rz|)rz

Combining the above we have for the acceleration of the point:

As = A + (5/2)(1/r)(rb X A)X Rz - (|Vs|^2/|Rz|)rz
<hr /></blockquote>A correction just for the record. Vs should be replaced with Vr in these particular equations, where:

Vr = W X Rz or W X Rs (equivalent)

So they should be:

W X dRz/dt = W X Vr = -(|Vr|^2/|Rz|)rz and

As = A + (5/2)(1/r)(rb X A)X Rz - (|Vr|^2/|Rz|)rz

Jim

Jal
05-22-2006, 04:37 PM
<blockquote><font class="small">Quote tailuge:</font><hr> Jim,

Currently I have a very simple bit of code
to maintain the rotation. It calculates how much
the ball rotates per video frame (from angularv)
and accumulates that rotation across frames.

I didnt implement your suggestion yet, because I
also think that there must be a better way. If
linear position can be calculated (pos(t)) simply from
the equations of motion then surely angular position
must have an analagous form. (But there my brain
quits)... <hr /></blockquote>Hi Tailuge,

Since this is essentially a cosmetic feature, perhaps your current way of tracking the angular position is the way to go. I've been working on another way. It involves the series for the position of a point on the surface:

Ps' = Ps + Vst + (1/2)Ast^2 + (1/(3*2))As't^3 + (1/(4*3*2))As''t^4 + ...

This is a Taylor series where the As factors are succesive derivatives with respect to time. The derivatives follow a pattern and each succeeding one is derivable from the preceeding ones, more or less.

When you observe it in action this function takes you on a wild ride. The components of the As''s commonly reach values of 10^200+. Yet, when multiplied by the t^n factors and divided by the every growing factorials, the thing does converge to apparently correct numbers. But to have it do so in a reasonable number of terms you have restrict the time interval t to a modest value of around .1 second.

For instance, if you try to figure a point's position at 1/2 second ahead with a fast spinning ball, and carry out the computation to 100 terms, the answer will be off by a mere 10^62 inches (or something like that). This is an unfathomably large number times the size of the observable universe! I guess if you're a stickler this might cause you a problem. It's not that it's wrong, you just need to continue the calculation to more terms. Eventually, it will converge. But if you lower your goal to .1 sec ahead, it can do it in about 80 terms.

In order to get a little more of a feel for how it performs, I tried shot speeds from 1 to 20 mph with full sidespin (rwz = v), and noted how long it took to complete the computations for a certain amount of shots. This was done with the goal of having it converge to within .0001" of the right number. (It doesn't take too much more to get it much closer than this, but I thought this might be a practical value.) On my machine (1.5 ghz slowed to about 1 ghz) it went as follows. The first column is the future time at which the new position is calculated. The second is the number of shots/sec the computer was able to complete.

.1, 2222 = .0045 real/virtual secs, (40 terms avg, 82 terms max)

.05, 8000 = .0025 r/v secs, (20.2/46)

.01, 40000 = .0025 r/v secs, (6.4/14)

.001, 66,667 = .015 r/v secs, (2.45/3)

Judging by the real vs virtual ratio, it looks like the optimal look ahead time is somewhere around .01 to .05 seconds, which is probably roughly your frame rate?

If you're interested and want to know how to calculate the acceleration factors, let me know. There may very well be a series which converges more quickly. Perhap's some math/modeling sophisticate will offer help. As I said, I can see why you might just want to stick with the way you're doing it now. But at least the above has the virtue of giving you any desired precision if you go out to enough terms and/or restrict the time interval sufficiently.

Congratulations on making the shot and I see you have some shading in there. It's looking good and thanks for the image. I'll address the z-axis thing in the reply to your next post (later).

Jim

Jal
05-22-2006, 10:10 PM
<blockquote><font class="small">Quote tailuge:</font><hr> Now that I have a good way to visualise the balls
rolling I've noticed that using the derivation of
throw given in your post (without second part) the
angular velocity of the balls after impact has
a component about the z axis. It all 'looks' very
convincing i.e. the balls pick up english as they
slide past and grip each other.
Was this your expectation? <hr /></blockquote>Definitely. If one ball has no sidespin to begin with, it will generally acquire some during the collision. If it has some sidespin, it may gain or lose some during the collision. If the other ball is at rest, it will almost always lose some of it.

<blockquote><font class="small">Quote tailuge:</font><hr>I've added the code but I think its correct.

I'm not sure how the equations of motion we have
used so far will behave with z component in angularv
non zero. Somehow this rotation does not decay to
zero at natural roll with the code I'm using, I think
this is the most obvious error.
<hr /></blockquote>It's not so much an error as something we haven't gotten to yet. The equations of motion predict the same trajectories with or without sidespin in the pre and post impact phases (we're ignoring a very small cloth effect discussed early on in this thread). This is because the torque applied to the ball by cloth friction is perpendicular to the z-axis. Any spin that developes is in the direction of the torque and that's also why it cannot lose sidespin (with the way things stand right now). The developing spin cannot cancel it out.

This seems to suggest that swerve will not occur, but not really. Swerve is due to a horizontal spin component.

To implement the decay, you can add a z-axis component to the angular acceleration. Some experimentation is called for but here is a guesstimate.

If you hit a ball hard with max sidespin head on into another ball, on pool cloth it'll spin for what, 10-20 seconds? Let's say 15 secs and that the cueball is traveling at 12 mph with a spin speed ratio of one (rwz/v = 1). That gives:

rwz = (12 mph) (1.467 ft/sec/mph)(12 in/ft) = 211 in/sec

Wz = 211/r = 211/1.125 = 188 radians/sec

So the angular acceleration is:

Awz = -188/15 sec = -12.5 radians/sec/sec

To implement it just multiply the unit vector in the direction of the z-spin component, which may be 1 or -1 in length, by this.

I looked over your code and everything seems to be okay. Onward and maybe even upward!

Jim

tailuge
05-23-2006, 04:51 AM
Jim,

I'll try the Z-axis angular acceleration,
guess this means I'll need another 'state'

sliding-&gt;rolling-&gt;spinning-&gt;atrest

and a simple function timeToStopSpinning.
As you know I'm excited to see the side
'bite' into the cushion soon.

As for the dot on the ball, I'll stick with my current
code (which is in essence an integration per frame)
just because its simple. I do however think there must
be a mapping (Wx,Wy,Wz) -&gt; local(x,y,z) that we are missing,
the equations of motion for
linear and angular systems are supposed to be analagous.
half way down the page.

T.

Jal
05-24-2006, 03:51 PM
Hi Tailuge,

The equations of motion are analagous to an extent. But unlike with linear motion, a point on the surface of a uniformly rotating sphere is undergoing a continuous acceleration in the direction of the sphere's center (centripetal). The direction of this acceleration is continuously changing, which complicates things further over the translational case. When you add a torque that changes the rotation, things get more complicated. Instead of just a few non-zero derivatives, you get an infinite series of them. For the case of constant rotation, all of these terms add up to a simple expression involving cos(wt) and sin(wt). When you have a torque acting on the sphere, extra terms are generated. It may be that they also add up to combinations of cos(wt) and sin(wt), or something like that, but with early terms of their series missing. I don't know, I'm looking at it. Since this could take a while, let's get on with the simpler treatment of the cushion. A number of things are ignored, but here it is.

We'll treat it just like we did throw. The cushion is the "object ball" and it has a position vector at the point of contact, Po. Using this, you do the same things to get the normal vector R, the "cut angle" theta, phi, etc.. Of course, it goes withou saying that to figure the surface speed Vsf, you only use the ball's speed and spin.

Let rho be the angle of inclination the nose of the cushion is above the center of the "cueball". On a pool table this is about 15.67 degrees. When you check to see if the ball will end up rolling across the cushion, you compare (u/cos(rho))(V)(1+e)cos(theta) to (2/7)Vsf, where u is approx. .2, and e is about .7 or .8 (coefficient of restitution of the cushion). Then either q=1, or q=(2/7)Vsf divided by the above expression, as before. (We use 2/7 here instead of 1/7 because the cushion can't move or spin to reduce the relative surface speed.)

The expression for Vcr when figuring Vc' is -[e/cos(rho)]Vcos(theta), to yield the post impact velocity in the normal (r) direction. In the parallel (tangential) direction, you use Vsin(theta)-[(1+e)/cos(rho)]Vcos(theta)tan(gamma). The (1+e) factor comes from the total impulse applied to the ball in the normal direction and we don't divide by 2 here.

The change in its spin is -(5/2)(1/r)(uq)[(1+e)/cos(rho)]Vcos(theta)rc X vsf, where rc is the unit vector pointing to the cushion nose and thus inclined upward from r by the angle rho. Its horizontal component is (rcos(rho))r and its vertical component is (rsin(rho))z.

The above ignores bed friction and the fact the ball can't move in the z-direction to reduce surface speed when figuring its "throw". We're also ignoring any cushion throwback (asymmetrical compression forces) and torsional friction with the cushion. I'm not sure that cushion throwback is a serious effect at this point. My more complete treatment also ignores these last two things and mainly deals with bed friction. This is a relatively small but not completely insignificant effect.

Sorry for the wait. I know you've been anxious to add this feature but the other problem now has my attention. Please don't ask what happens when a ball strikes a corner right now. /ccboard/images/graemlins/confused.gif

Jim

tailuge
05-25-2006, 01:21 AM
Jim,

Thanks once again. I've attached the code based on your
post here. It looks very convincing in action. I updated
my controls to allow me to strike the ball with lhs/rhs
aswell as follow/draw. I never expected this post would
lead me this far - its great fun. I'm adding some diamonds
on the rails to see if we can match one of the billiard
diamond systems (you know StartDiamond-First=Second ?)

To this end I might try and make the cuetip/ball interaction
a little more realistic, i.e. max spin/vel ratio, I'll post
and you can give it the once over.

As for the decay in Z axis spin, I have a slight problem
because zero Z spin might occur durring slide,roll or while
stationary. Thats something I'll have to think about in terms of my code.

T.

//
// CUSHION IMPACT
// given ball at impact with cushion
// calculate the velocity/spin after that impact
//
// cushPos is the relative position of the cushion to the ball
// at impact i.e. (1,0) rightcush.
//
// again based on post by Jim in billiards digest
//

ballEvent ballEvent::reflectBallInCushionThrow(const vector cushPos)
{

ballEvent r1= *this;

// to be consistent with Jim's posting extract vectors as follows
// Object ball placed to act as cushion

vector Pc = r1.pos;
vector Po = Pc + (cushPos*(bconst::radius+bconst::minsep));
vector Vc = r1.vel;
vector Wc = r1.angularv;
vector z = vector(0,0,-1);

// The radius vector from the center of the cueball to the contact point is then

vector R = (Po - Pc)*(1.0/2.0);
R.normalize();

vector V = Vc;
vector v = V/V.length();

double sinTheta = (v ^ r)*z;
double cosTheta = (v * r);

// The relative surface speed at the contact point is:

vector Vs = V + (Wc ^ R) + (Wo ^ R);

// Their relative surface speed in the tangential (friction) plane is

vector Vsf = r ^ Vs ^ r;
vector vsf = Vsf/Vsf.length();
if (Vsf.length()==0) vsf = vector(0,0,0);

// Let rho be the angle of inclination the nose of the cushion
// is above the center of the "cueball".

const double Rho = 15.67 * (bconst::pi / 180.0);
const double cosRho = cos(Rho);
const double sinRho = sin(Rho);

// The horizontal tangent line unit vector in the friction plane
// and the cosine of the angle the relative surface speed makes with it are

vector t = r ^ z;
double cosPhi = vsf * t;

// When you check to see if the ball will end up rolling across the cushion,
// you compare (u/cos(rho))(V)(1+e)cos(theta) to (2/7)Vsf,
// where u is approx. .2, and e is about .7 or .8 (coefficient of restitution of the cushion).
// Then either q=1, or q=(2/7)Vsf divided by the above expression, as before.
// (We use 2/7 here instead of 1/7 because the cushion can't
// move or spin to reduce the relative surface speed.)

double e = 0.75;
double u = 0.2;

double q = 1;
if ((u/cosRho)*V.length()*(1+e)*cosTheta &gt;= (2.0/7.0)*Vsf.length())
q = (2.0/7.0)*Vsf.length()/((u/cosRho)*V.length()*(1+e)*cosTheta);

// Tangent of throw angle

double tanGamma = u*q*cosPhi;

// Change in normal and tangential velocity

double Vcr = (-e/cosRho)*V.length()*cosTheta;
double Vct = V.length()*sinTheta - ((1+e)/cosRho)*V.length()*cosTheta*tanGamma;

// rc is the unit vector pointing to the cushion nose and thus inclined upward from r by the angle rho.
// Its horizontal component is (rcos(rho))r and its vertical component is (rsin(rho))z.

vector rc = r*cosRho + z*sinRho;

// Change in angular velocity
// The change in its spin is -(5/2)(1/r)(uq)[(1+e)/cos(rho)]Vcos(theta)rc X vsf

vector deltaWc = (rc ^ vsf) * (-5.0/2.0)*((1.0/bconst::radius)*u*q*((1+e)/cosRho)*V.length()*cosTheta);

// Post impact state

vector newVc = r*Vcr + t*Vct;
vector newWc = Wc + deltaWc;

r1.vel = newVc;
r1.angularv = newWc;
r1.desc = std::string("sliding");

return r1;
}

Jal
05-25-2006, 03:00 PM
Tailuge,

Thanks for the code. If I get up the gumption to learn your language, these examples will be an excellent aid.

I'm really amazed that you have balls moving around on screen and actually bouncing off of things so soon. When you first posted I assumed that this would be several months down the road, if you persued it at all.

Jim

tailuge
05-26-2006, 02:16 AM
Note:

In my posting of code:

vector z = vector(0,0,-1);
should probably be
vector z = vector(0,0,1);
But I cant actualy see any difference.

vector Vs = V + (Wc ^ R) + (Wo ^ R);
Should be:
vector Vs = V + (Wc ^ R);

If anyone out there is reading besides Jim and myself then
my program works on a basic data structure of a 'ball event':

class ballEvent
{
public:

int ballid;
double t;

// physical states

vector pos;
vector vel;
vector angularv;
vector angularpos;

std::string desc;
}

Each time a ball changes from 'sliding' to 'rolling'
to 'stationary' a new event is created in a list
of events. A ball can also obtain new event through
collision with cushion or other ball.

All these ball events are maintained in a list
in chronological order by the 'event manager'
this looks at each possible future event that
might occur for all the balls and chooses the
one that happens first and adds that to the list.
By repeating this step all state changes for a
shot can be calculated before any drawing on
screen occurs. The state expansion code is fairly
quick so that many possible outcomes can be compared
based on some variation in inputs.

class eventManager
{
public:

bool expandByOneEvent();
int expandAllEvents();

private:

eventListType nextEventsFromThisEvent(ballEvent be);

eventListType eventList;
}

A typicial list of expanded events looks like:
(ball id, time, pos,vel,angularvel,angularpos,state)

ballEvent be(0,0.000000,vector(0.000000,0.000000,0.000000), vector(0.142111,0.595493,0.000000), vector(0.000000,0.000000,0.000000),"sliding");
ballEvent be(1,0.000000,vector(0.000000,0.104850,0.000000), vector(0.000000,0.000000,0.000000), vector(0.000000,0.000000,0.000000),"atrest");
ballEvent be(2,0.000000,vector(0.000000,0.209700,0.000000), vector(0.000000,0.000000,0.000000), vector(0.000000,0.000000,0.000000),"atrest");
ballEvent be(0,0.063057,vector(0.008459,0.035444,0.000000), vector(0.160116,0.032347,0.000000), vector(-3.185792,1.334032,1.885110),"sliding");
ballEvent be(1,0.063057,vector(0.000000,0.104850,0.000000), vector(-0.033943,0.496360,0.000000), vector(1.591503,0.193959,1.885110),"sliding");
ballEvent be(0,0.099340,vector(0.013680,0.037027,0.000000), vector(0.127690,0.054917,0.000000), vector(-1.571310,3.653500,1.885110),"rolling");
ballEvent be(1,0.139986,vector(-0.002374,0.139821,0.000000), vector(-0.026420,0.011195,0.000000), vector(-2.862222,-0.298602,0.812407),"sliding");
ballEvent be(2,0.139986,vector(0.000000,0.209700,0.000000), vector(-0.001360,0.401625,0.000000), vector(1.521956,-0.051709,-1.072703),"sliding");
ballEvent be(1,0.163671,vector(-0.002946,0.140387,0.000000), vector(-0.021854,0.036577,0.000000), vector(-1.046565,-0.625280,0.812407),"rolling");
ballEvent be(2,0.259326,vector(-0.000170,0.249876,0.000000), vector(-0.001488,0.271677,0.000000), vector(-7.773316,-0.042562,-1.072703),"rolling");
ballEvent be(1,0.859320,vector(-0.010547,0.153109,0.000000), vector(0.000000,-0.000000,0.000000), vector(0.000000,-0.000000,0.812407),"atrest");
ballEvent be(0,2.368704,vector(0.158567,0.099341,0.000000), vector(0.000000,-0.000000,0.000000), vector(0.000000,0.000000,1.885110),"atrest");
ballEvent be(2,4.694942,vector(-0.003469,0.852404,0.000000), vector(0.000000,-0.000000,0.000000), vector(0.000000,0.000000,-1.072703),"atrest");

To visualise the shot before hitting, I draw straight lines
between the positions described by the events, but when
I animate the balls I interpolate between them using the
equations of motion described by the active events, splendid
curves abound.

Besides some 'refining' of roots there would be no
heavy iteration in this code and it could easily be
used in less efficient languages like flash/javascript.

T.

tailuge
05-27-2006, 09:53 PM
Without much thought I scanned the Shepard paper
and found this formula for spin ratio in terms
of offset of cuetip from dead centre. (level cue)

// p = power (0-1.0)
// dir = cue aim (x,y,0) (flat cue)
// facex,facey impact of face as viewed from cue (0-1.0)

void cue::setHitVel(double p)
{
hitvel = dir*(p*bconst::maxvel);
double bz = facex * bconst::maxcueoffset;
double by = facey * bconst::maxcueoffset;
hitangularv = (dir ^ bconst::vert)*Wy + vector(0,0,Wz);
}

With the model we currently have I think we could
actualy hit masse shots. If anyone has a formula
that can take a cue direction vector(x,y,z), cue pos
(any point in the cue) and ball position and can
tell me the balls initial spin I'll try it out.

T.

dr_dave
05-28-2006, 11:02 AM
<blockquote><font class="small">Quote tailuge:</font><hr>If anyone has a formula
that can take a cue direction vector(x,y,z), cue pos
(any point in the cue) and ball position and can
tell me the balls initial spin I'll try it out.<hr /></blockquote>
You should be able to get what you need from TP A.19 (http://www.engr.colostate.edu/~dga/pool/technical_proofs/new/TP_A-19.pdf). It is a complete analysis of masse shot physics and aiming methods.

Regards,
Dave

Jal
05-30-2006, 12:44 AM
<blockquote><font class="small">Quote tailuge:</font><hr> Without much thought I scanned the Shepard paper
and found this formula for spin ratio in terms
of offset of cuetip from dead centre. (level cue)

// p = power (0-1.0)
// dir = cue aim (x,y,0) (flat cue)
// facex,facey impact of face as viewed from cue (0-1.0)

void cue::setHitVel(double p)
{
hitvel = dir*(p*bconst::maxvel);
double bz = facex * bconst::maxcueoffset;
double by = facey * bconst::maxcueoffset;
hitangularv = (dir ^ bconst::vert)*Wy + vector(0,0,Wz);
}

With the model we currently have I think we could
actualy hit masse shots. If anyone has a formula
that can take a cue direction vector(x,y,z), cue pos
(any point in the cue) and ball position and can
tell me the balls initial spin I'll try it out.

T.
<hr /></blockquote>Tailuge,

If you haven't found what you're looking for, you can use this. (I hate to say it but the above code will not work.)

Let V be the speed of the cue stick with associated unit vector v.

Form the additional unit vectors, where z is the usual vertical unit vector:

z' = v X z X v/|v X z X v|

t' = v X z'

n' = z' X t'

n = z X t'

The vectors z' and t' define a plane perpendicular to the stick velocity V, with n' perpendicular to it and in the direction of V (ie, n' = v).

Suppose you want the tip offset to be Oe inches, centimeters, or whatever units you're using, off to the side, and Of inches above center, as you see it looking down the shaft of the stick. To be consistent with the vectors just defined, make Oe positive for right english, negative for left, and Of positive for follow, negative for draw.

The offset vector is then:

O = (Oe/R)t' + (Of/R)z'

Now let u = (5/2)|O|^2+1 and Mb/Ms be the ratio of the masses (or weights) of the ball to stick (about 1/3). Then the velocity of the cueball parallel to the table bed will be:

Vc = k[2/(u + Mb/Ms)](V * n)n

where k is a coefficient related to the elasticity of the collision (see below). Its spin will be:

Wc = (5/2)(1/R)O X Vc

Note that it's Vc and not V in the expression. The above should work for any orientation of V (ie, elevated cue).

You can set k=1 or use the following as you wish:

k = (1/2)[(2 - u + Mb/Ms)e + u + Mb/Ms]/[1 + Mb/Ms]

where u is as defined above and e is approximately .75. K has the interesting property of being equal to one at a certain tip offset (approx. .73R), no matter what the coefficient of restitution e happens to be. In other words, the collision behaves as if it were perfectly elastic even if it's perfectly inelastic! Sometimes if it looks like a duck, walks like a duck, and quacks like a duck, it's a caterpillar.

To get the position of the tip, actually the contact point, at the momemt of collision, set:

Ptn' = - Sqrt(R^2 - Oe^2 - Of^2)

Note the negative sign. Then its position will be:

Pt = P + O + (Ptn')n'

where P is the position of the center of the cueball.

Coriolis' amazing insight on how to determine the direction of the cueball at natural roll can be checked. Or perhaps more to the point, you can verify that your program is working properly by seeing if the cueball conforms to his prediction. To do this let Ptz be the vertical component of Pt, and vz the vertical component of v. Then let:

a = |(-R - Ptz)/vz| (absolute value is indicated), and

Pj = Pt + av

Pj is the position of the contact point projected onto the table bed in the direction V (if I did the math right). The vector from the cueball's contact point with the bed to Pj is:

C = Pj - (P - Rz)

The cueball at natural roll should run parallel to C. Although Coriolis' proposition applies to any shot, not just masses, you can try this only when there is considerable cue elevation because of the projected contact point's position.

Jim

Jal
05-30-2006, 02:03 PM
<blockquote><font class="small">Quote Jal:</font><hr>...Then the velocity of the cueball parallel to the table bed will be:

Vc = k[2/(u + Mb/Ms)](V * n)n

where k is a coefficient related to the elasticity of the collision (see below). Its spin will be:

Wc = (5/2)(1/R)O X Vc<hr /></blockquote>Whoops. The expression for the spin should use the cueball's velocity in the direction of the cue stick, as if the table bed wasn't there. So let:

Vcv = k[2/(u + Mb/Ms)]V

Then:

Wc = (5/2)(1/R)O X Vcv

Jim