Update to ReadUI method & examples

all examples may need a revision...
2010-07-14 13:09:18 +02:00 · 2010-07-14 13:09:18 +02:00 · f11b0e5ce4
parent 0d2811ff21
commit f11b0e5ce4
4 changed files with 55 additions and 42 deletions
--- a/matlab/ReadUI.m
+++ b/matlab/ReadUI.m
@ -1,12 +1,23 @@
-function UI = ReadUI(files, path)
-% function UI = ReadUI(files, path)
+function UI = ReadUI(files, path, freq)
+% function UI = ReadUI(files, path, freq)
 %
 % read current and voltages from multiple files found in path
 %
 % returns voltages/currents in time and frequency-domain
+% 
+% remarks on the frequency-domain:
+% - all signals are assumed to start at t=0
+% - currents that e.g. start at t = +delta_t/2 will be phase shifted by
+%    exp(-j*w*t(1))
+%
+% parameter:
+% freq (optional):
+%     frequency-domain values will be calculated according to 'freq'
+%     if 'freq' is not given, a FFT will be used
 %
 % e.g.
-% UI = ReadUI({'ut1_1','ut1_2','it1'},'tmp/');
+% U = ReadUI({'ut1_1','ut1_2'},'tmp' );
+% I = ReadUI('it1'            ,'tmp',[0.5e9 1e9 1.5e9]);
 % 
 % openEMS matlab interface
 % -----------------------
@ -32,5 +43,13 @@ for n=1:numel(filenames)
    UI.TD{n}.t = t;
    UI.TD{n}.val = val;
    
-    [UI.FD{n}.f,UI.FD{n}.val] = FFT_time2freq( t,val );
+    if (nargin<3)
+        [UI.FD{n}.f,UI.FD{n}.val] = FFT_time2freq( t,val );
+    else
+        UI.FD{n}.f = freq;
+        UI.FD{n}.val = DFT_time2freq( t, val, freq );
+    end
+    
+    %correct phase error for time-shifted signals
+    UI.FD{n}.val = UI.FD{n}.val .* exp(-1j*2*pi*UI.FD{n}.f * t(1)); 
 end
--- a/matlab/examples/Coax.m
+++ b/matlab/examples/Coax.m
@ -24,7 +24,7 @@ openEMS_opts = '';
 % openEMS_opts = [openEMS_opts ' --debug-operator'];

 % openEMS_opts = [openEMS_opts ' --disable-dumps --engine=fastest'];
-openEMS_opts = [openEMS_opts ' --engine=sse-compressed'];
+openEMS_opts = [openEMS_opts ' --engine=fastest'];

 Sim_Path = 'tmp';
 Sim_CSX = 'coax.xml';
@ -101,25 +101,18 @@ CSX = AddBox(CSX,'it1', 0 ,start,stop);
 WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX);

 %% cd to working dir and run openEMS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-savePath = pwd();
-cd(Sim_Path); %cd to working dir
-args = [Sim_CSX ' ' openEMS_opts];
-invoke_openEMS(args)
-cd(savePath);
+RunOpenEMS(Sim_Path,Sim_CSX,openEMS_opts);

 %% postproc & do the plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-UI = ReadUI({'ut1_1','ut1_2','it1'},'tmp/');
+UI = ReadUI({'ut1_1','ut1_2','it1'},Sim_Path);
 u_f = (UI.FD{1}.val + UI.FD{2}.val)/2;          %averaging voltages to fit current
 i_f = UI.FD{3}.val;

-delta_t = UI.TD{3}.t(1) - UI.TD{1}.t(1);        % half time-step (s) 
-i_f2 = i_f .* exp(-1i*2*pi*UI.FD{1}.f*delta_t); % compensate half time-step advance of H-field
-
 ZL = Z0/2/pi/sqrt(epsR)*log(coax_rad_ai/coax_rad_i); %analytic line-impedance of a coax
 plot(UI.FD{1}.f,ZL*ones(size(u_f)),'g');
 hold on;
 grid on;
-Z = u_f./i_f2;
+Z = u_f./i_f;
 plot(UI.FD{1}.f,real(Z),'Linewidth',2);
 plot(UI.FD{1}.f,imag(Z),'r','Linewidth',2);
 xlim([0 2*f0]);
--- a/matlab/examples/Coax_CylinderCoords.m
+++ b/matlab/examples/Coax_CylinderCoords.m
@ -10,8 +10,8 @@ Z0 = sqrt(MUE0/EPS0);

 f0 = 0.5e9;
 epsR = 1;
+kappa = 0;

-abs_length = 500;
 length = 3000;
 port_dist = 1500;
 rad_i  = 100;
@ -24,6 +24,7 @@ mesh_res = [max_mesh 2*pi/N_alpha max_mesh];
 %% define and openEMS options %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 openEMS_opts = '';
 openEMS_opts = [openEMS_opts ' --disable-dumps'];
+% openEMS_opts = [openEMS_opts ' --debug-operator'];
 % openEMS_opts = [openEMS_opts ' --debug-material'];

 Sim_Path = 'tmp';
@ -35,7 +36,7 @@ end
 mkdir(Sim_Path);

 %% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-FDTD = InitCylindricalFDTD(1e5,1e-5,'OverSampling',10);
+FDTD = InitCylindricalFDTD(1e5,1e-6,'OverSampling',10);
 FDTD = SetGaussExcite(FDTD,f0,f0);
 BC = [0 0 1 1 0 0];
 FDTD = SetBoundaryCond(FDTD,BC);
@ -49,34 +50,34 @@ mesh.z = 0 : mesh_res(3) : length;
 CSX = DefineRectGrid(CSX, 1e-3,mesh);

 %% fake pml %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-finalKappa = 0.3/abs_length^4;
