03-23-2025, 10:16 AM
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