%The electric field distribution around an antenna of infinite length
clear 

i=sqrt(-1);
%speed of light
c=3e8;
%c=1;
%pause
%amplitude of Ey
D1=1;
lamda = 1;
%prop. distance
Lz=10*lamda;
Lx=Lz;

%freespace wavenumber
k0 = 2*pi/lamda;

%ang. freq
om = k0*c;


%max. time
maxT = Lz/c
f0 = om/(2*pi)
%pause
delz = lamda/20;
delx = delz;
%number of z samples
Nz = floor(Lz/delz)
Nx = floor(Lx/delx)
%pause
delz = Lz/Nz;
delx = Lx/Nx;

T0 = maxT/8;
spread = T0/2;

delt = delz/(2*c);
facz= (c*delt/delz)^2;
facx= (c*delt/delx)^2;
Nt=floor(maxT/delt)
%Nt=100
%pause
zz=0:delz:Nz*delz;

zzp = -delz:delz:Nz*delz+delz;
xxp = -delx:delx:Nx*delx+delx;
zzp = zzp - Nz*delz/2;
xxp = xxp - Nz*delz/2;

size(zzp)
Nz+3

zs = Nz/2 + 1
xs = Nx/2 + 1
pause(2)

Ey = zeros(3,Nx+3,Nz+3);



figure(1)
m=1;
Nt = 200;
for n = 1:Nt
    %Y0 = zeros(1,Nx+1);
    %exac = zeros(1,Nx+1);
    Eyy = zeros(Nx+3,Nz+3);
    t = (n-1)*delt;
    n0=mod(n-1,3)+1;
    np=mod(n,3)+1;
    nm=mod(n-2,3)+1;
    for jx=2:Nz+2
    for jz=2:Nz+2
     z = (jz-2)*delz;
     if((jz == zs)&&(jx == xs))  
           Ey(n0,jx,jz) = Ey(n0,jx,jz)+D1*cos(om*(t-delt));
           vecx = Ey(n0,jx-1:jx+1,jz);
           vecz = Ey(n0,jx,jz-1:jz+1);
           term0 = (facx*DDf(vecx) + facz*DDf(vecz));
           Ey(np,jx,jz) =  D1*cos(om*t) + term0 + 2*Ey(n0,jx,jz) - Ey(nm,jx,jz);
            
%              Ey(np,jx,jz) = exp(-0.5*((t-T0)/spread)^2);
%              Ey(n0,jx,jz) = exp(-0.5*((t-delt-T0)/spread)^2);
     else  
         
         vecx = Ey(n0,jx-1:jx+1,jz);
         vecz = Ey(n0,jx,jz-1:jz+1);
         term0 = (facx*DDf(vecx) + facz*DDf(vecz));
         Ey(np,jx,jz) = term0 + 2*Ey(n0,jx,jz) - Ey(nm,jx,jz);

     end 
    end
    end
    Eyy(1:Nx+3,1:Nz+3) = Ey(np,:,:);
    Hp=surf(xxp,zzp,Eyy)
    view(30,37.5)
     axis([min(xxp) max(xxp) min(zzp) max(zzp) -1 1])
     xlabel('x (m)')
     ylabel('y (m)')
     zlabel('Ez (V/m)')
%     axis([-delx (Nx+2)*delx -delz (Nz+2)*delz -1 1])
    pause(0.2)
    %pause(0.5)
%     Eyy(1,:)=Ey(np,xs,:);
%     plot(zzp,Eyy)
%     axis([-delz (Nz+1)*delz -2 2])
%     pause(0.5)
    n
end
    
