function [wnrpm,V,Amat,list] = MDOFpulley(J,r,x,y,mass,EA,kT,Ldiv,Lom,ploton)
%list = {[lo1,lo2],[EA,kT],J,r,x,y} this gives all the important parameters
%in a large list

    %solve for Lengths
    L = [0,0,0,0,0,0,0];

    n = 1;
    while n <= 7
        if n==7;
            L(n) = ((x(1)-x(n))^2 + (y(1)-y(n))^2 - (r(1)-r(n))^2)^(1/2);
            n = 8;
        else
            L(n) = ((x(n+1)-x(n))^2 + (y(n+1)-y(n))^2 - (r(n+1)-r(n))^2)^(1/2);
            n = n+1;
        end
    end

    k = EA./L;

    %create stiffness matrix
    keq1 = [(k(1)+k(2))*r(2)^2, -k(2)*r(2)*r(3), 0,0,0,0,0];
    keq2 = [-k(2)*r(2)*r(3), (k(2)+k(3))*r(3)^2, -k(3)*r(3)*(-r(4)),0,0,0,0];
    keq3 = [0, -k(3)*r(3)*(-r(4)), (k(3)+k(4))*(-r(4))^2, -k(4)*(-r(4))*r(5),0,0,0];
    keq4 = [0,0,-k(4)*(-r(4))*r(5), (k(4)+k(5))*r(5)^2, -k(5)*r(5)*r(6),0,0];
    keq5 = [0,0,0, -k(5)*r(5)*r(6), (k(5)+k(6))*r(6)^2, -k(6)*r(6)*(-r(7)),(k(6)*r(6))*(x(6)-x(7))/L(6)];
    keq6 = [0,0,0,0,-k(6)*r(6)*(-r(7)),(k(6)+k(7))*(-r(7))^2 ,-(k(6)*(-r(7))*(x(6)-x(7))/L(6))-(k(7)*(-r(7))*(x(7)-x(1))/L(7))];
    keq7 = [0,0,0,0, k(6)*(x(6)-x(7))*r(6)/L(6), -(k(7)*(x(7)-x(1))*(-r(7))/L(7)) - (k(6)*(x(6)-x(7))*(-r(7))/L(6)), (k(7)*(x(7)-x(1))*(x(7)-x(1))/(L(7)^2)) + ((k(6)*(x(6)-x(7))^2)/(L(6)^2)) + kT];

    keq = [keq1;keq2;keq3;keq4;keq5;keq6;keq7];

    %adjust intertias, J = (mr^2)/2
    rorig = [81.25, 64.5, 70.60, -41.15, 30, 67.5, -38.1]/1000;
    for i = 1:6
        J(i) = (J(i)/((rorig(i)^4)))*(r(i)^4);
    end
    
    %create mass matrix   
    meq = [];
    for i = 1:6
        meq(i,i) = J(i+1);
    end
    meq(7,7) = mass;

    ceq = 0;

    [Q,lambda] = polyeig(keq,ceq,meq);

    %wn = imaginary component
    lambda2 = imag(lambda);

    %extract only positive lambdas
    lambda3 = [];
    n = 1;
    %zeros = 0
    for i = 1:length(lambda2)
        if lambda2(i) > 0
            lambda3(n) = lambda2(i);
            n = n+1;

        end
    end
    %natural freq in RPM
    wnrpm = lambda3./(2*pi)*60*(1/3);

    %extract mode shapes from eigenvectors
    Q2 = [];
    n = 1;
    for i = 1:(length(Q)/2)
        Q2(:,n) = Q(:,2*i-1);
        n = n+1;
    end
    Q3 = abs(Q2);
    V = Q3;

    for i = 1:length(Q3)
        maxQ3 = max(Q3(:,i));
        Q3(:,i) = Q3(:,i)./maxQ3;
    end
    %mode shapes for each pulley at each frequency
    V = Q3;

    %MDOF frequency response
    forces = [1;0;0;0;0;1;1];
    Amat = zeros(7,Lom/Ldiv+1);
    omegavector = zeros(1,Lom/Ldiv+1);

    n = 0;
    for omega = 0:Ldiv:Lom
        A = ((-(omega/60*(2*pi)*3).^2.*(meq) + (sqrt(-1))*(omega/60*(2*pi)*3).*(ceq) + (keq))^-1)*(forces);
        n = n+1;
        omegavector(n) = omega;
        A = abs(A);
        Amat(:,n) = A;
    end

    %full graph
    figure()
    clf
    hold on
    axis([0, Lom, 0, .025])
    [lo1,index1] = min(wnrpm);
    wnnew = [wnrpm(1:index1-1), wnrpm(index1+1:7)];
    [lo2,index2] = min(wnnew);
    if ploton >= 1
        plot(omegavector,Amat(1,:),'r','Linewidth',2)
        plot(omegavector,Amat(2,:),'g','Linewidth',2)
        plot(omegavector,Amat(3,:),'b','Linewidth',2)
        plot(omegavector,Amat(4,:),'c','Linewidth',2)
        plot(omegavector,Amat(5,:),'y','Linewidth',2)
        plot(omegavector,Amat(6,:),'k','Linewidth',2)
        plot(omegavector,Amat(7,:),'m','Linewidth',2)
%         text(1000,.014,'lo1-RPM=')
%         text(1400,.014,num2str(lo1))
%         text(1000,.012,'lo2-RPM=')
%         text(1400,.012,num2str(lo2))
    end
    
    list = {[lo1,lo2],[EA,kT],J,r,x,y};
    

    
end

