def Poisson_Solver_Neumann(u, v, Nx, Ny, dt, dx, dy, T, xp, yp):
"""Solves the 2D Poisson equation implicitly on a staggered grid
using Neumann Boundary Conditions

Params:
------
u, v           2D array of float, x and y velocities
Nx, Ny         float, Number of segments in x and y
dt             float, time step size
T              float, current time
X, Y           2D array of float, meshgrid

Returns:
-------
ANeum        2D array of float, A matrix with Neumann conditions
f_RHSn       1D array of float, f(x,y) for Neumann conditions
"""
#Building A
ANeum = np.zeros((Nx*Ny,Nx*Ny),dtype=float)

a = Nx * Ny
b = Nx * Ny
c_int = -4 * (Nx * Ny)
c_edge = -3 * (Nx * Ny)
c_corner = -2 * (Nx * Ny)
d = Nx * Ny
e = Nx * Ny

#Set corner points
ANeum[0,0] = c_corner
ANeum[-1,-1] = c_corner
ANeum[Nx*Ny-Ny,Nx*Ny-Ny] = c_corner
ANeum[Ny-1,Ny-1] = c_corner

#Set edges in first block
for j in range(1,Ny-1):
    ANeum[j,j] = c_edge
    j +=j

#Set edges in last block
for j in range((Nx*Ny)-Ny,(Nx*Ny)-2):
    ANeum[j+1,j+1] = c_edge
    j +=j

#Set edges along main diagonal except for first block
for j in range(Nx+1,Ny*Nx):
    if j % Nx ==0:
        ANeum[j-1,j-1] = c_edge
    j +=j

#Set edges on main diagonal except for last block
for j in range(Nx,(Ny*Nx)-Nx):
    if j % Nx ==0:
        ANeum[j,j] = c_edge
    j +=j

#Second diagonal above and below diagonal
for j in range(Ny,Ny*Nx):
    ANeum[j,j-Ny] = a
    ANeum[j-Nx,j] = e
    j +=j

#first diagonal below main diagonal
for j in range(1,Ny*Nx):
    if j % Ny ==0:
        ANeum[j,j-1] = 0
    else:
        ANeum[j,j-1] = b
    j +=j

#first diagonal above main diagonal
for j in range(0,Ny*Nx):
    if j % Nx ==0:
        ANeum[j-1,j] = 0
    else:
        ANeum[j-1,j] = d
    j +=j

#Main Diagonal
for j in range(0,Ny*Nx):
    if ANeum[j,j] ==0:
        ANeum[j,j] = c_int

#Setting random point to Dirichlet
ANeum[0,:] = 0
ANeum[0,0] = 1

return ANeum

for t in range(0,nt):
un = np.empty_like(u)
vn = np.empty_like(v)
#copy info from last loop, needed for current loop into variable_n
un = u.copy()
vn = v.copy()
#Pn = p.copy()
T = T_0 + t*dt

####STEP 1
### t -----> t + (1/3)*dt
un, vn = BCs(un, vn)
G1 = Build_G1(un, vn, dx, dy, nu)
G2 = Build_G2(un, vn, dx, dy, nu)
u1 = un + (dt*G1) #put 1/3 back in
v1 = vn + (dt*G2)
u1, v1 = BCs(u1, v1)
u_div = (u1[1:-1,2:-1] - u1[1:-1,1:-2])/dx
v_div = (v1[2:-1,1:-1] - v1[1:-2,1:-1])/dy
#Solve for Pressure Field
b1 = np.zeros((Ny+2,Nx+2), dtype=float)
b1[1:-1,1:-1] = (u1[1:-1,2:-1] - u1[1:-1,1:-2])/dx +
    (v1[2:-1,1:-1] - v1[1:-2,1:-1])/dy
b1 = np.reshape(b1[1:-1,1:-1], Ny*Nx)
A1 = Poisson_Solver_Neumann(u1, v1, Nx, Ny, dt, dx, dy, T, xp, yp)
b1 = b1[:]*(1.0/dt)
temp1 = np.linalg.solve(A1,b1)

p_star = np.reshape(temp1, (Ny,Nx))


F1_Pstar, F2_Pstar = Pressure_Terms(p_star, dx, dy)

#calc predictor velocities
u_star = u1.copy()
v_star = v1.copy()

u_star[1:-1,2:-2] = u1[1:-1,2:-2] - (dt*F1_Pstar)
v_star[2:-2,1:-1] = v1[2:-2,1:-1] - (dt*F2_Pstar)

u = u_star.copy()
v = v_star.copy()