From 09bbf49b1ad6349162021cdf5c98e8fec11e33b0 Mon Sep 17 00:00:00 2001
From: Sebastian Held <sebastian.held@uni-due.de>
Date: Wed, 31 Mar 2010 15:39:46 +0200
Subject: [PATCH] testsuite: new test: cavity.m (resonance frequency)

---
 TESTSUITE/combinedtests/cavity.m          | 124 ++++++++++++++++------
 TESTSUITE/helperscripts/check_frequency.m |  31 ++++++
 2 files changed, 123 insertions(+), 32 deletions(-)
 create mode 100644 TESTSUITE/helperscripts/check_frequency.m

diff --git a/TESTSUITE/combinedtests/cavity.m b/TESTSUITE/combinedtests/cavity.m
index 18310cb..39cbae9 100644
--- a/TESTSUITE/combinedtests/cavity.m
+++ b/TESTSUITE/combinedtests/cavity.m
@@ -7,9 +7,18 @@ ENABLE_PLOTS = 1;
 CLEANUP = 0;        % if enabled and result is PASS, remove simulation folder
 STOP_IF_FAILED = 1; % if enabled and result is FAILED, stop with error
 
-% LIMITS
-upper_error = 0.036; % max +3.6%
-lower_error = 0;     % max -0.0%
+% LIMITS - inside
+lower_rel_limit = 1.3e-3;    % -0.13%
+upper_rel_limit = 1.3e-3;    % +0.13%
+lower_rel_limit_TM = 2.5e-3; % -0.25%
+upper_rel_limit_TM = 0;      % +0%
+min_rel_amplitude = 0.6;     % 60%
+min_rel_amplitude_TM = 0.27; % 27%
+
+% LIMITS - outside
+outer_rel_limit = 0.02;
+max_rel_amplitude = 0.17;
+
 
 % structure
 a = 5e-2;
@@ -19,7 +28,7 @@ if ~((b<a) && (a<d))
     error 'correct the dimensions of the cavity'
 end
 
-f_start = 0;
+f_start = 1e9;
 f_stop = 10e9;
 
 Sim_Path = 'tmp';
@@ -35,23 +44,28 @@ FDTD = SetBoundaryCond(FDTD,BC);
 
 %setup CSXCAD geometry
 CSX = InitCSX();
-mesh.x = linspace(0,a,35);
-mesh.y = linspace(0,b,35);
-mesh.z = linspace(0,d,35);
+grid_res = 2e-3;
+mesh.x = 0:grid_res:a; %linspace(0,a,25);
+mesh.y = 0:grid_res:b; %linspace(0,b,25);
+mesh.z = 0:grid_res:d; %linspace(0,d,25);
 CSX = DefineRectGrid(CSX, 1,mesh);
 
 % excitation
-pos1 = [mesh.x(10) mesh.y(10) mesh.z(10)];
-pos2 = [mesh.x(10) mesh.y(11) mesh.z(10)];
-CSX = AddExcitation(CSX,'excite',1,[0 1 0]);
-CSX = AddBox(CSX, 'excite', 1, pos1, pos2);
+CSX = AddExcitation(CSX,'excite1',0,[1 1 1]);
+p(1,1) = mesh.x(floor(end*2/3));
+p(2,1) = mesh.y(floor(end*2/3));
+p(3,1) = mesh.z(floor(end*2/3));
+p(1,2) = mesh.x(floor(end*2/3)+1);
+p(2,2) = mesh.y(floor(end*2/3)+1);
+p(3,2) = mesh.z(floor(end*2/3)+1);
+CSX = AddCurve( CSX, 'excite1', 0, p );
  
-% %dump
+%dump
 % CSX = AddDump(CSX,'Et_',0,2);
-% pos1 = [mesh.x(1) mesh.y(10) mesh.z(1)];
-% pos2 = [mesh.x(end) mesh.y(10) mesh.z(end)];
+% pos1 = [mesh.x(1) mesh.y(1) mesh.z(1)];
+% pos2 = [mesh.x(end) mesh.y(end) mesh.z(end)];
 % CSX = AddBox(CSX,'Et_',0 , pos1,pos2);
-% 
+
 % %dump
 % CSX = AddDump(CSX,'Et2_',0,2);
 % pos1 = [mesh.x(1) mesh.y(1) mesh.z(1)];
@@ -65,13 +79,19 @@ CSX = AddBox(CSX, 'excite', 1, pos1, pos2);
 % CSX = AddBox(CSX,'Et3_',0 , pos1,pos2);
 
 %voltage calc
