ProcessFields now writes TD data into HDF5 group: /FieldData/TD + adapted matlab HDF5

Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>

Common/processfields.cpp

@@ -87,6 +87,8 @@ void ProcessFields::InitProcess()
 		delete group;
 		group = new H5::Group( file->createGroup( "/FieldData/FD" ));
 		delete group;
+		group = new H5::Group( file->createGroup( "/FieldData/TD" ));
+		delete group;
 		delete file;
 	}
 }
@@ -423,7 +425,7 @@ bool ProcessFields::DumpVectorArray2HDF5(string filename, string name, FDTD_FLOA
 
 	H5::H5File file( FILE_NAME, H5F_ACC_RDWR );
 
-	H5::Group group( file.openGroup( "/FieldData" ));
+	H5::Group group( file.openGroup( "/FieldData/TD" ));
 
 	hsize_t dimsf[4];              // dataset dimensions
 

matlab/GetField_Interpolation.m

@@ -49,11 +49,19 @@ mesh_i.lines{1} = x_i;
 mesh_i.lines{2} = y_i;
 mesh_i.lines{3} = z_i;
 
-NULL = zeros(numel(x_i),numel(y_i),numel(z_i),3);
+if (isfield(field,'TD'))
-for n=1:numel(field.values)
+    field_i.TD = interpolate_fields(field.TD,x,y,z, x_i, y_i, z_i);
-    field_i.values{n} = NULL;
+    field_i.TD.time = field.TD.time;
 end
-clear NULL;
+
+if (isfield(field,'FD'))
+    field_i.FD = interpolate_fields(field.FD,x,y,z, x_i, y_i, z_i);
+    field_i.FD.freq = field.FD.freq;
+end
+
+return
+
+function field_i = interpolate_fields(field, x,y,z, x_i, y_i, z_i)
 
 % matlab cannot handle 3D data to be 2D data, workaround for these cases
 if (numel(x)==1)

matlab/GetField_TD2FD.m

@@ -1,5 +1,5 @@
-function field_FD = GetField_TD2FD(field_TD, freq)
+function field = GetField_TD2FD(field, freq)
-% function field_FD = GetField_TD2FD(field_TD, freq)
+% function field = GetField_TD2FD(field, freq)
 %
 % Transforms time-domain field data into the frequency domain
 % Autocorrects the half-timestep offset of the H-field
@@ -15,20 +15,27 @@ function field_FD = GetField_TD2FD(field_TD, freq)
 %
 % See also ReadHDF5FieldData
 
-t = field_TD.time;
+if (~isfield(field,'TD'))
+    warning('openEMS:GetField_TD2FD','field has no time domain data... skipping FD transformation...');
+    return
+end
 
-field_FD.freq = freq;
+t = field.TD.time;
+
+clear field.FD
+
+field.FD.freq = freq;
 
 for nf = 1:numel(freq)
-    field_FD.values{nf} = 0;
+    field.FD.values{nf} = 0;
 end
     
-numTS = numel(field_TD.values);
+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} + 2/numTS * field.TD.values{n}.*exp(-1i*2*pi*f*t(n));
         % 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

matlab/PlotHDF5FieldData.m

@@ -27,9 +27,9 @@ fields = ReadHDF5FieldData(file);
 if (mesh.type==0)
     % cartesian mesh
     [X Y Z] = meshgrid(mesh.lines{1},mesh.lines{2},mesh.lines{3});
-    for n=1:numel(fields.values)
+    for n=1:numel(fields.TD.values)
         % since Matlab 7.1SP3 the field needs to be reordered
-        fields.values{n} = permute(fields.values{n},[2 1 3 4]); % reorder: y,x,z (or y,x)
+        fields.TD.values{n} = permute(fields.TD.values{n},[2 1 3 4]); % reorder: y,x,z (or y,x)
     end
 else
     disp(['PlotHDF5FieldData:: Error: unknown mesh type ' num2str(mesh.type)]);
