% 
% Copyright (c) 2008 Kristopher L. Kuhlman (klkuhlm at sandia dot gov)
% 
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
% copies of the Software, and to permit persons to whom the Software is
% furnished to do so, subject to the following conditions:
% 
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
% 
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
% THE SOFTWARE.
% 

% Kuhlman, K.L. and A.W. Warrick, 2008. Quasilinear infiltration from 
% an elliptical cavity, Advances in Water Resources, 31(8), 1057–1065. 
% http://dx.doi.org/10.1016/j.advwatres.2008.04.009
%
% this Matlab script computes the Mathieu functions required for the
% elliptical cutout problem

% $Id: elliptical_cutout.m,v 1.20 2008/04/07 22:05:07 kris Exp kris $

%% parameters to set
elA = 1.0;
ecc = 1.0;
N = 10;          % number of terms in generalized Fourier series expansion

%% computed parameters
F = elA*ecc;        % dimensionless semi-focal distance (*alpha/2)
e0 = acosh(1/ecc);      % elliptical radius of cutout
elB = F*sinh(e0);   % dim-less semi-minor dist

% compute solution on an elliptical mesh 
eta = linspace(e0,10.0,50);   % eta values to compute solution at
psi = linspace(-pi,pi,100);      % psi values "   "    "  "

ne = length(eta);
np = length(psi);
R = max(2*N+12,40);   % size of "infinite" matrix
q = (F*F)/4.0;        % Mathieu parameter (remove negative from definition)

%% setup and compute eigenvalue problem

