Visual Python

[agor95]

VPython code modified from the examples

Code: Select all

from visual import *

# pendulum

# The analysis is in terms of Lagrangian mechanics.
# The Lagrangian variables are angle of bar
# measured from the vertical.

scene.title = 'Pendulum'
scene.height = scene.width = 800

g = 9.8
M1 = 2.0
M2 = 1.0

d = 0.05 # thickness of each bar

gap = 2.*d # distance between two parts of upper, U-shaped assembly

L1 = 1.0 # physical length of upper assembly; distance between axles
L1display = L1+d # show upper assembly a bit longer than physical, to overlap axle

L2 = 1.0 # physical length of lower bar
L2display = L2+d/2. # show lower bar a bit longer than physical, to overlap axle

# Coefficients used in Lagrangian calculation
A = (1./4.)*M1*L1**2+(1./12.)*M1*L1**2+M2*L1**2
B = (1./2.)*M2*L1*L2
C = g*L1*(M1/2.+M2)
D = M2*L1*L2/2.
E = (1./12.)*M2*L2**2+(1./4.)*M2*L2**2
F = g*L2*M2/2.

hpedestal = 1*(L1+L2) # height of pedestal
wpedestal = 0.1 # width of pedestal

tbase = 0.05 # thickness of base
wbase = 8.*gap # width of base
offset = 2.*gap # from center of pedestal to center of U-shaped upper assembly

top = vector(0,0,0) # top of inner bar of U-shaped upper assembly

scene.center = top-vector(0,(L1+L2)/2.,0)

pedestal = box(pos=top-vector(0,hpedestal/2.,offset),
                 height=1.1*hpedestal, length=wpedestal, width=wpedestal,
                 color=(0.4,0.4,0.5))
base = box(pos=top-vector(0,hpedestal+tbase/2.,offset),
                 height=tbase, length=wbase, width=wbase,
                 color=pedestal.color)
axle = cylinder(pos=top-vector(0,0,gap/2.-d/4.), axis=(0,0,-offset), radius=d/4., color=color.yellow)

frame1 = frame(pos=top)
bar1 = box(frame=frame1, pos=(L1display/2.-d/2.,0,-(gap+d)/2.), size=(L1display,d,d), color=color.red)
bar1b = box(frame=frame1, pos=(L1display/2.-d/2.,0,(gap+d)/2.), size=(L1display,d,d), color=color.red)
axle1 = cylinder(frame=frame1, pos=(L1,0,-(gap+d)/2.), axis=(0,0,gap+d),
                 radius=axle.radius*20, color=axle.color)
frame1.axis = (0,-1,0)


theta1 = -2.001*pi/2. # initial upper angle (from vertical)
theta1dot = 0 # initial rate of change of theta1

frame1.rotate(axis=(0,0,1), angle=theta1)

scene.autoscale = 0

dt = 0.001
t = 0.

while 1:
    rate(1./dt)

    # Calculate accelerations of the Lagrangian coordinates:
    atheta1 = ((E*C/B)*sin(theta1))/(D-E*A/B)

    # Update velocities of the Lagrangian coordinates:
    theta1dot = theta1dot+atheta1*dt

    # Update Lagrangian coordinates:
    dtheta1 = theta1dot*dt

    theta1 = theta1+dtheta1
    
    frame1.rotate(axis=(0,0,1), angle=dtheta1)

    t = t+dt

[agor95]

This is the VPython Intro for you to see without installing the development application.

[ME]

Unfortunately I don't have vPython anymore..

But there's a slight improvement to make in the kinematics formula:

Code: Select all

    # Update Lagrangian coordinates: 
    dtheta1 = (1./2.)*atheta1*dt*dt + theta1dot*dt  

    # Update velocities of the Lagrangian coordinates:
    theta1dot = theta1dot+atheta1*dt 

    theta1 = theta1+dtheta1 

When you check this adjustment with a parabolic flight for example (constant acceleration), you'll see the velocities and positions match.

[agor95]

I appreciate your improvement to the maths.

I did have problems with the VPython of late.

I reinstalled Python-2.7.13.msi
and the

VPython-Win-32-Py2.7-6.01.exe

I have another version for Python 2.4; just in case.

[agor95]

Example of atoms in a box.

The number of atoms and their speeds are shown in a bar graph.

[agor95]

@ME

I am trying to reduce the formula from a double to a single pendulum.

Code: Select all

from __future__ import print_function, division
from visual import *
from visual.graph import *

# Double pendulum

# The analysis is in terms of Lagrangian mechanics.
# The Lagrangian variables are angle of upper bar and angle of lower bar,
# measured from the vertical.

