function Movie = infinite_basin_variable_h_and_rho(h,p,L,L_outer,F,G,time,accuracy)    %,initial_conditions)
% Movie = infinite_basin_variable_h_and_rho(h,p,L,L_outer,F,G,time,accuracy)
% This code was written to solve the PDE: \rho(x)\patrial_t^2\zeta = \partial_x(h(x)partial_x\zeta)
% over x in (-\inf,\inf) where h(x) and \rho(x) vary only over a finite region - from -L to L and
% are 1 outside this region.
%--------------------------------------------------------------------------
%-------------------INITIALISE VARIABLES-----------------------------------
%--------------------------------------------------------------------------
% if accuracry is set to zero (or not passed) then the solution is calulate
% to a lower accuracy which generally works fine.
%
% Default values
% accuracy = 0
% G = 0
% F = exp(-20*x.^2);
% time = 0:0.01:3;
% L = 1;
% L_outer = 2*L
% p = 1
% h = 1

if nargin == 0
    h = ones(1,1000);
    L = 1;
end
if nargin <= 1
    p = ones(size(h)); 
    L = 1;
end
if  nargin <= 2
    L = 1;
    L_outer = 2;
end
x = linspace(-L,L,length(h));
delta_x = x(2)-x(1);
if nargin <= 3
    L_outer = 2*L;
end

y1 = -L_outer:delta_x:-L-delta_x; % y1 and y2 are the regions of y on each side of of -L to L.
y2 = L+delta_x:delta_x:L_outer;
y = [y1,x,y2].'; % y contains the values of x for which the solution will be calculated

if nargin <= 4
    F = exp(-20*y.^2);
end
if nargin <= 5
    G = zeros(size(y));
end
if nargin <= 6
    time = 0:0.01:3;
end
if nargin <= 7
    accuracy = 0;
end
% Set N - the number of fourier functions to be used to approximate each
% eigenfunction and the number of eigenfunctions that will be used.
N = length(h)/(8*L);

% delta_x depends on N so that enough x points are used to numerically
% integrate the highest frequency fourier functions.


% Set the density \rho as a function of x for -L to L.
Density = diag(p); % a matrix used to weight the inner products
hdiff = gradient(h,x); % differentiate h
H = diag(h); % a matrix used to weight the integrals that give the matrix K

if accuracy == 0
    kmax = pi*N/20;  %kmax = 4*pi*N/5; gives exactly ten x points per wavelength but fewer k values is generally sufficient
    delta_k = 80/length(h); % separation between k values
else
    kmax = pi*N/5;  %kmax = 4*pi*N/5; gives exactly ten x points per wavelength but fewer k values is generally sufficient
    delta_k = 8/length(h); % separation between k values
end

k_vector = linspace(0,kmax,kmax/delta_k+1).';

% if nargin <= 5
%     % the initial conditions - zeta(x,0)=F(x) zeta_t(x,0) = G(x)
%     F = exp(-10/(L_outer-L)*(y+(L_outer+L)*0.5).^2).';
%     G = 20/(L_outer-L)*(y.'+(L_outer+L)*0.5).*F;
% end
% if nargin == 6
%     G = 0*F;
% end
%--------------------------------------------------------------------------
%----------------SIMPSON COEFFICIENT MATRICES------------------------------
%--------------------------------------------------------------------------
simp_x = 2*delta_x/3*ones(1,length(x)); simp_x(2:2:length(x)-1) = 4*delta_x/3; simp_x(1) = delta_x/3; simp_x(end) = delta_x/3;
simp_x = diag(simp_x); % Used to implement Simpson's rule for integrals wrt x over -L to L(the values of x)

simp_y = 2*delta_x/3*ones(1,length(y)); simp_y(2:2:length(y)-1) = 4*delta_x/3; simp_y(1) = delta_x/3; simp_y(end) = delta_x/3;
simp_y = diag(simp_y); % Used to implement Simpson's rule for integrals wrt x over -5L to 5L(the values of y)

simp_k = 4*delta_k/3*ones(1,round(kmax/delta_k)+1); simp_k(1:2:end) = 2*delta_k/3*ones(size(simp_k(1:2:end))); simp_k(1) = delta_k/3;simp_k(end) = delta_k/3;
simp_k = diag(simp_k); % Used to implement Simpson's rule for integrals wrt k over 0 to kmax(the values of k)
%--------------------------------------------------------------------------
%--EIGENVALUES/VECTORS FOR THE HOMOGENEOUS STURM LIOUVILLE PROBLEM FOR u,--
%----------------------WHERE \zeta = \zeta_p + u.--------------------------
%--------------------------------------------------------------------------

% The sinm functions are the fourier functions used to approximate the
% eigenfunctions
sinm = zeros(N,length(x));
for m=1:N
    sinm(m,:) = 1/sqrt(L)*sin(m*pi*0.5*(x/L + 1));
end
M = sinm*Density*simp_x*transpose(sinm); % M contains the weighted inner products of the sine functions

% the cosn functions are the derivatives of the fourier sine functions used
% to approxinate the eigenvectors of the Sturm Liouville problem
cosn = zeros(N,length(x));
for n=1:N
    cosn(n,:) = 1/sqrt(L)*n*pi*0.5*cos(n*pi*0.5*(x/L + 1))/L;
end
K=cosn*simp_x*H*transpose(cosn); % K is a matrix of integrals, satisfying Ka = lambdaMa


[A,D] = eig(K,M); % solves the generalised eigenvalue problem created by
% substituting the fourier approximation for an eigenfunction into the
% variational form of the Sturm Liouville problem and then differentiating
% wrt to the nth fourier coefficient - the derivative must be zero and the
% resulting expression can be written in matrix form as K*a = lambda*M*a
% where a is the vector of fourier coefficients for the eigenfunction
% substituted in. Hence A is a matrix of fourier cofficients for the
% eigenfunctions of the Sturm Liouville problem.

norms = sqrt(diag(transpose(A)*M*A)); % transpose(A)*M*A) is diagonal
normalise = ones(size(norms))./norms;
normalise = diag(normalise);

