%AQPMR(P,D,b_max,w_max,w_map,dc,e)
%    Quasi-Polynomial Mapping based Rootfinder utilizing Asymptotic features of the roots chains    
%    Finds all roots of the retarded quasi-polynomial QP(s) located to the left from b_max and within the strip given by the
%    boundaries [0,j*w_max]
%
%    QP(s)=(P(1,1)*s^n+...+P(1,n)*s+P(1,n+1))*exp(-D(1)*s)+
%     +(P(2,1)*s^n+...+P(2,n)*s+P(2,n+1))*exp(-D(2)*s)+
%                    ...
%     +(P(N,1)*s^n+...+P(N,n)*s+P(N,n+1))*exp(-D(N))*s)+
%     +(P(N+1,1)*s^n+...+P(N+1,n)*s+P(N+1,n+1))*exp(-D(N+1)*s)
%
%    located to the left from b_max and within the strip given by the
%    boundaries [0,j*w_max]
%
%    P - matrix of the polynomial coefficients
%    D - row matrix of the delays
%    b_max - upper real boundary of the region to be mapped 
%    w_max - upper imaginary boundary of the region to be mapped
%    w_map - upper imaginary boundary of the region to be mapped whole
%    w_map<w_max
%    dc - width of the strips close to the  to be mapped
%    e - required acuracy
%
%    The function AQPMR requires symbolic math toolbox
%
%    Example
%    find all the roots of the quasipolynomial
%    1.5*s^3+0.2*s^2+20.1-2.1*s*exp(-1.3*s)+3.2*s*exp(-3.5*s)+1.4*exp(-4.3*s)
%    located to the left from b_max=3 and j*w_max=100*j
%
%    P=[1.5 0.2 0 20.1;0 0 -2.1 0;0 0 3.2 0;0 0 0 1.4]
%    D=[0 1.3 3.5 4.3]
%    AQPMR(P,D,3,100,5,0.5,0.000001)
%
%    see also QPMR


function P=aqpmr(Pol,Del,b_max,w_max,w_map,dc,e)
close all
max_order=size(Pol);
max_order=max_order(2);
N_Del=length(Del);
dc=dc/2;

M='(';

for k=1:N_Del
    Mp='(';
    for p=1:max_order
        if Pol(k,p)~=0
            if Mp=='('
               Mp=strcat(Mp,num2str(Pol(k,p)),'*s^',num2str(max_order-p));
            else
               Mp=strcat(Mp,'+',num2str(Pol(k,p)),'*s^',num2str(max_order-p));
            end
        end
    end
       
    Mp=strcat(Mp,')');
    Mp=strrep(Mp,'+-','-');
    Mp=strrep(Mp,'s^1','s');
    Mp=strrep(Mp,'*s^0','');
    
    if M=='('
        M=strcat(Mp,'*exp(-',num2str(Del(k)),'*s)');
    else
        M=strcat(M,'+',Mp,'*exp(-',num2str(Del(k)),'*s)');   
    end
end
M=strrep(M,'*exp(-0*s)','');

%
%body distribucniho diagramu
for k=1:N_Del
    p=1;
    while Pol(k,p)==0
        p=p+1;
    end
    E_pol(1,k)=max_order-p;
    P_pol(1,k)=Pol(k,p);
end     

E_pol;
P_pol;

M_sub='(';

for k=1:N_Del
    if M_sub=='('
           M_sub=strcat(num2str(P_pol(k)),'*s^',num2str(E_pol(k)),'*exp(-',num2str(Del(k)),'*s)');
    else
           M_sub=strcat(M_sub,'+',num2str(P_pol(k)),'*s^',num2str(E_pol(k)),'*exp(-',num2str(Del(k)),'*s)');
    end
end
M_sub=strrep(M_sub,'+-','-');
M_sub=strrep(M_sub,'*exp(-0*s)','');
M_sub=strrep(M_sub,'s^1','s');
M_sub=strrep(M_sub,'*s^0','');

%preskupeni
for k=1:N_Del
    E_pr(N_Del-k+1)=E_pol(k);
    P_pr(N_Del-k+1)=P_pol(k);
end

Del_pr=zeros(1,N_Del);
for k=2:N_Del
    Del_pr(k)=Del(N_Del)-Del(N_Del-k+1);
end

figure(1)
plot(Del_pr,E_pr,'o')


p=1;
n_d=N_Del;
Del_r=Del_pr(p);
E_r=E_pr(1);
P_P=P_pr(1);
E_P=E_pr(1);
r_mi=1;
Del_P(1)=Del_r;
while p~=N_Del
    for k=p+1:n_d
        mi_semi(k-p)=(E_pol(N_Del-k+1)-E_r)/(Del_pr(k)-Del_r);
    end
        [mi(r_mi),b] = max(mi_semi);
        p=b+p;
        E_r=E_pr(p);
        E_P(r_mi+1)=E_r;
        Del_r=Del_pr(p);
        Del_P(r_mi+1)=Del_r;
        P_P(r_mi+1)=P_pr(p);
        r_mi=r_mi+1;
        mi_semi=0;
end

hold on
plot(Del_P,E_P,'k-')

title('Spectrum Distribution diagram')    
xlabel('Delay','FontSize',10)
ylabel('Polynomial degree','FontSize',10)

for k=1:r_mi
    P_Pp(k)=P_P(r_mi-k+1);
end


%vypocet c
k=1;
kp=1;
while k<r_mi
    clear Koef
    
    if k<r_mi-1
        podm=mi(k)~=mi(k+1);
    elseif r_mi==2
        podm=1;
    else
        podm=mi(k)~=mi(k-1);
    end
    
    if podm
        nP=E_P(k+1)-E_P(k);
        Koef(1)=P_P(k+1);
        Koef(nP+1)=P_P(k);
                
        w_P(1:nP)=roots(Koef);
        c(kp)=mi(k)*log(abs(w_P(1)));
        mi_bud(kp)=mi(k);
        k=k+1;
        kp=kp+1;
    else
        kk=k;
        i=k;
        opakuj=(mi(i)==mi(i+1));
        while opakuj
            kk=kk+1;
            if i==r_mi-1;
                opakuj=0;
            else
                opakuj=(mi(i)==mi(i+1));
                i=i+1;
            end
        end
        nP=0;
        Koef(1)=P_P(kk);
        for i=1:(kk-k)
            nP=nP+E_P(kk-i+1)-E_P(kk-i);
            Koef(nP+1)=P_P(kk-i);
        end 
        
        w_P(1:nP)=roots(Koef);
        
        c(kp)=mi(k)*log(abs(w_P(1)));
        mi_bud(kp)=mi(k);
        kp=kp+1;
        
        for i=2:nP
            if abs(w_P(i))~=abs(w_P(i-1))
                c(kp)=mi(k)*log(abs(w_P(i)));
                mi_bud(kp)=mi(k);
                kp=kp+1;
            end
        end
              
        k=kk;
    end   