# Bruce Sherwood
#
# Corrections to the Lagrangian calculations by Owen Long, UC. Riverside
#

scene.title = 'Double Pendulum'
scene.height = scene.width = 800

g = 9.8
M1 = 0.0000001
M2 = 2.0
L1 = 0.0000001 # physical length of upper assembly; distance between axles
L2 = 0.5 # physical length of lower bar

I1 = M1*L1**2/12 # moment of inertia of upper assembly

I2 = M2*L2**2/12 # moment of inertia of lower bar

d = 0.05 # thickness of each bar
gap = 2*d # distance between two parts of upper, U-shaped assembly
L1display = L1+d # show upper assembly a bit longer than physical, to overlap axle
L2display = L2+d/2 # show lower bar a bit longer than physical, to overlap axle

hpedestal = 1.3*(L1+L2) # height of pedestal
wpedestal = 0.1 # width of pedestal
tbase = 0.05 # thickness of base
wbase = 8*gap # width of base
offset = 2*gap # from center of pedestal to center of U-shaped upper assembly
top = vector(0,0,0) # top of inner bar of U-shaped upper assembly
scene.center = top-vector(0,(L1+L2)/2.,0)

theta1 = 0.0 # initial upper angle (from vertical)
theta2 = 2. # initial lower angle (from vertical)
theta1dot = 0 # initial rate of change of theta1
theta2dot = 0 # initial rate of change of theta2

pedestal = box(pos=top-vector(0,hpedestal/2.,offset),
                 height=1.1*hpedestal, length=wpedestal, width=wpedestal,
                 color=(0.4,0.4,0.5))
base = box(pos=top-vector(0,hpedestal+tbase/2.,offset),
                 height=tbase, length=wbase, width=wbase,
                 color=pedestal.color)
axle = cylinder(pos=top-vector(0,0,gap/2.-d/4.), axis=(0,0,-offset), radius=d/4., color=color.yellow)

frame1 = frame(pos=top)
bar1 = box(frame=frame1, pos=(L1display/2.-d/2.,0,-(gap+d)/2.), size=(L1display,d,d), color=color.red, opacity=0)
bar1b = box(frame=frame1, pos=(L1display/2.-d/2.,0,(gap+d)/2.), size=(L1display,d,d), color=color.red, opacity=0)
axle1 = cylinder(frame=frame1, pos=(L1,0,-(gap+d)/2.), axis=(0,0,gap+d),
                 radius=axle.radius, color=axle.color)
frame1.axis = (0,-1,0)
frame2 = frame(pos=frame1.axis*L1)
bar2 = box(frame=frame2, pos=(L2display/2.-d/2.,0,0), size=(L2display,d,d), color=color.green)
frame2.axis = (0,-1,0)
frame1.rotate(axis=(0,0,1), angle=theta1)
frame2.rotate(axis=(0,0,1), angle=theta2)
frame2.pos = top+frame1.axis*L1

dt = 0.00005
t = 0.

C11 = (0.25*M1+M2)*L1**2+I1
C22 = 0.25*M2*L2**2+I2

#### For energy check:
gdisplay(x=800)
gK = gcurve(color=color.yellow)
gU = gcurve(color=color.cyan)
gE = gcurve(color=color.red)

while t < 5&#58;
    rate&#40;1/dt&#41;
    
    # Calculate accelerations of the Lagrangian coordinates&#58;
    C12 = C21 = 0.5*M2*L1*L2*cos&#40;theta1-theta2&#41;
    Cdet = C11*C22-C12*C21
    a = .5*M2*L1*L2*sin&#40;theta1-theta2&#41;
    A = -&#40;.5*M1+M2&#41;*g*L1*sin&#40;theta1&#41;-a*theta2dot**2
    B = -.5*M2*g*L2*sin&#40;theta2&#41;+a*theta1dot**2
    atheta1 = &#40;C22*A-C12*B&#41;/Cdet
    atheta2 = &#40;-C21*A+C11*B&#41;/Cdet

    # Update velocities of the Lagrangian coordinates&#58;
    theta1dot += atheta1*dt
    theta2dot += atheta2*dt
    
    # Update Lagrangian coordinates&#58;
    dtheta1 = theta1dot*dt
    dtheta2 = theta2dot*dt
    theta1 += dtheta1
    theta2 += dtheta2

#    frame1.rotate&#40;axis=&#40;0,0,1&#41;, angle=dtheta1&#41;
#    frame2.pos = top+frame1.axis*L1
#    frame2.pos = top
    frame2.rotate&#40;axis=&#40;0,0,1&#41;, angle=dtheta2&#41;
    t = t+dt

