diff --git a/matlab/Dump2VTK.m b/matlab/Dump2VTK.m index fb8100c..5c01435 100644 --- a/matlab/Dump2VTK.m +++ b/matlab/Dump2VTK.m @@ -26,21 +26,21 @@ if (mesh.type==0) %write cartesian mesh to vtk fprintf(fid,'DIMENSIONS %d %d %d\n',numel(x),numel(y),numel(z)); - fprintf(fid,'X_COORDINATES %d float\n',numel(x)); + fprintf(fid,'X_COORDINATES %d double\n',numel(x)); fprintf(fid,'%e',x(1)); for n=2:numel(x) fprintf(fid,' %e',x(n)); end fprintf(fid,'\n'); - fprintf(fid,'Y_COORDINATES %d float\n',numel(y)); + fprintf(fid,'Y_COORDINATES %d double\n',numel(y)); fprintf(fid,'%e',y(1)); for n=2:numel(y) fprintf(fid,' %e',y(n)); end fprintf(fid,'\n'); - fprintf(fid,'Z_COORDINATES %d float\n',numel(z)); + fprintf(fid,'Z_COORDINATES %d double\n',numel(z)); fprintf(fid,'%e',z(1)); for n=2:numel(z) fprintf(fid,' %e',z(n)); @@ -54,7 +54,7 @@ elseif (mesh.type==1) %write cylindrical mesh to vtk fprintf(fid,'DIMENSIONS %d %d %d\n',numel(x),numel(y),numel(z)); - fprintf(fid,'POINTS %d float\n',numel(x)*numel(y)*numel(z)); + fprintf(fid,'POINTS %d double\n',numel(x)*numel(y)*numel(z)); for nz=1:numel(z) for ny=1:numel(y) @@ -82,9 +82,9 @@ fprintf(fid,'POINT_DATA %d\n',numel(x)*numel(y)*numel(z)); % dump vector field data if (ndims(fields)==4) if (nargin>3) - fprintf(fid,['VECTORS ' fieldname ' float\n']); + fprintf(fid,['VECTORS ' fieldname ' double\n']); else - fprintf(fid,'VECTORS field float\n'); + fprintf(fid,'VECTORS field double\n'); end fclose(fid); field_x = fields(:,:,:,1); @@ -101,9 +101,9 @@ end % dump scalar field data if (ndims(fields)==3) if (nargin>3) - fprintf(fid,['SCALARS ' fieldname ' float 1\nLOOKUP_TABLE default\n']); + fprintf(fid,['SCALARS ' fieldname ' double 1\nLOOKUP_TABLE default\n']); else - fprintf(fid,'SCALARS field float 1\nLOOKUP_TABLE default\n'); + fprintf(fid,'SCALARS field double 1\nLOOKUP_TABLE default\n'); end fclose(fid); dumpField = fields(:); diff --git a/matlab/GetField_Interpolation.m b/matlab/GetField_Interpolation.m index df0d4c1..787ef65 100644 --- a/matlab/GetField_Interpolation.m +++ b/matlab/GetField_Interpolation.m @@ -70,8 +70,10 @@ if (numel(x)==1) [Y_I Z_I] = ndgrid(y_i,z_i); for n=1:numel(field.values) field_i.values{n}(1,:,:,1) = interpn(Y,Z,squeeze(field.values{n}(1,:,:,1)),Y_I,Z_I); - field_i.values{n}(1,:,:,2) = interpn(Y,Z,squeeze(field.values{n}(1,:,:,2)),Y_I,Z_I); - field_i.values{n}(1,:,:,3) = interpn(Y,Z,squeeze(field.values{n}(1,:,:,3)),Y_I,Z_I); + if (size(field_i.values{n},4)>1) + field_i.values{n}(1,:,:,2) = interpn(Y,Z,squeeze(field.values{n}(1,:,:,2)),Y_I,Z_I); + field_i.values{n}(1,:,:,3) = interpn(Y,Z,squeeze(field.values{n}(1,:,:,3)),Y_I,Z_I); + end end return; end @@ -81,8 +83,10 @@ if (numel(y)==1) [X_I Z_I] = ndgrid(x_i,z_i); for n=1:numel(field.values) field_i.values{n}(:,1,:,1) = interpn(X,Z,squeeze(field.values{n}(:,1,:,1)),X_I,Z_I); - field_i.values{n}(:,1,:,2) = interpn(X,Z,squeeze(field.values{n}(:,1,:,2)),X_I,Z_I); - field_i.values{n}(:,1,:,3) = interpn(X,Z,squeeze(field.values{n}(:,1,:,3)),X_I,Z_I); + if (size(field_i.values{n},4)>1) + field_i.values{n}(:,1,:,2) = interpn(X,Z,squeeze(field.values{n}(:,1,:,2)),X_I,Z_I); + field_i.values{n}(:,1,:,3) = interpn(X,Z,squeeze(field.values{n}(:,1,:,3)),X_I,Z_I); + end end return; end @@ -92,8 +96,10 @@ if (numel(z)==1) [X_I Y_I] = ndgrid(x_i,y_i); for n=1:numel(field.values) field_i.values{n}(:,:,1,1) = interpn(X,Y,squeeze(field.values{n}(:,:,1,1)),X_I,Y_I); - field_i.values{n}(:,:,1,2) = interpn(X,Y,squeeze(field.values{n}(:,:,1,2)),X_I,Y_I); - field_i.values{n}(:,:,1,3) = interpn(X,Y,squeeze(field.values{n}(:,:,1,3)),X_I,Y_I); + if (size(field_i.values{n},4)>1) + field_i.values{n}(:,:,1,2) = interpn(X,Y,squeeze(field.values{n}(:,:,1,2)),X_I,Y_I); + field_i.values{n}(:,:,1,3) = interpn(X,Y,squeeze(field.values{n}(:,:,1,3)),X_I,Y_I); + end end return; end @@ -104,6 +110,8 @@ end [X_I Y_I Z_I] = ndgrid(x_i,y_i,z_i); for n=1:numel(field.values) field_i.values{n}(:,:,:,1) = interpn(X,Y,Z,field.values{n}(:,:,:,1),X_I,Y_I,Z_I); - field_i.values{n}(:,:,:,2) = interpn(X,Y,Z,field.values{n}(:,:,:,2),X_I,Y_I,Z_I); - field_i.values{n}(:,:,:,3) = interpn(X,Y,Z,field.values{n}(:,:,:,3),X_I,Y_I,Z_I); + if (size(field_i.values{n},4)>1) + field_i.values{n}(:,:,:,2) = interpn(X,Y,Z,field.values{n}(:,:,:,2),X_I,Y_I,Z_I); + field_i.values{n}(:,:,:,3) = interpn(X,Y,Z,field.values{n}(:,:,:,3),X_I,Y_I,Z_I); + end end diff --git a/matlab/GetField_Range.m b/matlab/GetField_Range.m new file mode 100644 index 0000000..973c4f6 --- /dev/null +++ b/matlab/GetField_Range.m @@ -0,0 +1,52 @@ +function [field_i mesh_i] = GetField_Range(field, mesh, range) +% [field_i mesh_i] = GetField_Range(field, mesh, range) +% +% Get a field dump subset within a given mesh range +% +% example: +% % specify a mesh range +% range{1} = [0 150] * 1e-3; +% range{2} = [0 0]; +% range{3} = [-850 -400] * 1e-3; +% +% % read hdf data +% [field mesh] = ReadHDF5Dump('Et.h5'); +% % extract a ranged subset +% [field_i mesh_i] = GetField_Range(field, mesh, range); +% +% %or both steps in one with the same result: +% [field_i mesh_i] = ReadHDF5Dump('Et.h5','Range', range); +% +% openEMS matlab interface +% ----------------------- +% author: Thorsten Liebig +% +% See also ReadHDF5Dump ReadHDF5FieldData ReadHDF5Mesh + +mesh_i = mesh; +for n=1:3 + ind_range{n} = find( mesh.lines{n}>=range{n}(1) & mesh.lines{n}<=range{n}(2)); + + if (isempty(ind_range{n})) + ind_range{n} = find( mesh.lines{n}>=range{n}(1) , 1, 'first'); + end + if (isempty(ind_range{n})) + ind_range{n} = find( mesh.lines{n}<=range{n}(2) , 1, 'last'); + end + + mesh_i.lines{n} = mesh.lines{n}(ind_range{n}); +end + +field_i = field; + +if (isfield(field,'FD')) + for n=1:numel(field.FD.values) + field_i.FD.values{n} = field.FD.values{n}(ind_range{1},ind_range{2},ind_range{3},:); + end +end + +if (isfield(field,'TD')) + for n=1:numel(field.TD.values) + field_i.TD.values{n} = field.TD.values{n}(ind_range{1},ind_range{2},ind_range{3},:); + end +end diff --git a/matlab/GetField_TD2FD.m b/matlab/GetField_TD2FD.m index 4bf6c91..e4dc60f 100644 --- a/matlab/GetField_TD2FD.m +++ b/matlab/GetField_TD2FD.m @@ -21,6 +21,7 @@ if (~isfield(field,'TD')) end t = field.TD.time; +dt = t(2)-t(1); clear field.FD @@ -35,7 +36,7 @@ numTS = numel(field.TD.values); for n=1:numTS for nf = 1:numel(freq) f = freq(nf); - field.FD.values{nf} = field.FD.values{nf} + 2/numTS * field.TD.values{n}.*exp(-1i*2*pi*f*t(n)); + field.FD.values{nf} = field.FD.values{nf} + field.TD.values{n}.*exp(-1i*2*pi*f*t(n)) * 2 * dt; % t(n) is absolute time and therefore the half-timestep offset of % the H-field is automatically compensated % openEMS output: E-fields start at t=0 diff --git a/matlab/ReadHDF5Dump.m b/matlab/ReadHDF5Dump.m index b263f95..365948f 100644 --- a/matlab/ReadHDF5Dump.m +++ b/matlab/ReadHDF5Dump.m @@ -5,12 +5,15 @@ function [field mesh] = ReadHDF5Dump(file, varargin) % transformation. % % possible arguments: +% 'Range' see GetField_Range % 'Interpolation' see GetField_Interpolation % 'Frequency' see GetField_TD2FD % % example: % [field mesh] = ReadHDF5Dump('Et.h5'); % or +% [field mesh] = ReadHDF5Dump('Et.h5','Range',{[0 100],[-20 20],[50 90]}); +% or % [field mesh] = ReadHDF5Dump('Et.h5','Interpolation',[21 1 101],'Frequency',300e6); % % openEMS matlab interface @@ -18,12 +21,18 @@ function [field mesh] = ReadHDF5Dump(file, varargin) % author: Thorsten Liebig % % See also ReadHDF5Mesh ReadHDF5FieldData GetField_Interpolation -% GetField_TD2FD +% GetField_TD2FD GetField_Range field = ReadHDF5FieldData(file); mesh = ReadHDF5Mesh(file); if (nargin>1) + for n=1:2:(nargin-1) + if (strcmp(varargin{n},'Range')==1); + [field mesh] = GetField_Range(field, mesh, varargin{n+1}) + end + end + for n=1:2:(nargin-1) if (strcmp(varargin{n},'Interpolation')==1); [field mesh] = GetField_Interpolation(field,mesh,varargin{n+1});