end

mi=mi_bud;


%samotny vypocet
syms s;
d=pi/Del(N_Del)/3/pi;
b_min=-1.1*(mi(1)*log(w_max));

M=M*exp(Del(N_Del)*s);


%kontrola poctu korenu v oblastech kde nemaji byt
figure(4)
hold on
axis([b_min b_max 0 w_max]);
xlabel('Re(\its \rm)','FontSize',10)
ylabel('Im(\its \rm)','FontSize',10)
figure(3)
axis([b_min b_max 0 w_max]);
hold on



%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
%dosazeny algoritmus z krivky_oblasti_2
k_mi=length(mi);


ds=0.3*d;

w=linspace(1e-10,w_max,1000);

for k=1:k_mi
    
    b=c(k)-mi(k)*log(w);figure(4),plot(b,w,'k-','LineWidth',1), figure(3),plot(b,w,'k-','LineWidth',1);
end

%nalezeni pruseciku s omega max

for k=1:k_mi
    prus_w_max_pom(k)=c(k)-mi(k)*log(w_max);
    prus_w_map_pom(k)=c(k)-mi(k)*log(w_map);
end

[bp_max i_b_max]=max(prus_w_max_pom);
[bp_map i_b_map]=max(prus_w_map_pom);
%serazeni krivek

b_max_k1=b_max;
for k=1:k_mi
    [prus_w_max(k) i_prus_w_max(k)]=min(prus_w_max_pom);
    prus_w_max_pom(i_prus_w_max(k))=bp_max+abs(bp_max);
    b_max=c(k)+3*dc-mi(k)*log(w_map);
    if b_max>b_max_k1
        b_max_k1=b_max;
    end    
end

b_max=b_max_k1;

b_min=c(i_prus_w_max(1))-3*dc-mi(i_prus_w_max(1))*log(w_max);
figure(3),axis([b_min b_max 0 w_max]);
figure(4),axis([b_min b_max 0 w_max]);
Oblasti=zeros(k_mi-1,k_mi+2);
smer_v_Oblasti=zeros(k_mi-1,k_mi+2);

