'' ============================================================
'' Atom Smash
''
'' Copyright (C) 2006 - Jeff Marshall (coder@execulink.com)
''
'' License: 1) Program is provided as-is.  
''          2) There is no warranty or guarantee.
''          3) There are no other conditions.
''
'' ============================================================
''
''
'' This particle system is not based on any real world model but
'' it does have some theory.  The math is simplified by having
'' all "atoms" of the same size and mass.  (There is no mass in
'' any of the calculations)  Also, certain physical effects are
'' simulated using factors and conditions to change the state of
'' the system.
''
'' The solid bodies are made up of "atoms". 
''
'' * Atoms have a spacing and a size (ATOMSPACE, ATOMSIZE)
'' * Atoms are affected by gravity (XGRAVITY, YGRAVITY)
'' * Atoms can bounce off the sides and bottom of the screen
''  and when they do, lose some momentum due to the
''  collision (BOUNDDECAY)
'' * Atoms have a max speed (MAXSPEED) to limit values and
''  make the simulation work better.
''
'' "Bonds" represent the relationship between two atoms.
'' When a bond is first created a distance between two atoms
'' is set to indicate how far the two atoms should be from
'' each other.  Too far or too close and the system will move
'' the atoms closer or further apart.  Let's call this the
'' balanced distance.
''
'' * The rate at which two atoms return to the balanced
''   distance is affected by stiffness (STIFFNESS).
'' * If two atoms are moving away from each other and the
''   instantaneous velocity is high enough the bond is
''   broken. (BREAKTHRESHOLD)
'' * If the magnitude of the relative velocity between two
''   atoms is high enough (DEFORMTHRESHOLD) the bond can be
''   deformed.  The amount of deformation is controlled
''   by a factor. (FDEFORM)
''
'' ============================================================

DECLARE SUB VSet (v AS ANY, x!, y!)
DECLARE SUB VNeg (v AS ANY, a AS ANY)
DECLARE SUB VSub2V (v AS ANY, a AS ANY)
DECLARE SUB VScale (v AS ANY, a AS ANY, k!)
DECLARE SUB VAddComp (v AS ANY, a AS ANY, x!, y!)
DECLARE SUB VSubComp (v AS ANY, a AS ANY, x!, y!)
DECLARE SUB VRot2V (v AS ANY, r!)
DECLARE SUB VAddComp2V (v AS ANY, x!, y!)
DECLARE SUB VNeg2V (v AS ANY)
DECLARE SUB VScale2V (v AS ANY, k!)
DECLARE SUB VAdd2V (v AS ANY, a AS ANY)
DECLARE SUB VLimit2V (v AS ANY, d!)
DECLARE SUB VSub (v AS ANY, a AS ANY, B AS ANY)
DECLARE SUB ApplyCollision (i%)
DECLARE SUB RotateAtoms (x!, y!, r!)
DECLARE SUB ApplyTension (i%)
DECLARE SUB ApplyBounds (i%)
DECLARE SUB DrawAtoms ()
DECLARE SUB SetBond (idx%, a1%, a2%, bal!, ten!, c%)
DECLARE SUB SetAtom (idx%, x1!, y1!, dx!, dy!, c%)
DECLARE SUB CreateBlock (x1!, y1!, xsiz%, ysiz%, c%)
DECLARE SUB CreateTriBlock(x1!, y1!, siz%, c%)
DECLARE SUB Initialize ()
DECLARE FUNCTION Dist! (a1%, a2%)
DECLARE FUNCTION GetFreeAtom% ()
DECLARE FUNCTION GetFreeBond% ()
declare Function FindAtom%(x!, y!)

CONST FALSE = 0, TRUE = NOT FALSE
CONST MAXATOMS = 1000
CONST MAXBONDS = 4000
CONST ATOMSIZE = 4
CONST YGRAVITY = -.004
CONST XGRAVITY = 0
CONST BOUNDDECAY = .99
CONST FDEFORM = .4
CONST ATOMSPACE = 10
CONST DEFORMTHRESHOLD = 1
CONST BREAKTHRESHOLD = 2
CONST STIFFNESS = 0.5
CONST MAXSPEED = 5

