update to matlab interface + helix example

matlab/ReadUI.m

@@ -0,0 +1,28 @@
+function UI = ReadUI(files, path)
+
+if (nargin<2)
+    path ='';
+end
+
+if (ischar(files))
+    filenames{1}=files;
+else
+    filenames=files;
+end
+
+for n=1:numel(filenames)
+    tmp = load([path filenames{n}]);
+    t = tmp(:,1)';
+    val = tmp(:,2)';
+
+    UI.TD{n}.t = t;
+    UI.TD{n}.val = val;
+    
+    dt=t(2)-t(1);
+    val = [val zeros(1,5000)];
+    L=numel(val);
+    UI.FD{n}.f = (0:L-1)/L/dt;
+    UI.FD{n}.f = UI.FD{n}.f(1:floor(L/2));
+    UI.FD{n}.val = fft(val)/L;
+    UI.FD{n}.val = UI.FD{n}.val(1:floor(L/2));
+end

matlab/examples/Helix.m

@@ -8,14 +8,18 @@ coil_rad = 10;
 coil_length = 50;
 coil_turns = 8;
 coil_res = 10;
-port_length = 10;
+port_length = coil_length;
 port_resist = 50;
 
+f_max = 50e6;
+f_excite = 1e9;
+mesh_size = wire_rad;
 
 openEMS_Path = [pwd() '/../../']
 openEMS_opts = '';
-% openEMS_opts = [openEMS_opts ' --disable-dumps'];
+openEMS_opts = [openEMS_opts ' --disable-dumps'];
 % openEMS_opts = [openEMS_opts ' --debug-material'];
+% openEMS_opts = [openEMS_opts ' --debug-operator'];
 
 Sim_Path = 'tmp';
 Sim_CSX = 'helix.xml';
@@ -24,15 +28,21 @@ mkdir(Sim_Path);
 
 %setup FDTD parameter
 FDTD = InitFDTD(5e5,1e-6);
-FDTD = SetGaussExcite(FDTD,0.5e9,0.5e9);
+FDTD = SetGaussExcite(FDTD,f_excite/2,f_excite/2);
 BC = [1 1 1 1 1 1];
 FDTD = SetBoundaryCond(FDTD,BC);
 
+add_Lines = mesh_size * 1.5.^(1:10)
+add_Lines = add_Lines(find(add_Lines<(3e8/f_max)*1e3))
+
 %setup CSXCAD geometry
 CSX = InitCSX();
-mesh.x = [-35,-25,-20,-15,-12,-11,-10.5,-10,-9.5,-9,-8.5,-8,-7.5,-7,-6.5,-6,-5.5,-5,-4.5,-4,-3.5,-3,-2.5,-2,-1.5,-1,-0.5,0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8,8.5,9,9.5,10,10.5,11,12,13.5,15,17,18,19,19.5,20,20.5,21,22,23,25,27.5,30,35,45];
+mesh.x = -coil_rad-mesh_size : mesh_size : coil_rad+mesh_size+feed_length;
-mesh.y = [-35,-25,-20,-15,-12,-11,-10.5,-10,-9.5,-9,-8.5,-8,-7.5,-7,-6.5,-6,-5.5,-5,-4.5,-4,-3.5,-3,-2.5,-2,-1.5,-1,-0.5,0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8,8.5,9,9.5,10,10.5,11,12,13,15,17.5,20,25,35];
+mesh.x = [mesh.x(1)-add_Lines mesh.x mesh.x(end)+add_Lines ];
-mesh.z = [-25,-15,-10,-5,-2,-1,-0.5,0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8,8.5,9,9.5,10,10.5,11,11.5,12,12.5,13,13.5,14,14.5,15,15.5,16,16.5,17,17.5,18,18.5,19,19.5,20,20.5,21,21.5,22,22.5,23,23.5,24,24.5,25,25.5,26,26.5,27,27.5,28,28.5,29,29.5,30,30.5,31,31.5,32,32.5,33,33.5,34,34.5,35,35.5,36,36.5,37,37.5,38,38.5,39,39.5,40,40.5,41,41.5,42,42.5,43,43.5,44,44.5,45,45.5,46,46.5,47,47.5,48,48.5,49,49.5,50,50.5,51,52,53,55,57.5,60,65,75];
+mesh.y = -coil_rad-mesh_size : mesh_size : coil_rad+mesh_size;
+mesh.y = [mesh.y(1)-add_Lines mesh.y mesh.y(end)+add_Lines ];
+mesh.z = -mesh_size : mesh_size : coil_length+mesh_size;
+mesh.z = [mesh.z(1)-add_Lines mesh.z mesh.z(end)+add_Lines ];
 CSX = DefineRectGrid(CSX, 1e-3,mesh);
 
 %create copper helix and feed lines...
