// Auto-construction of Magnetic Gear
// Mustafa Husain
// 2-10-2011 (2 years in Japan)

// Opening SciFemm Toolbox that connects FEMM and SciLab 
exec('C:\Program Files\femm42\scifemm\scifemm.sci',-1);
openfemm
// Creating a new Magnetostatics document to work on 
newdocument(0);
//                            
// * * * * * * * * * * * * * * GEAR PARAMETERS * * * * * * * * * * * * * * 
PH=4;                                          //Pole pairs of HSR
PL=10;                                        //Pole pairs of LSR
Ns=PH+PL;                                    //Poles of PP
La=26;                                      //Axial Thickness of Gear
Lm1=4.5;                                   //Thickness of HSR Magnet
Lm2=3;                                    //Thickness of LSR Magnet
Lair=5;                                    //Thickness of Air Region
apr1=0.5;                                //arc to pitch ratio of HSR
apr3=0.75;                              //arc to pitch ratio of LSR
apr2=0.5;                              //arc to pitch ratio of PP 
pitch1=360/(2*PH)*apr1*%pi/180;       //Pole pitch of HSR
pitch3=360/(2*PL)*apr3*%pi/180;      //Pole pitch of LSR
pitch2=360/Ns*apr2*%pi/180;         //Pole pitch of PP
Lc1=12.5;                          //Thickness of HSR steel core
Lc2=7;                            //Thickness of LSR steel core
Lp=10;                           //Radial Thickness of PP
g=0.5;                         //Air gap 
R1=7;                         //Shaft Radius
R2=R1+Lc1;                   //HSR Inner Radius
R3=R2+Lm1;                  //HSR Outer Radius
R4=R3+g;                   //Inner Radius of PP
R5=R4+Lp;                 //Outer Radius of PP
R6=R5+g;                 //Inner Radius of LSR
R7=R6+Lm2;              //Outer Radius of LSR (Inner Radius of LSR Core)
R8=R7+Lc2;             //Outer Radius of LSR Core
R9=R8+Lair;           //Radius of Air Region
// * * * * * * * * * * * * * * PROBLEM DEFINITION * * * * * * * * * * * * * 
Frequency=0;            //Frequency in Hz
Unit='millimeters';    //or 'centimeters' 'meters' 'inches''mils' 'micrometers'
Type='planar';        //Problem type 'planar' or 'axi'
Precision=1.e-8;     //Accuracy of solution
Depth=La;           //Axial length
minAngle=10;       //Minimum angle for solving
acSolver=0;       //AC solver method 0=Successive approximation 1=Newton 
mi_probdef(Frequency, Unit, Type,Precision ,Depth, minAngle,acSolver);

// * * * * * * * * * * * * *  * * DRAWING * * * * * * * * * * * * * * * * * 

// ................................SHAFT ..................................

mi_drawarc([-R1 0; R1 0], 180, 2.5);
mi_drawarc([R1 0; -R1 0], 180, 2.5);
mi_selectarcsegment(0,-R1)
mi_selectarcsegment(0,R1)
mi_setarcsegmentprop(2.5, '<none>', 0,0)

// .............................INNER ROTOR ...............................

// 1. Steel Core
mi_drawarc([-R2 0; R2 0], 180, 2.5);
mi_drawarc([R2 0; -R2 0], 180, 2.5);
// 2. Permanent Magnet
mi_drawarc([-R3 0; R3 0], 180, 2.5);
mi_drawarc([R3 0; -R3 0], 180, 2.5);
    //Grouping arcs to Group 1
    mi_selectarcsegment(0,-R2);
    mi_selectarcsegment(0,R2);
    mi_selectarcsegment(0,-R3);
    mi_selectarcsegment(0,R3);
    mi_setarcsegmentprop(2.5, '<none>', 0,1);
    //Drawing one pole
    mi_drawline(R2*cos(pitch1/2),R2*sin(pitch1/2),R3*cos(pitch1/2),R3*sin(pitch1/2));
    mi_drawline(R2*cos(-pitch1/2),R2*sin(-pitch1/2),R3*cos(-pitch1/2),R3*sin(-pitch1/2));
