      PROGRAM radiation
*     compute electric fields due to accelerating charge and 
*     plot electric field lines
      IMPLICIT NONE
      REAL*8 L, dt, dtheta, radius
      INTEGER IR_, isnap, lpc, nsnap
      CALL GWopen(IR_, 0)
      CALL initial(L,nsnap,radius,lpc,dt,dtheta)
      DO isnap = 1, nsnap
         CALL draw_field(L,isnap,radius,dt,dtheta)
      END DO
      CALL animate(L,nsnap)
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(L,nsnap,a,lpc,dt,dtheta)
      IMPLICIT NONE
      REAL*8 pi, L, a, dt, dtheta
      REAL xwin, ywin
      PARAMETER(pi = 3.141592653589793D0)
      INTEGER IR_, nsnap, lpc
      L = 20
      CALL compute_aspect_ratio(REAL(2*L),xwin,ywin)
      CALL GWindow(IR_, -xwin, -ywin, xwin, ywin)
      a = 0.3D0                             ! "radius" of visual image of charge
      lpc = 10                              ! number of field lines 
      dt = 0.4D0                            ! time between plots
      dtheta = 2*pi/lpc
      nsnap = 15                            ! number of snapshots
      END
C
      SUBROUTINE draw_field(L,isnap,a,dt,dtheta)
*     adopted from Program fieldline
      IMPLICIT NONE
      REAL*8 E(3), Emag, R(3), a_ret(3), ds, pi, r_ret(3), t, theta,
     +       v_ret(3), x, xline, y, yline, L, a, dt, dtheta, dotproduct
      PARAMETER(pi = 3.141592653589793D0)
      INTEGER IR_, isnap
      LOGICAL Ltmp1_
      DATA E, R, a_ret, r_ret, v_ret/15*0.0D0/
      CALL GWsetogn(IR_, 0, -isnap)
      CALL GWsetpen(IR_, -1, -1, -1, -1)
      ds = 1
      t = isnap*dt                          ! observation time
      R(1) = 0
      R(2) = 0
      CALL motion(R,t,r_ret,v_ret,a_ret)
*     find charge position at time t
      x = R(1) - r_ret(1)
      y = R(2) - r_ret(2)
      CALL GWellipse(IR_, REAL(x-a), REAL(y-a), REAL(x+a), REAL(y+a))
      DO theta = 0, 2*pi, dtheta
         xline = x + a*cos(theta)
         yline = y + a*sin(theta)
         CALL GWline2(IR_, REAL(xline), REAL(yline))
         Ltmp1_ = .TRUE.
         DO WHILE(Ltmp1_)
            R(1) = xline
            R(2) = yline
            CALL field(R,t,E)               ! compute fields 
            Emag = sqrt(dotproduct(E,E))
*           new position on field line
            xline = xline + ds*E(1)/Emag
            yline = yline + ds*E(2)/Emag
            CALL GWline2(IR_, REAL(xline), REAL(yline))
            Ltmp1_ = abs(xline) .LE. L .AND. abs(yline) .LE. L
         END DO
*        turn beam off
      END DO
      CALL GWrect(IR_, REAL(-L), REAL(-L), REAL(L), REAL(L))
      END
C
      SUBROUTINE motion(R,t_ret,r_ret,v_ret,a_ret)
*     compute motion of source
      IMPLICIT NONE
      REAL*8 R(3), t_ret, r_ret(3), v_ret(3), a_ret(3)
      r_ret(1) = R(1) - 0.2D0*cos(t_ret)    ! harmonic motion 
      r_ret(2) = R(2)
      r_ret(3) = R(3)
      v_ret(1) = -0.2D0*sin(t_ret)          ! particle velocity 
      v_ret(2) = 0
      v_ret(3) = 0
      a_ret(1) = -0.2D0*cos(t_ret)          ! particle acceleration 
      a_ret(2) = 0
      a_ret(3) = 0
      END
C
      SUBROUTINE field(R,t,E)
*     compute electric field vector
      IMPLICIT NONE
      REAL*8 E0, a_ret(3), dist_ret, r_ret(3), ru, t_ret, u_ret(3),
     +       uxa(3), v2, v_ret(3), w_ret(3), R(3), t, E(3), dotproduct
      INTEGER i
      CALL retarded_time(R,t,t_ret)         ! find retarded time
