# stacionarni reseni 1-D rovnice vedeni tepla
# (jednorozmerna tyc, nebo nekonecne siroka
# deska konecne tloustky)
#
import numpy as np
import scipy
import matplotlib.pyplot as plt

L= 1.0        # delka tyce
kmax= 1000    # diskretizace souradnice
h= L/kmax     # element delky
hh= 1.0/(h*h) # pre-faktor v matici 
x= np.linspace(0.0,L,kmax+1) # ekvidistantni diskretizace osy x

xk= np.zeros(kmax)   # pomocne souradnice (ve stredech)
kappa= np.ones(kmax) # koefficient vedeni tepla - konstantni hodnota
for k in range(kmax):
  xk[k]= 0.5*(x[k]+x[k+1])

# ruzne hodnoty kappa(x) - zde napr. tri vrstvy
# kappa[            0:int(  kmax/3)]= 1.0
# kappa[int(  kmax/3):int(2*kmax/3)]= 100.0
# kappa[int(2*kmax/3):kmax         ]= 1.0


T0= 20.0 # predepsana teplota v bode 0
T1= 0.0  # predepsana teplota  v bode L

T= np.zeros((kmax+1))   # vektor reseni
b= np.zeros(kmax+1)     # vektor prave strany
A= np.zeros((3,kmax+1)) # tri-diagonalni matice v uspornem ulozeni
# viz napoveda k scipy.linalg.solve_banded

# prvni radek matice - pouze diagonala
A[1,0]=    1.0
# posledni radek matice - pouze diagonala
A[1,kmax]= 1.0

# vnitrni radky matice podle usporneho schematu A[1+i-j]= a[i,j]
for k in range(1,kmax):
  A[0,k+1]= - hh* kappa[k]               # horni  diagonala
  A[1,k  ]= + hh*(kappa[k-1]+kappa[k])   # hlavni diagonala
  A[2,k-1]= - hh* kappa[k-1]             # spodni diagonala

# leva hranice
b[0]= T0
# prava hranice
b[kmax]= T1

# kontrolni vypis matice a prave strany
# print(A)
# print(b)


# reseni soustavy
T[:]= scipy.linalg.solve_banded((1,1),A,b,check_finite=True)

# kontrolni vypis reseni
# print(T)

# obrazek
plt.plot(x,T)
plt.plot(xk,kappa) # volitelne kresleni prubehu kapppa do stejneho obrazku
plt.title('Prubeh stacionarni teploty')
plt.xlabel('x')
plt.ylabel('T')
plt.show()