%
%   PROGRAM CALCULATES VALUES OF THE RLC PARAMETERS
%               FOR 
%   AN ULTRASONIC TRANSDUCER 
%   given transducer's frequency response
%   (reflection coefficients vs. frequency)
%

clear

%
%   WRITE NAME OF YOUR FILE HERE
%

filnam= input('WRITE NAME OF YOUR FILE      ','s')

%	Read data from text file

eval(['load ', filnam '.txt']);
Fdat = eval([filnam '(:,10) + j*' filnam '(:,11)']);  
freq = eval([filnam '(:,1)']); 

%   Convert reflection coefficients to admitance and impedance

Yin= (1 - Fdat)./(1 + Fdat)/50;
Zin = (1 + Fdat)./(1 - Fdat)*50;

% Parameter estimation
 
%Set number of points on either side of peak to find fit
N=5;
%Find the maximum in real(Yin)
[m,n]=max(real(Yin));

%Determine polynomial coefficients 
P=polyfit(freq(n-N:n+N),(Yin(n-N:n+N)),5);

%Extract region over which extrapolation will be done.
f=freq((n-N):(n+N));
y=polyval(P,f);



% Plot resonance approximation


figure(1)
clf
subplot(211)
plot(freq*1e-3,real(Yin))
hold
plot(f*1e-3,real(y),'r-.')
xlabel('Frequency (kHz)')
ylabel('Real(Yin)')
title( ['TRANSDUCER ' filnam ' - Estimated resonance'])

subplot(212)
plot(freq*1e-3,imag(Yin))
hold
plot(f*1e-3,imag(y),'r-.')

xlabel('Frequency (kHz)')
ylabel('Imag(Yin)')



%Determine Resonance Frequency
[m2,n2]=max(real(y));

f0e = f(n2);
w0 = 2*pi*f0e;

%Determine R1
R1=1/m2;
%Determine Co
Co=imag(y(n2))/(2*pi*f0e)


%Find X1
w= f*2*pi;
X1=imag(1./(y-j*w*Co));

%Determine C1
C1=0;
cnt=0;

for wi=1:length(X1)
if abs(X1(wi))>1e-9
val=((w(wi)/w0)^2 - 1)/(w(wi)*X1(wi));
C1=C1+val;
cnt=cnt+1;
end
end

C1=C1/cnt;

%Determine L1
L1=1/(C1*w0^2);

%Find admitance for RLC model

for f1i=1:length(freq)
   
   
   w=2*pi*freq(f1i);
   Zp1=R1 + j*w*L1 + 1/(j*w*C1);
   Yine(f1i)=(j*w*Co + 1/Zp1);
   Zine(f1i) =  1/Yine(f1i);
   
end

% Print estimated values

fprintf('\n Estimated RLC parameters');
fprintf('\n  fo = %0.3f kHz;',f0e/1000);
fprintf('\n  R1 = %0.2e ohm;',R1);
fprintf('\n  C0 = %0.4e nF;',Co*1e9);
fprintf('\n  C1 = %0.4e nF;',C1*1e9);
fprintf('\n  L1 = %0.4e mH;',L1*1e3);

fprintf('\n\n\n');

%   Plot admitance

figure(2)
clf
subplot(211)
plot(freq*1e-3,real(Yin)*1e3)
hold
plot(freq*1e-3,real(Yine)*1e3,'r')

legend('measured','estimated')
xlabel('Frequency (kHz)')
ylabel('Real(Y_i_n) [mS]')
title( ['TRANSDUCER ' filnam ' - Estimated admitance'])

subplot(212)
plot(freq*1e-3,imag(Yin)*1e3)
hold
plot(freq*1e-3,imag(Yine)*1e3,'r')

%grid
xlabel('Frequency (kHz)')
ylabel('Imag(Y_i_n)[mS]')

%   PLOT impedance

figure(3)
clf

subplot(211)

plot(freq*1e-3,real(Zin)/1e3)
hold
plot(freq*1e-3,real(Zine)/1e3,'r')
legend('measured','estimated')
xlabel('Frequency (kHz)')
ylabel('Real(Z_i_n) [k\Omega]')
title( ['TRANSDUCER ' filnam ' - Estimated impedance'])


subplot(212)

plot(freq*1e-3,imag(Zin)/1e3)
hold
plot(freq*1e-3,imag(Zine)/1e3,'r')
%grid
xlabel('Frequency (kHz)')
ylabel('Imag(Z_i_n) [k\Omega]')