# coding: utf-8

# In[21]:

# 3 Set 3
# 3.1 Eigenmodes of drums or membranes of dierent shapes


# In[45]:

import numpy as np
import scipy.sparse.linalg as sc
import scipy.linalg as scnorm
import matplotlib.pyplot as plt
size = 50
setting = 1

#1. create the matrix a
def matrix(size, setting):
    dx = size / (size - 1)
    B = np.identity(size ) * (-4)
    counter = 0
    sizeCounter = size
    for i in range(1, size):
        B[i, counter] = 1
        B[counter, i] = 1
        counter += 1
    # B[0,:]=0
    # B[:,0]=0
    # B[size-1,:]=0
    # B[:,size-1]=0

    I=np.identity(size)
    # I[0,:]=0
    # I[:,0]=0
    # I[size-1,:]=0
    # I[:,size-1]=0
    if setting == 1:
        middle = np.floor(len(B) / 2)
        for i in range(0, size):
            for j in range(0, size):
                if( np.sqrt((i - middle)**2 + (j - middle)**2) >  middle): B[i, j] = 0
                if( np.sqrt((i - middle)**2 + (j - middle)**2) >  middle): I[i, j] = 0
    else:
        middle = np.floor(len(B) / 2)
        L = 20
        N = 10
        for i in range(0, size):
            for j in range(0, size):
                if(i < L or i > size - L or j > size - N or j < N): B[i, j] = 0
                if(i < L or i > size - L or j > size - N or j < N): I[i, j] = 0
    return B,I

#used to find the middle of the matrix

#2. create the vector b with starting conditions
def vector(size):
    b = np.zeros([size, size])
    vertical = np.floor(len(b) / 2)
    # point06 = np.floor((len(b) - vertical) * 0.3)
    # point12 = np.floor((len(b) - vertical) * 0.6)
    # b[int(point06), int(point12)] = 10
    b[int(vertical * 1.2), int(vertical / 1.6)] = 1
    bcolumn = b.T.reshape(size**2, 1)

    return bcolumn


B, I = matrix(size, setting)
A = np.zeros((size, size))

M = []
for i in range(size):
    if i == 0:
        M.append(np.concatenate([ B, I, np.tile(A, (1, size-2)) ], axis = -1 ))
    elif i == size-1:
        M.append(np.concatenate([ np.tile(A, (1, size-2)), I, B ], axis = -1))
    else:
        M.append(np.concatenate([ np.tile(A, (1, i-1)), I, B, I, np.tile(A, (1, size-i-2))], axis = -1))
M = np.concatenate(M, axis = 0)

#circle domain
if setting == 1:
    middle = np.floor(len(B) / 2)
    for i in range(size):
        for j in range(size):
            if( np.sqrt((i - middle)**2 + (j - middle)**2) >  middle): B[i, j] = 0
            if( np.sqrt((i - middle)**2 + (j - middle)**2) >  middle): I[i, j] = 0
else:
    L = 20
    N = 10
    for i in range(size):
        for j in range(size):
            if(i < L or i > size - L or j > size - N or j < N): B[i, j] = 0
            if(i < L or i > size - L or j > size - N or j < N): I[i, j] = 0

# make the matrix not singular
for i in range(size**2):
    if M[i,i]==0: M[i,i]=1



# M = M / (1/size)**2
b = vector(size)
solution = sc.spsolve(M, b)
solution = solution.reshape(size, size)
middle = np.floor(len(solution) / 2)
if setting == 1:
    for i in range(size):
        for j in range(size):
            if( np.sqrt((i - middle)**2 + (j - middle)**2) >  middle): solution[i, j] = np.nan
else:
    L = 20
    N = 10
    for i in range(size):
        for j in range(size):
            if(i < L or i > size - L or j > size - N or j < N): solution[i, j] = np.nan

#print(M)
# plt.imshow(M,interpolation='none')
# plt.show()