function chi = sscptker(varargin)
% CHI = SSCPTKER(P,K,D,H1,ETA1,ETA2,R1,R2,RHO1,RHO2) Calculate the complex 
%       susceptibilities of supported polymeric bilyer film.
% 
% INPUT ARGUMENTS:
%       P (1/s)             : Laplace frequency
%       K (1/nm)            : in-plane wave vector
%       D (nm)              : total film thickness
%       H1 (nm)             : bottom layer thickness
%       ETA1 (kg/(s*nm))    : dyanmic viscosity of bottom layer
%       ETA2 (kg/(s*nm))    : dynamic viscosity of top layer
%       R1 (kg/s^2)         : interfacial tension
%       R2 (kg/s^2)         : surface tension
%       RHO1 (kg/nm^3)      : mass density of bottom layer
%       RHO2 (kg/nm^3)      : mass density of top layer
%
% OUTPUT ARGUMENTS:
%       CHI : 2x2 matrix, complex susceptibilities.
%
% Note:
%       1. Every parameter is in SI unit, except that the length is in nm.

if nargin ~=10
    error('Invalid input arguments!');
    return;
end

% --- get input arguments
w       = varargin{1}*i;
k       = varargin{2}; 
d       = varargin{3}; 
h1      = varargin{4}; 
eta1    = varargin{5}; 
eta2    = varargin{6}; 
r1      = varargin{7}; 
r2      = varargin{8}; 
rho1    = varargin{9}; 
rho2    = varargin{10};

% --- define some constant and matrices
g = 0*1e9;        % gravitational field is negligible for thin films
k1  = sqrt(k^2 - i*w*rho1/eta1);            % kappa1 in layer 1
k2  = sqrt(k^2 - i*w*rho2/eta2);            % kappa2 in layer 2
B   = [...
    k*sinh(k*h1)-k1*sinh(k1*h1),                            k1*cosh(k*h1)-k1*cosh(k1*h1);...
    -i*k*cosh(k*h1)+i*k*cosh(k1*h1),                        -i*k1*sinh(k*h1)+i*k*sinh(k1*h1);...
    w*sinh(k*h1),                                           w*(k1/k)*cosh(k*h1);...
    2*eta1*k^2*cosh(k*h1)-eta1*(k^2+k1^2)*cosh(k1*h1),      2*eta1*k*k1*sinh(k*h1)-eta1*(k^2+k1^2)*sinh(k1*h1);...
    ];                                      % B(4x2) at interface z=h1
AInv=[...
    0,                  i*eta2/(rho2*w)*(k^2+k2^2)*cosh(k*h1),   i*k/w*sinh(k*h1),     -k/(rho2*w)*cosh(k*h1);...
    0,                 -i*eta2/(rho2*w)*(k^2+k2^2)*sinh(k*h1),  -i*k/w*cosh(k*h1),      k/(rho2*w)*sinh(k*h1);...
    sinh(k2*h1)/k2,    -2*eta2/(rho2*w)*k*cosh(k2*h1),          -k/(k2*w)*sinh(k2*h1), -i*cosh(k2*h1)/(rho2*w);...
   -cosh(k2*h1)/k2,     2*eta2/(rho2*w)*k*sinh(k2*h1),           k/(k2*w)*cosh(k2*h1),  i*sinh(k2*h1)/(rho2*w);...
    ];                                      % inverse of A(4x4) at z=h1

%======================================================
% --- interfacial boundary condition at z=h1 ---
%======================================================
T   = [...
    1,                          0,                                      0,              0;...
    0,                          1,                                      0,              0;...
    (2*i*k/rho2)*(eta1-eta2),   -(i/(rho2*w))*(r1*k^2+g*(rho1-rho2)),   rho1/rho2,      0;...
    0,                          0,                                      0,              1;...
    ];                                      % T(4x4)
if r1 == 0
    T(3,2) = -(i/(rho2*w))*g*(rho1-rho2);
end
G   = T*B;                                  % G(4x2)
C   = AInv*G;                               % intermdiate matrix

%======================================================
% --- surface boudnary conditions at z=d ---
%======================================================
D   = [...
    2*k*cosh(k*d),...
    2*k*sinh(k*d),...
    i*(k^2+k2^2)*cosh(k2*d),...
    i*(k^2+k2^2)*sinh(k2*d);...
    ];                                      % D(1x4)
P   = D*C;                                  % P(1x2)
P0  = D*AInv(:,3);                          % P0(1x1)
E   = [...
    (-i*w/k+2*eta2/rho2*k)*sinh(k*d) + i/w*(r2*k^2/rho2+g)*cosh(k*d),...
    (-i*w/k+2*eta2/rho2*k)*cosh(k*d) + i/w*(r2*k^2/rho2+g)*sinh(k*d),...
    2*i*eta2/rho2*k*k2*sinh(k2*d)    - k/w*(r2*k^2/rho2+g)*cosh(k2*d),...
    2*i*eta2/rho2*k*k2*cosh(k2*d)    - k/w*(r2*k^2/rho2+g)*sinh(k2*d);...
    ];                                      % E(1x4)
Q   = E*C;                                  % Q(1x2)
Q0  = E*AInv(:,3);                          % Q0(1x1)

%======================================================
% --- construct other matrices ---
%======================================================
R   = [...
    cosh(k*d),...
    sinh(k*d),...
    i*k*cosh(k2*d),...
    i*k*sinh(k2*d);...
    ];                                      % R(1x4)
M   = R*C;                                  % M(1x2)
M0  = R*AInv(:,3);                          % M0(1x1)
N   = [...
    -i*(k*cosh(k*h1)-k*cosh(k1*h1)),...
    -i*(k1*sinh(k*h1)-k*sinh(k1*h1));...
    ];                                      % N(1x2)

%======================================================
% --- calculate susceptibility matrix ---
%======================================================
chi(1,1)    = P0*(Q(1)*N(2)-Q(2)*N(1)) - Q0*(P(1)*N(2)-P(2)*N(1));
chi(1,2)    = P(1)*N(2) -  P(2)*N(1);
chi(2,1)    = M0*(P(1)*Q(2)-P(2)*Q(1)) + M(1)*(P(2)*Q0-P0*Q(2)) + M(2)*(P0*Q(1)-P(1)*Q0);
chi(2,2)    = P(1)*M(2)-P(2)*M(1);
chi         = chi*(i/(rho2*w)) / (P(1)*Q(2)-P(2)*Q(1));
return;
% --- EOF ---