%This MATLAB code can be used to reproduce the data in Fig. 3-11
%	Courtesy of LARS ANDERSEN
%t=thickness of slab
%k0=free space propagation constant
%xa=location of the left computational domain endpoint
%alphaa=alpha coeff. to be used in the boundary condition at xa; see (3.43)
%betaa=betaa coeff. to be used in the boundary condition at xa; see (3.43)
%p=p coefficient appearing in the differential equation (3.2)
%epr=relative permittivity of the slab
%beta=as defined in Fig. 3-11 
%q=coeffient appearing in the differential equation (3.2)
%N=number of nodes (N-1=number of layers)
%R_FEM=computed reflection coeffient (to be plotted 

% Initialization 

clear;
%diary file1;
N=7;
t=0.25;
Dx=t/(N-2);
xa=t+Dx;
p=-1;
k0=2*pi;
f=0;
alphaa=j*k0;
betaa=2*j*k0*exp(j*k0*xa);
%values to generate plot in Fig. 3-11
Z=26;
%epsr=4-j*beta;
for z=1:Z
	beta=(z-1)/(Z-1)*5;
	epsr=4-j*beta

% Initialize global matrix and vector

for m=1:N,
  b(m,1)=0;
  for n=1:N,
    A(m,n)=0;
  end
end

% Compute element matrices and assemble global matrix

for n=1:N-1,
  xe1=(n-1)*xa/(N-1);
  xe2=n*xa/(N-1);
  if n==N-1,
    eps=1;
  else
    eps=epsr;
  end
  q=eps*k0^2;
  Ael(1,1)=p/abs(xe2-xe1)+q*abs(xe2-xe1)/3;
  Ael(2,2)=Ael(1,1);
  Ael(1,2)=-p/abs(xe2-xe1)+q*abs(xe2-xe1)/6;
  Ael(2,1)=Ael(1,2);
  bel(1)=f*abs(xe2-xe1)/2;
  bel(2)=bel(1);
  A(n,n)=A(n,n)+Ael(1,1);
  A(n,n+1)=A(n,n+1)+Ael(1,2);
  A(n+1,n)=A(n+1,n)+Ael(2,1);
  A(n+1,n+1)=A(n+1,n+1)+Ael(2,2);
  b(n,1)=b(n,1)+bel(1);
  b(n+1,1)=b(n+1,1)+bel(2);
end

% Enforce boundary conditions

A(1,1)=1;
for n=2:N,
  A(n,1)=0;
  A(1,n)=0;
end
A(N,N)=A(N,N)+alphaa*p;
b(N)=b(N)+betaa*p;

% Solve the resulting NxN system

x=inv(A)*b;

% Compute reflection coefficient

R_FEM=abs((x(N)-exp(j*k0*xa))/exp(j*k0*xa));
zeta0=377;
zeta1=sqrt(1/epsr)*zeta0;
k1=sqrt(epsr)*k0;
R_exact=abs((zeta1*tanh(j*k1*t)-zeta0)/(zeta1*tanh(j*k1*t)+zeta0));
Rf(z,1)=beta;
Rf(z,2)=R_FEM;
Re(z,1)=beta;
Re(z,2)=R_exact;

end

%Plot Reflection Coefficient

clf;
plot(Rf(:,1),Rf(:,2),'r');
hold;
plot(Re(:,1),Re(:,2));
xlabel('beta');
ylabel('|R|');

% End of program