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

% 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; % 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',1) call 
opendocument(myFile);
mi_setcurrent(myCoilName,1);
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=[];
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))];
        
	% 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')