% this code is going to simulate the solid-on-solid model.

function H=SOS(p_plus,p_minus,slope,T,L)

%% first define a matrix that will hold all of the height information. 
% the ith row corresponds to the matrix at time i
% the jth column corresponds to the height a distance j from the starting
% point

H=ones(T,L);

%% set the initial condition
H(1,2:2:end)=0;% set every other point to be -1 so it starts by looking like a saw blade

% now add in a slope
H(1,:)=H(1,:)+(1:L)*slope;

%% now run the growth equation in a for loop

for i=2:T
    H(i,:)=growth_equ_even(p_plus,p_minus,H(i-1,:));
    H(i,:)=growth_equ_odd(p_plus,p_minus,H(i,:));
end

%% make a plot
figure
imagesc(H)
axis xy
xlabel('position','FontSize',20)
ylabel('time','FontSize',20)


end% end the function

%% Now define the growth equation
function  H_=growth_equ_even(p_plus,p_minus,H_)

%% first find all of the maxima and minima
% Find the maxima and the minima using the second derivative test
dd_H=diff(diff(H_));

% put the entire thing on a torus
dd_H=[H_(2)-2*H_(1)+H_(end),dd_H,H_(1)-2*H_(end)+H_(end-1)];
%dd_H=[0,dd_H,0];
%% now decide which heights get changed

random_mat=rand(size(H_));

%% minima can get a block added to them
% find the minima for which the corresponding element in random_mat<p_minus
add_block=dd_H>0 & random_mat<p_plus; % this vector is 1 where ever the condition (dd_H>0 & random_mat<p_plus) is met and 0 elsewhere

H_(2:2:end)=H_(2:2:end)+2*add_block(2:2:end);

%% maxima can lose a block
% find the minima for which the corresponding element in random_mat<p_minus
rm_block=dd_H<0 & random_mat<p_minus; % this vector is 1 where ever the condition (dd_H<0 & random_mat<p_minus) is met and 0 elsewhere

H_(2:2:end)=H_(2:2:end)-2*rm_block(2:2:end);
end

%% Now define the growth equation
function  H_=growth_equ_odd(p_plus,p_minus,H_)

%% first find all of the maxima and minima
% Find the maxima and the minima using the second derivative test
dd_H=diff(diff(H_));

% put the entire thing on a torus
dd_H=[H_(2)-2*H_(1)+H_(end),dd_H,H_(1)-2*H_(end)+H_(end-1)];
%dd_H=[0,dd_H,0];
%% now decide which heights get changed

random_mat=rand(size(H_));

%% minima can get a block added to them
% find the minima for which the corresponding element in random_mat<p_minus
add_block=dd_H>0 & random_mat<p_plus; % this vector is 1 where ever the condition (dd_H>0 & random_mat<p_plus) is met and 0 elsewhere

H_(1:2:end)=H_(1:2:end)+2*add_block(1:2:end);

%% maxima can lose a block
% find the minima for which the corresponding element in random_mat<p_minus
rm_block=dd_H<0 & random_mat<p_minus; % this vector is 1 where ever the condition (dd_H<0 & random_mat<p_minus) is met and 0 elsewhere

H_(1:2:end)=H_(1:2:end)-2*rm_block(1:2:end);
end


