% Markov-chain simulation of soil textural profiles
% Budiman September 2000
% Based on:
% Li, Li, Shi (1999) Markov-chain simulation of soil textural profiles
% Geoderma 92,37-53

clear all;
% s = sand; sl=sandy loam; ll=light loam; ml =medium loam; cl=clay loam; c=clay
% A = IPD Initial Probability Density
%   a      b     c      d    e     f
A=[0.00  0.779 0.100 0.042 0.000 0.079];

% P = TPM Transition Probability Matrix (5 cm layers)
P=[
0.935	0.003	0.006	0.014	0.001	0.040
0.022	0.904	0.020	0.021	0.004	0.029
0.040	0.031	0.828	0.028	0.012	0.061
0.060	0.022	0.007	0.788	0.004	0.118
0.040	0.010	0.010	0.020	0.840	0.080
0.075	0.020	0.017	0.025	0.006	0.857];

% Initialize nlayer, nprofile
m=6;              % no. texture class
nlayer=40;        % no. layers in a profile
layer_thick=5;    % thickness of each layer
nprofile=100;      % no. profiles to simulate

% initialize matrix to store the results
texture=ones(nprofile,nlayer);
thick=zeros(nprofile,m);
depth=zeros(nprofile,m);

% randomize
rand('state',sum(100*clock));

% calculate cumulative probability CA & CP
CA=cumsum(A')';
CP=cumsum(P')';


for profile=1:nprofile,
   
   % simulate 1st layer (upper boundary condition)
	layer=1;
	% generate random number & compare with the probability matrix
    tex=min(find(CA>rand));
	texture(profile,layer)=tex;
	thick(profile,tex)=thick(profile,tex)+1;

	% Generate 2nd to n layers 

	for layer=2:nlayer,
        tex=min(find(CP(tex,:)>rand));
   	    texture(profile,layer)=tex;
        thick(profile,tex)=thick(profile,tex)+1;
        depth(profile,tex)=layer;      
	end;
   
end;

% display the matrix
pcolor(texture');

% calculate thickness & depth for each textural layer
depth=(depth-thick)*layer_thick;
thick=thick*layer_thick;
thick_mean=mean(thick)
thick_std=std(thick)
depth_mean=mean(depth)
depth_std=std(depth)

% plot histogram of depth
figure(1);
for tex=1: m
   subplot(2,3,tex); 
   hist(depth(:,tex),20);
end;


% calculate the simulated TPM
TPM=zeros(m,m);
for profile=1: nprofile,
   for layer=1:nlayer-1,
      i=texture(profile,layer);
      j=texture(profile,layer+1);
      TPM(i,j)=TPM(i,j)+1;
   end;  
end;
tsum=sum(TPM);
TPM=TPM./repmat(tsum,m,1)