-pos1 = [mesh.x(15) mesh.y(15) mesh.z(15)];
-pos2 = [mesh.x(15) mesh.y(16) mesh.z(15)];
+CSX = AddProbe(CSX,'ut1x',0);
+pos1 = [mesh.x(floor(end/4)) mesh.y(floor(end/2))   mesh.z(floor(end/5))];
+pos2 = [mesh.x(floor(end/4)+1) mesh.y(floor(end/2)) mesh.z(floor(end/5))];
+CSX = AddBox(CSX,'ut1x', 0 ,pos1,pos2);
+
 CSX = AddProbe(CSX,'ut1y',0);
+pos1 = [mesh.x(floor(end/4)) mesh.y(floor(end/2))   mesh.z(floor(end/5))];
+pos2 = [mesh.x(floor(end/4)) mesh.y(floor(end/2)+1) mesh.z(floor(end/5))];
 CSX = AddBox(CSX,'ut1y', 0 ,pos1,pos2);
-pos1 = [mesh.x(floor(end/2)) mesh.y(floor(end/2)) mesh.z(end)];
-pos2 = [mesh.x(floor(end/2)) mesh.y(floor(end/2)) mesh.z(end-1)];
+
 CSX = AddProbe(CSX,'ut1z',0);
+pos1 = [mesh.x(floor(end/2)) mesh.y(floor(end/2)) mesh.z(floor(end/5))];
+pos2 = [mesh.x(floor(end/2)) mesh.y(floor(end/2)) mesh.z(floor(end/5)+1)];
 CSX = AddBox(CSX,'ut1z', 0 ,pos1,pos2);
 
 %Write openEMS compatible xml-file
@@ -90,6 +110,16 @@ cd(savePath);
 % analysis
 %
 
+% remove excitation from time series
+t_start = 7e-10; % FIXME to be calculated
+t_idx_start = interp1( UI.TD{1}.t, 1:numel(UI.TD{1}.t), t_start, 'nearest' );
+for n=1:numel(UI.TD)
+    UI.TD{n}.t = UI.TD{n}.t(t_idx_start:end);
+    UI.TD{n}.val = UI.TD{n}.val(t_idx_start:end);
+    [UI.FD{n}.f,UI.FD{n}.val] = FFT_time2freq( UI.TD{n}.t, UI.TD{n}.val );
+end
+
+
 f = UI.FD{1}.f;
 ux = UI.FD{1}.val;
 uy = UI.FD{2}.val;
@@ -106,14 +136,31 @@ uz = uz(f_idx_start:f_idx_stop);
 k = @(m,n,l) sqrt( (m*pi/a)^2 + (n*pi/b)^2 + (l*pi/d)^2 );
 f_TE101 = c0/(2*pi) * k(1,0,1);
 f_TE102 = c0/(2*pi) * k(1,0,2);
-f_TE103 = c0/(2*pi) * k(1,0,3);
 f_TE201 = c0/(2*pi) * k(2,0,1);
 f_TE202 = c0/(2*pi) * k(2,0,2);
 f_TM110 = c0/(2*pi) * k(1,1,0);
-f_TM210 = c0/(2*pi) * k(2,1,0);
+f_TM111 = c0/(2*pi) * k(1,1,1);
+
+f_TE = [f_TE101 f_TE102 f_TE201 f_TE202];
+f_TM = [f_TM110 f_TM111];
+
+% calculate frequency limits
+temp = [f_start f_TE f_stop];
+f_outer1 = [];
+f_outer2 = [];
+for n=1:numel(temp)-1
+    f_outer1 = [f_outer1 temp(n) .* (1+outer_rel_limit)];
+    f_outer2 = [f_outer2 temp(n+1) .* (1-outer_rel_limit)];
+end
+
+temp = [f_start f_TM f_stop];
+f_outer1_TM = [];
+f_outer2_TM = [];
+for n=1:numel(temp)-1
+    f_outer1_TM = [f_outer1_TM temp(n) .* (1+outer_rel_limit)];
+    f_outer2_TM = [f_outer2_TM temp(n+1) .* (1-outer_rel_limit)];
+end
 
-f_TE = [f_TE101 f_TE102 f_TE103 f_TE201 f_TE202];
-f_TM = [f_TM110 f_TM210];
 
 if ENABLE_PLOTS
     figure
@@ -121,25 +168,38 @@ if ENABLE_PLOTS
     max1 = max(abs(uy));
     hold on
     plot( repmat(f_TE,2,1)/1e9, repmat([0; max1],1,numel(f_TE)), 'm-.', 'LineWidth', 2 )
