clear all
close all

% Annual deposition of salts:
% NaCl=24.105 g/m2/yr 
% CaSO4=7.014 g/m2/yr
dep_halite = 24.105;
dep_gypsum = 7.014;

% Equilibrium constant of salts:
% NaCl = 35.7g/100mL = 357 g/L = 357000 g/m3
% CaSO4 = 0.241g/100mL= 2.41 g/L = 2410 g/m3
ec_halite = 357;
ec_gypsum = 2.41;

% Enter volume of rainfall (in units cm/yr)
rainfall = 0.8;

% Soil Depth: [0:X]' where X=max depth, 0=surface layer
max_depth = 20;
depth = [0:2:max_depth]';
sample_pts = length(depth);
%%


% halite/gypsum = amount of salt in layer at begining of cycle. First layer =
% deposition amount
halite = zeros(sample_pts,1);
halite(1) = dep_halite;

gypsum = zeros(sample_pts,1);
gypsum(1) = dep_gypsum;

% Water holding capacity of each depth layer (0.1cm)
WHC = 0.2*ones(sample_pts,1);
WHC(1) = 0;

% soil_water is the amount of water in each depth layer
soil_water = zeros(sample_pts,1);

% free_water is the amount of water to be moved downward
free_water = zeros(sample_pts,1);
% free_water(1) always = to rainfall amount if soil completely dried out
% prior to rainfall event
free_water(1) = rainfall;

for i = 2:1:sample_pts;
    if free_water(i-1) > WHC(i);
        soil_water(i) = WHC(i);
    else 
        soil_water(i) = free_water(i-1);
    end
    
    free_water(i) = free_water(i-1) - soil_water(i);
    
end

% end salt deposition/water

%%

% start gypsum

% gypsum_soil is the amount of gypsum deposited in the soil at a given
% depth
gypsum_soil = zeros(sample_pts,1);
%gypsum_soil(1) = the amount of gypsm initially on the surface at beginning
%of cycle
gypsum_soil(1) = gypsum(1);



% gypsum_water is the amount of gypsum in the soil water at a given depth
gypsum_water = zeros(sample_pts,1);

% amount of gypsum on surface
if (((rainfall*10) - (WHC(1)*10)) * ec_gypsum) <= gypsum_soil(1);
    gypsum_water(1) = (((rainfall*10) - (WHC(1)*10)) * ec_gypsum);
else 
    gypsum_water(1) = gypsum_soil(1);
end

% amount of gypsum in first level (depth 2) of soil water
if ((rainfall*10) - ((soil_water(1))*10)) <= 0;
    gypsum_water(2) = 0;
else
    if (rainfall*10*ec_gypsum) >= gypsum(2);
        gypsum_water(2) = gypsum_water(1);
    else
        gypsum_water(2) = (gypsum_water(1) + (((rainfall*10) - ((soil_water(1))*10))* ec_gypsum));
    end
end



% amount of gypsum deposited in soil in first layer (depth = 2)
% need to calculate denominator of 2nd if statement first
if (rainfall - soil_water(1)) == 0;
   denominator_gypsum_soil(2) = soil_water(2);
else
   denominator_gypsum_soil(2) = (rainfall - soil_water(1));
end


% begin calculating 'gypsum_soil(2)' with 'denominator_gypsum_soil(2)'
% v_g_s_2 = 'value gypsum soil 2' = value needed to calculate
% 'gypsum_soil(2)'
if (((rainfall - soil_water(1)) * 10 * ec_gypsum) - gypsum_water(1)) >= gypsum(2);
   v_g_s_2 = 0;
else
   v_g_s_2 = (gypsum(2) - ((rainfall*10*ec_gypsum) - gypsum_water(1)));  
end

gypsum_soil(2) = v_g_s_2 + (gypsum_water(2)*(soil_water(2)/denominator_gypsum_soil(2)));


%%

% amount of gypsum in soil water for depth levels 3+
for j = 3:1:sample_pts;
    
if ((rainfall - (sum(soil_water(1:j-1))))*10) <= 0
   gypsum_water(j) = 0;
else
    if ((rainfall*10) - ((sum(soil_water(1:j-1)))*10)) * ec_gypsum - (gypsum_water(j-1) + gypsum_soil(j-1)) >= gypsum(j);
       gypsum_water(j) = (gypsum_water(j-1) - gypsum_soil(j-1) + gypsum(j));
    else
       gypsum_water(j) = (  (rainfall*10) - (sum(soil_water(1:j-1)))*10  )*ec_gypsum;
    end
end

% gypsum_soil(layers3+)
% amount of gypsum deposited in soil in second through sample_pts layers
% (layers 3+)
% need to calculate denominator first
if (rainfall - (sum(soil_water(1:j-1)))) == 0;
   denominator_gypsum_soil(j) = soil_water(j-1);
else
   denominator_gypsum_soil(j) = rainfall - (sum(soil_water(1:j-1)));
end

