function F=variable_beam(L,number,N,time,f,g,beta,gamma,avi)
%function F=variable_beam(L,number,N,time,f,g,beta,gamma,avi)
%
%This script file computes the motion of a non-uniform beam, generating an
%avi movie file if required. It requires beam_ev.m, beam_em.m, and beam_em_dd.m to run.
%Authors: Dr Mike Meylan, Michael Smith, University of Auckland.
%
%where  L=half the distance of the beam
%       number=number of modes required
%       N=number of subintervals is 2*N+1
%       time=time range of interest
%       f,g=initial conditions as on webpage
%       beta is the non-dimensionalised flexural rigidity function (D(x)=EI)
%       gamma is the non-dimensionalised linear mass density function (m(x)=rho A)
%       avi=create avi file (1=Yes,0=No)
%
%Note that: the beam is measured from -L<x<L
%           units of measure are not prescribed for D(x) or m(x) (MPa,GPa,...).
%           the number of subintervals needs to be considerably higher than the number
%           of modes chosen (for sufficient accuracy)
%           to preview the function, type: variable_beam() in the command window
%------ ------ ------ ------ ------ ------ ------ ------ ------ ------ ------ ------ ------ ------

%defaults
if nargin == 0, L=1/2;  end
if nargin <= 1, number=10;  end
if nargin <= 2, N = 100;   end
if nargin <= 3, time=0:0.01:0.5;  end
x=linspace(-L,L,2*N+1);
if nargin <= 5
    f=L*(x.^3/(6*L)-x.^2/2);
    g=0*x;
end
if nargin <=7
    beta=cat(2,5*ones(size(x(1:floor(end/3)))),50*ones(size(x(floor(end/3):floor(2*end/3)-1))),ones(size(x(floor(2*end/3)+1:end))));
    gamma=2*L/3*cat(2,5*ones(size(x(1:floor(end/3)))),50*ones(size(x(floor(end/3):floor(2*end/3)-1))),ones(size(x(floor(2*end/3)+1:end))));
end
if nargin <= 8, avi=0; end
%------ ------ ------ ------ ------ ------ ------ ------ ------ ------ ------ ------ ------ ------

%Rayleigh Ritz method for computing non-uniform eigenfunctions and non-uniform eigenvalues

eigvals = beam_ev(number,L);                            %vector of eigenvalues for uniform beam
eigmodes=beam_em(eigvals,L,x);                          %matrix of uniform eigenfunctions
eigmodes_dd=beam_em_dd(eigvals,L,x);                    %matrix of second derivatives of uniform eigenfunctions

weights=2*ones(1,length(x));                            %vector of weights (Composite Simpson's Rule)
weights(2:2:2*N) = 4;weights(1) = 1;weights(end) = 1;
weights = L/(3*N)*weights;

J_mat=diag(gamma.*weights);
M=eigmodes*J_mat*eigmodes';

H_mat=diag(beta.*weights);
K=eigmodes_dd*H_mat*eigmodes_dd';

[Ai_values,nonuniform_eigvals] = eig(K,M);              %recall though that mu=k_hat^4
nonuniform_eigvals=nonuniform_eigvals.^(1/4);
nonuniform_eigmodes=Ai_values'*eigmodes;                %this gets us a matrix of nonuniform eigenfunctions

%------ ------ ------ ------ ------ ------ ------ ------ ------ ------ ------ ------ ------ ------

%Working out coefficients A_n and B_n

A_n=(f*J_mat*nonuniform_eigmodes')./(diag(nonuniform_eigmodes*J_mat*nonuniform_eigmodes'))';        %contains A_1,A_2,...

nonzero_nu_eigvals=diag(nonuniform_eigvals(3:end,3:end))';                                          %obtaining nonzero nonuniform_eigenvalues
B_n=(g*J_mat*nonuniform_eigmodes')./(diag(nonuniform_eigmodes*J_mat*nonuniform_eigmodes'))';
B_n=B_n(3:end);                                                                                     % contains B_3,B_4,...
B_n=(1./nonzero_nu_eigvals.^2).*B_n;
B_n=cat(2,zeros(1,2),B_n);                                                                          %setting B_1 and B_2 to zero

%------ ------ ------ ------ ------ ------ ------ ------ ------ ------ ------ ------ ------ ------

% Generating image

%hf=figure('Position',[100,100,600,600]);
%rect = get(hf,'Position');
%rect(1:2) = [0 0]; % these commands ensure we capture the actual figure.
delete video.gif
for j = 1:length(time)
    out = A_n*diag(cos(time(j)*diag(nonuniform_eigvals).^2))*nonuniform_eigmodes+B_n*diag(sin(time(j)*diag(nonuniform_eigvals).^2))*nonuniform_eigmodes;
    refresh
    plot(x(1:ceil(end/3)),out(1:ceil(end/3)),'LineWidth',0.5)
    hold on
    plot(x(ceil(end/3)-1:2*end/3),out(ceil(end/3)-1:2*end/3),'LineWidth',2)
    plot(x(2*end/3:end),out(2*end/3:end),'LineWidth',0.5)
    hold off
    axis([-L,L,-L/2,L/2]);
    text(2*L/6,2*L/6,['t = ',num2str(time(j),'%1.2f')],'Fontsize',18)
    title('Variable Beam','Fontsize',18)
    set(gca,'FontSize',18)
    set(gca,'LineWidth',2);
    F(j) = getframe(gcf);
    %gif_add_frame(gca,'video.gif');
end

if avi==1
    if ismac==1
        movie2avi(F, 'variable_beam.avi', 'compression', 'none');
    elseif ispc==1
        movie2avi(F, 'variable_beam.avi', 'compression', 'Cinepak');
    else
        movie2avi(F, 'variable_beam.avi', 'compression', 'none');
    end
end