function [Ewave] =  sourceval(Ewave,media,n);
% [Ewave] =  sourceval(Ewave,media,n);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% edited 9/7/2011
% function changes source value for FDTD simulation.  Ewave and media are
% structures with simulation parameters and n is the time step
% edited 8/15/2012 to solve for spherical wave exactly. No more paraxial
% approximation for Gaussian

persistent dr drp drpd % pulse information that will only be calculated once

% create the actual source
if n == 0,
    % create source weighting
%     st = size(media.ior_matrix,2);
    Ewave.source.indicies = media.BC.num+1:size(media.ior_matrix,2)-media.BC.num;
    st = length(Ewave.source.indicies);
    x = (round(Ewave.source.xpos/media.delSpatial)-[1:st]+media.BC.num)*media.delSpatial; % distance from the source
    
    [Ewave] = wavetype(Ewave,media,x);
    wt = 0; % weighting
    dr = []; % resets persistent values
% multiply the source by a weighting based off time
elseif n < Ewave.source.delay, % source is off
    wt = 0;
elseif (n >= Ewave.source.delay)&(n < dr), % rampup
    wt = erf(2*n/(dr-Ewave.source.delay));
elseif (n >= dr)&(n < drp), % steady state
    wt = 1;
elseif (n >= drp)&(n < drpd), % rampdown
    wt = erf(2*(drpd-n)/(drpd-drp));
else, % shut off
    wt = 0;
end;

%% set time parameters for pulse
if isempty(dr),
    dr = Ewave.source.delay+Ewave.source.rampup;
    drp = Ewave.source.delay+Ewave.source.rampup+Ewave.source.pulsewidth;
    drpd = Ewave.source.delay+Ewave.source.rampup+Ewave.source.pulsewidth+Ewave.source.rampdwn;
    Ewave.source.fin = round(drpd/media.delTemporal); % timesteps for pulse to finish
    Ewave.w = 2*pi*Ewave.source.c0/Ewave.lambda0; % angular frequency in freespace
end;

% update source with phase and amplitude information
phs = exp(1i*Ewave.w*n);
Ewave.source.instant = real(phs*wt*Ewave.source.cnst); 


function [Ewave] = wavetype(Ewave,media,x);
% [Ewave] = wavetype(Ewave,x);
% wavetype is a subfunction of sourceval and produces a 'plane', 'gaussian'
% or 'play back' wave
if strcmpi(Ewave.source.type,'gaussian'),
        [Ewave.source.cnst] = advancedgaus(Ewave,media);
        
elseif strcmpi(Ewave.source.type,'plane'),
        Ewave.source.cnst = zeros(size(x));
        Ewave.source.cnst(abs(x) < Ewave.source.beamwidth/2) = 1;
        
elseif strcmpi(Ewave.source.type,'play back'),
    Ewave.source.cnst = Ewave.source.conj; % needs to be defined somewhere
end;
        
function [val] = advancedgaus(Ewave,media);
% function was written 10/4/2011
% simulates a gaussian wave focused to any depth within the media
% Ewave.source.beamwidth - width of focal spot

if ~isfield(Ewave.source,'guess_n'),
    Ewave.source.guess_n = 1.33;
end;
% guess at the wavelength in the medium
lam = Ewave.lambda0/Ewave.source.guess_n; 

%% setup grid area indicies
% [x] = size(media.ior_matrix,2);
x = length(Ewave.source.indicies);
xca = (1:x).*media.delSpatial-Ewave.source.xpos;
xax = (1:x).*media.delSpatial;

% find offsets
dx = Ewave.source.xpos-Ewave.source.xwaist;
dz = Ewave.source.dpos-Ewave.source.dwaist;

w0 = Ewave.source.beamwidth/2;  % Radius of beam waist 
b = pi*w0^2/lam; % Rayleigh range or confocal parameter (9.21)
k = 2*pi/lam;  % wavenumber (k)
zh = sqrt(dx^2+dz^2); % distance from center of the source to focus
zhs = sign(dz)*zh;
% theta = atan2(dx,dz); % account for scanning
rho = zhs+b.^2/zhs; % Local radius of curvature (9.9)
psi = atan(zhs/b);  % Gouy phase (9.22)
zcc = Ewave.source.dpos+rho*dz/zh;
xcc = Ewave.source.xpos+rho*dx/zh;

if zcc ~= zcc, zcc = Ewave.source.dpos; end;
if xcc ~= xcc, xcc = Ewave.source.xpos; end;

% delz = xca.*tan(theta);

%% create beamf
bp = (dz^2+b^2)/b;  %  local beam diameter (9.9)        
w = sqrt(2*bp/k); % local beam radius (9.10)
val = sqrt(2/pi./w^2).*exp(1i*k*sqrt((xax-xcc).^2+(zcc-Ewave.source.dpos).^2)).*exp(-xca.^2./w.^2).*exp(-1i*psi); % Gaussian beam (9.17)  -- actual!

val = (val)./max(abs(val(:))); % normalize and return real component
% % % % enforce zeros in boundaries - ensures stability
% val([1:media.BC.num,end-media.BC.num]) = 0;
        
        
        