Turtle Physics

globals[damping]

patches-own[density]

turtles-own [
  mass
  vx  ; Horizontal velocity
  vy  ; Vertical velocity
  ax  ; Horizontal acceleration
  ay  ; Vertical acceleration
  newxcor
  newycor
  radius
  restitution
]

to setup
  clear-all

  set damping 1 - damping-co  ; Velocity damping per tick (simulates air resistance)
  create-turtles ball-num [

    setxy random-xcor random-ycor
    set shape "circle"
    set radius 1.0 + random-float 1.5  ; Random radius between 0.5 and 2.0
    set size 2 * radius  ; Visual size based on collision radius

    set mass radius ^ 2
    set vx random-float 4 - 2  ; Random initial velocity
    set vy random-float 4 - 2
    set ax 0
    set ay 0
    set restitution 0.6  ; Bounciness between 0.7 and 1.0

  ]
  reset-ticks
end 

to go
  click
  apply-forces      ; Apply gravity
  update-velocities ; Integrate acceleration into velocity
  apply-damping     ; Reduce velocity due to damping
  detect-collisions ; Resolve collisions between turtles
  move-turtles      ; Update positions and handle edge bouncing
  tick
end 

to apply-forces
  ask turtles [
    set ay ay - gravity  ; Gravity affects vertical acceleration
  ]
end 

to update-velocities
  ask turtles [
    set vx vx + ax
    set vy vy + ay
    set ax 0  ; Reset acceleration
    set ay 0
  ]
end 

to apply-damping
  set damping 1 - damping-co  ; Velocity damping per tick (simulates air resistance)
  ask turtles [
    set vx vx * damping
    set vy vy * damping
  ]
end 

to detect-collisions
  ;; We only handle each pair once: for turtle i, we only compare with turtles j where j > i.
  ask turtles [
    let me self
    ask other turtles with [who > [who] of me] [

      ;; Compute the vector between me and the other turtle
      let detx (xcor - [xcor] of me)
      let dety (ycor - [ycor] of me)
      let dist distance me
      let sum-r (radius + [radius] of me)

      ;; Check overlap
      if dist  < sum-r [
        ;; If dist = 0, skip or nudge it to avoid divide-by-zero
        if dist = 0 [
          set dist 0.00000001
        ]

        ;;---------------------------------
        ;; (1) Separate the overlapping turtles
        ;;---------------------------------
        let overlap (sum-r - dist)
        let nx (detx / dist)
        let ny (dety / dist)


        let halfOverlap (overlap / 2)

        ;; Move this turtle away
        set newxcor xcor + (nx * halfOverlap)
        set newycor ycor + (ny * halfOverlap)




        ;; Move the other turtle away
        ask me [
          set newxcor xcor - (nx * halfOverlap)
          set newycor ycor - (ny * halfOverlap)
        ]

        ;;---------------------------------
        ;; (2) Compute new velocities (1D elastic collision)
        ;;---------------------------------
        let m1 mass
        let m2 [mass] of me

        let vx1 vx
        let vy1 vy
        let vx2 [vx] of me
        let vy2 [vy] of me

        set vx1 lim vx1
        set vy1 lim vy1
        set vx2 lim vx2
        set vy2 lim vy2


        ;; Compute normal velocity components
        let un1 (vx1 * nx + vy1 * ny)
        let un2 (vx2 * nx + vy2 * ny)

        ;; Compute tangential velocity components (unchanged by elastic collision)
        let ut1 (vx1 * -1 * ny + vy1 * nx)
        let ut2 (vx2 * -1 * ny + vy2 * nx)

        ;; Elastic collision formulas
        let un1Prime ((un1 * (m1 - m2) + 2 * m2 * un2) / (m1 + m2))

        let un2Prime ((un2 * (m2 - m1) + 2 * m1 * un1) / (m1 + m2))

        ;; Recombine into (vx, vy)
        set vx (un1Prime * nx + ut1 * -1 * ny)
        set vy (un1Prime * ny + ut1 *  nx)

        ask me [
          set vx (un2Prime * nx + ut2 * -1 * ny)
          set vy (un2Prime * ny + ut2 *  nx)
        ]
      ]
    ]
  ]
end 

to move-turtles
  ask turtles [

    ask patch-here[
      set pcolor [color]of myself
      set density 100
    ]

    ; Calculate new position
    set newxcor xcor + vx
    set newycor ycor + vy
    if 1 = 1[
    ; Bounce off horizontal edges
    if abs newxcor > max-pxcor - 1.5 [
      set vx (-1 * vx * restitution)
      set newxcor (sign newxcor) * (max-pxcor - 1.5)
    ]

    ; Bounce off vertical edges
    if abs newycor > max-pycor - 1.5[
      set vy (-1 * vy * restitution)
      set newycor (sign newycor) * (max-pycor - 1.5)
    ]
    ]
    setxy newxcor newycor
  ]

  ask patches[
    set density density * 0.9
    set pcolor scale-color pcolor density 0 100
  ]
end 

to-report sign [number]
  ifelse number > 0
  [ report 1
  ][
    ifelse number < 0[
      report -1
    ][
      report 0
    ]
  ]
end 

to-report lim [number]
  ifelse  number > 1E30[
    report 1E30][
    ifelse number < -1E30[
      report -1E30
    ][
     report number
    ]
  ]
end 

to up

  ask turtles[
    set ay ay + 1
  ]
end 

to click
  if mouse-down?[
    let x mouse-xcor
    let y mouse-ycor
    ask turtles[
      let d distancexy x y
      set ax ax + ( x - xcor ) / d * 0.05 / mass
      set ay ay + ( y - ycor ) / d * 0.05 / mass
    ]
  ]
end