Rectangle (simple physics)

[1micha.elok]

       
click the image to zoom in

The code is simply an attempt to simulate a rectangle falling from the top to the ground. Please be aware that there are still some bugs and errors, which do not accurately reflect real-world physics.

Code:

'===================================
'
' Physics Simulation of a Rectangle
'
'===================================

set window "Rectangle", 800, 600, false
set redraw off

' Initialize variables
x = 400          ' Center x-coordinate of the rectangle
y = 100          ' Center y-coordinate of the rectangle
width1 = 100     ' Width of the rectangle
height1 = 50     ' Height of the rectangle

' Physics
vx = 0           ' Horizontal velocity
vy = 0           ' Vertical velocity
gravity = 0.5    ' Gravity acceleration
elasticity = 0.7 ' Bounciness factor (0 to 1)
damping = 0.98   ' Damping factor to simulate softness
friction = 0.02  ' Friction coefficient
ground = 450     ' Y-coordinate of the ground
angle = 0        ' Initial angle of rotation (in degrees)
angularVelocity = 15 ' Angular velocity (rotation speed)
momentOfInertia = (width1^2 + height1^2) / 12 ' Moment of inertia for a rectangle
torqueFactor = 0.01 ' Factor for torque applied during collision

' miscellaneous
visible newX, newY
contact_x1 = 0
contact_x2 = 0

' Function to rotate a point around the origin
function rotatePoint(px, py, angleDeg)
    cosA = cos(rad(angleDeg))
    sinA = sin(rad(angleDeg))
    newX = px * cosA - py * sinA
    newY = px * sinA + py * cosA
endfunc

' Function to find the lowest value among four numbers
function lowestValue(a, b, c, d)
    smallest = a
    if b < smallest then smallest = b
    if c < smallest then smallest = c
    if d < smallest then smallest = d
    return smallest
endfunc

' Determine the sign of the input value vx
function sign(vx)
    if vx > 0 then
        return 1
    else if vx < 0 then
        return -1
    else
        return 0
    endif
endfunc


'-----------
' Main loop
'-----------
while not keydown(KEY_ESCAPE, true)
    ' Clear the screen
    set color 0, 0, 0
    cls
   
    ' Update velocities and positions
    vy = vy + gravity       ' Apply gravity
    y = y + vy              ' Update vertical position
    x = x + vx              ' Update horizontal position
    angle = angle + angularVelocity ' Update rotation angle
   
    ' Calculate the four corners of the rectangle after rotation
    cx1 = -width1/2
    cy1 = -height1/2
    cx2 = width1/2
    cy2 = -height1/2
    cx3 = width1/2
    cy3 = height1/2
    cx4 = -width1/2
    cy4 = height1/2
   
    ' Rotate the corners around the center of the rectangle
    rotatePoint(cx1, cy1, angle); cx1 = newX; cy1 = newY
    rotatePoint(cx2, cy2, angle); cx2 = newX; cy2 = newY
    rotatePoint(cx3, cy3, angle); cx3 = newX; cy3 = newY
    rotatePoint(cx4, cy4, angle); cx4 = newX; cy4 = newY
   
    ' Translate the corners to the rectangle's position
    x1 = x + cx1; y1 = y + cy1
    x2 = x + cx2; y2 = y + cy2
    x3 = x + cx3; y3 = y + cy3
    x4 = x + cx4; y4 = y + cy4
   
    ' Check for collision with the floor
    min_y = lowestValue(y1, y2, y3, y4) ' Find the lowest point of the rectangle
    if min_y >= ground then
        ' Prevent the rectangle from going below the ground
        penetration = min_y - ground
        y = y - penetration
       
        ' Reverse and reduce vertical velocity
        vy = -vy * elasticity
       
        ' Identify the contact points (corners touching the ground)
        contact_x1 = 0; contact_x2 = 0
        if abs(y1 - ground) < 1 then contact_x1 = x1
        if abs(y2 - ground) < 1 then contact_x2 = x2
        if abs(y3 - ground) < 1 then contact_x2 = x3
        if abs(y4 - ground) < 1 then contact_x1 = x4
       
        ' Ensure contact_x1 <= contact_x2
        if contact_x1 > contact_x2 then
            temp = contact_x1
            contact_x1 = contact_x2
            contact_x2 = temp
        endif
       
        ' Calculate torque based on contact points
        torque = 0
        if contact_x1 <> 0 then torque = torque + (contact_x1 - x) * torqueFactor
        if contact_x2 <> 0 then torque = torque + (contact_x2 - x) * torqueFactor
        angularVelocity = angularVelocity + torque / momentOfInertia
       
        ' Apply friction to horizontal velocity
        if abs(vx) > 0 then
            frictionForce = sign(vx) * friction
            vx = vx - frictionForce
        endif
       
        ' Gradually reduce angular velocity due to damping
        angularVelocity = angularVelocity * damping
    endif
   
    ' Simulate softness by gradually reducing velocities
    vy = vy * damping
    vx = vx * damping
   
    ' Stability check: Determine if the rectangle should topple
    com_x = x ' Center of mass x-coordinate
    com_y = y ' Center of mass y-coordinate
   
    ' Check if the center of mass is outside the base of support
    if com_x < contact_x1 or com_x > contact_x2 then
        ' Apply additional torque to simulate toppling
        angularVelocity = angularVelocity + sign(com_x - x) * 0.1
    endif
   
    ' Draw the rotated rectangle by connecting the corners
    set color 255, 255, 255
    draw poly [x1, y1, x2, y2, x3, y3, x4, y4]
   
    ' Refresh
    fwait 60
    redraw
