function [vp,vpfull] = powconv3(v,p)
% [vp,vpfull] = powconv3(v,p)
% computes the convolution power v^p for Fourier series in 3D 
% 
% The input v is in tensor form and should have odd length in all dimensions
% The output vp has the same size as v
% The output vpfull is the bigger tensor with all nonvanishing elements
%
% The function works both for floating point numbers and for intervals

cs = size(v);
N = (cs(1)-1)/2;
M = 2^nextpow2(2*p*N+1);

w = ILzeros([M,M,M],v(1));
s1 = M/2 + 1 - N;
s2 = M/2 + 1 + N;
w(s1:s2,s1:s2,s1:s2) = v;
ww = w([M/2+1:M,1:M/2],[M/2+1:M,1:M/2],[M/2+1:M,1:M/2]); 

u = ILfftn(ww);
up = reshape(u(:).^p,size(u));
v2 = ILifftn(up);

v3 = v2([M/2+1:M,1:M/2],[M/2+1:M,1:M/2],[M/2+1:M,1:M/2]);
vp = v3(s1:s2,s1:s2,s1:s2);

s1p = M/2 + 1 - p*N;
s2p = M/2 + 1 + p*N;
vpfull = v3(s1p:s2p,s1p:s2p,s1p:s2p);

return



