      PROGRAM magnetism
      IMPLICIT NONE
      REAL*8 a, dLx(50), dLy(50), dLz(50), delta, x(50), y(50),
     +       z(50)
      REAL xcursor, ycursor
      INTEGER IR_, N
      LOGICAL Ltmp1_
      CALL GWopen(IR_, 0)
      CALL wire(N,a,delta,x,y,z,dLx,dLy,dLz)
      CALL screen(a)
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         CALL GWcappnt(IR_, xcursor, ycursor, '')
         CALL draw(N,delta,x,y,z,dLx,dLy,dLz,xcursor,ycursor)
         Ltmp1_ = (abs(xcursor) .GE. .1) .OR. (abs(ycursor) .GE. .1)
      END DO
      CALL GWquit(IR_)
      END
C
      SUBROUTINE wire(N,a,delta,x,y,z,dLx,dLy,dLz)
*     read in current segments
*     assume current equals unity
      IMPLICIT NONE
      REAL*8 angle, dL, delta_angle, pi, a, delta, x(50), y(50),
     +       z(50), dLx(50), dLy(50), dLz(50)
      PARAMETER(pi = 3.141592653589793D0)
      INTEGER N, i
      WRITE(*,'(A,$)') 'radius of loop = '
      READ(*,*) a
      WRITE(*,'(A,$)') 'number of segments = '
      READ(*,*) N
*     magnitude of magnetic field line segment
      WRITE(*,'(A,$)') 'delta = '
      READ(*,*) delta
      delta_angle = 2*pi/N
      angle = 0
      dL = 2*pi*a/N
      DO i = 1, N
         x(i) = a*cos(angle)
         y(i) = 0
         z(i) = a*sin(angle)
*        direction of current segment at (x,y,z)
         dLx(i) = -dL*sin(angle)
         dLy(i) = 0
         dLz(i) = dL*cos(angle)
         angle = angle + delta_angle
      END DO
      END
C
      SUBROUTINE screen(a)
*     unit of distance is 1/10 screen width
      IMPLICIT NONE
      REAL*8 a
      REAL xwin, ywin
      INTEGER IR_
      CALL compute_aspect_ratio(REAL(3*a),xwin,ywin)
      CALL GWindow(IR_, -xwin, -ywin, xwin, ywin)
      CALL GWellipse(IR_, REAL(-a), -0.4, REAL(a), 0.4)
      CALL GWsetmrk(IR_, 1, 0.2, 16, -1, -1)
      CALL GWputmrk(IR_, 0.0, 0.0)
      WRITE(*,*) 'Click with mouse to start line'
      WRITE(*,*) 'Hit any key to stop line'
      WRITE(*,*) 'Click on filled box and hit any key to end'
      END
C
      SUBROUTINE draw(N,delta,x,y,z,dLx,dLy,dLz,xcursor,ycursor)
      IMPLICIT NONE
      REAL*8 B, B0, Bx, By, Bz, dx, dy, dz, r, rx, ry, rz, zcursor,
     +       delta, x(50), y(50), z(50), dLx(50), dLy(50), dLz(50)
      REAL xcursor, ycursor
      INTEGER IR_, N, i, kk
      LOGICAL TBkeyinput
      DATA rz, zcursor/2*0.0D0/
      rx = xcursor
      ry = ycursor
      CALL GWmove2(IR_, REAL(rx), REAL(ry));
      DO WHILE(.NOT.TBkeyinput())
         Bx = 0.0D0
         By = 0.0D0
         Bz = 0.0D0
         DO i = 1, N
*           compute displacement from position to current segment
            dx = rx - x(i)
            dy = ry - y(i)
            dz = zcursor - z(i)
            r = sqrt(dx*dx + dy*dy + dz*dz)
            B0 = 1/(r*r*r)                  ! assume current equals 1 
*           B proportional to dL x r
            Bx = Bx + B0*(dLy(i)*dz - dLz(i)*dy)
            By = By + B0*(dLz(i)*dx - dLx(i)*dz)
            Bz = Bz + B0*(dLx(i)*dy - dLy(i)*dx)
         END DO
         B = sqrt(Bx*Bx + By*By + Bz*Bz)
         rx = rx + delta*Bx/B               ! new position on flux line
         ry = ry + delta*By/B
         rz = rz + delta*Bz/B
         CALL GWline2(IR_, REAL(rx), REAL(ry))
      END DO
      CALL TBgetkey(kk)
      END
