function [Ewave,media] = execute_FDTD(filename,media,Ewave,timesteps,options);
% [Ewave,media] = execute_FDTD(filename,media,Ewave,timesteps,showresults,options);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% updated 9/27/2013 by Joseph Hollmann
% filename - name of file that will be created to store results
% media - object contains grid parameters
% Ewave - contains electromagnetic field parameters
% timesteps - Number of steps through FDTD code
% options - structure containing the following:
% options.truncate timesteps for future playback,
% options.callfnc - contains string of file to call. Allows for extra
% processing to occur without editing code. Requires the following:
% options.inputs.names - cell array that contains variables to be passed to
% callfnc
% options.pass - cell of other variables to be passed to the module.
% options.outputs - cell array of variable to be modified
% module will execute once before the FDTD loop and at every timestep after

%% check for existence of parameters
if ~exist('options','var'),options = []; end;
if ~isfield(options,'BC'),options.BC = [];end;

% set defaults for undefined parameters
[media.BC.type] = checkprop(options.BC,'type','cpml');
[showresults] = checkprop(options,'showoptions',1);
[numdisp] = checkprop(options,'numdisp',40);


%% setup boundary pixels and default parameters
if strcmpi(media.BC.type,'cpml'),
    [media.BC.num] = checkprop(options.BC,'num',10); % number of boundary pixels (for cpml only)
elseif strcmpi(media.BC.type,'mur'),
    media.BC.num = 5; % Mur's boundary condition Coefficient given in equation 15 EMC-23, No 4, Novemeber 1981
end;

%% insert extra space for BC's
% [y,x] = size(media.ior_matrix);
bc1 = repmat(media.ior_matrix(:,1),1,media.BC.num);
bc2 = repmat(media.ior_matrix(:,end),1,media.BC.num);
media.ior_matrix = [bc1,media.ior_matrix,bc2];

bc1 = repmat(media.ior_matrix(1,:),media.BC.num,1);
bc2 = repmat(media.ior_matrix(end,:),media.BC.num,1);
media.ior_matrix = [bc1;media.ior_matrix;bc2];

%% setup detector
if isfield(Ewave,'detector'),
    Ewave.detector.depind = round(Ewave.detector.dpos/media.delSpatial)+media.BC.num; % depth index location of detector
    strtpnt = floor(Ewave.detector.delay/media.delTemporal); % start point in time steps
    tm = timesteps-strtpnt;% utilize minimum value to preallocate timesteps
    Ewave.detector.timesig = zeros(timesteps-strtpnt,size(media.ior_matrix,2));
    detectflag = 1;
else,
    detectflag = 0;
    strtpnt = [];
    tm = [];
end;

%% constants
Ewave.source.c0 = 299792458; % speed of light in a vacuum (m/s) (updated)
media.mu0 = 4e-7*pi;
media.eps0 = 1/(media.mu0*Ewave.source.c0^2); % do not approximate

%% create source
[Ewave] =  sourceval(Ewave,media,0);

%% create Matricies for FDTD
[y,x] = size(media.ior_matrix); % update size
Ewave.Ez_field = zeros(y,x); Ewave.Hx_field = Ewave.Ez_field; Ewave.Hy_field = Ewave.Ez_field;

%% define constants for FDTD
media.hcnst = media.delTemporal/(media.mu0*media.delSpatial); % This is a constant for the FDTD (computed outside the loop for speed)
media.ecnst = media.delTemporal./(media.eps0*media.ior_matrix.^2.*media.delSpatial);

%% visualzation setup
findfig = findobj('name',2); if sum(findfig == 2) == 1, close(2); end;

if showresults == 1,
    w = media.delSpatial:Ewave.lambda0:media.width*1e6; % convert to \mum
    d = media.delSpatial:Ewave.lambda0:media.depth*1e6;
    figure(2); h = imagesc(w,d,Ewave.Ez_field); colorbar;
    ht = get(gca,'title');
end;

% allows for a seperate function to be called at startup and after every
% timestep - called right before the loop
if isfield(options,'callfcn'),
    if isfield(options.callfcn,'file'),
        ind = []; % declared just in case timesteps are needed by module
        eval(['[variables] = ',options.callfcn.file,'(',ind,');']);
        inputstr = variables{1};
        outputstr = variables{2};
        eval([outputstr,' = ',options.callfcn.file,'(',inputstr,');']);
    else,
        error('File name does not exist. This program will crash now');
    end;
end;


%% run FDTD code
for ind = 1:timesteps,
    % update source
    n = ind*media.delTemporal;
    [Ewave,media] = FDTD_2D(Ewave,media,n,ind);
    %% collect detector data
    if (ind > strtpnt) & (ind <= strtpnt+tm) & detectflag == 1,
        Ewave.detector.timesig(ind-strtpnt,:) = Ewave.Ez_field(Ewave.detector.depind,:);
    end;
    
    % visualize result
    if showresults == 1&ind/numdisp == round(ind/numdisp),
        set(h,'cdata',Ewave.Ez_field); %title(num2str(ind)); drawnow;
        set(ht,'string',num2str(ind));
        drawnow;
    end;
    
    % call auxilary fcn. Error checks should have executed earlier
    if isfield(options,'callfcn'),
        eval([outputstr,' = ',options.callfcn.file,'(',inputstr,');']);
    end;
end;

% trim BC pixels
[Ewave.Ez_field] = trimBCpixels(Ewave.Ez_field,media.BC.num,2);
[Ewave.Hx_field] = trimBCpixels(Ewave.Hx_field,media.BC.num,2);
[Ewave.Hy_field] = trimBCpixels(Ewave.Hy_field,media.BC.num,2);
if detectflag == 1,
    try 
        [Ewave.detector.timesig] = trimBCpixels(Ewave.detector.timesig,media.BC.num,1);
    catch
        warning('Problem triming detector signal');
    end;
end;
[media.ior_matrix] = trimBCpixels(media.ior_matrix,media.BC.num,2);

% save results
save(filename,'Ewave','media','-v7.3');


function [val] = checkprop(st,field,orig);
% [val] = checkprop(st,field,orig);
% Function was created to check for user dictated properties

if isfield(st,field),
    val = getfield(st,field);
else,
    val = orig;
end;


function [val] = trimBCpixels(val,numpixels,dim);
% [val] = trimBCpixels(val,numpixels)
% Function trims BC pixels from matricies

val = val(:,numpixels+1:end-numpixels);
if dim == 2,
    val = val(numpixels+1:end-numpixels,:);
end;




