' ga1lib.bas - jackord@kw.igs.net - updated 28 Feb 08 - Liberty Basic v4.02
' fieldlines and equipotentials for an equilateral mass array

' Initialize Window
  nomainwin
  WindowWidth=426    ' pixel scale 0-408
  WindowHeight=374   ' pixel scale 0-336
  UpperLeftX=100: UpperLeftY=50
  button#1, "Plot", [plot], UL, 6, 6, 40, 18
  open "Fieldlines and Equipotentials" for graphics_nsb as #1
  #1 "trapclose [quit]"

' Initialize and Plot Equilateral Array
  dim m(4): dim p(4): dim q(4)
  d=1: pi=4*atn(1)
  m(1)=1: m(2)=1: m(3)=1
  p(1)=121: q(1)=216: p(2)=204: q(2)=72: p(3)=287: q(3)=216
  #1 "down ; backcolor pink"
  for i=1 to 3
    #1 "place "; p(i); " "; q(i)
    #1 "circlefilled 12"
    #1 "place "; p(i)-6; " "; q(i)+5: #1 "\M"
  next i

[waitHere]
  wait

[plot]
' Plot Fieldlines
  #1 "color blue"
  for j=1 to 3
    for i=1 to 24
      t=(2*i-1)*pi/24: x=p(j)+12*cos(t): y=q(j)+12*sin(t)
      gosub [fdln]
    next i
  next j

' Plot Equipotentials
  #1 "color red"
  for i=1 to 11
    x=204: y=162+14*i
    gosub [eqpt]
  next i
goto [waitHere]

[quit]
  close#1
end

[fdln] ' Fieldline Subroutine
  #1 "place "; int(x); " "; int(y)
  tst=0: dx=0: dy=0: gosub [grav]
  dx=0-gx*d: dy=0-gy*d
  while tst=0
    dxo=dx: dyo=dy: gosub [grav]
    dx=0-gx*d: dy=0-gy*d: x=x+dx: y=y+dy
    #1 "goto "; int(x); " "; int(y)
    if x<0 or y<0 or x>408 or y>336 then tst=1
    if dxo*dx<0 and dyo*dy<0 then tst=1
  wend
return

[eqpt] ' Equipotential Subroutine
  #1 "place "; int(x); " "; int(y)
  dx=0: dy=0: gosub [grav]
  dx=gy*d: dy=0-gx*d: dxi=dx: xold=x
  while (xold-204)*dxi>=0 or (x-204)*dxi<0
    gosub [grav]
    dx=gy*d: dy=0-gx*d: xold=x: x=x+dx: y=y+dy
    #1 "goto "; int(x); " "; int(y)
  wend
return

[grav] ' Field Subroutine
  xx=x+dx/2: yy=y+dy/2: gx=0: gy=0
  for k=1 to 3
    r(k)=((xx-p(k))*(xx-p(k))+(yy-q(k))*(yy-q(k)))^1.5
    gx=gx+m(k)*(p(k)-xx)/r(k): gy=gy+m(k)*(q(k)-yy)/r(k)
  next k
  gg=(gx*gx+gy*gy)^.5: gx=gx/gg: gy=gy/gg
return