function generatorOSSRGFtransvection(l,M,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 allows to simulate different OSSRGFs satisfying the same scaling property 
%  X(a^Ex)=a^H X(x)
% for the transvection matrix E=[1, 0; l, 1]

% The program works as follows: 
%from the input, a 'canonical' OSSRGF is constructed (from a 'canonical' E-pseudo-norm constructed in the paper 'Explicit construction of operator scaling stable random gaussian fields' (M. Clausel and B. Vedel)) 
% 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)
%l : coeff of the matrix.
%H=Hurst coeff. chosen between 0 and 1 (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));
Y=(-2*2^M:2:2*2^M)/(2^(M+1));
X(2^M+1)=1/2^M;
Y(2^M+1)=1/2^M;




%CONSTRUCTION OF CANONICAL PSEUDO-NORM AND SPECTRAL DENSITY
A=abs(X);
B=X.*(log(abs(X)))/l;



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



%Pseudo Norm: abs(x)+abs(y-(xln(abs(x))/l)^(1/l)
rho=(L+abs(P-N)).^(1/l);
 
phi=rho.^(H+l);




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

Yfull=Y'*ones(1,2*2^M+1);


U=(rho.^(-l)).*Xfull;
V=(rho.^(-l)).*(Yfull-log(rho).*Xfull);


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



%EXEMPLE1: leads to the same field (test function)
%Geval=(abs(U)+abs(V-(U/l).*log(abs(U)))).^(1/l);

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

%EXEMPLE3:(marrant pour le champs)
Geval=(1+abs(U)).*(1+abs(V));

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

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

%...


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


%NEW ANISOTROPIC PSEUDO_NORM TAU AND NEW SPECTRAL DENISTY PSI

tau=Geval.*rho;
psi=tau.^(H+l);



%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)