Molecules

Copyright 2002, Dan Schroeder
Physics Department, Weber State University, Ogden, UT 84408-2508.  
http://physics.weber.edu/schroeder/

This program is free, open-source software.  You can redistribute it 
and/or modify it under the terms of the GNU General Public License 
as published by the Free Software Foundation, either version 2 of 
the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, 
but without any warranty; without even the implied warranty of 
merchantability or fitness for a particular purpose.  See the GNU 
General Public License for more details.

You should have received a copy of the GNU General Public License 
along with this program.  If not, contact the Free Software Foundation, 
Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA, 
www.gnu.org.

// Release notes and modification history

Version 1.0, 22 November 2002.  First public release.

There are still a couple of known bugs:  
  Shift key doesn't seem to work with little arrow controls on some Macs.
  Under Windows, full-screen save files don't open in full-screen mode.

Animation still isn't very smooth under OSX, but the problem is less
noticeable on newer, faster machines.

The Windows user interface could be improved, e.g., "ctrl-H" disappears
from menu after it's used once; menu bar appears in full-screen mode
only if main window has focus; "Quit" should really be called "Exit".

Code to compute heat capacity is disabled since it doesn't seem to be 
useful, though EkSquaredTotal is still written to save files in case it's 
needed in a later version.

Many more features would be nice to have:  recording/replaying movies;
plotting graphs of various quantities; manually setting velocities; 
permanent bonds to form diatomic molecules.

This is the source code for MacOS.  To compile for
Windows, you should check the GrowIcon box in the GWin 
properties window, and delete the ResizeWindow.  (The latter
is a kludge that's used on the Mac to prevent the grow icon 
from intruding into the window content while the simulation 
is running.)

//***********************************************************************

// App Properties:

AAattraction as single
AAForceCutoff as single
AAForceCutoff2 as single
Aax(999) as single
Aay(999) as single
ABattraction as single
ABForceCutoff as single
ABForceCutoff2 as single
ABrsum2 as single
Acolor as color
Acount as integer
Acountm1 as integer
Adiameter as single
Adiameter2 as single
Amass as single
Amassinv as single
Apixelwidth as single
Aradius as single
Ascreenx(999) as integer
Ascreeny(999) as integer
Avx(999) as single
Avy(999) as single
Ax(999) as single
Ay(999) as single
Bax(999) as single
Bay(999) as single
BBattraction as single
BBForceCutoff as single
BBForceCutoff2 as single
Bcolor as color
Bcount as integer
Bcountm1 as integer
Bdiameter as single
Bdiameter2 as single
BGcolor as color
Bmass as single
Bmassinv as single
Bpixelwidth as integer
Bradius as single
Bscreenx(999) as integer
Bscreeny(999) as integer
Bvx(999) as single
Bvy(999) as single
dt as single
dtover2 as single
dtsquaredover2 as single
Egravitational as single
EgTotal as single
Einteraction as single
EiTotal as single
Ekinetic as single
EkSquared as single
EkSquaredTotal as single
EkTotal as single
gravity as single
GWinPort as integer
GWinRegion as integer
iterations as integer
pressure as single
Ptotal as single
t as single
WallStiffness as single
WallStiffnessHalf as single
xmax as single
ymax as single

// Methods:

App.ChangeSpeeds:
Sub ChangeSpeeds(ratio as single)
  dim i as integer
  #pragma disableBackgroundTasks
  for i = 0 to Acountm1
    Avx(i) = ratio * Avx(i)
    Avy(i) = ratio * Avy(i)
  next
  for i = 0 to Bcountm1
    Bvx(i) = ratio * Bvx(i)
    Bvy(i) = ratio * Bvy(i)
  next
  ResetData  // if simulation isn't running, perhaps we should calculate new energy
End Sub

App.ColorsChanged:
Sub ColorsChanged()
  CWin.AcolorCircle.FillColor = Acolor
  CWin.BcolorCircle.FillColor = Bcolor
  CWin.BGcolorSquare.FillColor = BGcolor
  redraw
  ColorsRedBlueBeige.checked = false
  ColorsPurpleGreenSky.checked = false
  ColorsRedGreenBlue.checked = false
  ColorsSpring.checked = false
  ColorsSummer.checked = false
  ColorsFall.checked = false
  ColorsWinter.checked = false
  ColorsChristmas.checked = false
  ColorsIndependence.checked = false
  ColorsHalloween.checked = false
  ColorsFluorescent.checked = false
  ColorsEarthtones.checked = false
  ColorsBubbles.checked = false  
  ColorsGreens.checked = false
  ColorsBlues.checked = false
  ColorsPurples.checked = false
  ColorsReds.checked = false
  ColorsBrowns.checked = false
  ColorsGrays.checked = false
End Sub

App.ComputeAccelerations:
Sub ComputeAccelerations()
  // given molecules' positions, compute their accelerations; called by SingleStep
  // this is where the program spends most of its time (aside from drawing),
  // so speed is critical here
  dim i,j as integer
  dim thisx, thisy, deltax, deltay, dx2, dy2, temp, wallforce as single
  dim rsquared, r0overr2, attract, repel, foverr, fx, fy as single
  #pragma disableBackgroundTasks
  Einteraction = 0.0  // we'll also compute the potential energies
  Egravitational = 0.0
  pressure = 0.0  // and we'll compute the pressure on the walls
  // First handle bounces off walls (linear force):
  for i = 0 to Acountm1
    thisx = Ax(i)
    thisy = Ay(i)
    if gravity <> 0.0 then
      Egravitational = Egravitational + Amass * gravity * (ymax - thisy)
    end if
    if thisx < Aradius then
      temp = Aradius - thisx
      wallforce = WallStiffness * temp
      Aax(i) = wallforce
      pressure = pressure + wallforce
      Einteraction = Einteraction + WallStiffnessHalf * temp * temp
    else
      if thisx > xmax - Aradius then
        temp = xmax - Aradius - thisx
        wallforce = WallStiffness * temp
        Aax(i) = wallforce
        pressure = pressure - wallforce
        Einteraction = Einteraction + WallStiffnessHalf * temp * temp
      else
        Aax(i) = 0.0
      end if
    end if
    if thisy < Aradius then
      temp = Aradius - thisy
      Aay(i) = WallStiffness * temp
      Einteraction = Einteraction + WallStiffnessHalf * temp * temp
    else
      if thisy > ymax - Aradius then
        temp = ymax - Aradius - thisy
        Aay(i) = WallStiffness * temp
        Einteraction = Einteraction + WallStiffnessHalf * temp * temp
      else
        Aay(i) = 0.0
      end if
    end if
  next
  for i = 0 to Bcountm1
    thisx = Bx(i)
    thisy = By(i)
    if gravity <> 0.0 then
      Egravitational = Egravitational + Bmass * gravity * (ymax - thisy)
    end if
    if thisx < Bradius then
      temp = Bradius - thisx
      wallforce = WallStiffness * temp
      Bax(i) = wallforce
      pressure = pressure + wallforce
      Einteraction = Einteraction + WallStiffnessHalf * temp * temp
    else
      if thisx > xmax - Bradius then
        temp = xmax - Bradius - thisx
        wallforce = WallStiffness * temp
        Bax(i) = wallforce
        pressure = pressure - wallforce
        Einteraction = Einteraction + WallStiffnessHalf * temp * temp
      else
        Bax(i) = 0.0
      end if
    end if
    if thisy < Bradius then
      temp = Bradius - thisy
      Bay(i) = WallStiffness * temp
      Einteraction = Einteraction + WallStiffnessHalf * temp * temp
    else
      if thisy > ymax - Bradius then
        temp = ymax - Bradius - thisy
        Bay(i) = WallStiffness * temp
        Einteraction = Einteraction + WallStiffnessHalf * temp * temp
      else
        Bay(i) = 0.0
      end if
    end if
  next
  pressure = pressure / (2.0 * ymax)  // we're done computing the pressure
  // Now handle forces between molecules (Lennard-Jones interaction):
  for i = 0 to Acountm1
    thisx = Ax(i)
    thisy = Ay(i)
    for j = i+1 to Acountm1     // first do A-A interactions
      deltax = Ax(j) - thisx
      dx2 = deltax * deltax
      if dx2 < AAForceCutoff2 then
        deltay = Ay(j) - thisy
        dy2 = deltay * deltay
        if dy2 < AAForceCutoff2 then
          rsquared = dx2 + dy2
          if rsquared < AAForceCutoff2 then
            r0overr2 = Adiameter2 / rsquared
            attract = r0overr2 * r0overr2 * r0overr2
            repel = attract * attract
            foverr = 24.0 * (2.0 * repel - attract) / rsquared
            if foverr < 0.0 then
              foverr = foverr * AAattraction
              Einteraction = Einteraction + AAattraction * 4.0 * (repel - attract)
            else
              Einteraction = Einteraction + 4.0 * (repel - attract) + (1.0 - AAattraction)
            end if
            fx = foverr * deltax
            fy = foverr * deltay
            Aax(j) = Aax(j) + fx
            Aax(i) = Aax(i) - fx
            Aay(j) = Aay(j) + fy
            Aay(i) = Aay(i) - fy
          end if
        end if
      end if
    next
    for j = 0 to Bcountm1     // next do A-B interactions
      deltax = Bx(j) - thisx
      dx2 = deltax * deltax
      if dx2 < ABForceCutoff2 then
        deltay = By(j) - thisy
        dy2 = deltay * deltay
        if dy2 < ABForceCutoff2 then
          rsquared = dx2 + dy2
          if rsquared < ABForceCutoff2 then
            r0overr2 = ABrsum2 / rsquared
            attract = r0overr2 * r0overr2 * r0overr2
            repel = attract * attract
            foverr = 24.0 * (2.0 * repel - attract) / rsquared
            if foverr < 0.0 then
              foverr = foverr * ABattraction
              Einteraction = Einteraction + ABattraction * 4.0 * (repel - attract)
            else
              Einteraction = Einteraction + 4.0 * (repel - attract) + (1.0 - ABattraction)
            end if
            fx = foverr * deltax
            fy = foverr * deltay
            Bax(j) = Bax(j) + fx
            Aax(i) = Aax(i) - fx
            Bay(j) = Bay(j) + fy
            Aay(i) = Aay(i) - fy
          end if
        end if
      end if
    next
    Aax(i) = Aax(i) * Amassinv  // A's are done so divide by mass to get acceleration
    Aay(i) = Aay(i) * Amassinv + gravity
  next
  for i = 0 to Bcountm1   // now do B-B interactions
    thisx = Bx(i)
    thisy = By(i)
    for j = i+1 to Bcountm1
      deltax = Bx(j) - thisx
      dx2 = deltax * deltax
      if dx2 < BBForceCutoff2 then
        deltay = By(j) - thisy
        dy2 = deltay * deltay
        if dy2 < BBForceCutoff2 then
          rsquared = dx2 + dy2
          if rsquared < BBForceCutoff2 then
            r0overr2 = Bdiameter2 / rsquared
            attract = r0overr2 * r0overr2 * r0overr2
            repel = attract * attract
            foverr = 24.0 * (2.0 * repel - attract) / rsquared
            if foverr < 0.0 then
              foverr = foverr * BBattraction
              Einteraction = Einteraction + BBattraction * 4.0 * (repel - attract)
            else
              Einteraction = Einteraction + 4.0 * (repel - attract) + (1.0 - BBattraction)
            end if
            fx = foverr * deltax
            fy = foverr * deltay
            Bax(j) = Bax(j) + fx
            Bax(i) = Bax(i) - fx
            Bay(j) = Bay(j) + fy
            Bay(i) = Bay(i) - fy
          end if
        end if
      end if
    next
    Bax(i) = Bax(i) * Bmassinv  //  and now the B's are done
    Bay(i) = Bay(i) * Bmassinv + gravity
  next