#### Energy check&#58;
    K = .5*&#40;&#40;.25*M1+M2&#41;*L1**2+I1&#41;*theta1dot**2+.5*&#40;.25*M2*L2**2+I2&#41;*theta2dot**2+\
        .5*M2*L1*L2*cos&#40;theta1-theta2&#41;*theta1dot*theta2dot
    U = -&#40;.5*M1+M2&#41;*g*L1*cos&#40;theta1&#41;-.5*M2*g*L2*cos&#40;theta2&#41;
    gK.plot&#40;pos=&#40;t,K&#41;&#41;
    gU.plot&#40;pos=&#40;t,U&#41;&#41;
    gE.plot&#40;pos=&#40;t,K+U&#41;&#41;

This section is the hardest to untangle.
M1 & L1 are very small number

Code: Select all

    # Calculate accelerations of the Lagrangian coordinates:
    C12 = C21 = 0.5*M2*L1*L2*cos(theta1-theta2)
    Cdet = C11*C22-C12*C21
    a = .5*M2*L1*L2*sin(theta1-theta2)
    A = -(.5*M1+M2)*g*L1*sin(theta1)-a*theta2dot**2
    B = -.5*M2*g*L2*sin(theta2)+a*theta1dot**2
    atheta1 = (C22*A-C12*B)/Cdet
    atheta2 = (-C21*A+C11*B)/Cdet

[ME]

I can imagine... With Python it's possible to write as bad as in any other language.

Just replace all M1, and L1 with zero's, replace all multiplication with zero by a zero and then remove all variables containing a zero... I guess you tried that already as they (L1 and M1) are already very small.
The obvious problem is that (C12=C21=0) because it contains L1.
(C11) is also zero, also because of the same reason.
Hence (cdet) is zero, thus a division by zero....

So now you need to extract the formula's to understand why it calculates it this way.
That takes some time.

[agor95]

@ME

Thanks got there in the end comment out and revers if error.
I wish I knew what I was doing :-)

It was like burning a good book.

Code: Select all

    # Calculate accelerations of the Lagrangian coordinates:
    B = -.5*M2*g*L2*sin(theta2)
    atheta2 = B/C22
    # Update velocities of the Lagrangian coordinates:
    theta2dot += atheta2*dt
    
    # Update Lagrangian coordinates:
    dtheta2 = theta2dot*dt
    theta2 += dtheta2

At least this is understandable

[agor95]

Code: Select all

from __future__ import print_function, division
from visual import *
from visual.graph import *

# pendulum

# The analysis is in terms of Lagrangian mechanics.
# The Lagrangian variables are angle of bar
# measured from the vertical.

scene.title = 'Pendulum'
scene.height = scene.width = 800

g = 9.8          # meters per second  ** 2
M1 = 1.0         # 1 kg

L1 = 0.9         # physical length of assembly; distance between axles

d = 0.05         # thickness of each bar
gap = 2.*d       # distance between two parts, U-shaped assembly

L1display = (2.0-L1)+d # show assembly a bit longer than physical, to overlap axle
L2display = (L1)+d     # show assembly a bit longer than physical, to overlap axle

hpedestal = 1*L1+0.5 # height of pedestal
wpedestal = 0.1      # width of pedestal

tbase = 0.05    # thickness of base
wbase = 8.*gap  # width of base
offset = 2.*gap # from center of pedestal to center of U-shaped upper assembly

top = vector(0,0,0) # top of inner bar of U-shaped upper assembly

scene.center = top-vector(0,(L1)/2.,0)

pedestal = box(pos=top-vector(0,hpedestal/2.,offset),
                 height=1.1*hpedestal, length=wpedestal, width=wpedestal,
                 color=(0.4,0.4,0.5))
base = box(pos=top-vector(0,hpedestal+tbase/2.,offset),
                 height=tbase, length=wbase, width=wbase,
                 color=pedestal.color)
axle = cylinder(pos=top-vector(0,0,gap/2.-d/4.), axis=(0,0,-offset), radius=d/4., color=color.yellow)

frame1 = frame(pos=top)
bar1 = box(frame=frame1, pos=(L1display/2.-d/2.,0,-(gap+d)/2.), size=(L1display,d,d), color=color.red)
bar1b = box(frame=frame1, pos=(L1display/2.-d/2.,0,(gap+d)/2.), size=(L1display,d,d), color=color.red)

