function [ymu] = powermethod_LW(npixels,nangles)
%------------------------------------------------------------------
% Function: powermethod_LW(npixels,nangles)
%
% Inputs: npixels - Number of image/projection pixels
%         nangles - Number of projection angles
%
% Output:     ymu - Maximum eigenvalue for the Landweber algorithm
%------------------------------------------------------------------

x = ones(npixels,npixels);
xvec = lex_order(x,npixels,npixels);

maxiter = 100;
OK = true;
TOL = 10^-8;

% Step 1
k = 1;

% Step 2
lp = 1;
amax = abs(xvec(1));
for i = 2:npixels.^2;
    if abs(xvec(i)) > amax
        amax =abs(xvec(i));
        lp = i;
    end
end

% Step 3
xvec = xvec./amax;

% Step 4
while (k <= maxiter & OK == true)
    
    y = zeros(npixels,npixels);
    
    % Step 5
    x = unlex_order(xvec,npixels,npixels);
    temp = project2D(x,nangles);
    y = backproject2D(temp);
    
    % Step 6
    yvec = lex_order(y,npixels,npixels);
    ymu = yvec(lp);
    
    % Step 7
    lp = 1;
    amax = abs(yvec(1));
    for i = 2:npixels.^2;
        if abs(yvec(i)) > amax
          amax = abs(yvec(i));
          lp = i;
        end
    end
    
    % Step 8
    if amax <= 0 
        fprintf(1,'0 eigenvalue - select another ');
        fprintf(1,'initial vector and begin again\n');
        OK = false;
    else
        
        % Step 9
        err = 0;
        for i = 1:npixels.^2
            t = yvec(i)/y(lp);
            if abs(xvec(i)-t) > err 
                err = abs(xvec(i)-t);
            end
            xvec(i) = t;
        end
    
        % Step 10
        if err <= TOL 
            OK = false;
        end
    
        % Step 11
        k = k + 1;
    
    end
    
end

% Step 12
if k > maxiter
    fprintf(1,'Method did not converge within %d iterations\n',maxiter);
end

ymu = roundoff(2./ymu-0.005,3);
