dynamic_stiff(varargin)

function [D] = dynamic_stiff(varargin)
%dynamic_stiff   Calculate the dynamic stiffness of beam in space using
%spectral element
%
%   Function D = dynamic_stiff((w,L,E,rho,S,Izz,Iyy,Ixx,J,G,eta,ni)
% 
%   Function calculate the dynamic stiffness matrix for a beam in space
%   using the exact solution of the beam wave equations
%   D is the dynamic stiffness of the beam with dimension 12 x 12
%
%   Example
%       D = dynamic_stiff(w,L,E,rho,S,Izz,Iyy,Ixx,J,G,eta,ni)
%       where
%       w = angular frequency in rad/s
%       L = beam length
%       E = Young's Modulus
%       rho = density
%       S = cross sectional area
%       Izz = second moment of area over axis zz
%       Iyy = second moment of area over axis yy
%       Ixx = polar second moment of area (Ixx = Iyy + Izz)
%       J = torsional constant (equal to Ixx for beams with circular cross
%       sections
%       G = shear modulus
%       eta = loss factor
%       ni = Poisson's ratio
%
%   version 1.0. Developed by Paulo J. Paupitz Goncalves.

if length(varargin) < 12
    error('Not enough input arguments')
elseif length(varargin) > 13    
    error('Too many input arguments')
else
    w   = varargin{1}; 
    L   = varargin{2}; 
    E   = varargin{3}; 
    rho = varargin{4}; 
    S   = varargin{5}; 
    Izz = varargin{6}; 
    Iyy = varargin{7}; 
    Ixx = varargin{8}; 
    J   = varargin{9}; 
    G   = varargin{10}; 
    eta = varargin{11}; 
    ni  = varargin{12}; 
    
    D = zeros(12)     % Dynamic stiffness matrix
    E = E*(1+j*eta);  % Complex Young's Modulus
    G = G*(1+j*eta);  % Complex Shear Modulus

    gamay  = sqrt(w) * ((rho*S)/(E*Iyy))^(1/4) * L;
    gamaz  = sqrt(w) * ((rho*S)/(E*Izz))^(1/4) * L;
    psi    = w * L * sqrt(rho/E);
    lambda = w * L * sqrt(rho/G);

    deltay = cosh(gamay) * cos(gamay) - 1;
    F1y    = -gamay   * (sinh(gamay) - sin(gamay)) / deltay;
    F2y    = -gamay   * (cosh(gamay) * sin(gamay) - sinh(gamay) * cos(gamay)) / deltay;
    F3y    = -gamay^2 * (cosh(gamay) - cos(gamay)) / deltay;
    F4y    =  gamay^2 * (sinh(gamay) * sin(gamay)) / deltay;
    F5y    =  gamay^3 * (sinh(gamay) + sin(gamay)) / deltay;
    F6y    = -gamay^3 * (cosh(gamay) * sin(gamay) + sinh(gamay)*cos(gamay))/deltay;

    deltaz = cosh(gamaz) * cos(gamaz) - 1;
    F1z    = -gamaz   * (sinh(gamaz) - sin(gamaz)) / deltaz;
    F2z    = -gamaz   * (cosh(gamaz) * sin(gamaz) - sinh(gamaz) * cos(gamaz)) / deltaz;
    F3z    = -gamaz^2 * (cosh(gamaz) - cos(gamaz)) / deltaz;
    F4z    =  gamaz^2 * (sinh(gamaz) * sin(gamaz)) / deltaz;
    F5z    =  gamaz^3 * (sinh(gamaz) + sin(gamaz)) / deltaz;
    F6z    = -gamaz^3 * (cosh(gamaz) * sin(gamaz) + sinh(gamaz)*cos(gamaz))/deltaz;
    
    D(1,1)  =  (E*S)/L * psi * cot(psi);   D(7,7) =  D(1,1);
    D(1,7)  =  -(E*S)/L * psi * csc(psi);    D(7,1) =  D(1,7);
    D(2,2)   =  (E*Izz)/L^3 * F6z;
    D(2,6)   = -(E*Izz)/L^2 * F4z;     D(6,2)  = D(2,6);
    D(2,8)   =  (E*Izz)/L^3 * F5z;     D(8,2)  = D(2,8);
    D(2,12)  =  (E*Izz)/L^2 * F3z;     D(12,2) = D(2,12); 
    
    D(3,3)   =  (E*Iyy)/L^3 * F6y;
    D(3,5)   =  (E*Iyy)/L^2 * F4y;     D(5,3)  = D(3,5);
    D(3,9)   =  (E*Iyy)/L^3 * F5y;     D(9,3)  = D(3,9);
    D(3,11)  = -(E*Iyy)/L^2 * F3y;     D(11,3) = D(3,11);
    
    D(4,4)   =  (G*J)/L * lambda * cot(lambda);   D(10,10) =   D(4,4);
    D(4,10)  = -(G*J)/L * lambda * csc(lambda);     D(10,4) = D(4,10);
    
    D(5,5)   =  (E*Iyy)/L   * F2y;
    D(5,9)   =  (E*Iyy)/L^2 * F3y;     D(9,5)  = D(5,9);
    D(5,11)  =  (E*Iyy)/L   * F1y;     D(11,5) = D(5,11);
    
    D(6,6)   =  (E*Iyy)/L   * F2z;
    D(6,8)   = -(E*Iyy)/L^2 * F3z;     D(8,6)  = D(6,8);
    D(6,12)  =  (E*Iyy)/L   * F1z;     D(12,6) = D(6,12);
    
    
    D(8,8)   =  (E*Izz)/L^3 * F6z;
    D(8,12)  =  (E*Izz)/L^2 * F4z;     D(12,8) = D(8,12);
    
    D(9,9)   =  (E*Izz)/L^3 * F6y;
    D(9,11)  = -(E*Izz)/L^2 * F4y;     D(11,9) = D(9,11);
    
    D(11,11) =  (E*Iyy)/L   * F2y;    
    
    D(12,12) =  (E*Izz)/L   * F2z;    


end