function w_obs = beam_model(beamlength, rho, EI, cD, gamma, Kp, patch_loc, N, Delta_t, obs_point, voltage);

%
% function w_obs = beam_model(beamlength, rho, EI, cD, gamma, Kp, patch_loc, N, Delta_t, obs_point, voltage);
%
% This function models the movement of a beam forced by a PZT patch as a function
% of the beam properties and voltage input.  The movement of the beam is given at
% a user-specified point.  The model is run from time = 0 until 
% Delta_t * (length(voltage) - 1).  The beam is assumed to start at rest.
%
% Note that a variable time step solver is used, with linear interpolation
% between the force points.  If the force/voltage start of nearly constant,
% then it is possible the solver will skip over any regions of interest.
% In this case, the solver should be modified to set the 'MaxStep'
% parameter to some fraction of Delta_t.  However, doing so was found to 
% significantly increase computation time for some models, and is
% unnecessary if the force/voltage start off with a spike or similar.
% Thus, this parameter is not set in the current code.  Another result of
% this is that run-time will vary greatly depending on parameters being
% used.  Finally, linear interpolation is used between time samples of the
% voltage data.  This also slows the model, and thus it is disabled when
% and if it is no longer needed (i.e. when the voltage input is zero).  
% to take advantage of the speed increase, the voltage must be truly zero
% and not 1e-10 or similar.  
%
% Input Parameters (note: all units are assumed to be base metric, i.e. 
%        meters, kilograms, seconds, etc): 
%    beamlength -- length of the beam being modelled, in meters. 
%    rho -- 2 element vector containing linear densities, the first element
%      is the beamnad the second is the PZT patch.  The beam rho must be positive.  
%      The other may be positive or 0.
%    EI -- 2 element vector containing Young's modulus times moment of
%      inertia for the beam and patch, (in the same form as rho.
%      This may be called YI in the literature.
%    cD -- 2 element vector containing Kelvin-Voigt damping parameters for
%      the beam and the patch, respectively.  
%    gamma -- coefficient of air damping
%    Kp -- scaling factor for voltage input
%    patch_loc -- 2 element vector containing the left and right edge
%      locations of the PZT patch, where 0 is the location of the clamp
%    N -- number of finite elements to break the beam into
%    Delta_t -- time between samples for force and voltage.  This has no
%      effect on the differential equation solver.
%    obs_point -- location where beam movement is measured.  
%    voltage -- vector of voltage input terms. 
%
% Output:
%   vector containing displacement in meters at obs_point
%
% Tom Braun
% North Carolina State University
% Department of Mathematics/Center for Research in Scientific Computing
% April 2005
% Modified May 2005 for NCSU undergraduate workshop
% Based heavily on code taken from Ralph Smith's web site:
% http://www4.ncsu.edu/~rsmith/Smart_Material_Systems/Chapter8/Beam_Code/
%

[M_beam, K_beam, C_beam] = matrix_construct_beam(N, 0, beamlength, beamlength, rho(1), EI(1), cD(1), gamma);
[M_patch, K_patch, C_patch] = matrix_construct_beam(N, patch_loc(1), patch_loc(2), beamlength, rho(2), EI(2), cD(2), 0);
    
M = M_beam + M_patch;
K = K_beam + K_patch;
C = C_beam + C_patch;

B = force_beam(N, patch_loc(1), patch_loc(2), beamlength, Kp, 2);
tmp = find(voltage ~= 0);
tmp = tmp(end);
nointerp = Delta_t * (tmp + 1);

A_mat = [zeros(size(K)), eye(size(K)); -inv(M)*K, -inv(M)*C];
B_vec = [zeros(size(B)); inv(M)*B];
 
options = odeset('reltol', 1e-3, 'abstol', 1e-6, 'Jacobian', A_mat);
times = [0:Delta_t:(length(voltage)-1) * Delta_t];
[t,w] = ode15s(@yprime_beam, times, zeros(2*(N+1), 1), options, A_mat, B_vec, voltage, times, nointerp);

% this gives the value at whatever the observation point is
bvn = bevaluate_beam(N, 1, obs_point, beamlength./N); % cubic splines evaluated at x
w_obs = w(:, 1:N+1) * bvn';