function Movie = finite_basin_variable_h_and_rho(h,p,L,F,G,time,Problem)
% Movie = finite_basin_variable_h_and_rho(h,p,L,F,G,time,Problem)
% program to produce a movie of the solution to the equation 
% \rho \partial_t^2 u - \partial_x^2 (h  \partial_x^2 u) = 0.
%
% h is the variable depth, p is the variable density (\rho) and this is the
% region -L to L. 
% time is a vector of time values to calculate the solution,
% F if the initial value of u 
% G is the initial value of \partial_t u 
% Problem refers to the boundary conditions, which are either u(-L) = u(L)
% = 0 (Dirichlet) or \partial_x u(-L) = \partial_x u(L) = 0 (Neuman)
%
% Default values 
%
% if not specified we have
% Problem = 'a'
% G = 0
% F = exp(-20*x.^2);
% time = 0:0.01:3;
% L = 1;
% p = 1
% h = 1

% set defaults if needed 
if nargin == 0
    h = ones(1,100);h(26:75) =0.5;
    L = 1;
end
if nargin <= 1
    p = ones(size(h)); 
    L = 1;
end
if nargin <= 2
    L = 1 
end
x = linspace(-L,L,length(h));
if nargin <= 3
    F = 0.5*exp(-20*x.^2);
end
if nargin <= 4
    G = zeros(size(x));
end
if nargin <= 5
    time = 0:0.025:3;
end
if nargin <= 6
    Problem = 'a'; % put 'a' for Dirichlet boundary conditions and 'b' (or anything else for Neuman)
end








stepSize = x(2)-x(1);% h for Simpson's rule
N = length(h)/8; % the number of fourier sine/cosine terms used to approximate the eigenfunctions of u
H = diag(h); % a weighting matrix for where h is a product in an integral
P = diag(p); % a weighting matrix for where rho (p) is a product in an integral

W = zeros(N,length(h)); %initialises a matrix to contain fourier sine/cosine functions as row vectors
Wdiff = zeros(N,length(h)); % to contain the derivative of W w.r.t. x row-wise

if Problem == 'a' % fixed string BCs so eigenfunctions are sine waves
    for k = 1:N
        W(k,:) = (1/sqrt(L))*sin(k*pi*(x/(2*L) + 1/2));
        Wdiff(k,:) = (1/sqrt(L))*k*pi*cos(k*pi*(x/(2*L) + 1/2))/(2*L);
    end
else
    k=1; % normal derivative equals zero BCs so cosine waves are used
    W(k,:) = (1/sqrt(2*L))*ones(size(x)); %the constant fourier function
    Wdiff(k,:) = zeros(size(x)); %and its derivative
    for k = 2:N
        W(k,:) = (1/sqrt(L))*cos((k-1)*pi*(x/(2*L) + 1/2)); % fourier cosine functions
        Wdiff(k,:) = -(1/sqrt(L))*(k-1)*pi*sin((k-1)*pi*(x/(2*L) + 1/2))/(2*L); % their derivatives
    end
end


weight = 4*stepSize/3*ones(size(x)); weight(1:2:end) = 2*stepSize/3*ones(size(weight(1:2:end)));
weight(1) = stepSize/3; weight(end) = stepSize/3;
S = diag(weight); % Simpson's coefficient matrix for numerical integration

M = W*S*P*transpose(W); % a matrix of inner products of each of the fourier functions with each other
% - note they are not orthonormal for variable rho (contained in P) so M is not the identity

K = Wdiff*S*H*transpose(Wdiff); % K is a matrix of integrals, satisfying Ka = lambdaMa

[V,Dinitial] = eig(K,M); %solves the generalised eigenvalue/vector problem

norms = sqrt(diag(transpose(V)*M*V)); %transpose(V)*M*V) is diagonal 
normalise = ones(size(norms))./norms;% we normalise V and hence U so that the eigenfunctions are otrthonormal
normalise = diag(normalise);

U = normalise*transpose(V)*W; % U is a matrix of eigenfunctions as row vectors(x increasing along the rows)
A = U*S*P*transpose(F); % A is a vector of inner products of F with the eigenfunctions - ie. the fourier coefficients of the eigenfunctions
B = U*S*P*transpose(G); % B is a vector of inner products of G with the eigenfunctions - ie. the fourier coefficients of the eigenfunctions

for t = 1:length(time)
    for j = 1:N
        if Dinitial(j,j) ~= 0;
            D(j,t) = A(j)*cos(sqrt(Dinitial(j,j))*time(t)) ...
                + B(j)*sin(sqrt(Dinitial(j,j))*time(t))/sqrt(Dinitial(j,j)); % the time dependent part of the solution
        else
            D(j,t) = A(j)*cos(sqrt(Dinitial(j,j))*time(t)) ...
                + B(j);
        end
    end                                                      
end

f_time_estimate = transpose(U)*D; % multiplies each eigenfunction with its time varying coefficient vector D(j,:) and sums over j to give the time dependent solution


%plotting the solution

fontsize = 18;
linewidth = 2;


for count = 1:length(time)
    plot(x,f_time_estimate(:,count),'b','LineWidth',linewidth)
    hold on
    plot(x,-h,'g','LineWidth',linewidth)
    axis([-L L floor(min(min(f_time_estimate))) ceil(max(max(f_time_estimate)))]);
    text(L*0.65,0.8*ceil(max(max(real(f_time_estimate)))),['t = ',num2str(time(count),'%1.1f')],'Fontsize',fontsize)
    set(gca,'FontSize',fontsize)
    set(gca,'LineWidth',linewidth)
    xlabel('x','Fontsize',fontsize)
    ylabel('\zeta','Fontsize',fontsize)
    title('Variable Depth Finite Basin')
    hold off
    drawnow
    Movie(:,count)=getframe(gcf);
    
    % if you want to make a gif undelete the following   
%     gif_add_frame(gcf,'finite_string.gif');
    
end

return