function [Ewave,media] = FDTD_2D(Ewave,media,n,ind);
% [Ewave] = FDTD_2D(Ewave,media,n,x,y,BC);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Written by Joseph Hollmann for the Optical Science Lab at Northeastern
% Copyright holder is the Optical Science Lab, May 2011
% This file may be copied, used, or modified for educational and research purposes provided
% that this header information is not removed or altered. Other distribution is prohibited
% without permission. We assume no liability for use of the code and no obligation to provide
% support, Nevertheless we would like to improve this package, and would be happy to receive comments
% Email - Hollmann@ece.neu.edu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Propogates a wave from a source through a pre-defined media.
% Takes input from execute_FDTD
% Requires structures Ewave and media. x and y are
% predefined vectors that should be declared in program calling FDTD
% function. n - current timestep.  BC - boundary condition.
% Returns E-field in Ewave structure.

% declare variables to be kept around
persistent x y X XX Y YY bye cye byh cyh cxe bxe cxh bxh epsx epsy num rbci %xbcf ybcf

if ind == 1,
    if strcmpi(media.BC.type,'CPML'),
        %% setup FDTD matricies
        [y,x] = size(media.ior_matrix);
        X = 1:x-1;
        XX = 2:x-1;
        Y = 1:y-1;
        YY = 2:y-1;
        num = media.BC.num;
        rbci = 1:num;
        
        media.BC.val = 1;
        [bye,byh,cye,cyh,bxe,bxh,cxe,cxh,epsx,epsy] = setupparameters(media,num); % call setup for PML
        Ewave.BC.psiezy = zeros(num*2,x);
        Ewave.BC.psihxy = zeros(num*2,x);
        Ewave.BC.psiezx = zeros(y,num*2);
        Ewave.BC.psihyx = zeros(y,num*2);
    else,
        media.BC.val = 0;
    end;
    
    Ewave.source.dposi = round(Ewave.source.dpos/media.delSpatial)+num; % fixes an issue with previous versions
end;

if media.BC.val == 1,
    %% test to see if pulse is still active
    if ind < Ewave.source.fin,
        [Ewave] =  sourceval(Ewave,media,n);
         Ewave.Ez_field(Ewave.source.dposi,Ewave.source.indicies) = Ewave.Ez_field(Ewave.source.dposi,Ewave.source.indicies)+Ewave.source.instant;
