CalcLinEvolution(in_vec)

function out_vec = CalcLinEvolution(in_vec)
%  Version: 1.0, 19 March 2008
%
%   Calculates the predicted linear solution of a perturbed Kirchhoff vortex
%   See Love (1893), Guo et al. (2004), section IIC of paper

lambda_m            = in_vec(1);
mu_m_plus           = in_vec(2);
mu_m_neg            = in_vec(3);
m_1                 = in_vec(4);
sim_end_time        = in_vec(5);
frame_steps         = in_vec(6);
a                   = in_vec(7);
b                   = in_vec(8);
zeta                = in_vec(9);
n_ellipse_points    = in_vec(10);
A_m                 = in_vec(11);
B_m                 = in_vec(12);

s = linspace(0,2*pi,n_ellipse_points+1);
s = s(1:n_ellipse_points);
c = (a^2 - b^2)^0.5;
omega           = zeta * a*b/(a+b)^2;
u_0 = acosh(a/c);
J = c^2 / 2 * (cosh(2 * u_0) - cos(2 * s));         % = h_0^{-2} in Love notation           
        
for i_frame_step = 1:double(int16(sim_end_time / frame_steps - 0.01)) + 1        

    t                   = (i_frame_step-1) * frame_steps;           
    if mu_m_neg > 0                                                             % unstable
        q_vec =          (A_m * cosh(lambda_m * t) - mu_m_plus / lambda_m * B_m * sinh(lambda_m * t)) * cos(m_1 * s);
        q_vec = q_vec +  (B_m * cosh(lambda_m * t) - lambda_m / mu_m_plus * A_m * sinh(lambda_m * t)) * sin(m_1 * s);
    else                                                                        % stable
        q_vec =          (A_m * cos(lambda_m * t) - mu_m_plus / lambda_m * B_m * sin(lambda_m * t)) * cos(m_1 * s);
        q_vec = q_vec +  (B_m * cos(lambda_m * t) + lambda_m / mu_m_plus * A_m * sin(lambda_m * t)) * sin(m_1 * s);            
    end
    u_pert = q_vec ./ J; 

    if i_frame_step == 1
        T = [t];
        Y_Love_lin = [c * cosh(u_0 + u_pert) .* cos(s) + i * c * sinh(u_0 + u_pert) .* sin(s)];
    else
    T                   = [T.' t.'].';

    x_lin               = c * cosh(u_0 + u_pert) .* cos(s);                 % add rotation of ellipse at Kirchhoff frequency
    y_lin               = c * sinh(u_0 + u_pert) .* sin(s);
    ellipse_phase       = omega * t;
    x_lin_p             = x_lin * cos(ellipse_phase) - y_lin * sin(ellipse_phase);
    y_lin_p             = y_lin * cos(ellipse_phase) + x_lin * sin(ellipse_phase);           
    Y_Love_lin = [Y_Love_lin.', [x_lin_p + i * y_lin_p].'].';  
    end                        
end
Y_Love = Y_Love_lin;
out_vec = [T Y_Love];