%implicitni algoritmus

clear all;

h2=2;                        % delka dipolu (m)
a=0.002;                     % polomer dipolu (m)
Nz=40;                       % pocet segmentu (musi byt sude cislo)
LR=4e-7;                     % delka odezvy proudu (s)
N=Nz-1;

c=3e+8;                      % rychlost svetla
mi0=pi*4e-07;                % permeabilita vakua
ep0=8.85e-12;                % permitivita vakua

delz = h2/Nz;                % delka prostoroveho segmentu
delt = delz/c;               % delka casového segmentu = cas, za ktery

E0=120*pi;
T = 6e-9;                    % sirka Gaussova impulsu
t0 = 8e-9;                   % zpozdeni Gaussova impulsu
f0 = 0;                      % kmitocet stredu pasma analyzy  
t = 0:delt:LR;               % vzorkovaci okamziky
Nt=length(t);                % Nt = pocet casových kroku
cdelt=c*delt;

EG=E0*4/(c*T*sqrt(pi)).*exp(-(4/T)^2*(t-t0-delt/2).*(t-t0-delt/2)).*cos(2*pi*f0*(t-delt/2));;

kappa=zeros(Nz+1,1);
Rmn=zeros(Nz+1,1);           % priprava vypoctu kappa dle vztahu (5.1B.5)
for m=0:Nz                   % vztah (5.1B.13)
  hlp1=log((m+0.5)*delz + sqrt(((m+0.5)*delz)^2+a^2));
  hlp2=log((m-0.5)*delz + sqrt(((m-0.5)*delz)^2+a^2));
  kappa(m+1)=(hlp1-hlp2);        
 
    % vzdalenost diskret. elem. (5.1B.7) a (5.1B.12), mame moznost volby
    % pokud je splneno c*delt>>a je lepe pouzit (5.1B.7b) a (5.1B.12b) zanedbat
    % polomer vodice (viz Vrstva B)
  %R(m+1)=sqrt(a^2+(m*delz)^2);  %(5.1B.7a) a (5.1B.12a)
  R(m+1)= m*delz;               %(5.1B.7b) a (5.1B.12b) 

end                          
 
                          
kappaFI=-1/4/pi/ep0*kappa;    % kappa pro skalarni pot. (5.1B.13)
kappa=mi0/4/pi*kappa;         % kappa pro vektorovy pot. (5.1B.5)

A11C=zeros(N,N);              % koef. vekt. pot. dle (5.1B.35) a (5.1B.22)
A12C=zeros(N,N);

for m=1:N,  
    for n=1:N,
        mn=abs(m-n);          %uvazuj jen cleny v case tk>=tRk>tk-1
        if(R(mn+1)<cdelt)     %tj. vzdalenst musi byt mensi nez c*delt 
            A11C(m,n)=A11C(m,n)+kappa(mn+1)*(1-R(mn+1)/cdelt);
            A12C(m,n)=A12C(m,n)+kappa(mn+1)*R(mn+1)/cdelt;
        end;
    end;
end;

FI11C=zeros(N,N);           % koef. skal. pot. dle (5.1B.36) 
FI12C=zeros(N,N);           % pro ++ dle (5.1B.27b) a (5.1B.29); 
FI13C=zeros(N,N);           % pro ostatni kombinace analogicky 