//Grouping lines to Group 1
    mi_selectsegment((R2+R3)/2*cos(pitch1/2),(R2+R3)/2*sin(pitch1/2));
    mi_selectsegment((R2+R3)/2*cos(-pitch1/2),(R2+R3)/2*sin(-pitch1/2));
    mi_setsegmentprop( '<none>', 1, 0, 0, 1)
    //Deleting adjacent arcs
    mi_selectarcsegment(0,R3);
    mi_selectarcsegment(0,-R3);
    mi_deleteselectedarcsegments
    //Adding (2*PH-1) poles
    mi_seteditmode('group')
    mi_selectgroup(1)
    mi_copyrotate(0, 0, (pitch1/apr1*180/%pi), (2*PH-1) )

// ............................ POLE PIECES ...............................

//Drawing one pole
mi_drawarc([-R4 0; R4 0], 180, 2.5);
mi_drawarc([R4 0; -R4 0], 180, 2.5);
mi_drawarc([-R5 0; R5 0], 180, 2.5);
mi_drawarc([R5 0; -R5 0], 180, 2.5);
mi_selectarcsegment(0,-R4)
mi_selectarcsegment(0,R4)
mi_selectarcsegment(0,-R5)
mi_selectarcsegment(0,R5)
mi_setarcsegmentprop(2.5, '<none>', 0,2)
mi_drawline(R4*cos(pitch2/2),R4*sin(pitch2/2),R5*cos(pitch2/2),R5*sin(pitch2/2))
mi_drawline(R4*cos(-pitch2/2),R4*sin(-pitch2/2),R5*cos(-pitch2/2),R5*sin(-pitch2/2))
mi_selectsegment((R4+R5)/2*cos(pitch2/2),(R4+R5)/2*sin(pitch2/2));
mi_selectsegment((R4+R5)/2*cos(-pitch2/2),(R4+R5)/2*sin(-pitch2/2));
mi_setsegmentprop( '<none>', 1, 0, 0, 2)
mi_selectarcsegment(0,R4);
mi_selectarcsegment(0,R5);
mi_selectarcsegment(0,-R4);
mi_selectarcsegment(0,-R5);
mi_deleteselectedarcsegments
mi_seteditmode('group')
mi_selectgroup(2)
mi_copyrotate(0, 0, (pitch2/apr2*180/%pi), (Ns-1) )

// ...............................OUTER ROTOR .............................

// 1. Permanent Magnet
mi_drawarc([-R6 0; R6 0], 180, 2.5);
mi_drawarc([R6 0; -R6 0], 180, 2.5);
mi_selectarcsegment(0,-R6)
mi_selectarcsegment(0,R6)
mi_setarcsegmentprop(2.5, '<none>', 0,3)
mi_drawarc([-R7 0; R7 0], 180, 2.5);
mi_drawarc([R7 0; -R7 0], 180, 2.5);
mi_selectarcsegment(0,-R7)
mi_selectarcsegment(0,R7)
mi_setarcsegmentprop(2.5, '<none>', 0,3)
mi_drawline(R6*cos(pitch3/2),R6*sin(pitch3/2),R7*cos(pitch3/2),R7*sin(pitch3/2))
mi_drawline(R6*cos(-pitch3/2),R6*sin(-pitch3/2),R7*cos(-pitch3/2),R7*sin(-pitch3/2))
mi_selectsegment((R6+R7)/2*cos(pitch3/2),(R6+R7)/2*sin(pitch3/2));
mi_selectsegment((R6+R7)/2*cos(-pitch3/2),(R6+R7)/2*sin(-pitch3/2));
mi_setsegmentprop( '<none>', 1, 0, 0, 3)
mi_selectarcsegment(0,R6);
mi_selectarcsegment(0,-R6);
mi_deleteselectedarcsegments
mi_seteditmode('group')
mi_selectgroup(3)
mi_copyrotate(0, 0, (pitch3/apr3*180/%pi), (2*PL-1) )
// 2. Steel Core
mi_drawarc([-R8 0; R8 0], 180, 2.5);
mi_drawarc([R8 0; -R8 0], 180, 2.5);
mi_selectarcsegment(0,-R8)
mi_selectarcsegment(0,R8)
mi_setarcsegmentprop(2.5, '<none>', 0,3)