TYPE Vector
  x AS SINGLE
  y AS SINGLE
END TYPE

DECLARE FUNCTION VDist! (a AS ANY, b AS ANY)
DECLARE FUNCTION VMag! (v AS ANY)

TYPE Atom
  ex AS INTEGER
  p AS vector
  v AS vector
  flg as INTEGER
  c as INTEGER
END TYPE

TYPE Bond
  ex AS INTEGER
  a1 AS INTEGER
  a2 AS INTEGER
  ten AS SINGLE
  bal AS SINGLE
  c as INTEGER
END TYPE

DIM SHARED a(1 TO MAXATOMS) AS Atom
DIM SHARED B(1 TO MAXBONDS) AS Bond

SCREEN 18, 8, 2
WINDOW (-160, 0)-(160, 240)
screenset 1, 0

ReStart:

Randomize Timer
Initialize
CreateBlock 100, 150, 2, 5,14
RotateAtoms 0, 0, .1

'CreateBlock 0, 150, 2, 5,13
'RotateAtoms 0, 0, .1

CreateTriBlock 0, 100, 7, 15
RotateAtoms 0, 100, .5

CreateBlock -100, 150, 3, 3,12
RotateAtoms 0, 0, .1

CreateBlock 0, 2000, 5, 5,11

CreateBlock -100, 1000, 10, 3,10

CreateBlock -100, 500, 10, 8,9
RotateAtoms 0, 0, .1

CreateBlock -150, 3000, 5, 5,11
RotateAtoms -150, 3000, .5


grabatom% = 0

Dim tmx%,tmy%,mx!,my!,mw%,mb%

DO
  omb% = mb%
  GetMouse tmx%, tmy%, mw%, mb%
  mx! = (tmx% - 320) / 2
  my! = 240 - tmy% / 2

  if mb% <> 0 and omb = 0 then
    grabatom% = FindAtom%(mx!, my!)
  elseif mb% = 0 then
    grabatom% = 0
  end if

  if grabatom% > 0 then
    if a(grabatom%).ex = TRUE then
      a(grabatom%).p.x = mx!
      a(grabatom%).p.y = my!
      a(grabatom%).v.x = 0
      a(grabatom%).v.y = 0
    else
      grabatom% = 0
    end if
  end if


  FOR i% = 1 TO MAXBONDS
    IF B(i%).ex = TRUE THEN
      ApplyTension i%
    END IF
  NEXT i%
 
  FOR i% = 1 TO MAXATOMS
    IF a(i%).ex = TRUE THEN
      VAddComp2V a(i%).v, XGRAVITY, YGRAVITY  ' Apply Gravity
     
      ApplyBounds i%
     
      ApplyCollision i%
     
      VLimit2V a(i%).v, MAXSPEED               ' Limit Speed
      VAdd2V a(i%).p, a(i%).v           ' Apply Movement
      
    END IF
  NEXT i%
  
  CLS
  DrawAtoms
  'screensync
  flip

  k$ = INKEY$
  SELECT CASE k$
  CASE "R", "r"
    goto Restart
  CASE " "
    FOR i% = 1 TO MAXATOMS
      IF A(i%).ex = TRUE THEN
        VSet a(i%).v, (rnd - 0.5) * MAXSPEED * 2, (rnd - 0.5) * MAXSPEED * 2
      END IF
    NEXT i%
    
  END SELECT


LOOP UNTIL k$ = CHR$(27)


