function heat_mil
%PROGRAM DO ANALIZY STACJONARNEGO PRZEP£YWU CIEP£A MES
%AUTOR: S£AWOMIR MILEWSKI, PK WIL L5
%mail: slawek@L5.pk.edu.pl
%data utworzenia: paŸdziernik 2010
%data modyfikacji: styczeñ 2011
global stopien ns met;
global dt nt t0 T0 Tp Tn analiza sp_pocz wzor_pocz schemat tk czas cw;
close all
clc
cw = 1;
amp_zab = 0;
t0 = 0;
czas = 0;
T0 = [];
T = [];
R = [];
nt = 10;
dt = 0.02;
tk = 0;
met = 1;
sfor = 1;
sp_pocz = 2;
wzor_pocz = '0';
schemat = 2;
analiza = 1;
ksztalt = 2;
stopien = 1;
liczgwiaz = 9;
par_wyg = 0.01;
%par_wyg = 0;
wygladz = 0;
kx = 0;
ky = 0;
f = '0';
h = 1;
siat = 0;
nh = 1;
obl = 0;
naz_plik = '';
Kr = [];
WBQk = [];
WBTk = [];
Incs = [];
Xk = [];
Yk = [];
dane = 0;
zapis = 1;
if dane==1
    Xk = [0 2 5 0 2 5 0 2];
    Yk = [0 0 0 2 2 2 3 3];
    Kr = [1 2;2 3;3 6;6 5;5 8;8 7;7 4;4 1];
    WBTk = [1 10;2 10;7 10;8 10];
    WBQk = [4 -12;5 -12];
    siat = 1;
    Incs = [1 2 5 4; 4 5 8 7; 2 3 6 5];
    ksztalt=2;
end
nk = length(Xk);
ilk = size(Kr,1);
ilkQ = size(WBQk,1);
ilkT = size(WBTk,1);
Nes = size(Incs,1);
Xs = Xk;
Ys = Yk;
ns = nk;
nspop = 0;
lwe = 4;
if ksztalt==1
    lwe = 3;