@@ -43,7 +53,13 @@ CSX = SetMaterialProperty(CSX,'copper','Kappa',56e6);
 dt = 1.0/coil_res;
 height=0;
 wire.Vertex = {};
-count=0;
+p(1,1) = coil_rad + feed_length;
+p(2,1) = 0;
+p(3,1) = 0.5*(coil_length-port_length);
+p(1,2) = coil_rad + feed_length;
+p(2,2) = 0;
+p(3,2) = 0;
+count=2;
 for n=0:coil_turns-1
     for m=0:coil_res
         count = count + 1;
@@ -53,22 +69,21 @@ for n=0:coil_turns-1
     end
     height = height + coil_length/coil_turns;
 end
+p(1,count+1) = coil_rad + feed_length;
+p(2,count+1) = 0;
+p(3,count+1) = coil_length;
+p(1,count+2) = coil_rad + feed_length;
+p(2,count+2) = 0;
+p(3,count+2) = 0.5*(coil_length+port_length);
 CSX = AddWire(CSX, 'copper', 0, p, wire_rad);
 
-start = [coil_rad, 0 , 0];stop = [coil_rad+feed_length, 0 , 0];
-CSX = AddCylinder(CSX,'copper',0 ,start,stop,wire_rad);
-
-start(3)=coil_length;stop(3)=coil_length;
-CSX = AddCylinder(CSX,'copper',0 ,start,stop,wire_rad);
-
-start(3)=0;start(1)=coil_rad+feed_length;
-CSX = AddCylinder(CSX,'copper',0 ,start,stop,wire_rad);
-
 CSX = AddMaterial(CSX,'resist');
-kappa = port_length/port_resist/wire_rad^2/pi/1e-3
+kappa = port_length/port_resist/wire_rad^2/pi/1e-3;
 CSX = SetMaterialProperty(CSX,'resist','Kappa',kappa);
 
-start(3)=(coil_length-port_length)/2;stop(3)=(coil_length+port_length)/2;
+start=[coil_rad+feed_length 0 (coil_length-port_length)/2];
+stop=[coil_rad+feed_length 0 (coil_length+port_length)/2];
+%start(3)=(coil_length-port_length)/2;stop(3)=(coil_length+port_length)/2;
 CSX = AddCylinder(CSX,'resist',5 ,start,stop,wire_rad);
 
 %excitation
@@ -77,7 +92,7 @@ CSX = AddCylinder(CSX,'excite', 0 ,start,stop,wire_rad);
 
 %voltage calc
 CSX = AddProbe(CSX,'ut1',0);
-CSX = AddBox(CSX,'ut1', 0 ,start,stop);
+CSX = AddBox(CSX,'ut1', 0 ,stop,start);
 
 %current calc
 CSX = AddProbe(CSX,'it1',1);
@@ -86,6 +101,17 @@ start(1) = start(1)-2;start(2) = start(2)-2;
 stop(1) = stop(1)+2;stop(2) = stop(2)+2;
 CSX = AddBox(CSX,'it1', 0 ,start,stop);
 
+%dump
+CSX = AddDump(CSX,'Et_',0,0);
+start = [mesh.x(1) , 0 , mesh.z(1)];
+stop = [mesh.x(end) , 0 , mesh.z(end)];
+CSX = AddBox(CSX,'Et_',0 , start,stop);
+
+CSX = AddDump(CSX,'Ht_',1,0);
+start = [mesh.x(1) , 0 , mesh.z(1)];
+stop = [mesh.x(end) , 0 , mesh.z(end)];
+CSX = AddBox(CSX,'Ht_',0 , start,stop);
+
 %Write openEMS compatoble xml-file
 WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX);
 
@@ -97,3 +123,26 @@ disp(command);
 system(command)
 cd(savePath);
 
+%%%post-proc
+U = ReadUI('ut1','tmp/');
+I = ReadUI('it1','tmp/');
+
+Z = U.FD{1}.val./I.FD{1}.val;
+f = U.FD{1}.f;
+L = imag(Z)./(f*2*pi);
+R = real(Z);
+ind = find(f<f_max);
+
+subplot(2,1,1);
+plot(f(ind)*1e-6,L(ind)*1e9,'Linewidth',2);
+xlabel('frequency (MHz)');
+ylabel('coil inductance (nH)');
+grid on;
+subplot(2,1,2);
+plot(f(ind)*1e-6,R(ind),'Linewidth',2);
+xlabel('frequency (MHz)');
+ylabel('resistance (\Omega)');
+grid on;
+
+
+