%% Modified Halbach array by Jose B. Almeida, November 2017

%-----------------------------------------------------------------------------------------------
%% modify these values to suit your needs
units = 'millimeters'	% should be one of "inches","millimeters","centimeters","mils","meters","micrometers" 
forceUnits = 'Newtons'		% should be one of "lbf", "Newtons", "kgf"
magnetLength = 20.0
magnetThick = 5.0
magnetSep = 0.75  % one half of separation between magnet layers
magnetType= 'N52'  % has to be the EXACT name of a material in the materials library
steelThick = 1  % thickness of magnetic circuit steel
interMagnet = 0.7  % thickness of steel inter magnet pieces.
steelType = '1020 Steel'  % has to be the EXACT name of a material in the materials library
%-----------------------------------------------------------------------------------------------

addpath('C:\femm42\mfiles')
openfemm		% opens an OctaveFEMM session


newdocument(0)		% the 0 specifies a magnetics problem
mi_hidegrid()


%% define the problem type, equivalent to the top level menu    Problem
mi_probdef(0, units, "planar", 1E-8, magnetLength)

%% adds these materials from the Material Library to the project
mi_getmaterial(magnetType)	
mi_getmaterial(steelType)
mi_getmaterial('Air')


%% this bit is a bit of a pain. Although we have set the problem units above,
%% that has no effect on the calculation of c0 for the Asymptotic Boundary Condition
%% We therefore need to create a scaling factor for c0, dependant on the units
c0_scale=1.0;
if units=='micrometers' c0_scale= 1000000.0; 
elseif units=='millimeters' c0_scale= 1000.0; 
elseif units=='centimeters' c0_scale= 100.0;
elseif units=='meters'      c0_scale= 1.0;	
elseif units=='inches'      c0_scale= 1.0/0.0254; 
elseif units=='mils'        c0_scale= 1000.0/0.0254 end;


forceScale = 1.0;
if forceUnits=='Newtons' forceScale= 1.0;	%% the natural units of the weighted stress tensor
elseif forceUnits=='lbf'     forceScale= 0.2248089;
elseif forceUnits=='kgf'     forceScale= 0.1019716	end;


%% define the boundary
uo = 4 * pi * 10-7;  %% value of magnetic constant
mi_zoom(-2*magnetThick,-1.5*magnetThick,2*magnetThick,1.5*magnetThick)	  %% set the window to a nice size for the problem
c0= c0_scale/(2*uo*magnetThick);
mi_addboundprop("Asymptotic",0,0,0,0,0,0,c0,0,2)		%% create the Asymptotic Boundary Condition
mi_addboundprop("Periodic",0,0,0,0,0,0,0,0,4)		%% create the Periodic Boundary Condition
mi_addboundprop("Antiperiodic",0,0,0,0,0,0,0,0,5)		%% create the anti-Periodic Boundary Condition


mi_addnode(-1.5*magnetThick-2*interMagnet-0.1,2*magnetThick)
mi_addnode(-1.5*magnetThick-2*interMagnet-0.1,-2*magnetThick)
mi_addnode(2.5*magnetThick+2*interMagnet+0.1,2*magnetThick)
mi_addnode(2.5*magnetThick+2*interMagnet+0.1,-2*magnetThick)
mi_drawline(-1.5*magnetThick-2*interMagnet-0.1,2*magnetThick,-1.5*magnetThick-2*interMagnet-0.1,-2*magnetThick)
mi_drawline(2.5*magnetThick+2*interMagnet+0.1,2*magnetThick,2.5*magnetThick+2*interMagnet+0.1,-2*magnetThick)
mi_drawline(-1.5*magnetThick-2*interMagnet-0.1,2*magnetThick,2.5*magnetThick+2*interMagnet+0.1,2*magnetThick)
mi_drawline(-1.5*magnetThick-2*interMagnet-0.1,-2*magnetThick,2.5*magnetThick+2*interMagnet+0.1,-2*magnetThick)

%% set boundary properties

mi_selectsegment(0,2*magnetThick);
mi_selectsegment(0,-2*magnetThick);
mi_setsegmentprop("Asymptotic",0,0,0,0);
mi_clearselected()
mi_selectsegment(-1.5*magnetThick-2*interMagnet-0.1,0);
mi_selectsegment(2.5*magnetThick+2*interMagnet,0);
mi_setsegmentprop("Periodic",0,0,0,0);
mi_clearselected

%% draw magnets

