' gc3lib.bas - jackord@kw.igs.net - revised 28 Feb 08 - Liberty Basic v4.02
' magnetic fieldlines in the x-y plane due to
'   (1) a current-carrying loop centered in the x-z plane
'   (2) a north/south pole pair centered on the y-axis
' when "Far" is selected, the dimensions of the loop and pole pair
'   are reduced by a factor of ten below their plotted size

' Initialize Window
  nomainwin
  WindowWidth=426    ' pixel scale 0-408
  WindowHeight=374   ' pixel scale 0-336
  UpperLeftX=100: UpperLeftY=50
  button#1, "Plot", [plot], UL, 154, 5, 40, 18
  checkbox#1.ch, "Far", [waitHere], [waitHere], 214, 5, 40, 18
  open "Magnetic Dipole" for graphics_nsb as #1
  #1 "trapclose [quit]"
  d=.5: pi=4*atn(1)
  gosub [pin]

[waitHere]
  wait

[pin] ' Initialize Screen
  #1 "cls ; down ; color blue ; backcolor blue"
  #1 "line 204 156 204 180"
  #1 "place 204 156 ; circlefilled 3 ; place 204 180 ; circlefilled 3"
  #1 "color red ; backcolor red"
  #1 "place 192 167 ; boxfilled 216 170 ; backcolor white"
return

[plot] ' Plot Fieldlines
  gosub [pin]
  #1.ch "value?": input#1.ch, z$ 
  rb=6: if z$="set" then rb=rb/10
  for j=1 to 2
    for i=5 to 175 step 10
      z=i*pi/180: x=12*cos(z): y=12*sin(z)
      #1 "place "; 204+2*x; " "; 168-2*y
      xx=x: yy=y
      if j=1 then gosub [bfld] else gosub [qfld]
      dx=bx*d: dy=by*d
      while x*x+y*y>=143.9
        xx=x+dx/2: yy=y+dy/2
        if j=1 then gosub [bfld] else gosub [qfld]
        dx=bx*d: dy=by*d: x=x+dx: y=y+dy
        #1 "goto "; 204+2*x; " "; 168-2*y
      wend
      scan
    next i
    #1 "color blue"
  next j
goto [waitHere]

[quit]
  close#1
end

[bfld]                                      ' Biot-Savart law
  bx=0: by=0
  for k=1 to 60
    sn=sin(k*pi/30)
    r3=(xx*xx+yy*yy+rb*rb-2*xx*rb*sn)^1.5
    bx=bx+yy*sn/r3: by=by+(rb-xx*sn)/r3
  next k
  b=(bx*bx+by*by)^.5: bx=bx/b: by=by/b
return

[qfld]                                      ' Inverse-square law
  r=(xx*xx+(yy-rb)*(yy-rb))^1.5
  bx=xx/r: by=(yy-rb)/r
  r=(xx*xx+(yy+rb)*(yy+rb))^1.5
  bx=bx-xx/r: by=by-(yy+rb)/r
  b=sqr(bx*bx+by*by): bx=bx/b: by=by/b
return