function [soliton,acc_check] = Petviashvili_mod(beta,N,L,uinit)

% Petviashvili_mod.m - Computes the soliton by the modified Petviashvili iteration
% (Version with nonradial corrections)
%
% By Laurent Demanet, 2005

p = 2*beta + 1;
gamma = p/(p-1);
delta = -1/2;

h = L/N;
if mod(N,2),
    absc = -(L/2-h/2):h:(L/2-h/2);
else
    absc = -L/2:h:(L/2-h);
end
[X,Y,Z] = meshgrid(absc,absc,absc);
R = sqrt(X.^2 + Y.^2 + Z.^2);

xi = 1:N;
xifreq0 = floor(N/2+1);
xi = xi - xifreq0;
[XI1,XI2,XI3] = meshgrid(xi);
absXIsquared = XI1.^2 + XI2.^2 + XI3.^2;
Petv_mult = 1+4*pi^2*absXIsquared/L^2;
xi(1) = 0;  % this is only important for the computation of the gradient
[XI1,XI2,XI3] = meshgrid(xi);
grad1 = 2*pi*sqrt(-1)*XI1/L;
grad2 = 2*pi*sqrt(-1)*XI2/L;
grad3 = 2*pi*sqrt(-1)*XI3/L;

if nargin < 4,
    u = 4*exp(-R.^2);
else
    u = uinit;
end

prev_u = u;

Petv_consts = [];
err = 1;
k = 0;

disp('k | Petv_const | Petv_mod | Euler-Lagrange remainder');
while (err > 1e-11) & (k < 50)

    k = k+1;
    prevprev_u = prev_u;
    prev_u = u;
    up = abs(u).^p;
    U = fftshift(fftn(ifftshift(u)));
    UP = fftshift(fftn(ifftshift(up)));
    Petv_const = abs(sum(sum(sum(Petv_mult.*(U.*U))))/sum(sum(sum(U.*UP))));
    Petv_mod1 = real(sum(sum(sum(Petv_mult.*(U.*(grad1.*U)))))/sum(sum(sum(grad1.*(grad1.*(U.*UP))))));
    Petv_mod2 = real(sum(sum(sum(Petv_mult.*(U.*(grad2.*U)))))/sum(sum(sum(grad2.*(grad2.*(U.*UP))))));
    Petv_mod3 = real(sum(sum(sum(Petv_mult.*(U.*(grad3.*U)))))/sum(sum(sum(grad3.*(grad3.*(U.*UP))))));
    Petv_mod = abs(Petv_mod1) + abs(Petv_mod2) + abs(Petv_mod3);
    Ucorr = Petv_mod1*(grad1.*U) + Petv_mod2*(grad2.*U) + Petv_mod3*(grad3.*U);
    Unew = UP./Petv_mult * Petv_const^gamma + delta*Ucorr;
    u = abs(fftshift(ifftn(ifftshift(Unew))));
    
    % Soliton check
%    U = fftshift(fftn(ifftshift(u)));
%    lapl_U = - 4*pi^2*absXIsquared/L^2 .* U;
%    lapl_u = real(fftshift(ifftn(ifftshift(lapl_U))));
%    eulerlagrange = lapl_u - u + u.^p;
%    EL_test1 = sqrt(sum(sum(sum(abs(eulerlagrange).^2)))*h^3) / sqrt(sum(sum(sum(abs(u).^2)))*h^3);
    
    % Componentwise Aitken's acceleration + soliton check
    denom = (u - 2*prev_u + prevprev_u);
    if ~all(all(all(denom))),
        u1 = u;
        k = 50;
    else
        u1 = u - (u - prev_u).^2 ./ denom;
    end;
    U1 = fftshift(fftn(ifftshift(u1)));
    lapl_U1 = - 4*pi^2*absXIsquared/L^2 .* U1;
    lapl_u1 = real(fftshift(ifftn(ifftshift(lapl_U1))));
    eulerlagrange1 = lapl_u1 - u1 + (u1).^p;
    EL_test2 = sqrt(sum(sum(sum(abs(eulerlagrange1).^2)))*h^3) / sqrt(sum(sum(sum(abs(u1).^2)))*h^3);

    disp([num2str(k),' | ', num2str(Petv_const),' | ', num2str(Petv_mod),' | ', num2str(EL_test2)]);
    err = EL_test2;
    
end

acc_check = (err < 1e-11);
soliton = u1;
