setscreen ("graphics:max;400,offscreenonly")

type Node :
    record
        x, y, vx, vy, px, py : real
    end record

const g := 0
const wind := -0.0
const N := 100
const d := 0
const k := 0.5
const decay := 0.99
const wallFriction := 0.9
const radius := maxx div (N * 2)

var F : real
var mx, my, md, mxl, myl : int
var holding : int := -1
var x, y : array 1 .. N + 2 of int

var connected : array 1 .. N, 1 .. N of boolean
var nodes : array 1 .. N of Node

for i : 1 .. N
    for j : 1 .. N
        connected (i, j) := abs (i - j) = 1
    end for
    nodes (i).x := (maxx / (N - 1)) * (i - 1)
    nodes (i).px := (maxx / (N - 1)) * (i - 1)
    nodes (i).y := maxy / 2
    nodes (i).py := maxy / 2
    nodes (i).vx := 0
    nodes (i).vy := 0
end for


function whatAngle (dx, dy : real) : real
    var ratio, angle : real
    if abs (dx) > 0.0000001 then
        ratio := dy / dx
    else
        ratio := dy / 0.000001
    end if
    angle := arctand (abs (ratio))
    if dx < 0 then
        angle := 180 - angle
    end if
    if dy < 0 then
        angle := 360 - angle
    end if
    result angle
end whatAngle

loop
    mousewhere (mx, my, md)

    for i : 1 .. N
        nodes (i).vx := nodes (i).vx * decay + wind
        nodes (i).vy := nodes (i).vy * decay + g
        for j : -1 .. 1 by 2
            if i + j <= N & i + j >= 1 & connected (i, i + j) then
                F := -k * (Math.Distance (nodes (i + j).x, nodes (i + j).y, nodes (i).x, nodes (i).y) - d)
                nodes (i).vx += cosd (whatAngle (nodes (i).x - nodes (i + j).x, nodes (i).y - nodes (i + j).y)) * F
                nodes (i).vy += sind (whatAngle (nodes (i).x - nodes (i + j).x, nodes (i).y - nodes (i + j).y)) * F
            end if
        end for
        F := -k * (Math.Distance (nodes (i).px, nodes (i).py, nodes (i).x, nodes (i).y) - d) / 50
        nodes (i).vx += cosd (whatAngle (nodes (i).x - nodes (i).px, nodes (i).y - nodes (i).py)) * F
        nodes (i).vy += sind (whatAngle (nodes (i).x - nodes (i).px, nodes (i).y - nodes (i).py)) * F
    end for

    for i : 1 .. N
        x (i) := nodes (i).x div 1
        y (i) := nodes (i).y div 1
        if i ~= holding then
            if i ~= 1 & i ~= N then
                nodes (i).x += nodes (i).vx
            end if
            nodes (i).y += nodes (i).vy
            if (nodes (i).x < radius| nodes (i).x > maxx - radius) then
                nodes (i).vx *= -1
            end if
            if (nodes (i).y < radius| nodes (i).y > maxy - radius) then
                nodes (i).vy *= -1
            end if
            if (nodes (i).x < radius| nodes (i).x > maxx - radius)| (nodes (i).y < radius| nodes (i).y > maxy - radius) then
                nodes (i).vx *= wallFriction
                nodes (i).vy *= wallFriction
            end if
            nodes (i).x := min (max (radius, nodes (i).x), maxx - radius)
            nodes (i).y := min (max (radius, nodes (i).y), maxy - radius)
        end if
    end for
    x (N + 1) := maxx
    x (N + 2) := 0
    y (N + 1) := 0
    y (N + 2) := 0
    drawfillpolygon (x, y, N + 2, blue)

    if md > 0 then
        if holding = -1 then
            for i : 1 .. N
                if Math.Distance (mx, my, nodes (i).x, nodes (i).y) <= radius & i ~= 1 & i ~= N then
                    holding := i
                    exit
                end if
            end for
        end if
    else
        if holding ~= -1 then
            nodes (holding).vx := mx - mxl
            nodes (holding).vy := my - myl
        end if
        holding := -1
    end if
    if holding ~= -1 then
        nodes (holding).x := mx
        nodes (holding).y := my
    end if
    mxl := mx
    myl := my

    Time.DelaySinceLast (10)
    View.Update
    cls
end loop