From f11b0e5ce4d962c5ec52e858722b5ced7e00b796 Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Wed, 14 Jul 2010 13:09:18 +0200 Subject: [PATCH] Update to ReadUI method & examples all examples may need a revision... --- matlab/ReadUI.m | 27 +++++++++++--- matlab/examples/Coax.m | 15 +++----- matlab/examples/Coax_CylinderCoords.m | 52 ++++++++++++++------------- matlab/examples/Helix.m | 3 -- 4 files changed, 55 insertions(+), 42 deletions(-) diff --git a/matlab/ReadUI.m b/matlab/ReadUI.m index 97ec254..e519acb 100644 --- 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 diff --git a/matlab/examples/Coax.m b/matlab/examples/Coax.m index 417a5e1..629d063 100644 --- 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]); diff --git a/matlab/examples/Coax_CylinderCoords.m b/matlab/examples/Coax_CylinderCoords.m index 0588742..70f80f4 100644 --- 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]); diff --git a/matlab/examples/Helix.m b/matlab/examples/Helix.m index bd337ef..18e97c9 100644 --- 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);