##############################################################################
# Simulation of Plane-Polarized and Circular-Polarized Electromagnetic Waves
#
# Reinhard A. Schumacher (All rights reserved)
# Carnegie Mellon University
# May 2007
#
# Use MKS units for everything, except we use centimeters for length units
#
# Modifications:
#
from visual import *
from random import random, randrange
from visual.graph import *
print """
##############################################################################
Simulation of  Traveling Electromagnetic Wave: plane and circular polarization

 R. A. Schumacher
 Carnegie Mellon University
 All rights reserved.
------------------------------------------------------------------------------

 Two plane-polarized electromagnetic waves propagate along the same axis
 though space.  The red arrows represent the electric field strength along
 the white line.  The green arrows represent the magnetic field. To
 control the display:

 Click "Wave 1" ball to toggle the visibility of the first wave
 Click "Wave 2" ball to toggle the visibility of the second wave
 Click "Wave Total" ball to toggle the visibility of the sum of Waves 1 & 2
 Click anywhere else to pause or resume the amination.

 To start the electron moving, toggle one or the other of the waves off and on.
 The electron moves in response to the electic and magnetic field strengths it
 encounters along its way.

Version 1.5
###########################################################################
"""
twopi     = 2.* pi
dera      = pi/180.
###########################################################################
# USER ADJUSTABLE INPUTS GO HERE:
wlambda   = 100               #wavelength in centimeters
khat      = vector(1,0.5,0)   #propagation direction, not necessarily a unit vector
khat      = vector(1,0.0,0)   #propagation direction, not necessarily a unit vector
origin    = vector (0.,0.,0.) #opportunity to shift orgin of computation
clight    = 3E10              #speed of light in cm/sec
dt        = 2e-11             #time step in seconds
dr        = 3.                #distance step in centimeters
efieldmag = 95.               #scale factor for electric field in "centimeters"
escale    = 3.e8              #scale for electric field in Volts/centimeter
bscale    = 3.e-11            #scale for magnetic field in who- knows-what units
phi       = 90*dera           #phase angle beteen waves

speed     = 20                #maximum update speed per second: depends on CPU power
###########################################################################
thack     = vector (0,0,1)    #hack used to contruct a vector perpendicular to khat
#
scene1 = display(
    title = "Electromagnetic Waves" ,
    width=700, height=700,
    autoscale=1, background=(.9,.9,.9), #background=(.3,.4,.5), #background=color.black,  #
    x=0, y=0, forward=(-0.25,-0.25,-1))
def wirecube (s,x,y,z):
    c=curve (color=color.white, radius=0.5)
    pts = [(-s+x, -s+y, -s+z),(-s+x, -s+y, s+z), (-s+x, s+y, s+z),
           (-s+x, s+y, -s+z), (-s+x, -s+y, -s+z), (s+x, -s+y, -s+z),
           (s+x, s+y, -s+z), (-s+x, s+y, -s+z), (s+x, s+y, -s+z),
           (s+x, s+y, s+z), (-s+x, s+y, s+z), (s+x, s+y, s+z),
           (s+x, -s+y, s+z), (-s+x, -s+y, s+z), (s+x, -s+y, s+z),(s+x, -s+y, -s+z)]
    for pt in pts:
        c.append(pos=pt)
#
# Set up the visual scene
#
scene1.select()  #makes first scene the current default for new objects
wirecube(50.,0,0,0)  #2x50=100 cm on a side
base = box(pos=vector(0,-50,0),height=1.0, length=100.0, width=100.0,
           color=(0.7,0.7,0.9))  #color=(0.0,0.5,1.0))
ishow1 = 0
s_show1=sphere(pos=(-50,50,50), radius=2., color=color.green)
s_show1_label=label(pos=s_show1.pos, text="Wave 1", yoffset=+10,
                    height=12,opacity=0, box=0, line=0,color=color.blue)
ishow2 = 0
s_show2=sphere(pos=(50,50,50), radius=2., color=color.red)
s_show2_label=label(pos=s_show2.pos, text="Wave 2", yoffset=+10,
                    height=12,opacity=0, box=0, line=0,color=color.blue)
