%% derivfrac3D.m
%
% E. Koenig, 08/04/2010
%
% Pseudo-fractional derivation of a 3D field. 
%
% This operation is seen as the inverse of the pseudo-fractional
% integration. 
% It preserves the scale invariance of the initial 3D field.
% It is realized by a high-pass filter in k^H in the Fourier domain.
%
% The initial field is scale invariant and characterized by
%
% $$ \zeta(q) = \tau(q) + qH $$
%
% The pseudo-fractional derivation will create a rougher, more intermittent
% version of the initial field. This new field is scale invariant and 
% characterized by
%
% $$ \tau(q) $$
%
%
%% Use
%
% y = derivfrac3D(x, H, S, zero);
%
%
%% Inputs
%
% * *x*:  input 3D field
% * *H*: parameter of the high-pass filter
% * *S*: size of the initial filed
% * *zero*: 1 - management of k=0, 0 - else (default)
%
%
%% Outputs
%
% * *y*: derived 3D field.
%
%



function [y]= derivfrac3D(x, H, S, zero)

lx=size(x);

if nargin<4
    zero=0;    
    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 zero
if zero
    listnull=find(Kmod==0);
    Kmod(listnull)=min(Kmod(setdiff(1:length(Kmod(:)),listnull)))*ones(1,length(listnull));
    clear listnull
end

clear zero

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