MEv = diag((2*(0:R-1)').^2,0) + diag(ones(R-1,1)*q,-1) + diag(ones(R-1,1)*q,+1);
MEv(2,1) = 2.0*q;
[AEv,LEv] = eig(MEv);

MOd = diag((2*(0:R-1)' + 1).^2,0) + diag(ones(R-1,1)*q,-1) + diag(ones(R-1,1)*q,+1);
MOd(1,1) = 1.0 - q;
[AOd,LOd] = eig(MOd);

% set sign of eigenvalues
for col = 1:R
  if AEv(col,col) < 0.0
    AEv(:,col) = -AEv(:,col);
  end
  if AOd(col,col) < 0.0
    AOd(:,col) = -AOd(:,col);
  end
end

%% normalize eigenvectors (zero term included twice)
% other eigenvalues already normalized correctly by LAPACK
norm = 2*AEv(1,1)^2 + sum(AEv(2:R,1).^2);
AEv(1:R,1) = AEv(1:R,1)/sqrt(norm);

%% compute cosine-elliptic (ce), for each value of psi
ce = zeros(N,np);  se = zeros(N,np);
iv(1:R,1) = (-1.0).^(0:R-1).';

% MF are indexed by n, rather than 2n or 2n+1
% since only even ce and odd se are used
for p = 1:np
  for n = 1:N
    ce(n,p) = iv(n)*sum(iv(1:R).*AEv(1:R,n).*cos(2*(0:R-1)'*psi(p)));
    se(n,p) = iv(n)*sum(iv(1:R).*AOd(1:R,n).*sin((2*(0:R-1)' + 1)*psi(p)));
  end
end

%% compute Ke and Ko for each value of eta
Ke = zeros(N,ne);  Ko = zeros(N,ne);
besi = zeros(R+1,ne);   besk = zeros(R+1,ne);
v = zeros(2,ne);

expeta(1,1:ne) = exp(-eta(1:ne));
expeta(2,1:ne) = exp( eta(1:ne));
v(1,1:ne) = sqrt(q)*expeta(1,1:ne);
v(2,1:ne) = sqrt(q)*expeta(2,1:ne);

for e = 1:ne
  besi(1:R+1,e) = besseli(0:R, v(1,e))';
  besk(1:R+1,e) = besselk(0:R, v(2,e))';
end

% <<< these are not actually radial MF, the part that       >>>
% <<< cancels when forming the ratio of two functions       >>>
% <<< is left out, since they are only ever used as a ratio.>>>
for n = 1:N
  Ke(n,1:ne) = sum(AEv(1:R,n*ones(1,ne)).*besi(1:R,1:ne).*besk(1:R,1:ne),1);
  Ko(n,1:ne) = sum(AOd(1:R,n*ones(1,ne)).*(besi(1:R,1:ne).*besk(2:R+1,1:ne) + ...
                                           besi(2:R+1,1:ne).*besk(1:R,1:ne)),1);
end


%% compute Ke' for each value of eta
DKe = zeros(N,ne);   DKo = zeros(N,ne);
Dbesi = zeros(R+1,ne); Dbesk = zeros(R+1,ne);

% use recurrences to compute BF derivs
Dbesi(1,1) = besi(2,1);
Dbesk(1,1) = -besk(2,1);
Dbesi(2:R,1) =  0.5*(besi(1:R-1,1) + besi(3:R+1,1));
Dbesk(2:R,1) = -0.5*(besk(1:R-1,1) + besk(3:R+1,1));
Dbesi(R+1,1) =  besi(R,1) - (R+1)/v(1,1).*besi(R+1,1);
Dbesk(R+1,1) = -besk(R,1) - (R+1)/v(2,1).*besk(R+1,1);

% <<< similarly, these are not actually radial MF derivatives;   >>>
% <<< they are just what does not cancel when forming the ratio  >>>
for n = 1:N
  DKe(n,1) = sqrt(q)*sum(AEv(1:R,n).*...
    (besi(1:R,1).*Dbesk(1:R,1).*expeta(2*ones(R,1),1) - ...
    Dbesi(1:R,1).*besk(1:R,1).* expeta(1*ones(R,1),1)),1);
  DKo(n,1) = sqrt(q)*sum(AOd(1:R,n).*...
    (besi(1:R,1).*Dbesk(2:R+1,1).*expeta(2*ones(R,1),1) - ...
    Dbesi(1:R,1).*besk(2:R+1,1).*expeta(1*ones(R,1),1) + ...
    besi(2:R+1,1).*Dbesk(1:R,1).*expeta(2*ones(R,1),1) - ...
    Dbesi(2:R+1,1).*besk(1:R,1).*expeta(1*ones(R,1),1)),1);
end

%% compute H for all points on grid
besiREv(1:R,1) = besseli(2*(0:R-1),  elB).';
besiROd(1:R,1) = besseli(2*(0:R-1)+1,elB).';

coeffEv(1:N,1) =  2.0*iv(1:N).*sum(AEv(1:R,1:N).*besiREv(1:R,ones(1,N)),1)';
coeffOd(1:N,1) = -2.0*iv(1:N).*sum(AOd(1:R,1:N).*besiROd(1:R,ones(1,N)),1)';

H = zeros(ne,np);
for e=1:ne
  for p=1:np
    H(e,p) = sum(coeffEv(1:N).*ce(1:N,p).*Ke(1:N,e)./Ke(1:N,1) + ...
                 coeffOd(1:N).*se(1:N,p).*Ko(1:N,e)./Ko(1:N,1));
  end
end

% HCheck = zeros(N,np);
% for n=1:N
%   for p=1:np
%     HCheck(n,p) = sum(coeffEv(1:n).*ce(1:n,p) + coeffOd(1:n).*se(1:n,p));
%   end
% end


%% compute corresponding circular soluition for comparison
if ecc <= 0.1
  r(1:ne) = 0.5*F*exp(eta(1:ne)); % asymptotic relationship between eta and R
  HC = zeros(ne,np);
  
  beskEvC = zeros(N,ne); beskOdC = zeros(N,ne);
  
  besiEvRC(1:N,1) = besseli(2*(0:N-1),   r(1)).';
  besiOdRC(1:N,1) = besseli(2*(0:N-1)+1, r(1)).';
  beskEvC0(1:N,1) = besselk(2*(0:N-1),   r(1)).';
  beskOdC0(1:N,1) = besselk(2*(0:N-1)+1, r(1)).';
   
  half = ones(N,1);  half(1) = 0.5;
    
  CCoeffEv(1:N,1) =  2.0*iv(1:N).*besiEvRC(1:N).*half(1:N);
  CCoeffOd(1:N,1) = -2.0*iv(1:N).*besiOdRC(1:N);
  
  % circular Yukawan potential
  for e=1:ne
    beskEvC(1:N,e) = besselk(2*(0:N-1),r(e)).';
    beskOdC(1:N,e) = besselk(2*(0:N-1)+1,r(e)).';
    for p=1:np
      HC(e,p) = sum(...
        CCoeffEv(1:N).*beskEvC(1:N,e)./beskEvC0(1:N).*cos(2*(0:N-1)'*psi(p))) + sum(...
        CCoeffOd(1:N).*beskOdC(1:N,e)./beskOdC0(1:N).*sin((2*(0:N-1)'+1)*psi(p)));
    end
  end
   
%   HCCheck = zeros(2*N,np);
%   for n = 1:N
%      for p = 1:np
%       HCCheck(2*(n-1)+1,p) = sum(CCoeffEv(1:n).*cos(2*(0:n-1)'*psi(p))) ;
%       HCCheck(2*(n-1)+2,p) = sum(CCoeffOd(1:n).*sin((2*(0:n-1)'+1)*psi(p)));
%     end
%   end
  
end

%% compute dimensionless flux along the boundary of ellipse
% as a function of angular position
V0 = zeros(np,1);

for p=1:np
  flux = exp(elB*sin(psi(p)))*(...
      sum(coeffEv(1:N).*ce(1:N,p).*DKe(1:N,1)./Ke(1:N,1) + ...
          coeffOd(1:N).*se(1:N,p).*DKo(1:N,1)./Ko(1:N,1)) + ...
      elA*sin(psi(p))*sum(coeffEv(1:N).*ce(1:N,p) + coeffOd(1:N).*se(1:N,p)));
  V0(p) = 2.0*elA*sin(psi(p)) - flux;
end

%% compute things for plotting

X = F*cosh(eta(ones(np,1),1:ne))'.*cos(psi(ones(ne,1),1:np));
Z = F*sinh(eta(ones(np,1),1:ne))'.*sin(psi(ones(ne,1),1:np));  % z is positive down
phi = 0.5*log(H.*exp(Z));
PHI = phi - Z;
if F <= 0.1
  phiC = 0.5*log(HC.*exp(Z));
  PHIC = phiC - Z;
end

%% plot figures
range = [0,4,-1,6];
hcontours = -10:0.5:5;
pcontours = -10:0.25:5;

% hydraulic head contours
figure
[C,h] = contour(X,Z,PHI,hcontours);
set(gca,'YDir','reverse')
% xlabel('X','FontSize',15)
% ylabel('Z','FontSize',15)
daspect([1 1 1]);
axis(range);
hold on;
colormap(zeros(size(gray)))
plot(elA*cos(psi),elB*sin(psi),'-r','LineWidth',4);
clabel(C,h,'manual','FontSize',15,'Rotation',0)
h(1) = text('Position',[1.9,6.5],'fontsize',15, ...
  'interpreter','latex','string','$X$');
h(2) = text('Position',[-0.5,2.5],'fontsize',15, ...
  'interpreter','latex','string','$Z$');

% pressure head contours
figure
[C,h] = contour(X,Z,phi,pcontours);
set(gca,'YDir','reverse')
% xlabel('X','FontSize',15)
% ylabel('Z','FontSize',15)
daspect([1 1 1]);
axis(range);
hold on;
colormap(zeros(size(gray)))
plot(elA*cos(psi),elB*sin(psi),'-r','LineWidth',4);
clabel(C,h,'manual','FontSize',15,'Rotation',0)
h(1) = text('Position',[1.9,6.5],'fontsize',15, ...
  'interpreter','latex','string','$X$');
h(2) = text('Position',[-0.5,2.5],'fontsize',15, ...
  'interpreter','latex','string','$Z$');

% if ecc <= 0.1
%   figure
%   contourf(X,Z,PHIC,hcontours);
%   set(gca,'YDir','reverse')
%   daspect([1 1 1]);
%   axis(range);
%   title('circular solution PHIC')
%   colorbar;
%   hold on;
%   plot(elA*cos(psi),elB*sin(psi),'-k','LineWidth',3);
%   
%   figure
%   contourf(X,Z,phiC,pcontours);
%   set(gca,'YDir','reverse')
%   daspect([1 1 1]);
%   axis(range);
%   title('circular solution phiC')
%   colorbar;
%   hold on;
%   plot(elA*cos(psi),elB*sin(psi),'-k','LineWidth',3);
  
% % difference between circular and elliptical solutions
%   figure
%   contourf(X,Z,phi-phiC,20);
%   set(gca,'YDir','reverse')
%   daspect([1 1 1]);
%   axis([0,4,-4,6]);
%   xlabel('X')
%   ylabel('-Z')
%   colorbar;
%   hold on;
%   plot(elA*cos(psi),elB*sin(psi),'-k','LineWidth',3);
% end