ishowt = 0
s_showt=sphere(pos=(00,50,50), radius=2., color=color.red)
s_showt_label=label(pos=s_showt.pos, text="Wave Total", yoffset=+10,
                    height=12,opacity=0, box=0, line=0,color=color.blue)
ishowe = 0
s_showe=sphere(pos=(-50,-50,50), radius=2., color=color.red)
s_showe_label=label(pos=s_showe.pos, text="Reset Electron", yoffset=-10,
                    height=12,opacity=0, box=0, line=0,color=color.blue)
electron=sphere(pos=(0,0,0), radius=2., color=color.blue)
electron.charge = -1.6e-19 # Coulombs
electron.mass = 9.1e-31   # kilograms
electron.vel  = vector(0,0,0)      # velocity in cm/sec
electron.mom  = electron.vel*electron.mass
etrail        = curve(display=scene1,color=color.blue)  #follows the electron's path
electron.visible = 1
etrail.visble    = 1

#
khat     = norm(khat)
kvector  = khat*twopi/wlambda
omega    = mag(kvector)*clight
epsilon1 = norm(cross(khat,thack))
epsilon2 = cross(khat,epsilon1)     # epsilon1 
#epsilon2 = norm(epsilon1+epsilon2) # hack to put Wave2 at 45 degrees to Wave1
rstart = origin - 60*khat
rend   = origin + 60*khat
axis = curve(color=color.white,radius = 0.5)
zero = vector(0.,0.,0.)
pts = [rstart,rend]
for pt in pts:
    axis.append(pos=pt)

earrows1   = []
earrows2   = []
earrowstot = []
barrows1   = []
barrows2   = []
barrowstot = []
for i in range(-20,20):
    position   = origin+i*khat*dr
    evector1   = epsilon1*efieldmag*cos(dot(kvector,position))
    bvector1   = cross(khat,evector1)
    eposition1 = position - 0.5*evector1
    bposition1 = position - 0.5*bvector1
    earrows1 = earrows1+[arrow(pos=eposition1,axis=evector1,
                             shaftwidth=1.,color=color.red,visible=0)]
    barrows1 = barrows1+[arrow(pos=bposition1,axis=bvector1,
                             shaftwidth=1.,color=color.green,visible=0)]
    #
    evector2   = epsilon2*efieldmag*cos(dot(kvector,position)+phi)
    bvector2   = cross(khat,evector2)
    eposition2 = position - 0.5*evector2
    bposition2 = position - 0.5*bvector2
    earrows2 = earrows2+[arrow(pos=eposition2,axis=evector2,
                             shaftwidth=1.,color=color.red,visible=0)]
    barrows2 = barrows2+[arrow(pos=bposition2,axis=bvector2,
                             shaftwidth=1.,color=color.green,visible=0)]
    #
    evectortot   = evector1 + evector2
    bvectortot   = bvector2 + bvector2
    epositiontot = position - 0.5*evectortot
    bpositiontot = position - 0.5*bvectortot
    earrowstot = earrowstot + [arrow(pos=epositiontot,axis=evectortot,
                             shaftwidth=1.,color=color.red,visible=0)]
    barrowstot = barrowstot + [arrow(pos=bpositiontot,axis=bvectortot,
                             shaftwidth=1.,color=color.green,visible=0)]

tvalue = 0
#m = scene1.mouse.getevent()
scene1.autoscale = 0
for i in range(-20,20):
    earrows1[i].visible = 1
    barrows1[i].visible = 1
    earrows2[i].visible = 0
    barrows2[i].visible = 0
    earrowstot[i].visible = 0
    barrowstot[i].visible = 0
