emgsafwd.m

function [result] = emgsafwd(varargin)

global debug

% [result] = emgsafwd(config, layers, cparams);
% 
% emgsafwd.m is part of the CR1Dmod forward modeling package, and contains
% the code used to calculate the response of electrode arrays on the surface
% of a 1D layered half space.
%
% Input:
% config :    structure defining the configuation (generated by cr1dmod) 
% layers :    structure defining layer-parameters (generated by cr1dmod)
% cparams:    structure defining calculation options (generated by compute)
% name:       optional name to display in the waitbar window
%
% Output:
% result:     structure containing the calculated responses as well
%                as information on the configuration, such as the 
%                geometric factor.
%
% Written by:
% Thomas Ingeman-Nielsen
% The Arctic Technology Center, BYG
% Technical University of Denmark
% Email: tin@byg.dtu.dk

cpu_t=cputime;

global Rkern estimate
Rkern = [];
estimate = [];

if isappdata(0,'debug') && getappdata(0,'debug')
    dbstop if all error
    debug = 1;
else
    debug = 0;
end

if ( nargin < 3 )
    disp('A call to emgsafwd must provide three input arguments!');
    return
end

if ( nargin == 4 ) && ~isempty(varargin{4})
    modelname = [varargin{4} ':  '];
else 
    modelname = [];
end

% Set up structures of parameters...
ok = 0;
if ( isstruct(varargin{1})         )
    config = varargin{1};
    if ( isstruct(varargin{2})     )
        layers = varargin{2};
        if ( isstruct(varargin{3}) )
            cparams = varargin{3};                    
            ok = 1;
        end
    end
end
if ~ok
    disp('Non-structure inputs have not been implemented yet!');
    return
end

if ~isfield(cparams,'noPrint') || ~cparams.noPrint
    disp('Routine: EMgsafwd');
end    

% Set up vectors of layer parameters
eps0 = 8.85418782.*1e-12;
mu0 = 4e-7.*pi;

lparams.h = [layers.thickness].'; lparams.h(end) = []; % only include layers above halfspace

