' fcomlib.bas - jackord@kw.igs.net - revised 6 Apr 05 - Liberty Basic v4.01
' motion of two particles attracted to each other
'   by an inverse square force law
' revisions make graphics compatible with Windows XP

' Initialize Window
  nomainwin
  button#1, "LAB", [b1], UL, 10, 10, 30, 20
  button#1, "C M", [b2], UL, 50, 10, 30, 20
  WindowWidth=370     ' pixel scale 0-360
  WindowHeight=269   ' pixel scale 0-240
  UpperLeftX=10: UpperLeftY=10
  open "Two-Particle Motion" for graphics_nsb as #1   ' no slide bars
  #1 "trapclose [quit]"
  gosub [pin]

[waitHere]
  wait

[pin]                                                 ' initial screen
  #1 "rule xor ; backcolor white ; cls ; color blue ; down"
  #1 "place 12 210 ; goto 30 144 ; goto 32 150 ; goto 25 148 ; goto 30 144"
  #1 "place 72 210 ; goto 90 234 ; goto 84 233 ; goto 90 228 ; goto 90 234"
  #1 "color red ; backcolor red"
  #1 "place 12 210 ; circlefilled 5 ; circle 5"
  #1 "place 72 210 ; circlefilled 5 ; circle 5"
return

[b1]
  kk=1: goto [orbit]

[b2]
  kk=2: goto [orbit]

[orbit]
  gosub [pin]                                         ' show initial screen...
  st=time$("ms"): while (time$("ms")-st)<1000: wend   '   ...briefly
  #1 "backcolor white ; cls"                          '   ...then clear it
  k=60000: ma=2: mb=3: dt=.01: t=0: lp=0              ' masses, force constant
  vcx=3.6: vcy=0-2.4: xo=0: yo=190                    ' CofM velocity
  xa=12: ya=210: vax=0: vay=0-10.8                    ' initial positions
  xb=72: yb=210: vbx=0: vby=7.2                       '   and velocities
  if kk=1 then
    vax=vax+vcx: vay=vay+vcy: vbx=vbx+vcx: vby=vby+vcy
  else
    ya=85: yb=85                                      ' site CofM frame
    #1 "place 5 125 ; down; color black"
    #1 "\C M Frame"
  end if
  #1 "color red ; backcolor red"                      ' starting points
  #1 "place "; int(xa); " "; int(ya): #1 "circlefilled 5 ; circle 5"
  #1 "place "; int(xb); " "; int(yb): #1 "circlefilled 5 ; circle 5"
  dx=xb-xa: dy=yb-ya: r3=(dx*dx+dy*dy)^1.5
  aax=k/ma*dx/r3: aay=k/ma*dy/r3                      ' Newton's law
  abx=0-k/mb*dx/r3: aby=0-k/mb*dy/r3
  xa1=int(xa): ya1=int(ya): xb1=int(xb): yb1=int(yb)
  vax=vax+aax*dt/2: vay=vay+aay*dt/2                  ' Feynman half-step
  vbx=vbx+abx*dt/2: vby=vby+aby*dt/2
  tt=time$("ms")
  while t<75                                          ' main loop
    t=t+dt: lp=lp+1
    xa=xa+vax*dt: ya=ya+vay*dt                        ' update position
    xb=xb+vbx*dt: yb=yb+vby*dt
    xa2=int(xa): ya2=int(ya): xb2=int(xb): yb2=int(yb)
    if lp=6 then
      lp=0: scan
      #1 "place "; xa1; " "; ya1: #1 "circlefilled 5 ; circle 5"
      #1 "place "; xa2; " "; ya2: #1 "circlefilled 5 ; circle 5"
      #1 "place "; xb1; " "; yb1: #1 "circlefilled 5 ; circle 5"
      #1 "place "; xb2; " "; yb2: #1 "circlefilled 5 ; circle 5"
      #1 "color green ; rule over"
      #1 "line "; xa1; " "; ya1; " "; xa2; " "; ya2
      #1 "line "; xb1; " "; yb1; " "; xb2; " "; yb2
      xa1=xa2: ya1=ya2: xb1=xb2: yb1=yb2
      if kk=2 then
        #1 "rule over ; getbmp nb 0 60 80 50"
        xo=int(vcx*t): yo=190+int(vcy*t)
        #1 "drawbmp nb "; xo; " "; yo
      end if
      #1 "discard ; rule xor"
    end if
    dx=xb-xa: dy=yb-ya: r3=(dx*dx+dy*dy)^1.5
    aax=k/ma*dx/r3: aay=k/ma*dy/r3                    ' Newton's law
    abx=0-k/mb*dx/r3: aby=0-k/mb*dy/r3
    vax=vax+aax*dt: vay=vay+aay*dt                    ' update velocity
    vbx=vbx+abx*dt: vby=vby+aby*dt
  wend
  tt=time$("ms")-tt
  #1 "rule over ; color black ; backcolor white"
  #1 "place 5 235": #1 "\Run time(ms)="; tt
goto [waitHere]

[quit]
  close #1
end