em hw 4

import numpy
import pylab

#Defining constants#


N = 80

R = 0.025
rhonot = 1.0*10**10
wnot = 1.0*10**-3
mu0 = 4*numpy.pi*10**-7
e0 = 8.85*10**-12
c = 1/numpy.sqrt(mu0*e0)
L = 0.1
h = L/N
Toff = 8.0*10**-11
timer = 0.
dt = (0.25*h)/c
ntot = 30
ttot = ntot*dt

#Creating grid of data points#
grid = numpy.arange(0,2*L+.001,h)
x = numpy.zeros(grid.size)
y = numpy.zeros(grid.size)
z = numpy.zeros(grid.size)
for i in range(0,grid.size):
        x[i] = -L + i*h
        y[i] = -L + i*h
        z[i] = -L + i*h

#Creating 3d vectors of all fields#
ex = numpy.zeros([grid.size, grid.size, grid.size])
ey = numpy.zeros([grid.size, grid.size, grid.size])
ez = numpy.zeros([grid.size, grid.size, grid.size])
bx = numpy.zeros([grid.size, grid.size, grid.size])
by = numpy.zeros([grid.size, grid.size, grid.size])
bz = numpy.zeros([grid.size, grid.size, grid.size])
ehx = numpy.zeros([grid.size, grid.size, grid.size])
ehy = numpy.zeros([grid.size, grid.size, grid.size])
ehz = numpy.zeros([grid.size, grid.size, grid.size])
bhx = numpy.zeros([grid.size, grid.size, grid.size])
bhy = numpy.zeros([grid.size, grid.size, grid.size])
bhz = numpy.zeros([grid.size, grid.size, grid.size])
Jx = numpy.zeros([grid.size, grid.size, grid.size])
Jy = numpy.zeros([grid.size, grid.size, grid.size])
Jz = numpy.zeros([grid.size, grid.size, grid.size])
Jhx = numpy.zeros([grid.size, grid.size, grid.size])
Jhy = numpy.zeros([grid.size, grid.size, grid.size])
Jhz = numpy.zeros([grid.size, grid.size, grid.size])
estore = numpy.zeros([grid.size, grid.size, grid.size])
bstore = numpy.zeros([grid.size, grid.size, grid.size])

#Beginning loop#
while timer < ttot:
        print timer/ttot*100
        timer = timer + dt
     #  print timer
#Calculating n+1/2 Time step#
        for i in range(1,x.size-1):
                for j in range(1,x.size-1):
                        for k in range(1,x.size-1):
                                if (numpy.sqrt(x[i]**2 + y[j]**2 + z[k]**2) <= R) and (timer <= Toff):
                                        Jhx[i,j,k] = (1./2)*(rhonot*wnot)*(1 + numpy.cos(numpy.pi*numpy.sqrt(x[i]**2 + y[j]**2 + z[k]**2)/R))*(-y[j])
                                        Jhy[i,j,k] = (1./2)*(rhonot*wnot)*(1 + numpy.cos(numpy.pi*numpy.sqrt(x[i]**2 + y[j]**2 + z[k]**2)/R))*(x[i])
                                        Jhz[i,j,k] = 0

                                else:
                                        Jhx[i,j,k] = 0.
                                        Jhy[i,j,k] = 0.
                                        Jhz[i,j,k] = 0.

                                ehx[i,j,k] = ex[i,j,k] + (dt/(4*mu0*e0*h))*(bz[i,j+1,k] - bz[i,j-1,k] - by[i,j,k+1] + by[i,j,k-1]) - (dt/(2.*e0))*Jx[i,j,k]
                                ehy[i,j,k] = ey[i,j,k] + (dt/(4*mu0*e0*h))*(bx[i,j,k+1] - bx[i,j,k-1] - bz[i+1,j,k] + bz[i-1,j,k]) - (dt/(2.*e0))*Jy[i,j,k]
                                ehz[i,j,k] = ez[i,j,k] + (dt/(4*mu0*e0*h))*(by[i+1,j,k] - by[i-1,j,k] - bx[i,j+1,k] + bx[i,j-1,k]) - (dt/(2.*e0))*Jz[i,j,k]

                                bhx[i,j,k] = bx[i,j,k] - (dt/(4*h))*(ez[i,j+1,k] - ez[i,j-1,k] - ey[i,j,k+1] + ey[i,j,k-1])
                                bhy[i,j,k] = by[i,j,k] - (dt/(4*h))*(ex[i,j,k+1] - ex[i,j,k-1] - ez[i+1,j,k] + ez[i-1,j,k])
                                bhz[i,j,k] = bz[i,j,k] - (dt/(4*h))*(ey[i+1,j,k] - ey[i-1,j,k] - ex[i,j+1,k] + ex[i,j-1,k])

