% Do zadania przyjalem plyte z grupy B, czyli:
%
%              T_3                                   T_2
%                                      ******************************
%       *+----------------------------------------------------------+
%       *|                                                          |
%       *|                                                          |
%       *|                                                          |
%    T_1*|                          T_0                             |
%       *|                                                          |
%       *|                                                          |
%       *|                                                          |
%       *|                                                          |
%       *+----------------------------------------------------------+
%                                      ******************************
%                                                    T_2
%
% Temperatura otoczenia to T_3, gwiazdki oznaczają miejscie przyłożenia
% tempreatury.
%
% Pare założeń:
%  - wspolczynnik przenikania pomiedzy "segmentami" plyty - Kw
%  - wspolczynnik przenikania plyta->powietrze - Kz
%  - wspolczynnik przenikania pomiedzy plyta a punktem grzania - Kw
%
%
% A wiec, odrobina teorii
% Rozwazmy pojedynczy segment:
%
%              Ti,(j-1)
%            +----------+
%            |          |
%   T(i-1),j |   Ti,j   | T(i+1),j
%            |          |
%            +----------+
%              Ti,(j+1)
%
% Dodatkowo mamy jeszcze temperature "nad" i "pod" plyta.
%
% W ogolnosci, mamy bardzo proste r-nie rozniczkowe wiazace ze soba
% temperature oraz strumien ciepla:
%
%                       dC * Ti,j' = sum(q)
%
% gdzie q to sa przeplywy ciepla. Dla jasnych obliczen najlatwiej je
% ponazywac od stron geograficznych: polnocny, poludniowy, wschodni,
% zachodni oraz dolny i gorny, przyjmujac ze nazwa okresla w ktora strone
% mamy przeplyw. Zaleznosc wiazaca przeplyw z temparaturami okreslamy jako
%
%       q = k*S*(T_docelowa - T_zrodla)
%
% gdzie K to jakis wspolczynnik przeplywu ciepla, a S to powierzchnia przez
% ktora nastepuje wymiana ciepla.
%
% Tworzenie rownan: wektor zmiennych stanu
% Przyjmujemy ze zmiennymi stanu sa temperatury kolejnych segmentow:
% w przypadku zadania z laboratorium mamy 1250 zmiennych stanu. Trzeba je
% tylko jakos sensownie ulozyc - moja propozycja wyglada tak: pierwsze N
% zmiennych stanu to sa temperatury kolejnych segmentow w 1 wierszu plyty,
% kolejne N zmiennych to temperatury w 2 wierszu plyty, itp. Zapis wyglada
% mniej wiecej tak:
%
% X = [T1,1 T1,2 T1,3 ... T1,n T2,1 T2,2 T2,3 ... T2,n ... Tm,1 ... Tm,n]'
%
% gdzie n - liczba segmentow w 1 wierszu (czyli - liczba kolumn), a m to
% liczba wierszy.
%
% Trzeba ustalic zaleznosc pomiedzy i,j -> indeks w wektorze zmiennych
% stanu (pamietajac, ze MATLAB liczy od 1 - my tez tak bedziemy liczyc)
% indeks = (i-1)*n + j
% warto sobie ta zaleznosc zapisac w jakiejs funkcji - naprawde latwiej
% pozniej zapisac w kodzie funkcje, niz ta zaleznosc

clear;
clc;

l = 50;
h = 25;
dx = 1;
d = dx;

T_0 = 20;
T_1 = 10;
T_2 = 70;
T_3 = 30;

Kw = 1000;
Kz = 5;
C = 5e4;

Tsym = 20;
dT = 1e-3;

% Dla ulatwienia sobie funkcji dla ODE: przyjmujemy uklad typu
%
%      X' = A*X + B*U
%
% Gdzie U jest wektorem temperatur stalych (czyli otoczenie i temperatury
% ktore przylozylismy do plyty)
global U;
U = [T_1 T_2 T_3]';

% Liczymy sobie liczbe segmentow
global n m;
n = l / dx;
m = h / dx;

% Pole powierzchni segmentu
dS = dx*dx;

% Pojemnosc segmentu
dC = C/(n*m);

% Musimy sobie policzyc, od ktorego segmentu zaczynamy grzac temperatura
% T_2:
T2_start = round(n/2);

% Generowanie macierzy A i B
% Duzo prosciej jest wygenerowac od razu cala macierz A, niz najpierw
% srodek plyty, a pozniej krawedzie - dlaczego, zobaczycie w petli

global A B;
A = sparse(n*m, n*m);
B = sparse(n*m,3);