for m=1:N,
    for n=1:N,
        mn=abs(m-n);        % vypocet indexu vzdalenosti mezi elementy 
        mnPP=mn;            % P je +; M je - 
        mnPM=mn+1;
        mnMP=abs(mn-1);
        mnMM=mn;
        if(m>=n) 
            hlp=mnPM;
            mnPM=mnMP;
            mnMP=hlp;
        end;
                             % prispivaji jen cleny,pro ktere plati
                             % tRk>tk-1
        if(R(mnPP+1)<cdelt)
            FI11C(m,n)=FI11C(m,n)...             %(5.1B.29a)
                       -delt/delz/2*(1-R(mnPP+1)/cdelt)^2*kappaFI(mnPP+1);
            FI12C(m,n)=FI12C(m,n)...             %(5.1B.29b) 
                       -delt/delz/2*(1-(R(mnPP+1)/cdelt)^2)*kappaFI(mnPP+1);          
            FI13C(m,n)=FI13C(m,n)-1/delz*kappaFI(mnPP+1); %(5.1B.27b)
        end;
        if(R(mnPM+1)<cdelt)
            FI11C(m,n)=FI11C(m,n)+delt/delz/2*(1-R(mnPM+1)/cdelt)^2*kappaFI(mnPM+1);
            FI12C(m,n)=FI12C(m,n)+delt/delz/2*(1-(R(mnPM+1)/cdelt)^2)*kappaFI(mnPM+1);          
            FI13C(m,n)=FI13C(m,n)+1/delz*kappaFI(mnPM+1);
        end;
        if(R(mnMP+1)<cdelt)
            FI11C(m,n)=FI11C(m,n)+delt/delz/2*(1-R(mnMP+1)/cdelt)^2*kappaFI(mnMP+1);
            FI12C(m,n)=FI12C(m,n)+delt/delz/2*(1-(R(mnMP+1)/cdelt)^2)*kappaFI(mnMP+1);          
            FI13C(m,n)=FI13C(m,n)+1/delz*kappaFI(mnMP+1);
        end;
        if(R(mnMM+1)<cdelt)
            FI11C(m,n)=FI11C(m,n)-delt/delz/2*(1-R(mnMM+1)/cdelt)^2*kappaFI(mnMM+1);
            FI12C(m,n)=FI12C(m,n)-delt/delz/2*(1-(R(mnMM+1)/cdelt)^2)*kappaFI(mnMM+1);          
            FI13C(m,n)=FI13C(m,n)-1/delz*kappaFI(mnMM+1);
        end;
    end;
end;
                        % vypocet inverzni matice dle (5.1B.35)leva strana
invA11CFI11C=inv(A11C+delt/delz/2*FI11C); 

I=zeros(N,Nt);          %hledany proud
IS=zeros(N,Nt);         %integral proudu

A=zeros(N,1);          
FI=zeros(N,1);

E=zeros(N,1);           % vektor buzení
nap=(N+1)/2;            %dratovy dipol je napajen ve svem stredu

                        % CASOVA SMYCKA
