import numpy as np
import time

start = time.time()
t_steps = 100 # steps to take 
size_x = 20 # cells in x
size_y = 20 # cells in y

C = 343 #m/s = speed of sound
freq = C/(4*size_x) # standing wave for closed cylinder

# create epsilon
eps = np.ones((size_x, size_y))

#initialize 
Ex = np.zeros((size_x,size_y))
Ey = np.zeros((size_x, size_y))
Hz = np.zeros((size_x, size_y))

#spatial steps
dx = dy = 1

#time steps
dt = 1
c = 1

## courant constant
#cc = c*dt/min(dx, dy)
cc = 0.5
d = 0.99
for t in range(t_steps):
    #print('t= '+str(t))
    #update Hz, Ex, Ey
    # update H components
    deriv_y = np.zeros((size_x, size_y))
    deriv_x = np.zeros((size_x, size_y))
    for j in range(size_x):
        for k in range(size_y):
            indx = j + 1; indy = k + 1
            if (indx > size_x - 1):
                indx = 0
            if (indy > size_y - 1):
                indy = 0
            deriv_y[k, j] = (Ex[indy, j] - Ex[k, j])
            deriv_x[k, j] = (Ey[k, indx] - Ey[k, j])

    Hz -= (deriv_x - deriv_y)
    Hz[int(size_x/2), int(size_y/2)] -= np.sin(freq/10 * 2 * np.pi * t/30)#np.sin(2*np.pi*t/30)*(dt)

    #update E components
    Curl_x = np.zeros((size_x, size_y))
    Curl_y = np.zeros((size_x, size_y))
    for j in range(size_x): #y axis
        for k in range(size_y): #x axis
            indx = j - 1; indy = k - 1
            if(indx < 0):
                indx = size_x - 1
            if(indy < 0):
                indy = size_y - 1
            deriv_y = -(Hz[indy, j] - Hz[k, j])
            deriv_x = -(Hz[k, indx] - Hz[k, j])
            Curl_x[k, j]= (cc) * (deriv_x) #curl Hz, dy(Hz)
            Curl_y[k, j] = (cc) * (deriv_y)  #cul Hz dx(Hz)
    Ex += Curl_y
    Ey -= Curl_x

end = time.time()
print("time elapsed: ") 
print(end-start)