%bp4eig.m
%Hunter Gilbert
%Last Modified 9/15/08
%
%[X,E]=bp4eig(m,mtilde,ell,elltilde,alpha)
%Solves for the eigenvalues of a beaded string damped at a point
%
%m, mtilde - the masses on the left and on the right
%ell, elltilde - the lengths on the left and right, including the
%   segments that meet at the damper
%alpha - the damping coefficient (alpha>=0)
%X - The eigenvectors
%E - A diagonal matrix of eigenvalues corresponding to the eigenvectors in
%   X
%Note: the convention for indexing the masses and lengths is in the
%style of Boyko and Pivovarchik
%
function [X,E]=bp4eig(m,mtilde,ell,elltilde,alpha)
m=reshape(m,length(m),1);
mtilde=reshape(mtilde,length(mtilde),1);
ell=reshape(ell,length(ell),1);
elltilde=reshape(elltilde,length(elltilde),1);

if (alpha~=0)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Code 'A' Matrix with damping and a degree of freedom at the damper
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    A=zeros(2*(length(m)+length(mtilde))+1); %A will have size 2*n+1 where n is the number of masses

    %%% from eq. 4.1, but in the non-rotated coordinate system
    k=1;
    for j=1:2:2*length(m)-1
        A(j,j+1)=1; %Tautology, e.g. u1'=u1'
        A(j+1,j)=-1/m(k)*(1/ell(k)+1/ell(k+1)); %the u_k term
        A(j+1,j+2)=1/ell(k+1)/m(k); %the u_k+1 term
        if (k>1)
            A(j+1,j-2)=1/ell(k)/m(k); %the u_k-1 term
        end
        k=k+1;
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


    j=2*length(m)+1;

    %%% from eq. 4.4, but in the non-rotated coordinate system
    A(j,j)=-1/alpha*(1/ell(end)+1/elltilde(end));
    A(j,j-2)=1/alpha/ell(end);
    A(j,j+1)=1/alpha/elltilde(end);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%% from eq. 4.2, but in the non-rotated coordinate system
    j=j+1;
    k=1;
    elltilde = flipud(elltilde); %we are going through the right side 'backwards' (left to right)
    mtilde = flipud(mtilde); %we are going through the right side 'backwards' (left to right)
    
    for j=j:2:2*(length(m)+length(mtilde)),
        A(j,j+1)=1; %Tautology, e.g. u~1'=u~1'
        A(j+1,j)=-1/mtilde(k)*(1/elltilde(k)+1/elltilde(k+1)); %the u~_k term
        if (k~=1)
            A(j+1,j-2)=1/elltilde(k)/mtilde(k); %the u~_k+1 term if we are not at the bead closest to the center
        else
            A(j+1,j-1)=1/elltilde(k)/mtilde(k); %the u~_k+1 term if we are at the bead closest to the center
        end
        if (k<length(mtilde))
            A(j+1,j+2)=1/elltilde(k+1)/mtilde(k); %the u~_k-1 term
        end
        k=k+1;
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
else
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Code 'A' Matrix without damping and without a degree of freedom at the
    %damper
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    A=zeros(2*(length(m)+length(mtilde))); %A will have size 2*n where n is the number of masses

    %%% from eq. 4.1, but in the non-rotated coordinate system
    k=1;
    for j=1:2:2*length(m)-1
        A(j,j+1)=1; %Tautology, e.g. u1'=u1'
        if (k<length(m))
            A(j+1,j+2)=1/ell(k+1)/m(k); %the u_k+1 term
        else
            A(j+1,j+2)=1/(elltilde(end)+ell(end))/m(k); % the u_k+1 term
        end
        if (k>1)
            A(j+1,j-2)=1/ell(k)/m(k); %the u_k-1 term
        end
        if (k~=length(m))
            A(j+1,j)=-1/m(k)*(1/ell(k)+1/ell(k+1)); %the u_k term
        else
            A(j+1,j)=-1/m(k)*(1/ell(k)+1/(ell(k+1)+elltilde(end))); %the u_k term
        end
        k=k+1;
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%% from eq. 4.2, but in the non-rotated coordinate system
    j=2*length(m)+1;
    k=1;
    
    mtilde = flipud(mtilde); % we are going through the right side 'backwards' (left to right)
    elltilde = flipud(elltilde); % we are going through the right side 'backwards' (left to right)
    
    for j=j:2:2*(length(m)+length(mtilde)),
        A(j,j+1)=1; %Tautology, e.g. u~1'=u~1'
        if (k>1)
            A(j+1,j)=-1/mtilde(k)*(1/elltilde(k)+1/elltilde(k+1)); %the u~_k term
        else
            A(j+1,j)=-1/mtilde(k)*(1/(elltilde(k)+ell(end))+1/elltilde(k+1)); % the u~_k term
        end
        if (k~=1)
            A(j+1,j-2)=1/elltilde(k)/mtilde(k); %the u~_k+1 term if we are not at the bead closest to the center
        else
            A(j+1,j-2)=1/(ell(end)+elltilde(k))/mtilde(k); %the u~_k+1 term if we are at the bead closest to the center
        end
        if (k<length(mtilde))
            A(j+1,j+2)=1/elltilde(k+1)/mtilde(k); %the u~_k-1 term
        end
        k=k+1;
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
[X,E]=eig(A);