SUB ApplyBounds (i%)
 
  IF a(i%).p.y < 0 THEN
    a(i%).p.y = -a(i%).p.y
    a(i%).v.y = -a(i%).v.y * BOUNDDECAY
    a(i%).v.x = a(i%).v.x * BOUNDDECAY
  END IF
 
  IF a(i%).p.x > 160 THEN
    a(i%).p.x = 160 - (a(i%).p.x - 160)
    a(i%).v.x = -a(i%).v.x * BOUNDDECAY
  END IF

  IF a(i%).p.x < -160 THEN
    a(i%).p.x = -160 + (-a(i%).p.x - 160)
    a(i%).v.x = -a(i%).v.x * BOUNDDECAY
  END IF

END SUB

SUB ApplyCollision (i%)
 
  DIM cij AS vector
  DIM cji AS vector
  DIM rij AS vector
  DIM rji AS vector
 
  ' Check Atom i% to see if it has hit any other

  FOR j% = 1 TO MAXATOMS  ' i% - 1
    IF j% <> i% THEN
    IF a(j%).ex = TRUE THEN
     
      VSub cij, a(j%).p, a(i%).p
      VNeg cji, cij
      d! = VMag!(cij)
      IF d! < ATOMSIZE * 2 THEN
        ' We have a collision
        ' Is I moving towards J

        k! = (cij.x * cij.x + cij.y * cij.y)
       
        ' Project i.v on cij => rij
        k1! = (a(i%).v.x * cij.x + a(i%).v.y * cij.y) / k!

        IF k1! >= 0 THEN

        VScale rij, cij, k1!
       
        ' Project j.v on cji => rji
        k2! = (a(j%).v.x * cji.x + a(j%).v.y * cji.y) / k!
        VScale rji, cji, k2!

        ' Apply Results
        VAdd2V a(i%).v, rji
        VSub2V a(i%).v, rij
       
        VAdd2V a(j%).v, rij
        VSub2V a(j%).v, rji
        END IF
      END IF
    END IF
    END IF
  NEXT j%
END SUB

SUB ApplyTension (i%)
 
  DIM v AS vector
 
  a1% = B(i%).a1
  a2% = B(i%).a2

  VSub v, a(a1%).p, a(a2%).p

  d! = VMag!(v)

  IF d! / ATOMSPACE > BREAKTHRESHOLD THEN
    B(i%).ex = FALSE
    EXIT SUB
  END IF
 
  diff! = (d! - B(i%).bal) / 2 * STIFFNESS
 
  VScale2V v, diff! / B(i%).bal
  IF ABS(d! / ATOMSPACE) > DEFORMTHRESHOLD THEN
    B(i%).bal = B(i%).bal + diff! * FDEFORM
  END IF
  
  VAdd2V a(a2%).v, v
 
  VNeg2V v
  VAdd2V a(a1%).v, v

END SUB

Function FindAtom%(x!, y!)
  function = 0
  FOR i% = 1 TO MAXATOMS
    IF a(i%).ex = TRUE THEN
      if Abs(x! - a(i%).p.x) < ATOMSPACE then
        if Abs(y! - a(i%).p.y) < ATOMSPACE then
          function = i%
          exit for
        end if
      end if
    Endif
  NEXT i%
End Function

SUB CreateTriBlock(x1!, y1!, siz%, c%)
  xspc! = ATOMSPACE
  yspc! = SQR(0.75 * ATOMSPACE * ATOMSPACE)
  bal! = ATOMSPACE
  ten! = ATOMSPACE
  REDIM aidx%(0 TO siz% - 1, 0 TO siz% - 1)
  ' Create Atoms
  for x% = 0 to siz% - 1
    for y% = 0 to siz% - x% - 1
      aidx%(x%, y%) = GetFreeAtom%
      SetAtom aidx%(x%, y%), x1! + (x% + y% / 2) * xspc!, y1! + y% * yspc!, 0, 0, c%
    next y%
  next x%
  ' Create Diagonal Bonds
  for x% = 0 to siz% - 2
    for y% = 0 to siz% - x% - 2
      bidx% = GetFreeBond%
      SetBond bidx%, aidx%(x%, y%), aidx%(x%, y% + 1), ATOMSPACE, ten!, c%
      bidx% = GetFreeBond%
      SetBond bidx%, aidx%(x%, y%), aidx%(x% + 1, y%), ATOMSPACE, ten!, c%
      bidx% = GetFreeBond%
      SetBond bidx%, aidx%(x% + 1, y%), aidx%(x%, y% + 1), ATOMSPACE, ten!, c%
    NEXT y%
  NEXT x%
