Folk_S_Classification(Sand,Silt,Clay)

function handles = Folk_S_Classification(Sand,Silt,Clay)
% This script makes a ternary diagram for sand, silt and clay.  
% The plot is based on Folk (1954), reference below.
%
% This script is a modification of the TERNPLOT 
% script by Carl Sandrock (20020827),
% and the 'Schlee.m' written by Matt Arsenault.
%
% NOTE: The plotted result may be screen dependent.  If this happens,
% simply use the magnify tool to zoom in on the plot.
%         Sand
%         / \
%        /   \
%     Clay--- Silt
%
% Modified by: Matt Arsenault 20060222
%
%REFERENCE:
%Folk, R.L., 1954. The distinction between grain size and mineral 
%composition in sedimentary rock nomenclature. Journal of 
%Geology 62 (4), 344-359.
%
%TO RUN:  Folk_S_Classification(Sand,Silt,Clay)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Intersection with pure sand
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sand Left
[g_left(1), g_left(2)] = frac2xy(.10,.90); 
% Sand Right
[g_right(1), g_right(2)] = frac2xy(0,.90);
% Sand Left Center 
[g_cent(1), g_cent(2)] = frac2xy(.06,.90);
% Sand Right Center 
[g_frac(1),g_frac(2)] = frac2xy(.03,.90);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Mid-plot Intersection at 50%, moving left to right
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[g_leftthirty(1),g_leftthirty(2)] = frac2xy(.5,.5);
[g_rightthirty(1),g_rightthirty(2)] = frac2xy(0,.5);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Mid-plot Intersection at 10%, moving left to right
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[g_leftfive(1), g_leftfive(2)] = frac2xy(.90,.10);
[g_rightfive(1), g_rightfive(2)] = frac2xy(0,.10);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Mid-lines
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Center Left, 2:1 line
[mud_onetone(1), mud_onetone(2)] = frac2xy(.66,.0);
%Center Right, 1:2 line
[mud_ninetone(1), mud_ninetone(2)] = frac2xy(.33,.0);
% Done picking boundary points

A = Clay;
B = Silt;
C = Sand;

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);
	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 );

%Draw all the major horizontal lines
plot([g_left(1) g_right(1)],[g_left(2) g_right(2)], 'color', 'k', 'linewidth',1,...
                   'handlevisibility','off');
plot([g_leftthirty(1) g_rightthirty(1)],[g_leftthirty(2) g_rightthirty(2)], 'color', 'k', 'linewidth',1,...
                   'handlevisibility','off');
plot([g_leftfive(1) g_rightfive(1)],[g_leftfive(2) g_rightfive(2)], 'color', 'k', 'linewidth',1,...
                   'handlevisibility','off');
% The two vertical lines
plot([mud_onetone(1) g_cent(1)],[mud_onetone(2) g_cent(2)], 'color', 'k', 'linewidth',1,...
                   'handlevisibility','off');
