close all
clear all
clc
% Cubic grid Parameters
Nx=50;
Ny=50;
Nz=50;
dx=1;    % Lattice space incremetn
dy=1;    % Lattice space incremetn
dz=1;    % Lattice space incremetn
% Source location
xs=25,ys=25,zs=1
%%
% determine true permittivity and permeability matrices from supplied relative ones
ep0 = 8.8541878176e-12;             % dielectric permittivity of free space (Frads/m)
mu0 =4*pi*10^(-7);                 % magnetic permeability of free space 1.2566370614e-6 (henry/meter)
c0=1/sqrt(ep0*mu0);                 % Speed of light (m/s)
ep=1;
mu=1;
ep = ep*ep0;                        % true permittivity matrix
mu = mu*mu0;                        % true permeability matrix
%%
% determine time stepping parameters from supplied time vector
dt = 0.98/(c0*sqrt( (1/dx)^2 + (1/dy)^2 ));                                 % temporal discretization
numit = 500;                              % number of time steps% number of iterations
t=(0:1:numit-1).*dt;

tw=dt;
t0 = 4.0*tw;
T = (0:1:numit-1)'.*dt;
sourceb = -2.0*((T-t0)./tw).*exp(-1.0*((T-t0)./tw).^2.0); % Impolse electric source

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

DIELC = ones(Nx, Ny, Nz).*ep; %3D Permittivity array on half grid (cube centers)
MAGNC = ones(Nx+1, Ny+1, Nz+1).*mu; %3D Permeability array on integer grid (cube nodes)
RHOMC = zeros(Nx+1,Ny+1, Nz+1).*5; %3D Magnetic res. array on integer grid (cube nodes)Ohms/m) magnetic loose
SIGEC = zeros(Nx, Ny, Nz); %3D Electric conductivity array on half grid (cube centers)(siemens/m)
SIGEC(:,:,10:20)=0.5;


%% correcting for the center of the cell
nx = 2:Nx; ny = 2:Ny; nz =2:Nz;
Dtemp = (DIELC(nx,ny,:)+DIELC(nx-1,ny,:)+DIELC(nx,ny-1,:)+DIELC(nx-1,ny-1,:))/4;
Stemp = (SIGEC(nx,ny,:)+SIGEC(nx-1,ny,:)+SIGEC(nx,ny-1,:)+SIGEC(nx-1,ny-1,:))/4;
%%
Stemp = Stemp + 1e-9;
Ez= exp(-dt*Stemp./Dtemp);
Ez1 = Ez;
Ex1=Ez(:,:,nz); Ex1(nx,:,:)=Ex1;
Ey1=Ez(:,:,nz);Ey1(:,ny,:)=Ey1;

Ez =(1 - exp(-dt*Stemp./Dtemp))./(dx*Stemp);
Ez2 =Ez;
Ex=(1 - exp(-dt*Stemp./Dtemp))./(dy*Stemp);
Ex2 = Ex(:, :, nz);Ex2(nx,:,:)=Ex2;
Ey2=Ex(:,:,nz);Ey2(:,ny,:)=Ey2;

%
nx=1:Nx; ny = 1:Ny; nz = 1:Nz;
Mtemp = (MAGNC(:,ny,nz)+MAGNC(:,ny+1,nz)+MAGNC(:,ny,nz+1)+MAGNC(:,ny+1,nz+1))/4;
Rtemp = (RHOMC(:,ny,nz)+RHOMC(:,ny+1,nz)+RHOMC(:,ny,nz+1)+RHOMC(:,ny+1,nz+1))/4;
Rtemp = Rtemp + 1e-9;
%%

Hx = exp(-dt*Rtemp./Mtemp);
Hx1=Hx;
Hy1=Hx(nx,:,:); Hy1(:,ny+1,:)=Hy1;
Hz1=Hx(nx,:,:);Hz1(:,:,nz+1)=Hz1;

Hx2 = (1 - exp(-dt*Rtemp./Mtemp))./(dx*Rtemp + eps);
Hy= (1 - exp(-dt*Rtemp./Mtemp))./(dy*Rtemp + eps);
Hy2=Hy(nx,:,:);Hy2(:,ny+1,:)=Hy2;
Hz2=Hy(nx,:,:);Hz2(:,:,nz+1)=Hz2;

% Allocate field matrices
ElxO = zeros(Nx , Ny+1, Nz+1);
ElyO = zeros(Nx+1, Ny , Nz+1);
ElzO = zeros(Nx+1, Ny+1, Nz );
BxO = zeros(Nx+1, Ny , Nz );
ByO = zeros(Nx , Ny+1, Nz );
BzO = zeros(Nx , Ny , Nz+1);

ELx = zeros(Nx , Ny+1, Nz+1);
ELy = zeros(Nx+1, Ny , Nz+1);
ELz = zeros(Nx+1, Ny+1, Nz );
Bx = zeros(Nx+1, Ny , Nz );
By = zeros(Nx , Ny+1, Nz );
Bz = zeros(Nx , Ny , Nz+1);