END SUB
SUB CreateBlock (x1!, y1!, xsiz%, ysiz%, c%)
  xspc! = ATOMSPACE
  yspc! = ATOMSPACE
  bal! = ATOMSPACE
  ten! = ATOMSPACE
  REDIM aidx%(0 TO xsiz% - 1, 0 TO ysiz% - 1)
  ' Create Atoms
  FOR x% = 0 TO xsiz% - 1
    FOR y% = 0 TO ysiz% - 1
      aidx%(x%, y%) = GetFreeAtom%
      SetAtom aidx%(x%, y%), x1! + x% * xspc!, y1! + y% * yspc!, 0, 0, c%
    NEXT y%
  NEXT x%
  ' Create Horizontal Bonds
  FOR y% = 0 TO ysiz% - 1
    FOR x% = 0 TO xsiz% - 2
      bidx% = GetFreeBond%
      SetBond bidx%, aidx%(x%, y%), aidx%(x% + 1, y%), bal!, ten!, c%
    NEXT x%
  NEXT y%
  ' Create Vertical Bonds
  FOR x% = 0 TO xsiz% - 1
    FOR y% = 0 TO ysiz% - 2
      bidx% = GetFreeBond%
      SetBond bidx%, aidx%(x%, y%), aidx%(x%, y% + 1), bal!, ten!, c%
    NEXT y%
  NEXT x%
  ' Create Diagonal Bonds
  FOR x% = 0 TO xsiz% - 2
    FOR y% = 0 TO ysiz% - 2
      bidx% = GetFreeBond%
      SetBond bidx%, aidx%(x%, y%), aidx%(x% + 1, y% + 1), SQR(bal! * bal! * 2), ten!, c%
    NEXT y%
  NEXT x%
  ' Create Diagonal Bonds
  FOR x% = 1 TO xsiz% - 1
    FOR y% = 0 TO ysiz% - 2
      bidx% = GetFreeBond%
      SetBond bidx%, aidx%(x%, y%), aidx%(x% - 1, y% + 1), SQR(bal! * bal! * 2), ten!, c%
    NEXT y%
  NEXT x%
END SUB
SUB DrawAtoms
  ' Draw all the atoms
  FOR i% = 1 TO MAXATOMS
    IF a(i%).ex = TRUE THEN
      CIRCLE (a(i%).p.x, a(i%).p.y), ATOMSIZE, a(i%).c
    END IF
  NEXT i%
  ' Draw All the Bonds
  FOR i% = 1 TO MAXBONDS
    IF B(i%).ex = TRUE THEN
      LINE (a(B(i%).a1).p.x, a(B(i%).a1).p.y)-(a(B(i%).a2).p.x, a(B(i%).a2).p.y), B(i%).c
    END IF
  NEXT i%
END SUB
FUNCTION GetFreeAtom%
  DIM i AS INTEGER
  GetFreeAtom% = 0
  FOR i = 1 TO MAXATOMS
    IF a(i).ex = FALSE THEN
      GetFreeAtom% = i
      EXIT FOR
    END IF
  NEXT i
END FUNCTION
FUNCTION GetFreeBond%
  DIM i AS INTEGER
  GetFreeBond% = 0
  FOR i = 1 TO MAXBONDS
    IF B(i).ex = FALSE THEN
      GetFreeBond% = i
      EXIT FOR
    END IF
  NEXT i
END FUNCTION
SUB Initialize
  FOR i = 1 TO MAXATOMS
    a(i).ex = FALSE
  NEXT i
  FOR i = 1 TO MAXBONDS
    B(i).ex = FALSE
  NEXT i