End Sub

App.createA:
Sub createA()
  dim x, y as single
  dim trial as integer
  dim successful as boolean
  if Acount < 1000 then
    trial = 0 
    do
      x = (rnd * (xmax - Adiameter)) + Aradius
      y = (rnd * (ymax - Adiameter)) + Aradius
      trial = trial + 1
      successful = SpaceAvailable(x,y,Aradius)
    loop until successful or (trial > 99)
    if successful then
      Ax(Acount) = x
      Ay(Acount) = y
      Avx(Acount) = rnd - 0.5
      Avy(Acount) = rnd - 0.5
      Aax(Acount) = 0.0
      Aay(Acount) = 0.0
      Ascreenx(Acount) = -100
      drawA(Acount)
      Acount = Acount + 1
      Acountm1 = Acount - 1
      ResetData
    end if
    CWin.ANEdit.text = Format(Acount,"#")
  end if
End Sub

App.createB:
Sub createB()
  dim x, y as single
  dim trial as integer
  dim successful as boolean
  if Bcount < 1000 then
    trial = 0 
    do
      x = (rnd * (xmax - Bdiameter)) + Bradius
      y = (rnd * (ymax - Bdiameter)) + Bradius
      trial = trial + 1
      successful = SpaceAvailable(x,y,Bradius)
    loop until successful or (trial > 99)
    if successful then
      Bx(Bcount) = x
      By(Bcount) = y
      Bvx(Bcount) = rnd - 0.5
      Bvy(Bcount) = rnd - 0.5
      Bax(Bcount) = 0.0
      Bay(Bcount) = 0.0
      Bscreenx(Bcount) = -100
      drawB(Bcount)
      Bcount = Bcount + 1
      Bcountm1 = Bcount - 1
      ResetData
    end if
    CWin.BNEdit.text = Format(Bcount,"#")
  end if
End Sub

App.deleteA:
Sub deleteA(index as integer)
  if index >= 0 then
    Ax(index) = -100
    drawA(index)  // that should erase it
    if index < Acount-1 then   // if it's not the last one in the list...
      Ax(index) = Ax(Acount-1)  // move the last one into this slot
      Ay(index) = Ay(Acount-1)
      Avx(index) = Avx(Acount-1)
      Avy(index) = Avy(Acount-1)
      Aax(index) = Aax(Acount-1)
      Aay(index) = Aay(Acount-1)
      Ascreenx(index) = Ascreenx(Acount-1)
      Ascreeny(index) = Ascreeny(Acount-1)
    end if
    Acount = Acount - 1
    Acountm1 = Acount - 1
    CWin.ANEdit.text = Format(Acount,"#")
    ResetData
  end if
End Sub

App.deleteB:
Sub deleteB(index as integer)
  if index >= 0 then
    Bx(index) = -100
    drawB(index)  // that should erase it
    if index < Bcount-1 then   // if it's not the last one in the list...
      Bx(index) = Bx(Bcount-1)  // move the last one into this slot
      By(index) = By(Bcount-1)
      Bvx(index) = Bvx(Bcount-1)
      Bvy(index) = Bvy(Bcount-1)
      Bax(index) = Bax(Bcount-1)
      Bay(index) = Bay(Bcount-1)
      Bscreenx(index) = Bscreenx(Bcount-1)
      Bscreeny(index) = Bscreeny(Bcount-1)
    end if
    Bcount = Bcount - 1
    Bcountm1 = Bcount - 1
    CWin.BNEdit.text = Format(Bcount,"#")
    ResetData
  end if
End Sub

App.DisplayData:
Sub DisplayData()
  dim EkAv, EkSquaredAv as single
  CWin.Etotal.text = Format(Ekinetic+Einteraction+Egravitational,"-#.00")
  CWin.DataBox.list(0) = Format(t,"#.00")
  if CWin.AvButton.value and (iterations > 0) then
    CWin.DataBox.list(1) = Format(EkTotal/iterations,"#.00")
    CWin.DataBox.list(2) = Format(EiTotal/iterations,"-#.00")
    CWin.DataBox.list(3) = Format(EgTotal/iterations,"#.00")
    if Acount + Bcount <> 0 then
      EkAv = EkTotal / iterations
      //EkSquaredAv = EkSquaredTotal / iterations
      CWin.DataBox.list(4) = Format(EkAv/(Acount + Bcount),"#.0000")
      //CWin.DataBox.list(6) = Format( 1 / (1+ (1/(Acount+Bcount)) - (EkSquaredAv/(EkAv*EkAv)))  ,"#.00")
        // heat capacity formula from Ray and Graben, Mol Phys 43, 1293 (1981)
    end if
    CWin.DataBox.list(5) = Format(Ptotal/iterations, "#.0000")
  else
    CWin.DataBox.list(1) = Format(Ekinetic,"#.00")
    CWin.DataBox.list(2) = Format(Einteraction,"-#.00")
    CWin.DataBox.list(3) = Format(Egravitational,"#.00")
    if Acount + Bcount <> 0 then
      CWin.DataBox.list(4) = Format(Ekinetic / (Acount + Bcount),"#.0000")
    end if
    CWin.DataBox.list(5) = Format(pressure, "#.0000")
    //CWin.DataBox.list(6) = "?.??"  // can't compute instantaneous heat capacity
  end if
End Sub

App.drawA:
Sub drawA(index as integer)
  dim x,y,oldx,oldy as integer
  oldx = Ascreenx(index)
  oldy = Ascreeny(index)
  x = (Ax(index) - Aradius) * 10.0 + 0.5
  y = (Ay(index) - Aradius) * 10.0 + 0.5
  if (x <> oldx) or (y <> oldy) then
    GWin.graphics.forecolor = BGcolor
    GWin.graphics.filloval oldx,oldy,Apixelwidth,Apixelwidth
    GWin.graphics.forecolor = Acolor
    GWin.graphics.filloval x,y,Apixelwidth,Apixelwidth
    Ascreenx(index) = x
    Ascreeny(index) = y
  end if
End Sub

App.drawB:
Sub drawB(index as integer)
  dim x,y,oldx,oldy as integer
  oldx = Bscreenx(index)
  oldy = Bscreeny(index)
  x = (Bx(index) - Bradius) * 10.0 + 0.5
  y = (By(index) - Bradius) * 10.0 + 0.5
  if (x <> oldx) or (y <> oldy) then
    GWin.graphics.forecolor = BGcolor
    GWin.graphics.filloval oldx,oldy,Bpixelwidth,Bpixelwidth
    GWin.graphics.forecolor = Bcolor
    GWin.graphics.filloval x,y,Bpixelwidth,Bpixelwidth
    Bscreenx(index) = x
    Bscreeny(index) = y
  end if
End Sub

