' frdslib.bas - jackord@kw.igs.net - revised 20 Jan 04 - Liberty Basic v3.01
' gravitational force between identical rods of length D:
'   rod A lies along the x-axis with its center at (D/2, 0)
'   rod B is perpendicular to the x-axis with its center at (2*D, D/2)
' each rod is made up of N spheres of diameter D/N
' in the program, D=128, and N=8 or N=16
' the rods are released from rest and their motion is followed
'   until they collide
' the program takes a long time to execute and only the collision is plotted

' Initialize Window
  nomainwin
  button#1, "Run", [strt], UL, 10, 10, 40, 20
  button#1, "8/16", [setn], UL, 60, 10, 40, 20
  WindowWidth=310     ' pixel scale 0-300
  WindowHeight=229    ' pixel scale 0-200
  UpperLeftX=100: UpperLeftY=100
  open "Gravitational Attraction Between Rods" for graphics_nsb as #1
  #1 "trapclose [quit]"
  dim xa(17): dim ya(17): dim xb(17): dim yb(17)
  kk=0: n=8: d=128: mr=128
goto [rods]

[waitHere]
  wait

[setn]                                                ' Initialize & set n
  kk=0: n=128/n
goto [rods]

[strt]                                                ' Run if initialized
  kk=kk+1
goto [rods]

[rods]
  r=d/n/2: lim=8*r*r*r: rm=2*lim: m=mr/n: mm=m*m: dt=.1
  pi=4*atn(1): aa=0: ab=pi/2
  if kk=0 then
    #1 "backcolor white ; cls ; color red ; backcolor red ; down"
    for i=1 to n                                      ' Define and plot rods
      xa(i)=20+r*(2*i-1): ya(i)=30+128
      xb(i)=20+256: yb(i)=30+r*(2*i-1)
      #1 "place "; xa(i); " "; ya(i): #1 "circlefilled "; r
      #1 "place "; xb(i); " "; yb(i): #1 "circlefilled "; r
    next i
    #1 "line 175 126 186 126 ; line 180 121 180 132"
    Ic=.4*mr*r*r: xac=(xa(1)+xa(n))/2                 ' Moment of Inertia
    for i=1 to n/2
      Ic=Ic+2*m*(xa(i)-xac)*(xa(i)-xac)
    next i
  end if
  if kk=1 then
    gosub [force]                                     ' Feynman
    vax=fax/mr*dt/2: vay=fay/mr*dt/2                  '   half
    wa=Ta/Ic*dt/2: wb=Tb/Ic*dt/2                      '   step
    while rm>lim                                      ' Loop until collision
      xac=xac+vax*dt: yac=yac+vay*dt                  ' Translation
      xbc=xbc-vax*dt: ybc=ybc-vay*dt
      aa=aa+wa*dt: ab=ab+wb*dt                        ' Rotation
      for i=1 to n                                    ' New rod positions
        xa(i)=xac+(2*i-1-n)*r*cos(aa)
        ya(i)=yac+(2*i-1-n)*r*sin(aa)
        xb(i)=xbc+(2*i-1-n)*r*cos(ab)
        yb(i)=ybc+(2*i-1-n)*r*sin(ab)
      next i
      gosub [force]                                   ' Newton's law
      vax=vax+fax/mr*dt: vay=vay+fay/mr*dt
      wa=wa+Ta/Ic*dt: wb=wb+Tb/Ic*dt
      scan
    wend
    aa=int(aa*18000/pi)/100: ab=int(ab*18000/pi-9000)/100
    #1 "color blue ; backcolor blue"
    for i=1 to n                                      ' Rods after collision
      #1 "place "; xa(i); " "; ya(i): #1 "circlefilled "; r
      #1 "place "; xb(i); " "; yb(i): #1 "circlefilled "; r
    next i
    #1 "color black ; backcolor white"
    #1 "place 20 120": #1 "\Rotates "; 0-aa; " degrees"
    #1 "place 110 40": #1 "\Rotates "; ab; " degrees"
  end if
goto [waitHere]

[force]
  xac=(xa(1)+xa(n))/2: yac=(ya(1)+ya(n))/2            ' C of M
  xbc=(xb(1)+xb(n))/2: ybc=(yb(1)+yb(n))/2
  fax=0: fay=0: Ta=0: Tb=0
  for i=1 to n                                        ' Rod A
    fx=0: fy=0
    for j=1 to n                                      ' Rod B
      dx=xb(j)-xa(i): dy=yb(j)-ya(i)
      rcube=(dx*dx+dy*dy)^1.5
      dfx=mm*dx/rcube                                 ' Force between spheres
      dfy=mm*dy/rcube
      fx=fx+dfx: fy=fy+dfy
      Tb=Tb-(xb(j)-xbc)*dfy+(yb(j)-ybc)*dfx           ' Torque on Rod B
      if rcube<rm then rm=rcube
    next j
    fax=fax+fx: fay=fay+fy                            ' Force on Rod A
    Ta=Ta+(xa(i)-xac)*fy-(ya(i)-yac)*fx               ' Torque on Rod A
  next i
return

[quit]
  close #1
end