function [eps_strat,strat,dust] = generate_strat(stratlen, separation, max_thickness, max_dust, min_dust)
%Generates a single stratigraphy based on input strat length, mean
%layer separation distance, max dusty layer thickness, and max/min layer dust content 

layercheck          = 1/separation; %rough period at which there should be dusty layers
eps_ice             = matzler(4e6,175);   %ice permittivity at a given temp
eps_dust            = 8.8 - 1i*((2e-3)*8.8); %dust permittivity
var_ice             = separation/3; %variance in separation distance
var_dusty           = max_thickness/3; %variance for layer thickness

%create first layer
if rand<=0.5 %decide whether first layer is ice or dust
    dust = 0;
    eps_strat = eps_ice; %if it's ice, set ice permmitivity
    strat = (var_ice .* randn) + separation; %randomly decide thickness from normal dist
else
    dust = ((max_dust-min_dust).*rand + min_dust);
    eps_strat = maxwell_garnett(dust,eps_ice,eps_dust); %randomly set dusty layer composition
    %strat = (var_dusty .* randn) + av_thickness; %randomly decide thickness from normal dist
    strat = ( (max_thickness-1) .* rand) + 1;
end

ii=1;
stratcheck = sum(strat);
while stratcheck < stratlen %while the total thickness is less than the desired total
    
    if rand <= layercheck %randomly decide if there is a dusty layer or not
        thick_ice = ii + (rand) - 0.5; %add (-0.5:0.5) to icy layer thickness just so it isn't a whole integer
        %thick_dust = (var_dusty .* randn) + av_thickness; %randomly decide dusty layer thickness from normal dist
        thick_dust = ( (max_thickness-1) .* rand) + 1;
        dustnew =  ((max_dust-min_dust).*rand + min_dust);
        eps_dustlayer = maxwell_garnett(dustnew,eps_ice,eps_dust); %randomly set dusty layer composition
        
        strat = [strat; thick_ice; thick_dust]; %add ice and dusty layers to strat
        eps_strat = [eps_strat; eps_ice; eps_dustlayer]; %add corresponding permittivities to eps_strat
        dust  = [dust; 0; dustnew];

        stratcheck = sum(strat); %new total thickness
        ii=1;
    else
        stratcheck = sum(strat) + ii; %new total thickness
        ii=ii+1; %increment icy layer thickness
    end
end

if sum(strat) < stratlen
    strat = [strat; stratlen-sum(strat)];
    eps_strat = [eps_strat; eps_ice];
else
    strat = stratlen .* (strat./sum(strat));
end

end