*     dynamical variables of moving charge
      CALL motion(R,t_ret,r_ret,v_ret,a_ret)
      v2 = dotproduct(v_ret,v_ret)
      dist_ret = sqrt(dotproduct(r_ret,r_ret))
      DO i = 1, 3
         u_ret(i) = r_ret(i)/dist_ret - v_ret(i)
      END DO
      CALL crossproduct(u_ret,a_ret,uxa)
      CALL crossproduct(r_ret,uxa,w_ret)
      ru = dotproduct(r_ret,u_ret)
      E0 = dist_ret/ru**3
      DO i = 1, 3
         E(i) = E0*(u_ret(i)*(1 - v2) + w_ret(i))
      END DO
      END
C
      SUBROUTINE retarded_time(R,t,t_ret)
      IMPLICIT NONE
      REAL*8 dfdt, f, fa, fb, ta, tb, R(3), t, t_ret
      LOGICAL Ltmp1_
      tb = t                                ! upper guess for retarded time
      CALL fanddf(fb,dfdt,R,t,tb)
      ta = -1/1.6D0                         ! lower guess for retarded time
*     insure that f(ta) > 0 and f(tb) < 0
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         ta = ta*1.6D0
         CALL fanddf(fa,dfdt,R,t,ta)
         Ltmp1_ = fa .LE. 0.0D0
      END DO
      t_ret = 0.5D0*(ta + tb)
      CALL fanddf(f,dfdt,R,t,t_ret)
      CALL zero_of_f(f,dfdt,R,t,t_ret,ta,tb)
      END
C
      SUBROUTINE fanddf(f,dfdt,R,t,t_ret)
*     calculate f and df/dt_r
      IMPLICIT NONE
      REAL*8 a_ret(3), dist_ret, r_ret(3), v_ret(3), f, dfdt, R(3), t,
     +       t_ret, dotproduct
      CALL motion(R,t_ret,r_ret,v_ret,a_ret)
      dist_ret = sqrt(dotproduct(r_ret,r_ret))
      f = t - t_ret - dist_ret
*     derivative evaluated at retarded time
      dfdt = -1 + dotproduct(r_ret,v_ret)/dist_ret
      END
C
      SUBROUTINE zero_of_f(f,dfdt,R,t,t_ret,ta,tb)
*     do no more than 100 iterations to find the value of the reduced
*     time t_ret such that f = t - t_ret - dist_ret = 0
      IMPLICIT NONE
      REAL*8 dt, dt_old, dt_r, eps, j, sgn, f, dfdt, R(3), t, t_ret,
     +       ta, tb
      INTEGER IR_
      eps = 1d-6
      dt_r = tb - ta
      DO j = 1, 100
*        Newton difference between successive t_ret estimates
         dt = f/dfdt
         sgn = ((t_ret - ta) - dt)*((t_ret - tb) - dt)
         dt_old = dt_r
         IF(sgn .GE. 0 .OR. abs(2*dt) .GT. abs(dt_old)) THEN
*           use bisection method if next Newton iteration is not 
*           between ta and tb, or if dt is not less than half
*           of the old value of dt_r
            dt_r = 0.5D0*(ta - tb)
            t_ret = tb + dt_r
         ELSE                               ! use Newton's method
            dt_r = dt
            t_ret = t_ret - dt_r
         END IF
         IF(abs(dt_r) .LT. eps) RETURN      ! convergence test
         CALL fanddf(f,dfdt,R,t,t_ret)
         IF(f .LT. 0) THEN
            tb = t_ret
         ELSE
            ta = t_ret
         END IF
      END DO
      WRITE(*,*) 'too many iterations'
      CALL GWquit(IR_)
      STOP
      END
C
      SUBROUTINE animate(L,nsnap)
      IMPLICIT NONE
      REAL*8 L
      INTEGER IR_, isnap, nsnap, lsnap, i, loop
      DATA loop/10/
      isnap = 1
      CALL GWflush(IR_, isnap)
      DO i = 2, nsnap*loop                  ! show motion picture
         CALL csleep(100)
         lsnap = isnap
         isnap = isnap + 1
         if(isnap .GT. nsnap) isnap = 1
         CALL GWsetogn(IR_, lsnap, -lsnap)
         CALL GWflush(IR_, isnap)
      END DO
      END
C
      SUBROUTINE crossproduct(a,b,c)
      IMPLICIT NONE
      REAL*8 a(3), b(3), c(3)
      c(1) = a(2)*b(3) - a(3)*b(2)
      c(2) = a(3)*b(1) - a(1)*b(3)
      c(3) = a(1)*b(2) - a(2)*b(1)
      END
C
      FUNCTION dotproduct(a,b)
      IMPLICIT NONE
      REAL*8 dotproduct, a(3), b(3)
      dotproduct = a(1)*b(1) + a(2)*b(2) + a(3)*b(3)
      END