#
# Top of main loop
#
while 1:
    rate(speed)
    tvalue = tvalue+ dt
    for i in range(-20,20):
        position   = origin+i*khat*dr
        kayare     = dot(kvector,position)
        omegat     = omega*tvalue
        evector1   = epsilon1*efieldmag*cos(kayare-omegat)
        bvector1   = cross(khat,evector1)
        eposition1 = position - 0.5*evector1
        bposition1 = position - 0.5*bvector1
        earrows1[i].pos  = eposition1
        earrows1[i].axis = evector1
        barrows1[i].pos  = bposition1
        barrows1[i].axis = bvector1
        #
        evector2   = epsilon2*efieldmag*cos(kayare-omegat+phi)
        bvector2   = cross(khat,evector2)
        eposition2 = position - 0.5*evector2
        bposition2 = position - 0.5*bvector2
        earrows2[i].pos  = eposition2
        earrows2[i].axis = evector2
        barrows2[i].pos  = bposition2
        barrows2[i].axis = bvector2
        #     
        evectortot   = evector1 + evector2
        bvectortot   = bvector1 + bvector2
        epositiontot = position - 0.5*evectortot
        bpositiontot = position - 0.5*bvectortot
        earrowstot[i].pos  = epositiontot
        earrowstot[i].axis = evectortot
        barrowstot[i].pos  = bpositiontot
        barrowstot[i].axis = bvectortot
#
    position   = electron.pos
    kayare     = dot(kvector,position)
    omegat     = omega*tvalue
    evector1   = epsilon1*escale*cos(kayare-omegat)
    bvector1   = bscale*cross(khat,evector1)
#
    evector2   = epsilon2*escale*cos(kayare-omegat+phi)
    bvector2   = bscale*cross(khat,evector2)
#
    eforce = zero
    mforce = zero
    if ishow1==1 or ishowt==1:
        eforce = eforce+ electron.charge*evector1
        mforce = mforce+ electron.charge*cross(electron.vel,bvector1)
    if ishow2==1 or ishowt==1:
        eforce = eforce+ electron.charge*evector2
        mforce = mforce+ electron.charge*cross(electron.vel,bvector2)
    force  = eforce + mforce # force on electron due to E and B fields
    electron.mom = electron.mom + force*dt
    electron.vel = electron.mom/electron.mass #non-relativistic
    electron.pos = electron.pos + electron.vel*dt
    if ishowe == 0:
        etrail.append(pos=electron.pos)
    ishowe = 0
#    print "Forces",eforce, mforce
#    print "velocity,position",electron.vel,electron.pos
#    m = scene1.mouse.getevent()

#
#   Mousing around...
#
    if scene1.mouse.events:
        m = scene1.mouse.getevent()
#
#       Click on a labeled sphere to get the stated action.
#
        if m.click:
            if m.pick is s_show1:
                if ishow1 == 0:
                    print "Show first plane wave..."
                    ishow1 = 1
                    s_show1.color = color.green
                    for i in range(-20,20):
                        earrows1[i].visible = 1
                        barrows1[i].visible = 1
                else:
                    print "Hide first plane wave..."
                    ishow1 = 0
                    s_show1.color = color.red
                    for i in range(-20,20):
                        earrows1[i].visible = 0
                        barrows1[i].visible = 0
            if m.pick is s_show2:
                if ishow2 == 0:
                    print "Show second plane wave..."
                    ishow2 = 1
                    s_show2.color = color.green
                    for i in range(-20,20):
                        earrows2[i].visible = 1
                        barrows2[i].visible = 1
                else:
                    print "Hide second plane wave..."
                    ishow2 = 0
                    s_show2.color = color.red
                    for i in range(-20,20):
                        earrows2[i].visible = 0
                        barrows2[i].visible = 0
            if m.pick is s_showt:
                if ishowt == 0:
                    print "Show total wave..."
                    ishowt = 1
                    s_showt.color = color.green
                    for i in range(-20,20):
                        earrowstot[i].visible = 1
                        barrowstot[i].visible = 1
                else:
                    print "Hide total wave..."
                    s_showt.color = color.red
                    ishowt = 0
                    for i in range(-20,20):
                        earrowstot[i].visible = 0
                        barrowstot[i].visible = 0
            if m.pick is s_showe:
                    print "Move electron back to origin..."
                    ishowe = 1
                    #s_showe.color = color.green
                    electron.pos = zero
                    electron.mom = zero
                    etrail = curve(display=scene1,color=color.blue)  #reset curve

            else:
                print "Pausing  animation"
                scene1.mouse.getclick() #inhibits animation until the mouse is clicked
                print "Resuming animation"


print "All done."
#  END OF CODE