A = A*normalise; % A now contains the fourier coefficients of the normalised eigenfunctions
u = transpose(sinm)*A;
% u is a matrix of normalised eigenfunctions as column vectors. Self
% adjointness of the spatial operator implies the eigenfunctions are
% orthogonal, ie. u.'*Density*simp_x*u = eye(N), we can write this as
% transpose(A)*sinm*Density*simp_x*transpose(sinm)*A = eye(N), or:
% transpose(A)*M*A = eye(N) - hence we normalised A so that this would be
% true.

%--------------------------------------------------------------------------
%---WE NOW USE THE EIGENFUNCTIONS TO SOLVE THE INHOMOGENEOUS PROBLEM-------
%--------------------------------------------------------------------------

X_left = zeros(round(kmax/delta_k)+1,length(y1)); % will contain exp(i*k*y1)+(R1)*exp(-i*k*y1)
Y_left = zeros(round(kmax/delta_k)+1,length(x)); % will contain the eigenfunctions
% \zeta_left(x,k) of the PDE, where the boundary conditionsa are determined
% by a wave coming in from the left.
Z_left = zeros(round(kmax/delta_k)+1,length(y2)); % will contain (T1)*exp(i*k*y2);

X_right = zeros(round(kmax/delta_k)+1,length(y1)); % (T2)*exp(-i*k*y1);
Y_right = zeros(round(kmax/delta_k)+1,length(x)); % will contain the eigenfunctions
% \zeta_left(x,k) of the PDE, where the boundary conditionsa are determined
% by a wave coming in from the right.
Z_right = zeros(round(kmax/delta_k)+1,length(y2)); % exp(-i*k*y2)+(R2)*exp(i*k*y2);

for l= 1:length(k_vector) %loops through k values from 0 to 16 with separation delta_k
    k = k_vector(l);

    S = zeros(2);
    % initialises a matrix S to map [a; b] to [\partial_x \zeta |_{x=-L} ;
    % \partial_x \zeta |_{x=L} ]

    C = zeros(N,2); % will contain the fourier coefficients of u(x,k) for the BVP a=1,b=0 in the first column and a=0,b=1 in the second

    % We now consider the solution for the case when the lefthand side is one
    % and the righthand side is zero and vice versa.
    for j=1:2

        a = mod(j,2); b = mod(j+1,2);% sets a=1 and b=0 for j=1 and a=0 and b=1 for j=2

        zeta_p = (b-a)*0.5*(x/L + 1)+a; % a straight line particular solution.

        % f is the LHS of the inhomogeneous Sturm-Louiville problem for u
        f = -k^2*zeta_p*Density-(b-a)*hdiff/(2*L);


        if rank(D-k*eye(N)) < N % We ensure k is not an eigenvalue of D to avoid division by zero later
            k=k+1.e-6;
        end
        if k==0
            k=1.e-6;
        end

        B = f*simp_x*u; % B is a row vector, initially containing the
        % integrals of f times each of the eigenfunctions from the
        % homogeneous Sturm Liouville problem

        for n=1:N
            B(n) = B(n)/((k^2-D(n,n))); % B is now a vector of coefficients
            % of the eigenvectors from the homogeneous Sturm Liouville problem
        end

        C(:,j)=A*transpose(B); % c_n = (sum over m) a_(n,m)*b_m

        % Here we evaluate the slope of zeta(x,k)
        % at L and -L. The cosn(:,1).' etc comes from evaluating the derivatives
        % of the fourier functions at L and -L.

        S(1,j) = (b-a)/(2*L) + cosn(:,1).'*C(:,j); % slope at x=-L
        S(2,j) = (b-a)/(2*L) + cosn(:,length(x)).'*C(:,j); % slope at x=L
    end

    % We create and solve a linear system for T and R derived by writing
    % [\partial_x \zeta |_{x=-L}; \partial_x \zeta |_{x=L} ] = S[a,b] and
    % rearranging to solve for R,T.

    %WAVE COMING IN FROM THE LEFT
    U = [S(1,1)+i*k, S(1,2); S(2,1), S(2,2)-i*k]*exp(i*k*L);
    V = [(-S(1,1)+i*k)*exp(-i*k*L); -S(2,1)*exp(-i*k*L)];
    W = linsolve(U,V);
    R1 = W(1); T1 = W(2);

    % Evaluate a and b from T and R
    a_left = exp(-i*k*L)+(R1)*exp(i*k*L); b_left = (T1)*exp(i*k*L);

    % Using the linearity of the problem to add the a_left*(solution for
    % a=1,b=0) plus b_left*(solution for b=1,a=0).
    C_left = a_left*C(:,1)+b_left*C(:,2);

    % Evaluating the calculated soltion in three regions
    X_left(l,:) = exp(i*k*y1)+(R1)*exp(-i*k*y1);
    Y_left(l,:) = (b_left-a_left)*0.5*(x/L + 1)+ a_left +transpose(C_left)*sinm;
    Z_left(l,:) = (T1)*exp(i*k*y2);

    %WAVE COMING IN FROM THE RIGHT
    U = [S(1,2), S(1,1)+i*k; S(2,2)-i*k, S(2,1)]*exp(i*k*L);
    V = [-S(1,2)*exp(-i*k*L); (-S(2,2)-i*k)*exp(-i*k*L)];
    W = linsolve(U,V);
    R2 = W(1); T2 = W(2);

    % Evaluate a and b from T and R
    a_right = (T2)*exp(i*k*L); b_right = exp(-i*k*L)+(R2)*exp(i*k*L);

    % Using the linearity of the problem to add the a_left*(solution for
    % a=1,b=0) plus b_left*(solution for b=1,a=0).
    C_right = a_right*C(:,1)+b_right*C(:,2);

    % Evaluating the calculated soltion in three regions
    X_right(l,:) = (T2)*exp(-i*k*y1);
    Y_right(l,:) = (b_right-a_right)*0.5*(x/L + 1)+a_right+transpose(C_right)*sinm;
    Z_right(l,:) = exp(-i*k*y2)+(R2)*exp(i*k*y2);