for p=1:k_mi-1
    Oblasti(p,1)=i_prus_w_max(p);
    hledej=1;
    akt_krivka=i_prus_w_max(p);
    smer=-1;
    prus_krivka_k1=w_max;
    prus_krivka_k1_bud=w_max;
    dif_prus_krivka=w_max;
    smer_v_Oblasti(p,1)=smer;
    i_k=2;
    
    while hledej
        for k=1:k_mi
            if k~=akt_krivka
                prus_krivka=exp((c(akt_krivka)-c(k))/(mi(akt_krivka)-mi(k)));
                if (smer==-1)&((prus_krivka_k1-prus_krivka)<dif_prus_krivka)&(prus_krivka<prus_krivka_k1);
                    dif_prus_krivka=prus_krivka_k1-prus_krivka;
                    i_prus_krivka=k;
                    prus_krivka_k1_bud=prus_krivka;
                end
                if (smer==1)&((prus_krivka-prus_krivka_k1)<dif_prus_krivka)&(prus_krivka>prus_krivka_k1);
                    dif_prus_krivka=prus_krivka-prus_krivka_k1;
                    i_prus_krivka=k;
                    prus_krivka_k1_bud=prus_krivka;
                end
            end
        end
        
        %kdyz se neprotne se zadnou krivkou na intervalu nebo se protne
        %pod w_map
    
        if (prus_krivka_k1_bud==prus_krivka_k1)|prus_krivka_k1_bud<w_map
             Oblasti(p,i_k)=k_mi+1;
             smer=1;
             i_k=i_k+1;
             prus_krivka_k1=2*abs(bp_map);
             prus_krivka_akt=c(akt_krivka)-mi(akt_krivka)*log(w_map);
                 for k=1:k_mi
                      if k~=akt_krivka
                          prus_krivka=c(k)-mi(k)*log(w_map);
                          if (prus_krivka<prus_krivka_k1)&(prus_krivka>prus_krivka_akt)
                                i_prus_krivka=k;
                                prus_krivka_k1=prus_krivka;
                          end
                      end
                 end
             prus_krivka_k1=w_map;
             Oblasti(p,i_k)=i_prus_krivka;
             smer_v_Oblasti(p,i_k)=1;
             akt_krivka=i_prus_krivka;
             dif_prus_krivka=w_max-w_map; 
        else
        
       Oblasti(p,i_k)=i_prus_krivka;
       akt_krivka=i_prus_krivka;
       prus_krivka_k1=prus_krivka_k1_bud;
       
       if find(i_prus_w_max==i_prus_krivka)<p
           smer=-1;
           dif_prus_krivka=prus_krivka_k1_bud;
       end
       
       if find(i_prus_w_max==i_prus_krivka)>p
           smer=1;
           dif_prus_krivka=w_max-prus_krivka_k1_bud;
       end
       
       smer_v_Oblasti(p,i_k)=smer;
       
     end
     
       if i_prus_krivka==i_prus_w_max(p+1);
           hledej=0;
           posl_krivka(p)=i_k;
       end
       
       i_k=i_k+1;
           
    end
    end
    
  %uprava oblasti
  
    Oblasti_uprav=zeros(k_mi+1,k_mi+2);
    smer_v_Oblasti_uprav=zeros(k_mi+1,k_mi+2);
    
    if k_mi==1
        Oblasti_uprav(:,1:2)=[1 2;1 2];
        smer_v_Oblasti_uprav(:,1)=[-1;-1]
    else
      
    for obl=1:k_mi-1
           
    if c(Oblasti(obl,1))+dc+3*d-mi(Oblasti(obl,1))*log(w_max)<c(Oblasti(obl,posl_krivka(obl)))-dc-3*d-mi(Oblasti(obl,posl_krivka(obl)))*log(w_max)
    krivka_uprav=1;
    Oblasti_uprav(obl,krivka_uprav)=Oblasti(obl,1);
    smer_v_Oblasti_uprav(obl,krivka_uprav)=smer_v_Oblasti(obl,1);
    krivka=1;
    %b=c(Oblasti(obl,krivka))+dc-mi(Oblasti(obl,krivka))*log(w);plot(b,w,'r-','LineWidth',1);
    
    while krivka~=posl_krivka(obl)
    k=krivka+1;
    if smer_v_Oblasti(obl,krivka)==-1
        prus_k1=0;
    else
        prus_k1=w_max;
    end
      
    if Oblasti(obl,krivka)==k_mi+1
        Oblasti_uprav(obl,krivka)=k_mi+1;
        Oblasti_uprav(obl,krivka+1)=Oblasti(obl,krivka+1);
        smer_v_Oblasti_uprav(obl,krivka)=0;
        smer_v_Oblasti_uprav(obl,krivka+1)=1;
        krivka=krivka+1;            
        krivka_uprav=krivka_uprav+1;
    else
        
    while Oblasti(obl,k)~=0
        
        if Oblasti(obl,k)==k_mi+1;
            prus=w_map;
        else
            prus=exp((c(Oblasti(obl,krivka))-smer_v_Oblasti(obl,krivka)*dc-c(Oblasti(obl,k))+smer_v_Oblasti(obl,k)*dc)/(mi(Oblasti(obl,krivka))-mi(Oblasti(obl,k))));
            %b=c(Oblasti(obl,k))-dc*smer_v_Oblasti(obl,k)-mi(Oblasti(obl,k))*log(w);plot(b,w,'r-','LineWidth',1);
        end
        if (prus>prus_k1)&(prus<=w_max)&smer_v_Oblasti(obl,krivka)==-1
            nasl_krivka=k;
            prus_k1=prus;
        end  
        if (prus<prus_k1)&(prus<=w_max)&smer_v_Oblasti(obl,krivka)==1
            nasl_krivka=k;
            prus_k1=prus;
        end 
        k=k+1;
    end
    krivka_uprav=krivka_uprav+1;
    
    if prus_k1>=w_map
        krivka=nasl_krivka;
        Oblasti_uprav(obl,krivka_uprav)=Oblasti(obl,krivka);
        smer_v_Oblasti_uprav(obl,krivka_uprav)=smer_v_Oblasti(obl,krivka);
    else
        Oblasti_uprav(obl,krivka_uprav)=k_mi+1;
        
        k=krivka+1;
        prus_k1=b_max;
        while Oblasti(obl,k)~=0 
            prus=c(Oblasti(obl,k))-smer_v_Oblasti(obl,k)*dc-mi(Oblasti(obl,k))*log(w_map);
            if (prus<prus_k1)&(smer_v_Oblasti(obl,k)==1)
            nasl_krivka=k;
            prus_k1=prus;
            end         
        k=k+1;
        end
        krivka=nasl_krivka;
        krivka_uprav=krivka_uprav+1;
        Oblasti_uprav(obl,krivka_uprav)=Oblasti(obl,krivka);
        smer_v_Oblasti_uprav(obl,krivka_uprav)=smer_v_Oblasti(obl,krivka);
    end
    end
    end
    end
    end
    
    %oblasti nejvic nalevo a napravo
    for kk=1:2
    obl=obl+1;
    i_k=1;
    
    if kk==2
        akt_krivka=i_prus_w_max(1);
    else
        akt_krivka=i_prus_w_max(length(i_prus_w_max));
    end
    
    Oblasti_uprav(obl,i_k)=akt_krivka;
        
    prus_krivka_k1=w_max;
    dif_prus_krivka=w_max;
    smer_v_Oblasti_uprav(obl,i_k)=-1;
    hledej=1;
    
    while hledej
        hledej=0;
        for k=1:k_mi
            if k~=akt_krivka
                prus_krivka=exp((c(akt_krivka)-c(k))/(mi(akt_krivka)-mi(k)));
                if ((prus_krivka_k1-prus_krivka)<dif_prus_krivka)&(prus_krivka_k1-prus_krivka)>0
                    dif_prus_krivka=prus_krivka_k1-prus_krivka;
                    i_prus_krivka=k;
                    prus_krivka_k1_bud=prus_krivka;
                    hledej=1;
                end          
            end
        end
       
        prus_krivka_k1=prus_krivka_k1_bud;
        akt_krivka=i_prus_krivka;
        i_k=i_k+1;
        
       if  (prus_krivka_k1>w_map)&(hledej==1)
           Oblasti_uprav(obl,i_k)=i_prus_krivka;
           smer_v_Oblasti_uprav(obl,i_k)=-1;
       else
           Oblasti_uprav(obl,i_k)=k_mi+1;
           hledej=0;
       end
       dif_prus_krivka=w_max;
    end
    end
    end
    
    %uzavreni oblasti kde nemaji byt koreny
    hran_b=zeros(k_mi+1,2*k_mi+4);
    hran_w=zeros(k_mi+1,2*k_mi+4);
    
    for obl=1:k_mi
    if Oblasti_uprav(obl,1)~=0
    clear hranice
    k=2;
    prus_w_k1=w_max;
    prus_b_k1=c(Oblasti_uprav(obl,1))-smer_v_Oblasti_uprav(obl,1)*dc-mi(Oblasti_uprav(obl,1))*log(w_max);
    hran_w(obl+1,1)=w_max;
    hran_b(obl+1,1)=prus_b_k1;
    n_hr_k1=1;
    n_hr=1;
    
    while Oblasti_uprav(obl,k)~=0
    
    if Oblasti_uprav(obl,k)==k_mi+1
        prus_w=w_map;
        prus_b=c(Oblasti_uprav(obl,k-1))-smer_v_Oblasti_uprav(obl,k-1)*dc-mi(Oblasti_uprav(obl,k-1))*log(w_map);
            if abs(prus_w-prus_w_k1)>abs(prus_b-prus_b_k1)
                Ns=round(abs(prus_w-prus_w_k1)/ds);
                n_hr=n_hr+Ns;
                hr_w=linspace(prus_w_k1,prus_w,Ns);
                hr_b=c(Oblasti_uprav(obl,k-1))-smer_v_Oblasti_uprav(obl,k-1)*dc-mi(Oblasti_uprav(obl,k-1))*log(hr_w);
            else
                Ns=round(abs(prus_b-prus_b_k1)/ds);
                n_hr=n_hr+Ns;
                hr_b=linspace(prus_b_k1,prus_b,Ns);
                hr_w=exp((c(Oblasti_uprav(obl,k-1))-smer_v_Oblasti_uprav(obl,k-1)*dc-hr_b)/mi(Oblasti_uprav(obl,k-1)));
            end
    elseif Oblasti_uprav(obl,k-1)==k_mi+1
        prus_w=w_map;
        prus_b=c(Oblasti_uprav(obl,k))-smer_v_Oblasti_uprav(obl,k)*dc-mi(Oblasti_uprav(obl,k))*log(w_map);
        Ns=round(abs(prus_b-prus_b_k1)/ds);
        n_hr=n_hr+Ns;
        hr_b=linspace(prus_b_k1,prus_b,Ns);
        hr_w=w_map;
    else
        prus_w=exp((c(Oblasti_uprav(obl,k))-smer_v_Oblasti_uprav(obl,k)*dc-c(Oblasti_uprav(obl,k-1))+smer_v_Oblasti_uprav(obl,k-1)*dc)/(mi(Oblasti_uprav(obl,k))-mi(Oblasti_uprav(obl,k-1))));
        prus_b=c(Oblasti_uprav(obl,k))-smer_v_Oblasti_uprav(obl,k)*dc-mi(Oblasti_uprav(obl,k))*log(prus_w);
            if abs(prus_w-prus_w_k1)>abs(prus_b-prus_b_k1);
                Ns=round(abs(prus_w-prus_w_k1)/ds);
                n_hr=n_hr+Ns;
                hr_w=linspace(prus_w_k1,prus_w,Ns);
                hr_b=c(Oblasti_uprav(obl,k-1))-smer_v_Oblasti_uprav(obl,k-1)*dc-mi(Oblasti_uprav(obl,k-1))*log(hr_w);
            else
                Ns=round(abs(prus_b-prus_b_k1)/ds);
                n_hr=n_hr+Ns;
                hr_b=linspace(prus_b_k1,prus_b,Ns);
                hr_w=exp((c(Oblasti_uprav(obl,k-1))-smer_v_Oblasti_uprav(obl,k-1)*dc-hr_b)/mi(Oblasti_uprav(obl,k-1)));
            end
    end
    
    hran_w(obl+1,k)=prus_w;
    hran_b(obl+1,k)=prus_b;
    
    hranice(1,n_hr_k1:n_hr-1)=hr_b;
    hranice(2,n_hr_k1:n_hr-1)=hr_w;
    
    k=k+1;
    n_hr_k1=n_hr;
    prus_b_k1=prus_b;
    prus_w_k1=prus_w;
    end
    
    if obl==k_mi %oblast vpravo
            prus_w=w_map;
            prus_b=b_max;
            Ns=round(abs(prus_b-prus_b_k1)/ds);
            n_hr=n_hr+Ns;
            hranice(1,n_hr_k1:n_hr-1)=linspace(prus_b_k1,prus_b,Ns);
            hranice(2,n_hr_k1:n_hr-1)=w_map;
            hran_w(obl+1,k)=prus_w;
            hran_b(obl+1,k)=prus_b;
            n_hr_k1=n_hr;
            
            prus_w_k1=prus_w;
            prus_w=w_max;
            Ns=round(abs(prus_w-prus_w_k1)/ds);
            n_hr=n_hr+Ns;
            hranice(1,n_hr_k1:n_hr-1)=b_max;
            hranice(2,n_hr_k1:n_hr-1)=linspace(prus_w_k1,prus_w,Ns);
            hran_w(obl+1,k+1)=prus_w;
            hran_b(obl+1,k+1)=b_max;
            n_hr_k1=n_hr;
            
            prus_b_k1=prus_b;
            prus_b=c(Oblasti_uprav(obl,1))-smer_v_Oblasti_uprav(obl,1)*dc-mi(Oblasti_uprav(obl,1))*log(w_max);
            Ns=round(abs(prus_b-prus_b_k1)/ds);
            n_hr=n_hr+Ns;
            hr_b=linspace(prus_b_k1,prus_b,Ns);
            hranice(1,n_hr_k1:n_hr-1)=linspace(prus_b_k1,prus_b,Ns);
            hranice(2,n_hr_k1:n_hr-1)=w_max;
            hran_w(obl+1,k+2)=w_max;
            hran_b(obl+1,k+2)=prus_b;
            figure(4)
            d_arg=argOblast(M,hranice);
            text((c(Oblasti_uprav(obl,1))-mi(Oblasti_uprav(obl,1))*log(w_max)+b_max)/2,0.95*w_max,num2str(d_arg));
            
              
    else
       
    prus_w=w_max;
    prus_b=c(Oblasti_uprav(obl,k-1))-smer_v_Oblasti_uprav(obl,k-1)*dc-mi(Oblasti_uprav(obl,k-1))*log(w_max);
    Ns=round(abs(prus_w-prus_w_k1)/ds);
    n_hr=n_hr+Ns;
    hr_w=linspace(prus_w_k1,prus_w,Ns);
    hranice(1,n_hr_k1:n_hr-1)=c(Oblasti_uprav(obl,k-1))-smer_v_Oblasti_uprav(obl,k-1)*dc-mi(Oblasti_uprav(obl,k-1))*log(hr_w);
    hranice(2,n_hr_k1:n_hr-1)=hr_w;
    hran_w(obl+1,k)=prus_w;
    hran_b(obl+1,k)=prus_b;
    
    
    n_hr_k1=n_hr;
    prus_b_k1=prus_b;
    prus_b=c(Oblasti_uprav(obl,1))-smer_v_Oblasti_uprav(obl,1)*dc-mi(Oblasti_uprav(obl,1))*log(w_max);
    Ns=round(abs(prus_b-prus_b_k1)/ds);
    n_hr=n_hr+Ns;
    hranice(1,n_hr_k1:n_hr-1)=linspace(prus_b_k1,prus_b,Ns);
    hranice(2,n_hr_k1:n_hr-1)=w_max;
    hran_w(obl+1,k+1)=w_max;
    hran_b(obl+1,k+1)=prus_b;
    
    
    d_arg=argOblast(M,hranice);
    figure(4)
    text((c(Oblasti_uprav(obl,1))-mi(Oblasti_uprav(obl,1))*log(w_max)+c(Oblasti(obl,posl_krivka(obl)))-mi(Oblasti(obl,posl_krivka(obl)))*log(w_max))/2,0.95*w_max,num2str(d_arg));
    end
    figure(4)
    plot(hranice(1,:),hranice(2,:),'k--')
    figure(3)
    plot(hranice(1,:),hranice(2,:),'g') 
    end
    end
    
    %oblast vlevo
    
    obl=obl+1;
    clear hranice
    k=2;
    prus_w_k1=w_max;
    prus_b_k1=c(Oblasti_uprav(obl,1))-dc-mi(Oblasti_uprav(obl,1))*log(w_max);
    n_hr_k1=1;
    n_hr=1;
    hran_w(1,1)=w_max;
    hran_b(1,1)=prus_b_k1;
    
    while Oblasti_uprav(obl,k)~=0
        if Oblasti_uprav(obl,k)==k_mi+1
            prus_w=w_map;
            prus_b=c(Oblasti_uprav(obl,k-1))-dc-mi(Oblasti_uprav(obl,k-1))*log(prus_w);
            prus_b_dole=prus_b;
        else
           prus_w=exp((c(Oblasti_uprav(obl,k))-dc-c(Oblasti_uprav(obl,k-1))+dc)/(mi(Oblasti_uprav(obl,k))-mi(Oblasti_uprav(obl,k-1))));
           prus_b=c(Oblasti_uprav(obl,k))-dc-mi(Oblasti_uprav(obl,k))*log(prus_w);
        end
        
            if abs(prus_w-prus_w_k1)>abs(prus_b-prus_b_k1);
                Ns=round(abs(prus_w-prus_w_k1)/ds);
                n_hr=n_hr+Ns;
                hr_w=linspace(prus_w_k1,prus_w,Ns);
                hr_b=c(Oblasti_uprav(obl,k-1))-dc-mi(Oblasti_uprav(obl,k-1))*log(hr_w);
            else
                Ns=round(abs(prus_b-prus_b_k1)/ds);
                n_hr=n_hr+Ns;
                hr_b=linspace(prus_b_k1,prus_b,Ns);
                hr_w=exp((c(Oblasti_uprav(obl,k-1))-dc-hr_b)/mi(Oblasti_uprav(obl,k-1)));
            end
                     
           hranice(1,n_hr_k1:n_hr-1)=hr_b;
           hranice(2,n_hr_k1:n_hr-1)=hr_w;
           hran_w(1,k)=prus_w;
           hran_b(1,k)=prus_b;
           k=k+1;
           n_hr_k1=n_hr;
           prus_b_k1=prus_b;
           prus_w_k1=prus_w;
    end
                
            prus_w=w_map;
            prus_b=b_min;
            Ns=round(abs(prus_b-prus_b_k1)/ds);
            n_hr=n_hr+Ns;
            hranice(1,n_hr_k1:n_hr-1)=linspace(prus_b_k1,prus_b,Ns);
            hranice(2,n_hr_k1:n_hr-1)=w_map;
            hran_w(1,k)=prus_w;
            hran_b(1,k)=prus_b;
            n_hr_k1=n_hr;
            
            prus_w_k1=prus_w;
            prus_w=w_max;
            Ns=round(abs(prus_w-prus_w_k1)/ds);
            n_hr=n_hr+Ns;
            hranice(1,n_hr_k1:n_hr-1)=b_min;
            hranice(2,n_hr_k1:n_hr-1)=linspace(prus_w_k1,prus_w,Ns);
            hran_w(1,k+1)=prus_w;
            hran_b(1,k+1)=b_min;
            n_hr_k1=n_hr;
            
            prus_b_k1=prus_b;
            prus_b=c(Oblasti_uprav(obl,1))-dc-mi(Oblasti_uprav(obl,1))*log(w_max);
            Ns=round(abs(prus_b-prus_b_k1)/ds);
            n_hr=n_hr+Ns;
            hr_b=linspace(prus_b_k1,prus_b,Ns);
            hranice(1,n_hr_k1:n_hr-1)=linspace(prus_b_k1,prus_b,Ns);
            hranice(2,n_hr_k1:n_hr-1)=w_max;
            hran_w(1,k+2)=w_max;
            hran_b(1,k+2)=prus_b;
            
            d_arg=argOblast(M,hranice);
            figure(4)
            text((c(Oblasti_uprav(obl,1))-mi(Oblasti_uprav(obl,1))*log(w_max)+b_min)/2,0.95*w_max,num2str(-1*d_arg));
                figure(4)
                plot(hranice(1,:),hranice(2,:),'k--')
                figure(3)
                plot(hranice(1,:),hranice(2,:),'g')
            
 
   