omega = cparams.freq.*2.*pi;
sigma = 1./Z_CR(omega, [layers.rho].', [layers.tau].', [layers.c].', ...
    [layers.m].');

if layers(1).mu == 0
    cparams.nonmag = 1;
else
    cparams.nonmag = 0;
end

if ~isfield(cparams,'noDialog') || ~cparams.noDialog
    h1 = waitbar(0,[modelname 'Please wait...']);
    noDialog = 0;
else
    noDialog = 1;
end

% setup config structure for dipole-dipole calculations
if strcmpi(config(1).type, 'Dipole-Dipole')
    for k = 1:length(config(1,1).Aspac)
        for m = 1:length(config(1,1).Nspac)
            config(k,m) = config(1,1);
            config(k,m).C1 = [-config(k,m).Aspac(k)-...
                    0.5*config(k,m).Aspac(k)*config(k,m).Nspac(m) 0 0]; % transmitter dipole            C1  a  C2     P2  a  P1
            config(k,m).C2 = [config(k,m).C1(1)+...
                    config(k,m).Aspac(1) 0 0];                          % transmitter dipole            o------o      o------o    
            config(k,m).P1 = [config(k,m).C2(1)+...                     % receiver dipole                 Tx     n*a    Rx
                    config(k,m).Nspac(m)*config(k,m).Aspac(k) 0 0];   
            config(k,m).P2 = [config(k,m).P1(1)+...                     % receiver dipole
                    config(k,m).Aspac(k) 0 0];

            % Here we interchange C1 and C2, remember then to also change
            % cosine calculation
            C2 = config(k,m).C2;
            config(k,m).C2 = config(k,m).C1;
            config(k,m).C1 = C2;
            
            config(k,m).Cwire = [config(k,m).C1;config(k,m).C2];
            % Calculate length of each wire element
            config(k,m).S = sqrt(sum(diff(config(k,m).Cwire).^2,2));    
            
            config(k,m).Pwire = [config(k,m).P1;config(k,m).P2];
            % Calculate length of each wire element
            config(k,m).T = sqrt(sum(diff(config(k,m).Pwire).^2,2));
            
            %config(k,m).cosines = 1;
            config(k,m).cosines = (diff(config(k,m).Cwire)*diff(config(k,m).Pwire).')       ...
                ./(abs(config(k,m).S)*abs(config(k,m).T).');  % calculate cosines of angles between wire elements
            
            config(k,m).R = [];

            for s = 1:length(config(k,m).S)
                for t = 1:length(config(k,m).T)
                    config(k,m).R = [config(k,m).R; sqrt(sum(        ...
                            (config(k,m).Cwire(s:s+1,:)-             ...
                            config(k,m).Pwire(t:t+1,:)).^2,2))];
                    config(k,m).R = [config(k,m).R; sqrt(sum(        ...
                            (config(k,m).Cwire(s:s+1,:)-             ...
                            config(k,m).Pwire([t+1,t],:)).^2,2))];
                end
            end
            config(k,m).Rmin = min(config(k,m).R);
            config(k,m).Rmax = max(config(k,m).R);
            
            % calculate distance between electrodes
            config(k,m).r = sqrt(sum(([...
                    config(k,m).P2-config(k,m).C2;...           
                    config(k,m).P2-config(k,m).C1;...
                    config(k,m).P1-config(k,m).C2;...
                    config(k,m).P1-config(k,m).C1]).^2,2));
           
            % prepare result structure
            result(k,m).Aspac = config(k,m).Aspac(k);
            result(k,m).Nspac = config(k,m).Nspac(m);
            result(k,m).G_factor = 2*pi*(...
                1/config(k,m).r(4)-      ...
                1/config(k,m).r(3)-      ...
                1/config(k,m).r(2)+      ...
                1/config(k,m).r(1))^(-1);
            
            result(k,m).Z = omega.';  %allocate space for Z
        end % m
    end % k
    cparams.Rmin = min(min([config.Rmin]));
    cparams.Rmax = max(max([config.Rmax]));
    
else % if it is a general surface array
    for k = 1:length(config(:))
        config(k).Cwire = [config(k).C1;config(k).Cwire;config(k).C2];
        config(k).S = sqrt(sum(diff(config(k).Cwire).^2,2));     % Calculate length of each wire element
        
        config(k).Pwire = [config(k).P1;config(k).Pwire;config(k).P2];
        config(k).T = sqrt(sum(diff(config(k).Pwire).^2,2));     % Calculate length of each wire element
        
        config(k).cosines = (diff(config(k).Cwire)*diff(config(k).Pwire).')       ...
            ./(abs(config(k).S)*abs(config(k).T).');  % calculate cosines of angles between wire elements
        
        config(k).R = [];
        for s = 1:length(config(k).S)
            for t = 1:length(config(k).T)
                if config(k).cosines(s,t) ~=0
                    config(k).R = [config(k).R; sqrt(sum((config(k).Cwire(s:s+1,:)...
                            -config(k).Pwire(t:t+1,:)).^2,2))];
                    config(k).R = [config(k).R; sqrt(sum((config(k).Cwire(s:s+1,:)...
                            -config(k).Pwire([t+1,t],:)).^2,2))];
                end % if   
            end % t
        end % s
        
        % calculate distance between electrodes    
        config(k).r = sqrt(sum(([config(k).P2-config(k).C2; config(k).P2-config(k).C1;     ...
                config(k).P1-config(k).C2; config(k).P1-config(k).C1]).^2,2));
        
        % prepare result structure
        result(k).G_factor = 2*pi*(1/config(k).r(1)-1/config(k).r(2)-1/config(k).r(3)+  ...
            1/config(k).r(4))^(-1);
        result(k).Z = omega.';  %allocate space for Z
    end
    cparams.Rmin = min(min([config.R]));
    cparams.Rmax = max(max([config.R]));    
end % if

if ( ~cparams.Rspline )
    % In the case no spline interpolation is wanted
    %     return handle to FD calculation routine
    estimate = [];                             
    switch cparams.domain
        case {'FD'}
                    P_handle = @P_fun;
                    Q_handle = @Q_fun;
        case {'TD'}
            disp('Not implemented yet!');
    end % switch
elseif ~isempty(cparams.Rmin) % If there is coupling
    % If spline interpolation is wanted
    % setup some parameters and return handle to spline evaluation 

    Rlim = (floor(log2([cparams.Rmin cparams.Rmax]).*cparams.Rsp_NDEC+100)-[102 97])./cparams.Rsp_NDEC;
    cparams.R = 2.^(Rlim(1):1./cparams.Rsp_NDEC:Rlim(2));
    
    if debug
        disp('Evaluation distances for the creation of the spline function:');
        disp(cparams.R);
        assignin('base','R_values',cparams.R);
    end
    
    P_handle = @eval_Rspline;
    Q_handle = @Q_fun;
else
    Q_handle = @Q_fun;
end

original_Quad_tol = cparams.Quad_tol;

lo = length(omega);  % number of frequencies
lc = prod(size(config)); % number of configurations
lt = lo.*lc; % total number of calculations
for k = 1:lo        % cycle over frequencies 
    lparams.z = i.*omega(k).*mu0.*([layers.mu].'+1);
    if strcmp(cparams.calc_type,'Quasi')
        % setup free-space parameters
        lparams.z0 = i.*omega(k).*mu0;
        lparams.y0 = 0;
        lparams.k0_sq = 0;
        % Permittivity is neglected in the quasi-static half-space response
        lparams.y = sigma(:,k);
        % Non-magnetic first layer is assumed
        lparams.z(1,:) = lparams.z0(1,:);
    else
        % Use the full solution
        lparams.z0 = i.*omega(k).*mu0;
        lparams.y0 = i.*omega(k).*eps0;
        lparams.k0_sq = -lparams.z0.*lparams.y0;
        lparams.y = sigma(:,k) + i.*eps0.*[layers.eps_r].'*omega(k);
    end
    lparams.k_sq = -lparams.z.*lparams.y;
    
    if ( cparams.Rspline && isfield(cparams, 'R'))
        gsa_driver(cparams, omega(k), lparams);
    end
    
    
    for m = 1:lc    % cycle over configurations

        r_un = unique(config(m).r);
               
        drawnow;

        P_term = 0;
        
        if ~all(isequal(config(m).cosines, 0))
            % tolerance for each segment integration should add up to Quad_tol
            cparams.Quad_tol = original_Quad_tol/                         ...
                sum(sum(config(m).cosines&1));            

            total_num = length(config(m).S).*length(config(m).T);
            for s = 1:length(config(m).S)
                for t = 1:length(config(m).T)
                    if config(m).cosines(s,t) ~=0
                        if debug
                            disp(['Integration pass:  ' ...
                                    num2str((s-1)*length(config(m).T)+t) ...
                                    ' of ' num2str(total_num)]);
                        end
                        P_term = P_term + ...
                            lparams.z0./4./pi.*            ... % removed from the P-function and added here!    
                            dblquad(P_handle, 0, config(m).T(t), 0,      ...
                            config(m).S(s), cparams.Quad_tol, @quad,     ...
                            config(m).Pwire(t:t+1,:),                    ...
                            config(m).Cwire(s:s+1,:),                    ...
                            config(m).T(t), config(m).S(s),              ... 
                            omega(k), lparams, cparams).*                ...
                            config(m).cosines(s,t);           
                    end
                end
            end
        end
        if debug
            disp('Evaluating Q-function');
        end
        Q = feval(Q_handle, r_un, omega(k), lparams, cparams);   
        
        result(m).Z(k) = (P_term +                                 ...
            Q(config(m).r(1)==r_un)-Q(config(m).r(2)==r_un)-         ...
            Q(config(m).r(3)==r_un)+Q(config(m).r(4)==r_un));

        if ~noDialog
            titlestr = [modelname 'Just did: \omega = ' num2str(round(omega(k)     ...
                    *1000)/1000) '  Config: ' num2str(m)];  % 'Hz, Aspac = ' ...
            %            num2str(Aspac(m)) 'm, Nspac = ' sprintf('%4.2f',Nspac(n))];
            
            waitbar(((k-1)*lc+m)/lt, h1, titlestr);  
        end
    end % m    (config)
end % k        (omega)

if ~noDialog
    close(h1);
end

if ~isfield(cparams,'noPrint') || ~cparams.noPrint
    disp(['Time spent: ' num2str(cputime-cpu_t) ' s']);
end

%****************************************************
% Subfunctions
%****************************************************

% -------------------------------------------------------------------------
function gsa_driver(cparams, omega, lparams)

global Rkern debug

if debug
    disp(['GSA driver: ' num2str(omega)])
end

switch cparams.domain
    case {'FD'}
        if strcmp(cparams.calc_type,'Quasi')
            % Quasi static Halfspace response assumes also nonmagnetic first layer
            %disp('quasi+nonmag')
            r = cparams.R;
            kern = 2./(lparams.k_sq(1).*r.^3).*(                ...
                (i.*lparams.k_sq(1).^0.5.*r+1).*                ...
                exp(-i.*lparams.k_sq(1).^0.5.*r) - 1);
            if ~isempty(lparams.h)
                switch cparams.hank_type
                    case {'NHT'} 
                        for k = 1:length(r)
                            kern(k) = kern(k) + ...% j.*omega.*4.*pi.*1e-7./4./pi.*...
                                NJCST('J0', @P_kernel_quasinonmag,      ...
                                cparams.R(k), cparams.Seg_tol,          ...
                                cparams.NHT_tol, cparams.Max_seg,       ...
                                omega, lparams);
                        end
                    case {'FHT'}
                        for k = 1:length(r)
                            kern(k) = kern(k) + ...% j.*omega.*4.*pi.*1e-7./4./pi.*...
                                FJCST('J04', cparams.R(k), -1, 0, 1,    ...
                                cparams.FHT_err, @P_kernel_quasinonmag, ...
                                omega, lparams);
                        end
                end
            end
        elseif cparams.nonmag
            % nonmagnetic layered halfspace
            %disp('nonmag');
            r = cparams.R;
            kern = 2./((lparams.k_sq(1)-lparams.k0_sq).*r.^3).* ...
                ((i.*lparams.k_sq(1).^0.5.*r+1).*               ...
                exp(-i.*lparams.k_sq(1).^0.5.*r) -              ...
                (i.*lparams.k0_sq.^0.5.*r+1).*                  ...
                exp(-i.*lparams.k0_sq.^0.5.*r));
            if ~isempty(lparams.h)
                switch cparams.hank_type
                    case {'NHT'}
                        for k = 1:length(r)
                            kern(k) = kern(k) + ...% j.*omega.*4.*pi.*1e-7./4./pi.*...
                                NJCST('J0', @P_kernel_nonmag,           ...
                                cparams.R(k), cparams.Seg_tol,          ...
                                cparams.NHT_tol, cparams.Max_seg,       ...
                                omega, lparams);
                        end
                    case {'FHT'}
                        for k = 1:length(r)
                            kern(k) = kern(k) + ...% j.*omega.*4.*pi.*1e-7./4./pi.*...
                                FJCST('J04', cparams.R(k), -1, 0, 1,    ...
                                cparams.FHT_err, @P_kernel_nonmag,      ...
                                omega, lparams);
                        end
                end
            end                    
        else
            % full solution
            %disp('full solution');
            switch cparams.hank_type
                case {'NHT'}
                    for k = 1:length(cparams.R)
                        kern(k) = ... % j.*omega.*4.*pi.*1e-7./4./pi.*...
                            NJCST('J0', @P_kernel_full, cparams.R(k),   ...
                            cparams.Seg_tol, cparams.NHT_tol,           ...
                            cparams.Max_seg, omega, lparams);
                    end
                case {'FHT'}
                    for k = 1:length(cparams.R)
                        kern(k) = ... % j.*omega.*4.*pi.*1e-7./4./pi.*...
                            FJCST('J04', cparams.R(k), -1, 0, 1,        ...
                            cparams.FHT_err, @P_kernel_full, omega,     ...
                            lparams);
                    end
            end
        end
    case {'TD'}
        disp('Not implemented yet!');
end % switch

strct = spline(cparams.R,kern);

strct.R = cparams.R;
strct.kern = kern;
strct.estimate = [];

if ~isempty(Rkern)
%    Rkern(end+1) = strct;    % This line is used to save the splines for
%    later review
    Rkern = strct;
else
    Rkern = strct;    
end
if debug
    disp('End of GSA driver');
end
% -------------------------------------------------------------------------
function kern = eval_Rspline(y, x, Pwire, Cwire, T, S, omega,           ...
    lparams, cparams)

global Rkern

P = repmat(diff(Pwire)./T,length(y),1)...
    .*repmat(y.', 1, 3) + repmat(Pwire(1,:),length(y),1);  % Pwire-element coordinates (will be a vector)
C = diff(Cwire)./S.*x + Cwire(1,:);                       % Cwire-element coordinates (will be a scalar)

r = sqrt(sum((P-repmat(C,length(y),1)).^2,2));            % distances are calculated

kern = ppval(r, Rkern(end));
%Rkern(end).estimate = [Rkern(end).estimate; r]; % use this line to save
%r-values of estimation points.

% -------------------------------------------------------------------------
function P = P_fun(y, x, Pwire, Cwire, T, S, omega, lparams, cparams)

P = repmat(diff(Pwire)./T,length(y),1)...
    .*repmat(y.', 1, 3) + repmat(Pwire(1,:),length(y),1);  % Pwire-element coordinates (will be a vector)
C = diff(Cwire)./S.*x + Cwire(1,:);                       % Cwire-element coordinates (will be a scalar)

r = sqrt(sum((P-repmat(C,length(y),1)).^2,2));            % distances are calculated

% cparams.calc_type = 'Quasi';

if strcmp(cparams.calc_type,'Quasi')
    % Quasi static Halfspace response assumes also nonmagnetic first layer
    %disp('quasi+nonmag')
    P = 2./(lparams.k_sq(1).*r.^3).*((i.*lparams.k_sq(1).^0.5.*r+1).*   ...
        exp(-i.*lparams.k_sq(1).^0.5.*r) - 1);
    if ~isempty(lparams.h)
        switch cparams.hank_type
            case {'FHT'}
                for k = 1:length(r)
                    P(k) = P(k) +... % j.*omega.*4.*pi.*1e-7./4./pi.*...
                        FJCST('J04', r(k), -1, 0, 1, cparams.FHT_err, ...
                        @P_kernel_quasinonmag, omega, lparams);
                end
            case {'NHT'}
                for k = 1:length(r)
                    P(k) = P(k) +... % j.*omega.*4.*pi.*1e-7./4./pi.*...
                        NJCST('J0', @P_kernel_quasinonmag, r(k),         ...
                        cparams.Seg_tol, cparams.NHT_tol,               ...
                        cparams.Max_seg, omega, lparams);
                end
        end
    end
elseif cparams.nonmag
    % nonmagnetic layered halfspace
    %disp('nonmag');
    P = 2./((lparams.k_sq(1)-lparams.k0_sq).*r.^3).*(                   ...
        (i.*lparams.k_sq(1).^0.5.*r+1).*exp(-i.*lparams.k_sq(1).^0.5.*r)...
        - (i.*lparams.k0_sq.^0.5.*r+1).*exp(-i.*lparams.k0_sq.^0.5.*r));
    if ~isempty(lparams.h)
        switch cparams.hank_type
            case {'FHT'}
                for k = 1:length(r)
                    P(k) = P(k) +... % j.*omega.*4.*pi.*1e-7./4./pi.*...
                        FJCST('J04', r(k), -1, 0, 1, cparams.FHT_err, ...
                        @P_kernel_full, omega, lparams);
                end
            case {'NHT'}
                for k = 1:length(r)
                    P(k) = P(k) +... % j.*omega.*4.*pi.*1e-7./4./pi.*...
                        NJCST('J0', @P_kernel_nonmag, r(k),              ...
                        cparams.Seg_tol, cparams.NHT_tol,               ...
                        cparams.Max_seg, omega, lparams);
                end
        end
    end
else
    % full solution
    %disp('fullsol');
    switch cparams.hank_type
        case {'FHT'}
            for k = 1:length(r)
                P(k) = ... % j.*omega.*4.*pi.*1e-7./4./pi.*...
                    FJCST('J04', r(k), -1, 0, 1, cparams.FHT_err,     ...
                    @P_kernel_full, omega, lparams);
            end
        case {'NHT'}
            for k = 1:length(r)
                P(k) = ... % j.*omega.*4.*pi.*1e-7./4./pi.*...
                    NJCST('J0', @P_kernel_full, r(k), cparams.Seg_tol,   ...
                    cparams.NHT_tol, cparams.Max_seg, omega, lparams);
            end
    end
end

%estimate = [estimate; r, P];

% -------------------------------------------------------------------------
function Q = Q_fun(r, omega, lparams, cparams)

if strcmp(cparams.calc_type,'Quasi')
    % Quasi static Halfspace response assumes also nonmagnetic first layer
    %disp('quasi+nonmag')
    Q = -lparams.z0./(2.*pi.*lparams.k_sq(1,:).*r);
    if ~isempty(lparams.h)
        switch cparams.hank_type
            case {'FHT'}
                for k = 1:length(r)
                    Q(k) = -1./4./pi.*... 
                        (FJCST('J04', r(k), -1, 0, 1, cparams.FHT_err,        ...
                        @Q_kernel_quasinonmag, omega, lparams)) + Q(k);
                end
            case {'NHT'}
                for k = 1:length(r)
                    Q(k) = -1./4./pi.*... 
                        (NJCST('J0', @Q_kernel_quasinonmag, r(k),                ...
                        cparams.Seg_tol, cparams.NHT_tol, cparams.Max_seg,      ...
                        omega, lparams)) + Q(k);
                end
        end
    end
else
    % full solution
    %disp('full');
    switch cparams.hank_type
        case {'FHT'}
            for k = 1:length(r)
                Q(k) = 1./4./pi.*... 
                    (FJCST('J04', r(k), -1, 0, 1, cparams.FHT_err,            ...
                    @Q_kernel_full, omega, lparams));
            end
        case {'NHT'}
            for k = 1:length(r)
                Q(k) = 1./4./pi.*... 
                    (NJCST('J0', @Q_kernel_full, r(k), cparams.Seg_tol,          ...
                    cparams.NHT_tol, cparams.Max_seg, omega, lparams));
            end
    end
end


% -------------------------------------------------------------------------
function kern = P_kernel_full(lambda, varargin)

% kern = P_kernel(lambda, omega, eps, mu, sigma, h)
%
% Kernelfunction:       (1+R_TE).*lambda./u0
%
% This version is for the full solution!

%omega =   varargin{1};
lparams = varargin{2};

lambda_sq = lambda.^2;
u0 = sqrt(lambda_sq-lparams.k0_sq);
u = sqrt(repmat(lambda_sq,size(lparams.k_sq,1),1)-                      ...
    repmat(lparams.k_sq, 1, size(lambda_sq,2)));
Y = u./repmat(lparams.z,1,size(lambda,2));
Y0 = u0./lparams.z0;

gamma = 0;
if ~isempty(lparams.h)
    expuh2 = exp(-2.*u(1:end-1,:).*repmat(lparams.h,1,length(lambda)));
    phi = (Y(1:end-1,:)-Y(2:end,:))./(Y(1:end-1,:)+Y(2:end,:)); %  = (Yn-Yn+1)/(Yn+Yn+1)
    
    for m = length(lparams.h):-1:1
        gamma = expuh2(m,:).*(gamma+phi(m,:))./(gamma.*phi(m,:)+1);
    end
end

phi1 = (Y0-Y(1,:))./(Y0+Y(1,:));    %  as needed in the numerator
phi1mod = (2./lparams.z0)./(Y0+Y(1,:));     %  as needed in the denominator

kern = ((gamma+1).*phi1mod.*lambda)./(gamma.*phi1+1); % = ((R_TE+1)/u0) *lambda


% -------------------------------------------------------------------------
function kern = P_kernel_nonmag(lambda, varargin)

% kern = P_kernel(lambda, omega, eps, sigma, h)
%
% Kernelfunction:       (1+R_TE).*lambda./u0
%
% This version is for the nonmagnetic ground case!

%omega =   varargin{1};
lparams = varargin{2};

lambda_sq = lambda.^2;
u0 = sqrt(lambda_sq-lparams.k0_sq);
u = sqrt(repmat(lambda_sq,size(lparams.k_sq,1),1)-                      ...
    repmat(lparams.k_sq, 1, size(lambda_sq,2)));
Y = u./repmat(lparams.z,1,size(lambda,2));
%Y0 = u0./lparams.z0;

gamma = 0;
if ~isempty(lparams.h)
    expuh2 = exp(-2.*u(1:end-1,:).*repmat(lparams.h,1,length(lambda)));
    phi = (Y(1:end-1,:)-Y(2:end,:))./(Y(1:end-1,:)+Y(2:end,:)); %  = (Yn-Yn+1)/(Yn+Yn+1)
    
    for m = length(lparams.h):-1:1
        gamma = expuh2(m,:).*(gamma+phi(m,:))./(gamma.*phi(m,:)+1);
    end
end

kern = (4.*gamma.*u(1,:).*lambda)./(gamma.*(lparams.k_sq(1,:)-          ...
    lparams.k0_sq)+(u0+u(1,:)).^2);
% = layered-homhalf

% -------------------------------------------------------------------------
function kern = P_kernel_quasinonmag(lambda, varargin)

% kern = P_kernel(lambda, omega, eps, sigma, h)
%
% Kernelfunction:       (1+R_TE).*lambda./u0
%
% This version is for the quasistatic nonmagnetic first layer case!

%omega =   varargin{1};
lparams = varargin{2};

lambda_sq = lambda.^2;
u = sqrt(repmat(lambda_sq,size(lparams.k_sq,1),1)-                      ...
    repmat(lparams.k_sq, 1, size(lambda_sq,2)));
Y = u./repmat(lparams.z,1,size(lambda,2));
%Y0 = lambda./lparams.z0;

gamma = 0;
if ~isempty(lparams.h)
    expuh2 = exp(-2.*u(1:end-1,:).*repmat(lparams.h,1,length(lambda)));
    phi = (Y(1:end-1,:)-Y(2:end,:))./(Y(1:end-1,:)+Y(2:end,:)); %  = (Yn-Yn+1)/(Yn+Yn+1)
    for m = length(lparams.h):-1:1
        gamma = expuh2(m,:).*(gamma+phi(m,:))./(gamma.*phi(m,:)+1);
    end
end

kern = (4.*gamma.*u(1,:).*lambda)./(gamma.*lparams.k_sq(1,:)+           ...
    (lambda+u(1,:)).^2);
% = layered-homhalf

% -------------------------------------------------------------------------
function kern = Q_kernel_full(lambda, varargin)

% kern = Q_kernel(lambda, omega, eps, mu, sigma, h)
%
% Kernelfunction:       ((1-R_TM)*u0/y0 - (1+r_TE)*z0/u0)/lambda

%omega =   varargin{1};  
lparams = varargin{2};

lambda_sq = lambda.^2;
u0 = sqrt(lambda_sq-lparams.k0_sq);
u = sqrt(repmat(lambda_sq,size(lparams.k_sq,1),1)-                      ...
    repmat(lparams.k_sq, 1, size(lambda_sq,2)));
Y = u./repmat(lparams.z,1,size(lambda,2));
Z = u./repmat(lparams.y,1,size(lambda,2));
Y0 = u0./lparams.z0;
Z0 = u0./lparams.y0;

gamma_TE = 0;
gamma_TM = 0;

if ~isempty(lparams.h)
    expuh2 = exp(-2.*u(1:end-1,:).*repmat(lparams.h,1,length(lambda)));
    phi_TE = (Y(1:end-1,:)-Y(2:end,:))./(Y(1:end-1,:)+Y(2:end,:)); %  = (Yn-Yn+1)/(Yn+Yn+1)
    phi_TM = (Z(1:end-1,:)-Z(2:end,:))./(Z(1:end-1,:)+Z(2:end,:)); %  = (Zn-Zn+1)/(Zn+Zn+1)
    for m = length(lparams.h):-1:1
        gamma_TE = expuh2(m,:).*(gamma_TE+phi_TE(m,:))./(gamma_TE.*phi_TE(m,:)+1);
        gamma_TM = expuh2(m,:).*(gamma_TM+phi_TM(m,:))./(gamma_TM.*phi_TM(m,:)+1);    
    end
end

phi1_TE = (Y0-Y(1,:))./(Y0+Y(1,:)); % as needed in the numerator
phi1mod_TE = 2./(Y0+Y(1,:));        % as needed in the denominator

phi1_TM = (Z0-Z(1,:))./(Z0+Z(1,:)); % as needed in the numerator
phi1mod_TM = -2.*Z(1,:).*Z0./(Z0+Z(1,:));        % as needed in the denominator

kern_TE = ((gamma_TE+1).*phi1mod_TE)./(gamma_TE.*phi1_TE + 1); 
kern_TM = ((gamma_TM-1).*phi1mod_TM)./(gamma_TM.*phi1_TM + 1);

kern = (kern_TM-kern_TE)./lambda;   


% -------------------------------------------------------------------------
function kern = Q_kernel_quasinonmag(lambda, varargin)

% kern = Q_kernel(lambda, omega, eps, mu, sigma, h)
%
% Kernelfunction:       ((1-R_TM)*u0/y0 - (1+r_TE)*z0/u0)/lambda

%omega =   varargin{1};  
lparams = varargin{2};

%lambda = lambda.*induct;
lambda_sq = lambda.^2;
u = sqrt(repmat(lambda_sq,size(lparams.k_sq,1),1)-                      ...
    repmat(lparams.k_sq, 1, size(lambda_sq,2)));
Y = u./repmat(lparams.z,1,size(lambda,2));
Z = u./repmat(lparams.y,1,size(lambda,2));
%Y0 = lambda./lparams.z0;

gamma_TE = 0;
gamma_TM = 0;

if ~isempty(lparams.h)
    expuh2 = exp(-2.*u(1:end-1,:).*repmat(lparams.h,1,length(lambda)));
    phi_TE = (Y(1:end-1,:)-Y(2:end,:))./(Y(1:end-1,:)+Y(2:end,:)); %  = (Yn-Yn+1)/(Yn+Yn+1)
    phi_TM = (Z(1:end-1,:)-Z(2:end,:))./(Z(1:end-1,:)+Z(2:end,:)); %  = (Zn-Zn+1)/(Zn+Zn+1)
    for m = length(lparams.h):-1:1
        gamma_TE = expuh2(m,:).*(gamma_TE+phi_TE(m,:))./(gamma_TE.*phi_TE(m,:)+1);
        gamma_TM = expuh2(m,:).*(gamma_TM+phi_TM(m,:))./(gamma_TM.*phi_TM(m,:)+1);    
    end
end

kern_TE = gamma_TE.*lparams.z0./(gamma_TE.*lparams.k_sq(1,:) +          ...
    (lambda + u(1,:)).^2);
kern_TM = gamma_TM./((gamma_TM+1).*lparams.y(1,:));

kern = (kern_TM+kern_TE).*4.*u(1,:)./lambda;