wend

[Marcus]

Looks promising

I wish I could help you out with these things. But I know absolutely nothing when it comes to physics programming - I just "cheat" and hope no one notices.

[1micha.elok]

Looks promising

I wish I could help you out with these things. But I know absolutely nothing when it comes to physics programming - I just "cheat" and hope no one notices.

Don't worry about physics simulation, it's just for fun ! Honestly, I even like simplicity more than real physics math.  I tried to follow the math in this program, but somewhere between the nested loops and physics, my brain blue-screened. I think I accidentally solved a physics equation, summoned a parallel universe, and now my coffee is floating 

           
click the image to zoom -in

Code:

' Physics Simulation: Stacking
set window "boxes", 400, 600
set redraw off

' Constants
constant gravity = 2         ' Gravity acceleration
constant friction = 0.08     ' Friction factor (velocity decay)
constant groundY = 500       ' Ground level (pixels)
constant tolerance = 0.1     ' Tolerance for resting detection

' Color definiion
black   = [0,0,0]
white   = [255,255,255]

' Box properties
box = []
box.num = 500       'Number of boxes
box.w   = 10        'Width of each box
box.h   = 10        'Height of each box

' Initialize boxes
for i = 0 to box.num - 1
    box[i] = []
    box[i].x = rnd(width()/box.w)*13' Random X position
    box[i].y = rnd(100)             ' Random Y position (above ground)
    box[i].vx = rnd(-0.1,0.1)       ' Random X velocity
    box[i].vy = 0                   ' Initial Y velocity
    box[i].angle = 0                ' Initial angle (no rotation)
    box[i].va = 50                  ' Angular velocity
next

' Function to rotate a point around the origin
function rotatePoint(px, py, angle)
    result = dim(2)
    result[0] = px * cos(angle) - py * sin(angle)
    result[1] = px * sin(angle) + py * cos(angle)
    return result
endfunc


'-----------
' Main loop
'-----------
while not keydown(KEY_ESCAPE,true)
    ' Clear the screen
    set color black
    cls
    set color white
   
    ' Draw the ground
    draw rect 0, groundY, width(), 20
   
    ' Update and draw each box
    for i = 0 to box.num - 1
        ' Apply gravity
        box[i].vy = box[i].vy + gravity
       
        ' Update position
        box[i].x = box[i].x + box[i].vx
        box[i].y = box[i].y + box[i].vy
       
        ' Update rotation
        box[i].angle = box[i].angle + box[i].va
       
        ' Check for collision with the ground or another box below
        resting = false
       
        ' Check for resting on the ground
        if box[i].y + box.h / 2 >= groundY then
            box[i].y = groundY - box.h/2   ' Ensure no overlap with ground
            box[i].vy = 0
            resting = true
        else
            ' Check for stacking on another box
            j = 0
            while j < box.num and not resting
                if i <> j and abs(box[i].x - box[j].x) <= box.w / 2 + tolerance and
                box[i].y + box.h / 2 >= box[j].y - box.h / 2 - tolerance and box[i].y < box[j].y then
                    ' Resting on another box
                    box[i].y = box[j].y - box.h ' Ensure no overlap with the box below
                    box[i].vy = 0
                    resting = true
                endif
                j = j + 1
            wend
        endif
       
        ' Resolve overlaps iteratively
        for k = 0 to box.num - 1
            if i <> k then
                ' Check for horizontal overlap
                if abs(box[i].x - box[k].x) <= box.w / 2 + tolerance then
                    ' Check for vertical overlap
                    if box[i].y + box.h / 2 > box[k].y - box.h / 2 and box[i].y < box[k].y then
                        box[i].y = box[k].y - box.h  ' Adjust position to avoid overlap
                        box[i].vy = 0
                    endif
                endif
            endif
        next
       
        ' Stabilize rotation when resting
        if resting then
            snapAngle = round(box[i].angle / (PI / 2)) * (PI / 2)  ' Snap to nearest multiple of 90 degrees
            box[i].angle = snapAngle
            box[i].va = box[i].va * 0.1  ' Reduce angular velocity significantly
           
            ' Apply damping to horizontal velocity
            box[i].vx = box[i].vx * 0.9
        endif
       
        ' Collision with walls
        if box[i].x - box.w / 2 <= 0 or box[i].x + box.w / 2 >= 800 then
            box[i].vx = box[i].vx * -0.5  ' Reverse X velocity with damping
            box[i].va = box[i].va * -0.5  ' Reverse angular velocity with damping
        endif
       
        ' Draw the box
        set color white
        halfW = box.w / 2
        halfH = box.h / 2
       
        ' Calculate vertices of the rotated box
        points = dim(8)
        p = rotatePoint(-halfW, -halfH, box[i].angle)
        points[0] = box[i].x + p[0]
        points[1] = box[i].y + p[1]
        p = rotatePoint(halfW, -halfH, box[i].angle)
        points[2] = box[i].x + p[0]
        points[3] = box[i].y + p[1]
        p = rotatePoint(halfW, halfH, box[i].angle)
        points[4] = box[i].x + p[0]
        points[5] = box[i].y + p[1]
        p = rotatePoint(-halfW, halfH, box[i].angle)
        points[6] = box[i].x + p[0]
        points[7] = box[i].y + p[1]
       
        draw poly [points[0], points[1], points[2], points[3],
        points[4], points[5], points[6], points[7]]
       
    next
   
    ' Refresh
    fwait 60
    redraw
wend

[johnno56]

Cool demo...