plot([g_frac(1) mud_ninetone(1)],[g_frac(2) mud_ninetone(2)], 'color', 'k', 'linewidth',1,...
                   'handlevisibility','off');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Write in the Sand, Silt and Clay  Axis labels
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Sand
sandlocation = [0.168581 0.481926 12.052];
sandlabel = text(sandlocation(1), sandlocation(2), 'Percent of Sand (%)');
set(sandlabel,'FontName','Times','FontSize',18,'Rotation',60,'HorizontalAlignment','center');
%Silt
siltlocation = [0.506579 -0.0776316 12.052];
siltlabel = text(siltlocation(1), siltlocation(2), 'Clay : Silt Ratio');
set(siltlabel,'FontName','Times','FontSize',18,'HorizontalAlignment','center');
%Ratios
twooneratiolocation = [0.33 -0.025 12.052];
twooneratiolabel = text(twooneratiolocation(1), twooneratiolocation(2), '2:1');
set(twooneratiolabel,'FontName','Times','FontSize',12,'HorizontalAlignment','center');
%
onetworatiolocation = [0.66 -0.025 12.052];
onetworatiolabel = text(onetworatiolocation(1), onetworatiolocation(2), '1:2');
set(onetworatiolabel,'FontName','Times','FontSize',12,'HorizontalAlignment','center');
%Clay
claylocation = [0.832 0.481926 12.052];
claylabel = text(claylocation(1), claylocation(2), '');
set(claylabel,'FontName','Times','FontSize',18,'Rotation',-60,'HorizontalAlignment','center');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Apex label: Clay
apex_clay_location = [0.056579 -0.0376316 12.052];
apex_clay_label = text(apex_clay_location(1), apex_clay_location(2), 'Clay');
set(apex_clay_label,'FontName','Times','FontSize',18,'HorizontalAlignment','center');
%Apex label: Silt
apex_silt_location = [0.945 -0.0376316 12.052];
apex_silt_label = text(apex_silt_location(1), apex_silt_location(2), 'Silt');
set(apex_silt_label,'FontName','Times','FontSize',18,'HorizontalAlignment','center');
%Apex label: Sand
apex_sand_location = [0.508581 0.901926 12.052];
apex_sand_label = text(apex_sand_location(1), apex_sand_location(2), 'Sand');
set(apex_sand_label,'FontName','Times','FontSize',18,'HorizontalAlignment','center');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Write in the individual unit names
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sandonlylocation = [.5 .81];
sandonly = text(sandonlylocation(1), sandonlylocation(2),'S');
set(sandonly,'FontName','Times','FontSize',10,'HorizontalAlignment','center');

sandyclaylocation = [.35 .5];
sactext(1) = {'cS'};
sandyclay = text(sandyclaylocation(1),sandyclaylocation(2),sactext);
set(sandyclay,'FontName','Times','FontSize',10,'HorizontalAlignment','center');

muddysandlocation = [.5 sandyclaylocation(2)];
mstext(1) = {'mS'};
muddysand = text(muddysandlocation(1), muddysandlocation(2),mstext);
set(muddysand,'FontName','Times','FontSize',10,'HorizontalAlignment','center');

zslocation = [.65 sandyclaylocation(2)];
sctext(1) = {'zS'};
zs = text(zslocation(1), zslocation(2),sctext);
set(zs,'FontName','Times','FontSize',10,'HorizontalAlignment','center');

sandyclaylocation = [.25 .25];
scatext(1) = {'sC'};
sandyclay = text(sandyclaylocation(1), sandyclaylocation(2), scatext);
set(sandyclay,'FontName','Times','FontSize',10,'HorizontalAlignment','center');

clayeysiltlocation = [.75 sandyclaylocation(2)];
cstext(1) = {'sZ'};
clayeysilt = text(clayeysiltlocation(1), clayeysiltlocation(2), cstext);
set(clayeysilt,'FontName','Times','FontSize',10,'HorizontalAlignment','center');

muddysandtext(5) = {'sM'};
muddysand = text(.5, .28, muddysandtext);
set(muddysand,'FontName','Times','FontSize',10,'HorizontalAlignment','center');

clayonlylocation = [.18 .04];
clayonly = text(clayonlylocation(1),clayonlylocation(2),'C');
set(clayonly,'FontName','Times','FontSize',10,'HorizontalAlignment','center');

siltysandlocation = [.5 .08];
ssatext(2) = {'M'};
siltysand = text(siltysandlocation(1), siltysandlocation(2),ssatext);
set(siltysand,'FontName','Times','FontSize',10,'HorizontalAlignment','center');

siltonlylocation = [.84 .04];
siltonly = text(siltonlylocation(1),siltonlylocation(2),'Z');
set(siltonly,'FontName','Times','FontSize',10,'HorizontalAlignment','center');

% Plot data

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(radianfive_onetnine(60));
x = 1 - fA - y*cot(radianfive_onetnine(60));
end

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