App.OpenTheFile:
Sub OpenTheFile(f as FolderItem)
  Dim stream as BinaryStream
  Dim firstbyte, secondbyte, version, subversion as Integer
  Dim newwidth, newheight as integer
  Dim r,g,b as Integer
  Dim i as Integer
  #pragma disableBackgroundTasks  //  so the for loops won't take forever
  If f <> Nil then
    if CWin.simTimer.mode = 2 then
      CWin.startButton.push  // if it's running, stop (prevents errors under Windows)
    end if
    stream = f.OpenAsBinaryFile(false)
    firstbyte = stream.ReadByte  // first two bytes must be ASCII "MD" for Molecular Dynamics
    secondbyte = stream.ReadByte
    version = stream.ReadByte
    subversion = stream.ReadByte
    If (firstbyte <> 77) or (secondbyte <> 68) or (version <> 1) or (subversion <> 0) then
      MsgBox "Sorry, can't read that file."
    Else
      Acount = 0
      Acountm1 = -1
      Bcount = 0
      Bcountm1 = -1
      xmax = stream.ReadSingle
      ymax = stream.ReadSingle
      newwidth = round(10 * xmax)
      newheight = round(10 * ymax)
      // the following isn't fool-proof, but should work fine most of the time on Macs;
      // for some reason, the attempt to enter full-screen mode fails on Windows
      if (newwidth >= screen(0).width) and (newheight >= screen(0).height) then
        WindowFullScreen.checked = true
        GWin.left = 0
        GWin.top = 0
        GWin.width = screen(0).width
        GWin.height = screen(0).height
        GWin.MenuBarVisible = false
        GWin.FullScreen = true
        CWin.left = screen(0).width - 208
      else
        WindowFullScreen.checked = false
        GWin.FullScreen = false
        GWin.MenuBarVisible = true
        if (GWin.top + newheight > screen(0).height) or (GWin.top < 50) then
          GWin.top = 50
        end if
        if GWin.left + newwidth > screen(0).width then
          GWin.left = 10
        end if
        GWin.width = 10 * xmax
        GWin.height = 10 * ymax
      end if
      setGravity(stream.ReadSingle)
      setdt(stream.ReadSingle)
      Acount = stream.ReadShort
      Acountm1 = Acount - 1
      CWin.ANEdit.text = Format(Acount,"#")
      setAdiameter(stream.ReadSingle)
      setAmass(stream.ReadSingle)
      setAAattraction(stream.ReadSingle)
      r = stream.ReadByte
      g = stream.ReadByte
      b = stream.ReadByte
      Acolor = RGB(r,g,b)
      Bcount = stream.ReadShort
      Bcountm1 = Bcount - 1
      CWin.BNEdit.text = Format(Bcount,"#")
      setBdiameter(stream.ReadSingle)
      setBmass(stream.ReadSingle)
      setBBattraction(stream.ReadSingle)
      r = stream.ReadByte
      g = stream.ReadByte
      b = stream.ReadByte
      Bcolor = RGB(r,g,b)
      setABattraction(stream.ReadSingle)
      setWallStiffness(stream.ReadSingle)
      r = stream.ReadByte
      g = stream.ReadByte
      b = stream.ReadByte
      BGcolor = RGB(r,g,b)
      ColorsChanged
      iterations = stream.ReadLong
      t = stream.ReadSingle
      Ekinetic = stream.ReadSingle
      Einteraction = stream.ReadSingle
      Egravitational = stream.ReadSingle
      EkTotal = stream.ReadSingle
      EkSquaredTotal = stream.ReadSingle
      EiTotal = stream.ReadSingle
      EgTotal = stream.ReadSingle
      pressure = stream.ReadSingle
      pTotal = stream.ReadSingle
      for i = 0 to Acountm1
        Ax(i) = stream.ReadSingle
        Ay(i) = stream.ReadSingle
        Avx(i) = stream.ReadSingle
        Avy(i) = stream.ReadSingle
        Aax(i) = stream.ReadSingle
        Aay(i) = stream.ReadSingle
        Ascreenx(i) = -100
      next
      for i = 0 to Bcountm1
        Bx(i) = stream.ReadSingle
        By(i) = stream.ReadSingle
        Bvx(i) = stream.ReadSingle
        Bvy(i) = stream.ReadSingle
        Bax(i) = stream.ReadSingle
        Bay(i) = stream.ReadSingle
        Bscreenx(i) = -100
      next  // should really check to make sure none are off the screen
      redraw
      DisplayData
    end if
    stream.close
  end if
End Sub

App.redraw:
Sub redraw()
  dim i as integer
  dim x, y as integer
  #pragma disableBackgroundTasks
  GWin.graphics.forecolor = BGcolor
  GWin.graphics.fillrect 0,0,GWin.width, GWin.height
  GWin.graphics.forecolor = Acolor
  for i = 0 to Acountm1
    x = (Ax(i) - Aradius) * 10.0 + 0.5
    y = (Ay(i) - Aradius) * 10.0 + 0.5
    GWin.graphics.filloval x,y,Apixelwidth,Apixelwidth
    Ascreenx(i) = x  // these 2 lines may not be necessary
    Ascreeny(i) = y
  next
  GWin.graphics.forecolor = Bcolor
  for i = 0 to Bcountm1
    x = (Bx(i) - Bradius) * 10.0 + 0.5
    y = (By(i) - Bradius) * 10.0 + 0.5
    GWin.graphics.filloval x,y,Bpixelwidth,Bpixelwidth
    Bscreenx(i) = x  // these 2 lines may not be necessary
    Bscreeny(i) = y
  next
End Sub

App.ResetData:
Sub ResetData()
  iterations = 0
  t = 0.0
  EkTotal = 0.0
  EkSquaredTotal = 0.0
  EiTotal = 0.0
  EgTotal = 0.0
  Ptotal = 0.0
  DisplayData
End Sub

App.RunSimulation:
Sub RunSimulation()
  // this is the main simulation loop; it's called by the timer in the control window
  #if TargetCarbon
    declare sub QDSetDirtyRegion lib "CarbonLib" (port as integer, region as integer)
    declare sub LockPortBits lib "CarbonLib" (port as integer)
    declare sub UnlockPortBits lib "CarbonLib" (port as integer)
  #endif
  dim starttime, flushtime as double
  dim time1, time2 as double // for diagnosis only--delete later
  starttime = microseconds
  flushtime = starttime
  do
    #if TargetCarbon
      LockPortBits GWinPort  // tricks to speed up graphics on OS X
      //QDSetDirtyRegion GWinPort, GWinRegion  // this seems to slow it down under most conditions
    #endif
    do
      SingleStep
    loop until (microseconds - flushtime) > 30000  // update screen about 30 times per second (OSX only)
    #if TargetCarbon
      UnlockPortBits GWinPort 
    #endif
    GWin.UpdateNow
    flushtime = microseconds
    if Ekinetic > 0.01 * (Acount + Bcount) / dtsquaredover2 then // since each v should be < .1/dt
      CWin.startButton.push   // stop running
      MsgBox "I think the calculation has become unstable."
      exit
    end if
  loop until (microseconds - starttime) > 200000  // run a total of 0.2 seconds
  DisplayData
End Sub

App.setAAattraction:
Sub setAAattraction(newAA as single)
  if newAA < 0.0 then
    AAattraction = 0.0
  else
    if newAA > 2.0 then
      AAattraction = 2.0
    else
      AAattraction = newAA
    end if
  end if
  CWin.AAEdit.text = Format(AAattraction,"#.0")
  setForceCutoffs
  ResetData
End Sub

App.setABattraction:
Sub setABattraction(newAB as single)
  if newAB < 0.0 then
    ABattraction = 0.0
  else
    if newAB > 2.0 then
      ABattraction = 2.0
    else
      ABattraction = newAB
    end if
  end if
  CWin.ABEdit.text = Format(ABattraction,"#.0")
  setForceCutoffs
  ResetData
End Sub

App.setAdiameter:
Sub setAdiameter(newd as single)
  if newd < 1.0 then
    Adiameter = 1.0
  else
    if newd > 5.0 then
      Adiameter = 5.0
    else
      Adiameter = newd
    end if
  end if
  Apixelwidth = 10 * Adiameter
  Adiameter2 = Adiameter * Adiameter
  Aradius = Adiameter / 2.0
  Adiameter2 = Adiameter * Adiameter
  ABrsum2 = (Aradius + Bradius) * (Aradius + Bradius)
  setForceCutoffs
  CWin.ADEdit.text = Format(Adiameter,"#.0")
  redraw
  ResetData
End Sub

App.setAmass:
Sub setAmass(newM as single)
  if newM < 1.0 then
    Amass = 1.0
  else
    if newM > 10.0 then
      Amass = 10.0
    else
      Amass = newM
    end if
  end if
  Amassinv = 1 / Amass
  CWin.AMEdit.text = Format(Amass,"#.0")
  ResetData
End Sub

App.setBBattraction:
Sub setBBattraction(newBB as single)
  if newBB < 0.0 then
    BBattraction = 0.0
  else
    if newBB > 2.0 then
      BBattraction = 2.0
    else
      BBattraction = newBB
    end if
  end if
  CWin.BBEdit.text = Format(BBattraction,"#.0")
  setForceCutoffs
  ResetData
End Sub

App.setBdiameter:
Sub setBdiameter(newd as single)
  if newd < 1.0 then
    Bdiameter = 1.0
  else
    if newd > 5.0 then
      Bdiameter = 5.0
    else
      Bdiameter = newd
    end if
  end if
  Bpixelwidth = 10 * Bdiameter
  Bdiameter2 = Bdiameter * Bdiameter
  Bradius = Bdiameter / 2.0
  Bdiameter2 = Bdiameter * Bdiameter
  ABrsum2 = (Aradius + Bradius) * (Aradius + Bradius)
  setForceCutoffs
  CWin.BDEdit.text = Format(Bdiameter,"#.0")
  redraw
  ResetData
