REM gc3bbc.bas - jackord@kw.igs.net - 27 Feb 08 - BBC BASIC (RTR)
      REM magnetic fieldlines in the x-y plane due to
      REM   (1) a current-carrying loop centered in the x-z plane
      REM   (2) a north/south pole pair centered on the y-axis
      REM when "Far" is selected, the dimensions of the loop and pole pair
      REM   are reduced by a factor of ten below their plotted size
      
      REM Use 64-bit floats
      *FLOAT 64
      
      REM Install libraries
      INSTALL @lib$+"WINLIB5"
      _BM_GETCHECK = 240
      
      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", 154, 5, 40, 18, FN_setproc(PROCplot), 0)
      chkbox% = FN_button("Far", 214, 5, 40, 18, 0, 3)
      SYS "SetWindowText", @hwnd%, "Magnetic Dipole"
      d=.5
      PROCpin
      
      REPEAT
        WAIT 1
      UNTIL FALSE
      
      DEF PROCpin : REM Initialize Screen
      CLS
      GCOL 12
      LINE 408, 312, 408, 360
      CIRCLE FILL 410, 312, 8
      CIRCLE FILL 410, 360, 8
      GCOL 9
      RECTANGLE FILL 384, 334, 48, 6
      ENDPROC
      
      DEF PROCplot : REM Plot Fieldlines
      PROCpin
      SYS "SendMessage", chkbox%, _BM_GETCHECK, 0, 0 TO z%
      rb=6: IF z% 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)
          MOVE 408+4*x, 336-4*y
          xx=x: yy=y
          IF j=1 THEN PROCbfld ELSE PROCqfld
          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 PROCbfld ELSE PROCqfld
            dx=bx*d: dy=by*d: x=x+dx: y=y+dy
            DRAW 408+4*x, 336-4*y
          ENDWHILE
        NEXT i
        GCOL 12
      NEXT j
      ENDPROC
      
      DEF PROCbfld : REM 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
      ENDPROC
      
      DEF PROCqfld : REM 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
      ENDPROC