%Oblast vlevo dole
obl=obl+1;
    clear hranice
    n_hr_k1=1;
    n_hr=1;
    prus_w_k1=w_map;
    prus_b_k1=b_min;
    
    prus_w=w_map;
    prus_b=prus_b_dole;
    Ns=round(abs(prus_b-prus_b_k1)/ds);
    n_hr=n_hr+Ns;
    hranice(1,n_hr_k1:n_hr-1)=linspace(prus_b_k1,prus_b,Ns);
    hranice(2,n_hr_k1:n_hr-1)=w_map;
    n_hr_k1=n_hr;
    
    prus_w_k1=prus_w;
    prus_w=-d;
    Ns=round(abs(prus_w-prus_w_k1)/ds);
    n_hr=n_hr+Ns;
    hranice(1,n_hr_k1:n_hr-1)=prus_b_dole;
    hranice(2,n_hr_k1:n_hr-1)=linspace(prus_w_k1,prus_w,Ns);
    n_hr_k1=n_hr;
    
    prus_b_k1=prus_b_dole;
    prus_b=b_min;
    Ns=round(abs(prus_b-prus_b_k1)/ds);
    n_hr=n_hr+Ns;
    hranice(1,n_hr_k1:n_hr-1)=linspace(prus_b_k1,prus_b,Ns);
    hranice(2,n_hr_k1:n_hr-1)=-d;
    n_hr_k1=n_hr;
    
    prus_w_k1=-d;
    prus_w=w_map;
    Ns=round(abs(prus_w-prus_w_k1)/ds);
    n_hr=n_hr+Ns;
    hranice(1,n_hr_k1:n_hr-1)=b_min;
    hranice(2,n_hr_k1:n_hr-1)=linspace(prus_w_k1,prus_w,Ns);
    n_hr_k1=n_hr;
    
    d_arg=argOblast(M,hranice);
    figure(4)
    text((b_min+prus_b_dole)/2,0.8*w_map,num2str(-1*d_arg))
    figure(4)
    plot(hranice(1,:),hranice(2,:),'k--')


    %XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX


%Mapovaci algoritmus

pokracovat=input('for mapping press y   ','s');

if pokracovat=='y'
%vypocet
    figure(3)
    plot(hranice(1,:),hranice(2,:),'g')

b_map=prus_b_dole;

kb_max=floor((b_max+3*d-b_min+3*d)/d)+1;
kb_map=floor((b_map-b_min)/d+1);
kw_min=floor((w_map)/d)+3;
kw_max=floor((w_max+3*d)/d)+1;

Mm=zeros(kw_max,kb_max);
%Maska_P=zeros(kw_max,kb_max);

Mm(1:(kw_min),kb_map:kb_max)=subs(M,s,ones(1,kw_min)'*linspace(b_map-3*d,b_max+3*d,kb_max-kb_map+1)+j*((0:kw_min-1)'*d)*ones(1,kb_max-kb_map+1));
%Maska_P(1:kw_min,kb_map:kb_max)=ones(kw_min,kb_max-kb_map+1);

%Presun oblasti a hranic tak aby sly za sebou

pom1=Oblasti_uprav(k_mi+1,:);
pom2=smer_v_Oblasti_uprav(k_mi+1,:);
for k=1:k_mi
    Oblasti_uprav(k_mi+2-k,:)=Oblasti_uprav(k_mi+1-k,:);
    smer_v_Oblasti_uprav(k_mi+2-k,:)=smer_v_Oblasti_uprav(k_mi+1-k,:);
