%% Analysis of speaker with shorting ring
%% David Meeker
%% 23Nov2015

% Base problem name and configuration 
myFile='SSF-082.fem';
myCoilName='Icoil';
myCoilGroup=1;

% Displacement range under consideration 
xmin=-15;
xmax=15;
dx=5;

% Frequency range under consideration.
wmin = 0; % 10^0 = 1 Hz 
wmax = 4; % 10^4 = 10000 Hz 
dw = 1; % 5 points per decade on a log scale 

% Open up base problem and save as the DC operating point with no coil current 
openfemm;
opendocument(myFile);
mi_setcurrent(myCoilName,0);
mi_saveas('DCProblem.fem');

% Open up base problem and save as the incremental AC problem with a coil current of 1A 
% Incremental AC problem points back to the DC solution via the mi_setprevious('DCProblem.ans') call 
opendocument(myFile);
mi_setcurrent(myCoilName,1);
mi_setprevious('DCProblem.ans');
mi_saveas('ACProblem.fem');

% move the coil to the xmin position, assuming that the geometry is drawn so that the coil is nominally at x=0 
mi_setfocus('DCProblem.fem');
mi_selectgroup(myCoilGroup);
mi_movetranslate(0,xmin);
mi_setfocus('ACProblem.fem');
mi_selectgroup(myCoilGroup);
mi_movetranslate(0,xmin);

% Step through the full range of coil positions 
Lx=[];
Vx=[];
Bl=[];
Lo=[];
Xx=xmin:dx:xmax;
Ww=10.^(wmin:dw:wmax);

for x=xmin:dx:xmax
	Lfd=[];
	Vfd=[];
	mi_setfocus('DCProblem.fem');
	mi_analyze();
    
    % Get inductance in the limiting case as frequency goes to zero
    mi_setfocus('ACProblem.fem');
	mi_probdef(10^(-6));
	mi_analyze();
	mi_loadsolution();
    u=mo_getcircuitproperties(myCoilName);
	Lo=[Lo, abs(u(3))];
    Rdc=abs(u(2)/u(2));
        
	% At each coil position, evaluate Bl and coil inductance over a range of frequencies 
	for k=wmin:dw:wmax
		mi_setfocus('ACProblem.fem');
		mi_probdef(10^k);
		mi_analyze();
		mi_loadsolution();
		% Get Bl 
		if (k==0)
			mo_groupselectblock(1);
			f=abs(mo_blockintegral(29));
			Bl=[Bl, f];
		end;
		% Get Inductance and Voltage 
		u=mo_getcircuitproperties(myCoilName);
		Lfd=[Lfd, u(3)];
		Vfd=[Vfd, u(2)];
		fprintf('%g	%g %s\n',x,10^k,num2str(u(3)));
	end;
	Lx=[Lx; Lfd];
	Vx=[Vx; Vfd];
	if (x<xmax)
		mi_setfocus('DCProblem.fem');
		mi_selectgroup(myCoilGroup);
		mi_movetranslate(0,dx);
		mi_setfocus('ACProblem.fem');
		mi_selectgroup(myCoilGroup);
		mi_movetranslate(0,dx);
	end;
end;
closefemm;

% Plot up some results
figure(1);
loglog(Ww,abs(Vx));
xlabel('Frequency, Hz');
ylabel('Impedance, Ohms');
title('Impedance vs. Frequency and Position');

figure(2);
loglog(Ww,abs(1000*Lx));
xlabel('Frequency, Hz');
ylabel('Inductance, mH');
title('Inductance vs. Frequency and Position');

figure(3);
plot(Xx,Bl);
xlabel('Frequency, Hz')
ylabel('B*l, N/A, Ohms')
title('B*l vs. Frequency and Position')

% Compute network parameters
Ls=zeros(7,5);
Rs=zeros(7,4);
for k=1:7
    parms=GetParms(Ww,Lx(k,:),Lo(k));
    disp(parms);
    Ls(k,1:5)=parms(1:5);
    Rs(k,1:4)=parms(6:9);
end

disp('Inductance values for each position step:');
disp(Ls); disp(' ');
disp('Resistance values for each position step:');
disp(Rs); disp(' ');

% Regression of coefficients needed to describe inductances vs. position
Xmat=[ones(length(Xx),1), Xx', (Xx.^2)', (Xx.^3)'];
Lfit=zeros(5,4);
for k=1:5
    Lfit(k,:)=((Xmat'*Xmat)\(Xmat'*log(Ls(:,k))))';
end
disp('Coefficients needed to describe inductances vs. position:');
disp(Lfit); disp(' ');

% Regression of coefficients needed to describe resistances vs. position
Rfit=zeros(4,4);
for k=1:4
    Rfit(k,:)=((Xmat'*Xmat)\(Xmat'*log(Rs(:,k))))';
end
disp('Coefficients needed to describe resistances vs. position:');
disp(Rfit); disp(' ');

% Regression of coefficients needed to describe Bl curve vs. position
Blfit = ((Xmat'*Xmat)\(Xmat'*Bl'))';
disp('Coefficients needed to describe Bl curve vs. position:');
disp(Blfit);

% Save parameters for use in transient Simulink model
save('LoudspeakerParameters.mat','Lfit','Rfit','Blfit','Rdc');