PROGRAM magnetism
LIBRARY "mygraphics"
DIM x(50),y(50),z(50),dLx(50),dLy(50),dLz(50)
CALL wire(N,a,delta,x,y,z,dLx,dLy,dLz)
CALL screen(a)
DO
   GET POINT xcursor,ycursor
   CALL draw(N,delta,x,y,z,dLx,dLy,dLz,xcursor,ycursor)
LOOP until (abs(xcursor) < .1) and (abs(ycursor) < .1)
END

SUB wire(N,a,delta,x(),y(),z(),dLx(),dLy(),dLz())
    ! read in current segments
    ! assume current equals unity
    INPUT prompt "radius of loop = ": a
    INPUT prompt "number of segments = ": N
    ! magnitude of magnetic field line segment
    INPUT prompt "delta = ": delta
    LET delta_angle = 2*pi/N
    LET angle = 0
    LET dL = 2*pi*a/N
    FOR i = 1 to N
        LET x(i) = a*cos(angle)
        LET y(i) = 0
        LET z(i) = a*sin(angle)
        ! direction of current segment at (x,y,z)
        LET dLx(i) = -dL*sin(angle)
        LET dLy(i) = 0
        LET dLz(i) = dL*cos(angle)
        LET angle = angle + delta_angle
    NEXT i
END SUB

SUB screen(a)
    ! unit of distance is 1/10 screen width
    CALL compute_aspect_ratio(3*a,xwin,ywin)
    SET WINDOW -xwin,xwin,-ywin,ywin
    BOX CIRCLE -a,a,-0.4,0.4
    BOX AREA -0.1,0.1,-0.1,0.1
    PRINT "Click with mouse to start line"
    PRINT "Hit any key to stop line"
    PRINT "Click on filled box and hit any key to end"
END SUB

SUB draw(N,delta,x(),y(),z(),dLx(),dLy(),dLz(),xcursor,ycursor)
    LET rx = xcursor
    LET ry = ycursor
    DO
       LET Bx = 0.0
       LET By = 0.0
       LET Bz = 0.0
       FOR i = 1 to N
           ! compute displacement from position to current segment
           LET dx = rx - x(i)
           LET dy = ry - y(i)
           LET dz = zcursor - z(i)
           LET r = sqr(dx*dx + dy*dy + dz*dz)
           LET B0 = 1/(r*r*r)     ! assume current equals 1 
           ! B proportional to dL x r
           LET Bx = Bx + B0*(dLy(i)*dz - dLz(i)*dy)
           LET By = By + B0*(dLz(i)*dx - dLx(i)*dz)
           LET Bz = Bz + B0*(dLx(i)*dy - dLy(i)*dx)
       NEXT i
       LET B = sqr(Bx*Bx + By*By + Bz*Bz)
       LET rx = rx + delta*Bx/B   ! new position on flux line
       LET ry = ry + delta*By/B
       LET rz = rz + delta*Bz/B
       PLOT rx,ry;
    LOOP until key input
    GET KEY kk                    ! reset key input
    PLOT
END SUB