%% integfrac3D.m
%
% E. Koenig, 08/04/2010
%
% Pseudo-fractional integration of a 3D scale invariant field. 
%
% This operation preserves the scale invariance of the initial 3D field. 
% It is realized by a low-pass filter in 1/k^H in the Fourier domain.
% The value at k=0 is set to the minimum k value which is not null.
%
% The initial field is scale invariant and characterized by 
%
% $$\tau(q)$$
%
% and is very rough, intermittent.
%
% The pseudo-fractional integration will create a smoother version of the
% initial field. This new field is still scale invariant and characterized 
% by
%
% $$ \zeta(q) = \tau(q) + qH $$
%
%
%% Use
%
% y = integfrac3D(x, H, S);
%
%
%% Inputs
%
% * *x*: input 3D field
% * *H*: parameter of the low-pass filter.
% * *S*: size of the initial filed.
%
%
%% Outputs
%
% * *y*: integrated 3D field.
%
%



function [y]= integfrac3D(x, H, S)

lx=size(x);

if nargin<4   
    if nargin<3
        S = 2.^floor(log2(lx));
        if S ~= lx
            disp('x should be of size 2^N')
        end
    end
end

clear lx

x=x(1:S(1),1:S(2),1:S(3));

% Fourier Transform
X=fftshift(fftn(x));
clear x
kx=-S(1)/2:S(1)/2-1;
ky=-S(2)/2:S(2)/2-1;
kz=-S(3)/2:S(3)/2-1;
[Kx,Ky,Kz]=meshgrid(kx,ky,kz);
clear kx ky kz S


Kmod = sqrt(Kx.^2+Ky.^2+Kz.^2);
clear Kx Ky Kz

% Management of the value on (0,0,0)
listnull=find(Kmod==0);
Kmod(listnull)=min(Kmod(setdiff(1:length(Kmod(:)),listnull)))*ones(1,length(listnull));
clear listnull

% Low-pass filtering
Hfilt=1./(Kmod.^(H));
clear H Kmod
Y=Hfilt.*X;
clear X Hfilt
y=ifftn(ifftshift(Y));
clear Y



