REM gc2bbc.bas - jackord@kw.igs.net - revised 27 Feb 08 - BBC BASIC (RTR)
      REM magnetic fieldlines in the x-y plane due to a current-carrying loop
      REM   centered in the x-z plane
      
      REM Use 64-bit floats
      *FLOAT 64
      
      REM Install libraries
      INSTALL @lib$+"WINLIB5"
      
      REM Initialize Window
      WindowWidth=408    : REM pixel scale 0-408
      WindowHeight=336   : REM pixel scale 0-336
      VDU 23,22,WindowWidth;WindowHeight;8,15,16,128
      button% = FN_button("Plot", 5, 5, 40, 18, FN_setproc(PROCplot), 0)
      SYS "SetWindowText", @hwnd%, "Field of a Current Loop"
      
      n=200: r=6: d=.1: dt=2*PI/n
      PROCpin
      
      REPEAT
        WAIT 1
      UNTIL FALSE
      
      DEF PROCpin : REM Initialize Screen
      CLS
      GCOL 9
      RECTANGLE FILL 264,330,290,6
      ENDPROC
      
      DEF PROCplot : REM Plot Fieldlines
      PROCpin
      GCOL 12
      FOR i=-13 TO 13 STEP 2
        x=i: y=0-14: phase=1: PROCfdln
      NEXT i
      FOR i=0-14 TO 0-2 STEP 4
        x=0-17: y=i: phase=2: PROCfdln
      NEXT i
      FOR i=0-15 TO 0-7 STEP 4
        x=i: y=0: phase=4: PROCfdln
      NEXT i
      FOR i=0-14 TO 0-2 STEP 4
        x=17: y=i: phase=3: PROCfdln
      NEXT i
      FOR i=7 TO 15 STEP 4
        x=i: y=0: phase=4: PROCfdln
      NEXT i
      ENDPROC
      
      DEF PROCfdln : REM Fieldline Subroutine
      MOVE 24*(x+17), 24*(14-y)
      z=x: dx=0: dy=0: test=0
      PROCmfld: dx=Bx*d: dy=By*d: dyi=dy
      WHILE test=0
        PROCmfld: dx=Bx*d: dy=By*d: z=y: x=x+dx: y=y+dy
        DRAW 24*(x+17), 24*(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
      ENDWHILE
      ENDPROC
      
      DEF PROCmfld : REM 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
      ENDPROC