%implicit algorithm

clear all;

h2=2;                        % length of dipole (m)
a=0.002;                     % radius of dipole (m)
Nz=40;                       % number of segments (even number)
LR=4e-7;                     % length of response (s)
N=Nz-1;

c=3e+8;                      % speed of light
mi0=pi*4e-07;                % permeability of vacuum
ep0=8.85e-12;                % permittivity of vacuum

delz = h2/Nz;                % length of segment
delt = delz/c;               % length of time step

E0=120*pi;
T = 6e-9;                    % width of Gaussian pulse
t0 = 8e-9;                   % delay Gaussina pulse
f0 = 0;                      % center frequency of band  
t = 0:delt:LR;               % time
Nt=length(t);                % Number of time steps
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);           % kappa according to (5.1B.5)
for m=0:Nz                   % (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);        
 			       
  % distance of discretization elements (5.1B.7) a (5.1B.12), optional choice
    % if c*delt>>a is met, it is better to use(5.1B.7b) and (5.1B.12b) neglecting
    % of radius of dipole(layer B)
  %R(m+1)=sqrt(a^2+(m*delz)^2); %(5.1B.7a) and (5.1B.12a)
  R(m+1)= m*delz;              %(5.1B.7b) and (5.1B.12b) 
end                          
 
                          
kappaFI=-1/4/pi/ep0*kappa;    % kappa for scalar potential (5.1B.13)
kappa=mi0/4/pi*kappa;         % kappa for vector potential (5.1B.5)

A11C=zeros(N,N);             % coefficients of vector potential (5.1B.35) and (5.1B.22)
A12C=zeros(N,N);

for m=1:N,  
    for n=1:N,
        mn=abs(m-n);          %consider terms at time tk>=tRk>tk-1
        if(R(mn+1)<cdelt)     % distance has to be smaller than 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);           % coefficients of scalar potential (5.1B.36) 
FI12C=zeros(N,N);           % for ++ according to (5.1B.27b) and (5.1B.29) 
FI13C=zeros(N,N);           % for other combinations analogically 


for m=1:N,
    for n=1:N,
        mn=abs(m-n);        % computing of distance indexes between elements 
        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;
                             % terms contribute if the following 
                             % condition is met: 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;
                        % inverse matrix (5.1B.35)
invA11CFI11C=inv(A11C+delt/delz/2*FI11C); 

I=zeros(N,Nt);          %current
IS=zeros(N,Nt);         %integral of current

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

E=zeros(N,1);           % excitation vector
nap=(N+1)/2;            % dipole is excited at its center

                        % TIME LOOP
for k=3:Nt,
   A2=zeros(N,1);       % (5.1B.19b)
   FI2=zeros(N,1);      % (5.1B.36d), (5.1B.25b) and further analogically for +-,-+,-- 
   for m=1:N,
      for n=1:N,
                        % computing of distance indexes between elements
         mn=abs(m-n);   % vector potential
         mnPP=mn;       % scalar potential
         mnPM=mn+1;     % P is +; M is -
         mnMP=abs(mn-1);
         mnMM=mn;
         if(m>=n) 
            hlp=mnPM;
            mnPM=mnMP;
            mnMP=hlp;
         end;
                        % further the terms are computed
                        % if: tRk<=tk-1
      	 if(R(mn+1)>=cdelt) % vector potential (5.1B.19b) 
            tk=k-R(mn+1)/cdelt;
            if(tk>1)
               ftk=floor(tk); %linear approximation of current
               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)   % scalar potential ++; (5.1B.36d), (5.1B.25b) 
            tk=k-R(mnPP+1)/cdelt;
            if(tk>1)
               ftk=floor(tk);   %linear approximation of integral of current
               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;
                      % scalar potential +-; (5.1B.36d),similarly like (5.1B.25b)
         if(R(mnPM+1)>=cdelt)   
            tk=k-R(mnPM+1)/cdelt;
            if(tk>1)
               ftk=floor(tk);    %linear approximation of integral of current
               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;
                      % scalar potential +-; (5.1B.36d),similarly like (5.1B.25b)
         if(R(mnMP+1)>=cdelt)
            tk=k-R(mnMP+1)/cdelt;
            if(tk>1)
               ftk=floor(tk);    %linear approximation of integral of current
               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;
                      % scalar potential +-; (5.1B.36d),similarly like (5.1B.25b)
         if(R(mnMM+1)>=cdelt)
            tk=k-R(mnMM+1)/cdelt;
            if(tk>1)
               ftk=floor(tk);    %linear approximation of integral of current
               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);    % excitation vector is non-zero at the location of feeding (5.1B.7)
                    % (B.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));
                  
                    % integration of current with respect to time (trapezoid rule) (5.1B.15)
 IS(:,k)=(I(:,k)+I(:,k-1))*delt/2+IS(:,k-1);
                   
                    % preparation for the next step 
 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)
                     % considering (5.1B.26)-(5.1B.29)
end;

%computing the excitation pulse once more, since the original one is 
% about delt/2 delayed
E=E0*4/(c*T*sqrt(pi)).*exp(-(4/T)^2*(t-t0).*(t-t0)).*cos(2*pi*f0*t);

%display of excitation pulse
t=t*1e9; 
figure(1);hold on;plot(t,E,'b');
xlabel('Time [ns]');ylabel('Intensity of el. field [V/m]');
box on;

%display of current response 
                           % time is in ns
It=I(nap, :);                       % current response at feeding location
figure(2);hold on;plot(t,It,'r');
xlabel('Time [ns]');ylabel('Current [I]');
box on;


%computing frequency axis and input impedance
deltf=1/LR;                         %length of frequency step in spectrum
Fmax=400e+6;                        % max frequency

f=(deltf:deltf:Fmax)/1e6;           % computing of frequency in MHz
Nf=length(f);                       % number of frequencies 
Ut=delz*E;                          % voltage at antenna feeding 
Uf=fft(Ut);                         % Fourier's transformation of excitation
If=fft(It);                         % and current response
Zf=Uf(2:Nf+1)./If(2:Nf+1);          % computing of input impedance

% display of input impedance
figure(3);hold on;plot(f,real(Zf(1:Nf)),'r');
xlabel('Frequency [MHz]');ylabel('Real part of input impedance [\Omega]');
box on;

figure(4);hold on;plot(f,imag(Zf(1:Nf)),'r');
xlabel('Frequency [MHz]');ylabel('Imaginary part of input impedance [\Omega]');
box on;