axle1 = cylinder(frame=frame1, pos=((2.0-L1),0,-(gap+d)/2.), axis=(0,0,gap+d),
                 radius=axle.radius*L1*5, color=axle.color)

bar2 = box(frame=frame1, pos=(-L2display/2.+d/2.,0,-(gap+d)/2.), size=(L2display,d,d), color=color.red)
bar2b = box(frame=frame1, pos=(-L2display/2.+d/2.,0,(gap+d)/2.), size=(L2display,d,d), color=color.red)

axle2 = cylinder(frame=frame1, pos=(-L1,0,-(gap+d)/2.), axis=(0,0,gap+d),
                 radius=axle.radius*L1*5, color=axle.color)

frame1.axis = (0,-1,0)

theta1 = 1.57         # initial upper angle (from vertical)

theta1dot = 0.     # initial rate of change of theta1

frame1.rotate(axis=(0,0,1), angle=theta1)

scene.autoscale = 0

dt = 0.000025
t = 0.

I1  = (M1*(2.0-L1)**2/12)+(M1*(L1)**2/12)      # moment of inertia
C22 = (0.25*M1*(2.0-L1)**2)+(0.25*M1*L1**2)+I1

#### For energy check:
#gdisplay(x=800)
#gK = gcurve(color=color.yellow)
#gU = gcurve(color=color.cyan)
#gE = gcurve(color=color.red)

while 1:
    rate(1.0/dt)

    # Calculate accelerations of the Lagrangian coordinates:
    atheta1 = (-0.5*M1*g*(2.0-L1)*sin(theta1)+0.5*M1*g*L1*sin(theta1))/C22

    # Update velocities of the Lagrangian coordinates:
    theta1dot = theta1dot+atheta1*dt

    # Update Lagrangian coordinates:
    dtheta1 = theta1dot*dt

    theta1 = theta1+dtheta1
   
    frame1.rotate(axis=(0,0,1), angle=dtheta1)

    t = t+dt

#### Energy check:
# Yellow
    K = +(.5*(.25*M1*(2.0-L1)**2.+I1)*theta1dot**2+.5*(.25*M1*L1**2.+I1)*theta1dot**2.)
# Cyan
    U = -0.5*M1*g*(2.0-L1)*cos(theta1)+0.5*M1*g*L1*cos(theta1)
    
#    gK.plot(pos=(t,K))
#    gU.plot(pos=(t,U))
#    gE.plot(pos=(t,K+U))

[ME]

