% create known analytic array and slopes
% create random array by smoothing a random array with a power spectrum -1 
% calculate slopes from random array, and then add 50%  noise to the slopes

Qx = 3.7; Qy = 2.4; 							% periods of analytic array
phi = cos (Qy* (0.5:Ny-0.5)*pi/Ny)'*cos (Qx* (0.5:Nx-0.5)*pi/Nx) ;	% create array
phi= 1e2 * (phi + 0.2 * oneY' * X1 / Nx - 0.4 * Y1' * oneX  / Ny);	% adds small linear drift
Px = S; Py = S;   P = phi;    						% slopes
Px (1:Ny, 1:Nx - 1) = phi (1:Ny, 2:Nx) - phi (1:Ny, 1:Nx - 1);
Py (1:Ny - 1, 1:Nx) = phi (2:Ny, 1:Nx) - phi (1:Ny - 1, 1:Nx);

% random function, smoothed
l = 2 * max (Nx, Ny); lh = l/2+1; x = -l/2:l/2-1;  xh = 1:l; 		%l/4:3*l/4; % sizes
[XX, YY]  = meshgrid (x, x); R2 = XX.^2+YY.^2; R2 (lh, lh) = 1; 	% create r^2 array
XX  = real (ifft2 (fft2 ((rand(l) - 0.5)) .* fftshift (R2.^-1)));  	% smoothed random function 
XX = 1e4 * (XX  + .03*ones(l,1)*x/l+ .02*x'*ones(1,l)/l); % add small slope
R2 = R2<950;      					% create circle
R2 = 1;                 				% comment this line for circle case
YY = XX .* R2;     					% circle or full array
if R2 (1) ~= 1, R2 = R2 (Y1-Ny/2+l/2, X1-Nx/2+l/2); end % centre arrays
phi = YY (Y1-Ny/2+l/2, X1-Nx/2+l/2) .* R2;             	% chop out relevant central part

Ax = S; Ay = S;  A = phi; 				% slopes
Ax (1:Ny, 1:Nx-1) = phi (1:Ny, 2:Nx) - phi (1:Ny, 1:Nx-1);
Ay (1:Ny-1, 1:Nx) = phi (2:Ny, 1:Nx) - phi (1:Ny-1, 1:Nx);

% add noise to slopes
avg_abs_Sx= sqrt (trace(Ax*Ax')/ Nx );  avg_abs_Sy = sqrt (trace(Ay*Ay')/ Nx );
Noise = 50;    						 %noise strength, percent
Bx = Ax + (avg_abs_Sx/8*sqrt(Noise)/100)*(poissrnd(Noise, Ny, Nx) - Noise);      % add noise
By = Ay + (avg_abs_Sy/8*sqrt(Noise)/100)*(poissrnd(Noise, Ny, Nx) - Noise);      % Strength/100 noise