function [t,tt,ts,s] = mytrend(x,start,options)

freq = datfreq(start);
if islogical(options.season)
    if options.season
        s = freq;
    else
        s = [];
    end
else
    s = options.season;
end
if isnumericscalar(s) && (~any(s == [2,4,6,12]))
    s = [];
end

% Break points.
nper = size(x,1);
bp = round(options.break - start + 1);
bp = bp(:).';
bp(bp < 1 | bp > nper) = [];
bp = sort(bp);

% Logarithm requested by the user.
if options.log
    x = log(x);
end

if options.diff
    [tt,ts] = difftrend(x,bp,s,options);
    t = tt + ts;
else
    [tt,ts] = leveltrend(x,bp,s);
    t = tt + ts;
end

% Delogarithmise afterwards.
if options.log
    t = exp(t);
    if nargout > 1
        tt = exp(tt);
        ts = exp(ts);
    end
end

end

% Subfunctions.

%**************************************************************************
function [tt,ts] = leveltrend(xdata,bp,s)

[nper,nx] = size(xdata);
tt = nan([nper,nx]);
ts = nan([nper,nx]);
% Time line with breaks.
M = timeline(nper,bp,false);
% Matrix of seaonal factors.
S = season(nper,s);
X = [M,S];
nm = size(M,2);
ns = size(S,2);
for i = 1 : nx
    sample = getsample(xdata(:,i));
    b = X(sample,:) \ xdata(sample,i);
    if any(isnan(b))
        continue
    end
    tt(sample,i) = M(sample,:) * b(1:nm);
    if ns > 0
        ts(sample,i) = S(sample,:) * b(nm+1:end);
    else
        ts(sample,i) = 0;
    end
end

end
% leveltrend().

%**************************************************************************
function [tt,ts] = difftrend(x,bp,s,options)

dx = x(2:end,:) - x(1:end-1,:);
[nper,nx] = size(x);
tt = nan([nper,nx]);
ts = nan([nper,nx]);
% Time line with breaks.
M = timeline(nper-1,bp,true);
% Matrix of seaonal factors.
S = season(nper-1,s);
X = [M,S];
nm = size(M,2);
ns = size(S,2);
for i = 1 : nx
    sample = getsample(dx(:,i));
    first = find(sample,1);
    last = find(sample,1,'last');
    b = X(sample,:) \ dx(sample,i);
    if any(isnan(b))
        continue
    end
    dtt = M(sample,:) * b(1:nm);
    if ns > 0
        dts = S(sample,:) * b(nm+1:end);
    else
        dts = zeros([nper-1,1]);
    end
    tt(first,i) = x(first,i);
    tt(first+1:last+1,i) = dtt;
    tt(first:last+1,i) = cumsum(tt(first:last+1,i));
    if ~options.connect
        tt(first:last+1,i) = tt(first:last+1,i) + ...
            mean(x(first:last+1,i) - tt(first:last+1,i));
    end
    if ns > 0
        ts(first,i) = 0;
        ts(first+1:last+1,i) = cumsum(dts);
        if ~options.connect
            ts(first:last+1,i) = ts(first:last+1,i) ...
                - mean(ts(first:last+1,i));
        end
    else
        ts(first:last+1,i) = 0;
    end
end

end
% End of subfunction difftrend().

%**************************************************************************
function S = season(nper,s)

if ~isempty(s)
    S = zeros(nper,s-1);
    for i = 1 : s-1
        S(i:s:end,i) = 1;
    end
    S(s:s:end,:) = -1;
else
    S = zeros(nper,0);
end

end
% season().

%**************************************************************************
function M = timeline(nper,bp,diff)
% Time line with break points.

nbp = length(bp);
if diff
    x = ones([nper,1]);
    M = [x,zeros([nper,nbp])];
    k = 1;
else
    x = (0 : nper-1).';
    M = [ones([nper,1]),x,zeros([nper,nbp])];
    k = 2;
end
for i = 1 : nbp
    M(bp(i):end,k+i) = x(1 : nper-bp(i)+1);
end

end
% timeline().