I suspect your MoI formula is wrong (don't know if that solves your issue):
MoI (a beam from pivot to radius r) = (m·r²)/3

Perhaps I could create a function (Could take a while, I'm not used to Python)

[agor95]

@ME

The improvement in the formula is fine.

I can implement

Just done a quick learn.

The MOI was using the rotation around the center of mass.

Replaced it with the rotation from one end.

The mass distribution in the pendulum is mostly in the mass 1kg.
I can have a bar component using the formula.

I am not sure what the (m·r²)/4 is for at this time.

This variance in total energy did not happen during the single pendulum tests.

[ME]

Now I know again why I don't like Python...
Sorry for the large amounts of remarks, I think they're necessary.

Code: Select all

#######################
# Physics Formulas
#######################

def pf_MoI( m, r):
    # Calculate Moment of Inertia for a single point
    # Output: real  MoI [kg m^2]
    # Input:
    # m : real      Mass[kg]
    # r : real      Radius [m]
    return m * r**2

def pf_MoI_beam( m, r):
    # Calculate Moment of Inertia for a beam 
    # Output: real  MoI [kg m^2]
    # Input:
    # m : real      Mass[kg]
    # r : real      Radius [m]
    return m * r**2 / 3

def pf_MoI_centerbeam( m, r):
    # Calculate Moment of Inertia for a beam (at its center)
    # Output: real  MoI [kg m^2]
    # Input:
    # m : real      Mass[kg]
    # r : real      Radius [m]
    return m * r**2 / 12

def pf_E_RK( I, w):
    # Calculate Energy (Rotational Kinetic)
    # Output: real  E_RK [J]=[kg (m/s)^2]
    # Input:
    # I : real      Inertia[kg m^2]
    # w : real      Angular speed [rad / s^2]
    return 0.5 * I * w**2

#----------------------------------
# check some forumula's before use
#----------------------------------
def pf_chk_MoI_beam( m, r1, r2):
    # Evaluation formula
    # Calculate Moment of Inertia for a beam going from radii r1 to r2 
    # Output: real  MoI [kg m^2]
    # Input:
    # m : real      Mass[kg]
    # r1 : real     Radius 1 [m]
    # r2 : real     Radius 2 [m]

    # Calculate this in parts
    parts = 100
    m_part= m/parts
    length= (r2-r1)
    partlength = length/parts
    MoI   = 0
    
    for part in range(0,parts):
        r_part = r1 + partlength*(part+0.5) # +0.5 because Center of gravity of that part
        MoI = MoI + m_part * r_part**2
    return MoI

#--------------------------------
def pf_check():
    m = 7
    r1 = 0
    r2 = 4
    print("Physics Formula Check:")
    print("m=", m," r1=", r1, " r2=", r2, " Len=", r2-r1 ) 
    print("MoI_beam       = ", pf_MoI_beam(m,r2) )
    print("MoI_centerbeam = ", pf_MoI_centerbeam(m,r2) )
    print("MoI_chk (",r1," .. ",r2,") = ", pf_chk_MoI_beam( m,r1,r2) )
    r1=r1-1
    r2=r2-1
    print("MoI_chk (",r1," .. ",r2,") = ", pf_chk_MoI_beam( m,r1,r2) )
    r1=r1-1
    r2=r2-1
    print("MoI_chk (",r1," .. ",r2,") = ", pf_chk_MoI_beam( m,r1,r2) )
#--------------------------------

pf_check()

With pf_check() it outputs this (in Python 3.5):

Code: Select all

Physics Formula Check:
m= 7  r1= 0  r2= 4  Len= 4
MoI_beam       =  37.333333333333336
MoI_centerbeam =  9.333333333333334
MoI_chk ( 0  ..  4 ) =  37.332400000000014
MoI_chk ( -1  ..  3 ) =  16.332400000000007
MoI_chk ( -2  ..  2 ) =  9.332400000000003

Thus (if I'm correct here) your MoI formula could be replaced by this:

Code: Select all

L2 = 2 - L1
. . .
# I1  = (M1*(2.0-L1)**2/12)+(M1*(L1)**2/12)      # moment of inertia 

I1 = pf_MoI_beam(M1,L2) + pf_MoI_beam(M1,L1)

I don't know about the rest of those formula's

[agor95]

It is more important to enjoy what you do than the unpleasant presentation process.

For example; I learned a lot from seeing the process of creating the proof for MOI than just using them blind.

If you use pseudo code making sure the SI units are declared etc.

That would allow use to convert to whatever Computer Language we choose.

That also frees you up with the boring syntax of the implementation.

Thank you for you efforts.

You are helping to create a good foundation in the Forum.

[ME]

If you use pseudo code making sure the SI units are declared etc.
That would allow use to convert to whatever Computer Language we choose.
That also frees you up with the boring syntax of the implementation.

The main reason why many people just use simulators.

I learned a lot from seeing the process of creating the proof for MOI than just using them blind

That's why it's interesting to create your own. :-)

Thus (if I'm correct here) your MoI formula could be replaced by this:

Hmm...glanced over it again: apparently not, I didn't check as I warned myself....

For that moment I assumed a left beam and a right beam with each their weight (that's what I would have done)
It actually is a single beam with one weight but shifted, just as my actual MoI-check-function did.

So we need another formula for that beam MoI.
I found it to be this one: MoI = m·(s² + L·s + L²/3)
Where L is the beam-length, and the 's' is the shift towards the center.

Code: Select all

def pf_MoI_balancebeam( m, Length, shift):
    # Calculate Moment of Inertia for a beam (shifted from the center)
    # note: shift should be Negative for a balanced shift
    #       shift= 0, acts like a normal beam
    #       shift= -Length/2, acts like a centered beam
    # Output: real    MoI [kg m^2]
    # Input:
    # m : real        Mass[kg]
    # Length : real   Length [m]
    # shift: real     shifted [m]
    # assert shift <= 0, 'The shift should be negative'
    return m * &#40;&#40;shift + Length&#41;*shift + Length**2 / 3&#41;

Then use it like:

Code: Select all

	I1 = pf_MoI_balancebeam( M1, 2, -L1)

	# that should give about the same result as:
	I1_check = pf_chk_MoI_beam( m, -0.9, 2-0.9)

That should be (more :-) correct.

This variance in total energy did not happen during the single pendulum tests.

Your last program seems to calculate the effects of an out-of-balanced beam, which acts like a pendulum.

[agor95]

To assist; I am going to remove the mass 1kg and just use a beam/bar.

We can then use a standard

MoI = m.L²/3 formula

The summation of forces as proof for the formula
was helpful.

Note. I found the thread for that documentation.
The report states that the simulator used appeared to
have error correction sub-code.