%          Ewave.Ez_field(Ewave.source.dposi,Ewave.source.indicies) = Ewave.source.instant;
    end;

    % CPML - top and bottom
    Ewave.BC.psiezy(1:num-1,XX) = bye(2:num,XX).*Ewave.BC.psiezy(1:num-1,XX)+cye(2:num,XX).*(Ewave.Hx_field(YY(1:num-1),XX)-Ewave.Hx_field(YY(1:num-1)-1,XX))/media.delSpatial;
    Ewave.BC.psiezy(num+1:end,XX) = bye(num+1:end,XX).*Ewave.BC.psiezy(num+1:end,XX)+cye(num+1:end,XX).*(Ewave.Hx_field(YY(end-num+1:end)+1,XX)-Ewave.Hx_field(YY(end-num+1:end),XX))/media.delSpatial;
    Ewave.BC.psihxy(rbci,X) = byh(1:num,X).*Ewave.BC.psihxy(rbci,X)+cyh(rbci,X).*(Ewave.Ez_field(Y(rbci)+1,X)-Ewave.Ez_field(Y(rbci),X))/media.delSpatial;
    Ewave.BC.psihxy(num+1:end,X) = byh(num+1:end-1,X).*Ewave.BC.psihxy(num+1:end,X)+cyh(num+1:end,X).*(Ewave.Ez_field(Y(end-num+1:end)+1,X)-Ewave.Ez_field(Y(end-num+1:end),X))/media.delSpatial;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CPML - Side to side
    Ewave.BC.psiezx(YY,1:num-1) = bxe(YY,2:num).*Ewave.BC.psiezx(YY,1:num-1)+cxe(YY,2:num).*(Ewave.Hy_field(YY,XX(1:num-1))-Ewave.Hy_field(YY,XX(1:num-1)-1))/media.delSpatial;
    Ewave.BC.psiezx(YY,num+1:end) = bxe(YY,num+1:end).*Ewave.BC.psiezx(YY,num+1:end)+cxe(YY,num+1:end).*(Ewave.Hy_field(YY,XX(end-num+1:end)+1)-Ewave.Hy_field(YY,XX(end-num+1:end)))/media.delSpatial;
    Ewave.BC.psihyx(Y,rbci) = bxh(Y,rbci).*Ewave.BC.psihyx(Y,rbci)+cxh(Y,rbci).*(Ewave.Ez_field(Y,X(rbci)+1)-Ewave.Ez_field(Y,X(rbci)))/media.delSpatial;
    Ewave.BC.psihyx(Y,num+1:end) = bxh(Y,num+1:end-1).*Ewave.BC.psihyx(Y,num+1:end)+cxh(Y,num+1:end).*(Ewave.Ez_field(Y,X(end-num+1:end)+1)-Ewave.Ez_field(Y,X(end-num+1:end)))/media.delSpatial;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    if sum(sum(Ewave.BC.psiezx)) ~= 0,
        dude = 1;
    end;
    
    %% actual FDTD computation
    Ewave.Hx_field(Y,X) = Ewave.Hx_field(Y,X)-media.hcnst.*(Ewave.Ez_field(Y+1,X)-Ewave.Ez_field(Y,X));
    Ewave.Hy_field(Y,X) = Ewave.Hy_field(Y,X)+media.hcnst.*(Ewave.Ez_field(Y,X+1)-Ewave.Ez_field(Y,X));
    Ewave.Ez_field(YY,XX) = Ewave.Ez_field(YY,XX)+media.ecnst(YY,XX).*(Ewave.Hy_field(YY,XX)-Ewave.Hy_field(YY,XX-1)+Ewave.Hx_field(YY-1,XX)-Ewave.Hx_field(YY,XX));
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    Ewave.Ez_field(YY(rbci),XX) = Ewave.Ez_field(YY(rbci),XX)-media.delTemporal/media.eps0./epsy(rbci,XX).*Ewave.BC.psiezy(rbci,XX);
    Ewave.Ez_field(YY(end-num+1:end)+1,XX) = Ewave.Ez_field(YY(end-num+1:end)+1,XX)-media.delTemporal/media.eps0./epsy((num+1:end),XX).*Ewave.BC.psiezy(num+1:end,XX);
    Ewave.Hx_field(Y(rbci),X) = Ewave.Hx_field(Y(rbci),X)-media.delTemporal/media.mu0*Ewave.BC.psihxy(rbci,X);
    Ewave.Hx_field(Y(end-num+1:end),X) = Ewave.Hx_field(Y(end-num+1:end),X)-media.delTemporal/media.mu0*Ewave.BC.psihxy(num+1:end,X);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    Ewave.Ez_field(YY,XX(rbci)) = Ewave.Ez_field(YY,XX(rbci))+media.delTemporal/media.eps0./epsx(YY,rbci).*Ewave.BC.psiezx(YY,rbci);
    Ewave.Ez_field(YY,XX(end-num+1:end)+1) = Ewave.Ez_field(YY,XX(end-num+1:end)+1)+media.delTemporal/media.eps0./epsx(YY,num+1:end).*Ewave.BC.psiezx(YY,num+1:end);
    Ewave.Hy_field(Y,X(rbci)) = Ewave.Hy_field(Y,X(rbci))+media.delTemporal/media.mu0*Ewave.BC.psihyx(Y,rbci);
    Ewave.Hy_field(Y,X(end-num+1:end)) = Ewave.Hy_field(Y,X(end-num+1:end))+media.delTemporal/media.mu0*Ewave.BC.psihyx(Y,num+1:end);
else,
    [Ewave] = FDTDmur(media,Ewave,ind,n);
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [bye,byh,cye,cyh,bxe,bxh,cxe,cxh,epsx,epsy] = setupparameters(media,num);
% sets up parameters for CPML boundary conditions
%% set parameters
eta0 = sqrt(media.mu0/media.eps0);
dih = [0:num-1];
die = [0.5:num];
m = 3;
alpham = 1;
alphamax = 0.05;
kappa_e = 1+(1.1-1)*(die/num).^m;
alpha_e = alphamax*((num-die)/num).^(alpham);
sigma_e = 0.8*(m+1)/eta0/media.delSpatial*((die)/num).^m;
kappa_h = 1+(1.1-1)*(dih/num).^m;
alpha_h = alphamax*((num-dih)/num).^(alpham);
sigma_h = 0.8*(m+1)/eta0/media.delSpatial*((dih)/num).^m;
b_e = exp(-(sigma_e./kappa_e+alpha_e)*media.delTemporal/media.eps0);
b_h = exp(-(sigma_h./kappa_h+alpha_h)*media.delTemporal/media.eps0);
ce = sigma_e./(sigma_e.*kappa_e+kappa_e.^2.*alpha_e).*(b_e-1);
ch = sigma_h./(sigma_h.*kappa_h+kappa_h.^2.*alpha_h).*(b_h-1);
ce(ce~=ce) = 0;
ch(ch~=ch) = 0;

