function handles = Schlee(Sand,Mud,Gravel)
%
% This script makes a ternary diagram for sand, mud and gravel, assuming
% that mud=(silt+clay).  The plot is based on Schlee (1973), a modification
% on the Shepard (1954) classification system.  The Schlee modification
% utilizes the same ternary plotting style, but adjusts to include the
% larger phi sized portion of sediments.
%
% This script is a modification of the TERNPLOT script by Carl Sandrock (20020827).
% This modification simply adds the labels required for the Schlee diagram.
%
% NOTE: The plotted result may be screen dependent.  If this happens,
% simply use the magnify tool to zoom in on the plot.
%
%        Gravel
%         / \
%        /   \
%    Sand --- Mud
%
%REFRENCES:
%Shepard, F.P., 1954, Nomenclature based on 
%sand-silt-clay ratios: Journal Sedimentary 
%Petrology, v. 24, p. 151-158.
%
%Schlee, John, 1973, Atlantic Continental Shelf and 
%Slope of the United States sediment texture of the northeastern part: 
%U.S. Geological Survey Professional Paper 529-L, 64 p.
%
% Author: Carl Sandrock 20020827
% Modified by: Matthew Arsenault 20060222
%
%To run:  Schlee(Sand,Mud,Gravel)

A = Sand
B = Mud
C = Gravel

xoffset = 0.04;
yoffset = 0.02;

majors = 4;

% Normalization of data
Total = (A+B+C);
fA = A./Total;
fB = B./Total;
fC = 1-(fA+fB);

[x, y] = frac2xy(fA, fC);

% Sort data points in x order
[x, i] = sort(x);
y = y(i);

% get hold state
cax = newplot;
next = lower(get(cax,'NextPlot'));
hold_state = ishold;

% get x-axis text color so grid is in same color
tc = get(cax,'xcolor');
ls = get(cax,'gridlinestyle');

% Hold on to current Text defaults, reset them to the
% Axes' font attributes so tick marks use them.
fAngle  = get(cax, 'DefaultTextFontAngle');
fName   = get(cax, 'DefaultTextFontName');
fSize   = get(cax, 'DefaultTextFontSize');
fWeight = get(cax, 'DefaultTextFontWeight');
fUnits  = get(cax, 'DefaultTextUnits');

set(cax, 'DefaultTextFontAngle',  get(cax, 'FontAngle'), ...
    'DefaultTextFontName',   get(cax, 'FontName'), ...
    'DefaultTextFontSize',   get(cax, 'FontSize'), ...
    'DefaultTextFontWeight', get(cax, 'FontWeight'), ...
    'DefaultTextUnits','data')

