      PROGRAM fieldline
*     2 次元の電気力線を描く
      IMPLICIT NONE
      REAL*8 a, q(10), qtotal
      REAL x(10), y(10)
      INTEGER IR_, N
      DATA qtotal/0.0D0/
      CALL GWopen(IR_, 0)
      CALL screen(a)
      CALL charges(N,x,y,q,qtotal,a)
*     電気力線の描画
      CALL draw_lines(N,x,y,q,qtotal,a)
      CALL GWquit(IR_)
      END
C
      SUBROUTINE screen(a)
      IMPLICIT NONE
      REAL*8 a
      REAL xwin, ywin
      INTEGER IR_, L
      L = 10
      CALL compute_aspect_ratio(REAL(L),xwin,ywin)
      CALL GWindow(IR_, -xwin, -ywin, xwin, ywin)
      a = 0.2D0                             ! 画面上の電荷の"半径"
      CALL GWsetmrk(IR_, 6, REAL(a*2), -1, -1, -1)
      END
C
      SUBROUTINE charges(N,x,y,q,qtotal,a)
*     電荷の大きさと位置を入力
      IMPLICIT NONE
      CHARACTER*80 Stmp1_
      REAL*8 charge, q(10), qtotal, a
      REAL x(10), y(10)
      INTEGER IR_, N
      LOGICAL Ltmp1_
      N = 0                                 ! 点電荷の数
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         CALL draw_charges(N,x,y,q,a)
         CALL GWinput(IR_, 'charge (0 to exit input mode) = ', Stmp1_)
         IF(IR_ .EQ. 0) THEN
            charge = 0
         ELSE
            READ(Stmp1_,*) charge
         ENDIF
         IF(charge .NE. 0.0D0) THEN
            CALL GWsetmsg(IR_, 
     +           'place mouse at charge location and click.')
            N = N + 1
            CALL GWcappnt(IR_, x(N), y(N), '')  ! 電荷の位置
            q(N) = charge
            qtotal = qtotal + charge
         END IF
         Ltmp1_ = charge .NE. 0.0D0
      END DO
*     再度電荷を描画
      CALL draw_charges(N,x,y,q,a)
      END
C
      SUBROUTINE draw_charges(N,x,y,q,a)
      IMPLICIT NONE
      REAL*8 q(10), a
      REAL x(10), y(10)
      INTEGER IR_, i, N
      DO i = 1, N
         IF(q(i) .GT. 0) THEN
            CALL GWsetmrk(IR_, -1, -1.0, 16, -1, -1)
         ELSE
            CALL GWsetmrk(IR_, -1, -1.0, 13, -1, -1)
         END IF
         CALL GWputmrk(IR_, x(i), y(i))
      END DO
      END
C
      SUBROUTINE draw_lines(N,x,y,q,qtotal,a)
*     正電荷から出る単位電荷当たりの電気力線の本数
      IMPLICIT NONE
      REAL*8 E, E0, Emin, Ex, Ey, ds_big, ds_small, dtheta, dx, dy,
     +       pi, r, sgn, theta, xline, yline, q(10), qtotal, a
      PARAMETER(pi = 3.141592653589793D0)
      REAL x(10), y(10)
      INTEGER IR_, i, j, key, lpc, N
      CHARACTER stop_plot*80
      LOGICAL Ltmp1_, TBkeyinput
      DATA r/0.0D0/
      lpc = 8                               ! 単位電荷当たりの電気力線の本数
      Emin = 0.01D0                         ! E < Emin なら描画を停止する
      sgn = sign(1.0D0,qtotal)
      IF(sgn .EQ. 0) sgn = 1
      ds_small = 0.01D0*sgn
      ds_big = sgn
      DO i = 1, N
         IF(q(i)*sgn .GT. 0) THEN
            dtheta = 2*pi/(lpc*abs(q(i)))
            DO theta = 0, 2*pi, dtheta
               CALL GWmove2(IR_, x(i), y(i))
               stop_plot = 'no'
               xline = x(i) + a*cos(theta)
               yline = y(i) + a*sin(theta)
               Ltmp1_ = .TRUE.
               DO WHILE(Ltmp1_)
                  Ex = 0
                  Ey = 0
                  DO j = 1, N
*                    点から電荷jまでのx方向の距離
                     dx = xline - x(j)
*                    点から電荷jまでのy方向の距離
                     dy = yline - y(j)
                     r = sqrt(dx*dx + dy*dy)
                     IF(r .GT. 0.9D0*a) THEN
                        E0 = q(j)/r**3
                        Ex = Ex + E0*dx
                        Ey = Ey + E0*dy
                     ELSE
*                       電気力線が他の電荷の位置に到達
                        stop_plot = 'yes'
                     END IF
                  END DO
                  E = sqrt(Ex*Ex + Ey*Ey)
                  IF(E .GT. Emin .OR. r .LT. 20) THEN
*                    電気力線上の新しい位置
                     xline = xline + ds_small*Ex/E
                     yline = yline + ds_small*Ey/E
                  ELSE
*                    電気力線上の新しい位置
                     xline = xline + ds_big*Ex/E
                     yline = yline + ds_big*Ey/E
                  END IF
                  CALL GWline2(IR_, REAL(xline), REAL(yline))
                  IF(TBkeyinput()) THEN
*                    電気力線の描画を停止することができる
                     CALL TBgetkey(key)
                     stop_plot = 'yes'
                  END IF
                  Ltmp1_ = stop_plot .NE. 'yes'
               END DO
            END DO
         END IF
      END DO
      END
