Hi all!

I'm adding this post for those who are curious about how I handled the motion of the pendulum in the 1st video. It goes to the limits of my abilities in several ways.

I had an opportunity while working 3rd shift for a couple of years to finally

take on that long overdue review of classical mechanics. I had (typically) 6 hours per shift of free time, and it was nice being paid $22/hour to do it.

I obtained and studied from: Introduction to Classical Mechanics by David Morin. 1st edition (in case there's a 2nd by now).

My educational background: BS in physics long ago.

I'm using skills which would have faded to nothing by now if not for my usage in this hobby.

My goal was to make it through chapter 6 on the Lagrangian method. It's an alternative formulation (to vector based Newtonian mechanics) based on scalar quantities (just a number) such as kinetic and potential energy and is well suited to problems where motion is limited by a constraint, such as being limited to motion on the surface of a sphere in 3d space like this pendulum.

I highly recommend the textbook. It is as hardcore as anyone could want and it's loaded with fully worked solutions for a portion of the very full range of problems in each chapter. It's also loaded with clever little physics limerics. I'll post 1 or 2 here later.

This text has a place of honor in my work here. It can be seen floating in midair in the room in the 2nd video.

Please find attached type persPendulum. It's virtual update() has a bit more than usual going on, which gets right away to an issue with the enormous time step between frames, and to why type vec3d became involved.

My simplistic numerical integration scheme is this:

step 1: obtain the acceleration to be applied this frame.

In the case of the shots the acceleration is constant = (0,g,0).

In the case of the pendulum it's far from that.

step 2:

update the position:

position += velocity*dt + 0.5f*acceleration*dt*dt;

This is from Taylors series to 2 terms.

Then update the velocity

velocity += acceleration*dt;

The analysis will be in spherical coordinates.

So far these methods are persPendulum methods only, but I think they also

should be static methods available in vec3d?

// base vectors in spherical coordinates

vec3d get_rHat()const{ return sin(th)*( cos(phi)*Xu + sin(phi)*Yu ) - cos(th)*Zu; }

vec3d get_thHat()const{ return cos(th)*( cos(phi)*Xu + sin(phi)*Yu ) + sin(th)*Zu; }

vec3d get_phiHat()const{ return cos(phi)*Yu - sin(phi)*Xu; }

Where the vectors Xu, Yu, Zu are the cartesian basis.

I'll also attach my solution on paper soon, but I want it to be much cleaner than the working version I have.

The identifiers for the two angles theta (polar) and phi (azimuthal) and their derivatives are:

// polar angle and 1st and 2nd time derivatives

double th, th_1, th_2;

// azimuthal angle and 1st and 2nd time derivatives

double phi, phi_1, phi_2;

The 2 equations of motion obtained by the Lagrangian method are:

// step 1: find acceleration

th_2 = phi_1*phi_1*sinTh*cosTh - GravY*sinTh/L;// for theta

phi_2 = -2.0*th_1*phi_1*cosTh/sinTh;// for phi

Here, sinTh, cosTh are saved values of sin(th), cos(th), GravY = gravity.y and L = length of support. Note that the mass isn't a factor.

An additional condition may also be applied, perhaps to stabilize the integration?

Because the Lagrangian = kinetic - potential energy doesn't depend explicitly on the angle phi, there is a constant of the motion. Here it's the z component of the angular momentum so:

phi_1 = Lz/( sin(th)*sin(th) );// based on Lz conservation

where the value of Lz is determined by initial conditions. The quantity here has been scaled from the actual Lz: Scaled this is actually Lz /= m*L*L.

So the units of measure ought to balance above.

Integration at the frame rate ( dt = 20ms ) does not give a stable result.

It became stable at 10 updates per frame, so I tripled that and called it good.

The pendulum in the video is being updated 30 times per frame.

Once all of these micro updates are done I assign bobs position and velocity.

I'm still struggling with the float/double interface here.

// bobs new velocity

vec3d thHat = cosTh*( cosPhi*Xu + sinPhi*Yu ) + sinTh*Zu;// current thHat

vec3d phiHat = cosPhi*Yu - sinPhi*Xu;// current phiHat

Vel = L*( th_1*thHat + sinTh*phi_1*phiHat );

bob.vel = Vel.to_vec3f();// bob.vel update

vec3d rHat = sinTh*( cosPhi*Xu + sinPhi*Yu ) - cosTh*Zu;// updated rHat

// and position

vec3d Pos = L*rHat;

Pos.x += static_cast<double>(topQ.pos.x);

Pos.y += static_cast<double>(topQ.pos.y);

Pos.z += static_cast<double>(topQ.pos.z);

bob.setPosition( Pos.to_vec3f() );// bob update

I'll leave off here for now. But I must deal with the sudden change in th_1 and phi_1 which occur when poor bob gets shot.

I'm back to wrap up.

As for handling bobs hits:

I'm cheating a bit. I call hitFree() and do not account for the constraint in the collision,

so the momentum transferred isn't accurate, but it seems to do fine for effect. This serves the present purpose well enough for me. Next problem please.

Update the values for th_1 and phi_1

bobs velocity in polar coordinates is:

vec3d persPendulum::get_vel()const

{ return vec3d( L*( th_1*get_thHat() + sin(th)*phi_1*get_phiHat() ) ); }

Let vec3d Vel = bob.vel following the collision, then

Taking the dot product (of the return value above) with thHat (unit vector in theta direction) gives

Vel.dot( thHat ) = L*th_1

solving for th_1 permits the update

th_1 = Vel.dot(thHat)/L;

Take the dot product with phiHat instead to obtain

phi_1 = Vel.dot( phiHat )/( L*sin(th) );

following this I update Vel and then bob.vel

Vel = get_vel();

bob.vel = Vel.to_vec3f();

Just call

bool hit( persBall& PB, float cR, float dt );

where Cr = coefficient of restitution allows partially elastic collisions.

Why dt? It won't be moved here, just in update(). The collision test is ray based so dt is used to obtain the position in the previous frame = 2 points for a ray.

Finally, I'll admit that the integration details are well beyond my knowledge.

I think I should be using th_1 = Lz/( sin(th)*sin(th) ); instead of the 2nd acceleration equation, which actually follows from above. I'm presently experimenting with using

phi_1 = Lz/( sinTh*sinTh );

phi += phi_1*dT;

in place of

phi_2 = -2.0*th_1*phi_1*cosTh/sinTh;

phi += phi_1*dT + 0.5*phi_2*dT*dT;

phi_1 += phi_2*dT;

and it's giving a stable result for lower # of micro updates.

However, I'm once again getting runaway cases...

It's a fight to maintain stability and I don't know enough right here.

Perhaps someone who knows this part better may advise? Please?

I've also now attached my "on paper" solution.

And that's all I've got so I'll call it a wrap.

Any questions and/or suggestions are welcome.