% i - kolumna
% j - wiersz
for i=1:m
    for j=1:n
        % Dla ulatwienia: bedziemy generowac kolejne "wiersze" macierzy A
        A_ = sparse(1,n*m);
        B_ = sparse(1, 3);

        if i == 1
            if j < T2_start
                B_(3) = B_(3) + Kz*dS/dC;
                A_(indeks(i,j)) = A_(indeks(i,j)) - Kz*dS/dC;
            else
                B_(2) = B_(2) + Kw*dS/dC;
                A_(indeks(i,j)) = A_(indeks(i,j)) - Kw*dS/dC;
            end
        elseif ((i>5 && i<20) && (j>10 && j<40))  %TUTAJ
            A_(indeks(i,j)) = 0;
        elseif (i == 20 && (j>10 && j<40))
            A_(indeks(i,j)) = A_(indeks(i,j)) - Kz*dS/dC;
            B_(3) = B_(3) + Kz*dS/dC;

        else
            A_(indeks(i-1,j)) = A_(indeks(i-1,j)) + Kw*dS/dC;
            A_(indeks(i,j)) = A_(indeks(i,j)) - Kw*dS/dC;
        end

        if i == m
            if j < T2_start
                B_(3) = B_(3) + Kz*dS/dC;
                A_(indeks(i,j)) = A_(indeks(i,j)) - Kz*dS/dC;
            else
                B_(2) = B_(2) + Kw*dS/dC;
                A_(indeks(i,j)) = A_(indeks(i,j)) - Kw*dS/dC;
            end
        elseif (i == 5 && (j>10 && j<40))
            A_(indeks(i,j)) = A_(indeks(i,j)) - Kz*dS/dC;
        elseif ((i>5 && i<20) && (j>10 && j<40))  %TUTAJ
            A_(indeks(i,j)) = 0;
        else
            A_(indeks(i+1,j)) = A_(indeks(i+1,j)) + Kw*dS/dC;
            A_(indeks(i,j)) = A_(indeks(i,j)) - Kw*dS/dC;
        end

        if j == n
            B_(3) = B_(3) + Kz*dS/dC;
            A_(indeks(i,j)) = A_(indeks(i,j)) + Kz*dS/dC;
        elseif (j == 10 && (i>50 && i<20))
            A_(indeks(i,j)) = A_(indeks(i,j)) - Kz*dS/dC;
            B_(3) = B_(3) + Kz*dS/dC;
        elseif ((i>5 && i<20) && (j>10 && j<40))  %TUTAJ
            A_(indeks(i,j)) = 0;
        else
            A_(indeks(i,j+1)) = A_(indeks(i,j+1)) + Kw*dS/dC;
            A_(indeks(i,j)) = A_(indeks(i,j)) - Kw*dS/dC;
        end

        if j == 1
            B_(1) = B_(1) + Kw*dS/dC;
            A_(indeks(i,j)) = A_(indeks(i,j)) - Kw*dS/dC;
        elseif ((i>5 && i<20) && (j>10 && j<40))  %TUTAJ
            A_(indeks(i,j)) = 0;
        elseif (j == 40 && (i>5 && i<20))
            A_(indeks(i,j)) = A_(indeks(i,j)) - Kz*dS/dC;
            B_(3) = B_(3) + Kz*dS/dC;
        else
            A_(indeks(i,j-1)) = A_(indeks(i,j-1)) + Kw*dS/dC;
            A_(indeks(i,j)) = A_(indeks(i,j)) - Kw*dS/dC;
        end

        if ((i>5 && i<20) && (j>10 && j<40)) %TUTAJ
             A_(indeks(i,j)) = 0;
        else
             A_(indeks(i,j)) = A_(indeks(i,j)) - 2*Kz*dS/dC;
             B_(3) = B_(3) + 2*Kz*dS/dC;
        end;

        A(indeks(i,j),:) = A_;
        B(indeks(i,j),:) = B_;
    end
end

% Yay, mamy macierze - symulujemy? no jasne!

% Jeszcze tylko warunki poczatkowe
X0 = T_0*ones(n*m,1);

% Jakies parametry symulacji tez by bylo fajnie
TT = [0 Tsym];
opts = odeset('MaxStep', dT);

[X, T] = ode45('grzanie_plyty', TT, X0);

% Ok, no to co - tworzymy animacje ;)

% Tworzymy rysunek poczatkowy
% Teraz dla odmiany przekształcamy wektor zmiennych stanu na macierz, w
% ktorej odpowiednie wspolrzedne odpowiadaja odpowiedniemu fragmentowy
% plyty
for i = 1:m
    for j = 1:n
        ZZ(i,j) = T(1,indeks(i,j));
    end
end;
h1 = mesh(ZZ);
axis([0 50 0 25 0 70]);
% No i nastepne
pause;
for k = 1:10:length(T)
    for i = 1:m
        for j = 1:n
            ZZ(i,j) = T(k,indeks(i,j));
        end
    end;

    % Podmieniamy nasze temperatury
    set(h1,'zdata',ZZ);
    pause(0.01);
end;