#Calculating n+1 Time step#
        for i in range(1,x.size-1):
                for j in range(1,x.size-1):
                        for k in range(1,x.size-1):
                                if (numpy.sqrt(x[i]**2 + y[j]**2 + z[k]**2) <= R) and (timer <= Toff):
                                        Jx[i,j,k] = (1./2)*(rhonot*wnot)*(1 + numpy.cos(numpy.pi*numpy.sqrt(x[i]**2 + y[j]**2 + z[k]**2)/R))*(-y[j])
                                        Jy[i,j,k] = (1./2)*(rhonot*wnot)*(1 + numpy.cos(numpy.pi*numpy.sqrt(x[i]**2 + y[j]**2 + z[k]**2)/R))*(x[i])
                                        Jz[i,j,k] = 0

                                else:
                                        Jx[i,j,k] = 0.
                                        Jy[i,j,k] = 0.
                                        Jz[i,j,k] = 0.
                                ex[i,j,k] = ex[i,j,k] + (dt/(2*mu0*e0*h))*(bhz[i,j+1,k] - bhz[i,j-1,k] - bhy[i,j,k+1] + bhy[i,j,k-1]) - (dt/e0)*Jhx[i,j,k]
                                ey[i,j,k] = ey[i,j,k] + (dt/(2*mu0*e0*h))*(bhx[i,j,k+1] - bhx[i,j,k-1] - bhz[i+1,j,k] + bhz[i-1,j,k]) - (dt/e0)*Jhy[i,j,k]
                                ez[i,j,k] = ez[i,j,k] + (dt/(2*mu0*e0*h))*(bhy[i+1,j,k] - bhy[i-1,j,k] - bhx[i,j+1,k] + bhx[i,j-1,k]) - (dt/e0)*Jhz[i,j,k]

                                bx[i,j,k] = bx[i,j,k] - (dt/(2*h))*(ehz[i,j+1,k] - ehz[i,j-1,k] - ehy[i,j,k+1] + ehy[i,j,k-1])
                                by[i,j,k] = by[i,j,k] - (dt/(2*h))*(ehx[i,j,k+1] - ehx[i,j,k-1] - ehz[i+1,j,k] + ehz[i-1,j,k])
                                bz[i,j,k] = bz[i,j,k] - (dt/(2*h))*(ehy[i+1,j,k] - ehy[i-1,j,k] - ehx[i,j+1,k] + ehx[i,j-1,k])


#Plot#
half = numpy.around(x.size/1.9999)-1
estore = numpy.sqrt(e0)*ey[:,half,half]
bstore = bz[:,half,half]/numpy.sqrt(mu0)
pylab.figure(str(N))
pylab.title('Fields on the x-axis at t=8.34*10^-11s')
pylab.ylabel('Fields')
pylab.xlabel('x')
#pylab.xlim(-.10,.10)
#pylab.ylim(-0.6,1.0)
pylab.plot(x,estore,'b.-',x,bstore,'r.-')
pylab.savefig(str(N))