function drawMagnet (centerx,centery,magnetThick)
  mi_addnode(centerx-magnetThick/2,centery+magnetThick/2)
  mi_addnode(centerx+magnetThick/2,centery+magnetThick/2)
  mi_addnode(centerx-magnetThick/2,centery-magnetThick/2)
  mi_addnode(centerx+magnetThick/2,centery-magnetThick/2)
  mi_drawline(centerx-magnetThick/2,centery+magnetThick/2,centerx+magnetThick/2,centery+magnetThick/2)
  mi_drawline(centerx+magnetThick/2,centery+magnetThick/2,centerx+magnetThick/2,centery-magnetThick/2)
  mi_drawline(centerx+magnetThick/2,centery-magnetThick/2,centerx-magnetThick/2,centery-magnetThick/2)
  mi_drawline(centerx-magnetThick/2,centery-magnetThick/2,centerx-magnetThick/2,centery+magnetThick/2)
endfunction

drawMagnet(-magnetThick-interMagnet,magnetThick/2+magnetSep,magnetThick)
drawMagnet(0,magnetThick/2+magnetSep,magnetThick)
drawMagnet(magnetThick+interMagnet,magnetThick/2+magnetSep,magnetThick)
drawMagnet(2*magnetThick+2*interMagnet,magnetThick/2+magnetSep,magnetThick)
mi_addnode(-1.5*magnetThick - 2*interMagnet,magnetSep + magnetThick + steelThick)
mi_addnode(-1.5*magnetThick - 2*interMagnet,magnetSep)
mi_drawline(-1.5*magnetThick - 2*interMagnet,magnetSep + magnetThick + steelThick,-1.5*magnetThick - 2*interMagnet,magnetSep)


drawMagnet(-magnetThick-interMagnet,-magnetThick/2-magnetSep,magnetThick)
drawMagnet(0,-magnetThick/2-magnetSep,magnetThick)
drawMagnet(magnetThick+interMagnet,-magnetThick/2-magnetSep,magnetThick)
drawMagnet(2*magnetThick+2*interMagnet,-magnetThick/2-magnetSep,magnetThick)
mi_addnode(-1.5*magnetThick - 2*interMagnet,-magnetSep - magnetThick - steelThick)
mi_addnode(-1.5*magnetThick - 2*interMagnet,-magnetSep)
mi_drawline(-1.5*magnetThick - 2*interMagnet,-magnetSep - magnetThick - steelThick,-1.5*magnetThick - 2*interMagnet,-magnetSep)

% now close intermagnet spaces
mi_drawline(-1.5*magnetThick - 2*interMagnet,magnetSep + magnetThick,
  2.5*magnetThick+2*interMagnet,magnetThick+magnetSep)
mi_drawline(-1.5*magnetThick-2*interMagnet,magnetSep,2.5*magnetThick+2*interMagnet,magnetSep)

mi_drawline(-1.5*magnetThick-2*interMagnet,-magnetThick-magnetSep,
  2.5*magnetThick+2*interMagnet,-magnetThick-magnetSep)
mi_drawline(-1.5*magnetThick-2*interMagnet,-magnetSep,2.5*magnetThick+2*interMagnet,-magnetSep)

% draw steel sheet

mi_drawline(-1.5*magnetThick-2*interMagnet,magnetSep+magnetThick+steelThick,
2.5*magnetThick+2*interMagnet,magnetSep+magnetThick+steelThick)
  
  
mi_drawline(2.5*magnetThick+2*interMagnet,magnetSep+magnetThick+steelThick,
 2.5*magnetThick+2*interMagnet,magnetSep+magnetThick)
mi_drawline(-1.5*magnetThick-2*interMagnet,-magnetSep-magnetThick-steelThick,
2.5*magnetThick+2*interMagnet,-magnetSep-magnetThick-steelThick)  
    
mi_drawline(2.5*magnetThick+2*interMagnet,-magnetSep-magnetThick-steelThick,
  2.5*magnetThick+2*interMagnet,-magnetSep-magnetThick)

%% set magnetic properties
mi_clearselected
mi_addblocklabel(2*magnetThick+2*magnetSep,magnetSep+magnetThick/2);
mi_selectlabel(2*magnetThick+2*magnetSep,magnetSep+magnetThick/2);
mi_setblockprop(magnetType, 0, 0, "", 90, 0, 0);
mi_clearselected()
mi_addblocklabel(-magnetThick-interMagnet,magnetSep+magnetThick/2);
mi_selectlabel(-magnetThick-interMagnet,magnetSep+magnetThick/2);
mi_setblockprop(magnetType, 0, 0, "", 0, 0, 0);
mi_clearselected
mi_addblocklabel(magnetThick+interMagnet,magnetSep+magnetThick/2);
mi_selectlabel(magnetThick+interMagnet,magnetSep+magnetThick/2);
mi_setblockprop(magnetType, 0, 0, "", 180, 0, 0);
mi_clearselected
mi_addblocklabel(-0.5,magnetSep+magnetThick/2);
mi_selectlabel(-0.5,magnetSep+magnetThick/2);
mi_setblockprop(magnetType, 0, 0, "", -90, 0, 0);
mi_clearselected()

