clc; close all; clear;

openfemm;
newdocument(0);  

%Design Variables 
A=47.5;	%mm
B=57;	%mm
C=38;	%mm
D=4;	%mm
E=38;	%mm
F=28.5;	%mm
U=9.5;	%mm
V=9.5;	%mm
W=19;	%mm
X=9.5;	%mm
Y=28.5;	%mm
Clearance=1;	%mm
Current=0.2;	%Amps
turns=1000;

% Problem Definition
mi_probdef(0,'millimeters','planar', 1e-8, (B/2), (30), (0));  
mi_saveas('EI-57.fem');

% Materials & Current Definition (*mi_modifymaterial("BlockName",propnum,value))
mi_getmaterial('Air');
mi_getmaterial('1018 Steel');
mi_getmaterial('30 AWG');
mi_modifymaterial('1018 Steel',9,0);	% 9: LamType- 0:In plane
mi_modifymaterial('1018 Steel',6,0.5);	% 6: Lamination thickness(mm)
mi_modifymaterial('1018 Steel',8,0.85);	% 8: LammFill(fill fraction)

mi_addcircprop('Coil', Current, 1);

% Core design (I)
mi_drawpolygon([-B/2,Y/2; -B/2,A/2; B/2,A/2; B/2,Y/2]); 

mi_addblocklabel(0,E/2);
mi_selectlabel(0,E/2);
mi_setblockprop('1018 Steel',0,0,0,0,0,0);
mi_clearselected;

% Core desing (E)
mi_drawpolygon([-B/2,Y/2; -B/2,-A/2; B/2,-A/2; B/2,Y/2; B/2-V,Y/2; B/2-V,-Y/2; W/2,-Y/2; W/2,Y/2; -W/2,Y/2; -W/2,-Y/2; -B/2+V,-Y/2; -B/2+V,Y/2]); 

mi_addblocklabel(0,0);
mi_selectlabel(0,0);
mi_setblockprop('1018 Steel',0,0,0,0,0,0);
mi_clearselected;

% Air gap (inner openning)
mi_addblocklabel(-W/2-X/2,0);
mi_selectlabel(-W/2-X/2,0);
mi_setblockprop('Air',0,0,0,0,0,0);
mi_clearselected;

mi_addblocklabel(W/2+X/2,0);
mi_selectlabel(W/2+X/2,0);
mi_setblockprop('Air',0,0,0,0,0,0);
mi_clearselected;

% Holes design (*using mi_addarc(x1,y1,x2,y2,angle,maxseg))
% Left-up
mi_addnode(-F/2,E/2+D/2);
mi_addnode(-F/2,E/2-D/2);
mi_addarc(-F/2,E/2+D/2,-F/2,E/2-D/2,180,2);
mi_addarc(-F/2,E/2-D/2,-F/2,E/2+D/2,180,2);

mi_addblocklabel(-F/2,E/2);
mi_selectlabel(-F/2,E/2);
mi_setblockprop('Air',0,0,0,0,0,0);
mi_clearselected;

%Right-up
mi_addnode(F/2,E/2+D/2);
mi_addnode(F/2,E/2-D/2);
mi_addarc(F/2,E/2+D/2,F/2,E/2-D/2,180,2);
mi_addarc(F/2,E/2-D/2,F/2,E/2+D/2,180,2);

mi_addblocklabel(F/2,E/2);
mi_selectlabel(F/2,E/2);
mi_setblockprop('Air',0,0,0,0,0,0);
mi_clearselected;

%Left-down
mi_addnode(-F/2,-E/2+D/2);
mi_addnode(-F/2,-E/2-D/2);
mi_addarc(-F/2,-E/2+D/2,-F/2,-E/2-D/2,180,2);
mi_addarc(-F/2,-E/2-D/2,-F/2,-E/2+D/2,180,2);

mi_addblocklabel(-F/2,-E/2);
mi_selectlabel(-F/2,-E/2);
mi_setblockprop('Air',0,0,0,0,0,0);
mi_clearselected;

%Right-down
mi_addnode(F/2,-E/2+D/2);
mi_addnode(F/2,-E/2-D/2);
mi_addarc(F/2,-E/2+D/2,F/2,-E/2-D/2,180,2);
mi_addarc(F/2,-E/2-D/2,F/2,-E/2+D/2,180,2);

mi_addblocklabel(F/2,-E/2);
mi_selectlabel(F/2,-E/2);
mi_setblockprop('Air',0,0,0,0,0,0);
mi_clearselected;

%Left Coil design
mi_drawpolygon([-W/2,-Clearance; -W/2-X+Clearance,-Clearance; -W/2-X+Clearance,-Y/2+Clearance; -W/2,-Y/2+Clearance]);