% begin calculating 'gypsum_soil(j)' with 'denominator_gypsum_soil(3+)'
% where j = depth layers 2 through sample_pts
% v_g_s_j = 'value gypsum soil j+' = value needed to calculate
% 'gypsum_soil(j)'
if ((rainfall - (sum(soil_water(1:j-1)))) * 10 * ec_gypsum) - (gypsum_water(j-2) + gypsum_soil(j-2)) >= gypsum(j-1);
   v_g_s_j = 0;
else
   v_g_s_j = gypsum(j-1) - (((rainfall * 10) - (sum(soil_water(1:j-1))) * 10) * ec_gypsum);
end

gypsum_soil(j) = v_g_s_j + ((gypsum_water(j)) * (soil_water(j) / denominator_gypsum_soil(j)));

end


% end gypsum


%%
% start halite

% halite_soil is the amount of halite deposited in the soil at a given
% depth
halite_soil = zeros(sample_pts,1);
%halite_soil(1) = the amount of halite initially on the surface at beginning
%of cycle
halite_soil(1) = halite(1);



% halite_water is the amount of halite in the soil water at a given depth
halite_water = zeros(sample_pts,1);

% amount of halite on surface
if (((rainfall*10) - (WHC(1)*10)) * ec_halite) <= halite_soil(1);
    halite_water(1) = (((rainfall*10) - (WHC(1)*10)) * ec_halite);
else 
    halite_water(1) = (halite_soil(1));
end

% amount of halite in first level (depth 2) of soil water
% h1 = 0
if ((rainfall*10) - ((soil_water(1))*10)) <= 0;
    halite_water(2) = 0;
else
    if (rainfall*10*ec_halite) >= halite(2);
        halite_water(2) = halite_water(1);
    else
        halite_water(2) = (halite_water(1) + (((rainfall*10) - ((soil_water(1))*10))* ec_halite));
    end
end



% amount of halite deposited in soil in first layer (depth = 2)
% need to calculate denominator of 2nd if statement first
if (rainfall - soil_water(1)) == 0;
   denominator_halite_soil(2) = soil_water(2);
else
   denominator_halite_soil(2) = (rainfall - soil_water(1));
end


% begin calculating 'halite_soil(2)' with 'denominator_halite_soil(2)'
% v_h_s_2 = 'value halite soil 2' = value needed to calculate
% 'halite_soil(2)'
if (((rainfall - soil_water(1)) * 10 * ec_halite) - halite_water(1)) >= halite(2);
   v_h_s_2 = 0;
else
   v_h_s_2 = (halite(2) - ((rainfall*10*ec_halite) - halite_water(1)));  
end

halite_soil(2) = v_h_s_2 + (halite_water(2)*(soil_water(2)/denominator_halite_soil(2)));


%%

% amount of halite in soil water for depth levels 3+
for k = 3:1:sample_pts;
    
if ((rainfall - (sum(soil_water(1:k-1))))*10) <= 0
   halite_water(k) = 0;
else
    if ((rainfall*10) - ((sum(soil_water(1:k-1)))*10)) * ec_halite - (halite_water(k-1) + halite_soil(k-1)) >= halite(k);
       halite_water(k) = (halite_water(k-1) - halite_soil(k-1) + halite(k));
    else
       halite_water(k) = (  (rainfall*10) - (sum(soil_water(1:k-1)))*10  )*ec_halite;
    end
end

% halite_soil(layers3+)
% amount of halite deposited in soil in second through sample_pts layers
% (layers 3+)
% need to calculate denominator first
if (rainfall - (sum(soil_water(1:k-1)))) == 0;
   denominator_halite_soil(k) = soil_water(k-1);
else
   denominator_halite_soil(k) = rainfall - (sum(soil_water(1:k-1)));
end

% begin calculating 'halite_soil(k)' with 'denominator_halite_soil(3+)'
% where k = depth layers 2 through sample_pts
% v_h_s_k = 'value halite soil k+' = value needed to calculate
% 'halite_soil(k)'
if ((rainfall - (sum(soil_water(1:k-1)))) * 10 * ec_halite) - (halite_water(k-2) + halite_soil(k-2)) >= halite(k-1);
   v_h_s_k = 0;
else
   v_h_s_k = halite(k-1) - (((rainfall * 10) - (sum(soil_water(1:k-1))) * 10) * ec_halite);
end

halite_soil(k) = v_h_s_k + ((halite_water(k)) * (soil_water(k) / denominator_halite_soil(k)));

end


% end halite
%%

% % Show individual vectors
% soil_water
% free_water
% gypsum_water
% gypsum_soil
% halite_water
% halite_soil



% % Plot of deposited salts
% plot(halite_soil, max_depth - depth, 'r-', gypsum_soil, max_depth-depth,'b-')
% title('Halite (red) and Gypsum (blue) with depth','FontSize',18);
% xlabel('salt concentration (in grams)', 'FontSize',14);
% ylabel('depth','FontSize',14);