mi_addblocklabel(2*magnetThick+2*magnetSep,-magnetSep-magnetThick/2);
mi_selectlabel(2*magnetThick+2*magnetSep,-magnetSep-magnetThick/2);
mi_setblockprop(magnetType, 0, 0, "", 90, 0, 0);
mi_clearselected()
mi_addblocklabel(-magnetThick-interMagnet,-magnetSep-magnetThick/2);
mi_selectlabel(-magnetThick-interMagnet,-magnetSep-magnetThick/2);
mi_setblockprop(magnetType, 0, 0, "", 180, 0, 0);
mi_clearselected
mi_addblocklabel(magnetThick+interMagnet,-magnetSep-magnetThick/2);
mi_selectlabel(magnetThick+interMagnet,-magnetSep-magnetThick/2);
mi_setblockprop(magnetType, 0, 0, "", 0, 0, 0);
mi_clearselected
mi_addblocklabel(0,-magnetSep-magnetThick/2);
mi_selectlabel(0,-magnetSep-magnetThick/2);
mi_setblockprop(magnetType, 0, 0, "", -90, 0, 0);
mi_clearselected()

% set block properties

mi_addblocklabel(0,magnetSep/2 + 1.8*magnetThick);
mi_selectlabel(0,magnetSep/2 + 1.8*magnetThick);
mi_setblockprop("Air", 1, 0, "", 0, 1, 0);
mi_clearselected();

if interMagnet > 0
  mi_addblocklabel(-1.5*magnetThick-1.5*interMagnet,magnetSep+magnetThick/2);
  mi_addblocklabel(-0.5*magnetThick-0.5*interMagnet,magnetSep+magnetThick/2);
  mi_addblocklabel(1.5*magnetThick+1.5*interMagnet,magnetSep+magnetThick/2);
  mi_addblocklabel(0.5*magnetThick+0.5*interMagnet,magnetSep+magnetThick/2);
  mi_selectlabel(-1.5*magnetThick-1.5*interMagnet,magnetSep+magnetThick/2);
  mi_selectlabel(-0.5*magnetThick-0.5*interMagnet,magnetSep+magnetThick/2);
  mi_selectlabel(1.5*magnetThick+1.5*interMagnet,magnetSep+magnetThick/2);
  mi_selectlabel(0.5*magnetThick+0.5*interMagnet,magnetSep+magnetThick/2);
  
  mi_addblocklabel(-1.5*magnetThick-1.5*interMagnet,-magnetSep-magnetThick/2);
  mi_addblocklabel(-0.5*magnetThick-0.5*interMagnet,-magnetSep-magnetThick/2);
  mi_addblocklabel(1.5*magnetThick+1.5*interMagnet,-magnetSep-magnetThick/2);
  mi_addblocklabel(0.5*magnetThick+0.5*interMagnet,-magnetSep-magnetThick/2);
  mi_selectlabel(-1.5*magnetThick-1.5*interMagnet,-magnetSep-magnetThick/2);
  mi_selectlabel(-0.5*magnetThick-0.5*interMagnet,-magnetSep-magnetThick/2);
  mi_selectlabel(1.5*magnetThick+1.5*interMagnet,-magnetSep-magnetThick/2);
  mi_selectlabel(0.5*magnetThick+0.5*interMagnet,-magnetSep-magnetThick/2);
  
  mi_setblockprop("Air", 1, 0, "", 0, 1, 0);
  mi_clearselected();
endif

if steelThick > 0
  mi_addblocklabel(0,magnetSep + magnetThick + 0.5*steelThick);
  mi_selectlabel(0,magnetSep + magnetThick + 0.5*steelThick);
  mi_addblocklabel(0,-magnetSep - magnetThick - 0.5*steelThick);
  mi_selectlabel(0,-magnetSep - magnetThick - 0.5*steelThick);
  mi_setblockprop(steelType, 1, 0, "", 0, 1, 0);
  mi_clearselected();
endif


mi_refreshview()


%%----------------------------------------------------------------------------------------------
mi_saveas('modHalbach-test.FEM')



%%----------------------------------------------------------------------------------------------

% start processing

mi_analyze;
mi_loadsolution; % opens the solution window
mo_showdensityplot(1,0,2.0,0.0,'mag') %%Shows the flux density

mo_addcontour(-1.5*magnetThick-2*interMagnet,0)
mo_addcontour(2.5*magnetThick+2*interMagnet,0)
mo_makeplot(2,200)