1 function [perm,var_lnk_actual]=random_perm(var_lnk,corr_lenx,corr_leny,Nx,Ny,Lx,Ly)
   2 
   3 % Wave numbers
   4 kx = (2*pi/Lx)*[0:(Nx/2-1) (-Nx/2):(-1)]';
   5 ky = (2*pi/Ly)*[0:(Ny/2-1) (-Ny/2):(-1)]';
   6 [KX,KY]= meshgrid(kx,ky);
   7 KX2= KX.^2;
   8 KY2= KY.^2;
   9 dkx=sqrt(KX2(1,2)); dky=sqrt(KY2(2,1));
  10 
  11 
  12 %%%%Spectral density functions%%%%
  13 
  14 %Whittle-A spectrum, isotropic medium
  15 %a=pi/(4*sqrt(corr_lenx*corr_leny));
  16 %S=(2*var_lnk*a*a/pi)*( (KX2+KY2) ./ ((a*a+KX2+KY2).^3) );  
  17 
  18 %Gaussian spectrum, Random Field Generator, Ruan and McLaughlin, Adv Water Res. 1998,
  19 %S=(1/2/pi)*var_lnk*corr_lenx*corr_leny*exp(-0.5*(corr_lenx^2)*(KX2) - 0.5*(corr_leny^2)*(KY2));
  20 
  21 %Spectrum of modified exp autocovariance, Gelhar & Axness, WRR 1983
  22 S=(2/pi)*var_lnk*corr_lenx*corr_leny;
  23 S=S./( (1+(corr_lenx^2).*KX2+(corr_leny^2).*KY2).^3 ); 
  24 
  25 
  26 H=S.^(0.5); 
  27 theta=2*pi*rand(Ny,Nx); %random phase angle
  28 
  29 dZ=H.*exp(1i*theta).*sqrt(dkx*dky)*Nx*Ny;
  30 lnK=ifft2(dZ); %log(perm)
  31 K1=sqrt(2)*real(lnK); %K2=sqrt(2)*imag(lnK); 
  32 var_lnk_actual=var(reshape(K1,Ny*Nx,1));
  33 
  34 perm=exp(K1);