% Plot WKB approximation to bending of an obliquely upward-propagating
% light wave phi(x,z,t) = Re f(z)*exp(i(kx-om*t)) by a layer in
% which refractive index n(z) increases with z.  Local dispersion relation is
% om^2 = c0^2*(k^2 + m^2(z))/n^2(z).

% Assumptions:
% (1) n(z) = n1 + (n2-n1)*(1+erf(z))/2  where z is nondimensional height.
%     Thus n -> n1 as z->-inf, n -> n2 as z->inf.  Here n1 = 1, n2 = 2.
% (2) Nondimensionalized light speed in vacuum c0 = 1
% (3) Incoming wave of unit amplitude is at 45 degrees to the vertical. Hence
%     incoming vertical wavenumber m1 = hor. wavenumber k.
% (4) For WKB to hold. need m1 >>1. In practice we take m1 = 2*pi
%
% Plot is for t = 0

%------User-specified parameters------------------------
  c0 = 1; % Light speed in vacuum
  m1 = 2*pi;% Incoming vertical wavenumber
  k = m1; % Hor. wavenumber
  n1 = 1; % Refractive index below gradient layer.
  n2 = 2; % Refractive index above gradient layer (n2 > n1).
  f1 = 1; % Incoming wave amplitude
  zmin = -2;  % Min, max and increment in z for plots
  zmax = 2;
  dz = 0.05;
  xmin = 0;   % Min, max and increment in x for plots
  xmax = 4.4;
  dx = 0.1;

%------Calculations--------------------------------------
  om = sqrt(c0^2*(k^2 + m1^2)/n1^2); % Wave frequency
  x = xmin:dx:xmax;
  z = zmin:dz:zmax;
  nz = length(z);
  n = n1 + (n2-n1)*(1+erf(z))/2;
  m = sqrt(n.^2*om^2/c0^2 - k^2);
  
  % Find mint = int_zmin^z m(zeta) dzeta at each gridpoint using
  % trapezoidal rule. mint(k) = dz*(0.5*m(1)+m(2)+...+m(k-1)+0.5*m(k))
  
  mint = 0.5*dz*[0 cumsum(m(1:nz-1)+m(2:nz))]; 
  f = f1*sqrt(m1./m).*exp(1i*mint);
  [X F] = meshgrid(x,f);
  PHI = real(exp(1i*k*X).*F);

  % Contour plot
  
  contourf(x,z,PHI)
  colormap gray
  colorbar
  axis equal
  xlabel('x')
  ylabel('z')
  hold on
  plot([xmin xmax],[-1 -1],'r--',[xmin xmax],[1 1],'r--' ,'LineWidth',2)
  title('WKB approximation to refraction through a density gradient')
  hold off