' gc2lib.bas - jackord@kw.igs.net - revised 28 Feb 08 - Liberty Basic v4.02
' magnetic fieldlines in the x-y plane due to a current-carrying loop
'   centered in the x-z plane

' Initialize Window
  nomainwin
  WindowWidth=426    ' pixel scale 0-408
  WindowHeight=374   ' pixel scale 0-336
  UpperLeftX=100: UpperLeftY=50
  button#1, "Plot", [plot], UL, 5, 5, 40, 18
  open "Field of a Current Loop" for graphics_nsb as #1
  #1 "trapclose [quit]"
  n=200: r=6: d=.1: pi=4*atn(1): dt=2*pi/n
  gosub [pin]

[waitHere]
  wait

[pin] ' Initialize Screen
  #1 "cls ; down ; color red ; backcolor red"
  #1 "place 132 167 ; boxfilled 277 170 ; backcolor white"
return

[plot] ' Plot Fieldlines
  gosub [pin]
  #1 "color blue"
  for i=-13 to 13 step 2
    x=i: y=0-14: phase=1: gosub [fdln]
  next i
  for i=0-14 to 0-2 step 4
    x=0-17: y=i: phase=2: gosub [fdln]
  next i
  for i=0-15 to 0-7 step 4
    x=i: y=0: phase=4: gosub [fdln]
  next i
  for i=0-14 to 0-2 step 4
    x=17: y=i: phase=3: gosub [fdln]
  next i
  for i=7 to 15 step 4
    x=i: y=0: phase=4: gosub [fdln]
  next i
goto [waitHere]

[quit]
  close#1
end

[fdln] ' Fieldline Subroutine
  #1 "place "; int(12*(x+17)); " "; int(12*(14-y))
  z=x: dx=0: dy=0: test=0
  scan: gosub [mfld]: dx=Bx*d: dy=By*d: dyi=dy
  while test=0
    gosub [mfld]: dx=Bx*d: dy=By*d: z=y: x=x+dx: y=y+dy
    #1 "goto "; int(12*(x+17)); " "; int(12*(14-y))
    if phase=1 and y>14 then test=1
    if phase=2 and x<0-17 then test=1
    if phase=3 and x>17 then test=1
    if phase=4 and z*dyi<0 and y*dyi>0 then test=1
  wend
return

[mfld] ' Field Subroutine
  xx=x+dx/2: yy=y+dy/2: Bx=0: By=0
  for k=1 to n
    sn=sin(k*dt)
    r3=(xx*xx+yy*yy+r*r-2*xx*r*sn)^1.5
    Bx=Bx+yy*sn/r3: By=By+(r-xx*sn)/r3
  next k
  b=(Bx*Bx+By*By)^.5: Bx=Bx/b: By=By/b
return