% ========================================================================
% Beam example with boundary control
%
% written by: Binh Lieu
% ========================================================================

%% System parameters:
omval = logspace(5,7,100);  % temporal frequencies
 
% beam parameters
alpha = 5e-8; mu = 1.88e-7; ell = 240e-6; EI = 7.55e-12;
 
% convenience variables
a4 = EI/(ell^4); a3 = EI/(ell^3);
 
dom = domain(0,1);  % domain, y \in [0, 1]
y = chebfun('y',dom);
D1 = diff(dom,1); % first-order differentiation operator

% fone(y) = 1; fzero(y) = 0;
fone = chebfun(1,dom); fzero = chebfun(0,dom);
 
%% Boundary conditions on the right end:
Wa = [1, 0, 0, 0; 0, 1, 0, 0; zeros(2,4)];
r = [0; 0; 0; (1/a3)]; % non-homogeneous boundary conditions
 
%% Computing frequency responses
H1 = zeros(length(omval),1); 
KE = zeros(length(omval),1); 
PE = zeros(length(omval),1);
for indom = 1:length(omval)
    
    om = omval(indom);
    
    % boundary conditions on the left end
    Wb = [zeros(2,4); ...
        0, 0, 1, 0; ...
        0, 0, 0, (1 + alpha*1i*om)];
 
    % system coefficients
    P = [-mu*om*om*fone, fzero, fzero, fzero, a4*(1 + alpha*1i*om)*fone];
    
    % solving integral boundary value problem
    [f,fder] = intbvp(P,fzero,Wa,Wb,r);
    
    % Frequency responses
    % with the output g(t) = phi(1,t)
    H1(indom) = f(1);
    
    % with the output f(y,t) = \phi(y,t)
    % kinetic energy: KE = 0.5*integrate_{0}^{1} (f''(y))^2 dy
    KE(indom) = 0.5*norm(D1*D1*f)^2;
    % potential energy: PE = 0.5*om*om*integrate_{0}^{1} f(y)^2 dy
    PE(indom) = 0.5*norm(om*f)^2;
    
    [indom]
    
end


%% Plotting the magnitude and phase of H1
loglog(omval,abs(H1),'b-','LineWidth',1.1)
xlab = xlabel('\omega', 'interpreter', 'tex');
set(xlab, 'FontName', 'cmmi10', 'FontSize', 20);
h = get(gcf,'CurrentAxes');
set(h, 'FontName', 'cmr10', 'FontSize', 15);
 
semilogx(omval,unwrap(angle(H1)),'b-','LineWidth',1.1)
xlab = xlabel('\omega', 'interpreter', 'tex');
set(xlab, 'FontName', 'cmmi10', 'FontSize', 20);
h = get(gcf,'CurrentAxes');
set(h, 'FontName', 'cmr10', 'FontSize', 15);

% Plotting the kinetic and potential energy of the beam

loglog(omval,KE,'b-',omval,PE,'r-','LineWidth',1.1)
xlab = xlabel('\omega', 'interpreter', 'tex');
set(xlab, 'FontName', 'cmmi10', 'FontSize', 20);
h = get(gcf,'CurrentAxes');
set(h, 'FontName', 'cmr10', 'FontSize', 15);