End Sub

App.setBmass:
Sub setBmass(newM as single)
  if newM < 1.0 then
    Bmass = 1.0
  else
    if newM > 10.0 then
      Bmass = 10.0
    else
      Bmass = newM
    end if
  end if
  Bmassinv = 1 / Bmass
  CWin.BMEdit.text = Format(Bmass,"#.0")
  ResetData
End Sub

App.setdt:
Sub setdt(newdt as single)
  if newdt > 0.1 then
    dt = 0.1
  else
    if newdt < 0.0001 then
      dt = 0.0001
    else
      dt = newdt
    end if
  end if
  dtover2 = dt / 2.0
  dtsquaredover2 = dt * dt / 2.0
  CWin.dtEdit.text = Format(dt,"#.0000")
End Sub

App.setForceCutoffs:
Sub setForceCutoffs()
  // based on sizes and attraction strengths, sets cutoffs beyond
  // which we won't even compute the force
  if AAattraction < 0.001 then
    AAForceCutoff = Adiameter * pow(2,0.333333)
  else
    AAForceCutoff = Adiameter * 3.0
  end if
  AAForceCutoff2 = AAForceCutoff * AAForceCutoff
  if BBattraction < 0.001 then
    BBForceCutoff = Bdiameter * pow(2,0.333333)
  else
    BBForceCutoff = Bdiameter * 3.0
  end if
  BBForceCutoff2 = BBForceCutoff * BBForceCutoff
  if ABattraction < 0.001 then
    ABForceCutoff = (Aradius + Bradius) * pow(2,0.333333)
  else
    ABForceCutoff = (Aradius + Bradius) * 3.0
  end if
  ABForceCutoff2 = ABForceCutoff * ABForceCutoff
End Sub

App.setGravity:
Sub setGravity(newg as single)
  if newg > 1.00 then
    gravity = 1.00
  else
    if newg < 0.0 then
      gravity = 0.0
    else
      gravity = newg
    end if
  end if
  CWin.gEdit.text = Format(gravity,"#.000")
  ResetData
End Sub

App.SetWallStiffness:
Sub SetWallStiffness(newStiffness as single)
  if newStiffness < 1.0 then
    WallStiffness = 1.0
  else
    if newStiffness > 500.0 then
      WallStiffness = 500.0
    else
      WallStiffness = newStiffness
    end if
  end if
  WallStiffnessHalf = WallStiffness / 2.0
  CWin.WSEdit.text = Format(WallStiffness,"#")
End Sub

App.SingleStep:
Sub SingleStep()
  // carries out a single time step of the simulation
  dim i as integer
  dim tempvx, tempvy as single
  #pragma disableBackgroundTasks
  t = t + dt
  iterations = iterations + 1
  Ekinetic = 0.0
  for i = 0 to Acountm1
    Ax(i) = Ax(i) + Avx(i) * dt + Aax(i)*dtsquaredover2  // compute new position
    Ay(i) = Ay(i) + Avy(i) * dt + Aay(i)*dtsquaredover2
    Avx(i) = Avx(i) + Aax(i)*dtover2  // update velocity halfway
    Avy(i) = Avy(i) + Aay(i)*dtover2
  next
  for i = 0 to Bcountm1
    Bx(i) = Bx(i) + Bvx(i) * dt + Bax(i)*dtsquaredover2
    By(i) = By(i) + Bvy(i) * dt + Bay(i)*dtsquaredover2
    Bvx(i) = Bvx(i) + Bax(i)*dtover2
    Bvy(i) = Bvy(i) + Bay(i)*dtover2
  next
  ComputeAccelerations
  for i = 0 to Acountm1
    tempvx = Avx(i) + Aax(i)*dtover2  // finish updating v with new acceleration
    tempvy = Avy(i) + Aay(i)*dtover2
    Avx(i) = tempvx  // (temp variables eliminate multiple array lookups)
    Avy(i) = tempvy
    drawA(i)
    Ekinetic = Ekinetic + Amass * (tempvx*tempvx + tempvy*tempvy)  // we'll divide by 2 later
  next
  for i = 0 to Bcountm1
    tempvx = Bvx(i) + Bax(i)*dtover2  // finish updating v with new acceleration
    tempvy = Bvy(i) + Bay(i)*dtover2
    Bvx(i) = tempvx
    Bvy(i) = tempvy
    drawB(i)
    Ekinetic = Ekinetic + Bmass * (tempvx*tempvx + tempvy*tempvy)
  next
  Ekinetic = Ekinetic / 2.0
  EkSquared = Ekinetic * Ekinetic
  EkTotal = EkTotal + Ekinetic
  EkSquaredTotal = EkSquaredTotal + EkSquared
  EiTotal = EiTotal + Einteraction
  EgTotal = EgTotal + Egravitational
  Ptotal = Ptotal + pressure
End Sub

App.SpaceAvailable:
Function SpaceAvailable(x as single, y as single, radius as single) As boolean
  // determines whether there's space available to put a new molecule at (x,y)
  dim i as integer
  dim dx, dy as single
  dim r2limit as single
  dim tooclose as boolean
  #pragma disableBackgroundTasks
  tooclose = false
  i = 0
  r2limit = (radius + Aradius) * (radius + Aradius) * 1.26
  while i < Acount and not tooclose
    dx = Ax(i) - x
    dy = Ay(i) - y
    if dx*dx + dy*dy < r2limit then
      tooclose = true
    end if
    i = i + 1
  wend
  i = 0
  r2limit = (radius + Bradius) * (radius + Bradius) * 1.26
  while i < Bcount and not tooclose
    dx = Bx(i) - x
    dy = By(i) - y
    if dx*dx + dy*dy < r2limit then
      tooclose = true
    end if
    i = i + 1
  wend
  return not tooclose
End Function

// Event handlers:

App.Activate:
Sub Activate()
  if WindowHideControls.text = "Hide Controls" then  // if CWin is supposed to be visible...
    CWin.show
  end if
End Sub

App.Deactivate:
Sub Deactivate()
  CWin.hide
End Sub

App.EnableMenuItems:
Sub EnableMenuItems()
  AppleAboutMolecules.enable
  FileNew.enable
  FileOpen.enable
  FileSave.enable
  FileSavePicture.enable
  If CWin.DataBox.SelCount > 0 then
    EditCopy.Enable
  End if
  If CWin.DataBox.SelCount < 6 then
    EditSelectAll.Enable
  End if
  ColorsRedBlueBeige.enable
  ColorsPurpleGreenSky.enable
  ColorsRedGreenBlue.enable
  ColorsSpring.enable
  ColorsSummer.enable
  ColorsFall.enable
  ColorsWinter.enable
  ColorsChristmas.enable
  ColorsIndependence.enable
  ColorsHalloween.enable
  ColorsEarthtones.enable
  ColorsBubbles.enable
  ColorsFluorescent.enable
  ColorsGreens.enable
  ColorsBlues.enable
  ColorsPurples.enable
  ColorsReds.enable
  ColorsBrowns.enable
  ColorsGrays.enable
  WindowFullScreen.enable
  WindowHideControls.enable
End Sub

App.Open:
Sub Open()
  #if TargetCarbon
    declare function GetWindowPort lib "CarbonLib" (window as WindowPtr) as integer
    declare function NewRgn lib "CarbonLib" as integer
    declare sub SetRectRgn lib "CarbonLib" (region as integer, left as short, top as short, right as short, bottom as short)
  #endif
  dim i as integer
  #pragma disableAutoWaitCursor  // affects the whole program
  #pragma disableBackgroundTasks  // just for loops in this subroutine
  GWin.left = screen(0).width/2 - 306
  GWin.top = 50
  CWin.left = GWin.left + 412
  CWin.top = GWin.top
  #if TargetPPC
    CWin.HasBackColor = True  // gray looks better in OS9
    CWin.Backcolor = RGB(238,238,238)
  #endif
  CWin.show
  xmax = GWin.width / 10.0
  ymax = GWin.height / 10.0
  Acount = 0
  Acountm1 = Acount - 1
  Bcount = 0
  Bcountm1 = Bcount - 1
  Acolor = RGB(130,0,200)
  Bcolor = RGB(0,160,100)
  BGcolor = RGB(245,245,255)
  ColorsChanged
  ColorsPurpleGreenSky.checked = true
  setAdiameter(1.0)
  setAmass(1.0)
  setAAattraction(1.0)
  setBdiameter(1.0)
  setBmass(1.0)
  setBBattraction(1.0)
  setWallStiffness(50.0)
  setABattraction(1.0)
  for i = 0 to 999
    Ascreenx(i) = -100  // put all molecules way off screen so they'll be drawn the first time
    Bscreenx(i) = -100
  next
  GWin.show
  #if TargetCarbon
    GWinPort = GetWindowPort(GWin)
    GWinRegion = NewRgn()
    SetRectRgn GWinRegion,0,0,GWin.width,GWin.height
  #endif
  for i = 1 to 50
    createA
  next
  setGravity(0.0)
  setdt(0.025)
  ResetData
End Sub

App.OpenDocument:
Sub OpenDocument(item As FolderItem)
  OpenTheFile(item)
End Sub

// Menu handlers:

App.AppleAboutMolecules:
Function Action as Boolean
  AboutWindow.BackColor = BGcolor
  AboutWindow.StaticText1.TextColor = Acolor
  AboutWindow.StaticText2.TextColor = Bcolor
  AboutWindow.StaticText3.TextColor = Acolor
  AboutWindow.StaticText4.TextColor = Acolor
  AboutWindow.StaticText5.TextColor = Acolor
  AboutWindow.StaticText6.TextColor = Acolor
  AboutWindow.StaticText7.TextColor = Acolor
  AboutWindow.StaticText8.TextColor = Acolor
  AboutWindow.Show
  return true
End Function

App.ColorsBlues:
Function Action as Boolean
  Acolor = RGB(0,0,110)
  Bcolor = RGB(100,75,230)
  BGcolor = RGB(232,232,255)
  ColorsChanged
  ColorsBlues.checked = true
End Function

App.ColorsBrowns:
Function Action as Boolean
  Acolor = RGB(115,30,0)
  Bcolor = RGB(180,125,60)
  BGcolor = RGB(245,235,200)
  ColorsChanged
  ColorsBrowns.checked = true
End Function

App.ColorsBubbles:
Function Action as Boolean
  Acolor = RGB(130,255,255)
  Bcolor = RGB(255,200,255)
  BGcolor = RGB(0,0,210)
  ColorsChanged
  ColorsBubbles.checked = true
End Function

App.ColorsChristmas:
Function Action as Boolean
  Acolor = RGB(0,190,0)
  Bcolor = RGB(240,0,0)
  BGcolor = RGB(255,255,255)
  ColorsChanged
  ColorsChristmas.checked = true
End Function

App.ColorsEarthtones:
Function Action as Boolean
  Acolor = RGB(170,100,0)
  Bcolor = RGB(30,170,50)
  BGcolor = RGB(255,255,230)
  ColorsChanged
  ColorsEarthtones.checked = true
End Function

App.ColorsFall:
Function Action as Boolean
  Acolor = RGB(255,100,0)
  Bcolor = RGB(255,230,0)
  BGcolor = RGB(65,65,30)
  ColorsChanged
  ColorsFall.checked = true
End Function

App.ColorsFluorescent:
Function Action as Boolean
  Acolor = RGB(255,0,255)
  Bcolor = RGB(0,240,200)
  BGcolor = RGB(0,0,40)
  ColorsChanged
  ColorsFluorescent.checked = true
End Function

App.ColorsGrays:
Function Action as Boolean
  Acolor = RGB(35,35,35)
  Bcolor = RGB(115,115,115)
  BGcolor = RGB(225,225,225)
  ColorsChanged
  ColorsGrays.checked = true
End Function

App.ColorsGreens:
Function Action as Boolean
  Acolor = RGB(20,100,0)
  Bcolor = RGB(20,210,20)
  BGcolor = RGB(232,255,232)
  ColorsChanged
  ColorsGreens.checked = true
End Function

App.ColorsHalloween:
Function Action as Boolean
  Acolor = RGB(255,190,0)
  Bcolor = RGB(120,255,0)
  BGcolor = RGB(0,0,0)
  ColorsChanged
  ColorsHalloween.checked = true
End Function

App.ColorsIndependence:
Function Action as Boolean
  Acolor = RGB(240,0,0)
  Bcolor = RGB(0,0,210)
  BGcolor = RGB(255,255,255)
  ColorsChanged
  ColorsIndependence.checked = true
End Function

App.ColorsPurpleGreenSky:
Function Action as Boolean
  Acolor = RGB(130,0,200)
  Bcolor = RGB(0,160,100)
  BGcolor = RGB(240,240,255)
  ColorsChanged
  ColorsPurpleGreenSky.checked = true
End Function

App.ColorsPurples:
Function Action as Boolean
  Acolor = RGB(90,0,130)
  Bcolor = RGB(175,80,185)
  BGcolor = RGB(245,230,255)
  ColorsChanged
  ColorsPurples.checked = true
End Function

App.ColorsRedBlueBeige:
Function Action as Boolean
  Acolor = RGB(200,50,50)
  Bcolor = RGB(50,50,200)
  BGcolor = RGB(255,255,230)
  ColorsChanged
  ColorsRedBlueBeige.checked = true
End Function

App.ColorsRedGreenBlue:
Function Action as Boolean
  Acolor = RGB(255,0,0)
  Bcolor = RGB(0,255,0)
  BGcolor = RGB(0,0,130)
  ColorsChanged
  ColorsRedGreenBlue.checked = true
End Function

App.ColorsReds:
Function Action as Boolean
  Acolor = RGB(130,5,5)
  Bcolor = RGB(230,0,0)
  BGcolor = RGB(255,232,232)
  ColorsChanged
  ColorsReds.checked = true
End Function

App.ColorsSpring:
Function Action as Boolean
  Acolor = RGB(255,190,255)
  Bcolor = RGB(255,255,0)
  BGcolor = RGB(0,200,0)
  ColorsChanged
  ColorsSpring.checked = true
End Function

App.ColorsSummer:
Function Action as Boolean
  Acolor = RGB(80,0,255)
  Bcolor = RGB(0,180,0)
  BGcolor = RGB(220,232,255)
  ColorsChanged
  ColorsSummer.checked = true
End Function

App.ColorsWinter:
Function Action as Boolean
  Acolor = RGB(115,160,255)
  Bcolor = RGB(150,50,255)
  BGcolor = RGB(255,255,255)
  ColorsChanged
  ColorsWinter.checked = true
End Function

App.EditCopy:
Function Action as Boolean
  dim RowNum as integer
  dim c as Clipboard
  dim copyList as String
  c = new Clipboard
  for RowNum = 0 to 5
    if CWin.DataBox.Selected(RowNum) then
      copyList = copyList + CWin.DataBox.List(RowNum)+Chr(13)
    end if
  next
  c.SetText copyList
  c.Close
  return true
End Function

App.EditSelectAll:
Function Action as Boolean
  dim RowNum as integer
  for RowNum = 0 to 5
    CWin.DataBox.Selected(RowNum) = True
  next
  return true
End Function

App.FileNew:
Function Action as Boolean
  // delete all molecules and recreate them with random positions and velocities
  // (especially useful in starting over after calculation becomes unstable)
  dim i, tempAcount, tempBcount as integer
  tempAcount = Acount
  tempBcount = Bcount
  Acount = 0
  Acountm1 = -1
  Bcount = 0
  Bcountm1 = -1
  redraw
  for i = 0 to tempAcount - 1
    createA
  next
  for i = 0 to tempBcount - 1
    createB
  next
  return true
End Function

App.FileOpen:
Function Action as Boolean
  Dim f as FolderItem
  #if TargetWin32
    CWin.Hide
    f = GetOpenFolderItem("Molecules Save File;All Types (*.*)")
    if WindowHideControls.text = "Hide Controls" then
    CWin.show
    end if
  #else
    f = GetOpenFolderItem("Molecules Save File")
  #endif
  If f <> Nil then
    OpenTheFile(f)
  End if
  return true
End Function

App.FileSave:
Function Action as Boolean
  Dim f as FolderItem
  Dim stream as BinaryStream
  Dim i as Integer
  #pragma disableBackgroundTasks  //  so the for loops won't take forever
  #if TargetWin32
    CWin.Hide
    f = GetSaveFolderItem("Molecules Save File","Untitled.mol2")
    if WindowHideControls.text = "Hide Controls" then
    CWin.show
    end if
  #else
    f = GetSaveFolderItem("Molecules Save File","Untitled")
  #endif
  If f <> Nil then
    stream = f.CreateBinaryFile("Molecules Save File")
    stream.WriteByte(77)  // ASCII "M" to identify MD program
    stream.WriteByte(68)  // ASCII "D" to identify MD program
    stream.WriteByte(1)  // main version number
    stream.WriteByte(0)  // sub-version number
    stream.WriteSingle(xmax)
    stream.WriteSingle(ymax)
    stream.WriteSingle(gravity)
    stream.WriteSingle(dt)
    stream.WriteShort(Acount)
    stream.WriteSingle(Adiameter)
    stream.WriteSingle(Amass)
    stream.WriteSingle(AAattraction)
    stream.WriteByte(Acolor.red)
    stream.WriteByte(Acolor.green)
    stream.WriteByte(Acolor.blue)
    stream.WriteShort(Bcount)
    stream.WriteSingle(Bdiameter)
    stream.WriteSingle(Bmass)
    stream.WriteSingle(BBattraction)
    stream.WriteByte(Bcolor.red)
    stream.WriteByte(Bcolor.green)
    stream.WriteByte(Bcolor.blue)
    stream.WriteSingle(ABattraction)
    stream.WriteSingle(WallStiffness)
    stream.WriteByte(BGcolor.red)
    stream.WriteByte(BGcolor.green)
    stream.WriteByte(BGcolor.blue)
    stream.WriteLong(iterations)
    stream.WriteSingle(t)
    stream.WriteSingle(Ekinetic)
    stream.WriteSingle(Einteraction)
    stream.WriteSingle(Egravitational)
    stream.WriteSingle(EkTotal)
    stream.WriteSingle(EkSquaredTotal)
    stream.WriteSingle(EiTotal)
    stream.WriteSingle(EgTotal)
    stream.WriteSingle(pressure)
    stream.WriteSingle(pTotal)
    for i = 0 to Acountm1
      stream.WriteSingle(Ax(i))
      stream.WriteSingle(Ay(i))
      stream.WriteSingle(Avx(i))
      stream.WriteSingle(Avy(i))
      stream.WriteSingle(Aax(i))
      stream.WriteSingle(Aay(i))
    next
    for i = 0 to Bcountm1
      stream.WriteSingle(Bx(i))
      stream.WriteSingle(By(i))
      stream.WriteSingle(Bvx(i))
      stream.WriteSingle(Bvy(i))
      stream.WriteSingle(Bax(i))
      stream.WriteSingle(Bay(i))
    next
    stream.close
  end if
  return true
