clc; clear all; close all

disp('Projekt "Simulation und Entwicklung eines induktiven Sensors zur Abstandsmessung"');
disp(' ');
disp('In diesem Projekt wird ein induktiver Sensor nach dem Wirbelstromprinzip zur Messung');
disp('des Abstandes zwischen dem Messobjekt aus Aluminium und der Sensorfl?che entwickelt.');
disp(' ');

% Zuerst muss eine Schnittstelle zwischen MATLAB und FEMM hergestellt werden
addpath('C:\femm42\mfiles');

% FEMM wird geöffnet
openfemm;

% Ein neues Dokument vom Typ "Magnetostatic" wird geöffnet
newdocument(0);

% Dimensionierung der einzelnen Elemente
disp('Dimensionierung der einzelnen Elemente');
disp(' ');

disp('1. Sensorgehäuse:');

disp('Innendurchmesser Ds_i des Sensorgeh?uses (in mm):');
Ds_i = 8
Rs_i = Ds_i / 2;

disp('Au?endurchmesser Ds_a des Sensorgeh?uses (in mm):');
Ds_a = 12
Rs_a = Ds_a / 2;

disp('L?nge Ls des Sensorgeh?uses (in mm):');
Ls = 60

disp('Material des Sensorgeh?uses: Edelstahl');
disp(' ');
disp(' ');

disp('2. Schalenkern:');

disp('Durchmesser Kernloch D_sk vom Schalenkern (in mm):');
D_sk = 1.05
R_sk = D_sk / 2;

disp('Durchmesser Bolzen D_sb vom Schalenkern (in mm):');
D_sb = 3
R_sb = D_sb / 2;

disp('Innendurchmesser der Schale D_ssi vom Schalenkern (in mm):');
D_ssi = 5.8
R_ssi = D_ssi / 2;

disp('Au?endurchmesser der Schale D_ssa vom Schalenkern (in mm):');
D_ssa = 7.35
R_ssa = D_ssa / 2;

disp('Gesamtl?nge L_sg vom Schalenkern (in mm):');
L_sg = 3.6

disp('Tiefe der Aussparung L_sa vom Schalenkern (in mm):');
L_sa = 2.8

disp('Material des Schalenkerns: M33');
disp('');
disp('');

disp('3. Spule:');

disp('Drahtdurchmesser D_sd der Spule (in mm):');
D_sd = 0.5

disp('Anzahl der max. Wicklungen N pro Lage von der Spule:');
N_max = floor(L_sa / D_sd)
N = 1:1:N_max;

disp('Anzahl der max. Lagen M von der Spule:');
M_max = floor((R_ssi - R_sb)/D_sd)
M = 1:1:M_max;

disp('Anzahl der max. Wicklungen N_ges von der Spule:');
N_ges = N_max * M_max

disp('Material der Spule: 24 AWG (Kupfer)');
disp('');
disp('');

disp('4. Messobjekt:');
disp('Durchmesser D_m des Messobjektes (in mm):');
D_m = 12
R_m = D_m / 2;

disp('L?nge L_m des Messobjektes (in mm):');
a = 1;
L_m = D_m * a

disp('Material des Messobjektes: Aluminium 6061-T6');
disp('');

disp('maximaler Abstand s_max zwischen Messobjekt und Sensorfl?che (in mm):');
s_max = 1.1

disp('minimalster Abstand s_min zwischen Messobjekt und Sensorfl?che (in mm):');
s_min = 0.1

disp('Verschiebung des Messobjektes S_Schritt (in mm) jeweils um:');
s_Schritt = 0.1

disp('Stromfluss durch die Spule (in A):');
i = 0.1

% Definierung des Problemtypen
% (0=Magnetostatic; Einheit in Millimeter; achsensymmetrisch; Aufl?sung; Tiefe der Bauteile im Dokument; Winkel)
disp('Frequenz an der Spule (in Hz):');
f=10^6;

mi_probdef(f, 'millimeters', 'axi', 1.e-8, 0, 30);


% Sensorgeh?use wird gezeichnet
mi_drawrectangle(Rs_i, 0, Rs_a, -Ls);

% Schalenkern wird gezeichnet
mi_drawline([R_sk 0; R_sb 0])
mi_drawline([R_sb 0; R_sb -L_sa])
mi_drawline([R_sb -L_sa; R_ssi -L_sa])
mi_drawline([R_ssi -L_sa; R_ssi 0])
mi_drawline([R_ssi 0; R_ssa 0])
mi_drawline([R_ssa 0; R_ssa -L_sg])
mi_drawline([R_ssa -L_sg; R_sk -L_sg])
mi_drawline([R_sk -L_sg; R_sk 0])

% Spule wird gezeichnet
Sx = M_max*D_sd + R_sb;
Sy = N_max*D_sd - L_sa;
mi_drawrectangle([R_sb -L_sa; Sx Sy]);

% Messobjekt wird gezeichnet
My = L_m + s_max;
mi_drawrectangle([0 s_max; R_m My]);

% Kreis f?r die Umgebung wird gezeichnet
mi_drawarc([0 -100; 0 100], 180, 2.5);
mi_addsegment([0 -100; 0 100]);

% F?gt Knotenpunkte an den einzelnen Elementen zur Bestimmung des Materials hinzu
mi_addblocklabel(20,0); %Luft
mi_addblocklabel(5,-2); %Sensorgeh?use
mi_addblocklabel(2,-3); %Schalenkern
mi_addblocklabel(1.6,-2.7); %Spule
M_ly = s_max + 0.2;
mi_addblocklabel(2,M_ly); %Messobjekt