end
    Oblasti_uprav(1,:)=pom1;
    smer_v_Oblasti_uprav(1,:)=pom2;
      
%vyber max min v oblasti

hran_b_hran=zeros(k_mi+1,3);
i_hran_b_hran=zeros(k_mi+1,3);

       
if k_mi==1
        hran_b_hran(1,1)=b_min;
        hran_b_hran(1,2)=hran_b(1,1);
        hran_b_hran(1,3)=hran_b(1,2);
    
        hran_b_hran(2,1)=hran_b(2,1);
        hran_b_hran(2,2)=hran_b(2,2);
        hran_b_hran(2,3)=b_max;
else
    
        hran_b_hran(1,1)=b_min;
        hran_b_hran(1,2)=hran_b(1,1);
        hran_b_hran(1,3)=max(hran_b(1,1:min(find(hran_b(1,:)==0))-1));

for k=2:k_mi
    if hran_b(k,1)~=0
         i_min_b=min(find(hran_b(k,:)==0));
         hran_b_hran(k,1)=min(hran_b(k,1:i_min_b-1));
         hran_b_hran(k,2)=c(Oblasti(k-1,posl_krivka(k-1)))-dc-mi((Oblasti(k-1,posl_krivka(k-1))))*log(w_max);
         hran_b_hran(k,3)=max(hran_b(k,1:min(find(hran_b(k,:)==0))-1));
    end
end
end

    if hran_b(k+1,1)~=0
        i_min_b=min(find(hran_b(k+1,:)==0));
        hran_b_hran(k+1,1)=min(hran_b(k+1,1:i_min_b-1));
        hran_b_hran(k+1,2)=b_max;
        hran_b_hran(k+1,3)=b_max;      
    end
    
%urceni indexu hranicnich bodu
i_bmax_prus=zeros(k_mi+1,1);
i_bmax_prus(1)=Oblasti_uprav(1,find(Oblasti_uprav(1,:)==k_mi+1)-1);
for k=2:k_mi
    if smer_v_Oblasti_uprav(k,1)~=0
        i_bmax_prus(k)=Oblasti_uprav(k,min(find(smer_v_Oblasti_uprav(k,:)==1))-1);
    end
end