% only do grids if hold is off
if ~hold_state
	%plot axis lines
	hold on;
	
    plot ([0 1 0.5 0],[0 0 sin(1/3*pi) 0], 'color', tc, 'linewidth',1,...
                   'handlevisibility','off');
	set(gca, 'visible', 'off');
    

    % plot background if necessary
    if ~isstr(get(cax,'color')),
       patch('xdata', [0 1 0.5 0], 'ydata', [0 0 sin(1/3*pi) 0], ...
             'edgecolor',tc,'facecolor',get(gca,'color'),...
             'handlevisibility','off');
    end
    
	% Generate labels
	majorticks = linspace(0, 1, majors + 1);
    majorticks = majorticks(3:end-2);
	labels = num2str(majorticks'*100);
	
	zerocomp = zeros(1, length(majorticks));
	
	% Plot left labels
	[lxa, lya] = frac2xy(majorticks, 1-majorticks);
	text(lxa-xoffset, lya, labels);
	
end;

% Reset defaults
set(cax, 'DefaultTextFontAngle', fAngle , ...
    'DefaultTextFontName',   fName , ...
    'DefaultTextFontSize',   fSize, ...
    'DefaultTextFontWeight', fWeight, ...
    'DefaultTextUnits',fUnits );

% Set up the intersection points to form the interior boundaries of the
% Shclee plot, as modified from Shepard
%
% Boundary points for the pure phase regions
[c_sac(1), c_sac(2)] = frac2xy(.125,.75);
[c_sc(1), c_sc(2)] = frac2xy(0,.75);
[sa_csa(1), sa_csa(2)] = frac2xy(.0,.0);
[sa_ssa(1), sa_ssa(2)] = frac2xy(0,0);
[s_cs(1), s_cs(2)] = frac2xy(0,0);
[s_sas(1), s_sas(2)] = frac2xy(0,0);
% Boundary points for the intersections with the pure phase regions
[c(1), c(2)] = frac2xy(.125,.75);
[sa(1), sa(2)] = frac2xy(.0,0.0);
[s(1), s(2)] = frac2xy(.0,.0);
% Boundary Points for homogeneous central region
[sasc_c(1), sasc_c(2)] = frac2xy(.2,.6);
[sasc_csa(1), sasc_csa(2)] = frac2xy(.4,.4);
[sasc_sa(1), sasc_sa(2)] = frac2xy(.4,.4);
[sasc_sas(1),sasc_sas(2)] = frac2xy(.4,.2);
[sasc_s(1), sasc_s(2)] = frac2xy(.2,.2);
[sasc_cs(1), sasc_cs(2)] = frac2xy(0,0);
% Boundary points on the midpoints of the triangle sides
[sac_csa(1),sac_csa(2)] = frac2xy(.5,.5);
[ssa_sas(1),ssa_sas(2)] = frac2xy(0,0);
[cs_sc(1),cs_sc(2)] = frac2xy(0,0);
% Boundary points on the bottom 10 percent
[gx(1),gx(2)] = frac2xy(.9,.1);
[gy(1),gy(2)] = frac2xy(.0,.1);
% Done picking boundary points

%Draw the separator lines on the Schlee Diagram

                %Plot Bottom Horizontal line at <10% Gravel
plot([gx(1) gy(1)],[gx(2) gy(2)], 'color', 'k', 'linewidth',1,...
                   'handlevisibility','off');
               %Plot Bottom Horizontal line below pure Gravel
plot([c_sac(1) c_sc(1)],[c_sac(2) c_sc(2)], 'color', 'k', 'linewidth',1,...
                   'handlevisibility','off');
               %Plot upper completely vertical line
plot([sasc_c(1) c(1)],[sasc_c(2) c(2)], 'color', 'k', 'linewidth',1,...
                   'handlevisibility','off');
               %Plot lower left of gravel portion
plot([sac_csa(1) sasc_csa(1)],[sac_csa(2) sasc_csa(2)], 'color', 'k', 'linewidth',1,...
                   'handlevisibility','off');
               %Plot angled line of lower gravel
plot([sasc_c(1) sasc_sa(1)],[sasc_c(2) sasc_sa(2)], 'color', 'k', 'linewidth',1,...
                   'handlevisibility','off');               

% Write in labels

%Axis label: Gravel
gravellocation = [0.168581 0.481926 12.052];
gravellabel = text(gravellocation(1), gravellocation(2), 'Gravel Content (%)');
set(gravellabel,'FontName','Times','FontSize',18,'Rotation',60,'HorizontalAlignment','center');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Apex label: Sand
siltlocation = [0.056579 -0.0376316 12.052];
siltlabel = text(siltlocation(1), siltlocation(2), 'Sand');
set(siltlabel,'FontName','Times','FontSize',18,'HorizontalAlignment','center');
%
%Apex label: Mud
siltlocation = [0.945 -0.0376316 12.052];
siltlabel = text(siltlocation(1), siltlocation(2), 'Mud');
set(siltlabel,'FontName','Times','FontSize',18,'HorizontalAlignment','center');
%
%Apex label: Gravel
gravellocation = [0.508581 0.901926 12.052];
gravellabel = text(gravellocation(1), gravellocation(2), 'Gravel');
set(gravellabel,'FontName','Times','FontSize',18,'HorizontalAlignment','center');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%Inner label: Gravel
gravelonlylocation = [.5 .73];
gravelonly = text(gravelonlylocation(1), gravelonlylocation(2),'Gravel');
set(gravelonly,'FontName','Times','FontSize',10,'HorizontalAlignment','center');
%
% Center Label of 'Gravelly Sediment'
homogenouslocation = [gravelonlylocation(1) .3];
homtext(4) = {'Gravelly'};
homtext(5) = {'Sediment'};
homogeneous = text(homogenouslocation(1), homogenouslocation(2), homtext);
set(homogeneous,'FontName','Times','FontSize',10,'HorizontalAlignment','center');
%
% Label for lower 10%
sandonlylocation = [.50 .045];
sandonly = text(sandonlylocation(1),sandonlylocation(2),'Sand,Silt & Clay (Gravel <10%)');
set(sandonly,'FontName','Times','FontSize',10,'HorizontalAlignment','center');

q = plot(x, y,'bo');
%
if nargout > 0
    hpol = q;
end
if ~hold_state
    set(gca,'dataaspectratio',[1 1 1]), axis off; set(cax,'NextPlot',next);
end

end
%-----------------------------------------%

% Support Functions frac2xy and radians

function [x, y] = frac2xy(fA, fC);
y = fC*sin(radians(60));
x = 1 - fA - y*cot(radians(60));
end

function radians = radians(degrees)
% RADIANS (DEGREES)
% Conversion function from Radians to Degrees.
% Richard Medlock 12-03-2002
radians = ((2*pi)/360)*degrees;
end