response(sys)

function [To W dW Z dZ] = response(sys)
%RESPONSE  Computes the midpoint displacement and velocity of a given OCCD
%system.
%    [t W dW Z dZ] = RESPONSE(SYS) where SYS is the structure array as given by
%    OCSET, returns:
%       t is the n-by-1 vector of time.
%		W is the n-by-1 vector of total midpoint displacement.
%		dW is the n-by-1 vector of total midpoint velocity.
%		Z is the n-by-m vector of midpoint displacement for m modes.
%		dZ is the n-by-m vector of midpoint velocity for m modes.


if sys.T>0
    % Call delayed RESPONSE function
    sys.To = sys.To(end); % Only need terminal time (not vector)
    [To W dW Z dZ] = del_response(sys);
    return
end
% Extract parameters
Nt = length(sys.To); % Number of discrete time terms
Nx = length(sys.xo); % Number of initial modes specified
%%Adjust inputs to size Nt-by-Nx
n = 1:Nx;
To = repmat(sys.To(:), 1, Nx);
xo = repmat(sys.xo(:)', Nt, 1);
vo = repmat(sys.vo(:)', Nt, 1);
a = repmat(sys.a(:)', Nt, 1);
b = repmat(sys.b(:)', Nt, 1);
N = repmat(n, Nt, 1);
cp = sys.cp; cd = sys.cd;

% Compute dependent open-loop relations
[R u Ro uo] = ComputeRu;


%Compute modal displacement Z = zh + zq
Z = ((vo - R.*xo)./u.*sin(u.*To) + xo.*cos(u.*To)).*exp(R.*To) + ...
    exp(Ro.*To)./2.*(-b.*cos(uo.*To).*To.*uo + b.*sin(uo.*To) + a.*sin(uo.*To).*To.*uo.^2)./uo.^3;

%Compute modal velocity dZ = dzh + dzq
dZ = ((-xo.*u + (-xo.*R.^2 + vo.*R)./u).*sin(u.*To) + vo.*cos(u.*To)).*exp(R.*To) + ...
    ((-1./2.*Ro.*b.*To./uo.^2+1./2.*a.*To).*cos(uo.*To)+(1./2.*Ro.*(b+a.*To.*uo.^2)./uo.^3+1./2.*(b.*uo.^2.*To+a.*uo.^2)./uo.^3).*sin(uo.*To)).*exp(Ro.*To);

w = sin(n'*pi/2);%Generate Eigenfunctions

W = Z*w; % Calculate total displacement
dW = dZ*w; % Calculate total velocity
To = To(:,1);

    function [R u Ro uo] = ComputeRu
        % Compute dependent open-loop relations
        wo = (N*pi).^2*sys.Alpha;
        R = -sys.Gamma*wo - cd/2;
        u = sqrt(wo.^2 + cp - R.^2);
        % Compute open loop relations
        if sys.IOL
            Ro = -sys.Gamma*wo;
            uo = sqrt(wo.^2 - Ro.^2);
        else
            Ro = R;
            uo = u;
        end
    end
end


%%% DELAYED SUBFUNCTION

function [T W dW Z dZ] = del_response(sys)
n = sys.M;
N = 1:length(sys.xo);
wo = (N*pi).^2*sys.Alpha;
R = -sys.Gamma*wo - sys.cd/2;
u = sqrt(wo.^2 + sys.cp - R.^2);
lambda = -2*sys.Gamma*wo;

% Compute open loop relations
if sys.IOL
    Ro = -sys.Gamma*wo;
    uo = sqrt(wo.^2 - Ro.^2);
else
    Ro = R;
    uo = u;
end
% opt = ddeset('RelTol',1e-20,'initialstep',sys.T/5,'NormControl','on','AbsTol',1e-20);
T = linspace(0,sys.To,n)';
for n = 1:length(sys.xo)
    sol = dde23(@(T, X, Z)ddex1de(T,X,Z, sys, n), sys.T, [sys.xo(n); sys.vo(n)], [0 sys.To]);
    dsoln = deval(sol, T);
    Z(:,n) = dsoln(1, :)';
    dZ(:,n) = dsoln(2, :)';
end

w = sin((1:length(sys.xo))'*pi/2);%Generate Eigenfunctions
W = Z*w;
dW = dZ*w;%Calculate total disp. and vel.

    function dydt = ddex1de(t,x,Z, sys, n)
        % Differential equations function for DDEX1.
        q = exp(Ro(n)*t)*(sys.a(n)*cos(uo(n)*t) + sys.b(n)/uo(n)*sin(uo(n)*t));
        dydt = [ x(2);
            q-2*sys.Gamma*wo(n)*x(2)-wo(n).^2*x(1)-sys.cp*Z(1)-sys.cd*Z(2)];
    end
end