mi_addblocklabel(-W/2-X/2,-Y/3);
mi_selectlabel(-W/2-X/2,-Y/3);
mi_setblockprop('30 AWG',0,0,'Coil',0,0,turns); 
mi_clearselected;

%Right Coil design
mi_drawpolygon([W/2,-Clearance; W/2+X-Clearance,-Clearance; W/2+X-Clearance,-Y/2+Clearance; W/2,-Y/2+Clearance]);

mi_addblocklabel(W/2+X/2,-Y/3);
mi_selectlabel(+W/2+X/2,-Y/3);
mi_setblockprop('30 AWG',0,0,'Coil',0,0,-turns);
mi_clearselected;

% Boundary
mi_addnode(0,5*max(A,B)/2);
mi_addnode(0,-5*max(A,B)/2);
mi_addarc(0,5*max(A,B)/2,0,-5*max(A,B)/2,180,2);
mi_addarc(0,-5*max(A,B)/2,0,5*max(A,B)/2,180,2);

mi_addblocklabel(0,5*max(A,B)/3);
mi_selectlabel(0,5*max(A,B)/3);
mi_setblockprop('Air',0,0,0,0,0,0);
mi_clearselected;

% Mesh
mi_createmesh;
mi_zoom(-B/2-5,-A/2-5,B/2+5,A/2+5);
mi_savemetafile('EI-57_mesh.emf'); 

% Magnetic field lines
mi_analyze;
mi_loadsolution;
mo_zoom(-B/2-5,-A/2-5,B/2+5,A/2+5);
mo_savemetafile('EI-57_2a.emf'); 

% Magnetic flux density
mo_showdensityplot(1,0,0,2,'mag')
mo_zoom(-B/2-5,-A/2-5,B/2+5,A/2+5);
mo_savemetafile('EI-57_2b.emf');

% Max and min values of magnetic flux density
B_tot=zeros(101,101);
max_b=-10;
min_b=10;
i=1;
for xi=-B/2:B/100:B/2
	j=1;
    for yj=-A/2:A/100:A/2
		B_tot(i,j)=norm(mo_getb(xi,yj));
		if max_b<=B_tot(i,j)
			max_b=B_tot(i,j);
		end
		if min_b>B_tot(i,j)
			min_b=B_tot(i,j);
        end
     j=j+1;
    end
    i=i+1;
end

% Pcu and Wc calculation
mo_selectblock(-W/2-X/2,-Y/3);
mo_selectblock(W/2+X/2,-Y/3); 
Pcu=mo_blockintegral(4); 	% 4:Resistive losses (Watt)
Wc=mo_blockintegral(2); 	% 2:Magnetic field energy (Joule)
mo_clearblock;
 
% Calculation of amplitude of the magnetic flux density B, along the horizontal axis
B_horiz=zeros(1,1000);
x_axis=zeros(1,1000);

for i=1:1000
	x_init=-B/2-B/1000;
	x_calc=-B/2+i*B/1000;
	B_horiz(1,i)=norm(mo_getb(x_calc,0));
	x_axis(1,i)=x_calc;
end

f1=figure;
plot(x_axis,B_horiz);
title('Amplitude of Magnetic Flux Density (Horizontal)');
xlabel('x (mm)');
ylabel('|B| (T)');
saveas(f1,'EI-57_2d_h.emf');

% Calculation of amplitude of the magnetic flux density B, along the vertical axis
B_vert=zeros(1,1000);
y_axis=zeros(1,1000);

for i=1:1000
	y_init=-A/2-A/1000;
	y_calc=-A/2+i*A/1000;
	B_vert(1,i)=norm(mo_getb(0,y_calc));
	y_axis(1,i)=y_calc;
end

f2=figure;
plot(B_vert,y_axis);
title('Amplitude of Magnetic Flux Density (Vertical)');
xlabel('|B| (T)');
ylabel('y (mm)');
saveas(f2,'EI-57_2d_v.emf');

% Calculation of magnetic flux
% Phi
mo_addcontour(-W/2,0);
mo_addcontour(W/2,0);
Phi_vec=mo_lineintegral(0);	%0:total B.n
Phi=Phi_vec(1); %Weber
mo_clearcontour;

% Phi1
mo_addcontour(-B/2,0);
mo_addcontour(-W/2-X,0);
Phi_vec=mo_lineintegral(0);
Phi1=Phi_vec(1);    %Weber
mo_clearcontour;

% Phi2
mo_addcontour(W/2+X,0);
mo_addcontour(B/2,0);
Phi_vec=mo_lineintegral(0);
Phi2=Phi_vec(1);    %Weber
mo_clearcontour;