End Function

App.FileSavePicture:
Function Action as Boolean
  dim f as FolderItem
  dim offscreenimage as Picture  // for some reason we can't just use GWin.Graphics
  dim i,j,x,y as Integer
  #pragma disableBackgroundTasks
  #if TargetMacOS
    f = GetSaveFolderItem("image/x-pict","Untitled")
  #else
    CWin.hide
    f = GetSaveFolderItem("image/x-bmp","Untitled.bmp")
    if WindowHideControls.text = "Hide Controls" then
    CWin.show
    end if
  #endif
  if f <> nil then
    offscreenimage = NewPicture(GWin.width,GWin.height,16)  // bitdepth is 16 to save space
    offscreenimage.graphics.ForeColor = BGcolor
    offscreenimage.graphics.FillRect 0,0,GWin.width,GWin.height
    offscreenimage.graphics.ForeColor = Acolor
    for i = 0 to Acountm1
      offscreenimage.graphics.FillOval Ascreenx(i), Ascreeny(i), Apixelwidth, Apixelwidth
    next
    offscreenimage.graphics.ForeColor = Bcolor
    for i = 0 to Bcountm1
      offscreenimage.graphics.FillOval Bscreenx(i), Bscreeny(i), Bpixelwidth, Bpixelwidth
    next
    f.SaveAsPicture offscreenimage,4  // the 4 gives raster PICT on Mac, BMP on Windows
  end if
  return true
End Function

App.WindowFullScreen:
Function Action as Boolean
  if WindowFullScreen.checked then
    WindowFullScreen.checked = false
    GWin.FullScreen = false
    GWin.MenuBarVisible = true
    GWin.top = 50
    GWin.left = 10
    GWin.width = screen(0).width - 20
    GWin.height = screen(0).height - 80
  else
    WindowFullScreen.checked = true
    GWin.left = 0
    GWin.top = 0
    GWin.width = screen(0).width
    GWin.height = screen(0).height
    GWin.MenuBarVisible = false
    GWin.FullScreen = true
    CWin.left = screen(0).width - 208
  end if
  return true
End Function

App.WindowHideControls:
Function Action as Boolean
  if WindowHideControls.text = "Hide Controls" then
    WindowHideControls.text = "Show Controls"
    CWin.hide
  else
    WindowHideControls.text = "Hide Controls"
    CWin.show
  end if
  return true
End Function

//***********************************************************************

// GWin Properties:

DraggedAtom as Integer

// Event handlers:

GWin.Close:
Sub Close()
  quit // under Windows it seems i can't get rid of the close box,
  // so we'll just quit if they click it
End Sub

GWin.KeyDown:
Function KeyDown(Key As String) As Boolean
  select case key
  case "S"
    CWin.StartButton.Push
    return true
  case "H"
    if Keyboard.ShiftKey then
      App.ChangeSpeeds(2.0)
    else
      App.ChangeSpeeds(1.0/0.95)
    end if
    return true
  case "C"
    if Keyboard.ShiftKey then
      App.ChangeSpeeds(0.5)
    else
      App.ChangeSpeeds(0.95)
    end if
    return true
  end select
End Function

GWin.MouseDown:
Function MouseDown(X As Integer, Y As Integer) As Boolean
  dim i as integer
  dim pixelr, pixelr2 as integer
  dim dx, dy as integer
  pixelr = App.Apixelwidth / 2
  pixelr2 = pixelr * pixelr
  DraggedAtom = -1
  i = 0
  while (DraggedAtom = -1) and (i < App.Acount)
    dx = x - (App.Ascreenx(i) + pixelr)
    dy = y - (App.Ascreeny(i) + pixelr)
    if dx*dx + dy*dy < pixelr2 then
      DraggedAtom = i
    end if
    i = i + 1
  wend
  pixelr = App.Bpixelwidth / 2
  pixelr2 = pixelr * pixelr
  i = 0
  while (DraggedAtom = -1) and (i < App.Bcount)
    dx = x - (App.Bscreenx(i) + pixelr)
    dy = y - (App.Bscreeny(i) + pixelr)
    if dx*dx + dy*dy < pixelr2 then
      DraggedAtom = i + 1000  // kludge:  add 1000 to denote type B atom
    end if
    i = i + 1
  wend
  if DraggedAtom = -1 then
    return false
  else
    return true
  end if
End Function

GWin.MouseDrag:
Sub MouseDrag(X As Integer, Y As Integer)
  if DraggedAtom < 1000 then
    App.Ax(DraggedAtom) = x / 10.0
    App.Ay(DraggedAtom) = y / 10.0
    App.drawA(DraggedAtom)
  else
    App.Bx(DraggedAtom-1000) = x / 10.0
    App.By(DraggedAtom-1000) = y / 10.0
    App.drawB(DraggedAtom-1000)
  end if
End Sub

GWin.MouseMove:
Sub MouseMove(X As Integer, Y As Integer)
  #if TargetMacOS
    if (x > me.width - 16) and (y > me.height - 16) and (not me.FullScreen) then  // if it's near the bottom-right corner...
    ResizeWindow.left = me.left  // show an indentical-looking window with a grow box (quite a kludge, eh?)
    ResizeWindow.top = me.top
    ResizeWindow.width = me.width
    ResizeWindow.height = me.height
    if CWin.simTimer.mode = 2 then
      CWin.startButton.push
    end if
    ResizeWindow.show
    ResizeWindow.Timer1.mode = 1
    else
    if (not me.MenuBarVisible) and (y < 10) then
      me.MenuBarVisible = true
    else
      if me.FullScreen and me.MenuBarVisible and (y > 10) then
        me.MenuBarVisible = false
      end if
    end if
    end if
  #else  // no such kludge needed in Windows, since grow box doesn't intrude into window content
    if (not me.MenuBarVisible) and (y < 10) then
    me.MenuBarVisible = true
    else
    if me.FullScreen and me.MenuBarVisible and (y > 10) then
      me.MenuBarVisible = false
    end if
    end if
  #endif
End Sub

GWin.MouseUp:
Sub MouseUp(X As Integer, Y As Integer)
  // fortunately, this event will occur even if the mouse has left the window
  dim pixelr as integer
  if DraggedAtom < 1000 then
    pixelr = App.Apixelwidth / 2
    if (x < pixelr) or (x > me.width - pixelr) or (y < pixelr) or (y > me.height - pixelr) then
      App.deleteA(DraggedAtom)
    end if
  else
    pixelr = App.Bpixelwidth / 2
    if (x < pixelr) or (x > me.width - pixelr) or (y < pixelr) or (y > me.height - pixelr) then
      App.deleteB(DraggedAtom-1000)
    end if
  end if
  App.ResetData
End Sub

GWin.Moved:
Sub Moved()
  if me.top >= 40 then
    CWin.top = me.top
  end if
  if (me.width + 212 < screen(0).width) and (not me.Fullscreen) then
    if me.left + (me.width/2) <= screen(0).width/2 then
      CWin.left = me.left + me.width + 12
    else
      CWin.left = me.left - 212
    end if
  end if
End Sub

GWin.Paint:
Sub Paint(g As Graphics)
  App.redraw
End Sub

GWin.Resized:
Sub Resized()
  #if TargetCarbon
    declare sub SetRectRgn lib "CarbonLib" (region as integer, left as short, top as short, right as short, bottom as short)
  #endif
  dim i as integer
  App.xmax = me.width / 10.0
  CWin.WidthEdit.text = Format(App.xmax,"#.0")
  App.ymax = me.height / 10.0
  CWin.HeightEdit.text = Format(App.ymax,"#.0")
  #if TargetCarbon
    SetRectRgn App.GWinRegion, 0, 0, me.width, me.height
  #endif
  for i = App.Acount-1 downto 0  // gotta go backwards because "deleteA" moves 'em down the list
    if (App.Ax(i) > App.xmax - App.Aradius) or (App.Ay(i) > App.ymax - App.Aradius) then
      App.deleteA(i)
    end if
  next
  for i = App.Bcount-1 downto 0
    if (App.Bx(i) > App.xmax - App.Bradius) or (App.By(i) > App.ymax - App.Bradius) then
      App.deleteB(i)
    end if
  next
End Sub

//***********************************************************************

// CWin event handlers:

CWin.Close:
Sub Close()
  WindowHideControls.text = "Show Controls"
End Sub

CWin.KeyDown:
Function KeyDown(Key As String) As Boolean
  select case key
  case "S"
    CWin.StartButton.Push
    return true
  case "H"
    if Keyboard.ShiftKey then
      App.ChangeSpeeds(2.0)
    else
      App.ChangeSpeeds(1.0/0.95)
    end if
    return true
  case "C"
    if Keyboard.ShiftKey then
      App.ChangeSpeeds(0.5)
    else
      App.ChangeSpeeds(0.95)
    end if
    return true
  end select
End Function

CWin.MouseDown:
Function MouseDown(X As Integer, Y As Integer) As Boolean
  me.SetFocus  // remove focus from controls
End Function

// Controls:

CWin.StartButton.Action:
Sub Action()
  if self.simTimer.mode = 0 then
    self.simTimer.mode = 2
    me.caption = "Stop"
  else
    self.simTimer.mode = 0
    me.caption = "Start"
  end if
End Sub

CWin.Tarrows.Up:
Sub Up()
  if Keyboard.AsyncShiftKey then
    App.ChangeSpeeds(2.0)
  else
    App.ChangeSpeeds(1.0/0.95)
  end if
End Sub

CWin.Tarrows.Down:
Sub Down()
  if Keyboard.AsyncShiftKey then
    App.ChangeSpeeds(0.5)
  else
    App.ChangeSpeeds(0.95)
  end if
End Sub

CWin.TabPanel1.MouseDown:
Function MouseDown(X As Integer, Y As Integer) As Boolean
  dim i as integer
  for i = 0 to 6
    self.DataBox.Selected(i) = false
  next
End Function

CWin.gEdit.LostFocus:
Sub LostFocus()
  App.setGravity(round(1000*CDbl(me.text))/1000)
End Sub

CWin.gEdit.KeyDown:
Function KeyDown(Key As String) As Boolean
  if key = chr(13) then
    self.setFocus
  end if
End Function

CWin.Garrows.Down:
Sub Down()
  if Keyboard.AsyncShiftKey then
    App.setGravity(App.gravity - 0.01)
  else
    App.setGravity(App.gravity - 0.001)
  end if
End Sub

CWin.Garrows.Up:
Sub Up()
  if Keyboard.AsyncShiftKey then
    App.setGravity(App.gravity + 0.01)
  else
    App.setGravity(App.gravity + 0.001)
  end if
End Sub

CWin.dtEdit.KeyDown:
Function KeyDown(Key As String) As Boolean
  if key = chr(13) then
    self.setFocus
  end if
End Function

CWin.dtEdit.LostFocus:
Sub LostFocus()
  // i may want to go to a fourth decimal place here...
  App.setdt(round(10000.0*CDbl(me.text))/10000.0)
End Sub

CWin.dtArrows.Down:
Sub Down()
  if Keyboard.AsyncShiftKey then
    App.setdt(App.dt - 0.01)
  else
    App.setdt(App.dt - 0.001)
  end if
End Sub

CWin.dtArrows.Up:
Sub Up()
  if Keyboard.AsyncShiftKey then
    App.setdt(App.dt + 0.01)
  else
    App.setdt(App.dt + 0.001)
  end if
End Sub

CWin.ReverseButton.Action:
Sub Action()
  dim i as integer
  #pragma disableBackgroundTasks
  for i = 0 to App.Acount - 1
    App.Avx(i) = -App.Avx(i)
    App.Avy(i) = -App.Avy(i)
  next
End Sub

CWin.ADEdit.LostFocus:
Sub LostFocus()
  App.setAdiameter(round(10.0*CDbl(me.text))/10.0)
End Sub

CWin.ADEdit.KeyDown:
Function KeyDown(Key As String) As Boolean
  if key = chr(13) then
    self.setFocus
  end if
End Function

CWin.ADarrows.Down:
Sub Down()
  if Keyboard.AsyncShiftKey then
    App.setAdiameter(App.Adiameter - 1.0)
  else
    App.setAdiameter(App.Adiameter - 0.1)
  end if
End Sub

CWin.ADarrows.Up:
Sub Up()
  App.setAdiameter(App.Adiameter + 0.1)
End Sub

CWin.AMEdit.LostFocus:
Sub LostFocus()
  App.setAmass(round(10.0*CDbl(me.text))/10.0)
End Sub

CWin.AMEdit.KeyDown:
Function KeyDown(Key As String) As Boolean
  if key = chr(13) then
    self.setFocus
  end if
End Function

CWin.AMarrows.Up:
Sub Up()
  if Keyboard.AsyncShiftKey then
    App.SetAmass(App.Amass + 1.0)
  else
    App.setAmass(App.Amass + 0.1)
  end if
End Sub

CWin.AMarrows.Down:
Sub Down()
  if Keyboard.AsyncShiftKey then
    App.SetAmass(App.Amass - 1.0)
  else
    App.setAmass(App.Amass - 0.1)
  end if
End Sub

CWin.AAEdit.LostFocus:
Sub LostFocus()
  App.setAAattraction(round(10.0*CDbl(me.text))/10.0)
End Sub

CWin.AAEdit.KeyDown:
Function KeyDown(Key As String) As Boolean
  if key = chr(13) then
    self.setFocus
  end if
End Function

CWin.AAarrows.Up:
Sub Up()
  if Keyboard.AsyncShiftKey then
    App.setAAattraction(App.AAattraction + 1.0)
  else
    App.setAAattraction(App.AAattraction + 0.1)
  end if
End Sub

CWin.AAarrows.Down:
Sub Down()
  if Keyboard.AsyncShiftKey then
    App.setAAattraction(App.AAattraction - 1.0)
  else
    App.setAAattraction(App.AAattraction - 0.1)
  end if
End Sub

CWin.AcolorButton.Action:
Sub Action()
  Dim theColor as Color
  Dim b as Boolean
  theColor = App.Acolor
  b = SelectColor(theColor,"Select a color")
  App.Acolor = theColor
  App.ColorsChanged
End Sub

CWin.simTimer.Action:
Sub Action()
  App.RunSimulation
End Sub

CWin.ANEdit.LostFocus:
Sub LostFocus()
  dim newN, i as integer
  newN = round(CDbl(me.text))
  if newN < 0 then
    newN = 0
  end if
  if newN > 1000 then
    newN = 1000
  end if
  if newN > App.Acount then
    for i = App.Acount + 1 to newN
      App.createA
    next
  else
    if newN < App.Acount then
      do
        App.deleteA(App.Acount-1)
      loop until App.Acount = newN
    end if
  end if
End Sub

CWin.ANEdit.KeyDown:
Function KeyDown(Key As String) As Boolean
  if key = chr(13) then
    self.setFocus
  end if
End Function

CWin.ANarrows.Up:
Sub Up()
  dim i as integer
  if Keyboard.AsyncShiftkey then
    for i = 1 to 10
      App.createA
    next
  else
    App.createA
  end if
End Sub

CWin.ANarrows.Down:
Sub Down()
  dim i as integer
  if Keyboard.AsyncShiftKey then
    for i = 1 to 10
      App.deleteA(App.Acount-1)
    next
  else
    App.deleteA(App.Acount-1)
  end if
End Sub

CWin.BNEdit.LostFocus:
Sub LostFocus()
  dim newN, i as integer
  newN = round(CDbl(me.text))
  if newN < 0 then
    newN = 0
  end if
  if newN > 1000 then
    newN = 1000
  end if
  if newN > App.Bcount then
    for i = App.Bcount + 1 to newN
      App.createB
    next
  else
    if newN < App.Bcount then
      do
        App.deleteB(App.Bcount-1)
      loop until App.Bcount = newN
    end if
  end if
End Sub

CWin.BNEdit.KeyDown:
Function KeyDown(Key As String) As Boolean
  if key = chr(13) then
    self.setFocus
  end if
End Function

CWin.BNarrows.Up:
Sub Up()
  dim i as integer
  if Keyboard.AsyncShiftkey then
    for i = 1 to 10
      App.createB
    next
  else
    App.createB
  end if
End Sub

CWin.BNarrows.Down:
Sub Down()
  dim i as integer
  if Keyboard.AsyncShiftKey then
    for i = 1 to 10
      App.deleteB(App.Bcount-1)
    next
  else
    App.deleteB(App.Bcount-1)
  end if
End Sub

CWin.BDEdit.LostFocus:
Sub LostFocus()
  App.setBdiameter(round(10.0*CDbl(me.text))/10.0)
End Sub

CWin.BDEdit.KeyDown:
Function KeyDown(Key As String) As Boolean
  if key = chr(13) then
    self.setFocus
  end if
End Function

CWin.BDarrows.Up:
Sub Up()
  App.setBdiameter(App.Bdiameter + 0.1)
End Sub

CWin.BDarrows.Down:
Sub Down()
  if Keyboard.AsyncShiftKey then
    App.setBdiameter(App.Bdiameter - 1.0)
  else
    App.setBdiameter(App.Bdiameter - 0.1)
  end if
End Sub

CWin.BMEdit.LostFocus:
Sub LostFocus()
  App.setBmass(round(10.0*CDbl(me.text))/10.0)
End Sub

CWin.BMEdit.KeyDown:
Function KeyDown(Key As String) As Boolean
  if key = chr(13) then
    self.setFocus
  end if
End Function

CWin.BMarrows.Up:
Sub Up()
  if Keyboard.AsyncShiftKey then
    App.setBmass(App.Bmass + 1.0)
  else
    App.setBmass(App.Bmass + 0.1)
  end if
End Sub

CWin.BMarrows.Down:
Sub Down()
  if Keyboard.AsyncShiftKey then
    App.setBmass(App.Bmass - 1.0)
  else
    App.setBmass(App.Bmass - 0.1)
  end if
End Sub

CWin.BBEdit.LostFocus:
Sub LostFocus()
  App.setBBattraction(round(10.0*CDbl(me.text))/10.0)
End Sub

CWin.BBEdit.KeyDown:
Function KeyDown(Key As String) As Boolean
  if key = chr(13) then
    self.setFocus
  end if
End Function

CWin.BBarrows.Up:
Sub Up()
  if Keyboard.AsyncShiftKey then
    App.setBBattraction(App.BBattraction + 1.0)
  else
    App.setBBattraction(App.BBattraction + 0.1)
  end if
