ÿþfunction OSSRGF_diagonal(M,l1,l2,H)

%Authors : M. Clausel and B. Vedel, 2008
%clausel@univ-paris12.fr
%beatrice.vedel@ens-lyon.fr

%M. Clausel: Laboratoire d'Analyse Mathématique et Appliquée 
% UMR 8050 du CNRS 
% Université Paris Est, 
% 61 Avenue du Général de Gaulle
%94100 Créteil Cedex, France

%B. Vedel: Laboratoire de Physique
%CNRS UMR 5672
%ENS de Lyon
%46, allée d'Italie 
%69007 Lyon



%This program simulates a OSSRGF X satisfying the scaling property
% X(a^E x)= a^H X(x)
%for the diagonal matrix E=[l1,0;0,l2]


%INPUT
%M=scale of the cartesian grid (2^(M+1) x 2^(M+1) points)
%l1 and l2 : eigen values of the diagonal matrix (coeff of the matrix). They have to be strictly positive

%H=Hurst coeff. chosen between 0 and min(l1,l2) (strictly)

%OUTPUT: 
%figure(1)
%1/phi where phi is the spectral density associated to the canonical pseudo-norm 

%figure(2)
%simulation of the Fourier transform of the OSSRGF of spectral density phi

%figure(3)
%simulation of the OSSRGF


%COORDONATES


X=(-2*2^M:2:2*2^M)/(2^(M+1));
X(2^M+1)=1/2^M;
Y=(-2*2^M:2:2*2^M)/(2^(M+1));
Y(2^M+1)=1/2^M;
XX=X(ones(1,2*2^M+1),:);
YY=Y(ones(1,2*2^M+1),:)';
clear X Y


%CONSTRUCTION OF THE PSEUDO NORM AND THE SPECTRAL DENSITY

%rho is the classical pseudonorm associated to the diagonal matrix with eigenvalues  l1 et l2: 
%rho(x,y)=(abs(x)^(2/l1) + abs(y)^(2/l2) )^(1/2)

rho=sqrt(abs(XX).^(2/l1)+abs(YY).^(2/l2));
clear XX YY

%phi is the spectral density of the field obtained via rho

phi=rho.^(H+(l1+l2)/2);
clear rho


%CONSTRUCTION OF THE FOURIER TRANSFORM OF THE FIELD (WITHOUT RENORMALIZATION)
%W= Fourier transform of the OSSRGF 
Z=randn(2*2^M+1,2*2^M+1);
W=fftshift(fft2(Z))./phi;
clear Z

%figure(1)
%mesh(X,Y,phi)

%figure(2)
%mesh(X,Y,W)

%CONSTRUCTION OF THE FIELD (FOURIER INVERSE + RENORMALIZATION)
T=real(ifft2(ifftshift(W)));
Zp=T-T(2^M+1,2^M+1);
mini=min(min(Zp));
Maxi=max(max(Zp));


figure(3)
imagesc(Zp,[mini Maxi]);
%colormap(gray);