PROGRAM fieldline
! draw electric field lines in two dimensions
LIBRARY "csgraphics"
DIM x(10),y(10),q(10)
CALL screen(a,pos$,neg$)
CALL charges(N,x(),y(),q(),qtotal,a,pos$,neg$)
! draw field lines
CALL draw_lines(N,x(),y(),q(),qtotal,a)
END

SUB screen(a,pos$,neg$)
    LET L = 10
    CALL compute_aspect_ratio(L,xwin,ywin)
    SET WINDOW -xwin,xwin,-ywin,ywin
    LET a = 0.2              ! "radius" of visual image of charges
    SET COLOR "blue"
    BOX CIRCLE -a,a,-a,a
    FLOOD 0,0
    BOX KEEP -a,a,-a,a in pos$
    CLEAR
    SET COLOR "red"
    BOX CIRCLE -a,a,-a,a
    FLOOD 0,0
    BOX KEEP -a,a,-a,a in neg$
    CLEAR
END SUB

SUB charges(N,x(),y(),q(),qtotal,a,pos$,neg$)
    ! input charge values and location
    LET N = 0                ! # of point charges
    SET COLOR "black"
    DO
       CALL draw_charges(N,x(),y(),q(),a,pos$,neg$)
       INPUT prompt "charge (0 to exit input mode) = ": charge
       IF charge <> 0 then
          PRINT "place mouse at charge location and click."
          LET N = N + 1
          GET POINT x(N),y(N)     ! location of charge
          LET q(N) = charge
          LET qtotal = qtotal + charge
       END IF
       CLEAR
    LOOP until charge = 0
    ! redraw charges
    CALL draw_charges(N,x(),y(),q(),a,pos$,neg$)
END SUB

SUB draw_charges(N,x(),y(),q(),a,pos$,neg$)
    FOR i = 1 to N
        IF q(i) > 0 then
           BOX SHOW pos$ at x(i)-a,y(i)-a
        ELSE
           BOX SHOW neg$ at x(i)-a,y(i)-a
        END IF
    NEXT i
END SUB

SUB draw_lines(N,x(),y(),q(),qtotal,a)
    ! number of lines per unit charge leaving positive charges
    SET COLOR "black"
    LET lpc = 8              ! lines per charge
    LET Emin = 0.01          ! if E < Emin stop drawing fieldline
    LET sign = sgn(qtotal)
    IF sign = 0 then LET sign = 1
    LET ds_small = 0.01*sign
    LET ds_big = sign
    FOR i = 1 to N           ! loop over all charges
        IF q(i)*sign > 0 then     ! start fields at positive charges
           LET dtheta = 2*pi/(lpc*abs(q(i)))
           FOR theta = 0 to 2*pi step dtheta
               LET stop_plot$ = "no"
               LET xline = x(i) + a*cos(theta)
               LET yline = y(i) + a*sin(theta)
               DO
                  LET Ex = 0
                  LET Ey = 0
                  FOR j = 1 to N
                      ! x-distance from point to charge j
                      LET dx = xline - x(j)
                      ! y-distance from point to charge j
                      LET dy = yline - y(j)
                      LET r = sqr(dx*dx + dy*dy)
                      IF r > 0.9*a then
                         LET E0 = q(j)/(r*r*r)
                         LET Ex = Ex + E0*dx
                         LET Ey = Ey + E0*dy
                      ELSE
                         ! field line reached another charge
                         LET stop_plot$ = "yes"
                      END IF
                  NEXT j
                  LET E = sqr(Ex*Ex + Ey*Ey)
                  IF E > Emin or r < 20 then
                     ! new position on fieldline
                     LET xline = xline + ds_small*Ex/E
                     LET yline = yline + ds_small*Ey/E
                  ELSE
                     ! new position on fieldline
                     LET xline = xline + ds_big*Ex/E
                     LET yline = yline + ds_big*Ey/E
                  END IF
                  PLOT xline,yline;
                  IF key input then
                     ! user can stop drawing field line
                     GET KEY key
                     LET stop_plot$ = "yes"
                  END IF
               LOOP until stop_plot$ = "yes"
               PLOT          ! turn beam off
           NEXT theta
        END IF
    NEXT i
END SUB