End Sub

CWin.BBarrows.Down:
Sub Down()
  if Keyboard.AsyncShiftKey then
    App.setBBattraction(App.BBattraction - 1.0)
  else
    App.setBBattraction(App.BBattraction - 0.1)
  end if
End Sub

CWin.ABEdit.LostFocus:
Sub LostFocus()
  App.setABattraction(round(10.0*CDbl(me.text))/10.0)
End Sub

CWin.ABEdit.KeyDown:
Function KeyDown(Key As String) As Boolean
  if key = chr(13) then
    self.setFocus
  end if
End Function

CWin.ABarrows.Up:
Sub Up()
  if Keyboard.AsyncShiftKey then
    App.setABattraction(App.ABattraction + 1.0)
  else
    App.setABattraction(App.ABattraction + 0.1)
  end if
End Sub

CWin.ABarrows.Down:
Sub Down()
  if Keyboard.AsyncShiftKey then
    App.setABattraction(App.ABattraction - 1.0)
  else
    App.setABattraction(App.ABattraction - 0.1)
  end if
End Sub

CWin.BcolorButton.Action:
Sub Action()
  Dim theColor as Color
  Dim b as Boolean
  theColor = App.Bcolor
  b = SelectColor(theColor,"Select a color")
  App.Bcolor = theColor
  App.ColorsChanged
End Sub

CWin.BGcolorButton.Action:
Sub Action()
  Dim theColor as Color
  Dim b as Boolean
  theColor = App.BGcolor
  b = SelectColor(theColor,"Select a color")
  App.BGcolor = theColor
  App.ColorsChanged
End Sub

CWin.DataBox.Open:
Sub Open()
  dim i as integer
  for i = 0 to 5
    me.CellAlignment(i,0) = 3
  next
End Sub

CWin.InstButton.Action:
Sub Action()
  self.AvButton.value = not me.value
  App.DisplayData
End Sub

CWin.AvButton.Action:
Sub Action()
  self.InstButton.value = not me.value
  App.DisplayData
End Sub

CWin.ResetButton.Action:
Sub Action()
  App.ResetData
End Sub

CWin.WidthEdit.LostFocus:
Sub LostFocus()
  dim newwidth as single
  newwidth = CDbl(me.text)
  if newwidth < 5.0 then
    newwidth = 5.0
  else
    if newwidth > screen(0).width / 10.0 then
      newwidth = screen(0).width / 10.0
    end if
  end if
  GWin.fullScreen = false
  GWin.menuBarVisible = true
  WindowFullScreen.checked = false
  if GWin.top < 50 then
    GWin.top = 50
  end if
  GWin.width = round(10.0 * newwidth)
End Sub

CWin.WidthEdit.KeyDown:
Function KeyDown(Key As String) As Boolean
  if key = chr(13) then
    self.setFocus
  end if
End Function

CWin.HeightEdit.KeyDown:
Function KeyDown(Key As String) As Boolean
  if key = chr(13) then
    self.setFocus
  end if
End Function

CWin.HeightEdit.LostFocus:
Sub LostFocus()
  dim newheight as single
  newheight = CDbl(me.text)
  if newheight < 5.0 then
    newheight = 5.0
  else
    if newheight > screen(0).height / 10.0 then
      newheight = screen(0).height / 10.0
    end if
  end if
  GWin.fullScreen = false
  GWin.menuBarVisible = true
  WindowFullScreen.checked = false
  if GWin.top < 50 then
    GWin.top = 50
  end if
  GWin.height = round(10.0 * newheight)
End Sub

CWin.WSEdit.LostFocus:
Sub LostFocus()
  App.setWallStiffness(round(CDbl(me.text)))
End Sub

CWin.WSEdit.KeyDown:
Function KeyDown(Key As String) As Boolean
  if key = chr(13) then
    self.setFocus
  end if
End Function

CWin.WSArrows.Up:
Sub Up()
  If Keyboard.AsyncShiftkey then
    App.setWallStiffness(App.WallStiffness + 10.0)
  else
    App.setWallStiffness(App.WallStiffness + 1.0)
  end if
End Sub

CWin.WSArrows.Down:
Sub Down()
  If Keyboard.AsyncShiftkey then
    App.setWallStiffness(App.WallStiffness - 10.0)
  else
    App.setWallStiffness(App.WallStiffness - 1.0)
  end if
End Sub

CWin.StepButton.Action:
Sub Action()
  dim i as integer
  if Keyboard.ShiftKey then
    for i = 1 to 10
      App.SingleStep
    next
  else
    App.SingleStep
  end if
  App.DisplayData
End Sub

//***********************************************************************

// ResizeWindow note:

This window is a kludge to allow the graphics window (GWin) to be resized
without always showing the grow icon, which would intrude into the content.
On Windows OS this isn't necessary, so this window can be deleted from the 
source code; simply check the GrowIcon box in the GWin properties window.

// ResizeWindow event handlers:

ResizeWindow.MouseEnter:
Sub MouseEnter()
  me.hide
End Sub

ResizeWindow.Paint:
Sub Paint(g As Graphics)
  dim i as integer
  #pragma disableBackgroundTasks
  me.graphics.forecolor = App.BGcolor
  me.graphics.fillrect 0,0,me.width, me.height
  me.graphics.forecolor = App.Acolor
  for i = 0 to App.Acountm1
    me.graphics.filloval App.Ascreenx(i),App.Ascreeny(i),App.Apixelwidth,App.Apixelwidth
  next
  me.graphics.forecolor = App.Bcolor
  for i = 0 to App.Bcountm1
    me.graphics.filloval App.Bscreenx(i),App.Bscreeny(i),App.Bpixelwidth,App.Bpixelwidth
  next
End Sub

ResizeWindow.Resized:
Sub Resized()
  GWin.width = me.width
  GWin.height = me.height
End Sub

// Controls:

ResizeWindow.Timer1.Action:
Sub Action()
  self.hide
End Sub

//***********************************************************************

// AboutWindow text content:

Molecules    version 1.0

This program simulates the interactions of simple molecules (atoms) that 
attract each other when they are close but strongly repel when they touch.  
By changing the energy and density you can study the gas, liquid, and 
solid phases.

Copyright 2002 by Dan Schroeder, Weber State University, Ogden, UT, USA.  
http://physics.weber.edu/schroeder/.

Click the Next button below for more information about the simulation, 
using this program, and copying this program.

Each pair of molecules in this two-dimensional simulation interacts via 
the potential energy function

     u(r) = 4*((d/r)^12 - (d/r)^6),

where r is the distance between their centers and d is the sum of their 
two radii.  (This function is known as the Lennard-Jones potential.)  
The force in the radial direction is minus the r-derivative of this 
function; it is repulsive out to a distance slightly greater than r, 
and attractive beyond that.

When the "attraction" parameter is set to a value other than 1, the force 
is unchanged in the repulsive region but multiplied by the parameter in 
the attractive region.

This program uses natural units based on the default diameter, mass, 
and potential-well depth.  One unit of distance is represented by 10 pixels 
on the screen.  The unit of time is such that a unit-mass molecule moving 
one unit of distance per unit of time has a kinetic energy of 1/2.  

A unit of energy divided by a unit of distance gives the standard unit 
of force.  This then defines the units of the gravitational constant 
(force/mass), wall stiffness (force/distance), and pressure (also 
force/distance in two dimensions).  Temperature is measured in units 
of energy (i.e., Boltzmann's constant equals 1).  For this classical, 
two-dimensional system, the temperature is simply equal to the average 
kinetic energy per molecule.

Hints and Shortcuts

Use the S key to start or stop the simulation.  Use H and C to add and 
remove energy.  The shift key will augment the functions of most of the 
arrow controls.

If the time step is too large, the program will not treat collisions 
accurately and the total energy will increase, producing a runaway effect.  
If this happens, you can use the File-New command to start over.

You can copy and paste from the data list into a spreadsheet or word 
processor.  You can drag molecules around with the mouse, but try not 
to place them so they overlap.  Similarly, don't increase the molecular 
diameter too much at once.

This program is free, open-source software.  You can redistribute it 
and/or modify it under the terms of the GNU General Public License as 
published by the Free Software Foundation, either version 2 of the 
License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but 
without any warranty; without even the implied warranty of merchantability 
or fitness for a particular purpose.  See the GNU General Public License 
for more details.

You should have received a copy of the GNU General Public License along 
with this program.  If not, contact the Free Software Foundation, Inc., 
59 Temple Place, Suite 330, Boston, MA 02111-1307, USA, www.gnu.org.

The integration algorithm used by this program, as well as several other 
features, is taken from H. Gould and J. Tobochnik, An Introduction to 
Computer Simulation Methods (2nd edition, Addison-Wesley, 1996).  I highly 
recommend this book as a source of ideas for using and modifying this 
program.

Please feel free to download the source code to this program from the 
author's web site, http://physics.weber.edu/schroeder/.

Made with REALbasic

// Event handler:

AboutWindow.MouseDown:
Function MouseDown(X As Integer, Y As Integer) As Boolean
  me.hide
End Function

// Controls:

AboutWindow.CloseButton.Action:
Sub Action()
  self.hide
End Sub

AboutWindow.NextButton.Action:
Sub Action()
  if Self.PagePanel1.Value = 5 then
    Self.PagePanel1.Value = 0
  else
    Self.PagePanel1.Value = Self.PagePanel1.Value + 1
  end if
End Sub