f=0;
figure;
% Waitbar:
h=waitbar(1,['Time steping: Please wait...']);
for it=1:length(t),  
   
    %% E-field update (everywhere except on the boundary
    for i=1:Nx
        for j=1:Ny
            for k=1:Nz
                if ((j~=Ny)& (k~=Nz))
                    ElxO(i,j+1,k+1) = Ex1(i,j,k).*ELx(i,j+1,k+1)+Ex2(i,j,k).*((Bz(i,j+1,k+1)-Bz(i,j,k+1))-(By(i,j+1,k+1)-By(i,j+1,k)));
                end
                if ((i~=Nx)& (k~=Nz))
                    ElyO(i+1,j,k+1) = Ey1(i,j,k).*ELy(i+1,j,k+1)+Ey2(i,j,k).*((Bx(i+1,j,k+1)-Bx(i+1,j,k))-(Bz(i+1,j,k+1)-Bz(i,j,k+1)));
                end
                if ((i~=Nx)& (j~=Ny))
                    ElzO(i+1,j+1,k) = Ez1(i,j,k).*ELz(i+1,j+1,k)+Ez2(i,j,k).*((By(i+1,j+1,k)-By(i,j+1,k))-(Bx(i+1,j+1,k)-Bx(i+1,j,k)));
                end
             end
        end
    end    % Plane wave source
    [J]=source(Nx,Ny,Nz,dx,dy,dz,t(it));
    % ElzO=ElzO+J(:,:,1:Nz);
    ElzO(xs,ys,zs)=sourceb(it);
    
    
    %% H-field update (everywhere)
    for i=1:Nx+1
        for j=1:Ny+1
            for k=1:Nz+1
                %% H-field update (everywhere)
                if((j<=Ny)& (k<=Nz))
                    BxO(i,j,k) = Hx1(i,j,k).*Bx(i,j,k) + Hx2(i,j,k).*((ElyO(i,j,k+1)-ElyO(i,j,k))- (ElzO(i,j+1,k)-ElzO(i,j,k)));
                end
                if((i<=Nx)& (k<=Nz))
                    ByO(i,j,k) = Hy1(i,j,k).*By(i,j,k) + Hy2(i,j,k).*((ElzO(i+1,j,k)-ElzO(i,j,k))- (ElxO(i,j,k+1)-ElxO(i,j,k)));
                end
                if((i<=Nx)& (j<=Ny) & (k<=Nz))
                    BzO(i,j,k) = Hz1(i,j,k).*Bz(i,j,k) + Hz2(i,j,k).*((ElxO(i,j+1,k)-ElxO(i,j,k))- (ElyO(i+1,j,k)-ElyO(i,j,k)));
                end
            end
        end
    end
    ELx=ElxO;
    ELy=ElyO;
    ELz=ElzO;
    Bx=BxO;
    By=ByO;
    Bz=BzO;

    waitbar(it/length(t),h);
    [x,y,z] = meshgrid(1:Nx,1:Ny,1:Nz);
    %%
    subplot(2,2,1)
    % Impedance
    zzxy=ELy(1:Nx,1:Ny,1:Nz)./Bx(1:Nx,1:Ny,1:Nz);
    [CURLX, CURLY, CURLZ]=curl(x,y,z,Bx(1:Nx,1:Ny,1:Nz),By(1:Nx,1:Ny,1:Nz),Bz(1:Nx,1:Ny,1:Nz));
    [ECURLX, ECURLY, ECURLZ]=curl(x,y,z,ELx(1:Nx,1:Ny,1:Nz),ELy(1:Nx,1:Ny,1:Nz),ELz(1:Nx,1:Ny,1:Nz));
    div = divergence(x,y,z,ELx(1:Nx,1:Ny,1:Nz),ELy(1:Nx,1:Ny,1:Nz),ELz(1:Nx,1:Ny,1:Nz));
    
    %Visulization:
    
    xmin=min(x(:)); xmax=max(x(:));
    ymin=min(y(:)); ymax=max(y(:));
    zmin=min(z(:)); zmax=max(z(:));
    xlim([xmin xmax])
    ylim([ymin ymax])
    zlim([zmin zmax])
    
    hsurfaces=slice(x,y,z,SIGEC(1:Nx,1:Ny,1:Nz),25,25,25);
    set(hsurfaces,'EdgeColor','none',...
        'AmbientStrength',0.6)
    daspect([2,2,2]), view(30,40), axis tight;set(gca,'Zdir','reverse');
    xlabel('East(m)');ylabel('North(m)');zlabel('Elevation(m)');
    xlim([xmin xmax]);title('Conductivity');
    subplot(2,2,2)
    plot_fieldE = streamslice(x,y,z,ELx(1:Nx,1:Ny,1:Nz),ELy(1:Nx,1:Ny,1:Nz),ELz(1:Nx,1:Ny,1:Nz),25,25,25);
    set(plot_fieldE,'Color','black','LineWidth',1);
        
    %Magnetic field
    plot_fieldH = streamslice(x,y,z,BxO(1:Nx,1:Ny,1:Nz),ByO(1:Nx,1:Ny,1:Nz),BzO(1:Nx,1:Ny,1:Nz),25,25,25);
    set(plot_fieldH,'Color','Green','LineWidth',1);
    daspect([2,2,2]), view(30,40), axis tight;
    xlim([xmin xmax])
    ylim([ymin ymax])
    zlim([zmin zmax])
    xlabel('East(m)');ylabel('North (m)');zlabel('Elevation(m)');
    title('Electric and magnetic Field');
    set(gca,'Zdir','reverse')
    %streamline(x,y,z,ELx(1:Nx,1:Ny,1:Nz),ELy(1:Nx,1:Ny,1:Nz),ELz(1:Nx,1:Ny,1:Nz),[0:1:50],[0:1:50],[0:1:50])
    %%
    subplot(2,2,3);
    imagesc(squeeze(ELx(:,25,:))');title('Electric field');
    colorbar;      set(gca,'clim',[-1e-4 1e-4]);
    xlabel('East profile(m)');ylabel('Elevation(m)')
    %%
    subplot(2,2,4);
    imagesc(squeeze(Bx(:,25,:))');title('Magnetic field');
    colorbar;      set(gca,'clim',[-1e-10 1e-10]);
    xlabel('East profile(m)');ylabel('Elevation(m)');

    f=f+1;
    saveas(gcf, ['plane_',sprintf('%03d',f),'.png'])
    
end
close(h);