end
ilop = 16;
op = 1;
while (op<ilop)
    s1 = 'WARUNKI BRZEGOWE (TEMPERATURA, STRUMIEÑ CIEP£A)';
    if analiza==2
        s1 = 'WARUNKI BRZEGOWO-POCZ¥TKOWE (TEMPERATURA, STRUMIEÑ CIEP£A)';
    end
    op = menu(sprintf('\t H E A T M I L \n - ANALIZA MES i (B)MRS \n USTALONEGO I NIEUSTALONEGO PRZEP£YWU CIEP£A \n\n autor: S³awomir Milewski \n http://www.L5.pk.edu.pl/~slawek \n slawek@L5.pk.edu.pl'),...
        'DEFINIOWANIE KRAWÊDZI OBSZARU','DEFINIOWANIE MATERIA£U',...
        'GENERATOR CIEP£A W OBSZARZE',s1,'RODZAJ ANALIZY I PARAMETRY METOD MES / (B)MRS','GENERACJA RÊCZNA SIATKI MES',...
        'GENERACJA AUTOMATYCZNA SIATKI MES/ WÊZ£ÓW (B)MRS','RYSUNEK: SCHEMAT CIEPLNY KONSTRUKCJI','RYSUNEK: SIATKA MES / WÊZ£Y BMRS','OBLICZENIA',...
        'WYNIKI I ICH WIZUALIZACJA','ZAPIS DANYCH DO PLIKU','ZAPIS WYNIKÓW DO PLIKU','ODCZYT DANYCH Z PLIKU','WYCZYSZCZENIE BAZY DANYCH','WYJŒCIE');
    WBQ = zeros(ilkQ,3);
    WBT = zeros(ilkT,3);
    for i=1:ilkT
        nr = WBTk(i,1);
        n1 = Kr(nr,1);
        n2 = Kr(nr,2);
        WBT(i,:) = [n1 n2 WBTk(i,2)];
    end
    for i=1:ilkQ
        nr = WBQk(i,1);
        n1 = Kr(nr,1);
        n2 = Kr(nr,2);
        WBQ(i,:) = [n1 n2 WBQk(i,2)];
    end
    if (op==ilop | op==ilop-1) & (zapis==0)
        war = questdlg('NIEZAPISANE DANE I WYNIKI ZOSTAN¥ UTRACONE! OPUŒCIÆ PROGRAM?',...
            'PROŒBA O POTWIERDZENIE','TAK','NIE','NIE');
        switch war
            case 'NIE'
                op = 0;
        end
    end
    if op>=1 & op<=7
        zapis = 0;
    end
    if op==1
        siat = 0;
    end
    if op<=7
        close all;
        obl = 0;
    end
    if op<=7 | op>=10
        try
            close(msg);
        catch
        end
    end
    if op==15
        close all
        kx = 0;
        ky = 0;
        h = 1;
        Xs = [];
        Ys = [];
        ns = 0;
        nspop = 0;
        Xk = [];
        Yk = [];
        nk = 0;
        ilk = 0;
        Kr = [];
        Nes = 0;
        Incs = 0;
        f = 0;
        ilkQ = 0;
        ilkT = 0;
        WBQk = [];
        WBTk = [];
        siat = 0;
        obl = 0;
        nh = 1;
        schemat = 2;
        analiza = 1;
        ksztalt = 2;
        stopien = 1;
        met = 1;
        uiwait(msgbox('USUNIÊTO WSZYSTKIE DANE','INFO','help','modal'));
        zapis = 1;
    end
    if op==1
        close all;
        op1 = 1;
        ilop1 = 4;
        while op1<ilop1
            op1 = menu('DEFINIOWANIE KRAWÊDZI OBSZARU','PUNKTY','LINIE','USUÑ DANE GEOMETRYCZNE','POWRÓT');
            if op1==3
                Xk = [];
                Yk = [];
                ilk = 0;
                Kr = [];
                nk = 0;
                Xs = [];
                Ys = [];
                ns = 0;
                nspop = 0;
                WBQk = [];
                WBTk = [];
                ilkT = 0;
                ilkQ = 0;
                Incs = [];
                Nes = 0;
            end
            if op1==1
                op2 = 1;
                ilop2 = 5;
                while op2<ilop2
                    s = '';
                    for i=1:nk
                        t = ['i = ' num2str(i) ', x = ' num2str(Xk(i)) ', y = ' num2str(Yk(i))];
                        if i==1
                            s = str2mat(t);
                        else
                            s = str2mat(s,t);
                        end
                    end
                    rysunek(1,nk,Xk,Yk,ilk,Kr);
                    msg = msgbox(s,'PUNKTY KRAWÊDZI','help','replace');
                    op2 = menu('DEFINIOWANIE PUNKTÓW KRAWÊDZI','DODAJ PUNKT','USUÑ PUNKT','POPRAW PUNKT','USUÑ WSZYSTKIE PUNKTY','POWRÓT');
                    options.Resize='on';
                    options.WindowStyle='normal';
                    options.FontSize = 20;
                    if op2==2 | op2==4
                        if ilk>0
                            uiwait(msgbox('TA OPERACJA WYMAGA USUNIÊCIA LINII BRZEGOWEJ!','UWAGA','warn','modal'));
                            op2 = 5;
                        end
                    end
                    if op2>=1 & op2<=4
                        if Nes>0
                            uiwait(msgbox('TA OPERACJA WYMAGA USUNIÊCIA SIATKI MES / WÊZ£ÓW BMRS','UWAGA','warn','modal'));
                            op2 = 5;
                        end
                    end
                    if op2==1
                        pr = {['x(' num2str(nk+1) ') = '],['y(' num2str(nk+1) ') = ']};
                        odp = inputdlg(pr,'PUNKTY',1,{'',''},options);
                        if isempty(odp)==0
                            bl = 0;
                            try
                                [x y] = deal(odp{:});
                                x = eval(x);
                                y = eval(y);
                            catch
                                bl = 1;
                            end
                            if bl==0
                                nk
                                Xk
                                Yk
                                nk = nk + 1;
                                Xk1 = Xk;
                                Yk1 = Yk;
                                Xk = zeros(nk,1);
                                Yk = zeros(nk,1);
                                Xk(1:nk-1) = Xk1;
                                Yk(1:nk-1) = Yk1;
                                Xk(nk) = x;
                                Yk(nk) = y;
                            else
                                uiwait(msgbox('B£¥D - Z£A WARTOŒÆ','B£¥D','error','modal'));
                            end
                        end
                    end
                    if op2==2
                        pr = {['nr (1..' num2str(nk) ') = ']};
                        odp = inputdlg(pr,'PUNKTY',1,{''},options);
                        if isempty(odp)==0
                            x = deal(odp{:});
                            nr = 0;
                            bl = 0;
                            try
                                nr = eval(x);
                            catch
                                bl = 1;
                            end
                            if nr>0 & nr<=nk
                                for j=nr:nk-1
                                    Xk(j) = Xk(j+1);
                                    Yk(j) = Yk(j+1);
                                end
                                Xk(nk) = [];
                                Yk(nk) = [];
                                nk = nk-1;
                            else
                                uiwait(msgbox('B£¥D - Z£Y NUMER PUNKTU','B£¥D','error','modal'));
                            end
                        end
                    end
                    if op2==3
                        pr = {['nr (1..' num2str(nk) ') = ']};
                        odp = inputdlg(pr,'PUNKTY',1,{''},options);
                        if isempty(odp)==0
                            x = deal(odp{:});
                            bl = 0;
                            try
                                nr = eval(x);
                            catch
                                bl = 1;
                            end
                            if (nr>0 & nr<=nk) & bl==0
                                pr = {['x(' num2str(nr) ') = '],['y(' num2str(nr) ') = ']};
                                odp = inputdlg(pr,'PUNKTY',1,{num2str(Xk(nr)) num2str(Yk(nr))},options);
                                if isempty(odp)==0
                                    [x y] = deal(odp{:});
                                    bl = 0;
                                    try
                                        z = eval(x);
                                        z = eval(y);
                                    catch
                                        bl = 1;
                                    end
                                    if bl==0
                                        Xk(nr) = eval(x);
                                        Yk(nr) = eval(y);
                                    else
                                        uiwait(msgbox('B£¥D - Z£A WARTOŒÆ','B£¥D','error','modal'));
                                    end
                                end
                            else
                                uiwait(msgbox('B£¥D - Z£Y NUMER PUNKTU','B£¥D','error','modal'));
                            end
                        end
                    end
                    if op2==4
                        Xk = [];
                        Yk = [];
                        nk = 0;
                    end
                end
                close all;
                try
                    close(msg);
                catch
                end
            end
            if op1==2
                if nk<4
                    uiwait(msgbox('B£¥D - ZA MA£O PUNKTÓW','B£¥D','error','modal'));
                else
                    op2 = 1;
                    ilop2 = 5;
                    while op2<ilop2
                        rysunek(1,nk,Xk,Yk,ilk,Kr);
                        if op2==4
                            s = '';
                            for i=1:ilk
                                t = ['i = ' num2str(i) ', n1 = ' num2str((Kr(i,1))) ', n2 = ' num2str((Kr(i,2)))];
                                if i==1
                                    s = str2mat(t);
                                else
                                    s = str2mat(s,t);
                                end
                            end
                            uiwait(msgbox(s,'LINIE KRAWÊDZI','help','replace'));
                        end
                        op2 = menu('DEFINIOWANIE LINII BRZEGOWYCH','RÊCZNE TWORZENIE LINII BRZEGOWEJ','AUTOMATYCZNE TWORZENIE LINII BRZEGOWEJ',...
                            'USUNIÊCIE LINII BRZEGOWEJ','LISTA KRAWÊDZI','POWRÓT');
                        if op2==1
                            if ilk==0
                                ilk1 = 0;
                                koniec = 0;
                                Kr1 = zeros(nk,2);
                                while koniec==0
                                    if ilk1==0
                                        pr = {'PUNKT POCZ¥TKOWY'};
                                    else
                                        pr = {['KOLEJNY PUNKT (nr = ' num2str(ilk1+1) ', pop = ' num2str(nrp) ') = ']};
                                    end
                                    odp = inputdlg(pr,'LINIE',1,{''});
                                    if isempty(odp)==0
                                        bl = 0;
                                        try
                                            nr = eval(deal(odp{:}));
                                        catch
                                            bl = 1;
                                        end
                                        if bl==1 | (nr<1 | nr>nk)
                                            koniec = 2;
                                        else
                                            if ilk1==0
                                                nr0 = nr;
                                            end
                                            ilk1=ilk1+1;
                                            if ilk1>1
                                                Kr1(ilk1-1,:) = [nrp nr];
                                            end
                                            nrp = nr;
                                            if nr==nr0 & ilk1>1
                                                koniec = 1;
                                            end
                                        end
                                    else
                                        koniec = 2;
                                    end
                                end
                                if koniec==1
                                    ilk = ilk1-1;
                                    Kr = Kr1;
                                end
                            else
                                uiwait(msgbox('TA OPERACJA WYMAGA USUNIÊCIA LINII BRZEGOWEJ!','UWAGA','warn','modal'));
                            end
                        end
                        if op2==2
                            if ilk==0
                                ilk = nk;
                                Kr = zeros(ilk,2);
                                for i=1:ilk
                                    n1 = i;
                                    n2 = i+1;
                                    if i==ilk
                                        n2 = 1;
                                    end
                                    Kr(i,:)= [n1 n2];
                                end
                            else
                                uiwait(msgbox('TA OPERACJA WYMAGA USUNIÊCIA LINII BRZEGOWEJ!','UWAGA','warn','modal'));
                            end
                        end
                        if op2==3
                            if ilkQ>0 | ilkT>0
                                uiwait(msgbox('TA OPERACJA WYMAGA USUNIÊCIA WSZYSTKICH WARUNKÓW BRZEGOWYCH Z KRAWÊDZI!','UWAGA','warn','modal'));
                            else
                                ilk = 0;
                                Kr = [];
                            end
                        end
                    end
                    close all;
                end
            end
        end
    end
    if op==2
        op1 = 1;
        ilop1 = 3;
        while op1<ilop1
            op1 = menu('STA£E MATERIA£OWE','STA£A PRZEWODNICTWA k_x',...
                'STA£A PRZEWODNICTWA k_y','POWRÓT');
            if op1<ilop1
                if op1==1
                    naz = 'kx';
                    z = kx;
                end
                if op1==2
                    naz = 'ky';
                    z = ky;
                end
                pr = {[naz ' = ']};
                odp = inputdlg(pr,'STA£E',1,{num2str(z)});
                if isempty(odp)==0
                    bl = 0;
                    try
                        z = eval(deal(odp{:}));
                    catch
                        bl = 1;
                    end
                    if bl==0
                        if op1==1
                            kx = z;
                        end
                        if op1==2
                            ky = z;
                        end
                    else
                        uiwait(msgbox('B£¥D - Z£A WARTOŒÆ','B£¥D','error','modal'));
                    end
                end
            end
        end
    end
    if op==3
        pr = {'INTENSYWNOŒÆ STRUMIENIA f(x,y,t) = '};
        odp = inputdlg(pr,'GENERATOR',1,{f});
        if isempty(odp)==0
            bl = 0;
            try
                f0 = (deal(odp{:}));
                x = 1;
                y = 1;
                t = 1;
                f0w = eval(f0);
            catch
                bl = 1;
            end
            if bl==1
                uiwait(msgbox('B£¥D - Z£A WARTOŒÆ','B£¥D','error','modal'));
            end
            if bl==0
                f = f0;
                il = 0;
                for i=1:length(f0)
                    if f0(i)~=' '
                        il = il + 1;
                        f(il) = f0(i);
                    end
                end
                f(il+1:end) = [];
            end
        end
    end
    if op==4
        if ilk<4
            uiwait(msgbox('B£¥D - ZA MA£O KRAWÊDZI BRZEGOWYCH','B£¥D','error','modal'));
        else
            ilop1 = 9;
            op1 = 1;
            while op1<ilop1
                WBQ = zeros(ilkQ,3);
                WBT = zeros(ilkT,3);
                for i=1:ilkT
                    nr = WBTk(i,1);
                    n1 = Kr(nr,1);
                    n2 = Kr(nr,2);
                    WBT(i,:) = [n1 n2 WBTk(i,2)];
                end
                for i=1:ilkQ
                    nr = WBQk(i,1);
                    n1 = Kr(nr,1);
                    n2 = Kr(nr,2);
                    WBQ(i,:) = [n1 n2 WBQk(i,2)];
                end
                if op1==7
                    s = 'KRAWÊDZIE Z TEMPERATUR¥';
                    for i=1:ilkT
                        t = ['nr = ' num2str(WBTk(i,1)) ', T = ' num2str((WBTk(i,2)))];
                        s = str2mat(s,t);
                    end
                    s = str2mat(s,' ');
                    s = str2mat(s,'KRAWÊDZIE ZE STRUMIENIEM CIEP£A');
                    for i=1:ilkQ
                        t = ['nr = ' num2str(WBQk(i,1)) ', q = ' num2str((WBQk(i,2)))];
                        s = str2mat(s,t);
                    end
                    uiwait(msgbox(s,'WARUNKI BRZEGOWE','help','replace'));
                end
                rysunek(1,nk,Xk,Yk,ilk,Kr,1,WBTk,WBQk,Nes,Incs,Xs,Ys,lwe);
                op1 = menu('WARUNKI BRZEGOWE','NOWY WARUNEK NA KRAWÊDZI','KOREKTA WARUNKU Z WYBRANEJ KRAWÊDZI','USUNIÊCIE WARUNKU Z WYBRANEJ KRAWÊDZI',...
                    'USUNIÊCIE WSZYSTKICH TEMPERATUR',...
                    'USUNIÊCIE WSZYSTKICH STRUMIENI',...
                    'USUNIÊCIE WSZYSTKICH WARUNKÓW','LISTA WARUNKÓW','DEFINICJA I WIZUALIZACJA TEMPERATURY POCZ¥TKOWEJ','POWRÓT');
                if op1==8
                    if analiza==1
                        uiwait(msgbox('WARUNKI POCZ¥TKOWE DOTYCZ¥ TYLKO PRZEP£YWU NIEUSTALONEGO','B£¥D','error','modal'));
                    end
                    if analiza==2
                        il = 0;
                        if sp_pocz==1
                            bl = 0;
                            if Nes==1 | length(T)~=ns | isempty(T)==1
                                uiwait(msgbox('B£¥D - NIEAKTUALNE WYNIKI','B£¥D','error','modal'));
                                bl=1;
                            end
                            if bl==0
                                T0 = T;
                                il = 1;
                            end
                        end
                        if sp_pocz==2
                            il = 0;
                            pr = {'FUNKCJA TEMPERATURY POCZ¥TKOWEJ T0(x,y,t) = '};
                            odp = inputdlg(pr,'FUNKCJA TEMPERATURY POCZ¥TKOWEJ',1,{wzor_pocz});
                            if isempty(odp)==0
                                bl = 0;
                                try
                                    x = 0;
                                    y = 0;
                                    t = t0;
                                    z = eval(deal(odp{:}));
                                catch
                                    bl = 1;
                                end
                                if bl==1
                                    uiwait(msgbox('B£¥D - Z£A WARTOŒÆ','B£¥D','error','modal'));
                                else
                                    wzor_pocz0 = (deal(odp{:}));
                                    
                                    wzor_pocz = wzor_pocz0;
                                    il = 0;
                                    for i=1:length(wzor_pocz0)
                                        if wzor_pocz0(i)~=' '
                                            il = il + 1;
                                            wzor_pocz(il) = wzor_pocz0(i);
                                        end
                                    end
                                    wzor_pocz(il+1:end) = [];
                                    
                                    il = 1;
                                    if ilkT>0
                                        il = 0;
                                        dec = questdlg('WARTOŒCI TEMPERATUR NA KRAWÊDZIACH MOG¥ ULEC ZMIANIE! KONTYNUOWAÆ?','ZMIANA WARTOŒCI TEMPERATUR BRZEGOWYCH','TAK','NIE','TAK');
                                        if strcmp(dec,'TAK')==1
                                            il = 1;
                                        end
                                    end
                                    if il==1
                                        T0 = zeros(ns,1);
                                        for i=1:ns
                                            x = Xs(i);
                                            y = Ys(i);
                                            t = 0;
                                            T0(i) = eval(wzor_pocz);
                                        end
                                        for i=1:ilkT
                                            nr = WBTk(i,1);
                                            n1 = Kr(nr,1);
                                            n2 = Kr(nr,2);
                                            x = (Xs(n1) + Xs(n2))/2;
                                            y = (Ys(n1) + Ys(n2))/2;
                                            t = t0;
                                            WBTk(i,2) = eval(wzor_pocz);
                                        end
                                    end
                                end
                            end
                        end
                        if il==1
                            uiwait(msgbox('TEMPERATURA POCZ¥TKOWA ZOSTA£A PRZYPISANA','INFO','help','modal'));
                        end
                    end
                end
                if op1==3 | op1==2
                    pr = {['nr (1..' num2str(ilk) ') = ']};
                    odp = inputdlg(pr,'KRAWÊDZIE',1,{''});
                    if isempty(odp)==0
                        x = deal(odp{:});
                        blad = 0;
                        try
                            NUMER = eval(x);
                        catch
                            blad = 1;
                        end
                        if (NUMER>0 & NUMER<=ilk) & blad==0
                            nr = 0;
                            typ = 0;
                            for i=1:ilkT
                                if WBTk(i,1)==NUMER
                                    nr = i;
                                    typ = 1;
                                end
                            end
                            if op1==3
                                if nr>0
                                    for i=nr:ilkT-1
                                        WBTk(i,:) = WBTk(i+1,:);
                                    end
                                    WBTk(ilkT,:) = [];
                                    ilkT = ilkT-1;
                                end
                            end
                            if nr==0
                                for i=1:ilkQ
                                    if WBQk(i,1)==NUMER
                                        nr = i;
                                        typ = 2;
                                    end
                                end
                                if op1==3
                                    if nr>0
                                        for i=nr:ilkQ-1
                                            WBQk(i,:) = WBQk(i+1,:);
                                        end
                                        WBQk(ilkQ,:) = [];
                                        ilkQ = ilkQ-1;
                                    end
                                end
                            end
                            if op1==2
                                if typ>0
                                    if typ==1
                                        naz = 'TEMPERATURA T = ';
                                        dom = WBTk(nr,2);
                                    end
                                    if typ==2
                                        naz = 'STRUMIEÑ CIEP£A q = ';
                                        dom = WBQk(nr,2);
                                    end
                                    odp = inputdlg({naz},'WARUNKI BRZEGOWE - KOREKTA',1,{num2str(dom)});
                                    if isempty(odp)==0
                                        bl = 0;
                                        try
                                            [z] = deal(odp{:});
                                            z = eval(z);
                                        catch
                                            bl = 1;
                                        end
                                        if bl==0
                                            if typ==1
                                                WBTk(nr,2) = z;
                                            end
                                            if typ==2
                                                WBQk(nr,2) = z;
                                            end
                                        else
                                            uiwait(msgbox('B£¥D - Z£A WARTOŒÆ','B£¥D','error','modal'));
                                        end
                                    end
                                else
                                    uiwait(msgbox('B£¥D - TA KRAWÊD� NIE MA PRZYPISANEGO WARUNKU','B£¥D','error','modal'));
                                end
                            end
                        else
                            uiwait(msgbox('B£¥D - Z£A WARTOŒÆ','B£¥D','error','modal'));
                        end
                    end
                end
                if op1==1
                    pr = {['NUMER KRAWÊDZI (1..' num2str(ilk) ') = '],'TEMPERATURA T = ','STRUMIEÑ CIEP£A q = '};
                    odp = inputdlg(pr,'WARUNKI BRZEGOWE',1,{'',num2str(0),num2str(0)});
                    if isempty(odp)==0
                        bl = 0;
                        try
                            [z1 z2 z3] = deal(odp{:});
                            nr = eval(z1);
                            Te = eval(z2);
                            qe = eval(z3);
                        catch
                            bl = 1;
                        end
                        if bl==1 | (nr<=0 | nr>ilk)
                            uiwait(msgbox('B£¥D - Z£A WARTOŒÆ','B£¥D','error','modal'));
                        else
                            il = 0;
                            for i=1:ilkT
                                if nr==WBTk(i,1)
                                    il = 1;
                                    Tep = WBTk(i,2);
                                end
                            end
                            for i=1:ilkQ
                                if nr==WBQk(i,1)
                                    il = 2;
                                    qep = WBQk(i,2);
                                end
                            end
                            if Te~=0 & qe~=0
                                il = 3;
                            end
                            if il==0
                                if Te~=0 | qe==0
                                    ilkT = ilkT + 1;
                                    WBTk(ilkT,:) = [nr Te];
                                end
                                if qe~=0
                                    ilkQ = ilkQ + 1;
                                    WBQk(ilkQ,:) = [nr qe];
                                end
                            else
                                if il==1
                                    uiwait(msgbox(['B£¥D - DO TEJ KRAWÊDZI (nr = ' num2str(nr) ') JU¯ JEST PRZYPISANA TEMPERATURA T = ' num2str(Tep)],'B£¥D','error','modal'));
                                end
                                if il==2
                                    uiwait(msgbox(['B£¥D - DO TEJ KRAWÊDZI (nr = ' num2str(nr) ') JU¯ JEST PRZYPISANY STRUMIEÑ q = ' num2str(qep)],'B£¥D','error','modal'));
                                end
                                if il==3
                                    uiwait(msgbox('B£¥D - DO KRAWÊDZI MO¯E BYÆ PRZYPISANY TYLKO STRUMIEÑ LUB TYLKO TEMPERATURA','B£¥D','error','modal'));
                                end
                            end
                        end
                    end
                end
                if op1==4 | op1==6
                    ilkT = 0;
                    WBTk = [];
                end
                if op1==5 | op1==6
                    ilkQ = 0;
                    WBQk = [];
                end
            end
            close all;
        end
    end
    if op==5
        op1 = 1;
        ilop1 = 5;
        while op1<ilop1
            s11 = 'MES';
            if met==2
                s11 = 'BMRS';
            end
            s21 = 'PROSTOK¥TY';
            if ksztalt==1
                s21 = 'TRÓJK¥TY';
            end
            if met==1
                s31 = 'ELEMENTY BILINIOWE';
                if ksztalt==1 & stopien==1
                    s31 = 'ELEMENTY LINIOWE';
                end
                if stopien==2
                    s31 = 'ELEMENTY KWADRATOWE';
                    if ksztalt==2
                        s31 = 'ELEMENTY BIKWADRATOWE';
                    end
                end
            end
            if met==2
                s31 = 'LINIOWY';
                if stopien==2
                    s31 = 'KWADRATOWY';
                end
            end
            s41 = 'WARIACYJNE';
            if sfor==2
                s41 = 'LOKALNE';
            end
            s51 = num2str(liczgwiaz);
            s61 = 'USTALONY (t=const)';
            if analiza==2
                s61 = 'NIEUSTALONY';
            end
            if met==1
                s2 = ['KSZTA£T ELEMENTÓW : '  s21];
                s3 = ['RZ¥D INTERPOLACJI : '  s31];
            end
            if met==2
                s2 = ['CA£KOWANIE : '  s21];
                s3 = ['RZ¥D APROKSYMACJI MWLS : '  s31];
            end
            s4 = ['RODZAJ SFORMU£OWANIA ZADANIA :'  s41];
            s5 = ['LICZEBNOŒÆ GWIAZD RÓ¯NICOWYCH : '  s51];
            %op1 = menu('RODZAJ ANALIZY I PARAMETRY METODY',s1,s2,s3,...
            %   s4,s5,s6,'POWRÓT');
            s = str2mat('RODZAJ ANALIZY I PARAMETRY METODY');
            s = str2mat(s,'');
            s = str2mat(s,['METODA: ' num2str(s11)]);
            if met==1
                s = str2mat(s,s2);
                s = str2mat(s,s3);
            end
            if met==2
                s = str2mat(s,s2);
                s = str2mat(s,s3);
                s = str2mat(s,s4);
                s = str2mat(s,s5);
            end
            s = str2mat(s,'');
            s = str2mat(s,['PRZEP£YW: ' num2str(s61)]);
            if analiza==2
                sp1 = 'JAWNY';
                if schemat==2
                    sp1 = 'NIEJAWNY';
                end
                s = str2mat(s,['SCHEMAT CA£KOWANIA PO CZASIE: ' sp1]);
                s = str2mat(s,['KROK CZASOWY dt = ' num2str(dt)]);
                s = str2mat(s,['LICZBA KROKÓW CZASOWYCH nt = ' num2str(nt)]);
                sp2 = 'ROZWI¥ZANIE ZADANIA PRZEP£YWU USTALONEGO';
                if sp_pocz==2
                    sp2 = ['WG DANEJ ZALE¯NOŒCI T(x,y,t) = ' wzor_pocz];
                end
                sp3 = sprintf('�RÓD£O TEMPERATURY POCZ¥TKOWEJ: \n %s',sp2);
                s = str2mat(s,sp3);
            end
            op1 = menu(s,'WYBÓR METODY KOMPUTEROWEJ',['PARAMETRY ' s11],'CHARAKTER PRZEP£YWU','PARAMETRY PRZEP£YWU NIEUSTALONEGO','WYJŒCIE');
            if op1==4
                if analiza==1
                    uiwait(msgbox('PARAMETRY DOTYCZ¥ TYLKO PRZEP£YWU NIEUSTALONEGO','B£¥D','error','non-modal'));
                end
                if analiza==2
                    s1 = 'TEMPERATURA DLA PRZEP£YWU USTALONEGO';
                    s2 = 'ZADANIE FUNKCJI TEMPERATURY';
                    s = s2;
                    if sp_pocz == 1
                        s = s1;
                    end
                    dec = questdlg('PROSZÊ WYBRAÆ SPOSÓB GENERACJI TEMPERATURY POCZ¥TKOWEJ','TEMPERATURA POCZ¥TKOWA',s1,s2,'DALEJ',s);
                    il = 0;
                    if strcmp(dec,s1)==1
                        il = 1;
                        sp_pocz = 1;
                    end
                    if strcmp(dec,s2)==1
                        il = 1;
                        sp_pocz = 2;
                    end
                    if il==1
                        uiwait(msgbox(['USTAWIONO WARUNEK POCZ¥TKOWY: ' dec],'INFO','help','non-modal'));
                    end
                    pr = {'KROK CZASOWY dt = ','LICZBA KROKÓW CZASOWYCH n_t = ','CZAS KOÑCOWY t_k = '};
                    odp = inputdlg(pr,'CA£KOWANIE PO CZASIE',1,{num2str(dt),num2str(nt),num2str(tk)});
                    if isempty(odp)==0
                        bl = 0;
                        try
                            [z1 z2 z3] = deal(odp{:});
                            dt0 = eval(z1);
                            nt0 = eval(z2);
                            tk0 = eval(z3);
                        catch
                            bl = 1;
                        end
                        if bl==1 | (dt0<0 | nt0<0 | nt0/2~=round(nt0/2) | tk0<0) | (dt0==0 & nt0==0) | (nt0==0 & tk0==0) | (dt0==0 & tk0==0)
                            uiwait(msgbox('B£¥D - Z£A WARTOŒÆ','B£¥D','error','modal'));
                        else
                            if nt0~=0 & dt0~=0 & tk0~=0
                                tk0 = 0;
                            end
                            if nt0==0
                                nt = tk0/dt0;
                                tk = tk0;
                                dt = dt0;
                            end
                            if dt0==0
                                dt = tk0/nt0;
                                tk = tk0;
                                nt = nt0;
                            end
                            if tk0==0
                                tk0 = nt0*dt0;
                                dt = dt0;
                                nt = nt0;
                            end
                            uiwait(msgbox(['USTAWIONO KROK CZASOWY dt = ' num2str(dt) ' oraz LICZBÊ KROKÓW nt = ' num2str(nt)],'INFO','help','non-modal'));
                            s1 = 'SCHEMAT JAWNY (WARUNKOWO STABILNY)';
                            s2 = 'SCHEMAT NIEJAWNY (BEZWARUNKOWO STABILNY)';
                            s = s2;
                            if schemat == 1
                                s = s1;
                            end
                            dec = questdlg('PROSZÊ WYBRAÆ SCHEMAT CA£KOWANIA PO CZASIE','CA£KOWANIE PO CZASIE',s1,s2,'DALEJ',s);
                            il = 0;
                            if strcmp(dec,s1)==1
                                il = 1;
                                schemat = 1;
                            end
                            if strcmp(dec,s2)==1
                                il = 1;
                                schemat = 2;
                            end
                            if il==1
                                uiwait(msgbox(['USTAWIONO SCHEMAT CA£KOWANIA PO CZASIE: ' dec],'INFO','help','non-modal'));
                            end
                        end
                    end
                end
            end
            if op1==2
                s1 = 'TRÓJK¥TY';
                s2 = 'PROSTOK¥TY';
                s3 = 'DALEJ';
                s = s1;
                if ksztalt==2
                    s = s2;
                end
                if ksztalt==3
                    s = s3;
                end
                sp = ' CA£KOWANIA';
                if met==1
                    sp = '';
                end
                dec = questdlg(['PROSZÊ WYBRAÆ KSZTA£T ELEMENTÓW' sp], ...
                    ['KSZTA£T ELEMENTÓW' sp], ...
                    s1,s2,s3,s);
                if Nes>0
                    if strcmp(dec,s1)==1 | strcmp(dec,s2)==1
                        sp1 = 'SIATKI MES';
                        if met==2
                            sp1 = 'WÊZ£ÓW BMRS';
                        end
                        uiwait(msgbox(['TA OPERACJA JEST MO¯LIWA PO UPRZEDNIM USUNIÊCIU ' sp1],'B£¥D','error','non-modal'));
                    end
                else
                    il = 0;
                    if strcmp(dec,s1)==1
                        lwe = 3;
                        il = 1;
                        ksztalt=1;
                    end
                    if strcmp(dec,s2)==1
                        lwe = 4;
                        il = 1;
                        ksztalt=2;
                    end
                    if il==1
                        s = '';
                        if met==1
                            s = [' (LWE = ' num2str(lwe) ')'];
                        end
                        uiwait(msgbox(['WYBRANO KSZTA£T ELEMENTÓW ' sp ' = ' dec s],['KSZTA£T ELEMENTÓW ' sp],'help','non-modal'));
                    end
                end
                s1 = 'LINIOWY';
                s2 = 'KWADRATOWY';
                if ksztalt==2 & met==1
                    s1 = 'BILINIOWY';
                    s2 = 'BIKWADRATOWY';
                end
                s = s1;
                if stopien==2
                    s = s2;
                end
                sp = s11;
                dec = questdlg(['PROSZÊ WYBRAÆ RZ¥D APROKSYMACJI dla ' sp], ...
                    ['RZ¥D APROKSYMACJI dla ' sp], ...
                    s1,s2,'DALEJ',s);
                il = 0;
                if strcmp(dec,s1)==1
                    il = 1;
                    stopien=1;
                end
                if strcmp(dec,s2)==1
                    stopien=2;
                    il = 1;
                end
                if il==1
                    uiwait(msgbox(['WYBRANO RZ¥D APROKSYMACJI dla ' sp ' = ' dec],['RZ¥D APROKSYMACJI dla ' sp],'help','non-modal'));
                end
                if met==2
                    s1 = 'WARIACYJNE';
                    s2 = 'LOKALNE';
                    s = s1;
                    if sfor==2
                        s = s2;
                    end
                    dec = questdlg(['PROSZÊ WYBRAÆ SFORMU£OWANIE PROBLEMU BRZEGOWEGO dla ' sp], ...
                        ['SFORMU£OWANIE PROBLEMU BRZEGOWEGO dla ' sp], ...
                        s1,s2,'DALEJ',s);
                    il = 0;
                    if strcmp(dec,s1)==1
                        sfor=1;
                        il = 1;
                    end
                    if strcmp(dec,s2)==1
                        sfor=2;
                        il = 1;
                    end
                    if il==1
                        uiwait(msgbox(['WYBRANO SFORMU£OWANIE PROBLEMU BRZEGOWEGO dla ' sp ' = ' dec],['SFORMU£OWANIE PROBLEMU BRZEGOWEGO dla ' sp],'help','non-modal'));
                    end
                    LICZ = [9 16 24];
                    s1 = num2str(LICZ(1));
                    s2 = num2str(LICZ(2));
                    s3 = num2str(LICZ(3));
                    s = num2str(liczgwiaz);
                    dec = questdlg(['PROSZÊ WYBRAÆ LICZEBNOŒÆ GWIAZD RÓ¯NICOWYCH dla ' sp], ...
                        ['LICZEBNOŒÆ GWIAZD RÓ¯NICOWYCH dla ' sp], ...
                        s1,s2,s3,s);
                    if strcmp(dec,s1)==1
                        liczgwiaz = LICZ(1);
                    end
                    if strcmp(dec,s2)==1
                        liczgwiaz = LICZ(2);
                    end
                    if strcmp(dec,s3)==1
                        liczgwiaz = LICZ(3);
                    end
                    uiwait(msgbox(['WYBRANO LICZEBNOŒÆ GWIAZD RÓ¯NICOWYCH dla ' sp ' = ' dec],['LICZEBNOŒÆ GWIAZD RÓ¯NICOWYCH dla ' sp],'help','non-modal'));
                end
            end
            if op1==1
                s1 = 'METODA ELEMENTÓW SKOÑCZONYCH (MES)';
                s2 = 'BEZSIATKOWA METODA RÓ¯NIC SKOÑCZONYCH (BMRS)';
                s = s1;
                if met==2
                    s = s2;
                end
                dec = questdlg('PROSZÊ WYBRAÆ METODÊ OBLICZENIOW¥', ...
                    'METODA OBLICZENIOWA', ...
                    s1,s2,'POWRÓT',s);
                if strcmp(dec,s1)==1
                    met = 1;
                    sfor = 1;
                end
                if strcmp(dec,s2)==1
                    met = 2;
                end
                uiwait(msgbox(['WYBRANO METODÊ OBLICZENIOW¥ : ' dec],'METODA OBLICZENIOWA','help','non-modal'));
            end
            if op1==3
                s1 = 'USTALONY (t=const)';
                s2 = 'NIEUSTALONY';
                s = s1;
                if analiza==2
                    s = s2;
                end
                dec = questdlg('PROSZÊ WYBRAÆ CHARAKTER PRZEP£YWU', ...
                    'PRZEP£YW', ...
                    s1,s2,'POWRÓT',s);
                if strcmp(dec,s1)==1
                    analiza = 1;
                end
                if strcmp(dec,s2)==1
                    analiza = 2;
                end
            end
        end
    end
    if op==6
        if ilk<4
            uiwait(msgbox('B£¥D - ZA MA£O KRAWÊDZI BRZEGOWYCH','B£¥D','error','modal'));
        elseif met~=1
            uiwait(msgbox('B£¥D - EDYCJI PODLEGAJ¥ TYLKO ELEMENTY SKOÑCZONE','B£¥D','error','modal'));
        elseif Nes>50
            uiwait(msgbox('B£¥D - ZA DU¯O ELEMENTÓW DO EDYCJI','B£¥D','error','modal'));
        else
            op1 = 1;
            ilop1 = 8;
            if (isempty(Xs)==1 | Nes==0) %| siat==2 % | siat==0
                Xs = Xk;
                Ys = Yk;
                ns = nk;
                nspop = 0;
            end
            siat = 1;
            while op1<ilop1
                rysunek(1,nk,Xk,Yk,ilk,Kr,2,WBTk,WBQk,Nes,Incs,Xs,Ys,lwe);
                if op1==6
                    s = 'LISTA WÊZ£ÓW';
                    for i=1:ns
                        s = str2mat(s,['i = ' num2str(i) ', x = ' num2str(Xs(i)) ', y = ' num2str(Ys(i))]);
                    end
                    s = str2mat(s,' ');
                    s = str2mat(s,'LISTA ELEMENTÓW');
                    for i=1:Nes
                        s = str2mat(s,['i = ' num2str(i) ', Inc = ' num2str(Incs(i,:))]);
                    end
                    uiwait(msgbox(s,'SIATKA MES','help','non-modal'));
                end
                op1 = menu('RÊCZNA GENERACJA SIATKI','NOWY ELEMENT','KOREKTA ELEMENTU','USUNIÊCIE ELEMENTU',...
                    'USUNIÊCIE WSZYSTKICH ELEMENTÓW','USUNIÊCIE SIATKI MES','LISTA WÊZ£ÓW I ELEMENTÓW','PRZEKSZTA£CENIE SIATKI PROSTOK¥TNEJ NA TRÓJK¥TN¥','POWRÓT');
                NUMER = Nes+1;
                if op1==7
                    if ksztalt==1
                        uiwait(msgbox('ZMIANA MO¯LIWA JEST JEDYNIE Z SIATKI PROSTOK¥TNEJ NA TRÓJK¥TN¥','B£¥D','error','non-modal'));
                    end
                    if ksztalt==2
                        if Nes==0
                            uiwait(msgbox('BRAK SIATKI DO PRZEKSZTA£CENIA','B£¥D','error','non-modal'));
                        else
                            Nes1 = 0;
                            Incs1 = zeros(2*Nes,3);
                            for i=1:Nes
                                Nes1 = Nes1 + 1;
                                Incs1(Nes1,:) = Incs(i,[1 2 3]);
                                Nes1 = Nes1 + 1;
                                Incs1(Nes1,:) = Incs(i,[1 3 4]);
                            end
                            Incs = Incs1;
                            Nes = Nes1;
                            lwe = 3;
                            ksztalt = 1;
                        end
                    end
                end
                if op1==2 | op1==3
                    pr = {['nr (1..' num2str(Nes) ') = ']};
                    odp = inputdlg(pr,'ELEMENTY',1,{''});
                    if isempty(odp)==0
                        x = deal(odp{:});
                        blad = 0;
                        try
                            NUMER = eval(x);
                        catch
                            blad = 1;
                        end
                        if (NUMER>0 & NUMER<=Nes) & blad==0
                            if op1==3
                                for j=NUMER:Nes-1
                                    Incs(j,:) = Incs(j+1,:);
                                end
                                Nes = Nes-1;
                            end
                            if op1==2
                                op1 = 1;
                            end
                        else
                            uiwait(msgbox('B£¥D - Z£A WARTOŒÆ','B£¥D','error','modal'));
                        end
                    end
                end
                if op1==1
                    naz = ['ELEMENT nr = ' num2str(NUMER)];
                    if ksztalt==2
                        pr = {[naz ' : wêze³ 1 = ' ],[naz ' : wêze³ 2 = ' ],[naz ' : wêze³ 3 = ' ],[naz ' : wêze³ 4 = ' ]};
                    end
                    if ksztalt==1
                        pr = {[naz ' : wêze³ 1 = ' ],[naz ' : wêze³ 2 = ' ],[naz ' : wêze³ 3 = ' ]};
                    end
                    if NUMER==Nes+1
                        nr1 = '';
                        nr2 = '';
                        nr3 = '';
                        if ksztalt==2
                            nr4 = '';
                        end
                    else
                        nr1 = num2str(Incs(NUMER,1));
                        nr2 = num2str(Incs(NUMER,2));
                        nr3 = num2str(Incs(NUMER,3));
                        if ksztalt==2
                            nr4 = num2str(Incs(NUMER,4));
                        end
                    end
                    if ksztalt==2
                        s={nr1,nr2,nr3,nr4};
                    end
                    if ksztalt==1
                        s={nr1,nr2,nr3};
                    end
                    odp = inputdlg(pr,'ELEMENT SKOÑCZONY',1,s);
                    %                     nk
                    %                     Xk
                    %                     Yk
                    %                     ns
                    %                     Ys
                    %                     Xs
                    if isempty(odp)==0
                        bl = 0;
                        try
                            if ksztalt==2
                                [z1 z2 z3 z4] = deal(odp{:});
                                nr = [eval(z1) eval(z2) eval(z3) eval(z4)];
                            end
                            if ksztalt==1
                                [z1 z2 z3] = deal(odp{:});
                                nr = [eval(z1) eval(z2) eval(z3)];
                            end
                        catch
                            bl = 1;
                        end
                        if bl==0
                            for i=1:lwe
                                if nr(i)<1 | nr(i)>ns
                                    bl = 1;
                                end
                            end
                        end
                        if bl==1
                            uiwait(msgbox('B£¥D - Z£A WARTOŒÆ','B£¥D','error','modal'));
                        else
                            if NUMER==Nes+1
                                Nes = Nes + 1;
                            end
                            Incs(NUMER,:) = nr;
                        end
                    end
                end
                if op1==4 | op1==5
                    Nes = 0;
                    Incs = [];
                    WBs = [];
                end
                if op1==5
                    ns = nk;
                    nspop = 0;
                    Xs = Xk;
                    Ys = Yk;
                end
            end
            WBs = zeros(ns,3);
            for i=1:ns
                [WBs(i,1) nrkr] = pbrzeg(Xs(i),Ys(i),ilk,Kr,Xk,Yk);
                if WBs(i,1)==1
                    for j=1:length(nrkr)
                        WBs(i,1+j) = nrkr(j);
                    end
                end
            end
            close all;
            if Nes==0
                siat = 0;
            end
        end
    end
    if op==7
        if ilk<4
            uiwait(msgbox('B£¥D - ZA MA£O KRAWÊDZI BRZEGOWYCH','B£¥D','error','modal'));
        else
            op1 = 1;
            ilop1 = 6;
            while op1<ilop1
                op1 = menu('GENERATOR SIATKI MES / WÊZ£ÓW BMRS','LICZBA PODZIA£ÓW KRAWÊDZI','GENERACJA SIATKI MES / WÊZ£ÓW BMRS',...
                    'USUNIÊCIE SIATKI MES / WÊZ£ÓW BMRS','PRZEKSZTA£CENIE PROSTOK¥TÓW NA TRÓJK¥TY','PARAMET ZABURZENIA WÊZ£ÓW WEWNÊTRZNYCH','POWRÓT');
                if op1==4
                    if ksztalt==1
                        uiwait(msgbox('ZMIANA MO¯LIWA JEST JEDYNIE Z SIATKI (CA£KOWANIA) PROSTOK¥TNEJ NA TRÓJK¥TN¥','B£¥D','error','non-modal'));
                    end
                    if ksztalt==2
                        if Nes==0
                            uiwait(msgbox('BRAK ELEMENTÓW DO PRZEKSZTA£CENIA','B£¥D','error','non-modal'));
                        else
                            Nes1 = 0;
                            Incs1 = zeros(2*Nes,3);
                            for i=1:Nes
                                Nes1 = Nes1 + 1;
                                Incs1(Nes1,:) = Incs(i,[1 2 3]);
                                Nes1 = Nes1 + 1;
                                Incs1(Nes1,:) = Incs(i,[1 3 4]);
                            end
                            Incs = Incs1;
                            Nes = Nes1;
                            lwe = 3;
                            ksztalt = 1;
                            uiwait(msgbox('SIATKA (CA£KOWANIA) ZOSTA£A PRZEKSZTA£CONA NA SIATKÊ TRÓJK¥TN¥','INFORMACJA','help','non-modal'));
                        end
                    end
                end
                if op1==5
                    pr = {'AMPLITUDA ZABURZENIA = '};
                    odp = inputdlg(pr,'LOSOWE ZABURZENIE WÊZ£ÓW WEWNÊTRZNYCH',1,{num2str(amp_zab)});
                    bl = 0;
                    try
                        z = deal(odp{:});
                        nr = eval(z);
                    catch
                        bl = 1;
                    end
                    if nr<0
                        bl = 1;
                    end
                    if bl==1
                        uiwait(msgbox('B£¥D - Z£A WARTOŒÆ','B£¥D','error','modal'));
                    else
                        amp_zab = nr;
                    end
                end
                if op1==1
                    pr = {'LICZBA PODZIA£ÓW = '};
                    odp = inputdlg(pr,'GÊSTOŒÆ WÊZ£ÓW',1,{num2str(nh)});
                    bl = 0;
                    try
                        z = deal(odp{:});
                        nr = eval(z);
                    catch
                        bl = 1;
                    end
                    if nr<=0
                        bl = 1;
                    end
                    if bl==1
                        uiwait(msgbox('B£¥D - Z£A WARTOŒÆ','B£¥D','error','modal'));
                    else
                        nh = nr;
                    end
                end
                if op1==2
                    [ns_test,Xs_test,Ys_test,Nes_test,Incs_test,WBs_test] = siatka(nh,ilk,Kr,Xk,Yk,ksztalt);
                    if ns_test>0
                        siat = 2;
                        ns = ns_test;
                        Xs = Xs_test;
                        Ys = Ys_test;
                        Nes = Nes_test;
                        Incs = Incs_test;
                        WBs = WBs_test;
                        dlmax = 0;
                        for i=1:ilk
                            n1 = Kr(i,1);
                            n2 = Kr(i,2);
                            x1 = Xk(n1);
                            x2 = Xk(n2);
                            y1 = Yk(n1);
                            y2 = Yk(n2);
                            dl = sqrt( (x1-x2)^2 + (y1-y2)^2 );
                            if dl>dlmax
                                dlmax = dl;
                            end
                        end
                        for i=1:ns
                            if WBs(i,1)==0
                                Xs(i) = Xs(i) - (0.5 - rand(1,1))*amp_zab*Xs(i)/dlmax;
                                Ys(i) = Ys(i) - (0.5 - rand(1,1))*amp_zab*Ys(i)/dlmax;
                            end
                        end
                        uiwait(msgbox('SIATKA MES / WÊZ£Y BMRS ZOSTA£Y WYGENEROWANE','INFO','help','modal'));
                    else
                        uiwait(msgbox('PRÓBA GENERACJI ZBYT GÊSTEGO ZBIORU WÊZ£ÓW! PROSZÊ ZMNIEJSZYÆ JEGO GÊSTOŒÆ','B£¥D','error','modal'));
                    end
                end
                if op1==3
                    ns = nk;
                    nspop = 0;
                    Xs = Xk;
                    Ys = Yk;
                    Incs = [];
                    Nes = 0;
                    WBs = [];
                    uiwait(msgbox('SIATKA MES / WÊZ£Y BMRS ZOSTA£Y USUNIÊTE','INFO','help','modal'));
                    WBs = zeros(ns,3);
                    for i=1:ns
                        [WBs(i,1) nrkr] = pbrzeg(Xs(i),Ys(i),ilk,Kr,Xk,Yk);
                        if WBs(i,1)==1
                            for j=1:length(nrkr)
                                WBs(i,1+j) = nrkr(j);
                            end
                        end
                    end
                end
            end
        end
    end
    if op==5 | op==6 | op==7 | op==10
        lsse = lwe;
        if stopien>1 & met==1
            lsse = stopien*lwe;
        end
        if stopien==1 & Nes>0
            Incse = Incs;
            if nspop>0
                ns = nspop;
                nspop = 0;
            end
        end
        if stopien>1 & Nes>0 & met==1
            Incs1 = zeros(Nes,lsse);
            Xs1 = zeros(ns+4*ns,1);
            Ys1 = zeros(ns+4*ns,1);
            WBs1 = zeros(ns+4*ns,3);
            Xs1(1:ns) = Xs;
            Ys1(1:ns) = Ys;
            WBs1(1:ns,:) = WBs;
            ns1 = ns;
            w = zeros(lwe+1,1);
            for i=1:Nes
                w(1:lwe) = Incs(i,:);
                w(lwe+1) = Incs(i,1);
                for j=1:lwe
                    x = 0.5*(Xs(w(j))+Xs(w(j+1)));
                    y = 0.5*(Ys(w(j))+Ys(w(j+1)));
                    nr = 0;
                    for k=1:ns1
                        d = sqrt( (x-Xs1(k))^2 + (y-Ys1(k))^2 );
                        if d<eps
                            nr = k;
                        end
                    end
                    if nr==0
                        ns1 = ns1+1;
                        nr = ns1;
                        Xs1(ns1) = x;
                        Ys1(ns1) = y;
                        [WBs1(ns1,1) nrkr] = pbrzeg(x,y,ilk,Kr,Xk,Yk);
                        if WBs1(ns1,1)==1
                            for k=1:length(nrkr)
                                WBs1(ns1,1+k) = nrkr(k);
                            end
                        end
                    end
                    Incs1(i,2*j-1) = w(j);
                    Incs1(i,2*j) = nr;
                end
            end
            Xs = Xs1(1:ns1);
            Ys = Ys1(1:ns1);
            if nspop==0
                nspop = ns;
            end
            ns = ns1;
            WBs = WBs1(1:ns1,:);
            Incse = Incs1;
        end
    end
    if op==10
        il = 0;
        s = 'B£ÊDY W MODELU - SZCZEGÓ£Y PONI¯EJ';
        if h == 0
            s = str2mat(s,[' ' '- ZEROWA GRUBOŒÆ TARCZY h']);
            il = il + 1;
        end
        if kx == 0
            s = str2mat(s,[' ' '- ZEROWY WSPÓ£CZYNNIK PRZEWODNOŒCI CIEPLNEJ kx']);
            il = il + 1;
        end
        if ky == 0
            s = str2mat(s,[' ' '- ZEROWY WSPÓ£CZYNNIK PRZEWODNOŒCI CIEPLNEJ ky']);
            il = il + 1;
        end
        if analiza==2
            if schemat==1
                hmin = 10^10;
                for i=1:Nes
                    w = Incs(i,:);
                    w(lwe+1) = w(1);
                    for j=1:lwe
                        hx = abs(Xs(w(j+1)) - Xs(w(j)));
                        hy = abs(Ys(w(j+1)) - Ys(w(j)));                        
                        if hx<hmin & hx~=0
                            hmin = hx;
                        end
                        if hy<hmin & hy~=0
                            hmin = hy;
                        end
                    end
                end
                hmin
                wsp = hmin^2/max([kx ky]);                
                wsp
                if dt>=wsp
                    dts = wsp*0.25;
                    s1 = ['OBLICZENIA DLA NIESTABILNEGO dt = ' num2str(dt)];
                    s2 = ['ZMIANA KROKU dt = ' num2str(dts)];
                    s3 = ['ZMIANA NA SCHEMAT NIEJAWNY dla dt = ' num2str(dt)];
                    dec = questdlg('NIESTABILNY SCHEMAT CA£KOWANIA PO CZASIE','NIESTABILNY SCHEMAT CA£KOWANIA PO CZASIE',s1,s2,s3,s2);
                    if strcmp(dec,s2)
                        dt = dts;
                    end
                    if strcmp(dec,s3);
                        schemat = 2;
                    end
                end
            end
            if sp_pocz==2
                T0 = zeros(ns,1);
                for i=1:ns
                    x = Xs(i);
                    y = Ys(i);
                    t = t0;
                    T0(i) = eval(wzor_pocz);
                end
            end
            if sp_pocz==1
                if isempty(T)==0 & length(T)==ns
                    T0 = T;
                else
                    %s = str2mat(s,[' ' '- B£ÊDNA TEMPERATURA POCZ¥TKOWA LUB JEJ BRAK']);
                    %il = il + 1;
                    analiza = 1;
                    [T,R] = solver(ns,Xs,Ys,WBs,Incs,Incse,Nes,lwe,lsse,ksztalt,stopien,met,sfor,liczgwiaz,ilkT,WBTk,ilkQ,WBQk,kx,ky,f);
                    analiza = 2;
                    T0 = T;
                end
            end
        end
        if il>0
            uiwait(msgbox(s,'B£¥D','error','modal'));
        else
            if siat==0
                uiwait(msgbox('B£¥D - BRAK SIATKI MES / WÊZ£ÓW BMRS','B£¥D','error','modal'));
            else
                if analiza==1
                    mes = msgbox('ZADANIE JEST ROZWI¥ZYWANE, PROSZÊ CZEKAÆ...','INFORMACJA','help','replace');
                end
                if analiza==1
                    czas = t0;
                    [T,R] = solver(ns,Xs,Ys,WBs,Incs,Incse,Nes,lwe,lsse,ksztalt,stopien,met,sfor,liczgwiaz,ilkT,WBTk,ilkQ,WBQk,kx,ky,f);
                    [QG,PG] = postprocess(ns,Xs,Ys,Incs,Incse,Nes,lwe,lsse,ksztalt,met,sfor,stopien,liczgwiaz,kx,ky,T);
                end
                if analiza==2
                    %                     dt = 0.1;
                    %                     nt = 20;
                    %                     t0 = 0;
                    %                     T0 = 10;
                    %                     Tp = ones(ns,1)*T0;
                    Tp = T0;
                    Tnu = zeros(nt,ns);
                    for t=1:nt
                        czas = t0 + t*dt;
                        mes = msgbox(['ZADANIE JEST ROZWI¥ZYWANE, PROSZÊ CZEKAÆ... (KROK #' num2str(t) '/' num2str(nt) ' , t = ' num2str(czas) ')'],'INFORMACJA','help','replace');
                        [Tn,R] = solver(ns,Xs,Ys,WBs,Incs,Incse,Nes,lwe,lsse,ksztalt,stopien,met,sfor,liczgwiaz,ilkT,WBTk,ilkQ,WBQk,kx,ky,f);                        
                        if analiza==2 & schemat==1
                            Tn = Tp + dt*Tn;
                        end
                        if t==nt & met==1
                            [QG,PG] = postprocess(ns,Xs,Ys,Incs,Incse,Nes,lwe,lsse,ksztalt,met,sfor,stopien,liczgwiaz,kx,ky,Tn);
                        end
                        if t==1
                            Rmax = max(R);
                            Rmin = min(R);
                        end
                        if t>1
                            rmax = max(R);
                            if rmax>Rmax
                                Rmax = rmax;
                            end
                            rmin = min(R);
                            if rmin<Rmin
                                Rmin = rmin;
                            end
                        end
                        Tnu(t,:) = Tn';
                        Tp = Tn;
                        try
                            close(mes);
                        catch
                        end
                    end
                end
                try
                    close(mes);
                catch
                end
                close all;
                uiwait(msgbox('ZADANIE ZOSTA£O ROZWI¥ZANE','INFO','help','modal'));
                obl = 1;
            end
        end
    end
    if op==9
        rysunek(2,nk,Xk,Yk,ilk,Kr,2,WBTk,WBQk,Nes,Incs,Xs,Ys,lwe);
        for i=1:ns
            p = plot(Xs(i),Ys(i),'ro');
            set(p,'LineWidth',2);
            % if siat==2
            if WBs(i,1)==1
                set(p,'MarkerFaceColor','r');
            end
            % end
            %t = text(Xs(i)+0.01*(xmax+xmin),Ys(i)+0.01*(ymax+ymin),num2str(i));
            %set(t,'FontSize',8);
            %set(t,'Color','r');
        end
        fig = gcf;
        set(fig,'Position',get(0,'ScreenSize'));
    end
    if op==8
        rysunek(1,nk,Xk,Yk,ilk,Kr,1,WBTk,WBQk,Nes,Incs,Xs,Ys,lwe);
        fig = gcf;
        set(fig,'Position',get(0,'ScreenSize'));
    end
    if op==11
        if obl==1
            warwez = 0;
            xmax = max(Xk);
            xmin = min(Xk);
            ymax = max(Yk);
            ymin = min(Yk);
            ilop1 = 7;
            op1 = 1;
            while op1<ilop1
                s1 = 'W£¥CZONE';
                if wygladz == 0
                    s1 = 'WY£¥CZONE';
                end
                s2 = 'W£¥CZONE';
                if warwez == 0
                    s2 = 'WY£¥CZONE';
                end
                op1 = menu('WYNIKI - TEKSTOWO I GRAFICZNIE','ROZK£AD TEMPERATURY','SK£ADOWE STRUMIENIA CIEP£A','WYNIKI WÊZ£OWE - LISTA',...
                    'WYNIKI ELEMENTOWE - LISTA',['WYG£ADZANIE WYKRESÓW: ' s1],['WARTOŒCI WÊZ£OWE NA WYKRESACH: ' s2],'POWRÓT');
                if op1==5
                    if wygladz==1
                        wygladz=0;
                    else
                        wygladz=1;
                    end
                end
                if op1==6
                    if warwez==1
                        warwez=0;
                    else
                        warwez=1;
                    end
                end
                if op1==3
                    try
                        close(msg);
                    catch
                    end
                    s = 'TEMPERATURA I STRUMIEÑ CIEP£A WE WÊZ£ACH';
                    for i=1:ns
                        s = str2mat(s,['i = ' sprintf('%4d',i) ', x = ' sprintf('%-8.4f',Xs(i)) ', y = ' sprintf('%-8.4f',Ys(i)) ', T  = ' sprintf('%-10.4f',T(i)) ', q = ' sprintf('%-10.4f',R(i))]);
                    end
                    msg = msgbox(s,'WYNIKI WÊZ£OWE','help','modal');
                end
                if op1==4
                    try
                        close(msg);
                    catch
                    end
                    s = 'STRUMIEÑ CIEP£A W ELEMENTACH (WÊZ£Y I ŒRODEK)';
                    for i=1:Nes
                        %                         t1 = '';
                        %                         t2 = '';
                        %                         for j=1:5
                        %                             t1 = strcat(t1,['  q_x(' num2str(PG(i,j,1),2) ',' num2str(PG(i,j,2),2) ') = ' num2str(QG(i,1,j),4) ' ; ']);
                        %                             t2 = strcat(t2,['  q_y(' num2str(PG(i,j,1),2) ',' num2str(PG(i,j,2),2) ') = ' num2str(QG(i,2,j),4) ' ; ']);
                        %                         end
                        if ksztalt==2
                            t1 = ['(' sprintf('%10.5f',QG(i,1,4)) ' ; ' sprintf('%10.5f',QG(i,2,4)) ')' '                   ' '(' sprintf('%10.5f',QG(i,1,3)) ' ; ' sprintf('%10.5f',QG(i,2,3)) ')'];
                            t2 = ['                            ' '(' sprintf('%10.5f',QG(i,1,5)) ' ; ' sprintf('%10.5f',QG(i,2,5)) ')'];
                            t3 = ['(' sprintf('%10.5f',QG(i,1,1)) ' ; ' sprintf('%10.5f',QG(i,2,1)) ')' '                   ' '(' sprintf('%10.5f',QG(i,1,2)) ' ; ' sprintf('%10.5f',QG(i,2,2)) ')'];
                        end
                        if ksztalt==1
                            t1 = ['                            ' '(' sprintf('%10.5f',QG(i,1,1)) ' ; ' sprintf('%10.5f',QG(i,2,1)) ')'];
                            t2 = ['                            ' '(' sprintf('%10.5f',QG(i,1,4)) ' ; ' sprintf('%10.5f',QG(i,2,4)) ')'];
                            t3 = ['(' sprintf('%10.5f',QG(i,1,2)) ' ; ' sprintf('%10.5f',QG(i,2,2)) ')' '                   ' '(' sprintf('%10.5f',QG(i,1,3)) ' ; ' sprintf('%10.5f',QG(i,2,3)) ')'];
                        end
                        s = str2mat(s,['i = ' num2str(i)]);
                        s = str2mat(s,t1);
                        s = str2mat(s,t2);
                        s = str2mat(s,t3);
                        s = str2mat(s,' ');
                    end
                    msg = msgbox(s,'WYNIKI ELEMENTOWE','help','modal');
                    set(msg,'Resize','on')
                end
                if op1==1
                    fig = figure(3);
                    pocz = 1;
                    N = 1;
                    if analiza==2
                        N = nt;
                        s1 = 'ROZK£AD TEMPERATURY W CZASIE CA£EGO PROCESU';
                        s2 = 'KOÑCOWY ROZK£AD TEMPERATURY';
                        s = s1;
                        dec = questdlg('PREZENTACJA WYNIKÓW TEMPERATURY','TEMPERATURA',s1,s2,'POWRÓT',s);
                        il = 0;
                        if strcmp(dec,s1)==1
                            il = 1;
                            pocz = 1;
                        end
                        if strcmp(dec,s2)==1
                            il = 1;
                            pocz = nt;
                        end
                        if il==0
                            pocz = N+1;
                            close all;
                        end
                        Tmax = max(max(Tnu));
                        Tmin = min(min(Tnu));
                    end
                    for k=pocz:N
                        if analiza==1
                            Temp = T;
                        end
                        if analiza==2
                            Temp = Tnu(k,:)';
                        end
                        hold off
                        plot3(0,0,0,'w');
                        view(2);
                        hold on
                        grid on
                        set(fig,'Position',get(0,'ScreenSize'));
                        x = zeros(lwe,1);
                        y = zeros(lwe,1);
                        c = zeros(lwe,1);
                        for i=1:Nes
                            w = zeros(lwe,1);
                            w(1:lwe) = Incs(i,:);
                            %w(lwe+1) = Incs(i,1);                                
                            for j=1:lwe
                                x(j)=Xs(w(j));
                                y(j)=Ys(w(j));
                            end
                            if ksztalt==10                                                          
                                x_min = min(x);
                                x_max = max(x);
                                y_min = min(y);
                                y_max = max(y);
                                nr1 = 0;
                                nr2 = 0;
                                nr3 = 0;
                                nr4 = 0;
                                for j=1:lwe
                                    if x(j)==x_min & y(j)==y_min
                                        nr1 = w(j);
                                    end
                                    if x(j)==x_max & y(j)==y_min
                                        nr2 = w(j);
                                    end
                                    if x(j)==x_max & y(j)==y_max
                                        nr3 = w(j);
                                    end
                                    if x(j)==x_min & y(j)==y_max
                                        nr4 = w(j);
                                    end
                                end
                                w = [nr1 nr2 nr3 nr4];
                                x = [x_min x_max x_max x_min];    
                                y = [y_min y_min y_max y_max];
                            end
                            for j=1:lwe
                                if wygladz==0
                                    c(j) = Temp(w(j));
                                    %                                     if analiza==1
                                    %                                         c(j)=T(we(j));
                                    %                                     end
                                    %                                     if analiza==2
                                    %                                         c(j) = Tnu(k,we(j));
                                    %                                     end
                                end
                                if wygladz==1
                                    S = star(ns,x(j),y(j),Xs,Ys);
                                    g = 0;
                                    if wygladz==1
                                        g = par_wyg;
                                    end
                                    m = liczgwiaz;
                                    if m>ns
                                        m = ns;
                                    end
                                    M = mwls(x(j),y(j),Xs,Ys,S,m,stopien,g);
                                    c(j) = 0;
                                    for p=1:m
                                        c(j) = c(j) + M(1,p)*Temp(S(p,1));
                                    end
                                end   
                            end
                            %                         for j=2:2:8
                            %                             n1 = j-1;
                            %                             n2 = j+1;
                            %                             if(n2>8)
                            %                                 n2 = 1;
                            %                             end
                            %                             x(j)=(x(n1)+x(n2))/2;
                            %                             y(j)=(y(n1)+y(n2))/2;
                            %                             c(j)=(c(n1)+c(n2))/2;
                            %                         end
                            if ksztalt==2
                                s = patch(x(1:3),y(1:3),c(1:3),c(1:3));
                                set(s,'EdgeColor','interp');
                                s = patch(x([1 3 4]),y([1 3 4]),c([1 3 4]),c([1 3 4]));
                                set(s,'EdgeColor','interp');                               
                            end
                            if ksztalt==1
                                s = patch(x,y,c,c);
                                set(s,'EdgeColor','interp');
                            end
                        end
                        axis([xmin-0.1*(xmax+xmin) xmax+0.1*(xmax+xmin) ymin-0.1*(ymax+ymin) ymax+0.1*(ymax+ymin)]);
                        t = title('');
                        if analiza==1
                            set(t,'String',{'ROZK£AD TEMPERATURY',['T_m_i_n = ' num2str(min(T)) ', T_m_a_x = ' num2str(max(T))]});
                        end
                        if analiza==2
                            set(t,'String',{['ROZK£AD TEMPERATURY dla t = ' num2str(dt*k) ' / ' num2str(dt*N)],['T_m_i_n = ' num2str(min(Tnu(k,:))) ', T_m_a_x = ' num2str(max(Tnu(k,:)))]});
                        end
                        set(t,'FontSize',12);
                        colorbar();
                        if warwez==1
                            if analiza==1
                                Temp = T;
                            end
                            if analiza==2
                                Temp = Tnu(k,:);
                            end
                            for i=1:ns
                                kol = 'b';
                                if WBs(i,1)>0
                                    kol = 'r';
                                end
                                plot3([Xs(i) Xs(i)],[Ys(i) Ys(i)],[0 Temp(i)],'k-.');
                                plot3(Xs(i),Ys(i),0,'k.','MarkerFaceColor','k','MarkerSize',3);
                                plot3(Xs(i),Ys(i),Temp(i),[kol 's'],'MarkerFaceColor',kol,'MarkerSize',3);
                            end
                        end
                        view(3);
                        if analiza==2                            
                            Z = axis;
                            axis([Z(1:4) Tmin Tmax]);
                            pause(0.001);
                        end
                    end
                end
                if op1==2
                    pocz = 1;
                    N = 1;                                                     
                    s1 = 'SK£ADOWA Qx';
                    s2 = 'SK£ADOWA Qy';
                    s3 = 'WSZYSTKIE SK£ADOWE I WYPADKOWA';
                    s = s1;
                    dec = questdlg('PREZENTACJA WYNIKÓW STRUMIENIA CIEP£A','STRUMIEÑ CIEP£A',s1,s2,s3,s);                         
                    if strcmp(dec,s1)==1                           
                        poczk = 1;                          
                        konk = 1;
                    end
                    if strcmp(dec,s2)==1                           
                        poczk = 2;
                        konk = 2;                         
                    end                       
                    if strcmp(dec,s3)==1                           
                        poczk = 1;
                        konk = 3;                            
                    end 
                    il = konk-poczk+1;
                    if analiza==2
                        N = nt;
                        s1 = 'ROZK£AD STRUMIENIA W CZASIE CA£EGO PROCESU';
                        s2 = 'KOÑCOWY ROZK£AD STRUMIENIA';
                        s = s1;
                        dec = questdlg('PREZENTACJA WYNIKÓW STRUMIENIA CIEP£A','STRUMIEÑ CIEP£A',s1,s2,'POWRÓT',s);
                        il1 = 0;
                        if strcmp(dec,s1)==1
                            il1 = 1;
                            pocz = 1;
                        end
                        if strcmp(dec,s2)==1
                            il1 = 1;
                            pocz = nt;
                        end
                        if il1==0
                            pocz = N+1;
                            close all;
                        end                        
                    end
                    mes = msgbox('TRWAJ¥ OBLICZENIA, PROSZÊ CZEKAÆ...','INFORMACJA','help','replace');
                    D = diag([kx ky]);
                    close all
                    fig = figure(4);
                    set(fig,'Position',get(0,'ScreenSize'));                    
                    x = zeros(Nes,lwe);
                    y = zeros(Nes,lwe);                    
                    c = zeros(lwe,1);
                    w = zeros(lwe+1,1);
                   
                    Qw = zeros(N-pocz+1,ns,konk-poczk+1);
                    Qe = zeros(N-pocz+1,Nes,lwe,konk-poczk+1);    
                    Qs = zeros(N-pocz+1,Nes,konk-poczk+1);    
                    xs = zeros(Nes,1);
                    ys = zeros(Nes,1);
                    for h=pocz:N
                        if analiza==1
                            Temp = T;
                        end
                        if analiza==2
                            Temp = Tnu(h,:)';
                        end                        
                        if met==1
                            for i=1:Nes
                                w(1:lwe) = Incs(i,:);
                                w(lwe+1) = Incs(i,1);
                                we = Incse(i,:);
                                dlx = zeros(lwe,1);
                                dly = zeros(lwe,1);
                                x0 = 10^16;
                                y0 = 10^16;
                                for j=1:lwe
                                    dlx(j) = abs(Xs(w(j)) - Xs(w(j+1)));
                                    dly(j) = abs(Ys(w(j)) - Ys(w(j+1)));
                                    if Xs(w(j))<x0
                                        x0 = Xs(we(j));
                                    end
                                    if Ys(w(j))<y0
                                        y0 = Ys(we(j));
                                    end
                                end
                                a = max(dlx);
                                b = max(dly);
                                if ksztalt==1
                                    x0 = 0;
                                    y0 = 0;
                                end
                                for j=1:lwe
                                    x(i,j)=Xs(w(j));
                                    y(i,j)=Ys(w(j));
                                    %if wygladz==0
                                        qe = q(Xs(w(j))-x0,Ys(w(j))-y0,a,b,Temp,we(1:lsse),D,Xs,Ys,lsse);
                                        zes = [qe(1) qe(2) sqrt(qe(1)^2 + qe(2)^2)];
                                        if il==3
                                            Qe(h,i,j,1) = zes(1);
                                            Qe(h,i,j,2) = zes(2);
                                            Qe(h,i,j,3) = zes(3);
                                        else
                                            Qe(h,i,j,1) = zes(poczk);
                                        end
                                    %end
                                end
                                if wygladz==1
                                    if h==1
                                        xs(i) = sum(x(i,:))/lwe;
                                        ys(i) = sum(y(i,:))/lwe;
                                    end
                                    qe = q(xs(i)-x0,ys(i)-y0,a,b,Temp,we(1:lsse),D,Xs,Ys,lsse);
                                    zes = [qe(1) qe(2) sqrt(qe(1)^2 + qe(2)^2)];
                                    if il==3
                                        Qs(h,i,1) = zes(1);
                                        Qs(h,i,2) = zes(2);
                                        Qs(h,i,3) = zes(3);
                                    else
                                        Qs(h,i,1) = zes(poczk);
                                    end
                                        
                                end
                            end
                            if wygladz==1
                                for i=1:ns
                                    S = star(Nes,Xs(i),Ys(i),xs,ys);
                                    g = par_wyg;
                                    m = liczgwiaz;
                                    if m>Nes
                                        m = Nes;
                                    end
                                    M = mwls(Xs(i),Ys(i),xs,ys,S,m,1,g);
                                    zes = zeros(3,1);
                                    for p=1:m
                                        for k=poczk:konk
                                            zes(k-poczk+1) = zes(k-poczk+1) + M(1,p)*Qs(h,S(p,1),k-poczk+1);
                                        end                                        
                                    end          
                                    for k=poczk:konk
                                        Qw(h,i,k-poczk+1) = zes(k-poczk+1);
                                    end
                                end
                            end
                        end  
                        if met==2
                            for i=1:ns                                
                                S = star(ns,Xs(i),Ys(i),Xs,Ys);
                                g = 0;
                                if wygladz==1
                                    g = par_wyg;
                                end
                                m = liczgwiaz;
                                if m>ns
                                    m = ns;
                                end
                                M = mwls(Xs(i),Ys(i),Xs,Ys,S,m,1,g);
                                z1 = 0;
                                z2 = 0;
                                for p=1:m
                                    z1 = z1 + M(2,p)*Temp(S(p,1));
                                    z2 = z2 + M(3,p)*Temp(S(p,1));
                                end
                                zes = [z1 z2 sqrt(z1^2 + z2^2)]; 
                                if il==3
                                    Qw(h,i,1) = zes(1);
                                    Qw(h,i,2) = zes(2);
                                    Qw(h,i,3) = zes(3);
                                else
                                    Qw(h,i,1) = zes(poczk);
                                end
                            end
                        end
                    end                    
                    qmin = zeros(konk-poczk+1,1);
                    qmax = zeros(konk-poczk+1,1);
                    for k=poczk:konk
                        if met==1% & wygladz==0
                            qmin(k-poczk+1) = min(min(min(Qe(:,:,:,k-poczk+1))));
                            qmax(k-poczk+1) = max(max(max(Qe(:,:,:,k-poczk+1))));
                        end
                        if met==2% | wygladz==1
                            qmin(k-poczk+1) = min(min(Qw(:,:,k-poczk+1)));
                            qmax(k-poczk+1) = max(max(Qw(:,:,k-poczk+1)));
                        end
                    end
                    try
                        close(mes);
                    catch
                    end
                    for h=pocz:N      
                        for k=poczk:konk
                            if il==3
                                subplot(2,2,k);
                            end                           
                            hold off
                            plot3(0,0,0,'w');                            
                            hold on   
                            grid on
                            if analiza==1
                                axis([xmin-0.1*(xmax+xmin) xmax+0.1*(xmax+xmin) ymin-0.1*(ymax+ymin) ymax+0.1*(ymax+ymin)]);
                            end                                                        
                            if analiza==2
                                view(3);
                            end
                        end                       
                        for k=poczk:konk
                            if il==3
                                subplot(2,2,k);
                                hold on
                                grid on
                            end    
                            if analiza==2
                                view(3);
                            end
                            for i=1:Nes                                  
                                if met==1 & wygladz==0
                                    for p=1:lwe
                                        c(p) = Qe(h,i,p,k-poczk+1);
                                    end
                                end
                                if met==2 | wygladz==1
                                    w = Incs(i,:);
                                    x = zeros(lwe,1);
                                    y = zeros(lwe,1);
                                    for p=1:lwe
                                        c(p) = Qw(h,w(p),k-poczk+1);
                                        x(p) = Xs(w(p));
                                        y(p) = Ys(w(p));
                                    end
                                end
                               % c = Qe(h,i,:,k);
                               % x(i,:)
                               % y(i,:)
                                %c