END SUB
SUB RotateAtoms (x!, y!, r!)
  DIM v AS vector
  FOR i% = 1 to MAXATOMS
    IF a(i%).ex = TRUE THEN
    IF a(i%).flg = FALSE THEN
      VSubComp v, a(i%).p, x!, y!
      VRot2V v, r!
      VAddComp a(i%).p, v, x!, y!
      a(i%).flg = TRUE
    END IF
    END IF
  NEXT i%
END SUB
SUB SetAtom (idx%, x1!, y1!, dx!, dy!, c%)
  a(idx%).ex = TRUE
  a(idx%).p.x = x1!
  a(idx%).p.y = y1!
  a(idx%).v.x = dx!
  a(idx%).v.y = dy!
  a(idx%).flg = FALSE
  a(idx%).c = c%
END SUB
SUB SetBond (idx%, a1%, a2%, bal!, ten!, c%)
  B(idx%).ex = TRUE
  B(idx%).a1 = a1%
  B(idx%).a2 = a2%
  B(idx%).bal = bal!
  B(idx%).ten = ten!
  B(idx%).c = c%
END SUB
SUB VAdd (v AS Vector, a AS Vector, b AS Vector)
  v.x = a.x + b.x
  v.y = a.y + b.y
END SUB
SUB VAdd2V (v AS Vector, a AS Vector)
  v.x = v.x + a.x
  v.y = v.y + a.y
END SUB
SUB VAddComp (v AS Vector, a AS Vector, x!, y!)
  v.x = a.x + x!
  v.y = a.y + y!
END SUB
SUB VAddComp2V (v AS Vector, x!, y!)
  v.x = v.x + x!
  v.y = v.y + y!
END SUB
SUB VAddMul (v AS Vector, k1!, a AS Vector, k2!, b AS Vector)
  v.x = k1! * a.x + k2! * b.x
  v.y = k1! * a.y + k2! * b.y
END SUB
FUNCTION VDist! (a AS Vector, b AS Vector)
  DIM v AS Vector
  VSub v, a, b
  VDist! = VMag!(v)
END FUNCTION
SUB VLimit2V (v AS Vector, d!)
  IF v.x > d! THEN v.x = d!
  IF v.x < -d! THEN v.x = -d!
  IF v.y > d! THEN v.y = d!
  IF v.y < -d! THEN v.y = -d!
END SUB
FUNCTION VMag! (v AS Vector)
  VMag! = SQR(v.x * v.x + v.y * v.y)
END FUNCTION
SUB VNeg (v AS Vector, a AS Vector)
  VScale v, a, -1
END SUB
SUB VNeg2V (v AS Vector)
  VScale2V v, -1
END SUB
SUB VRot2V (v AS Vector, r!)
  DIM a AS Vector
  a.x = v.x * COS(r!) + v.y * SIN(r!)
  a.y = -v.x * SIN(r!) + v.y * COS(r!)
  v = a
END SUB
SUB VRot902V (v AS Vector)
  DIM a AS Vector
  a = v
  v.x = -a.y
  v.y = a.x
END SUB
SUB VScale (v AS Vector, a AS Vector, k!)
  v.x = k! * a.x
  v.y = k! * a.y
END SUB
SUB VScale2V (v AS Vector, k!)
  v.x = k! * v.x
  v.y = k! * v.y
END SUB
SUB VSet (v AS Vector, x!, y!)
  v.x = x!
  v.y = y!
END SUB
SUB VSub (v AS Vector, a AS Vector, b AS Vector)
  v.x = a.x - b.x
  v.y = a.y - b.y
END SUB
SUB VSub2V (v AS Vector, a AS Vector)
  v.x = v.x - a.x
  v.y = v.y - a.y
END SUB
SUB VSubComp (v AS Vector, a AS Vector, x!, y!)
  v.x = a.x - x!
  v.y = a.y - y!
END SUB