%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LRK Motor Flux Linkage Analysis        %
% (Octave Version)                       %
% dmeeker@ieee.org
% February 18, 2006                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

openfemm

% load geometry
try
  v=ver;
  opendocument([cd,'\lrk.fem']);
catch
  opendocument('lrk.fem');
end

% and save as a temporary file name so we don't 
% ovewrite the original geometry
mi_saveas('temp.fem');

% make sure that the currents are zeroed out
% we just want flux linkage due to the magnets
mi_modifycircprop('A',1,0);
mi_modifycircprop('B',1,0);
mi_modifycircprop('C',1,0);

% move the rotor through one pole pitch in small
% increments, recording the flux linkage of each
% phase at each rotor position
steps = 60;
data = zeros(steps+1,5);
for k = 1:(steps+1)
  disp(sprintf('%i/%i',k-1, steps));
  mi_analyze(1);
  mi_loadsolution;
  data(k,1) = (k-1)*360/14/steps;

  % collect the flux linkage for each phase
  v = mo_getcircuitproperties('A'); data(k,2)=v(3);
  v = mo_getcircuitproperties('B'); data(k,3)=v(3); 
  v = mo_getcircuitproperties('C'); data(k,4)=v(3);

  % compute the cogging torque
  mo_groupselectblock(2);
  data(k,5) = mo_blockintegral(22);

  % increment the rotor's position
  mi_seteditmode('group');
  mi_selectgroup(2);
  mi_moverotate(0, 0, 360/14/steps);
  mi_clearselected();
end

% plot cogging torque results...
plot(data(:,1),data(:,5));

% write results to disk for further analysis.
fp=fopen('flux-results.txt','w');
for k = 1:(steps+1)
	fprintf(fp,'%g	%g	%g	%g	%g\n',data(k,1),data(k,2),data(k,3),data(k,4),data(k,5));
end
fclose(fp);

closefemm