for k=3:Nt,
   A2=zeros(N,1);       % (5.1B.19b)
   FI2=zeros(N,1);      % (5.1B.36d), (5.1B.25b) a dale analog. pro +-,-+,-- 
   for m=1:N,
      for n=1:N,
                        % vypocet indexu vzdalenosti mezi elementy
         mn=abs(m-n);   % vektorovy potencial
         mnPP=mn;       % skalarni potencial
         mnPM=mn+1;     % P je +; M je -
         mnMP=abs(mn-1);
         mnMM=mn;
         if(m>=n) 
            hlp=mnPM;
            mnPM=mnMP;
            mnMP=hlp;
         end;
                        % dale jsou pocitany jen cleny,
                        % pro ktere plati: tRk<=tk-1
      	 if(R(mn+1)>=cdelt) % vektorovy potencial (5.1B.19b) 
            tk=k-R(mn+1)/cdelt;
            if(tk>1)
               ftk=floor(tk); %linearni aproximace proudu
               In=(I(n,ftk+1)-I(n,ftk))*(tk-ftk)+I(n,ftk);              
               A2(m,1)=A2(m,1)+In*kappa(mn+1);
            end;
         end;
      	
         if(R(mnPP+1)>=cdelt)   % skalarni potencial ++; (5.1B.36d), (5.1B.25b) 
            tk=k-R(mnPP+1)/cdelt;
            if(tk>1)
               ftk=floor(tk);   %linearni aproximace integralu proudu
               ISn=(IS(n,ftk+1)-IS(n,ftk))*(tk-ftk)+IS(n,ftk);              
               FI2(m,1)=FI2(m,1)-1/delz*ISn*kappaFI(mnPP+1);
            end;
         end;
                      % skalarni potencial +-; (5.1B.36d),podobne jako (5.1B.25b)
         if(R(mnPM+1)>=cdelt)   
            tk=k-R(mnPM+1)/cdelt;
            if(tk>1)
               ftk=floor(tk);    %linearni aproximace integralu proudu
               ISn=(IS(n,ftk+1)-IS(n,ftk))*(tk-ftk)+IS(n,ftk);              
               FI2(m,1)=FI2(m,1)+1/delz*ISn*kappaFI(mnPM+1);
            end;
         end;
                      % skalarni potencial +-; (5.1B.36d),podobne jako (5.1B.25b)
         if(R(mnMP+1)>=cdelt)
            tk=k-R(mnMP+1)/cdelt;
            if(tk>1)
               ftk=floor(tk);    %linearni aproximace integralu proudu
               ISn=(IS(n,ftk+1)-IS(n,ftk))*(tk-ftk)+IS(n,ftk);              
               FI2(m,1)=FI2(m,1)+1/delz*ISn*kappaFI(mnMP+1);
            end;
         end;
                      % skalarni potencial +-; (B.36d),podobne jako (5.1B.25b)
         if(R(mnMM+1)>=cdelt)
            tk=k-R(mnMM+1)/cdelt;
            if(tk>1)
               ftk=floor(tk);    %linearni aproximace integralu proudu
               ISn=(IS(n,ftk+1)-IS(n,ftk))*(tk-ftk)+IS(n,ftk);              
               FI2(m,1)=FI2(m,1)-1/delz*ISn*kappaFI(mnMM+1);
            end;
         end;
     end;     
 end;
 
    
 E(nap,1)=EG(k);    % vektor buzení je nenulovy jen v miste napajeni vztah (5.1B.7)
                    % (5.1B.35)                 
 I(:,k)=invA11CFI11C*(delt*E+A-delt/2/delz*FI-(A12C*I(:,k-1)+A2)...
                -delt/2/delz*(FI12C*I(:,k-1)+FI13C*IS(:,k-1)+FI2));
                  
                    % integrace proudu dle lichobeznikoveho pravidla (5.1B.15)
 IS(:,k)=(I(:,k)+I(:,k-1))*delt/2+IS(:,k-1);
                   
                    % vypocet pro dalsi krok 
 A=A11C*I(:,k)+A12C*I(:,k-1)+A2;    % (5.1B.18) a (5.1B.21)   
 FI=FI11C*I(:,k)+FI12C*I(:,k-1)+FI13C*IS(:,k-1)+FI2; %(5.1B.30)-(5.1B.33)
                     % s uvazenim (5.1B.26)-(5.1B.29)
end;

%jeste jednou vypoctu prubeh budiciho pulsu, protoze puvodni je o delt/2
%zpozden
E=E0*4/(c*T*sqrt(pi)).*exp(-(4/T)^2*(t-t0).*(t-t0)).*cos(2*pi*f0*t);

%zobrazeni budiciho pulsu
t=t*1e9;
figure(1);hold on;plot(t,E,'b');
xlabel('Èas [ns]');ylabel('Intenzita el. pole [V/m]');
box on;

%zobrazeni proudove odezvy
                            %cas je v ns
It=I(nap, :);                       % odezva proudu v miste napajeni
figure(2);hold on;plot(t,It,'r');
xlabel('Èas [ns]');ylabel('Proud [I]');
box on;


%vypocet deleni frekvencni osy a vstupni impedance
deltf=1/LR;                         %velikost frekvencniho kroku ve spektru
Fmax=400e+6;                        %maximalni kmitoce spektra

f=(deltf:deltf:Fmax)/1e6;           % vypocet frekvence c MHz
Nf=length(f);                       % pocet kmitoctu 
Ut=delz*E;                          % napeti na svorkach anteny
Uf=fft(Ut);                         % Fourierova transformace budiciho napeti
If=fft(It);                         % a odezvy proudu
Zf=Uf(2:Nf+1)./If(2:Nf+1);          % vypocet vstupni impedance

% zobrazeni vstupni impedance
figure(3);hold on;plot(f,real(Zf(1:Nf)),'r');
xlabel('Frekvence [MHz]');ylabel('Reálná èást vstupní impedance [\Omega]');
box on;

figure(4);hold on;plot(f,imag(Zf(1:Nf)),'r');
xlabel('Frekvence [MHz]');ylabel('Imaginární èást vstupní impedance [\Omega]');
box on;