// ............................AIR REGION .................................

mi_drawarc([-R9 0; R9 0], 180, 2.5);
mi_drawarc([R9 0; -R9 0], 180, 2.5);
mi_selectarcsegment(0,-R9)
mi_selectarcsegment(0,R9)

// * * * * * * * * * * * * * * MATERIAL LIBRARY * * * * * * * * * * * * * *
//NdFeB Permanent Magnet NMX-38SR of 38MGOe 
//NdFeB Permanent Magnet N42SH of 42MGOe
//JFE Steel Silicon-Steel Core 50JN400
mi_addmaterial('Air', 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0);
mi_addmaterial('NMX-38SR', 1.0718, 1.0718, 891000, 0, 0.625, 0, 0, 1, 0, 0, 0);
mi_addmaterial('N42SH',1.0305,1.0305,1000000, 0, 0.7142857, 0, 0, 1, 0, 0, 0);
mi_addmaterial('50JN400',1,1,0, 0, 2.12766, 0.5, 0, 0.96, 0, 0, 0);
bhcurve = [0.,1.31,1.45,1.49,1.53,1.543,1.57,1.577,1.59,1.61,1.62,1.63,1.6424,1.651,1.6634,1.6714,1.6812,1.689,1.7046,1.7124,1.7182,1.725,1.7305, 1.7366,   1.7423,1.7493,1.755,1.7607,1.766, 1.7736, 1.7793, 1.7853,  1.7888, 1.7945,    1.7986,1.8028,1.8081,1.81285,1.8134,1.817,1.8216,1.8287,1.8353,1.842,1.8487,  1.8555,1.8623,1.869,1.8755,1.8823,1.889,1.8977,1.9063,1.9148,1.9235,1.932,    1.94357, 1.9555,1.9673,1.9793,1.991,2.00657,2.0225,2.0383,2.054,2.07052,2.0873,2.104; 0., 250., 500., 750., 1000.,1250.,1500.,1750.,2000.,2250., 2500., 2750., 3000., 3250., 3500., 3750., 4000., 4250., 4500., 4750., 5000., 5250., 5500., 5750., 6000., 6250., 6500., 6750., 7000., 7250., 7500., 7750., 8000.,  8250., 8500., 8750., 9000., 9250., 9500., 9750., 10000., 10500.,11000.,11500.,12000.,12500., 13000.,13500.,14000.,14500.,15000.,16000.,17000.,18000.,19000.,20000., 22000.,24000.,26000.,28000.,30000.,35000.,40000.,45000.,50000.,60000.,70000.,80000.]';
mi_addbhpoints('50JN400', bhcurve);
// ........................................................................
// * * * * * * * * * *  APPLYING MATERIAL TO BLOCKS  * * * * * * * * * * *
//mi_setblockprop(’blockname’,automesh,meshsize,’incircuit',magdir,
// group,turns)

// ................................SHAFT..................................
mi_addblocklabel(0,0) 
mi_selectlabel(0,0);
mi_setblockprop('Air', 0, 1, '<None>', 0, 0, 0);
mi_clearselected

