% Calculates the relative permittivity of a PCB substrate
% using lab measurements of the resonances of the ring resonator
% with known geometry.

clear all;
close all;

c0 = 2.998e8;
B = 1.51;          % Substrate thickness in mm
T = 0.0175/B;   % Ratio of conductor thickness to substrate thicness
U = 1.5/B;      % Ratio of trace width to substrate thickness
radius = 293;   % Radius of ring
% Frequencies of measured minimum in S21:
freq = [0.8877867, 1.7878, 2.6753, 3.5627, 4.4377, 5.3127, 6.1877, 7.0377, 7.9127, 8.7752, 9.6252, 10.4876, 11.2876, 12.1251, 12.9626];

U1 = U + T.*log(1 + 4.*exp(1)./T.*(tanh(sqrt(6.51.*U))).^2)./pi;
Ur = @(er) U + (U1 - U).*(1 + 1./(cosh(sqrt(er-1))))./2;

Au = @(er) 1 + log((Ur(er).^4 + Ur(er).^2./2704)./(Ur(er).^4 + 0.432))./49 ...
    + log((Ur(er)./18.1).^3 + 1)./18.7;

Ber = @(er) 0.564.*((er-0.9)./(er+3)).^0.053;

Y = @(er) (er + 1)./2 + (er - 1)./2.*(1+10./Ur(er)).^(-Au(er).*Ber(er));

Z01 = @(x) 376.73.*log((6+(2.*pi-6).*exp(-(30.666./x).^0.7528))./x ...
    + sqrt(4./x.^2 + 1))./(2.*pi);

Z0 = @(er) Z01(Ur(er))./sqrt(Y(er,Ur(er)));

e_eff0 = @(er) Y(er).*(Z01(U1)./Z01(Ur(er))).^2; % Effective perm. at DC

P1 = @(F) 0.27488 + U.*(0.6315+0.525*(0.0157.*F.*B+1).^(-20)) - 0.065683.*exp(-8.7513.*U);
P2 = @(er) 0.33622.*(1-exp(-0.03442.*er));
P3 = @(F) 0.0363.*exp(-4.6.*U).*(1-exp(-(F.*B/38.7).^4.97));
P4 = @(er) 2.751.*(1-exp(-(er./15.916).^8)) + 1;
P = @(er, F) P1(F).*P2(er).*(F.*B.*(0.1844+P3(F).*P4(er))).^1.5763;

e_eff = @(er, F) er - (er-e_eff0(er))./(1+P(er, F));    % Effective perm. at frequency F
Z0F = @(er, F) Z0(er).*sqrt(e_eff0(er)./e_eff(er,F)).*(e_eff(er,F)-1)./(e_eff0(er)-1); % Characteristic impedance at F

v = @(F,rad,k) 2.*pi.*rad.*F./k;
e_eff_ring = @(F, rad, k) (c0./v(F.*1e9,rad.*1e-4,k)).^2;   % Effective perm. of resonator

% Calculate relative permittivity using meausured resonant frequencies
er_calc = @(F,k) fzero(@(er) e_eff(er,F) - e_eff_ring(F,radius,k),4.4);
for i=1:length(freq)
    er_calc2(i) = er_calc(freq(i), i);
end

figure();
hold on;
xlabel('Frequency (GHz)');
ylabel('\epsilon_r');
title(sprintf('B = %f mm', B))
plot(freq, er_calc2,'o');
grid on;