%                                 pause
                                if met==1 & wygladz==0
                                    if ksztalt==1 | ksztalt==2
                                        s = patch(x(i,:),y(i,:),c,c);
                                        set(s,'EdgeColor','interp');        
                                    end
                                    if ksztalt==2                                  
                                        s = patch(x(i,[1 2 3])',y(i,[1 2 3])',c([1 2 3]),c([1 2 3]));
                                        set(s,'EdgeColor','interp');  
                                        s = patch(x(i,[1 3 4])',y(i,[1 3 4])',c([1 3 4]),c([1 3 4]));
                                        set(s,'EdgeColor','interp');   
                                    end
                                end
                                if met==2 | wygladz==1
                                    if ksztalt==1
                                        s = patch(x,y,c,c);
                                        set(s,'EdgeColor','interp');        
                                    end
                                    if ksztalt==2
                                        s = patch(x(1:3),y(1:3),c(1:3),c(1:3));
                                        set(s,'EdgeColor','interp');
                                        s = patch(x([1 3 4]),y([1 3 4]),c([1 3 4]),c([1 3 4]));
                                        set(s,'EdgeColor','interp');
                                    end
                                end                                                          
                            end
                            t = title('');
                            if k==1
                                naz = 'SK£ADOWA Q_x';
                            end
                            if k==2
                                naz = 'SK£ADOWA Q_y';
                            end
                            if k==3
                                naz = 'WARTOŒÆ ŒREDNIA Q = (Q_x^2 + Q_y^2)^0^.^5';
                            end
                            sp = '';
                            if analiza==2
                                sp = [' dla t = ' num2str(dt*h) ' / ' num2str(dt*N)];
                            end
                            set(t,'String',{['ROZK£AD WEKTORA STRUMIENIA CIEP£A ' sp],' ',naz,['Q_m_i_n = ' num2str(qmin(k-poczk+1)) ', Q_m_a_x = ' num2str(qmax(k-poczk+1))]});
                            set(t,'FontSize',12);
                            colorbar();
                        end                       
                        for k=poczk:konk        
                            if il==3
                                subplot(2,2,k);
                            end    
                            if analiza==2                                
                                Z = axis;
                                axis([Z(1:4) qmin(k-poczk+1) qmax(k-poczk+1)]);                                
                                view(3);
                                pause(0.1);                                                               
                            end           
                        end                     
                    end
                        
                                
                            %                             for i=1:Nes
                            %
                            %                                 for j=1:lsse
                            %                                     x(j)=Xs(we(j));
                            %                                     y(j)=Ys(we(j));
                            %                                     if wygladz==0 & met==1
                            %                                         qe = q(Xs(we(j))-x0,Ys(we(j))-y0,a,b,Temp,we(1:lsse),D,Xs,Ys,lsse);
                            %                                         if k==1
                            %                                             z = qe(1);
                            %                                         end
                            %                                         if k==2
                            %                                             z = qe(2);
                            %                                         end
                            %                                         if k==3
                            %                                             z = sqrt(qe(2)^2 + qe(1)^2);
                            %                                         end
                            %                                         c(j)=z;
                            %                                     end
                            %                                     if wygladz==1 | met==2
                            %                                         S = star(ns,x(j),y(j),Xs,Ys);
                            %                                         g = 0;
                            %                                         if wygladz==1
                            %                                             g = par_wyg;
                            %                                         end
                            %                                         m = liczgwiaz;
                            %                                         if m>ns
                            %                                             m = ns;
                            %                                         end
                            %                                         M = mwls(x(j),y(j),Xs,Ys,S,m,1,g);
                            %                                         z1 = 0;
                            %                                         z2 = 0;
                            %                                         for p=1:m
                            %                                             z1 = z1 + M(2,p)*Temp(S(p,1));
                            %                                             z2 = z2 + M(3,p)*Temp(S(p,1));
                            %                                         end
                            %                                         if k==1
                            %                                             c(j) = z1;
                            %                                         end
                            %                                         if k==2
                            %                                             c(j) = z2;
                            %                                         end
                            %                                         if k==3
                            %                                             c(j) = sqrt(z1^2 + z2^2);
                            %                                         end
                            %                                     end
                            %                                 end                                                                                                                                    
                end
            end
        else
            uiwait(msgbox('ZADANIE WYMAGA ROZWI¥ZANIA!','B£¥D','error','modal'));
        end
    end
    if op==13
        if obl==1
            pr = {'NAZWA PLIKU (BEZ ROZSZERZENIA) : '};
            odp = inputdlg(pr,'ZAPIS',1,{''});
            if isempty(odp)==0
                naz = deal(odp{:});
                il = fopen([naz '.txt']);
                dec = 'Tak';
                if il>0
                    dec = questdlg('Plik ju¿ istnieje! Nadpisaæ?', ...
                        'PYTANIE', ...
                        'Tak','Nie','Nie');
                end
                if strcmp(dec,'Tak')==1
                    plik = fopen([naz '.txt'],'wt');
                    fprintf(plik,'%s\n','HEATMIL - WYNIKI ANALIZY MES i BMRS STACJONARNEGO I NIESTACJONARNEGO PRZEP£YWU  CIEP£A');
                    fprintf(plik,'%s\n','AUTOR PROGRAMU: S£AWOMIR MILEWSKI, PK WIL L5, slawek@L5.pk.edu.pl');
                    fprintf(plik,'%s\n','-----------------------------------------------------------------');
                    fprintf(plik,'%s\n',['Data: ' num2str(date)]);
                    fprintf(plik,'%s\n','-----------------------------------------------------------------');
                    fprintf(plik,'%s\n','DANE OBLICZENIOWE: ');
                    %         fprintf(plik,'%s\n',['GRUBOŒÆ TARCZY h = ' num2str(h) ' [m]']);
                    fprintf(plik,'%s\n',['PRZEWODNOŒÆ CIEPLNA kx = ' num2str(kx) ' [J/Cms]']);
                    fprintf(plik,'%s\n',['PRZEWODNOŒÆ CIEPLNA ky = ' num2str(ky) ' [J/Cms]']);
                    fprintf(plik,'%s\n',['GENERATOR CIEP£A W OBSZARZE f(x,y,t) = ' f ' [J/m2s]']);
                    fprintf(plik,'%s\n','-----------------------------------------------------------------');
                    if met==1
                        fprintf(plik,'%s\n\n','WYNIKI ANALIZY MES: ');
                    end
                    if met==2
                        fprintf(plik,'%s\n\n','WYNIKI ANALIZY BMRS: ');
                    end
                    fprintf(plik,'%s\n','WYNIKI WÊZ£OWE (TEMPERATURA I STRUMIEÑ CIEP£A WE WÊZ£ACH POD KONIEC PROCESU CIEPLNEGO)');
                    if analiza==1
                        Temp = T;
                    end
                    if analiza==2
                        Temp = Tn;
                    end
                    for i=1:ns
                        fprintf(plik,'%s\n',['wêze³ = ' sprintf('%4d',i) ', x = ' sprintf('%-8.4f',Xs(i)) ', y = ' sprintf('%-8.4f',Ys(i)) ', T  = ' sprintf('%-10.4f',Temp(i)) ', q = ' sprintf('%-10.4f',R(i))]);
                    end
                    if met==1
                        fprintf(plik,'\n%s\n','WYNIKI ELEMENTOWE (STRUMIEÑ CIEP£A WE WÊZ£ACH ELEMENTU ORAZ W JEGO ŒRODKU CIÊ¯KOŒCI POD KONIEC PROCESU CIEPLNEGO)');
                        for i=1:Nes
                            fprintf(plik,'%s\n',['element = ' num2str(i)]);
                            %                         t1 = ['(' sprintf('%10.5f',QG(i,1,4)) ' ; ' sprintf('%10.5f',QG(i,2,4)) ')' '                   ' '(' sprintf('%10.5f',QG(i,1,3)) ' ; ' sprintf('%10.5f',QG(i,2,3)) ')'];
                            %                         t2 = ['                      ' '(' sprintf('%10.5f',QG(i,1,5)) ' ; ' sprintf('%10.5f',QG(i,2,5)) ')'];
                            %                         t3 = ['(' sprintf('%10.5f',QG(i,1,1)) ' ; ' sprintf('%10.5f',QG(i,2,1)) ')' '                   ' '(' sprintf('%10.5f',QG(i,1,2)) ' ; ' sprintf('%10.5f',QG(i,2,2)) ')'];
                            if ksztalt==2
                                t1 = ['(' sprintf('%10.5f',QG(i,1,4)) ' ; ' sprintf('%10.5f',QG(i,2,4)) ')' '                   ' '(' sprintf('%10.5f',QG(i,1,3)) ' ; ' sprintf('%10.5f',QG(i,2,3)) ')'];
                                t2 = ['                            ' '(' sprintf('%10.5f',QG(i,1,5)) ' ; ' sprintf('%10.5f',QG(i,2,5)) ')'];
                                t3 = ['(' sprintf('%10.5f',QG(i,1,1)) ' ; ' sprintf('%10.5f',QG(i,2,1)) ')' '                   ' '(' sprintf('%10.5f',QG(i,1,2)) ' ; ' sprintf('%10.5f',QG(i,2,2)) ')'];
                            end
                            if ksztalt==1
                                t1 = ['                            ' '(' sprintf('%10.5f',QG(i,1,1)) ' ; ' sprintf('%10.5f',QG(i,2,1)) ')'];
                                t2 = ['                            ' '(' sprintf('%10.5f',QG(i,1,4)) ' ; ' sprintf('%10.5f',QG(i,2,4)) ')'];
                                t3 = ['(' sprintf('%10.5f',QG(i,1,2)) ' ; ' sprintf('%10.5f',QG(i,2,2)) ')' '                   ' '(' sprintf('%10.5f',QG(i,1,3)) ' ; ' sprintf('%10.5f',QG(i,2,3)) ')'];
                            end
                            fprintf(plik,'%s\n',t1);
                            fprintf(plik,'%s\n',t2);
                            fprintf(plik,'%s\n\n',t3);
                        end
                    end
                    fclose(plik);
                end
            end
        else
            uiwait(msgbox('ZADANIE WYMAGA ROZWI¥ZANIA!','B£¥D','error','modal'));
        end
    end
    if op==12
        pr = {'NAZWA PLIKU (BEZ ROZSZERZENIA) : '};
        odp = inputdlg(pr,'ZAPIS',1,{naz_plik});
        if isempty(odp)==0
            naz = deal(odp{:});
            il = fopen([naz '.txt']);
            dec = 'Tak';
            if il>0
                dec = questdlg('Plik ju¿ istnieje! Nadpisaæ?', ...
                    'PYTANIE', ...
                    'Tak', 'Nie','Nie');
            end
            if strcmp(dec,'Tak')==1
                plik = fopen([naz '.txt'],'wt');
                naz_plik = naz;
                fprintf(plik,'%d\n',2011);
                fprintf(plik,'%d\n',nk);
                fprintf(plik,'%d\n',ilk);
                fprintf(plik,'%d\n',ilkT);
                fprintf(plik,'%d\n',ilkQ);
                fprintf(plik,'%d\n',ns);
                fprintf(plik,'%d\n',Nes);
                fprintf(plik,'%.16f\n',h);
                fprintf(plik,'%.16f\n',kx);
                fprintf(plik,'%.16f\n',ky);
                fprintf(plik,'%s\n',f);
                fprintf(plik,'%d\n',nh);
                fprintf(plik,'%d\n',siat);
                if nk>0
                    fprintf(plik,'%.16f\n',Xk);
                    fprintf(plik,'%.16f\n',Yk);
                end
                if ilk>0
                    fprintf(plik,'%d\n',Kr(1:ilk,:)');
                end
                if ilkT>0
                    fprintf(plik,'%.16f\n',WBTk');
                end
                if ilkQ>0
                    fprintf(plik,'%.16f\n',WBQk');
                end
                if ns>0
                    fprintf(plik,'%.16f\n',Xs);
                    fprintf(plik,'%.16f\n',Ys);
                    %if siat==2
                    fprintf(plik,'%d\n',WBs');
                    %end
                end
                fprintf(plik,'%d\n%d\n%d\n%d\n%d\n%d\n%d\n%d\n%d\n',met,analiza,schemat,sfor,liczgwiaz,sp_pocz,stopien,ksztalt,wygladz);
                fprintf(plik,'%.16f\n',par_wyg);
                %fprintf(plik,'%d\n%d\n',lwe,lsse);
                fprintf(plik,'%d\n',lwe);
                if Nes>0
                    fprintf(plik,'%d\n',Incs');
                    %fprintf(plik,'%d\n',Incse');
                end
                fprintf(plik,'%s\n',wzor_pocz);
                fprintf(plik,'%.16f\n%.16f\n%.16f\n%d\n%.16f\n%.16f\n',t0,dt,tk,nt,cw,amp_zab);
                fclose(plik);
                uiwait(msgbox(['DOKONANO ZAPISU DO PLIKU "' naz '.txt"'],'INFO','help','modal'));
                zapis = 1;
            end
        end
    end
    if op==14
        d = dir;
        str = {d.name};
        PLIK = struct([]);
        ilplik = 0;
        for i=1:length(str)
            naz = char(str(i));
            if length(naz)>4
                if strcmp(naz(end-3:end),'.txt')==1
                    ilplik = ilplik + 1;
                    PLIK(ilplik).naz = naz;
                end
            end
        end
        if ilplik==0
            uiwait(msgbox('BRAK DOSTÊPNYCH PLIKÓW','INFO','help','modal'));
        else
            [s,v] = listdlg('PromptString','Wybierz plik:',...
                'SelectionMode','single',...
                'ListString',{PLIK.naz});
            if isempty(s)==0
                bl = 0;
                try
                    Xk = [];
                    Yk = [];
                    Xs = [];
                    Ys = [];
                    Incs = [];
                    WBTk = [];
                    WBQk = [];
                    Kr = [];
                    naz = PLIK(s).naz;
                    naz_plik = naz(1:(end-4));
                    plik = fopen(naz,'rt');
                    wersja = fscanf(plik,'%d',[1 1]);
                    if wersja~=2011
                        nk = wersja;
                    end
                    if wersja==2011
                        nk = fscanf(plik,'%d',[1 1]);
                    end
                    ilk = fscanf(plik,'%d',[1 1]);
                    ilkT = fscanf(plik,'%d',[1 1]);
                    ilkQ = fscanf(plik,'%d',[1 1]);
                    ns = fscanf(plik,'%d',[1 1]);
                    Nes = fscanf(plik,'%d',[1 1]);
                    h = fscanf(plik,'%f',[1 1]);
                    kx = fscanf(plik,'%f',[1 1]);
                    ky = fscanf(plik,'%f',[1 1]);
                    if wersja==2011
                        f = fscanf(plik,'%s',[1 1]);
                    end
                    if wersja~=2011
                        z = fscanf(plik,'%f',[1 1]);
                        f = num2str(z);
                    end
                    nh = fscanf(plik,'%d',[1 1]);
                    siat = fscanf(plik,'%d',[1 1]);
                    if nk>0
                        Xk = fscanf(plik,'%f',[nk 1]);
                        Yk = fscanf(plik,'%f',[nk 1]);
                    end
                    if ilk>0
                        Kr = zeros(ilk,2);
                        for i=1:ilk
                            n1 = fscanf(plik,'%d',[1 1]);
                            n2 = fscanf(plik,'%d',[1 1]);
                            Kr(i,:) = [n1 n2];
                        end
                    end
                    if ilkT>0
                        WBTk = zeros(ilkT,2);
                        for i=1:ilkT
                            nr = round(fscanf(plik,'%f',[1 1]));
                            te = fscanf(plik,'%f',[1 1]);
                            WBTk(i,:) = [nr te];
                        end
                    end
                    if ilkQ>0
                        WBQk = zeros(ilkQ,2);
                        for i=1:ilkQ
                            nr = round(fscanf(plik,'%f',[1 1]));
                            qe = fscanf(plik,'%f',[1 1]);
                            WBQk(i,:) = [nr qe];
                        end
                    end
                    if ns>0
                        Xs = fscanf(plik,'%f',[ns 1]);
                        Ys = fscanf(plik,'%f',[ns 1]);
                        %if siat==2
                        WBs = zeros(ns,3);
                        for i=1:ns
                            n1 = fscanf(plik,'%d',[1 1]);
                            n2 = fscanf(plik,'%d',[1 1]);
                            n3 = fscanf(plik,'%d',[1 1]);
                            WBs(i,:) = [n1 n2 n3];
                        end
                        %end
                    end
                    
                    if wersja==2011
                        met = fscanf(plik,'%d',[1 1]);
                        analiza = fscanf(plik,'%d',[1 1]);
                        schemat = fscanf(plik,'%d',[1 1]);
                        sfor = fscanf(plik,'%d',[1 1]);
                        liczgwiaz = fscanf(plik,'%d',[1 1]);
                        sp_pocz = fscanf(plik,'%d',[1 1]);
                        stopien = fscanf(plik,'%d',[1 1]);
                        ksztalt = fscanf(plik,'%d',[1 1]);
                        wygladz = fscanf(plik,'%d',[1 1]);
                        par_wyg = fscanf(plik,'%f',[1 1]);
                        lwe = fscanf(plik,'%d',[1 1]);
                    end
                    if wersja~=2011
                        met = 1;
                        analiza = 1;
                        schemat = 2;
                        sfor = 1;
                        liczgwiaz = 16;
                        sp_pocz = 2;
                        stopien = 1;
                        ksztalt = 2;
                        wygladz = 0;
                        par_wyg = 0.5;
                        lwe = 4;
                    end
                    
                    if Nes>0
                        Incs = zeros(Nes,lwe);
                        for i=1:Nes
                            n1 = fscanf(plik,'%d',[1 1]);
                            n2 = fscanf(plik,'%d',[1 1]);
                            n3 = fscanf(plik,'%d',[1 1]);
                            if ksztalt==2
                                n4 = fscanf(plik,'%d',[1 1]);
                                Incs(i,:) = [n1 n2 n3 n4];
                            end
                            if ksztalt==1
                                Incs(i,:) = [n1 n2 n3];
                            end
                        end
                    end
                    if wersja==2011
                        wzor_pocz = fscanf(plik,'%s',[1 1]);
                        t0 = fscanf(plik,'%f',[1 1]);
                        dt = fscanf(plik,'%f',[1 1]);
                        tk = fscanf(plik,'%f',[1 1]);
                        nt = fscanf(plik,'%d',[1 1]);
                        cw = fscanf(plik,'%f',[1 1]);
                        if cw==0
                            cw = 1;
                        end
                        amp_zab = fscanf(plik,'%f',[1 1]);
                    end
                    if wersja~=2011
                        t0 = 0;
                        dt = 0.02;
                        nt = 10;
                        tk = dt*nt;
                        cw = 1;
                        wzor_pocz = '0';
                        amp_zab = 0;
                    end
                    fclose(plik);
                catch
                    bl = 1;
                end
                if bl==0
                    uiwait(msgbox(['DOKONANO ODCZYTU Z PLIKU "' naz '"'],'INFO','help','modal'));
                    nspop = 0;
                else
                    uiwait(msgbox(['PLIK "' naz '" MA NIEW£AŒCIW¥ STRUKTURÊ'],'B£¥D','error','modal'));
                    siat = 0;
                    Xk = [];
                    Yk = [];
                    Xs = [];
                    Ys = [];
                    Incs = [];
                    WBTk = [];
                    WBQk = [];
                    Kr = [];
                    ilkT = 0;
                    ilkQ = 0;
                    Nes = 0;
                    ns = 0;
                    nk = 0;
                    nspop = 0;
                    ilk = 0;
                end
                close all;
                obl = 0;
            end
        end
    end
end
try
    close(msg);
catch
end
close all;

function [QG,PG] = postprocess(ns,Xs,Ys,Incs,Incse,Nes,lwe,lsse,ksztalt,met,sfor,stopien,liczgwiaz,kx,ky,T)
D = diag([kx ky]);
QG = zeros(Nes,2,lwe+1);
PG = zeros(Nes,2,lwe+1);
for i=1:Nes
    w = Incs(i,:);
    w(lwe+1) = Incs(i,1);
    we = Incse(i,:);
    we(lsse+1) = Incse(i,1);
    dlx = zeros(lwe,1);
    dly = zeros(lwe,1);
    for j=1:lwe
        dlx(j) = abs(Xs(w(j)) - Xs(w(j+1)));
        dly(j) = abs(Ys(w(j)) - Ys(w(j+1)));
    end
    a = max(dlx);
    b = max(dly);
    %i
    Q = zeros(2,lwe+1);
    for j=1:lwe+1
        if ksztalt==2
            if j==1
                x = 0;
                y = 0;
            end
            if j==2
                x = a;
                y = 0;
            end
            if j==3
                x = a;
                y = b;
            end
            if j==4
                x = 0;
                y = b;
            end
            if j==5
                x = a/2;
                y = b/2;
            end
        end
        if ksztalt==1
            if j<=3
                x = Xs(w(j));
                y = Ys(w(j));
            end
            if j==4
                x = (Xs(w(1)) + Xs(w(2)) + Xs(w(3)))/3;
                y = (Ys(w(1)) + Ys(w(2)) + Ys(w(3)))/3;
            end
            if ksztalt==2
                x = x - min(Xs);
                y = y - min(Ys);
            end
        end
        qe = q(x,y,a,b,T,we,D,Xs,Ys,lsse);
        Q(:,j) = qe;
        PG(i,j,1) = x;
        PG(i,j,2) = y;
    end
    for j=1:lwe+1
        for k=1:2
            QG(i,k,j) = Q(k,j);
        end
    end
    %Q
end

function [T,R] = solver(ns,Xs,Ys,WBs,Incs,Incse,Nes,lwe,lsse,ksztalt,stopien,met,sfor,liczgwiaz,ilkT,WBTk,ilkQ,WBQk,kx,ky,f)
global dt nt t0 T0 Tp Tn analiza;
global sp_pocz wzor_pocz schemat tk czas cw;
D = diag([kx ky]);
KG = zeros(ns,ns);
FG = zeros(ns,1);
xt = zeros(lwe,1);
yt = zeros(lwe,1);

for i=1:Nes
    w = Incs(i,:);
    w(lwe+1) = Incs(i,1);
    we = Incse(i,:);
    we(lsse+1) = Incse(i,1);
    dlx = zeros(lwe,1);
    dly = zeros(lwe,1);
    for j=1:lwe
        dlx(j) = abs(Xs(w(j)) - Xs(w(j+1)));
        dly(j) = abs(Ys(w(j)) - Ys(w(j+1)));
        xt(j) = Xs(w(j));
        yt(j) = Ys(w(j));
    end
    a = max(dlx);
    b = max(dly);
    if ksztalt==2
        Ng = 4;
        [PGx,PGy,WAG,Jx,Jy] = cal_pros(Ng,0,a,0,b);
        NoGp = Ng*Ng;
        J = Jx*Jy;
    end
    if ksztalt==1
        Ng = 7;
        [PGx,PGy,WAG] =  cal_troj(Ng,Xs,Ys,w(1:3));
        wek = [Xs(w(1)) - Xs(w(2)), Ys(w(1)) - Ys(w(2)); Xs(w(2)) - Xs(w(3)), Ys(w(2)) - Ys(w(3)); Xs(w(3))-Xs(w(1)), Ys(w(3))-Ys(w(1))];
        dl = [sqrt(wek(1,1)^2 + wek(1,2)^2), sqrt(wek(2,1)^2 + wek(2,2)^2), sqrt(wek(3,1)^2 + wek(3,2)^2)];
        p = 0.5*(dl(1) + dl(2) + dl(3));
        J = sqrt(p*(p-dl(1))*(p-dl(2))*(p-dl(3))); % Jakobian - pole trójk¹ta
        NoGp = Ng;
    end
    if met==1
        ke = zeros(lsse,lsse);
        fe = zeros(lsse,1);
        if analiza==2
            Te = zeros(lsse,1);
            for j=1:lsse
                Te(j) = Tp(we(j));
            end
        end
        Tx = 0;
        Ty = 0; 
        for p=1:NoGp
            xg = PGx(p);
            yg = PGy(p);
            x = xg;
            y = yg;
            t = czas;
            f_war = eval(f);
            %                         if ksztalt==1
            %                             xg = xg - min(Xs);
            %                             yg = yg - min(Ys);
            %                         end
            Ne = N(xg,yg,a,b,we,Xs,Ys,lsse);            
            Be = B(xg,yg,a,b,we,Xs,Ys,lsse);
            if analiza==2
                Tx = Tx + Be(1,:)*Te;
                Ty = Ty + Be(2,:)*Te;
            end
            if analiza==1 | (analiza==2 & schemat==2)
                ke = ke + J*WAG(p)*Be'*D*Be;
                if analiza==2
                    ke = ke + (cw/dt)*J*WAG(p)*Ne'*Ne;
                end
                fe = fe + J*WAG(p)*Ne'*f_war;
                if analiza==2
                    te = Ne*Te;
                    fe = fe + (cw/dt)*J*WAG(p)*Ne'*te;
                end
            end
            if analiza==2 & schemat==1
                ke = ke + (cw/dt)*J*WAG(p)*Ne'*Ne;
                Nex =  N(xg,yg,a,b,we,Xs,Ys,lsse,1);
                Ney =  N(xg,yg,a,b,we,Xs,Ys,lsse,2);     
                te = Ne*Te;
                fe = fe + J*WAG(p)*( -(kx*Nex'*Tx + ky*Ney'*Ty) + (cw/dt)*Ne'*te + Ne'*f_war);                               
            end
        end        
        %i
        %ke
        for j=1:lsse
            for k=1:lsse
                KG(we(j),we(k)) = KG(we(j),we(k)) + ke(j,k);
            end
            FG(we(j)) = FG(we(j)) + fe(j);
        end
    end
    if met==2 & sfor==1
        m = liczgwiaz;
        if m>ns
            m = ns;
        end
        for p=1:NoGp % loop over Gauss points in triangle
            xg = PGx(p); % x-coordinate of the Gauss point
            yg = PGy(p); % y-coordinate of the Gauss point
            if ksztalt==2
                xg = xg + min(xt);
                yg = yg + min(yt);
            end
            S = star(ns,xg,yg,Xs,Ys); % MFD star generation
            % MFD formulas generation (MWLS) for trial function (m star nodes, p=1)
            MF = mwls(xg,yg,Xs,Ys,S,m,stopien);
            % MFD formulas generation (MWLS) for test function (3 star nodes, p=1)
            stv = 1;
            if ksztalt==2
                stv = 1.5;
            end
            Mv = mwls(xg,yg,Xs,Ys,w(1:lwe)',lwe,stv);
            x = xg;
            y = yg;
            t = czas;
            f_war = eval(f);
            if analiza==2
                te = 0;
                for j=1:m
                    nr = S(j,1);
                    te = te + MF(1,j)*Tp(nr);
                end
            end
            for k=1:lwe % loop over d.o.f. of test function
                for j=1:m % loop over d.o.f. of trial function
                    nr = S(j,1);
                    KG(w(k),nr) = KG(w(k),nr) + J*WAG(p)*(kx*MF(2,j)*Mv(2,k) + ...
                        ky*MF(3,j)*Mv(3,k)); % aggregation of the global coefficient matrix
                    if analiza==2
                        KG(w(k),nr) = KG(w(k),nr) + (cw/dt)*J*WAG(p)*MF(1,j)*Mv(1,k);
                    end                            
                end
                % aggregation of the global right hand side vector
                FG(w(k)) = FG(w(k)) + f_war*J*WAG(p)*Mv(1,k);
                if analiza==2
                    FG(w(k)) = FG(w(k)) + (cw/dt)*J*WAG(p)*Mv(1,k)*te;
                end
            end
        end
    end
    %                     i
    %                     ke
end

%if siat==2 | siat==1
wbT = zeros(ns,2);
for i=1:ns
    if WBs(i,1)==1
        for k=1:2
            nr1 = WBs(i,1+k);
            il = 0;
            Te = 0;
            for j=1:ilkT
                nr2 = WBTk(j,1);
                if nr1==nr2
                    il = 1;
                    Te = WBTk(j,2);
                end
            end
            if il==1
                if analiza==1 | (analiza==2 & sp_pocz==1)
                    wbT(i,1) = wbT(i,1) + 1;
                    wbT(i,2) = wbT(i,2) + Te;
                end
                if analiza==2 & sp_pocz==2                    
                    wbT(i,1) = 1;
                    if schemat==2% | schemat==1
                        x = Xs(i);
                        y = Ys(i);
                        t = czas;
                        wbT(i,2) = eval(wzor_pocz);
                    end
                    if schemat==1
                        x = Xs(i);
                        y = Ys(i);
                        t = czas;
                        T2 = eval(wzor_pocz);
                        t = czas-dt;
                        T1 = eval(wzor_pocz);
                        wbT(i,2) = (T2-T1)/dt;
                    end
                end
            end
        end
    end
end
%WBs
for i=1:Nes
    w = Incs(i,:);
    w(lwe+1) = Incs(i,1);
    dlx = zeros(lwe,1);
    dly = zeros(lwe,1);
    for j=1:lwe
        dlx(j) = abs(Xs(w(j)) - Xs(w(j+1)));
        dly(j) = abs(Ys(w(j)) - Ys(w(j+1)));
    end
    for p=1:lwe
        n1 = w(p);
        n2 = w(p+1);
        if WBs(n1)==1 & WBs(n2)==1
            il1p = zeros(2,1);
            qe1p = zeros(2,1);
            for k=1:2
                nr1 = WBs(n1,1+k);
                for j=1:ilkQ
                    nr2 = WBQk(j,1);
                    if nr1==nr2
                        il1p(k) = nr2;
                        qe1p(k) = WBQk(j,2);
                    end
                end
            end
            il2p = zeros(2,1);
            qe2p = zeros(2,1);
            for k=1:2
                nr1 = WBs(n2,1+k);
                for j=1:ilkQ
                    nr2 = WBQk(j,1);
                    if nr1==nr2
                        il2p(k) = nr2;
                        qe2p(k) = WBQk(j,2);
                    end
                end
            end
            il = 0;
            qe = 0;
            for k=1:2
                nr1 = il1p(k);
                for j=1:2
                    nr2 = il2p(j);
                    if nr1==nr2 & (nr1~=0) & (nr2~=0)
                        il = nr1;
                        qe = 0.5*(qe1p(k) + qe2p(j));
                    end
                end
            end
            if il>0
                dl = sqrt( (Xs(n1) - Xs(n2))^2 + (Ys(n1) - Ys(n2))^2 );
                if stopien==1 | met==2
                    FG(n1) = FG(n1) + qe*dl/2;
                    FG(n2) = FG(n2) + qe*dl/2;
                end
                if stopien==2 & met==1
                    FG(n1) = FG(n1) + qe*dl/6;
                    FG(n2) = FG(n2) + qe*dl/6;
                    n3 = 0;
                    k = 0;
                    x = 0.5*(Xs(n1)+Xs(n2));
                    y = 0.5*(Ys(n1)+Ys(n2));
                    while (n3==0)
                        k = k+1;
                        d = sqrt( (x-Xs(k))^2 + (y-Ys(k))^2 );
                        if d<eps
                            n3 = k;
                        end
                    end
                    FG(n3) = FG(n3) + 2*qe*dl/3;
                end
            end
        end
    end
end
%end

ilwb = 0;
wb = zeros(ns,2);
for i=1:ns
    if wbT(i,1)>0
        ilwb = ilwb + 1;
        wb(ilwb,:) = [i wbT(i,2)/wbT(i,1)];
    end
end
wb(ilwb+1:end,:) = [];
if ns<=30
    'MACIERZ SZTYWNOŒCI K'
    KG
    'WEKTOR OBCI¥¯EÑ'
    FG
end
[T,R] = solveq(KG,FG,wb);


function rysunek(fig,nk,Xk,Yk,ilk,Kr,opcja,WBTk,WBQk,Nes,Incs,Xs,Ys,lwe)
global stopien ns met;
global dt nt t0 T0 Tp Tn analiza sp_pocz wzor_pocz schemat tk;
mes = msgbox('RYSUNEK JEST ODŒWIE¯ANY, PROSZÊ CZEKAÆ...','INFORMACJA','help','replace');
if nargin<7
    opcja = 0;
    WBTk = [];
    WBQk = [];
    Nes = 0;
    Incs = [];
    Xs = [];
    Ys = [];
    lwe = 4;
end
xmax = max(Xk);
xmin = min(Xk);
ymax = max(Yk);
ymin = min(Yk);
figure(fig);
hold off
plot(0,0,'w.');
hold on
grid on
if analiza==2 & opcja==1
    if Nes>0 & isempty(T0)==0
        maxT0 = max(T0);
        minT0 = min(T0);
        if abs(maxT0-minT0)>eps
            x = zeros(lwe,1);
            y = zeros(lwe,1);
            c = zeros(lwe,1);
            for i=1:Nes
                w = Incs(i,:);
                w(lwe+1) = Incs(i,1);
                for j=1:lwe
                    x(j)=Xs(w(j));
                    y(j)=Ys(w(j));
                    c(j)=T0(w(j));
                end
                s = patch(x,y,c,c);
                set(s,'EdgeColor','interp');
            end
            %legend('TEMPERATURA POCZ¥TKOWA');
            colorbar();
        end
    end
end
if xmax~=xmin & ymax~=ymin
    %    axis([xmin-0.1*(xmax+xmin) xmax+0.1*(xmax+xmin) ymin-0.1*(ymax+ymin) ymax+0.1*(ymax+ymin)]);
end
for i=1:nk
    plot(Xk(i),Yk(i),'bo','LineWidth',2);
    t = text(Xk(i)+0.01*(xmax+xmin),Yk(i)+0.01*(ymax+ymin),num2str(i));
    set(t,'BackgroundColor','w');
    set(t,'FontSize',8);
    set(t,'Color','r');
end
for i=1:ilk
    n1 = Kr(i,1);
    n2 = Kr(i,2);
    plot([Xk(n1) Xk(n2)],[Yk(n1) Yk(n2)],'y-','LineWidth',2);
    xs = 0.5*(Xk(n1) + Xk(n2));
    ys = 0.5*(Yk(n1) + Yk(n2));
    t = text(xs,ys,num2str(i),'FontSize',10,'Color','k');
    set(t,'BackgroundColor','w');
end
if opcja==1
    for i=1:size(WBQk,1)
        %         n1 = WBQ(i,1);
        %         n2 = WBQ(i,2);
        %         qe = WBQ(i,3);
        nr = WBQk(i,1);
        n1 = Kr(nr,1);
        n2 = Kr(nr,2);
        qe = WBQk(i,2);
        plot([Xk(n1) Xk(n2)],[Yk(n1) Yk(n2)],'g-','LineWidth',2.5);
        xs = 0.5*(Xk(n1)+Xk(n2));
        ys = 0.5*(Yk(n1)+Yk(n2));
        s = num2str(nr);
        s1 = '';
        for j=1:length(s)
            s1(2*j-1) = '_';
            s1(2*j) = s(j);
        end
        t = text(xs,ys,['Q' s1 '=' num2str(qe)]);
        set(t,'BackgroundColor','w');
        set(t,'EdgeColor','g');
        set(t,'FontWeight','bold');
        set(t,'FontSize',12);
        set(t,'Color','g');
    end
    for i=1:size(WBTk,1)
        %         n1 = WBT(i,1);
        %         n2 = WBT(i,2);
        %         Te = WBT(i,3);
        nr = WBTk(i,1);
        n1 = Kr(nr,1);
        n2 = Kr(nr,2);
        Te = WBTk(i,2);
        plot([Xk(n1) Xk(n2)],[Yk(n1) Yk(n2)],'m-','LineWidth',2.5);
        xs = 0.5*(Xk(n1)+Xk(n2));
        ys = 0.5*(Yk(n1)+Yk(n2));
        s = num2str(nr);
        s1 = '';
        for j=1:length(s)
            s1(2*j-1) = '_';
            s1(2*j) = s(j);
        end
        t = text(xs,ys,['T' s1 '=' num2str(Te)]);
        set(t,'BackgroundColor','w');
        set(t,'EdgeColor','m');
        set(t,'FontWeight','bold');
        set(t,'FontSize',12);
        set(t,'Color','m');
    end
end
if opcja==2
    if met==1 | met==2
        w = zeros(1,lwe+1);
        for i=1:Nes
            w(1:lwe) = Incs(i,:);
            w(lwe+1) = Incs(i,1);
            xs = 0;
            ys = 0;
            for j=1:lwe
                styl = 'b-';
                if met==2
                    styl = 'g-.';
                end
                p = plot([Xs(w(j)) Xs(w(j+1))],[Ys(w(j)) Ys(w(j+1))],styl);
                if met==1
                    set(p,'LineWidth',1.5);                
                end
                xs = xs + Xs(w(j));
                ys = ys + Ys(w(j));
            end
            if Nes<=100 & met==1
                xs = xs/lwe;
                ys = ys/lwe;
                t = text(xs,ys,num2str(i));
                set(t,'FontSize',12);
            end
        end
    end
    %ns = length(Xs);
    for i=1:ns
        plot(Xs(i),Ys(i),'ko','LineWidth',2);
        if ns<=100
            t = text(Xs(i)+0.01*(xmax+xmin),Ys(i)+0.01*(ymax+ymin),num2str(i));
            set(t,'BackgroundColor','w');
            set(t,'FontSize',10);
            set(t,'Color','r');
        end
    end
end
dod = '';
if opcja==1
    dod = '+ WARUNKI BRZEGOWE';
end
if opcja==2
    dod = '+ DYSKRETYZACJA MES / BMRS';
end
title(['GEOMETRIA ZADANIA ' dod],'FontSize',12);
close(mes);

function [ns,Xs,Ys,Nes,Incs,WBs] = siatka(nh,ilk,Kr,Xk,Yk,ksztalt)
% if nargin<6
%     ksztalt=2;
% end
ilp = ilk;
Xp = zeros(ilp,1);
Yp = zeros(ilp,1);
for i=1:ilp
    nr = Kr(i,1);
    Xp(i) = Xk(nr);
    Yp(i) = Yk(nr);
end
Dlx = zeros(ilk,1);
Dly = zeros(ilk,1);
Dlxmin = 1;
Dlymin = 1;
eps = 10^-12;
ilx = 0;
ily = 0;
for i=1:ilk
    n1 = Kr(i,1);
    n2 = Kr(i,2);
    x1 = Xk(n1);
    y1 = Yk(n1);
    x2 = Xk(n2);
    y2 = Yk(n2);
    dlx = abs(x2-x1);
    dly = abs(y2-y1);
    if abs(dlx)>eps
        il = 0;
        for j=1:ilx
            if abs(Dlx(j,1)-dlx)<eps
                il = il + 1;
            end
        end
        if il==0
            ilx = ilx + 1;
            Dlx(ilx,1) = dlx;
            Dlxmin = Dlxmin*dlx;
        end
    end
    
    if abs(dly)>eps
        il = 0;
        for j=1:ily
            if abs(Dly(j,1)-dly)<eps
                il = il + 1;
            end
        end
        if il==0
            ily = ily + 1;
            Dly(ily,1) = dly;
            Dlymin = Dlymin*dly;
        end
    end
end
Dlx(ilx+1:end) = [];
Dly(ily+1:end) = [];


Dlxmin = Dlx(1);
for i=1:ilx-1
    if i==1
        a = Dlx(i,1);
    else
        a = Dlxmin;
    end
    b = Dlx(i+1,1);
    Dlxmin = NWD(a,b);
end

Dlymin = Dly(1);
for i=1:ily-1
    if i==1
        a = Dly(i,1);
    else
        a = Dlymin;
    end
    b = Dly(i+1,1);
    Dlymin = NWD(a,b);
end

xmax = max(Xp);
xmin = min(Xp);
ymax = max(Yp);
ymin = min(Yp);

nhc = nh;
hxc = Dlxmin/(nhc);
hyc = Dlymin/(nhc);
nxc = round((xmax - xmin)/hxc + 1);
nyc = round((ymax - ymin)/hyc + 1);
nc = nxc*nyc;

if nc>1000
    ns = 0;
    Xs = [];
    Ys = [];
    Nes = 0;
    Incs = [];
    WBs = [];
else
    Xc = zeros(nc,1);
    Yc = zeros(nc,1);
    x = xmin;
    y = ymin;
    j = 1;
    for i=1:nc
        Xc(i) = x;
        Yc(i) = y;
        x = x + hxc;
        j = j + 1;
        if j>nxc % przejœcie na jeden poziom wy¿ej
            j = 1;
            x = xmin;
            y = y + hyc;
        end
    end
    if ksztalt==2
        Pc = struct([]);
        il = 0;
        for j=1:nyc-1
            for i=1:nxc-1
                n1 = i + (j-1)*nxc;
                n2 = i + 1 + (j-1)*nxc;
                n3 = i + 1 + (j)*nxc;
                n4 = i + j*nxc;
                %             n1
                %             n2
                %             n3
                %             n4
                Pc(il+1).X = [Xc(n1) Xc(n2)];
                Pc(il+1).Y = [Yc(n1) Yc(n3)];
                Pc(il+1).w = [n1 n2 n3 n4];
                xs = 0.5*(sum(Pc(il+1).X));
                ys = 0.5*(sum(Pc(il+1).Y));
                Pc(il+1).c = pwew(xs,ys,Xp,Yp,ilp);
                if Pc(il+1).c==1
                    il = il + 1;
                end
            end
        end
        
        Xc1 = zeros(nc,1);
        Yc1 = zeros(nc,1);
        nc = il;
        eps = 10^-12;
        ilc = 0;
        for i=1:nc
            w = Pc(i).w;
            w1 = zeros(4,1);
            for j=1:4
                x = Xc(w(j));
                y = Yc(w(j));
                il = 0;
                for k=1:ilc
                    d = sqrt( (x - Xc1(k))^2 + (y - Yc1(k))^2 );
                    if d<eps
                        il = k;
                    end
                end
                if il==0
                    ilc = ilc + 1;
                    Xc1(ilc) = x;
                    Yc1(ilc) = y;
                    il = ilc;
                end
                w1(j) = il;
            end
            Pc(i).w = w1;
            Pc(i).X = [Xc1(w1(1)) Xc1(w1(2))];
            Pc(i).Y = [Yc1(w1(1)) Yc1(w1(3))];
        end
        Xc1(ilc+1:end) = [];
        Yc1(ilc+1:end) = [];
        Xc = Xc1;
        Yc = Yc1;
        clear Xc1 Yc1;
    end
    if ksztalt==1
        Xc1 = zeros(nc,1);
        Yc1 = zeros(nc,1);
        nc1 = 0;
        for i=1:nc
            c = pwew(Xc(i),Yc(i),Xp,Yp,ilp);
            c1 = pbrzeg(Xc(i),Yc(i),ilk,Kr,Xk,Yk);
            if c==1 | c1==1
                nc1 = nc1 + 1;
                Xc1(nc1) = Xc(i);
                Yc1(nc1) = Yc(i);
            end
        end
        Xc = Xc1(1:nc1);
        Yc = Yc1(1:nc1);
        ilc = nc1;
        T = delaunay(Xc,Yc);
        iltr = length(T);
        T1 = zeros(iltr,3);
        iltr1 = 0;
        for i=1:iltr
            x = zeros(3,1);
            y = zeros(3,1);
            for j=1:3
                nr = T(i,j);
                x(j) = Xc(nr);
                y(j) = Yc(nr);
            end
            xs = sum(x)/3;
            ys = sum(y)/3;
            c = pwew(xs,ys,Xp,Yp,ilp);
            if c==1
                iltr1 = iltr1+1;
                T1(iltr1,:) = T(i,:);
            end
        end
        iltr = iltr1;
        T = T1(1:iltr1,:);
    end
    WBc = zeros(ilc,3);
    for i=1:ilc
        [WBc(i,1) nrkr] = pbrzeg(Xc(i),Yc(i),ilk,Kr,Xk,Yk);
        if WBc(i,1)==1
            for j=1:length(nrkr)
                WBc(i,1+j) = nrkr(j);
            end
        end
    end
    Xs = Xc;
    Ys = Yc;
    WBs = WBc;
    ns = ilc;
    if ksztalt==2
        Nes = nc;
        Incs = zeros(Nes,4);
        for i=1:Nes
            Incs(i,:) = Pc(i).w;
        end
    end
    if ksztalt==1
        Nes = iltr;
        Incs = zeros(Nes,3);
        for i=1:Nes
            Incs(i,:) = T(i,:);
        end
    end
end


function [z] = NWD(a,b)
eps = 10^-12;
if abs(b)<eps
    z = a;
else
    z = NWD(b,mod(a,b));
end

function [c,nrkr] = pbrzeg(x,y,ilk,Kr,Xp,Yp)
wb = 0;
eps = 10^-12;
nrkr = zeros(ilk,1);
ilk1 = 0;
for k=1:ilk
    n1 = Kr(k,1);
    n2 = Kr(k,2);
    x1 = Xp(n1);
    y1 = Yp(n1);
    x2 = Xp(n2);
    y2 = Yp(n2);
    wb1 = 0;
    if abs(x2-x1)>eps
        a = (y2-y1)/(x2-x1);
        b = y2 - a*x2;
        if abs(a*x + b - y)<eps & x>=min([x1 x2])-eps & x<=max([x1 x2])+eps & y>=min([y1 y2])-eps & y<=max([y1 y2])+eps
            wb = 1;
            wb1 = 1;
        end
        if abs(a*x + b - y)<eps & wb == 0
            if abs(y2-y1)<eps
                if x>=min([x1 x2])-eps & x<=max([x1 x2])+eps
                    wb = 1;
                    wb1 = 1;
                end
            end
        end
    else
        if abs(x-x2)<eps & ( (y>=min([y1 y2]) & y<=max([y1 y2])) | abs(y-y1)<eps | abs(y-y2)<eps )
            wb = 1;
            wb1 = 1;
        end
    end
    if wb1==1
        ilk1 = ilk1 + 1;
        nrkr(ilk1) = k;
    end
end
nrkr(ilk1+1:end) = [];
c = wb;


function [c] = pwew(x,y,Xp,Yp,ilp)
c = 0;
p = ilp;
for k=1:ilp
    if ((Yp(k)>y) ~= (Yp(p)>y)) & (x<(Xp(p)-Xp(k))*(y-Yp(k))/(Yp(p)-Yp(k))+Xp(k))
        if c==0
            c = 1;
        else
            c = 0;
        end
    end
    p = k;
end

function [z] = B(x,y,a,b,w,Xs,Ys,lsse)
% dx = 0.0001;
% dy = 0.0001;
% Nx = (N(x+dx,y,a,b,w,Xs,Ys,lsse) - N(x-dx,y,a,b,w,Xs,Ys,lsse))/(2*dx);
% Ny = (N(x,y+dy,a,b,w,Xs,Ys,lsse) - N(x,y-dy,a,b,w,Xs,Ys,lsse))/(2*dy);
Nx = N(x,y,a,b,w,Xs,Ys,lsse,1);
Ny = N(x,y,a,b,w,Xs,Ys,lsse,2);
z = [Nx;Ny];

function [z] = N(x,y,a,b,w,Xs,Ys,lsse,poch)
if nargin<9
    poch = 0;
end
%z = [(x-a)*(y-b)/(a*b) -x*(y-b)/(a*b) x*y/(a*b) -(x-a)*y/(a*b)];
M = zeros(lsse,lsse);
X = zeros(lsse,1);
Y = zeros(lsse,1);
for i=1:lsse
    X(i) = Xs(w(i));
    Y(i) = Ys(w(i));
end
Xmin = min(X);
Ymin = min(Y);
for i=1:lsse
    xw = Xs(w(i));
    yw = Ys(w(i));
    if lsse==8 | lsse==4
        xw = Xs(w(i))-Xmin;
        yw = Ys(w(i))-Ymin;
    end
    if lsse==8
        M(i,:) = [1 xw yw xw^2 xw*yw yw^2 xw^2*yw xw*yw^2];
    end
    if lsse==6
        M(i,:) = [1 xw yw xw^2 xw*yw yw^2];
    end
    if lsse==4
        M(i,:) = [1 xw yw xw*yw];
    end
    if lsse==3
        M(i,:) = [1 xw yw];
    end
end
F = eye(lsse,lsse);
W = M\F;
z = zeros(1,lsse);
for i=1:lsse
    if poch==0
        z(i) = W(1,i) + (x)*W(2,i) + (y)*W(3,i);
        if lsse==4
            z(i) = z(i) + (x)*(y)*W(4,i);
        end
        if lsse==6 | lsse==8
            z(i) = z(i) + x^2*W(4,i) + x*y*W(5,i) + y^2*W(6,i);
            if lsse==8
                z(i) = z(i) + x^2*y*W(7,i) + x*y^2*W(8,i);
            end
        end
    end
    if poch==1
        z(i) = W(2,i);
        if lsse==4
            z(i) = z(i) + (y)*W(4,i);
        end
        if lsse==6 | lsse==8
            z(i) = z(i) + 2*x*W(4,i) + y*W(5,i);
            if lsse==8
                z(i) = z(i) + 2*x*y*W(7,i) + y^2*W(8,i);
            end
        end
    end
    if poch==2
        z(i) = W(3,i);
        if lsse==4
            z(i) = z(i) + (x)*W(4,i);
        end
        if lsse==6 | lsse==8
            z(i) = z(i) + x*W(5,i) + 2*y*W(6,i);
            if lsse==8
                z(i) = z(i) + x^2*W(7,i) + 2*x*y*W(8,i);
            end
        end
    end
end


function [z] = q(x,y,a,b,T,w,D,Xs,Ys,lsse)
Te = zeros(lsse,1);
for i=1:lsse
    Te(i) = T(w(i));
end
Be = B(x,y,a,b,w,Xs,Ys,lsse);
z = -(D*Be)*Te;

function [PG,wag] = cal(Ng)
if Ng==1
    wag = 2;
    PG = 0;
end
if Ng==2
    wag = [1 1];
    PG = [-1/sqrt(3) 1/sqrt(3)];
end
if Ng==3
    PG = [-1/sqrt(0.6) 0 1/sqrt(0.6)];
    wag = [ 5/9 8/9 5/9 ];
end
if Ng==4
    PG = [-0.861136311594053 -0.339981043584856 0.339981043584856 0.861136311594053];
    wag = [0.347854845137454 0.652145154862546 0.652145154862546 0.347854845137454];
end

function [PGx,PGy,WAG,Jx,Jy] = cal_pros(Ng,a,b,c,d)
[PG,wag] = cal(Ng);
k = 0;
Jx = 0.5*(b-a);
Jy = 0.5*(d-c);
WAG = zeros(Ng*Ng,1);
PGx = zeros(Ng*Ng,1);
PGy = zeros(Ng*Ng,1);
for i=1:Ng
    for j=1:Ng
        k = k + 1;
        WAG(k) = wag(i)*wag(j);
        PGx(k) = Jx*PG(i) + (b+a)/2;
        PGy(k) = Jy*PG(j) + (c+d)/2;
    end
end

function [PGx,PGy,WAG1] =  cal_troj(Ng1,X,Y,w) % funkcja przekazuj¹ca rozk³ad p.Gaussa w trójk¹cie i wartoœci wag
PGx = zeros(1,Ng1);
PGy = zeros(1,Ng1);
WAG1 = zeros(1,Ng1);
switch Ng1
    case 1
        WAG1 = [1];
    case 3
        WAG1 = [1/3 1/3 1/3];
    case 4
        WAG1 = [-27/48 25/48 25/48 25/48];
    case 7
        WAG1 = [0.225 0.1323941527 0.1323941527 0.1323941527 0.1259391805 0.1259391805 0.1259391805];
    otherwise
        disp('ZLA LICZBA PUNKTÓW GAUSSA!');
        pause
end
switch Ng1
    case 1
        PGx(1) = (1./3.) * (X(w(1)) + X(w(2)) + X(w(3)));
        PGy(1) = (1./3.) * (Y(w(1)) + Y(w(2)) + Y(w(3)));
    case 3
        PGx(1) = 0.5*(X(w(1)) + X(w(2)));
        PGy(1) = 0.5*(Y(w(1)) + Y(w(2)));
        PGx(2) = 0.5*(X(w(2)) + X(w(3)));
        PGy(2) = 0.5*(Y(w(2)) + Y(w(3)));
        PGx(3) = 0.5*(X(w(1)) + X(w(3)));
        PGy(3) = 0.5*(Y(w(1)) + Y(w(3)));
    case 4
        PGx(1) = (X(w(1)) + X(w(2)) + X(w(3)))/3;
        PGy(1) = (Y(w(1)) + Y(w(2)) + Y(w(3)))/3;
        PGx(2) = 0.6*X(w(1)) + 0.2*X(w(2)) + 0.2*X(w(3));
        PGy(2) = 0.6*Y(w(1)) + 0.2*Y(w(2)) + 0.2*Y(w(3));
        PGx(3) = 0.2*X(w(1)) + 0.6*X(w(2)) + 0.2*X(w(3));
        PGy(3) = 0.2*Y(w(1)) + 0.6*Y(w(2)) + 0.2*Y(w(3));
        PGx(4) = 0.2*X(w(1)) + 0.2*X(w(2)) + 0.6*X(w(3));
        PGy(4) = 0.2*Y(w(1)) + 0.2*Y(w(2)) + 0.6*Y(w(3));
    case 7
        a1 = 0.0597158717;
        b1 = 0.4701420641;
        a2 = 0.7974269853;
        b2 = 0.1012865073;
        PGx(1) = (1/3)*(X(w(1)) + X(w(2)) + X(w(3)));
        PGy(1) = (1/3)*(Y(w(1)) + Y(w(2)) + Y(w(3)));
        PGx(2) = a1 * X(w(1)) + b1 * X(w(2)) + b1 * X(w(3));
        PGy(2) = a1 * Y(w(1)) + b1 * Y(w(2)) + b1 * Y(w(3));
        PGx(3) = b1 * X(w(1)) + a1 * X(w(2)) + b1 * X(w(3));
        PGy(3) = b1 * Y(w(1)) + a1 * Y(w(2)) + b1 * Y(w(3));
        PGx(4) = b1 * X(w(1)) + b1 * X(w(2)) + a1 * X(w(3));
        PGy(4) = b1 * Y(w(1)) + b1 * Y(w(2)) + a1 * Y(w(3));
        PGx(5) = a2 * X(w(1)) + b2 * X(w(2)) + b2 * X(w(3));
        PGy(5) = a2 * Y(w(1)) + b2 * Y(w(2)) + b2 * Y(w(3));
        PGx(6) = b2 * X(w(1)) + a2 * X(w(2)) + b2 * X(w(3));
        PGy(6) = b2 * Y(w(1)) + a2 * Y(w(2)) + b2 * Y(w(3));
        PGx(7) = b2 * X(w(1)) + b2 * X(w(2)) + a2 * X(w(3));
        PGy(7) = b2 * Y(w(1)) + b2 * Y(w(2)) + a2 * Y(w(3));
end

function [d,Q]=solveq(K,f,bc)
% a=solveq(K,f)
% [a,Q]=solveq(K,f,bc)
%-------------------------------------------------------------
% PURPOSE
%  Solve static FE-equations considering boundary conditions.
%
% INPUT: K : global stiffness matrix, dim(K)= nd x nd
%        f : global load vector, dim(f)= nd x 1
%
%        bc : boundary condition matrix
%            dim(bc)= nbc x 2, nbc : number of b.c.'s
%
% OUTPUT:  a : solution including boundary values
%          Q : reaction force vector
%              dim(a)=dim(Q)= nd x 1, nd : number of dof's
%-------------------------------------------------------------

% LAST MODIFIED: M Ristinmaa   1993-10-06
% Copyright (c)  Division of Structural Mechanics and
%                Department of Solid Mechanics.
%                Lund Institute of Technology
%-------------------------------------------------------------
if nargin==2 ;
    d=K\f ;
elseif nargin==3;
    [nd,nd]=size(K);
    fdof=[1:nd]';
    %
    d=zeros(size(fdof));
    Q=zeros(size(fdof));
    %
    pdof=bc(:,1);
    dp=bc(:,2);
    fdof(pdof)=[];
    %
    s=K(fdof,fdof)\(f(fdof)-K(fdof,pdof)*dp);
    %
    d(pdof)=dp;
    d(fdof)=s;
end
Q=K*d-f;

%--------------------------end--------------------------------

function [S] = star(N,x,y,X,Y)
S = zeros(N,2);
for i=1:N
    S(i,1) = i;
    S(i,2) = sqrt( (x-X(i))^2 + (y-Y(i))^2 );
end
S = sortrows(S,2);

function [M] = mwls(x,y,X,Y,S,m,p,g)
if nargin<8
    g = 0;
end
if nargin<7
    p = 2;
end
s = 6;
if p==1
    s = 3;
end
if p==1.5
    s = 4;
end
P = zeros(m,s); % local interpolates matrix
W = zeros(m,m); % diagonal weights matrix
for i=1:m
    nr = S(i,1);
    h = X(nr)-x;
    k = Y(nr)-y;
    z = [1 h k 0.5*h^2 h*k 0.5*k^2]; % 2D interpolates up to and including second order
    P(i,:) = z(1:s);
    if p==1.5
        P(i,:) = [z(1:3) z(5)];
    end
    d = sqrt(h^2 + k^2); % distance from the star central point
    W(i,i) = 1/((d^2 + (g^4)/(d^2 + g^2 + 10^-12))^(p+1) + 10^-12); % singular (g=0) or non-singular (g>0) weights
end
M = (P'*W*P)^-1*P'*W; % matrix of the MFD formulas