function generatorOSSRGFdiag(lambda1,lambda2,M,s)


%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 allows to simulate different OSSRGFs satisfying the same scaling property 
%  X(a^Ex)=a^H X(x)
% for the diagonal matrix E=[lambda1  0;  0 lambda2]

% The program works as follows: 
%from the input, a canonical OSSRGF is constructed (from the canonical E-pseudo-norm) 
% The program propose to chose a function Geval, continuous and strictly positive : different examples are proposed in the code, but one can compute others.
% It has to be done inside the code (selecting the one or writing a new one)
%Selecting this function Geval allows to generate different OSSRGFs satisfying the same scaling property.



%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) 
%OSSRGF constructed from the canonical pseudo-norm

%figure(2)
%Other OSSRGF contructed via the function Geval


%COORDONNATES 

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;





%CONSTRUCTION OF THE USUAL PSEUDO-NORM AND THE USUAL SPECTRAL DENSITY
A=(abs(X)).^(2/lambda1);
B=(abs(Y)).^(2/lambda2);
%a=size(A)

N=A(ones(1,2*2^M+1),:);
Pa=B(ones(1,2*2^M+1),:);
P=Pa';


rho=(sqrt(N+P));
phi=(sqrt(N+P)).^(s+(lambda1+lambda2)/2);




%CREATION OF COORDONATES (U,V)=(rho(x,y)^(-lambda1)x, rho(x,y)^(-lambda2)y)
Xfull=(X'*ones(1,2*2^M+1))';
Yfull=Y'*ones(1,2*2^M+1);

Xdiag=Xfull.*eye(2*2^M+1);
Ydiag=Yfull.*eye(2*2^M+1);

U=(rho.^(-lambda1))*Xdiag;

V=((rho.^(-lambda2))'*Ydiag')';


%FUNCTION G EVALUATED IN (U,V) 
%G HAS TO BE STRICTLY POSITIVE AND CONTINUOUS OUTSIDE O
%SEVERAL EXAMPLES ARE PROPOSED...


%EXAMPLE1: leads to the same field (test function)
%Geval=sqrt((abs(U)).^(2/lambda1)+(abs(V)).^(2/lambda2));

%EXEMPLE2:
%Geval=1+abs(U)+abs(V);

%EXEMPLE3:
Geval=(1+abs(U)).*(1+abs(V));

%EXEMPLE4
%Geval=1+abs(U);

%EXEMPLE5
%Geval=(1+abs(U)).*exp(V);

%EXEMPLE6
%Geval6=exp(V).*exp(1+cos(U).^2);

%...

%NEW ANISOTROPIC PSEUDO_NORM TAU AND NEW SPECTRAL DENISTY PSI
tau=Geval.*rho;
psi=tau.^(s+(lambda1+lambda2)/2);


%CREATION OF THE OSSRGFs
Z=randn(2*2^M+1,2*2^M+1);

%Index 1= OSSRGF constructed from the canonical pseudo norm
%Index 2= OSSRGF constructed from the pseudo norm obtained via Geval

W1=fftshift(fft2(Z))./phi;

W2=fftshift(fft2(Z))./psi;



I1=real(ifft2(ifftshift(W1)));
I2=real(ifft2(ifftshift(W2)));

J1=I1-I1(2^M+1,2^M+1);
J2=I2-I2(2^M+1,2^M+1);



mini1=min(min(J1));
Maxi1=max(max(J1));
mini2=min(min(J2));
Maxi2=max(max(J2));



figure(1)
imagesc(J1,[mini1 Maxi1]);
colormap(gray)

figure(2)
imagesc(J2,[mini2 Maxi2]);
colormap(gray)