end
Depth =[-ones(1,length(y1)), -h, -ones(1,length(y2))]; % So the depth can be plotted
Density = [ones(1,length(y1)), p, ones(1,length(y2))]; % to weight inner products
Density = diag(Density);

%The eigenfunctions zeta(x,k) are stored in A_left and A_right
A_left = [X_left(:,1:length(y1)) Y_left Z_left(:,1:length(y2))];
A_right = [X_right(:,1:length(y1)) Y_right Z_right(:,1:length(y2))];


%--------------------------------------------------------------------------
%-----------------GENERALISED FOURIER TRANSFORM----------------------------
%--------------------------------------------------------------------------

% The inner products of F and G with the eigenfunctions
c1 = conj(A_left)*Density*simp_y*F/(2*pi);
d1 = conj(A_right)*Density*simp_y*F/(2*pi);
c2 = conj(A_left)*Density*simp_y*G/(2*pi);
d2 = conj(A_right)*Density*simp_y*G/(2*pi);

%THIS IMMEDATELY GIVES THE INITIAL SOLUTION : zeta0 =
%transpose(A_left)*simp_k*c1 + transpose(A_right)*simp_k*d1;

% plots the initial conditions
% zeta0 = transpose(A_left)*simp_k*c1 + transpose(A_right)*simp_k*d1;
% plot(y,real(zeta0),y,F,y,Depth)


%Calculating the time dependent solution 
% time = linspace(0,time,15*time);
Solution = zeros(length(time),length(y));
for s = 1:length(time)
    s
    c1t = c1.*cos(k_vector*time(s));
    d1t = d1.*cos(k_vector*time(s));
    c2t = c2.*sin(k_vector*time(s))./(k_vector + eps);
    d2t = d2.*sin(k_vector*time(s))./(k_vector + eps);
    ct = c1t+c2t; dt = d1t+d2t;
    
    Solution(s,:) = transpose(A_left)*simp_k*ct + transpose(A_right)*simp_k*dt; % The inverse 'generalised fourier transform'
end

hf=figure('Position',[100,100,600,600]);
rect = get(hf,'Position'); 
rect(1:2) = [0 0]; % these commands ensure we capture the actual figure. 

%plotting the solution

fontsize = 14;

for count = 1:length(time)
    plot(y,real(Solution(count,:)),'b','LineWidth',2)
    hold on
    plot(y,Depth,'g','LineWidth',1.5)
    axis([-L_outer, L_outer, floor(min(min(real(Solution)))), ceil(max(max(real(Solution))))]);
    text(L_outer*0.7,0.8*ceil(max(max(real(Solution)))),['t = ',num2str(time(count),'%1.1f')],'Fontsize',fontsize)
    set(gca,'FontSize',fontsize)
    xlabel('x','Fontsize',fontsize)
    ylabel('\zeta','Fontsize',fontsize)
    hold off
    drawnow
    Movie(:,count)=getframe(hf,rect);
end


% movie(Movie)

return