%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Analysis of a Woofer Motor           %%
%%                                      %%   
%% David Meeker                         %% 
%% dmeeker@ieee.org                     %%
%%                                      %%
%% 15 May2009                           %%
%%                                      %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Introduction
%
% The purpose of this script is to automate the analysis of a woofer motor.
% The notebook produces the static performance curves that designers often
% desire during the development of a new transducer: the actuator's "BL" 
% (the amount of force produced on the voice coil per amp of coil current); 
% the coil's self-inductance, and a plot of the magnetic field in which the 
% coil is immersed.
%
% Because there are so many possible variations in speaker motor geometry,
% the notebook does not create the speaker geometry within the notebook.  Rather, 
% it reads an existing geometry info FEMM and analyses it at a number of different
% coil locations. It is assumed that the speaker motor has been drawn with all
% elements in the coil belonging to group number 1, so that the voice coil can
% be easily selected and moved by the notebook.  It is also assumed that the .fem
% file describing the motor is located in the same directory as the notebook 
% (although the user could manually change the path to the .fem file).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Design - Specific Parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Model Name
ModelName = 'Woofer.fem';

% Maximum excursion + /- from the centered position:
Xlim = 5;

% Movement increments used during the analysis
dX = 1;

% Define the range over which to plot the flux density :
Xcenter = 17.15;    % vertical location of the center of the airgap 
Xspan = 10;         % field plotted between Xcenter-Xspan and Xcenter+Xspan

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Analysis Routines
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Open an instance of FEMM
openfemm;

% Analyze BL and incremental inductance at 1 mm steps between - Xlim and + Xlim
opendocument(ModelName);
mi_saveas('temp.fem');
mi_selectgroup(1);
mi_movetranslate(0, -Xlim);
mi_clearselected;

bl = [];
inductance = [];
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('% Speaker Motor Analysis Script');
disp('%');
disp('% David Meeker');
disp('% dmeeker@ieee.org');
disp('% 15May2009');
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');
disp('');
disp(sprintf('Disp(mm)\t\tBL(N/A)\t\tInductance(uH)'));
for k=-Xlim:dX:Xlim,
    mi_modifycircprop('icoil',1,1);
    mi_analyze(1);
    mi_loadsolution();
    mo_groupselectblock(1);
    fz = mo_blockintegral(12);
    bl = [bl; [k, fz]];
    R   = mo_getcircuitproperties('icoil')*[0;1;0];
    fl1 = mo_getcircuitproperties('icoil')*[0;0;1];
    mo_close;
    mi_modifycircprop('icoil', 1, 0);
    mi_analyze(1);
    mi_loadsolution;
    fl0 = mo_getcircuitproperties('icoil')*[0;0;1];
    L = (fl1 - fl0)*10^6;
    inductance = [inductance; [k, L]];
    disp(sprintf('%g\t\t%g\t\t%g',k, fz, L));
    mi_selectgroup(1);
    mi_movetranslate(0, 1);
end

mi_close;

% Evaluate flux density in the gap at the position of the coil
opendocument(ModelName)
mi_saveas('temp.fem');
mi_modifycircprop('icoil', 1, 0);
mi_analyze(1);
mi_loadsolution;

fluxdensity = [];
x=-Xspan;
for x=-Xspan:(Xspan/100):Xspan 
  b=mo_getb(16,Xcenter+x);
  fluxdensity=[fluxdensity;[x,abs(b(1))]];
end

closefemm;

delete('temp.fem');
delete('temp.ans');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Plot Analysis Results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1);
plot(bl(:,1),bl(:,2));
xlabel('Displacement, mm');
ylabel('BL, N/A');
title('BL vs. Displacement');

figure(2);
plot(inductance(:,1),inductance(:,2));
xlabel('Displacement, mm');
ylabel('Inductnace, uH');
title('Coil DC Self Inductance vs. Displacement');

figure(3);
plot(fluxdensity(:,1),fluxdensity(:,2));
xlabel('Vertical distance from center of gap, mm');
ylabel('Flux density, T');
title('Air Gap Radial Flux Density vs. Position');

% DC Coil Resistance in Ohms:
disp(sprintf('\nDC Coil Resistance = %g',R));

% Compute the nominal range of operation of the speaker.  It is assumed that the 
% operating range is the range in which BL is within 10% of the maximum value of BL

BLMax = max(bl(:,2));
BLMin = 0.9*BLMax;

for k=2:length(bl)
    if (bl(k,2)>=BLMin) 
        break;
    end      
end

Xmin = bl(k-1,1)+(bl(k,1)-bl(k-1,1))*(BLMin-bl(k-1,2))/(bl(k,2)-bl(k-1,2));

kk=k;
for k=kk:length(bl)
    if (bl(k,2)<=BLMin) 
        break;
    end      
end

Xmax = bl(k-1,1)+(bl(k,1)-bl(k-1,1))*(BLMin-bl(k-1,2))/(bl(k,2)-bl(k-1,2));

disp(sprintf('Xmin = %g',Xmin));
disp(sprintf('Xmax = %g',Xmax));