function [J, model, ac_mod] = costfunc(q, beam, input, data);

%
% function [J, model, ac_mod] = costfunc(q, beam, input, data);
%
% This file you may want to modify to play with different cost functionals.
% It sets up parameters for beam_model_demo, calls that code, and compares
% the results
%
% Inputs:
%   q -- vector of parameters, in the following order:
%     1.  gamma -- coefficient of air damping
%     2.  YI for beam -- Young's modulus times moment of intertia
%     3.  CI for beam -- internal damping
%     4.  f0 for beam 1, Kp for beam 2
%     5.  rho for patch -- linear density of patch
%     6.  YI for patch
%     7.  CI for patch
%     8.  rho for accelerometer -- beam 1 only
%     9.  YI for accelerometer -- beam 1 only
%     10. CI for accelerometer -- beam 1 only
%     NOTE: rho for beam has been made a constant, and all other parameters
%       are chosen in relation to it.  This removes a simple source of 
%       nonuniqueness in the parameter set.
%   beam -- 1 for the beam excited by the hammer, 2 for the beam excited by
%     the patch
%   input -- force input as measured by the hammer for beam 1, voltage
%     input for beam 2
%   data -- data being compared to the model -- this is acceleration for
%     beam 1 and displacement for beam 2
% 
% Outputs:
%   J -- cost for the given parameters 
%   model -- displacement results for the given parameters
%   ac_mod -- acceleration results for the given parametes (beam 1 only)
%
% Tom Braun
% North Carolina State University
% Department of Mathematics/Center for Research in Scientific Computing
% April 2005
%

% disallow bad parameters.  This shouldn't matter for DIRECT, but does for
% fminsearch, lsqnonlin
if length(q) < 10
    q(length(q)+1:10) = 0;
end
q(q < 0) = 0;
if (q(2) < 1e-8)
    q(2) = 1e-8;
end
if (q(4) < 1e-8);
    q(4) = 1e-8;
end

if beam == 1 % hammer impact beam with accelerometers -- has no voltage input
    rho = [0.071209; q(5); q(8)];
    gamma = q(1);
    YI = [q(2); q(6); q(9)];
    CI = [q(3); q(7); q(10)];
    beamlen = 0.287;
    f0 = q(4);
    kp = 0;
    patch_loc = [0.02; 0.046];
    accel_loc = [0.26; 0.275];
    force_point = [0.0979; 0.0981];
    Delta_t = 1/256;
    observe = sum(accel_loc) / 2;
    zeroforcingafter = 0.1;
    force = input;
    voltage = zeros(size(input));
else % patch excited beam -- has no force input and no accelerometer
    rho = [0.072184; q(5); 0];
    gamma = q(1);
    YI = [q(2); q(6); 0];
    CI = [q(3); q(7); 0];
    beamlen = 0.393;
    f0 = 0;
    kp = q(4);
    patch_loc = [0.041, 0.092];
    accel_loc = [0, beamlen]; % this is a dummy placeholder
    force_point = [0, beamlen]; % this is a dummy placeholder
    Delta_t = 0.001;
    observe = 0.128;
    zeroforcingafter = 0.01;
    force = zeros(size(input));
    voltage = input;
end

model = beam_model_demo(beamlen, rho, YI, CI, gamma, f0, kp, patch_loc, accel_loc, force_point, ...
    16, Delta_t, observe, force, voltage, zeroforcingafter);

fftsize = 2^(ceil(log2(length(data))));
if beam == 1
    ac_mod = (model(3:end) - 2*model(2:end-1) + model(1:end-2)) ./ Delta_t^2;
    J = 2.*sum(abs(data(230:end-1) - ac_mod(229:end)).^2) + sum(abs(abs(fft(data(230:end-1), fftsize)) - abs(fft(ac_mod(229:end), fftsize))).^2);
else 
    ac_mod = zeros(size(voltage)); % must be present, since this is an optional output.  It's meaningless for this model
    J = 2.*sum(abs(data(:) - model(:)).^2) + sum(abs(abs(fft(data(:), fftsize)) - abs(fft(model(:), fftsize))).^2);
end