for k_it=1:kb_max
 bs=b_min+k_it*d;
 prus_obl=zeros(k_mi+1,1);
 %plot([bs bs],[0 w_max],'b:')
 kk=1;
 
 %prvni oblast - ta nejvic nalevo
 
 for k=1:k_mi+1
     if bs>hran_b_hran(k,1)&bs<hran_b_hran(k,3)
         prus_obl(kk)=k;
         kk=kk+1;
     end
 end
 
 if bs<hran_b_hran(1,3)
 k=2;
 prus_k1=hran_b_hran(1,3)-bs;
 wb_min=exp((c(i_bmax_prus(prus_obl(1)))-dc-bs)/mi(i_bmax_prus(prus_obl(1))));
 if prus_obl(1)==1
     while hran_b(1,k)~=hran_b_hran(1,3)
         prus=hran_b(1,k)-bs;
         if (prus<prus_k1)&(prus>0)
             wb_min=exp((c(Oblasti_uprav(1,k-1))-dc-bs)/mi(Oblasti_uprav(1,k-1)));
             prus_k1=prus;
         end
         k=k+1;
     end
     if prus_obl(2)==0
         wb_max=w_max;
     else
         k=2;
         prus_k1=hran_b_hran(prus_obl(2),3)-bs;
         if i_bmax_prus(prus_obl(2))==k_mi+1|i_bmax_prus(prus_obl(2))==0
             wb_max=w_map;
         else
             wb_max=exp((c(i_bmax_prus(prus_obl(2)))+dc-bs)/mi(i_bmax_prus(prus_obl(2))));
         end
         
         while hran_b(prus_obl(2),k)~=hran_b_hran(prus_obl(2),3)
             prus=hran_b(prus_obl(2),k)-bs;
             if (prus<prus_k1)&(prus>0)
                 wb_max=exp((c(Oblasti_uprav(prus_obl(2),k-1))+dc-bs)/mi(Oblasti_uprav(prus_obl(2),k-1)));
                 prus_k1=prus;
             end
             k=k+1;
         end
        
     end
 end
 %plot([bs bs],[wb_min wb_max],'b-')
 
 ky_h=floor((wb_max)/d)+3;
 ky_d=floor((wb_min)/d)-2;    
 Mm(ky_d:ky_h,k_it+4)=subs(M,s,bs+j*((ky_d:ky_h)'*d));
 %Maska_P(ky_d+3:ky_h-4,k_it+4)=ones(ky_h-ky_d-6,1);
 
 if prus_obl(2)~=0 
 %nasledujici oblasti
 i_obl=2;
 w_b_max=0;
 k=1;
      while ~((hran_b_hran(prus_obl(i_obl),1)<bs)&(bs<hran_b_hran(prus_obl(i_obl),2)))&wb_max<w_max   
          
      while ~((hran_b(prus_obl(i_obl),k)>bs)&(bs>hran_b(prus_obl(i_obl),k+1)))
         k=k+1;
      end
      wb_min=exp((c(Oblasti_uprav(prus_obl(i_obl),k))-dc-bs)/mi(Oblasti_uprav(prus_obl(i_obl),k)));
      
      if prus_obl(i_obl+1)==0
          wb_max=w_max;
      else
          i_obl=i_obl+1;
          k=1;
          while ~((hran_b(prus_obl(i_obl),k)<bs)&(bs<hran_b(prus_obl(i_obl),k+1)))
                k=k+1;
          end
          wb_max=exp((c(Oblasti_uprav(prus_obl(i_obl),k))+dc-bs)/mi(Oblasti_uprav(prus_obl(i_obl),k)));
          if wb_max>w_max;
              wb_max=w_max;
          end
      end          
          
      %plot([bs bs],[wb_min wb_max],'b-') 
      
      ky_h=floor((wb_max)/d)+3;
      ky_d=floor((wb_min)/d)-2;    
      Mm(ky_d:ky_h,k_it+4)=subs(M,s,bs+j*((ky_d:ky_h)'*d));
      %Maska_P(ky_d+3:ky_h-4,k_it+4)=ones(ky_h-ky_d-6,1);
 end 
 end
 else
     
     wb_min=w_map;
     if prus_obl(1)==0
         if bs<b_max
         wb_max=w_max;
         %plot([bs bs],[wb_min wb_max],'b-')
         
         ky_h=floor((wb_max)/d)+3;
         ky_d=floor((wb_min)/d)-2;    
         Mm(ky_d:ky_h,k_it+4)=subs(M,s,bs+j*((ky_d:ky_h)'*d));
         %Maska_P(ky_d+3:ky_h-4,k_it+4)=ones(ky_h-ky_d-6,1);
         end
         
     else
         k=1;
         while ~((hran_b(prus_obl(1),k)<bs)&(bs<hran_b(prus_obl(1),k+1)))
              k=k+1;
         end
         if Oblasti_uprav(prus_obl(1),k)~=k_mi+1
              wb_max=exp((c(Oblasti_uprav(prus_obl(1),k))+dc-bs)/mi(Oblasti_uprav(prus_obl(1),k)));
              %plot([bs bs],[wb_min wb_max],'b-')
              
             ky_h=floor((wb_max)/d)+3;
             ky_d=floor((wb_min)/d)-2;    
             Mm(ky_d:ky_h,k_it+4)=subs(M,s,bs+j*((ky_d:ky_h)'*d));
             %Maska_P(ky_d+3:ky_h-4,k_it+4)=ones(ky_h-ky_d-6,1);
             
         end
     
     if bs>hran_b_hran(prus_obl(1),2)&prus_obl(2)==0
         k=1;
         while ~((hran_b(prus_obl(1),k)>bs)&(bs>hran_b(prus_obl(1),k+1)))
              k=k+1;
         end
         wb_min=exp((c(Oblasti_uprav(prus_obl(1),k))-dc-bs)/mi(Oblasti_uprav(prus_obl(1),k)));
         wb_max=w_max;
         %plot([bs bs],[wb_min wb_max],'b-')
         
         ky_h=floor((wb_max)/d)+3;
         ky_d=floor((wb_min)/d)-2;    
         Mm(ky_d:ky_h,k_it+4)=subs(M,s,bs+j*((ky_d:ky_h)'*d));
         %Maska_P(ky_d+3:ky_h-4,k_it+4)=ones(ky_h-ky_d-6,1);
         
     end 
     end
     
     
     %nasledujici oblasti
     if prus_obl(2)~=0 
     i_obl=1;
     wb_max=0;
     k=1;
      while ~((hran_b_hran(prus_obl(i_obl),1)<bs)&(bs<hran_b_hran(prus_obl(i_obl),2)))&wb_max<w_max   
          
      while ~((hran_b(prus_obl(i_obl),k)>bs)&(bs>hran_b(prus_obl(i_obl),k+1)))
         k=k+1;
      end
      wb_min=exp((c(Oblasti_uprav(prus_obl(i_obl),k))-dc-bs)/mi(Oblasti_uprav(prus_obl(i_obl),k)));
      
      if prus_obl(i_obl+1)==0
          wb_max=w_max;
      else
          i_obl=i_obl+1;
          k=1;
          while ~((hran_b(prus_obl(i_obl),k)<bs)&(bs<hran_b(prus_obl(i_obl),k+1)))
                k=k+1;
          end
          wb_max=exp((c(Oblasti_uprav(prus_obl(i_obl),k))+dc-bs)/mi(Oblasti_uprav(prus_obl(i_obl),k)));
          if wb_max>w_max;
              wb_max=w_max;
          end
      end          
      
      
      ky_h=floor((wb_max)/d)+3;
      ky_d=floor((wb_min)/d)-2;    
      Mm(ky_d:ky_h,k_it+4)=subs(M,s,bs+j*((ky_d:ky_h)'*d));
      %Maska_P(ky_d+3:ky_h-4,k_it+4)=ones(ky_h-ky_d-6,1);
      
      end      
          
      end     
 end
end

%otupeni hran
for k=2:k_mi
    if smer_v_Oblasti_uprav(k,min(find(smer_v_Oblasti_uprav(k,:)==1))-1)==-1
        distance=0;
        b_leva=hran_b_hran(k,3);
        b_prava=b_leva;
        i_leva=min(find(smer_v_Oblasti_uprav(k,:)==1))-1;
        i_prava=i_leva+1;
        leva=Oblasti_uprav(k,i_leva);
        prava=Oblasti_uprav(k,i_prava);
        i_hran_leva=find(hran_b(k,:)==hran_b_hran(k,3))-1;
        i_hran_prava=i_leva+2;
        wb_max=exp((c(leva)+dc-b_leva)/mi(leva))
                
        while distance<4*d
            wb_max=wb_max+d;
            b_prava=c(prava)-dc-mi(prava)*log(wb_max);
            b_leva=c(leva)+dc-mi(leva)*log(wb_max);
            distance=b_prava-b_leva;
            %figure(4)
            %plot([b_leva b_prava],[wb_max wb_max])
            
            kx_l=floor((b_leva-b_min)/d)+4;
            kx_p=floor((b_prava-b_min)/d)+4;
            ky_h=floor((wb_max)/d);
            Mm(ky_h,kx_l:kx_p)=subs(M,s,b_min+(kx_l:kx_p)*d-3*d+j*ky_h*d);
            %Maska_P(ky_h,kx_l:kx_p)=ones(1,kx_p-kx_l+1);
            
            
            if b_leva<hran_b(k,i_hran_leva)
                i_hran_leva=i_hran_leva-1;
                i_leva=i_leva-1;
                leva=Oblasti_uprav(k,i_leva);
            end
            if b_prava<hran_b(k,i_hran_prava)
                i_hran_prava=i_hran_prava+1;
                i_prava=i_prava+1;
                prava=Oblasti_uprav(k,i_prava);
            end            
        end           
    end
end


'Mapping the contours'

%zacatek vypoctu mapovacim algoritmem
MM=double(Mm);
R=real(MM);
I=imag(MM);
W=linspace(0,w_max+3*d,kw_max)';
B=linspace(b_min-3*d,b_max+3*d,kb_max)';
figure(3)
hold on
[Ci,Hi]=contour(B,W,I,[0 0],'r-');
[Cr,Hr]=contour(B,W,R,[0 0],'b-');

axis([b_min b_max 0 w_max]);
 

xlabel('Re(\its \rm)','FontSize',10)
ylabel('Im(\its \rm)','FontSize',10)
            
pCr=size(Cr);     %kontrola po sobe jdoucich hodnot
pCr=pCr(2);
Crr=Cr;
Cr=Cr(1,:)+Cr(2,:)*j;
IIc=imag(subs(M,s,Cr));

'Locating the intersection points'
np=1;
ch=1;
for no=2:pCr-2
    if ((IIc(no-1)*IIc(no+1))<0)&(ch==1)
        if (abs(Crr(1,no)-Crr(1,no+1))<2*d)&(abs(Crr(2,no)-Crr(2,no+1))<2*d)
               Pp(np)=Cr(no);
               np=np+1;
               ch=0;
        end
    else
        ch=1;
    end
end


%selekce polu ktere jsou mimo
kk=1;
for k=1:length(Pp)
    if imag(Pp(k))>w_map&floor((real(Pp(k)-b_min)/d))+5<kb_max
        %if Maska_P(floor(imag(Pp(k))/d),floor((real(Pp(k)-b_min)/d))+4)~=0
         %if Maska_P(floor(imag(Pp(k))/d),floor((real(Pp(k)-b_min)/d))+3)~=0&Maska_P(floor(imag(Pp(k))/d),floor((real(Pp(k)-b_min)/d))+5)~=0
          if (Mm(floor(imag(Pp(k))/d),floor((real(Pp(k)-b_min)/d))+3)~=0)&(Mm(floor(imag(Pp(k))/d),floor((real(Pp(k)-b_min)/d))+5)~=0)
          P(kk)=Pp(k);
            kk=kk+1;
         end
         
    else
            if (real(Pp(k))<=b_max)&(real(Pp(k))>=b_map)&imag(Pp(k))<=w_map
               P(kk)=Pp(k);
               kk=kk+1;
            end
    end
        
end   

'Increasing accuracy by the Newton method'
clear Pp;
Pp=P;
clear P;
%num. zpresneni polu

if e==-1
   e=0.001*d
end

   H=M/diff(M);
    n=size(Pp);
    n=n(2);
for o=1:n
    a=Pp(o);
    Newton_run=1;
while Newton_run
    a_k1=a;
    s=a;
    a=a-eval(H);
    Newton_run=(abs(a-a_k1)>e);
end
  if (~isreal(a))&(abs(imag(a))<0.1*e)  %orezani imag polu 
      a=real(a);
end
      
  P(o,1)=a;

end

Pmax=max(abs(P));
for k=1:n
    [Pm o]=min(abs(P));
    PM(k,1)=P(o);
    P(o)=2*Pmax;
end

P=PM;
    

figure(4)
plot(real(P),imag(P),'ko','MarkerFaceColor',[0 0 0],'MarkerSize',3)
figure(3)
plot(real(P),imag(P),'ko','MarkerFaceColor',[0 0 0],'MarkerSize',3)
else
  P=[];   
end

end

function Vysledek=argOblast(M,hranice)
syms s

B=hranice(1,:);
W=hranice(2,:);
N=length(hranice);

Q=subs(M,s,B'+j*W');
Arg_pi=angle(Q);
%figure(10)
% plot(1:N,Arg_pi)
Arg=0;

for k=1:N-1
    if Arg_pi(k)-Arg_pi(k+1)>pi/2;
       Arg=Arg+1;
   elseif Arg_pi(k)-Arg_pi(k+1)<-pi/2;
       Arg=Arg-1;
   end
end

Vysledek=Arg;
