% The following code generates a simplified model of the balanced armature
% receiver (small hearing aid speaker) from Kim et al. 2013. It does not
% include the semi-inductor or transmission line

NFT=1024; NF=1+NFT/2;%number of non-negative frequencies
Fmax=1e4; Fs=2*Fmax; Fmin=Fs/NFT;
f=0:Fmin:Fmax; %sweep frequency from 0 to Fmax [Hz] in steps of Fmin [Hz]
f = f(10:end); % we don't need the lowest frequencies

% define load condition
A_canal = pi*0.0075^2; % ear canal area m^2

% useful physical constants
rho = 1.2; % density of air kg/m^3
c = 342; % speed of sound m/s
P0 = 1E5; % atmospheric pressure
eta = 1.4; % ratio of specific heats

%Define basic parameters based on Kim and Allen
Re=195; 
Le=9e-3; %electrical section
Lem = 52e-3;
LeTot = Le+Lem;
G = 7.5; %Gryator Bo l: Electrical to Mechanical transformation
Cm=1.25e-3; Lm=4.3e-6; Rm=3e-3; %Mechanical section
A = 2.4e-6; %Mechanical to Acoustical transformation (earphone port radius ~ 1 mm)
Ca = 4.3E-15; % acoustic volume
z0tx = 1e9; ltx = 1e-4; % Kim et al 2013 transmission line parameters
Vtx = ltx*(rho*c/z0tx); % volume of transmission line tube
Ctx = Vtx/(eta*P0); % capacitance of transmission line volume
CaTot = Ca+Ctx;

for k=1:length(f); 
    s=2*pi*f(k)*1j; %Define complex frequency
      
    Te = [1 Re+s*LeTot;0 1]; %electrical
    Tem = [0 G; 1/G 0]; % gyrator electrical to mechanics
    Tm = [1 1./(s*Cm); 0 1]*[1 s*Lm; 0 1]*[1 Rm; 0 1]; % mechanical
    Tma = [A 0; 0 1/A]; % coupling of mechanics to acoustics
    Ta = [1 0; s*CaTot 1]; % acoustical

    % indices: 0= front end; 1 = left side of mechanical part (right of
    % gyrator); 2 = right side of mechanical part (left of transformer); 3
    % = left side of acoustical part (right side of transformer)
    T03(k,:,:)=Te*Tem*Tm*Tma*Ta; %avoids a nightmare of algebra
    T01(k,:,:)=Te*Tem;
    T02(k,:,:)=Te*Tem*Tm;
    T13(k,:,:)=Tm*Tma*Ta;
    T23(k,:,:)=Tma*Ta;
    
    % calculate Z matrices here to make calculations easier later
    Z03(k,1,1) = T03(k,1,1)/T03(k,2,1); % A/C
    Z03(k,1,2) = -1/T03(k,2,1); % det(T)/C where det(T) = -1 for anti-reciprocal system (gyrator)
    Z03(k,2,1) = 1/T03(k,2,1); % 1/C
    Z03(k,2,2) = T03(k,2,2)/T03(k,2,1); % D/C
    
    Z02(k,1,1) = T02(k,1,1)/T02(k,2,1); % A/C
    Z02(k,1,2) = -1/T02(k,2,1); % det(T)/C where det(T) = -1 for anti-reciprocal system (gyrator)
    Z02(k,2,1) = 1/T02(k,2,1); % 1/C
    Z02(k,2,2) = T02(k,2,2)/T02(k,2,1); % D/C
    
    Z13(k,1,1) = T13(k,1,1)/T13(k,2,1); % A/C
    Z13(k,1,2) = 1/T13(k,2,1); % det(T)/C where det(T) = 1 for recriprocal system
    Z13(k,2,1) = 1/T13(k,2,1); % 1/C
    Z13(k,2,2) = T13(k,2,2)/T13(k,2,1); % D/C
    
    Z23(k,1,1) = T23(k,1,1)/T23(k,2,1); % A/C
    Z23(k,1,2) = 1/T23(k,2,1); % det(T)/C where det(T) = 1 for reciprocal system
    Z23(k,2,1) = 1/T23(k,2,1); % 1/C
    Z23(k,2,2) = T23(k,2,2)/T23(k,2,1); % D/C
   
end %end for freq loop

% Solve for transfer functions: 
I0divV0_U10 = T01(:,2,1)./T01(:,1,1);
% TO DO: Solve for other transfer functions

% Solve for input impedance looking to the right of solution point. Use 
% ABCD matrix to model transmission between input (V=1) and solution point.
ZL = rho*c/A_canal; % ear canal
Zin2 = Z23(:,1,1) - Z23(:,1,2).*Z23(:,2,1)./(ZL + Z23(:,2,2));
Zin1 = Z13(:,1,1) - Z13(:,1,2).*Z13(:,2,1)./(ZL + Z13(:,2,2));
Zin0 = Z03(:,1,1) - Z03(:,1,2).*Z03(:,2,1)./(ZL + Z03(:,2,2));

% forces
V0 = ones(1,length(f)); % input 1 V
% TO DO: Calculate F1, F2, P3

% flows
I0 = 1./Zin0;
% TO DO: Calculate U1, U2, V3

% PLOT INPUT IMPEDANCE
figure('units','normalized','outerposition',[.05 .15 .7,.7])
subplot(121)
plot(f/1000, abs(Zin0)); hold on 
set(gca,'xscale','log','yscale','log','xlim',[.2 10]); 
title('|Z_{in}|')
xlabel('Frequency [kHz]')
subplot(122)
plot(f/1000, unwrap(angle(Zin0))/pi); hold on
set(gca,'xscale','log','xlim',[.2 10],'ylim',[-.6 .6]); 
title('\angle Z_{in} [rad/pi]')
xlabel('Frequency [kHz]')
% TO DO: Calculate Zin with other load conditions, compare

% PLOT TRANSFER FUNCTIONS
figure('units','normalized','outerposition',[.05 .15 .7,.7])
subplot(231) %(a)
plot(f/1000, abs(I0divV0_U10)); 
set(gca,'xscale','log','xlim',[.2 10]); title('(a) |I_0/V_0|_{U_1=0}')
xlabel('Frequency [kHz]')
% TO DO: add subplots for (b)-(f)

% PLOT FORCES AND FLOWS
[val,i1k] = min(abs(f-1000)); % find 1kHz index for normalizing curves
figure('units','normalized','outerposition',[.05 .15 .7,.7])
subplot(121)
plot(f/1000, abs(V0)/abs(V0(i1k)),'k','linewidth',2); hold on
legend('|V_0|')
xlabel('Frequency [kHz]'); set(gca,'xscale','log','xlim',[.2 10]); title('Generalized Forces')
% TO DO: Add F1, F2, P3 
subplot(122)
plot(f/1000, abs(I0)/abs(I0(i1k)),'k','linewidth',2); hold on
legend('|I_0|')
xlabel('Frequency [kHz]'); set(gca,'xscale','log','xlim',[.2 10]); title('Generalized Flows')
% TO DO: Add U1, U2, V3