// ..............................INNER ROTOR..............................
// 1. Core 
mi_addblocklabel(R1+((R2-R1)/2),0) 
mi_selectlabel(R1+((R2-R1)/2),0);
mi_setblockprop('50JN400', 0, 1, '<None>', 0, 1, 0);
mi_clearselected
// 2. Permanent Magnet
//Making 1st pole
mi_addblocklabel(R2+((R3-R2)/2),0);
mi_selectlabel(R2+((R3-R2)/2),0);
mi_setblockprop('NMX-38SR', 0, 1, '<None>', 180, 1, 0);
mi_clearselected;
//Making 2nd pole
mi_addblocklabel((R2+((R3-R2)/2))*cos(pitch1/apr1),(R2+((R3-R2)/2))*sin(pitch1/apr1));
mi_selectlabel((R2+((R3-R2)/2))*cos(pitch1/apr1),(R2+((R3-R2)/2))*sin(pitch1/apr1));
mi_setblockprop('NMX-38SR', 0, 1, '<None>', pitch1/apr1*180/%pi, 1, 0);
mi_clearselected;
mi_selectlabel(R2+((R3-R2)/2),0);
mi_selectlabel((R2+((R3-R2)/2))*cos(pitch1/apr1),(R2+((R3-R2)/2))*sin(pitch1/apr1));
//Copying one pole pair to make other poles
mi_copyrotate(0, 0, 2*(pitch1/apr1*180/%pi), PH-1 )
mi_clearselected

// ..............................POLE PIECES..............................

mi_addblocklabel(R4+((R5-R4)/2),0) 
mi_selectlabel(R4+((R5-R4)/2),0);
mi_setblockprop('50JN400', 0, 1, '<None>', 0, 2, 0);
mi_clearselected
mi_selectlabel(R4+((R5-R4)/2),0);
mi_copyrotate(0, 0, (pitch2/apr2*180/%pi), Ns-1 )
mi_clearselected

// ..............................OUTER ROTOR ..............................

// 1. Permanent Magnet
//Making 1st pole
mi_addblocklabel(R6+((R7-R6)/2),0);
mi_selectlabel(R6+((R7-R6)/2),0);
mi_setblockprop('N42SH', 0, 1, '<None>', 180, 3, 0);
mi_clearselected;
//Making 2nd pole
mi_addblocklabel((R6+((R7-R6)/2))*cos(pitch3/apr3),(R6+((R7-R6)/2))*sin(pitch3/apr3));
mi_selectlabel((R6+((R7-R6)/2))*cos(pitch3/apr3),(R6+((R7-R6)/2))*sin(pitch3/apr3));
mi_setblockprop('N42SH', 0, 1, '<None>', pitch3/apr3*180/%pi, 3, 0);
mi_clearselected;
mi_selectlabel(R6+((R7-R6)/2),0);
mi_selectlabel((R6+((R7-R6)/2))*cos(pitch3/apr3),(R6+((R7-R6)/2))*sin(pitch3/apr3));
//Copying one pole pair to make other poles
mi_copyrotate(0, 0, 2*(pitch3/apr3*180/%pi), PL-1 )
mi_clearselected

// 2. Core 
mi_addblocklabel(R7+((R8-R7)/2),0) 
mi_selectlabel(R7+((R8-R7)/2),0);
mi_setblockprop('50JN400', 0, 1, '<None>', 0, 3, 0);
mi_clearselected

// ...............................AIR REGION ..............................
//Air gap 
mi_addblocklabel(R3+((R4-R3)/2),0) 
mi_selectlabel(R3+((R4-R3)/2),0);
mi_setblockprop('Air', 0,R4-R3 , '<None>', 0, 0, 0);
mi_clearselected
//Outside air
mi_addblocklabel(R8+((R9-R8)/2),0) 
mi_selectlabel(R8+((R9-R8)/2),0);
mi_setblockprop('Air', 0, 1, '<None>', 0, 0, 0);
mi_clearselected

// * * * * * * * * * * * * * * SETTING BOUNDARY * * * * * * * * * * * * * *     

mi_addboundprop('A=0', 0, 0, 0, 0, 0, 0, 0, 0, 0);
mi_selectarcsegment(0,-R9);
mi_setarcsegmentprop(2.5, 'A=0', 0, 0);
mi_clearselected
mi_selectarcsegment(0,R9);
mi_setarcsegmentprop(2.5, 'A=0', 0, 0);
mi_clearselected
mi_zoomnatural
mi_saveas('D:\PhD\Simulations\Testing\test.fem');
