import numpy as np
import time
start = time.time()

t_steps = 100 # steps to take 
size_x = 1000 # cells in x
size_y = 1000 # cells in y


C = 343 #m/s = speed of sound
CA_DX = 0.001 #m = cell size
CA_DT = 1/C 

pos_x = int(size_x/2)
pos_y = int(size_y/2)

damping = 0.99
omega = 3 / (2 * np.pi)

freq = C/(4*size_x) # standing wave for closed cylinder
amplitude = 10

# update rules for pressure and velocity based on neighbouring cells
pressure = [[0.0 for x in range(size_x)] for y in range (size_y)]
V = [[[0.0, 0.0, 0.0, 0.0] for x in range(size_x)] for y in range (size_y)]
c = 1 / np.sqrt(2) 

def update_V():
    # V(x,dt) = V(x,t) - {P(x + dx, t) - P(x,t)} 
    for i in range(size_y):
        for j in range(size_x):
            #p = pressure;
            #V = velocity
            cell_pressure = pressure[i][j]
            V[i][j][0] = V[i][j][0] + cell_pressure - pressure[i - 1][j] if i > 0 else cell_pressure
            V[i][j][1] = V[i][j][1] + cell_pressure - pressure[i][j + 1] if j < size_x -1 else cell_pressure
            V[i][j][2] = V[i][j][2] + cell_pressure - pressure[i + 1][j] if i < size_y -1 else cell_pressure
            V[i][j][3] = V[i][j][3] + cell_pressure - pressure[i][j - 1] if j > 0 else cell_pressure
            
def update_P():
    for i in range(size_y):
        for j in range(size_x):
            pressure[i][j] -= c**2 * damping * sum(V[i][j])
            
def step(i):
    pressure[pos_y][pos_x] = np.sin(freq/10 * 2 * np.pi * i)
    update_V()
    update_P()

for i in range(t_steps):
    step(i)

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