function [ R, mr, vr, e, a0, c0 ] = doDisaggExperiment( A, g, v, m )

%arguments:

% A = input image, e.g. uniform random (must be square)
% g = Lagrange multiplier
% optional::
% v = the variance to scale to
% m = the mean to scale to

%return values:

% R  = the regularised matrix
% mr = mean of R (the regularised output)
% vr = variance of R
% e  = image entropy
% a0 = variogram range
% c0 = variogram sill

% if a0 and c0 are not specified as outputs then the variogram will not
% be calculated. This is VERY useful as the variogram can take a long
% time to generate

rescaleVariance=0 ;

% image size
s=size( A );
N=s( 1 );

% check matrix is square
if s( 1 ) ~= s( 2 )
	error 'Input matrix must be square'
end

if nargin >= 3
	rescaleVariance=1 ;
end
if nargin >= 4
	rescaleMean=1 ;
end

% apply the langrangian constraint =====================
R=reg2D( A, g );

r=reshape( R, N*N, 1 );
%a=reshape( A, N*N, 1 );

if rescaleVariance
	r = (v^0.5)*r/std(r) ;
	R=reshape( r, N, N );
end

if rescaleMean
	r = r - mean( r ) + m ;
	R=reshape( r, N, N );
end
    
mr=mean( r );
vr=var( r );

% calculate entropy ===================================
h=hist( r, N );
h=h/(N*N);
h(h==0)=eps;
h=h.*log(h);
h(isnan(h)) = 0;
e = -(sum(h)) ;

%optionally calculate the variogram ===================

if nargout == 6

	% set up x and y vectors for calculating the variogram
	y=[];
	for i=1:N
	 y=[y [1:N]];
	end
	y=y';

	x=[];
	for i=1:N
	 for j=1:N
	  x=[x i];
	 end
	end
	x=x';

	% calculate experimental variogram
	vgrm=variogram( [x y], r );
    
	% fit spherical model    
	p = fitSpherical( [10;2], vgrm );
    
	a0=p(1);
	c0=p(2);    
end  
%end