function [r, gd, gdd] = stackgddfast(ks, ds, n0, ns, dns, theta, pol)
% STACKGDDFAST  Calculate complex reflection, group delay and GDD of
% arbitrary dielectric structure using approximate analytic methods,
% including only first order material dispersion. Assumes two indices.
% This software is provided as an example. In practice, it is best to
% use one of the included MEX versions that are significantly faster.
%
% [r, gd, gdd] = stackgddfast(ks, ds, n0, ns, dns, ddns, theta, pol)
%
% Input:
% ks: wavenumber column vector (1/um)
% ds: thickness row vector (um)
% ns(length(ks),3): zeroth order index matrix (first layer, second, substrate)
% dns(length(ks),3): first order index (first layer, second, substrate)
% n0: initial index, assumed to be dispersionless
% theta: incidence angle (radians)
% pol: {'TE','TM'}
%
% Output:
% r(ks): complex reflection coefficient
% gd(ks): group delay (fs)
% gdd(ks): first order dispersion (fs^2)

% Implemenation Notes:
%
% Each T matrix goes from the end of the last layer through the Lth layer.
% Only the first two elements of the transfer matrix are kept
% at each point due to the symmetry of T matrices.
%
% All matrices are arranged with wavenumbers in columns.


%%% Initialization and input handling.
c = .2997924580;  % um/fs
if (nargin < 7) pol = 'TE'; end
if (nargin < 6) theta = 0; end
nk = length(ks);  % number of wavelengths
ds = [NaN, ds, 0];  % create dummy layers for outside and substrate
n = length(ds);  % number of layers
if (dns == 0)
    dns = zeros(nk, 3);
end

%%% Initialize loop variables.
T = [ones(nk,1), zeros(nk,1)];  % initial matrix
dT = zeros(nk,2);  % initial 1st derivative matrix
ddT = zeros(nk,2);  % initial 2nd derivative matrix

%%% Precalculate material variables.
% index variables: [n1, n2, nsub]
% p variables: [p01, p12, p21, p2sub]

% Effective indices.
neffs = sqrt(ns.^2 - (n0*sin(theta))^2);
neff0 = n0*cos(theta);
dneffs = ns.*dns./neffs;  % not really needed for nsub

% Calculate pTE (needed regardless of polarization).
pTEs(:,1) = neff0./neffs(:,1); % n1 from air
pTEs(:,2) = neffs(:,1)./neffs(:,2); % n2 from n1
pTEs(:,3) = 1./pTEs(:,2);  % n1 from n2
pTEs(:,4) = neffs(:,2)./neffs(:,3);  % nsub from n2
% Calculate p = pTM or p = pTE as needed.
if (pol == 'TM')
    % Calculate p = pTM from p0 and pTE.
    % p0 is just p (as above) for normal incidence.
    % p0 = ns(:,L - 1)./nL;
    % p = pTE./p0.^2;
    p0s(:,1) = n0./ns(:,1);
    p0s(:,2) = ns(:,1)./ns(:,2);
    p0s(:,3) = 1./p0s(:,2);
    p0s(:,4) = ns(:,2)./ns(:,3);
    ps = pTEs./p0s.^2;
else  % pol == 'TE'
    ps = pTEs;
end

% Common subexpressions.
dterms = (neffs + [ks,ks,ks].*dneffs);
dterms2 = dterms.^2;
kneffs = [ks,ks,ks].*neffs;
twoidneffs = 2i*dneffs;
pps = 1 + ps;
pms = 1 - ps;

%%% Step through structure, calculating all wavelengths in parallel.
for L = 2:n,
    % Initialize local variables.
    d = ds(L);
    id = 1i*d;  % layer thickness times i
    
    % Select appropriate material parameter indices
    if (L == n)  % last layer
        nx = 3; px = 4;
    elseif (L == 2)  % first layer
        nx = 1; px = 1;
    elseif (mod(L,2))  % even layer, n1->n2
        nx = 2; px = 2;
    else  % odd layer, n2->n1
        nx = 1; px = 3;
    end
    
    % Calculate new matrices for full structure up to L.
    % Only three matrix multiplys are required due to the approximation.
    ephi = exp(-id*kneffs(:,nx))/2;  % phasor of layer
    Dterm = [-id*dterms(:,nx)];
    D = [Dterm, Dterm];
    D2term = -d.*(twoidneffs(:,nx) + d*dterms2(:,nx));
    D2 = [D2term, D2term];
    TL1 = ephi.*pps(:,px); TL2 = ephi.*pms(:,px);  % current layer elements
    T = [TL1.*T(:,1) + TL2.*conj(T(:,2)), TL1.*T(:,2) + TL2.*conj(T(:,1))];  % new T
    TLdT = [TL1.*dT(:,1) + TL2.*conj(dT(:,2)), TL1.*dT(:,2) + TL2.*conj(dT(:,1))];
    dT = D.*T + TLdT;
    ddT = D2.*T + 2*D.*TLdT + ...
      [TL1.*ddT(:,1) + TL2.*conj(ddT(:,2)), TL1.*ddT(:,2) + TL2.*conj(ddT(:,1))];
end

% Calculate complex reflectance and group delay from total T matrices.
% [This is not as numerically stable as it could be?]
r = -T(:,2)./T(:,1);  % complex reflection coefficient
dr = (T(:,2).*dT(:,1) - dT(:,2).*T(:,1))./T(:,1).^2;
ddr = (T(:,2).*(-2*dT(:,1).^2 + T(:,1).*ddT(:,1)) + ...
    T(:,1).*(2*dT(:,1).*dT(:,2) - T(:,1).*ddT(:,2)))./T(:,1).^3;
dphi = (imag(dr).*real(r) - real(dr).*imag(r))./abs(r).^2;
dmag = (real(dr).*real(r) + imag(dr).*imag(r))./abs(r);
ddphi = (imag(ddr).*real(r) - real(ddr).*imag(r))./abs(r).^2 - ...
    2*dmag.*dphi./abs(r);
gd = dphi/c;
gdd = ddphi/c^2;