+    plot( (repmat(f_TE,2,1) .* repmat(1-lower_rel_limit,2,numel(f_TE)))/1e9, repmat([0; max1],1,numel(f_TE)), 'r-', 'LineWidth', 1 )
+    plot( (repmat(f_TE,2,1) .* repmat(1+upper_rel_limit,2,numel(f_TE)))/1e9, repmat([0; max1],1,numel(f_TE)), 'r-', 'LineWidth', 1 )
+    plot( (repmat(f_TE,2,1) .* repmat([1-outer_rel_limit;1+outer_rel_limit],1,numel(f_TE)))/1e9, repmat(max1*min_rel_amplitude,2,numel(f_TE)), 'r-', 'LineWidth', 1 ) % freq limits
+    plot( [f_outer1;f_outer2]/1e9, repmat(max1*max_rel_amplitude,2,numel(f_outer1)), 'g-', 'LineWidth', 1 ) % amplitude limits
     xlabel('Frequency (GHz)')
-    legend( {'u_y','theoretical'} );
+    legend( {'u_y','theoretical'} )
+    title( 'TE-modes' )
 
     figure
     plot(f/1e9,abs(uz))
     max1 = max(abs(uz));
     hold on
     plot( repmat(f_TM,2,1)/1e9, repmat([0; max1],1,numel(f_TM)), 'm-.', 'LineWidth', 2 )
+    plot( (repmat(f_TM,2,1) .* repmat(1-lower_rel_limit_TM,2,numel(f_TM)))/1e9, repmat([0; max1],1,numel(f_TM)), 'r-', 'LineWidth', 1 )
+    plot( (repmat(f_TM,2,1) .* repmat(1+upper_rel_limit_TM,2,numel(f_TM)))/1e9, repmat([0; max1],1,numel(f_TM)), 'r-', 'LineWidth', 1 )
+    plot( (repmat(f_TM,2,1) .* repmat([1-lower_rel_limit_TM;1+upper_rel_limit_TM],1,numel(f_TM)))/1e9, repmat(max1*min_rel_amplitude_TM,2,numel(f_TM)), 'r-', 'LineWidth', 1 ) % freq limits
+    plot( [f_outer1_TM;f_outer2_TM]/1e9, repmat(max1*max_rel_amplitude,2,numel(f_outer1_TM)), 'g-', 'LineWidth', 1 ) % amplitude limits
     xlabel('Frequency (GHz)')
-    legend( {'u_z','theoretical'} );
+    legend( {'u_z','theoretical'} )
+    title( 'TM-modes' )
 end
 
-pass = 1;
-% pass = check_limits( Z, upper_limit, lower_limit );
-% if pass
-%     disp( 'combinedtests/Coax.m (characteristic impedance):  pass' );
-% else
-%     disp( 'combinedtests/Coax.m (characteristic impedance):  * FAILED *' );
-% end
+pass1 = check_frequency( f, abs(uy), f_TE*(1+upper_rel_limit),    f_TE*(1-lower_rel_limit),    min_rel_amplitude, 'inside' );
+pass2 = check_frequency( f, abs(uz), f_TM*(1+upper_rel_limit_TM), f_TM*(1-lower_rel_limit_TM), min_rel_amplitude_TM, 'inside' );
+pass3 = check_frequency( f, abs(uy), f_outer2, f_outer1, max_rel_amplitude, 'outside' );
+pass4 = check_frequency( f, abs(uz), f_outer2_TM, f_outer1_TM, max_rel_amplitude, 'outside' );
+pass = pass1 && pass2 && pass3 && pass4;
+if pass
+    disp( 'combinedtests/cavity.m (resonance frequency):  pass' );
+else
+    disp( 'combinedtests/cavity.m (resonance frequency):  * FAILED *' );
+end
 
 
 
diff --git a/TESTSUITE/helperscripts/check_frequency.m b/TESTSUITE/helperscripts/check_frequency.m
new file mode 100644
index 0000000..33a7cca
--- /dev/null
+++ b/TESTSUITE/helperscripts/check_frequency.m
@@ -0,0 +1,31 @@
+function pass = check_frequency( f, val, f_upper, f_lower, rel_amplitude, type )
+
+pass = true;
+max1 = max(val);
+
+if numel(f_upper) ~= numel(f_lower)
+    error 'inconsistant vectors'
+end
+
+for n=1:numel(f_upper)
+    f1 = f_lower(n);
+    f2 = f_upper(n);
+    f1_idx = interp1( f, 1:numel(f), f1, 'nearest' );
+%    if f(f1_idx) < f1, f1_idx = f1_idx + 1; end
+    f2_idx = interp1( f, 1:numel(f), f2, 'nearest' );
+%    if f(f2_idx) > f2, f2_idx = f2_idx - 1; end
+
+    if strcmp( type, 'inside' )
+        if max( val(f1_idx:f2_idx) ) < max1 * rel_amplitude
+            pass = false;
+            return
+        end
+    elseif strcmp( type, 'outside' )
+        if max( val(f1_idx:f2_idx) ) > max1 * rel_amplitude
+            pass = false;
+            return
+        end
+    else
+        error 'unsupported operation'
+    end
+end