Rectangle (simple physics) - Printable Version

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

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