@@ -38,14 +38,14 @@ end
 max_amp = 0;
 
 if (component>0)
-    for n=1:numel(fields.values)
+    for n=1:numel(fields.TD.values)
-        Field{n} = fields.values{n}(:,:,:,component);
+        Field{n} = fields.TD.values{n}(:,:,:,component);
     end
 else 
-    for n=1:numel(fields.values)
+    for n=1:numel(fields.TD.values)
-        fx = fields.values{n}(:,:,:,1);
+        fx = fields.TD.values{n}(:,:,:,1);
-        fy = fields.values{n}(:,:,:,2);
+        fy = fields.TD.values{n}(:,:,:,2);
-        fz = fields.values{n}(:,:,:,3);
+        fz = fields.TD.values{n}(:,:,:,3);
         Field{n} = sqrt(fx.^2 + fy.^2 + fz.^2);
     end        
 end

matlab/ReadHDF5FieldData.m

@@ -25,21 +25,51 @@ if isOctave
 end
 
 info = hdf5info(file);
+TD.names = {};
+FD.names = {};
+hdf_fielddata = [];
 
 for n=1:numel(info.GroupHierarchy.Groups)
     if strcmp(info.GroupHierarchy.Groups(n).Name,'/FieldData')
-        for m=1:numel(info.GroupHierarchy.Groups(n).Datasets)
+        %found /FieldData, look for either TD or FD data
-            names{m} = info.GroupHierarchy.Groups(n).Datasets(m).Name;
+        for nGroup=1:numel(info.GroupHierarchy.Groups(n).Groups)
-            hdf_fielddata.time(m) = double(info.GroupHierarchy.Groups(n).Datasets(m).Attributes.Value);
+            %search and read TD data
+            if strcmp(info.GroupHierarchy.Groups(n).Groups(nGroup).Name,'/FieldData/TD')
+                for m=1:numel(info.GroupHierarchy.Groups(n).Groups(nGroup).Datasets)
+                    TD.names{m} = info.GroupHierarchy.Groups(n).Groups(nGroup).Datasets(m).Name;
+                    hdf_fielddata.TD.time(m) = double(info.GroupHierarchy.Groups(n).Groups(nGroup).Datasets(m).Attributes.Value);
+                end
+            end
+            %search and read FD data
+            if strcmp(info.GroupHierarchy.Groups(n).Groups(nGroup).Name,'/FieldData/FD')
+                for m=1:numel(info.GroupHierarchy.Groups(n).Groups(nGroup).Datasets)
+                    FD.names{m} = info.GroupHierarchy.Groups(n).Groups(nGroup).Datasets(m).Name;
+                    FD.freq(m) = double(info.GroupHierarchy.Groups(n).Groups(nGroup).Datasets(m).Attributes.Value);
+                end
+            end
         end
+        
     end
 end
 
-hdf_fielddata.names = names;
+if (numel(TD.names)>0)
-for n=1:numel(names)
+    for n=1:numel(TD.names)
-    hdf_fielddata.values{n} = double(hdf5read(file,names{n}));
+        hdf_fielddata.TD.values{n} = double(hdf5read(file,TD.names{n}));
+    end
 end
 
+if (numel(FD.names)>0)
+    Nr_freq = numel(FD.names)/2;
+    for n=1:Nr_freq
+        name = ['/FieldData/FD/f' int2str(n-1) '_real'];
+        ind = find(strcmp(FD.names,name));
+        hdf_fielddata.FD.values{n} = double(hdf5read(file,FD.names{ind}));
+        hdf_fielddata.FD.freq(n) = FD.freq(ind);
+        name = ['/FieldData/FD/f' int2str(n-1) '_imag'];
+        ind = find(strcmp(FD.names,name));
+        hdf_fielddata.FD.values{n} = hdf_fielddata.FD.values{n} + 1j*double(hdf5read(file,FD.names{ind}));
+    end
+end
 
 
 function hdf_fielddata = ReadHDF5FieldData_octave(file)