# 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 = 100

#1. create the matrix a
def matrix(size):
    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


    return B,I

#used to find the middle of the matrix

# def middle(B, dx):
#     B = B / dx**2
#     horizontal = np.floor(len(B) / 2)
#     vertical = horizontal
#     for i in range(len(B)):
#         for j in range(len(B)):
#             if( np.sqrt((i - horizontal)**2 + (j - vertical)**2) >  horizontal): B[i, j] = 0
#             if(i == j and B[i, j] == 0): B[i, j] = 1

#     return B, vertical, horizontal

#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)] = 1
    bcolumn = b.T.reshape(size**2, 1)

    return bcolumn
#
# #3. solve the problem
# solution = sc.spsolve(matrix(size), vector(size))
# solution.reshape(size, size)


B, I = matrix(size)
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)


# make the matrix not singular
for i in range(size**2):

    for j in range(size**2):

        if i==j:
            if M[i,j]==0 :

                M[i,j]=1

M = M / (1/size)**2

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