% create matricies for 
bye1D = repmat(b_e',1,size(media.ior_matrix,2));
byh1D = repmat(b_h',1,size(media.ior_matrix,2));
cye1D = repmat(ce',1,size(media.ior_matrix,2));
cyh1D = repmat(ch',1,size(media.ior_matrix,2));

bye = zeros(2*size(bye1D,1),size(bye1D,2));
cye = zeros(2*size(cye1D,1),size(cye1D,2));
bye(1:num,:) = bye1D(end:-1:1,:);
bye(num+1:end,:) = bye1D(1:end,:);
cye(1:num,:) = cye1D(end:-1:1,:);
cye(num+1:end-1,:) = cye1D(1:end-1,:);
byh = zeros(2*size(bye1D,1)+1,size(bye1D,2));
byh(1:num,:) = byh1D(end:-1:1,:);
byh(num+1:end-1,:) = byh1D;
cyh = zeros(2*size(cye1D,1),size(cye1D,2));
cyh(1:num,:) = cyh1D(end:-1:1,:);
cyh(num+1:end,:) = cyh1D;

%%%%%
bxe1D = repmat(b_e,size(media.ior_matrix,1),1);
bxh1D = repmat(b_h,size(media.ior_matrix,1),1);
cxe1D = repmat(ce,size(media.ior_matrix,1),1);
cxh1D = repmat(ch,size(media.ior_matrix,1),1);

%%%%%%%%%%%%%%%%%%%%%%%%
bxe = zeros(size(bxe1D,1),2*size(bxe1D,2));
cxe = zeros(size(cxe1D,1),2*size(cxe1D,2));
bxh = zeros(size(bxe1D,1),2*size(bxe1D,2)+1);
cxh = zeros(size(cxe1D,1),2*size(cxe1D,2));
bxe(:,1:num) = bxe1D(:,end:-1:1);
bxe(:,num+1:end) = bxe1D(:,1:end);
cxe(:,1:num) = cxe1D(:,end:-1:1);
cxe(:,num+1:end-1) = cxe1D(:,1:end-1);
bxh(:,1:num) = bxh1D(:,end:-1:1);
bxh(:,num+1:end-1) = bxh1D;
cxh(:,1:num) = cxh1D(:,end:-1:1);
cxh(:,num+1:end) = cxh1D;

% set up indicies of refraction
epsx = ([repmat(media.ior_matrix(:,1),1,num),repmat(media.ior_matrix(:,end),1,num)]).^2;
epsy = ([repmat(media.ior_matrix(1,:),num,1);repmat(media.ior_matrix(end,:),num,1)]).^2;

% solves FDTD equations using Mur BC's. In here for historical reasons
function [Ewave] = FDTDmur(media,Ewave,ind,n);
persistent x y X XX Y YY;

if abs(n-media.delTemporal) < ind,
    %% create Matricies for FDTD
    [y,x] = size(media.ior_matrix);
    
    % set up index vectors
    X = 2:x-1;
    XX = 3:x-2;
    Y = 2:y-1;
    YY = 3:y-2;
    
    media.murr = (Ewave.source.c0*media.delTemporal-media.delSpatial)/(Ewave.source.c0*media.delTemporal+media.delSpatial);
    Ewave.source.dposi = round(Ewave.source.dpos/media.delSpatial)+5; % fixes an issue with previous versions
end;

% test to see if pulse is still active
if ind < Ewave.source.fin,
    [Ewave] =  sourceval(Ewave,media,n);
    Ewave.Ez_field(Ewave.source.dposi,Ewave.source.indicies) = Ewave.Ez_field(Ewave.source.dposi,Ewave.source.indicies)+Ewave.source.instant;  
end;

% store Ez fields for BC's
Eholdx = [Ewave.Ez_field(3,X);Ewave.Ez_field(end-2,X)];
Eholdy = [Ewave.Ez_field(Y,3),Ewave.Ez_field(Y,end-2)]'; % store field

% calculate magnetic field
Ewave.Hx_field(Y,X) = Ewave.Hx_field(Y,X)-media.hcnst.*(Ewave.Ez_field(Y+1,X)-Ewave.Ez_field(Y,X));
Ewave.Hy_field(Y,X) = Ewave.Hy_field(Y,X)+media.hcnst.*(Ewave.Ez_field(Y,X+1)-Ewave.Ez_field(Y,X));

Ewave.Ez_field(YY,XX) = Ewave.Ez_field(YY,XX)+media.ecnst(YY,XX).*(Ewave.Hy_field(YY,XX)...
    -Ewave.Hy_field(YY,XX-1)+Ewave.Hx_field(YY-1,XX)-Ewave.Hx_field(YY,XX));


Ewave.Ez_field(2,X) = Eholdx(1,:)+media.murr.*(Ewave.Ez_field(3,X)-Ewave.Ez_field(2,X));
Ewave.Ez_field(end-1,X) = Eholdx(2,:)+media.murr.*(Ewave.Ez_field(end-2,X)-Ewave.Ez_field(end-1,X));

Ewave.Ez_field(Y,2) = Eholdy(1,:).'+media.murr.*(Ewave.Ez_field(Y,3)-Ewave.Ez_field(Y,2));
Ewave.Ez_field(Y,end-1) = Eholdy(2,:).'+media.murr.*(Ewave.Ez_field(Y,end-2)-Ewave.Ez_field(Y,end-1));

