function [Ewave] = storecx(varargin);
% [Ewave] = storeLast20(Ewave,timesteps,ind,passdata);
% modified 10/15/2012 by Joseph Hollmann
% Function now tells the program what options it requires
% modified by Chuck DiMarzio Sat Oct 29 14:48:52 2011 from storeLastX
% written by Joseph Hollmann, 10/22/2011
% Module written to work with execute_FDTD
% function takes Ewave structure in and creates a new substructure in
% modules.storelast. This stores Ez_field from the last X number of runs, 
% where X is defined in passdata
% Ex: options.inputs.pass = {20}; options.outputs = ['Ewave'];
% options.inputs.names = {'Ewave','timesteps','ind'}; options.callfcn = 'storeLastX';


persistent ti Eztemp Hxtemp Hytemp ran;

% input check
if isempty(varargin),
    variables{1} = 'Ewave,media,timesteps,ind,options';
    variables{2} = 'Ewave';
    ind = [];
    Ewave = variables; % just used this once to return variable names
else,
    Ewave = varargin{1};
    media = varargin{2};
    timesteps = varargin{3};
    ind = varargin{4};
    options = varargin{5};
    
    % test to see if within loop
    if isempty(ind) == 1,
        ti = 0; % reset persistent variable
        if ~isfield(options.callfcn,'cxrange'),
            ran = 40;
        else,
            ran = options.callfcn.cxrange;
        end;
        
    elseif ind == timesteps-ran,
        % preallocate memory - but not until needed
        [y,x] = size(Ewave.Ez_field);
        Ewave.module.firsttimestep=ind;
        Ewave.module.lasttimestep=timesteps;
        Ewave.module.Ezcomplex = zeros(y,x);
        Eztemp=zeros(y,x,ran);
        Hxtemp=zeros(y,x,ran);
        Hytemp=zeros(y,x,ran);
        Ewave.module.Hxcomplex = zeros(y,x);
        Ewave.module.Hycomplex = zeros(y,x);
        ti = 1;
    elseif ind >= timesteps-ran+1,
        % store last X number of runs
        % if not, preallocate memory (important all modules do this)
        Eztemp(:,:,ti)=Ewave.Ez_field;
        Hxtemp(:,:,ti)=Ewave.Hx_field;
        Hytemp(:,:,ti)=Ewave.Hy_field;
        ti = ti+1;
    end;
    if ind == timesteps;
        Ewave.module.Ezcomplex(:,:) = cxcalc(Eztemp,Ewave,media);
        Ewave.module.Hxcomplex(:,:) = cxcalc(Hxtemp,Ewave,media);
        Ewave.module.Hycomplex(:,:) = cxcalc(Hytemp,Ewave,media);
        [Ewave.module.Ezcomplex] = trimBCpixels(Ewave.module.Ezcomplex,media.BC.num,2);
        [Ewave.module.Hycomplex] = trimBCpixels(Ewave.module.Hycomplex,media.BC.num,2);
        [Ewave.module.Hxcomplex] = trimBCpixels(Ewave.module.Hxcomplex,media.BC.num,2);
    end;
end;

function cxfield=cxcalc(field,Ewave,media);
%
%  Computes the conplex field from the last X field values in FDTD.
%  see saveLastx.m and execute_FDTD.m
%  by Chuck DiMarzio, NU/Caltech, October 2011
%
[tt,ff]=meshgrid([-1:1]*3e8/Ewave.lambda0,...
    [1:size(field,3)]*media.delTemporal);
m=exp(1i*2*pi*tt.*ff);
mi=pinv(m);
cxfield=zeros(size(field,1),size(field,2));
for nnn=1:size(field,3);
    cxfield=cxfield+mi(3,nnn)*field(:,:,nnn);
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;