% Definiere die Eigenschaften der achsensymmetrischen Umgebung
muo = pi*4.e-7;
mi_addboundprop('Asymptotic', 0, 0, 0, 0, 0, 0, 1/(muo*0.2), 0, 2);

% Anwendung der Eigenschaften der achsensymmetrischen Umgebung auf den Umgebungskreis
mi_selectarcsegment(100,0);
mi_setarcsegmentprop(2.5, 'Asymptotic', 0, 0);

load('HBDiagramm_M33.mat')
% Add some materials properties
%mi addmaterial(?matname?, mu x, mu y, H c, J, Cduct, Lam d, Phi hmax, lam fill,
%LamType, Phi hx, Phi hy, nstr, dwire) adds a new material with called ?matname?
%with the material properties:
mi_addmaterial('Air', 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0);
mi_addmaterial('24 AWG', 1, 1, 0, 0, 58, 0, 0, 1, 0, 0, 0);
mi_addmaterial('M33', 1, 1, 80, 0, 0, 0, 0, 1, 0, 0, 0);
bhcurve = HBDiagramm_M33;
mi_addbhpoints('M33', bhcurve);
mi_addmaterial('Aluminium', 1, 1, 0, 0, 24.59, 0, 0, 1, 0, 0, 0);
mi_addmaterial('StainlessSteel304', 1, 1, 0, 0, 1.45, 0, 0, 1, 0, 0, 0);


% Add a "circuit property" so that we can calculate the properties of the
% coil as seen from the terminals.

mi_addcircprop('Spule', i, 1);

% Apply the materials to the appropriate block labels

mi_selectlabel(20,0);
mi_setblockprop('Air', 1, 0, '<None>', 0, 0, 0);
mi_clearselected

mi_selectlabel(1.6,-2.7);
mi_setblockprop('24 AWG', 1, 0, 'Spule', 0, 0, N_ges);
mi_clearselected

mi_selectlabel(2,-3);
mi_setblockprop('M33', 1, 0, '<None>', 90, 0, 0);
mi_clearselected

mi_selectlabel(5,-2);
mi_setblockprop('StainlessSteel304', 1, 0, '<None>', 0, 0, 0);
mi_clearselected

mi_selectlabel(2,M_ly);
mi_setblockprop('Aluminium', 1, 0, '<None>', 0, 99, 0);
mi_clearselected

mi_selectnode(0, s_max);
mi_selectnode(0,My);
mi_selectnode(R_m,s_max);
mi_selectnode(R_m,My);
mi_setnodeprop('Aluminium',99)
mi_clearselected

% Now, the finished input geometry can be displayed.
mi_zoom(-4, -9, 4,9)

% We have to give the geometry a name before we can analyze it.
mi_saveas('coil.fem');

% Now,analyze the problem and load the solution when the analysis is finished
ergebnis = []

%Schleife
for z= s_max:-s_Schritt:s_min
	mi_analyze()
	mi_loadsolution()
    mo_zoom(-4, -9, 4,9)
	%mo_groupselectblock(99)
    t = mo_getcircuitproperties('Spule')
    ergebnis = [ergebnis;t]
	if (z>s_min) 
		mi_selectgroup(99)
    	mi_movetranslate(0,-s_Schritt)
	end
end

% {i, v, \[Phi]} = MOGetCircuitProperties["icoil"]
% If we were interested in inductance, it could be obtained by
% dividing flux linkage by current
T=table
T.x=(s_max:-s_Schritt:s_min)'
T.current=ergebnis(:,1)
T.volt=ergebnis(:,2)
T.flux=ergebnis(:,3)
subplot(4,1,1);
plot(T.x,abs(T.volt))
grid on
ylabel('U [V]\rightarrow')
xlabel('Abstand s [mm]\rightarrow')
title('Spannungsabfall')
subplot(4,1,2);
plot(T.x,imag(T.flux))
grid on
ylabel('\Phi [Weber]\rightarrow')
xlabel('Abstand s [mm]\rightarrow')
title('Flussverkettung an der Spule')
subplot(4,1,3);
plot(T.x,imag(T.flux)/i)
grid on
ylabel('L [H]\rightarrow')
xlabel('Abstand s [mm]\rightarrow')
title('Selbstinduktion an der Spule')
subplot(4,1,4);
plot(T.x,real(T.flux)/i)
grid on
ylabel('R [Ohm]\rightarrow')
xlabel('Abstand s [mm]\rightarrow')
title('Widerstand an der Spule')

% Netliste für die Schaltung im LTspice
%* D:\THGA\Simulation elektrischer Systeme\FEMM\Induktiver Sensor.asc
%L1 N001 N002 {L}
%C1 N002 N001 {C}
%R2 N002 N001 {R}
%XU1 N003 N004 N007 N005 N006 N009 N008 +10 0 LTC1043
%XU2 +5 N001 +10 0 N002 LT1013
%XU3 +5 N002 +10 0 N008 LT1013
%R3 N001 N008 20k
%R4 N003 N002 10k
%R5 N005 +5 10k
%XU4 N006 N003 N011 N004 N005 N012 N008 +10 0 LTC1043
%C2 N007 N009 100n
%C3 N011 N012 100n
%V1 N010 0 12
%R6 N010 +10 20 tol=1 pwr=0.1
%R7 +10 +5 50 tol=1 pwr=0.1
%R8 +5 0 50 tol=1 pwr=0.1
%.tran 0 200m 0 1u startup
%.options gmin=1e-10
%.options method=gear
%.params L = 10m C = 1.6u R = 10k
%.lib LTC.lib
%.backanno
%.end
netlist = ['Induktiver Sensor.net']