clear;
% close all;
%%This script follows the steps in Gilad and Hardenburgh, 2006 for
%%calculating quantites determined by convolution integrals with
%%space-variant kernels. In our case, the flux goes as
%%int(E(x')f(x-x',x')dx', integrated from zero to x. A true convolution
%%integral would not have a spatial dependence on x'. To get away from
%%this, Gilad and Hardenburg present a method that approximates the a true
%%kernel p(x-x') from a series of kernels with the same form as f(x-x')
%%multiplied by a weighting function alpha. This essentially takes the
%%spatial dependence out of the kernel and places it on the weighting
%%function alpha. This then allows us to redefine the kernel as the
%%approximated kernel, K(x-x'), and another function g(x')=E(x')alpha(x').
%%The brute force method of calculating the flux requires N^2 calculations
%%at least. This method requires N(l)Nlog(N) calculations where N(l) is
%%the number of functions combined to make up the approximating kernel, and
%%N is the number of nodes on the domain. So long as N(l)<0)=0;
%Calculate phi values which are the kernel for specific decay rates
S=linspace(-1.15,0,2);
Sp=linspace(0,1.15,2);
phiP=(1/L0)*(2*Sc./(Sc-S)-1);%create vector of decay rates (disentrainment rates)
phiN=(1/L0)*(2*Sc./(Sc+S)-1);%create vector of decay rates for negative flux if needed.
phiP=1./phiP;
phiN=1./phiN;
% phiL=linspace(min(phiP),max(phiP),16);%create vector of decay rates used to approximate the flux
phiTemp=min(phiP);
phiTempN=max(phiN);
phiL=0;
i=1;
while phiL <= max(phiP);
phiTemp=phiTemp+0.5*phiTemp;
phiLN(i)=phiTempN -0.5*phiTempN;
phiL(i)=phiTemp;
i=i+1;
end
%%initialize matrices for speed
alpha=zeros(length(phiL),lZ);%weights that will vary as a function of position that weight the different kernels
alpha_1=zeros(length(phiL),lZ);%A temporary matrix used in determining alpha
alpha_2=zeros(length(phiL));%A temporary matrix used in determining alpha
A=zeros(1,length(phiL));%A temporary matrix used in determining alpha
kernTempP=zeros(length(phiL),lZ);%A matrix that contains the kernels of all approximating functions at different positions
g=zeros(size(alpha));%a matrix that becomes one of the convolved functions
% phiL=1./phiL;
t=0;
%Calculate kernels for different decay rates (phiL) at all positions.
for k=1:lZ;
for i=1:length(phiL);
kernTempP(i,k)=exp(-(k*dx - dx)/phiL(i));
end
end
fftKern=fft(kernTempP,[],2);
SComp=linspace(0,-1.15,11501);%Create vector of slopes from zero to near the critical slope
phiComp=((1/L0)*(2*Sc./(Sc-SComp)-1));%Calculate the decay rate for all slopes
phiCompN=((1/L0)*(2*Sc./(Sc+SComp)-1));
phiComp=1./phiComp;%take one over the value to make it consistent with decay rates.
phiCompN=1./phiCompN;
%%This loop attempts to solve a system of linear equations to determine values for
%%alpha at every position x.
for k=1:length(phiComp);
% Temp=zeros(length(phiL),1);
for j=1:length(phiL);
b(j,1)=1*phiComp(k)*phiL(j)/(phiComp(k) + phiL(j));
for i=1:length(phiL);
M(j,i)=2*phiL(i)*phiL(j)/(phiL(i) + phiL(j));
end
end
Temp=M\b;
alpha(:,k)=Temp(:);
end
%%
%%Time loop
tic
while t<1000;
t=t+1;
%%Calculate slope
slp=H*z;
slp(lZ)=slp(lZ-1);
slp(slp>0)=0;
E=(E0 + E1*abs(slp).^2);
PhiInd=round(abs(slp*(10^4)))+1;
%This turns the slope value into an
%index such that it immediately tells where to go in a matrix of
%weights, alpha, that weight the respective kernels. For example, this
%turns a slope of -0.7 into an index of 7001, which then can be used to
%get the alpha value of each kernel at each position.
for k=1:lZ;
g(:,k)=E(k)*alpha(:,PhiInd(k));%determine one of the functions to be convolved
end
fftG=fft(g,[],2);%take its fft
fftD=(fftKern.*fftG);%multiply the two fourier transformed functions
Q=ifft(sum(fftD,1));%sum the approximating kernels and take the inverse fft
Qreal=real(Q);%the flux is the real part (there only is a real part)
dz=-H2*Qreal';%Divergence of the flux
z=z+dz;%change in elevation
if rem(t,100)==0;
plot(z,'-k');
hold on
plot(zInit);
axis image
drawnow
hold off
end
end
toc