function OSSRGF_spirale(M,alpha,beta,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^Ex)=a^H X(x)
% for the spirale matrix [alpha, beta; - beta, alpha]

%ENTRIES
%M=scale of the cartesian grid (2^(M+1) x 2^(M+1) points)
%alhpa beta : coefficients of the spirale matrix E (matrix of anisotropy or enxponent), 
% alpha has to be stricly positive, beta real

%H=Hurst coeff, has to be chosen between 0 and alpha (srictly)


%OUTPUT: 
%figure(1)
%Simulation of the field



%COORDINATES

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

NormCarre=XX.^2+YY.^2;

%CONSTRUCTION OF THE SPECTRAL DENSITY PHI OF THE OSSRGF

phi=(abs(XX.*cos(beta/(2*alpha)*log(NormCarre))-YY.*sin(beta/(2*alpha)*log(NormCarre)))./(NormCarre)).^((H+2*alpha)/alpha);
%norm(phi)

clear XX YY NormCarre

%CONSTRUCTION OF THE FOURIER TRANSFORM W OF THE OSSRGF (WITHOUT RENORMALIZATION) 

%W= la transformee de Fourier du champs gaussien d'anisotropie E
Z=randn(2*2^M+1,2*2^M+1);
W=fftshift(fft2(Z)).*phi;
clear Z

%CONSTRUCTION OF THE FIELD 
T=real(ifft2(ifftshift(W)));
%clear W

Zp=T-T(2^M+1,2^M+1);
mini=min(min(Zp));
Maxi=max(max(Zp));

figure(1)
imagesc(Zp,[mini,Maxi])
colormap(gray)