# -*- coding: utf-8 -*-
"""
Created on Sun Oct 29 15:35:23 2017

@author: jared's account
"""

# point_charge .py - Iterative solution of 2D Poisson equation
import numpy as np
import matplotlib
import matplotlib . pyplot as plt
# System size and resolution
L = 2.0
N = 21
h = L/(N-1)
# Make the charge density matrix
q = 10
rho0 = q/h**2
rho = np.zeros((N,N))
rho[N/4 , N/2] = rho0 # Intentional integer division
rho[3*N/4 , N/2] = -1*rho0 # Intentional integer division

# Make the initial guess for solution matrix
V = np.zeros((N,N))
V[:,N-1]=1
V[:0]=-1
# Solver
iterations = 0
eps = 1e-8 # Convergence threshold
error = 2*eps # Dummy error to enter the first loop
e=8.854e-12
while iterations < 1e4 and error > eps :
    V_temp = np.copy ( V )
    error = 0.0 # We make this accumulate in the loop
    for j in range (1 ,N -1): # Avoid updating the boundaries
        for i in range (1 ,N -1):
            V [i , j ] = 0.25*( V_temp[i+1,j] + V_temp[i-1,j] +
            V_temp[i,j-1] + V_temp [i,j+1] + rho[i, j ]* h**2/(4*e))
            error += np.abs(V[i,j] - V_temp[i,j])
    iterations += 1
    error /= float(N**2)
    if iterations%200 == 0:
        print("iterations = " , iterations)
   
# Define arrays used for plotting
x = np.linspace (0 ,L , N )
y = np.linspace (0 ,L , N )
# meshgrid () by default orders the arrays differently from what
# I prefer . This way , (X[i,j] ,Y[i,j]) equals (x[i] ,y[j])
Y , X = np.meshgrid(y , x )
# Plot the potential V
#matplotlib.rcParams['ytick.direction'] = ’out’
CS = plt.contour(X ,Y ,V ,500) # Make a contour plot
# plt. clabel (CS , inline =1 , fontsize =10) # Contour labels
plt.title('Point charge in the cavity of a grounded box')
CB = plt.colorbar(CS, shrink=0.8 , extend = 'both')
plt.show()



h=L/(N-1)
Exx,Eyy = np.gradient(V,h/2,h/2)
plt.quiver(X,Y,-1*Exx,-1*Eyy)
plt.title('Fields lines of Electric Field')
plt.show()
#fig = plt.figure(11)
#
#ax = fig.gca()
#ax.quiver(X,Y, -1*Exx, -1*Eyy, color='r',
#          angles='xy', scale_units='inches', linewidths=2)