PROGRAM radiation
! compute electric fields due to accelerating charge and 
! plot electric field lines
LIBRARY "csgraphics"
DIM snapshot$(100)
CALL initial(L,nsnap,radius,lpc,dt,dtheta)
FOR isnap = 1 to nsnap
    CALL draw_field(L,isnap,radius,dt,dtheta,snapshot$())
NEXT isnap
CALL animate(L,nsnap,snapshot$())
END

SUB initial(L,nsnap,a,lpc,dt,dtheta)
    LET L = 20
    CALL compute_aspect_ratio(2*L,xwin,ywin)
    SET WINDOW -xwin,xwin,-ywin,ywin
    LET a = 0.3                   ! "radius" of visual image of charge
    LET lpc = 10                  ! number of field lines 
    LET dt = 0.4                  ! time between plots
    LET dtheta = 2*pi/lpc
    LET nsnap = 15                ! number of snapshots
END SUB

SUB draw_field(L,isnap,a,dt,dtheta,snapshot$())
    ! adopted from Program fieldline
    DIM E(3),R(3),r_ret(3),v_ret(3),a_ret(3)
    DECLARE DEF dotproduct
    LET ds = 1
    LET t = isnap*dt              ! observation time
    LET R(1) = 0
    LET R(2) = 0
    CALL motion(R(),t,r_ret(),v_ret(),a_ret())
    ! find charge position at time t
    LET x = R(1) - r_ret(1)
    LET y = R(2) - r_ret(2)
    BOX CIRCLE x-a,x+a,y-a,y+a
    FOR theta = 0 to 2*pi step dtheta
        LET xline = x + a*cos(theta)
        LET yline = y + a*sin(theta)
        PLOT xline,yline;
        DO
           LET R(1) = xline
           LET R(2) = yline
           CALL field(R(),t,E())  ! compute fields 
           LET Emag = sqr(dotproduct(E(),E()))
           ! new position on field line
           LET xline = xline + ds*E(1)/Emag
           LET yline = yline + ds*E(2)/Emag
           PLOT xline,yline;
        LOOP until abs(xline) > L or abs(yline) > L
        PLOT                      ! turn beam off
    NEXT theta
    BOX LINES -L,L,-L,L
    BOX KEEP -L,L,-L,L in snapshot$(isnap)
    CLEAR
END SUB

SUB motion(R(),t_ret,r_ret(),v_ret(),a_ret())
    ! compute motion of source
    LET r_ret(1) = R(1) - 0.2*cos(t_ret)    ! harmonic motion 
    LET r_ret(2) = R(2)
    LET r_ret(3) = R(3)
    LET v_ret(1) = -0.2*sin(t_ret)     ! particle velocity 
    LET v_ret(2) = 0
    LET v_ret(3) = 0
    LET a_ret(1) = -0.2*cos(t_ret)     ! particle acceleration 
    LET a_ret(2) = 0
    LET a_ret(3) = 0
END SUB

SUB field(R(),t,E())
    ! compute electric field vector
    DECLARE DEF dotproduct
    DIM r_ret(3),v_ret(3),a_ret(3),u_ret(3),uxa(3),w_ret(3)
    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())
    LET v2 = dotproduct(v_ret(),v_ret())
    LET dist_ret = sqr(dotproduct(r_ret(),r_ret()))
    FOR i = 1 to 3
        LET u_ret(i) = r_ret(i)/dist_ret - v_ret(i)
    NEXT i
    CALL crossproduct(u_ret(),a_ret(),uxa())
    CALL crossproduct(r_ret(),uxa(),w_ret())
    LET ru = dotproduct(r_ret(),u_ret())
    LET E0 = dist_ret/ru^3
    FOR i = 1 to 3
        LET E(i) = E0*(u_ret(i)*(1 - v2) + w_ret(i))
    NEXT i
END SUB

SUB retarded_time(R(),t,t_ret)
    LET tb = t                    ! upper guess for retarded time
    CALL fanddf(fb,dfdt,R(),t,tb)
    LET ta = -1/1.6               ! lower guess for retarded time
    ! insure that f(ta) > 0 and f(tb) < 0
    DO
       LET ta = ta*1.6
       CALL fanddf(fa,dfdt,R(),t,ta)
    LOOP until fa > 0
    LET t_ret = 0.5*(ta + tb)
    CALL fanddf(f,dfdt,R(),t,t_ret)
    CALL zero_of_f(f,dfdt,R(),t,t_ret,ta,tb)
END SUB

SUB fanddf(f,dfdt,R(),t,t_ret)
    ! calculate f and df/dt_r
    DIM r_ret(3),v_ret(3),a_ret(3)
    DECLARE DEF dotproduct
    CALL motion(R(),t_ret,r_ret(),v_ret(),a_ret())
    LET dist_ret = sqr(dotproduct(r_ret(),r_ret()))
    LET f = t - t_ret - dist_ret
    ! derivative evaluated at retarded time
    LET dfdt = -1 + dotproduct(r_ret(),v_ret())/dist_ret
END SUB

SUB 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
    LET eps = 1e-6
    LET dt_r = tb - ta
    FOR j = 1 to 100
        ! Newton difference between successive t_ret estimates
        LET dt = f/dfdt
        LET sign = ((t_ret - ta) - dt)*((t_ret - tb) - dt)
        LET dt_old = dt_r
        IF sign >= 0 or abs(2*dt) > 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
           LET dt_r = 0.5*(ta - tb)
           LET t_ret = tb + dt_r
        ELSE                      ! use Newton's method
           LET dt_r = dt
           LET t_ret = t_ret - dt_r
        END IF
        IF abs(dt_r) < eps then EXIT SUB    ! convergence test
        CALL fanddf(f,dfdt,R(),t,t_ret)
        IF f < 0 then
           LET tb = t_ret
        ELSE
           LET ta = t_ret
        END IF
    NEXT j
    PRINT "too many iterations"
    STOP
END SUB

SUB animate(L,nsnap,snapshot$())
    FOR isnap = 1 to nsnap        ! show motion picture
        BOX SHOW snapshot$(isnap) at -L,-L
        PAUSE 1
    NEXT isnap
END SUB

SUB crossproduct(a(),b(),c())
    LET c(1) = a(2)*b(3) - a(3)*b(2)
    LET c(2) = a(3)*b(1) - a(1)*b(3)
    LET c(3) = a(1)*b(2) - a(2)*b(1)
END SUB

DEF dotproduct(a(),b())
    LET dotproduct = a(1)*b(1) + a(2)*b(2) + a(3)*b(3)
END DEF