%% Analysis of speaker blocked coil impedance vs. position
%% David Meeker
%% 01Nov2015

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

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

% Frequency under consideration in Hz
w = 1000;

% 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',1) call 
opendocument(myFile);
mi_setcurrent(myCoilName,1);
mi_probdef(w);
mi_setprevious('DCProblem.ans',1);
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=[];
Xx=xmin:dx:xmax;

for x=xmin:dx:xmax
	mi_setfocus('DCProblem.fem');
	mi_analyze();
	% At each coil position, evaluate Bl and coil inductance over a range of frequencies 
	
	mi_setfocus('ACProblem.fem');
	mi_analyze();
	mi_loadsolution();
    
    % Evaluate Bl curve
    mo_groupselectblock(1);
	f=abs(mo_blockintegral(29));
    Bl=[Bl, f];
            
	% Get Inductance and Voltage 
	u=mo_getcircuitproperties(myCoilName);
	Lx=[Lx, u(3)];
	Vx=[Vx, u(2)];
	fprintf('%g %s\n',x,num2str(u(3)));
    
    % Move coil to next position
	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);
plot(Xx,abs(Vx));
xlabel('Displacement, mm');
ylabel('Impedance, Ohms');
title('abs(Impedance) vs. Position at 1kHz');

figure(2);
plot(Xx,abs(1000*Lx));
xlabel('Displacement, mm');
ylabel('abs(Inductance), mH');
title('abs(Inductance) vs. Position at 1kHz');

figure(3);
plot(Xx,real(1000*Lx));
xlabel('Displacement, mm');
ylabel('real(Inductance), mH')
title('real(Impedance) vs. Frequency and Position at 1kHz');

figure(4);
plot(Xx,real(2j*w*Lx));
xlabel('Displacement, mm');
ylabel('omega*imag(inductance), Ohms')
title('real(jw*Inductance) vs. Frequency and Position');

figure(5);
plot(Xx,Bl);
xlabel('Displacement, mm');
ylabel('B*l, N/A')
title('B*l vs. Position at 1kHz');