+abs_length = 30*(mesh.z(2)-mesh.z(1))
+finalKappa = 0.3;
 finalSigma = finalKappa*MUE0/EPS0/epsR;
 CSX = AddMaterial(CSX,'pml');
-CSX = SetMaterialProperty(CSX,'pml','Kappa',finalKappa,'Epsilon',epsR);
+CSX = SetMaterialProperty(CSX,'pml','Kappa',finalKappa+kappa,'Epsilon',epsR);
 CSX = SetMaterialProperty(CSX,'pml','Sigma',finalSigma);
-CSX = SetMaterialWeight(CSX,'pml','Kappa',['pow(abs(z)-' num2str(length-abs_length) ',4)']);
-CSX = SetMaterialWeight(CSX,'pml','Sigma',['pow(abs(z)-' num2str(length-abs_length) ',4)']);
+CSX = SetMaterialWeight(CSX,'pml','Kappa',['pow(abs(z)-' num2str(length-abs_length) ',4)/' num2str(abs_length^4)]);
+CSX = SetMaterialWeight(CSX,'pml','Sigma',['pow(abs(z)-' num2str(length-abs_length) ',4)/' num2str(abs_length^4)]);

 start = [rad_i mesh.y(1) length-abs_length];
 stop = [rad_a mesh.y(end) length];
 CSX = AddBox(CSX,'pml',0 ,start,stop);

-
+%% material
 CSX = AddMaterial(CSX,'fill');
-CSX = SetMaterialProperty(CSX,'fill','Epsilon',epsR);
+CSX = SetMaterialProperty(CSX,'fill','Epsilon',epsR,'Kappa',kappa);
 start = [mesh.x(1) mesh.y(1) 0];
 stop = [mesh.x(end) mesh.y(end) length];
 CSX = AddBox(CSX,'fill',0 ,start,stop);

-start = [rad_i mesh.y(1) 0];
-stop = [rad_a mesh.y(end) 0];
-
 %% apply the excitation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 CSX = AddExcitation(CSX,'excite',0,[1 0 0]);
 weight{1} = '1/rho';
 weight{2} = 0;
 weight{3} = 0;
 CSX = SetExcitationWeight(CSX, 'excite', weight );
+start = [rad_i mesh.y(1) 0];
+stop = [rad_a mesh.y(end) 0];
 CSX = AddBox(CSX,'excite',0 ,start,stop);
 
 %% define dump boxes... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@ -97,35 +98,38 @@ CSX = AddProbe(CSX,'ut1_2',0);
 start = [ rad_i 0 port_dist+mesh_res(3) ];stop = [ rad_a 0 port_dist+mesh_res(3) ];
 CSX = AddBox(CSX,'ut1_2', 0 ,start,stop);

+
+CSX = AddProbe(CSX,'ut_ex',0);
+start = [ rad_i 0 0 ];stop = [ rad_a 0 0 ];
+CSX = AddBox(CSX,'ut_ex', 0 ,start,stop);
+
 % current calc
 CSX = AddProbe(CSX,'it1',1);
 mid = 0.5*(rad_i+rad_a);
 start = [ 0 mesh.y(1) port_dist ];stop = [ mid mesh.y(end) port_dist ];
 CSX = AddBox(CSX,'it1', 0 ,start,stop);

-%Write openEMS compatoble xml-file
+%% run openEMS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX);
-
-%% cd to working dir and run openEMS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-savePath = pwd();
-cd(Sim_Path); %cd to working dir
-args = [Sim_CSX ' ' openEMS_opts];
-invoke_openEMS(args);
-cd(savePath);
+RunOpenEMS(Sim_Path,Sim_CSX,openEMS_opts);

 %% postproc & do the plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-UI = ReadUI({'ut1_1','ut1_2','it1'},'tmp/');
+UI = ReadUI({'ut1_1','ut1_2','it1'},Sim_Path);
+plot(UI.TD{1}.t,UI.TD{1}.val)
+
+UI_ex = ReadUI({'ut_ex'},'tmp/');
+hold on;
+plot(UI_ex.TD{1}.t,UI_ex.TD{1}.val,'r--');
+
 u_f = (UI.FD{1}.val + UI.FD{2}.val)/2;          %averaging voltages to fit current
 i_f = UI.FD{3}.val;

-delta_t = UI.TD{3}.t(1) - UI.TD{1}.t(1);        % half time-step (s) 
-i_f2 = i_f .* exp(-1i*2*pi*UI.FD{1}.f*delta_t); % compensate half time-step advance of H-field
-
+figure
 ZL = Z0/2/pi/sqrt(epsR)*log(rad_a/rad_i); %analytic line-impedance of a coax
 plot(UI.FD{1}.f,ZL*ones(size(u_f)),'g');
 hold on;
 grid on;
-Z = u_f./i_f2;
+Z = u_f./i_f;
 plot(UI.FD{1}.f,real(Z),'Linewidth',2);
 plot(UI.FD{1}.f,imag(Z),'r','Linewidth',2);
 xlim([0 2*f0]);
--- a/matlab/examples/Helix.m
+++ b/matlab/examples/Helix.m
@ -127,9 +127,6 @@ RunOpenEMS(Sim_Path, Sim_CSX, openEMS_opts);
 U = ReadUI('ut1','tmp/');
 I = ReadUI('it1','tmp/');

-delta_t_2 = I.TD{1}.t(1) - U.TD{1}.t(1);                        % half time-step (s) 
-I.FD{1}.val = I.FD{1}.val .* exp(-1i*2*pi*I.FD{1}.f*delta_t_2); % compensate half time-step advance of H-field
-
 Z = U.FD{1}.val./I.FD{1}.val;
 f = U.FD{1}.f;
 L = imag(Z)./(f*2*pi);