summaryrefslogtreecommitdiff
path: root/openEMS
diff options
context:
space:
mode:
Diffstat (limited to 'openEMS')
-rw-r--r--openEMS/CMakeLists.txt17
-rw-r--r--openEMS/Common/processmodematch.cpp12
-rw-r--r--openEMS/FDTD/excitation.cpp55
-rw-r--r--openEMS/FDTD/extensions/engine_ext_excitation.cpp12
-rw-r--r--openEMS/FDTD/extensions/engine_ext_tfsf.cpp4
-rw-r--r--openEMS/FDTD/openems_fdtd_mpi.cpp7
-rw-r--r--openEMS/main.cpp8
-rw-r--r--openEMS/matlab/DelayFidelity.m93
-rw-r--r--openEMS/matlab/RunOpenEMS.m2
-rw-r--r--openEMS/matlab/RunOpenEMS_MPI.m2
-rw-r--r--openEMS/matlab/Tutorials/CylindricalWave_CC.m2
-rw-r--r--openEMS/matlab/Tutorials/RadarUWBTutorial.m236
-rw-r--r--openEMS/matlab/Tutorials/Simple_Patch_Antenna.m53
-rw-r--r--openEMS/matlab/Tutorials/StripLine2MSL.m133
-rwxr-xr-xopenEMS/matlab/h5readatt_octave.cc21
-rw-r--r--openEMS/matlab/plotRefl.m145
-rw-r--r--openEMS/matlab/setup.m20
-rw-r--r--openEMS/nf2ff/CMakeLists.txt3
-rw-r--r--openEMS/openEMS_MPI.pro5
-rw-r--r--openEMS/openems.cpp9
-rw-r--r--openEMS/tools/hdf5_file_reader.cpp4
-rw-r--r--openEMS/tools/hdf5_file_writer.cpp14
22 files changed, 731 insertions, 126 deletions
diff --git a/openEMS/CMakeLists.txt b/openEMS/CMakeLists.txt
index 5434138..898f280 100644
--- a/openEMS/CMakeLists.txt
+++ b/openEMS/CMakeLists.txt
@@ -6,13 +6,13 @@ ELSE()
SET( CMAKE_BUILD_TYPE Release CACHE STRING "Set to either \"Release\" or \"Debug\"" )
ENDIF()
-PROJECT(openEMS CXX)
+PROJECT(openEMS CXX C)
cmake_minimum_required(VERSION 2.8)
# default
set(LIB_VERSION_MAJOR 0)
set(LIB_VERSION_MINOR 0)
-set(LIB_VERSION_PATCH 34)
+set(LIB_VERSION_PATCH 35)
set(LIB_VERSION_STRING ${LIB_VERSION_MAJOR}.${LIB_VERSION_MINOR}.${LIB_VERSION_PATCH})
set(VERSION "v${LIB_VERSION_STRING}")
@@ -21,7 +21,7 @@ IF(EXISTS ${PROJECT_SOURCE_DIR}/localConfig.cmake)
include(${PROJECT_SOURCE_DIR}/localConfig.cmake)
ENDIF()
-set(VERSION "v0.0.34")
+set(VERSION "v0.0.35")
# add git revision
IF(EXISTS ${PROJECT_SOURCE_DIR}/.git )
@@ -109,11 +109,8 @@ ADD_DEFINITIONS( -DTIXML_USE_STL )
# hdf5
find_package(HDF5 1.8 COMPONENTS C HL REQUIRED)
-INCLUDE_DIRECTORIES (${HDF5_INCLUDE_DIR})
-link_directories(${HDF5_LIBRARY_DIRS})
-
-# hdf5 compat
-ADD_DEFINITIONS( -DH5_USE_16_API )
+INCLUDE_DIRECTORIES (${HDF5_INCLUDE_DIRS})
+link_directories(${HDF5_LIBRARIES})
# boost
find_package(Boost 1.46 COMPONENTS
@@ -158,6 +155,8 @@ set(SOURCES
openems.cpp
)
+set(PUB_HEADERS openems.h)
+
# libs
ADD_SUBDIRECTORY( tools )
ADD_SUBDIRECTORY( FDTD )
@@ -181,6 +180,7 @@ TARGET_LINK_LIBRARIES( openEMS
${fparser_LIBRARIES}
tinyxml
${HDF5_LIBRARIES}
+ ${HDF5_HL_LIBRARIES}
${Boost_LIBRARIES}
${vtk_LIBS}
${MPI_LIBRARIES}
@@ -204,5 +204,6 @@ if (UNIX)
PERMISSIONS OWNER_WRITE OWNER_READ GROUP_READ WORLD_READ OWNER_EXECUTE GROUP_EXECUTE WORLD_EXECUTE)
endif()
endif ()
+INSTALL(FILES ${PUB_HEADERS} DESTINATION include/openEMS)
INSTALL( DIRECTORY matlab DESTINATION share/openEMS )
# TODO mpi, tarball, debug, release
diff --git a/openEMS/Common/processmodematch.cpp b/openEMS/Common/processmodematch.cpp
index 9c71c7a..620257d 100644
--- a/openEMS/Common/processmodematch.cpp
+++ b/openEMS/Common/processmodematch.cpp
@@ -135,22 +135,23 @@ void ProcessModeMatch::InitProcess()
m_ModeDist[n] = Create2DArray<double>(m_numLines);
}
+ bool dualMesh = m_ModeFieldType==1;
unsigned int pos[3] = {0,0,0};
double discLine[3] = {0,0,0};
double gridDelta = 1; // 1 -> mode-matching function is definied in drawing units...
double var[7];
pos[m_ny] = start[m_ny];
- discLine[m_ny] = Op->GetDiscLine(m_ny,pos[m_ny],m_dualMesh);
+ discLine[m_ny] = Op->GetDiscLine(m_ny,pos[m_ny],dualMesh);
double norm = 0;
double area = 0;
for (unsigned int posP = 0; posP<m_numLines[0]; ++posP)
{
pos[nP] = start[nP] + posP;
- discLine[nP] = Op->GetDiscLine(nP,pos[nP],m_dualMesh);
+ discLine[nP] = Op->GetDiscLine(nP,pos[nP],dualMesh);
for (unsigned int posPP = 0; posPP<m_numLines[1]; ++posPP)
{
pos[nPP] = start[nPP] + posPP;
- discLine[nPP] = Op->GetDiscLine(nPP,pos[nPP],m_dualMesh);
+ discLine[nPP] = Op->GetDiscLine(nPP,pos[nPP],dualMesh);
var[0] = discLine[0] * gridDelta; // x
var[1] = discLine[1] * gridDelta; // y
@@ -169,7 +170,7 @@ void ProcessModeMatch::InitProcess()
var[5] = sqrt(pow(discLine[0],2)+pow(discLine[2],2)) * gridDelta; // r
var[6] = asin(1)-atan(var[2]/var[3]); //theta (t)
}
- area = Op->GetNodeArea(m_ny,pos,m_dualMesh);
+ area = Op->GetNodeArea(m_ny,pos,dualMesh);
for (int n=0; n<2; ++n)
{
m_ModeDist[n][posP][posPP] = m_ModeParser[n]->Eval(var); //calc mode template
@@ -227,6 +228,7 @@ double* ProcessModeMatch::CalcMultipleIntegrals()
double field = 0;
double purity = 0;
double area = 0;
+ bool dualMesh = m_ModeFieldType==1;
int nP = (m_ny+1)%3;
int nPP = (m_ny+2)%3;
@@ -242,7 +244,7 @@ double* ProcessModeMatch::CalcMultipleIntegrals()
for (unsigned int posPP = 0; posPP<m_numLines[1]; ++posPP)
{
pos[nPP] = start[nPP] + posPP;
- area = Op->GetNodeArea(m_ny,pos,m_dualMesh);
+ area = Op->GetNodeArea(m_ny,pos,dualMesh);
if (m_ModeFieldType==0)
m_Eng_Interface->GetEField(pos,out);
if (m_ModeFieldType==1)
diff --git a/openEMS/FDTD/excitation.cpp b/openEMS/FDTD/excitation.cpp
index f095b53..4e1ba85 100644
--- a/openEMS/FDTD/excitation.cpp
+++ b/openEMS/FDTD/excitation.cpp
@@ -137,7 +137,7 @@ unsigned int Excitation::GetMaxExcitationTimestep() const
{
FDTD_FLOAT maxAmp=0;
unsigned int maxStep=0;
- for (unsigned int n=1; n<Length+1; ++n)
+ for (unsigned int n=0; n<Length; ++n)
{
if (fabs(Signal_volt[n])>maxAmp)
{
@@ -152,7 +152,7 @@ void Excitation::CalcGaussianPulsExcitation(double f0, double fc, int nTS)
{
if (dT==0) return;
- Length = (unsigned int)(2.0 * 9.0/(2.0*PI*fc) / dT);
+ Length = (unsigned int)ceil(2.0 * 9.0/(2.0*PI*fc) / dT);
if (Length>(unsigned int)nTS)
{
cerr << "Operator::CalcGaussianPulsExcitation: Requested excitation pusle would be " << Length << " timesteps or " << Length * dT << " s long. Cutting to max number of timesteps!" << endl;
@@ -160,13 +160,11 @@ void Excitation::CalcGaussianPulsExcitation(double f0, double fc, int nTS)
}
delete[] Signal_volt;
delete[] Signal_curr;
- Signal_volt = new FDTD_FLOAT[Length+1];
- Signal_curr = new FDTD_FLOAT[Length+1];
- Signal_volt[0]=0.0;
- Signal_curr[0]=0.0;
- for (unsigned int n=1; n<Length+1; ++n)
+ Signal_volt = new FDTD_FLOAT[Length];
+ Signal_curr = new FDTD_FLOAT[Length];
+ for (unsigned int n=0; n<Length; ++n)
{
- double t = (n-1)*dT;
+ double t = n*dT;
Signal_volt[n] = cos(2.0*PI*f0*(t-9.0/(2.0*PI*fc)))*exp(-1*pow(2.0*PI*fc*t/3.0-3,2));
t += 0.5*dT;
Signal_curr[n] = cos(2.0*PI*f0*(t-9.0/(2.0*PI*fc)))*exp(-1*pow(2.0*PI*fc*t/3.0-3,2));
@@ -182,12 +180,12 @@ void Excitation::CalcDiracPulsExcitation()
{
if (dT==0) return;
- Length = 1;
+ Length = 2;
// cerr << "Operator::CalcDiracPulsExcitation: Length of the excite signal: " << ExciteLength << " timesteps" << endl;
delete[] Signal_volt;
delete[] Signal_curr;
- Signal_volt = new FDTD_FLOAT[Length+1];
- Signal_curr = new FDTD_FLOAT[Length+1];
+ Signal_volt = new FDTD_FLOAT[Length];
+ Signal_curr = new FDTD_FLOAT[Length];
Signal_volt[0]=0.0;
Signal_volt[1]=1.0;
Signal_curr[0]=0.0;
@@ -203,11 +201,11 @@ void Excitation::CalcStepExcitation()
{
if (dT==0) return;
- Length = 1;
+ Length = 2;
delete[] Signal_volt;
delete[] Signal_curr;
- Signal_volt = new FDTD_FLOAT[Length+1];
- Signal_curr = new FDTD_FLOAT[Length+1];
+ Signal_volt = new FDTD_FLOAT[Length];
+ Signal_curr = new FDTD_FLOAT[Length];
Signal_volt[0]=1.0;
Signal_volt[1]=1.0;
Signal_curr[0]=1.0;
@@ -228,10 +226,8 @@ void Excitation::CalcCustomExcitation(double f0, int nTS, string signal)
// cerr << "Operator::CalcSinusExcitation: Length of the excite signal: " << ExciteLength << " timesteps" << endl;
delete[] Signal_volt;
delete[] Signal_curr;
- Signal_volt = new FDTD_FLOAT[Length+1];
- Signal_curr = new FDTD_FLOAT[Length+1];
- Signal_volt[0]=0.0;
- Signal_curr[0]=0.0;
+ Signal_volt = new FDTD_FLOAT[Length];
+ Signal_curr = new FDTD_FLOAT[Length];
FunctionParser fParse;
fParse.AddConstant("pi", 3.14159265358979323846);
fParse.AddConstant("e", 2.71828182845904523536);
@@ -242,9 +238,9 @@ void Excitation::CalcCustomExcitation(double f0, int nTS, string signal)
exit(1);
}
double vars[1];
- for (unsigned int n=1; n<Length+1; ++n)
+ for (unsigned int n=0; n<Length; ++n)
{
- vars[0] = (n-1)*dT;
+ vars[0] = n*dT;
Signal_volt[n] = fParse.Eval(vars);
vars[0] += 0.5*dT;
Signal_curr[n] = fParse.Eval(vars);
@@ -260,17 +256,16 @@ void Excitation::CalcSinusExcitation(double f0, int nTS)
if (dT==0) return;
if (nTS<=0) return;
- Length = (unsigned int)(2.0/f0/dT);
- //cerr << "Operator::CalcSinusExcitation: Length of the excite signal: " << Length << " timesteps " << Length*dT << "s" << endl;
+ Length = (unsigned int)round(2.0/f0/dT);
delete[] Signal_volt;
delete[] Signal_curr;
- Signal_volt = new FDTD_FLOAT[Length+1];
- Signal_curr = new FDTD_FLOAT[Length+1];
+ Signal_volt = new FDTD_FLOAT[Length];
+ Signal_curr = new FDTD_FLOAT[Length];
Signal_volt[0]=0.0;
Signal_curr[0]=0.0;
- for (unsigned int n=1; n<Length+1; ++n)
+ for (unsigned int n=1; n<Length; ++n)
{
- double t = (n-1)*dT;
+ double t = n*dT;
Signal_volt[n] = sin(2.0*PI*f0*t);
t += 0.5*dT;
Signal_curr[n] = sin(2.0*PI*f0*t);
@@ -286,8 +281,8 @@ void Excitation::DumpVoltageExcite(string filename)
file.open( filename.c_str() );
if (file.fail())
return;
- for (unsigned int n=1; n<Length+1; ++n)
- file << (n-1)*dT << "\t" << Signal_volt[n] << "\n";
+ for (unsigned int n=0; n<Length; ++n)
+ file << n*dT << "\t" << Signal_volt[n] << "\n";
file.close();
}
@@ -297,8 +292,8 @@ void Excitation::DumpCurrentExcite(string filename)
file.open( filename.c_str() );
if (file.fail())
return;
- for (unsigned int n=1; n<Length+1; ++n)
- file << (n-1)*dT + 0.5*dT << "\t" << Signal_curr[n] << "\n";
+ for (unsigned int n=0; n<Length; ++n)
+ file << n*dT + 0.5*dT << "\t" << Signal_curr[n] << "\n";
file.close();
}
diff --git a/openEMS/FDTD/extensions/engine_ext_excitation.cpp b/openEMS/FDTD/extensions/engine_ext_excitation.cpp
index b5bd8bd..402193a 100644
--- a/openEMS/FDTD/extensions/engine_ext_excitation.cpp
+++ b/openEMS/FDTD/extensions/engine_ext_excitation.cpp
@@ -54,7 +54,7 @@ void Engine_Ext_Excitation::Apply2Voltages()
exc_pos = numTS - (int)m_Op_Exc->Volt_delay[n];
exc_pos *= (exc_pos>0);
exc_pos %= p;
- exc_pos *= (exc_pos<=(int)length);
+ exc_pos *= (exc_pos<(int)length);
ny = m_Op_Exc->Volt_dir[n];
pos[0]=m_Op_Exc->Volt_index[0][n];
pos[1]=m_Op_Exc->Volt_index[1][n];
@@ -71,7 +71,7 @@ void Engine_Ext_Excitation::Apply2Voltages()
exc_pos = numTS - (int)m_Op_Exc->Volt_delay[n];
exc_pos *= (exc_pos>0);
exc_pos %= p;
- exc_pos *= (exc_pos<=(int)length);
+ exc_pos *= (exc_pos<(int)length);
ny = m_Op_Exc->Volt_dir[n];
pos[0]=m_Op_Exc->Volt_index[0][n];
pos[1]=m_Op_Exc->Volt_index[1][n];
@@ -87,7 +87,7 @@ void Engine_Ext_Excitation::Apply2Voltages()
exc_pos = numTS - (int)m_Op_Exc->Volt_delay[n];
exc_pos *= (exc_pos>0);
exc_pos %= p;
- exc_pos *= (exc_pos<=(int)length);
+ exc_pos *= (exc_pos<(int)length);
ny = m_Op_Exc->Volt_dir[n];
pos[0]=m_Op_Exc->Volt_index[0][n];
pos[1]=m_Op_Exc->Volt_index[1][n];
@@ -124,7 +124,7 @@ void Engine_Ext_Excitation::Apply2Current()
exc_pos = numTS - (int)m_Op_Exc->Curr_delay[n];
exc_pos *= (exc_pos>0);
exc_pos %= p;
- exc_pos *= (exc_pos<=(int)length);
+ exc_pos *= (exc_pos<(int)length);
ny = m_Op_Exc->Curr_dir[n];
pos[0]=m_Op_Exc->Curr_index[0][n];
pos[1]=m_Op_Exc->Curr_index[1][n];
@@ -141,7 +141,7 @@ void Engine_Ext_Excitation::Apply2Current()
exc_pos = numTS - (int)m_Op_Exc->Curr_delay[n];
exc_pos *= (exc_pos>0);
exc_pos %= p;
- exc_pos *= (exc_pos<=(int)length);
+ exc_pos *= (exc_pos<(int)length);
ny = m_Op_Exc->Curr_dir[n];
pos[0]=m_Op_Exc->Curr_index[0][n];
pos[1]=m_Op_Exc->Curr_index[1][n];
@@ -157,7 +157,7 @@ void Engine_Ext_Excitation::Apply2Current()
exc_pos = numTS - (int)m_Op_Exc->Curr_delay[n];
exc_pos *= (exc_pos>0);
exc_pos %= p;
- exc_pos *= (exc_pos<=(int)length);
+ exc_pos *= (exc_pos<(int)length);
ny = m_Op_Exc->Curr_dir[n];
pos[0]=m_Op_Exc->Curr_index[0][n];
pos[1]=m_Op_Exc->Curr_index[1][n];
diff --git a/openEMS/FDTD/extensions/engine_ext_tfsf.cpp b/openEMS/FDTD/extensions/engine_ext_tfsf.cpp
index cd3b107..9942656 100644
--- a/openEMS/FDTD/extensions/engine_ext_tfsf.cpp
+++ b/openEMS/FDTD/extensions/engine_ext_tfsf.cpp
@@ -44,7 +44,7 @@ void Engine_Ext_TFSF::DoPostVoltageUpdates()
{
if ( numTS < n )
m_DelayLookup[n]=0;
- else if ((numTS-n > length) && (p==0))
+ else if ((numTS-n >= length) && (p==0))
m_DelayLookup[n]=0;
else
m_DelayLookup[n] = numTS - n;
@@ -134,7 +134,7 @@ void Engine_Ext_TFSF::DoPostCurrentUpdates()
{
if ( numTS < n )
m_DelayLookup[n]=0;
- else if ((numTS-n > length) && (p==0))
+ else if ((numTS-n >= length) && (p==0))
m_DelayLookup[n]=0;
else
m_DelayLookup[n] = numTS - n;
diff --git a/openEMS/FDTD/openems_fdtd_mpi.cpp b/openEMS/FDTD/openems_fdtd_mpi.cpp
index eaab59a..964bdd4 100644
--- a/openEMS/FDTD/openems_fdtd_mpi.cpp
+++ b/openEMS/FDTD/openems_fdtd_mpi.cpp
@@ -210,7 +210,7 @@ bool openEMS_FDTD_MPI::SetupMPI()
if (numProcs!=m_NumProc)
{
if (m_MyID==0)
- cerr << "openEMS_FDTD_MPI::SetupMPI: Error: Requested splits require " << numProcs << " processes, but only " << m_NumProc << " were found! Exit! " << endl;
+ cerr << "openEMS_FDTD_MPI::SetupMPI: Error: Requested splits require " << numProcs << " processes, but " << m_NumProc << " were found! Exit! " << endl;
exit(10);
}
@@ -256,8 +256,8 @@ bool openEMS_FDTD_MPI::SetupMPI()
grid->AddDiscLine(2, m_Original_Grid->GetLine(2,n) );
m_MPI_Op->SetSplitPos(0,m_SplitNumber[0].at(i));
- m_MPI_Op->SetSplitPos(1,m_SplitNumber[1].at(i));
- m_MPI_Op->SetSplitPos(2,m_SplitNumber[2].at(i));
+ m_MPI_Op->SetSplitPos(1,m_SplitNumber[1].at(j));
+ m_MPI_Op->SetSplitPos(2,m_SplitNumber[2].at(k));
if (i>0)
m_MPI_Op->SetNeighborDown(0,procTable[i-1][j][k]);
@@ -398,6 +398,7 @@ bool openEMS_FDTD_MPI::SetupProcessing()
{
//type is integral processing --> disable! Needs to be fixed!
cerr << "openEMS_FDTD_MPI::SetupProcessing(): Warning: Processing: " << proc->GetName() << " occures multiple times and is being deactivated..." << endl;
+ cerr << "openEMS_FDTD_MPI::SetupProcessing(): Note: Processing: Make sure that there are no splits inside probes or sources." << endl;
deactivate = true;
rename = false;
}
diff --git a/openEMS/main.cpp b/openEMS/main.cpp
index 1d03481..8a8000f 100644
--- a/openEMS/main.cpp
+++ b/openEMS/main.cpp
@@ -68,8 +68,12 @@ int main(int argc, char *argv[])
}
int EC = FDTD.ParseFDTDSetup(argv[1]);
+ if(!EC) {
+ cerr << "openEMS - ParseFDTDSetup failed." << endl;
+ exit(1);
+ }
EC = FDTD.SetupFDTD();
- if (EC) return EC;
+ if (EC) exit(EC);
FDTD.RunFDTD();
#ifdef MPI_SUPPORT
@@ -77,5 +81,5 @@ int main(int argc, char *argv[])
MPI::Finalize();
#endif
- return 0;
+ exit(0);
}
diff --git a/openEMS/matlab/DelayFidelity.m b/openEMS/matlab/DelayFidelity.m
new file mode 100644
index 0000000..3cc0ae2
--- /dev/null
+++ b/openEMS/matlab/DelayFidelity.m
@@ -0,0 +1,93 @@
+function [delay, fidelity, nf2ff_out] = DelayFidelity(nf2ff, port, path, weight_theta, weight_phi, theta, phi, f_0, f_c, varargin)
+% [delay, fidelity] = DelayFidelity(nf2ff, port, path, theta, phi, f_lo, f_hi, varargin)
+%
+%
+% This function calculates the time delay from the source port to the phase center of the antenna and the fidelity.
+% The fidelity is the similarity between the excitation pulse and the radiated pulse (normalized scalar product).
+% The resolution of the delay will be equal to or better than ((f_0 + f_c)*Oversampling)^-1 when using Gaussian excitation.
+% Oversampling is an input parameter to InitFDTD. The rows of delay and fidelity correspond to theta and the columns to phi.
+%
+% input:
+% nf2ff: return value of CreateNF2FFBox.
+% port: return value of AddLumpedPort
+% path: path of the simulation results.
+% weight_theta: weight if the E_theta component
+% weight_phi: eight of the E_phi component
+% -> with both (possibly complex) parameters any polarization can be examined
+% theta: theta values to be simulated
+% phi: phi values to be simulated
+% f_0: center frequency of SetGaussExcite
+% f_c: cutoff frequency of SetGaussExcite
+%
+% variable input:
+% 'Center': phase center of the antenna for CalcNF2FF
+% 'Radius': radius for CalcNF2FF
+% 'Mode': mode CalcNF2FF
+%
+% example:
+% theta = [-180:10:180] * pi / 180;
+% phi = [0, 90] * pi / 180;
+% [delay, fidelity] = DelayFidelity2(nf2ff, port, Sim_Path, sin(tilt), cos(tilt), theta, phi, f_0, f_c, 'Mode', 1);
+% figure
+% polar(theta.', delay(:,1) * 3e11); % delay in mm
+% figure
+% polar(theta', (fidelity(:,1)-0.95)/0.05); % last 5 percent of fidelity
+%
+% Author: Georg Michel
+
+C0 = 299792458;
+center = [0, 0, 0];
+radius = 1;
+nf2ff_mode = 0;
+
+for n=1:2:numel(varargin)
+ if (strcmp(varargin{n},'Center')==1);
+ center = varargin{n+1};
+ elseif (strcmp(varargin{n},'Radius')==1);
+ radius = varargin{n+1};
+ elseif (strcmp(varargin{n},'Mode')==1);
+ nf2ff_mode = varargin{n+1};
+ end
+end
+
+
+port_ut = load(fullfile(path, port.U_filename));
+port_it = load(fullfile(path, port.I_filename));
+dt = port_ut(2,1) - port_ut(1,1);
+fftsize = 2^(nextpow2(size(port_ut)(1)) + 1);
+df = 1 / (dt * fftsize);
+uport = fft(port_ut(:, 2), fftsize)(1:fftsize/2+1);
+iport = fft(port_it(:, 2), fftsize)(1:fftsize/2+1);
+fport = df * (0:fftsize/2);
+f_ind = find(fport > (f_0 - f_c ) & fport < (f_0 + f_c));
+disp(["frequencies: ", num2str(numel(f_ind))]);
+exc_f = uport.' + iport.' * port.Feed_R; %excitation in freq domain
+exc_f(!f_ind) = 0;
+exc_f /= sqrt(exc_f * exc_f'); % normalization (transposing also conjugates)
+
+nf2ff = CalcNF2FF(nf2ff, path, fport(f_ind), theta, phi, ...
+ 'Center', center, 'Radius', radius, 'Mode', nf2ff_mode);
+radfield = weight_theta * cell2mat(nf2ff.E_theta) + weight_phi * cell2mat(nf2ff.E_phi); % rows: theta(f1), columns: phi(f1), phi(f2), ...phi(fn)
+radfield = reshape(radfield, [length(nf2ff.theta), length(nf2ff.phi), length(nf2ff.freq)]);
+correction = reshape(exp(-2i*pi*nf2ff.r/C0*nf2ff.freq), 1,1,numel(nf2ff.freq)); %dimensions: theta, phi, frequencies
+radfield = radfield./correction; % correct for radius delay
+% normalize radfield
+radnorm = sqrt(dot(radfield, radfield, 3));
+radfield ./= radnorm;
+
+%initialize radiated field in fully populated frequency domain
+rad_f = zeros([numel(nf2ff.theta), numel(nf2ff.phi), numel(fport)]);
+rad_f(:, :, f_ind) = radfield; % assign selected frequencies
+exc_f = reshape(exc_f, [1,1,numel(exc_f)]); %make exc_f confomant with rad_f
+
+cr_f = rad_f .* conj(exc_f); % calculate cross correlation
+% calculate the cross correlation in time domain (analytic signal)
+cr = ifft(cr_f(:, :, 1:end-1), [], 3) * (numel(fport) -1); % twice the FFT normalization (sqrt^2) because product of two normalized functions
+%search for the maxiumum of the envelope
+[fidelity, delay_ind] = max(abs(cr), [], 3);
+delay = (delay_ind - 1) * dt * 2; % double time step because of single-sided FFT
+nf2ff_out = nf2ff; %possibly needed for plotting the far field and other things
+disp(["DelayFidelity: delay resolution = ", num2str(dt*2e9), "ns"]);
+return;
+
+
diff --git a/openEMS/matlab/RunOpenEMS.m b/openEMS/matlab/RunOpenEMS.m
index a4947e3..3894adc 100644
--- a/openEMS/matlab/RunOpenEMS.m
+++ b/openEMS/matlab/RunOpenEMS.m
@@ -21,7 +21,7 @@ function RunOpenEMS(Sim_Path, Sim_File, opts, Settings)
% --engine=fastest fastest available engine (default)
% --engine=basic basic FDTD engine
% --engine=sse engine using sse vector extensions
-% --engine=sse_compressed engine using compressed operator + sse vector extensions
+% --engine=sse-compressed engine using compressed operator + sse vector extensions
% --engine=MPI engine using compressed operator + sse vector extensions + MPI parallel processing
% --engine=multithreaded engine using compressed operator + sse vector extensions + MPI + multithreading
% --numThreads=<n> Force use n threads for multithreaded engine
diff --git a/openEMS/matlab/RunOpenEMS_MPI.m b/openEMS/matlab/RunOpenEMS_MPI.m
index 495f7e5..779b038 100644
--- a/openEMS/matlab/RunOpenEMS_MPI.m
+++ b/openEMS/matlab/RunOpenEMS_MPI.m
@@ -87,7 +87,7 @@ end
if isfield(Settings.MPI,'Hosts')
disp(['Running remote openEMS_MPI in working dir: ' work_path]);
- [status] = system(['mpiexec -host ' HostList ' -n ' int2str(NrProc) ' -wdir ' work_path ' ' Settings.MPI.Binary ' ' Sim_File ' ' opts ' ' append_unix]);
+ [status] = system(['mpiexec ' Settings.MPI.GlobalArgs ' -host ' HostList ' -n ' int2str(NrProc) ' -wdir ' work_path ' ' Settings.MPI.Binary ' ' Sim_File ' ' opts ' ' append_unix]);
else
disp('Running local openEMS_MPI');
[status] = system(['mpiexec ' Settings.MPI.GlobalArgs ' -n ' int2str(NrProc) ' ' Settings.MPI.Binary ' ' Sim_File ' ' opts ' ' append_unix]);
diff --git a/openEMS/matlab/Tutorials/CylindricalWave_CC.m b/openEMS/matlab/Tutorials/CylindricalWave_CC.m
index d55c470..b81d2d4 100644
--- a/openEMS/matlab/Tutorials/CylindricalWave_CC.m
+++ b/openEMS/matlab/Tutorials/CylindricalWave_CC.m
@@ -28,7 +28,7 @@ exite_offset = 1300;
excite_angle = 45;
%% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-FDTD = InitFDTD(100000,1e-4,'CoordSystem',1,'MultiGrid',split);
+FDTD = InitFDTD('NrTS', 100000, 'EndCriteria', 1e-4, 'CoordSystem', 1, 'MultiGrid', split);
FDTD = SetGaussExcite(FDTD,f0,f0/2);
BC = [0 3 0 0 0 0]; % pml in positive r-direction
FDTD = SetBoundaryCond(FDTD,BC);
diff --git a/openEMS/matlab/Tutorials/RadarUWBTutorial.m b/openEMS/matlab/Tutorials/RadarUWBTutorial.m
new file mode 100644
index 0000000..6fd2f5b
--- /dev/null
+++ b/openEMS/matlab/Tutorials/RadarUWBTutorial.m
@@ -0,0 +1,236 @@
+% Tutorial on time delay and signal integrity for radar
+% and UWB applications
+%
+% Tested with
+% - Octave 4.0
+% - openEMS v0.0.35
+%
+% Author: Georg Michel, 2016
+
+ clear;
+ close all;
+
+physical_constants;
+
+% --- start of configuration section ---
+
+% In radar and ultrawideband applications it is important to know the
+% delay and fidelity of RF pulses. The delay is the retardation of the
+% signal from the source to the phase center of the antenna. It is
+% composed out of linear delay, dispersion and minimum-phase
+% delay. Dispersion due to waveguides or frequency-dependent
+% permittivity and minimum-phase delay due to resonances will degrade
+% the fidelity which is the normalized similarity between excitation and
+% radiated signal. In this tutorial you can examine the performance of a
+% simple ultrawideband (UWB) monopole. The delay and fidelity of this
+% antenna are calculated and plotted. You can compare these properties
+% in different channels.
+%
+% The Gaussian excitation is set to the same 3dB bandwidth as the
+% channels of the IEEE 802.15.4 UWB PHY. One exeption is channel4twice
+% which has the double bandwidth of channel 4. It can be seen that the
+% delay is larger and the fidelity is smaller in the vicinity of the
+% (undesired) resonances of the antenna. Note that for a real UWB system
+% the total delay and fidelity result from both the transmitting and
+% receiving antenna or twice the delay and the square of the fidelity
+% for monostatic radar.
+%
+% The resolution of the delay will depend on the 'Oversampling'
+% parameter to InitFDTD. See the description of DelayFidelity
+%
+% In the configuration section below you can uncomment the respective
+% parameter settings. As an exercise, you can examine how the permittivity
+% of the substrate influences gain, delay and fidelity.
+
+
+%suffix = "channel1";
+%f_0 = 3.5e9; % center frequency of the channel
+%f_c = 0.25e9 / 0.3925; % 3dB bandwidth is 0.3925 times 20dB bandwidth for Gaussian excitation
+
+%suffix = "channel2";
+%f_0 = 4.0e9; % center frequency of the channel
+%f_c = 0.25e9 / 0.3925;
+
+%suffix = "channel3";
+%f_0 = 4.5e9; % center frequency of the channel
+%f_c = 0.25e9 / 0.3925;
+
+suffix = "channel4";
+f_0 = 4.0e9; % center frequency of the channel
+f_c = 0.5e9 / 0.3925;
+
+%suffix = "channel5";
+%f_0 = 6.5e9; % center frequency of the channel
+%f_c = 0.25e9 / 0.3925;
+
+%suffix = "channel7";
+%f_0 = 6.5e9; % center frequency of the channel
+%f_c = 0.5e9 / 0.3925;
+
+%suffix = "channel4twice"; % this is just to demonstrate the degradation of the fidelity with increasing bandwidth
+%f_0 = 4.0e9; % center frequency of the channel
+%f_c = 1e9 / 0.3925;
+
+tilt = 45 * pi / 180; % polarization tilt angle against co-polarization (90DEG is cross polarized)
+
+% --- end of configuration section ---
+
+% path and filename setup
+Sim_Path = 'tmp';
+Sim_CSX = 'uwb.xml';
+
+% properties of the substrate
+substrate.epsR = 4; % FR4
+substrate.height = 0.707;
+substrate.cells = 3; % thickness in cells
+
+% size of the monopole and the gap to the ground plane
+gap = 0.62; % 0.5
+patchsize = 14;
+
+% we will use millimeters
+unit = 1e-3;
+
+% set the resolution for the finer structures, e.g. the antenna gap
+fineResolution = C0 / (f_0 + f_c) / sqrt(substrate.epsR) / unit / 40;
+% set the resolution for the coarser structures, e.g. the surrounding air
+coarseResolution = C0/(f_0 + f_c) / unit / 20;
+
+
+% initialize the CSX structure
+CSX = InitCSX();
+
+% add the properties which are used to model the antenna
+CSX = AddMetal(CSX, 'Ground' );
+CSX = AddMetal(CSX, 'Patch');
+CSX = AddMetal(CSX, 'Line');
+CSX = AddMaterial(CSX, 'Substrate' );
+CSX = SetMaterialProperty(CSX, 'Substrate', 'Epsilon', substrate.epsR);
+
+% define the supstrate and sheet-like primitives for the properties
+CSX = AddBox(CSX, 'Substrate', 1, [-16, -16, -substrate.height], [16, 18, 0]);
+CSX = AddBox(CSX, 'Ground', 2, [-16, -16, -substrate.height], [16, 0, -substrate.height]);
+CSX = AddBox(CSX, 'Line', 2, [-1.15, -16, 0], [1.15, gap, 0]);
+CSX = AddBox(CSX, 'Patch', 2, [-patchsize/2, gap, 0], [patchsize/2, gap + patchsize, 0]);
+
+% setup a mesh
+mesh.x = [];
+mesh.y = [];
+
+% two mesh lines for the metal coatings of teh substrate
+mesh.z = linspace(-substrate.height, 0, substrate.cells +1);
+
+% find optimal mesh lines for the patch and ground, not yes the microstrip line
+mesh = DetectEdges(CSX, mesh, 'SetProperty',{'Patch', 'Ground'}, '2D_Metal_Edge_Res', fineResolution/2);
+
+%replace gap mesh lines which are too close by a single mesh line
+tooclose = find (diff(mesh.y) < fineResolution/4);
+if ~isempty(tooclose)
+ mesh.y(tooclose) = (mesh.y(tooclose) + mesh.y(tooclose+1))/2;
+ mesh.y(tooclose + 1) = [];
+endif
+
+% store the microstrip edges in a temporary variable
+meshline = DetectEdges(CSX, [], 'SetProperty', 'Line', '2D_Metal_Edge_Res', fineResolution/2);
+% as well as the edges of the substrate (without 1/3 - 2/3 rule)
+meshsubstrate = DetectEdges(CSX, [], 'SetProperty', 'Substrate');
+% add only the x mesh lines of the microstrip
+mesh.x = [mesh.x meshline.x];
+% and only the top of the substrate, the other edges are covered by the ground plane
+mesh.y = [mesh.y, meshsubstrate.y(end)]; % top of substrate
+
+% for now we have only the edges, now calculate mesh lines inbetween
+mesh = SmoothMesh(mesh, fineResolution);
+
+% add the outer boundary
+mesh.x = [mesh.x -60, 60];
+mesh.y = [mesh.y, -60, 65];
+mesh.z = [mesh.z, -46, 45];
+
+% add coarse mesh lines for the free space
+mesh = SmoothMesh(mesh, coarseResolution);
+
+% define the grid
+CSX = DefineRectGrid( CSX, unit, mesh);
+% and the feeding port
+[CSX, port] = AddLumpedPort( CSX, 999, 1, 50, [-1.15, meshline.y(2), -substrate.height], [1.15, meshline.y(2), 0], [0 0 1], true);
+
+%setup a NF2FF box for the calculation of the far field
+start = [mesh.x(10) mesh.y(10) mesh.z(10)];
+stop = [mesh.x(end-9) mesh.y(end-9) mesh.z(end-9)];
+[CSX nf2ff] = CreateNF2FFBox(CSX, 'nf2ff', start, stop);
+
+% initialize the FDTD structure with excitation and open boundary conditions
+FDTD = InitFDTD( 'NrTs', 30000, 'EndCriteria', 1e-5, 'OverSampling', 20);
+FDTD = SetGaussExcite(FDTD, f_0, f_c );
+BC = {'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8'};
+FDTD = SetBoundaryCond(FDTD, BC );
+
+
+% remove old data, show structure, calculate new data
+[status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory
+[status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder
+
+% write the data to the working directory
+WriteOpenEMS([Sim_Path '/' Sim_CSX], FDTD, CSX);
+% show the geometry for checking
+CSXGeomPlot([Sim_Path '/' Sim_CSX]);
+% run the simulation
+RunOpenEMS( Sim_Path, Sim_CSX);
+
+% plot amplitude and phase of the reflection coefficient
+freq = linspace(f_0-f_c, f_0+f_c, 200);
+port = calcPort(port, Sim_Path, freq);
+s11 = port.uf.ref ./ port.uf.inc;
+s11phase = unwrap(arg(s11));
+figure %("visible", "off"); % use this to plot only into files at the end of this script
+ax = plotyy( freq/1e6, 20*log10(abs(s11)), freq/1e6, s11phase);
+grid on
+title( ['reflection coefficient ', suffix, ' S_{11}']);
+xlabel( 'frequency f / MHz' );
+ylabel( ax(1), 'reflection coefficient |S_{11}|' );
+ylabel(ax(2), 'S_{11} phase (rad)');
+
+% define an azimuthal trace around the monopole
+phi = [0] * pi / 180;
+theta = [-180:10:180] * pi / 180;
+
+% calculate the delay, the fidelity and the farfield
+[delay, fidelity, nf2ff] = DelayFidelity(nf2ff, port, Sim_Path, sin(tilt), cos(tilt), theta, phi, f_0, f_c, 'Mode', 1);
+
+%plot the gain at (close to) f_0
+f_0_nearest_ind = nthargout(2, @min, abs(nf2ff.freq -f_0));
+%turn directivity into gain
+nf2ff.Dmax(f_0_nearest_ind) *= nf2ff.Prad(f_0_nearest_ind) / calcPort(port, Sim_Path, nf2ff.freq(f_0_nearest_ind)).P_inc;
+figure %("visible", "off");
+polarFF(nf2ff, 'xaxis', 'theta', 'freq_index', f_0_nearest_ind, 'logscale', [-4, 4]);
+title(["gain ", suffix, " / dBi"]);
+
+
+% We trick polarFF into plotting the delay in mm because
+% the axes of the vanilla polar plot can not be scaled.
+plotvar = delay * C0 * 1000;
+maxplot = 80;
+minplot = 30;
+nf2ff.Dmax(1) = 10^(max(plotvar)/10);
+nf2ff.E_norm{1} = 10.^(plotvar/20);
+figure %("visible", "off");
+polarFF(nf2ff, 'xaxis', 'theta', 'logscale', [minplot, maxplot]);
+title(["delay ", suffix, " / mm"]);
+
+% The same for the fidelity in percent.
+plotvar = fidelity * 100;
+maxplot = 100;
+minplot = 98;
+nf2ff.Dmax(1) = 10^(max(plotvar)/10);
+nf2ff.E_norm{1} = 10.^(plotvar/20);
+figure %("visible", "off");
+polarFF(nf2ff, 'xaxis', 'theta', 'logscale', [minplot, maxplot]);
+title(["fidelity ", suffix, " / %"]);
+
+% save the plots in order to compare them afer simulating the different channels
+print(1, ["s11_", suffix, ".png"]);
+print(2, ["farfield_", suffix, ".png"]);
+print(3, ["delay_mm_", suffix, ".png"]);
+print(4, ["fidelity_", suffix, ".png"]);
+return; \ No newline at end of file
diff --git a/openEMS/matlab/Tutorials/Simple_Patch_Antenna.m b/openEMS/matlab/Tutorials/Simple_Patch_Antenna.m
index ee2a8f0..1e24c4c 100644
--- a/openEMS/matlab/Tutorials/Simple_Patch_Antenna.m
+++ b/openEMS/matlab/Tutorials/Simple_Patch_Antenna.m
@@ -1,20 +1,20 @@
-%
-% Tutorials / simple patch antenna
+%% Simple Patch Antenna Tutorial
%
% Describtion at:
-% http://openems.de/index.php/Tutorial:_Simple_Patch_Antenna
+% <http://openems.de/index.php/Tutorial:_Simple_Patch_Antenna>
%
% Tested with
% - Matlab 2013a / Octave 4.0
-% - openEMS v0.0.33
+% - openEMS v0.0.35
%
-% (C) 2010-2015 Thorsten Liebig <thorsten.liebig@uni-due.de>
+% (C) 2010-2017 Thorsten Liebig <thorsten.liebig@uni-due.de>
+%%
close all
clear
clc
-%% setup the simulation
+%% Setup the Simulation
physical_constants;
unit = 1e-3; % all length in mm
@@ -38,7 +38,7 @@ feed.R = 50; %feed resistance
% size of the simulation box
SimBox = [200 200 150];
-%% setup FDTD parameter & excitation function
+%% Setup FDTD Parameter & Excitation Function
f0 = 2e9; % center frequency
fc = 1e9; % 20 dB corner frequency
FDTD = InitFDTD( 'NrTs', 30000 );
@@ -46,7 +46,7 @@ FDTD = SetGaussExcite( FDTD, f0, fc );
BC = {'MUR' 'MUR' 'MUR' 'MUR' 'MUR' 'MUR'}; % boundary conditions
FDTD = SetBoundaryCond( FDTD, BC );
-%% setup CSXCAD geometry & mesh
+%% Setup CSXCAD Geometry & Mesh
CSX = InitCSX();
%initialize the mesh with the "air-box" dimensions
@@ -54,13 +54,13 @@ mesh.x = [-SimBox(1)/2 SimBox(1)/2];
mesh.y = [-SimBox(2)/2 SimBox(2)/2];
mesh.z = [-SimBox(3)/3 SimBox(3)*2/3];
-%% create patch
+% Create Patch
CSX = AddMetal( CSX, 'patch' ); % create a perfect electric conductor (PEC)
start = [-patch.width/2 -patch.length/2 substrate.thickness];
stop = [ patch.width/2 patch.length/2 substrate.thickness];
CSX = AddBox(CSX,'patch',10,start,stop); % add a box-primitive to the metal property 'patch'
-%% create substrate
+% Create Substrate
CSX = AddMaterial( CSX, 'substrate' );
CSX = SetMaterialProperty( CSX, 'substrate', 'Epsilon', substrate.epsR, 'Kappa', substrate.kappa );
start = [-substrate.width/2 -substrate.length/2 0];
@@ -70,18 +70,18 @@ CSX = AddBox( CSX, 'substrate', 0, start, stop );
% add extra cells to discretize the substrate thickness
mesh.z = [linspace(0,substrate.thickness,substrate.cells+1) mesh.z];
-%% create ground (same size as substrate)
+% Create Ground same size as substrate
CSX = AddMetal( CSX, 'gnd' ); % create a perfect electric conductor (PEC)
start(3)=0;
stop(3) =0;
CSX = AddBox(CSX,'gnd',10,start,stop);
-%% apply the excitation & resist as a current source
+% Apply the Excitation & Resist as a Current Source
start = [feed.pos 0 0];
stop = [feed.pos 0 substrate.thickness];
[CSX port] = AddLumpedPort(CSX, 5 ,1 ,feed.R, start, stop, [0 0 1], true);
-%% finalize the mesh
+% Finalize the Mesh
% detect all edges except of the patch
mesh = DetectEdges(CSX, mesh,'ExcludeProperty','patch');
% detect and set a special 2D metal edge mesh for the patch
@@ -90,35 +90,41 @@ mesh = DetectEdges(CSX, mesh,'SetProperty','patch','2D_Metal_Edge_Res', c0/(f0+f
mesh = SmoothMesh(mesh, c0/(f0+fc)/unit/20);
CSX = DefineRectGrid(CSX, unit, mesh);
-%% add a nf2ff calc box; size is 3 cells away from MUR boundary condition
+CSX = AddDump(CSX,'Hf', 'DumpType', 11, 'Frequency',[2.4e9]);
+CSX = AddBox(CSX,'Hf',10,[-substrate.width -substrate.length -10*substrate.thickness],[substrate.width +substrate.length 10*substrate.thickness]); %assign box
+
+% add a nf2ff calc box; size is 3 cells away from MUR boundary condition
start = [mesh.x(4) mesh.y(4) mesh.z(4)];
stop = [mesh.x(end-3) mesh.y(end-3) mesh.z(end-3)];
[CSX nf2ff] = CreateNF2FFBox(CSX, 'nf2ff', start, stop);
-%% prepare simulation folder
+%% Prepare and Run Simulation
Sim_Path = 'tmp_Patch_Ant';
Sim_CSX = 'patch_ant.xml';
+% create an empty working directory
[status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory
[status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder
-%% write openEMS compatible xml-file
+% write openEMS compatible xml-file
WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX );
-%% show the structure
+% show the structure
CSXGeomPlot( [Sim_Path '/' Sim_CSX] );
-%% run openEMS
+% run openEMS
RunOpenEMS( Sim_Path, Sim_CSX);
-%% postprocessing & do the plots
+%% Postprocessing & Plots
freq = linspace( max([1e9,f0-fc]), f0+fc, 501 );
port = calcPort(port, Sim_Path, freq);
-Zin = port.uf.tot ./ port.if.tot;
-s11 = port.uf.ref ./ port.uf.inc;
+%% Smith chart port reflection
+plotRefl(port, 'threshold', -10)
+title( 'reflection coefficient' );
% plot feed point impedance
+Zin = port.uf.tot ./ port.if.tot;
figure
plot( freq/1e6, real(Zin), 'k-', 'Linewidth', 2 );
hold on
@@ -130,6 +136,7 @@ ylabel( 'impedance Z_{in} / Ohm' );
legend( 'real', 'imag' );
% plot reflection coefficient S11
+s11 = port.uf.ref ./ port.uf.inc;
figure
plot( freq/1e6, 20*log10(abs(s11)), 'k-', 'Linewidth', 2 );
grid on
@@ -139,7 +146,7 @@ ylabel( 'reflection coefficient |S_{11}|' );
drawnow
-%% NFFF contour plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% NFFF Plots
%find resonance frequncy from s11
f_res_ind = find(s11==min(s11));
f_res = freq(f_res_ind);
@@ -166,7 +173,7 @@ plotFFdB(nf2ff,'xaxis','theta','param',[1 2])
drawnow
-%%
+% Show 3D pattern
disp( 'calculating 3D far field pattern and dumping to vtk (use Paraview to visualize)...' );
thetaRange = (0:2:180);
phiRange = (0:2:360) - 180;
diff --git a/openEMS/matlab/Tutorials/StripLine2MSL.m b/openEMS/matlab/Tutorials/StripLine2MSL.m
new file mode 100644
index 0000000..64eeaa1
--- /dev/null
+++ b/openEMS/matlab/Tutorials/StripLine2MSL.m
@@ -0,0 +1,133 @@
+%
+% Stripline to Microstrip Line Transition
+%
+% Describtion at:
+% <http://openems.de/index.php/Tutorial:_Stripline_to_MSL_Transition>
+%
+% Tested with
+% - Octave 4.0
+% - openEMS v0.0.35
+%
+% (C) 2017 Thorsten Liebig <thorsten.liebig@gmx.de>
+
+close all
+clear
+clc
+
+%% Setup the Simulation
+physical_constants;
+unit = 1e-6; % specify everything in um
+
+line_length = 15000; % line length of strip line and microstrip line
+substrate_width = 6000;
+air_spacer = 4000; % air spacer above the substrate
+
+msl_width = 500;
+msl_substrate_thickness = 254;
+
+strip_width = 500;
+strip_substrate_thickness = 512;
+
+connect_via_rad = 500/2;
+connect_via_gap = 1250/2;
+
+substrate_epr = 3.66;
+substrate_kappa = 1e-3 * 2*pi*2.45e9 * EPS0*substrate_epr; % substrate losses
+
+f_max = 10e9;
+resolution = 250;
+edge_res = 25;
+feed_shift = 2500;
+meas_shift = 5000;
+
+%% Setup FDTD Parameters & Excitation Function
+FDTD = InitFDTD();
+FDTD = SetGaussExcite( FDTD, f_max/2, f_max/2);
+BC = {'PML_8' 'PML_8' 'PEC' 'PEC' 'PEC' 'MUR'};
+FDTD = SetBoundaryCond( FDTD, BC );
+
+%% Setup CSXCAD Geometry & Mesh
+CSX = InitCSX();
+edge_mesh = [-1/3 2/3]*edge_res; % 1/3 - 2/3 rule for 2D metal edges
+
+mesh.x = SmoothMeshLines( [-connect_via_gap 0 connect_via_gap], 2*edge_res, 1.5 );
+mesh.x = SmoothMeshLines( [-line_length mesh.x line_length], resolution, 1.5);
+mesh.y = SmoothMeshLines( [0 msl_width/2+edge_mesh substrate_width/2], resolution/4 , 1.5);
+mesh.y = sort(unique([-mesh.y mesh.y]));
+mesh.z = SmoothMeshLines( [linspace(-strip_substrate_thickness,0,5) linspace(0,strip_substrate_thickness,5) linspace(strip_substrate_thickness,msl_substrate_thickness+strip_substrate_thickness,5) 2*strip_substrate_thickness+air_spacer] , resolution );
+CSX = DefineRectGrid( CSX, unit, mesh );
+
+% Create Substrate
+CSX = AddMaterial( CSX, 'RO4350B' );
+CSX = SetMaterialProperty( CSX, 'RO4350B', 'Epsilon', substrate_epr, 'Kappa', substrate_kappa );
+start = [mesh.x(1), mesh.y(1), -strip_substrate_thickness];
+stop = [mesh.x(end), mesh.y(end), +strip_substrate_thickness+msl_substrate_thickness];
+CSX = AddBox( CSX, 'RO4350B', 0, start, stop );
+
+% Create a PEC called 'metal' and 'gnd'
+CSX = AddMetal( CSX, 'gnd' );
+CSX = AddMetal( CSX, 'metal' );
+
+% Create strip line port (incl. metal stip line)
+start = [-line_length -strip_width/2 0];
+stop = [0 +strip_width/2 0];
+[CSX,port{1}] = AddStripLinePort( CSX, 100, 1, 'metal', start, stop, strip_substrate_thickness, 'x', [0 0 -1], 'ExcitePort', true, 'FeedShift', feed_shift, 'MeasPlaneShift', meas_shift );
+
+% Create MSL port on top
+start = [line_length -strip_width/2 strip_substrate_thickness+msl_substrate_thickness];
+stop = [0 +strip_width/2 strip_substrate_thickness];
+[CSX,port{2}] = AddMSLPort( CSX, 100, 2, 'metal', start, stop, 'x', [0 0 -1], 'MeasPlaneShift', meas_shift );
+
+% transitional via
+start = [0, 0, 0];
+stop = [0, 0, strip_substrate_thickness+msl_substrate_thickness];
+CSX = AddCylinder(CSX, 'metal', 100, start, stop, connect_via_rad);
+
+% metal plane between strip line and MSL, including hole for transition
+p(1,1) = mesh.x(1);
+p(2,1) = mesh.y(1);
+p(1,2) = 0;
+p(2,2) = mesh.y(1);
+for a = linspace(-pi, pi, 21)
+ p(1,end+1) = connect_via_gap*sin(a);
+ p(2,end) = connect_via_gap*cos(a);
+endfor
+p(1,end+1) = 0;
+p(2,end ) = mesh.y(1);
+p(1,end+1) = mesh.x(end);
+p(2,end ) = mesh.y(1);
+p(1,end+1) = mesh.x(end);
+p(2,end ) = mesh.y(end);
+p(1,end+1) = mesh.x(1);
+p(2,end ) = mesh.y(end);
+CSX = AddPolygon( CSX, 'gnd', 1, 'z', strip_substrate_thickness, p);
+
+%% Write/Show/Run the openEMS compatible xml-file
+Sim_Path = 'tmp';
+Sim_CSX = 'strip2msl.xml';
+
+[status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory
+[status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder
+
+WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX );
+CSXGeomPlot( [Sim_Path '/' Sim_CSX] );
+RunOpenEMS( Sim_Path, Sim_CSX );
+
+%% Post-Processing
+close all
+f = linspace( 0, f_max, 1601 );
+port = calcPort( port, Sim_Path, f, 'RefImpedance', 50);
+
+s11 = port{1}.uf.ref./ port{1}.uf.inc;
+s21 = port{2}.uf.ref./ port{1}.uf.inc;
+
+plot(f/1e9,20*log10(abs(s11)),'k-','LineWidth',2);
+hold on;
+grid on;
+plot(f/1e9,20*log10(abs(s21)),'r--','LineWidth',2);
+legend('S_{11}','S_{21}');
+ylabel('S-Parameter (dB)','FontSize',12);
+xlabel('frequency (GHz) \rightarrow','FontSize',12);
+ylim([-40 2]);
+
+
diff --git a/openEMS/matlab/h5readatt_octave.cc b/openEMS/matlab/h5readatt_octave.cc
index 8bd58d0..13e1765 100755
--- a/openEMS/matlab/h5readatt_octave.cc
+++ b/openEMS/matlab/h5readatt_octave.cc
@@ -5,13 +5,7 @@
// this special treatment is necessary because Win32-Octave ships with a very old hdf5 version (1.6.10)
void CloseH5Object(hid_t obj)
{
-#if ((H5_VERS_MAJOR == 1) && (H5_VERS_MINOR == 6))
- // try group close, than Dataset close
- if (H5Gclose(obj)<0)
- H5Dclose(obj);
-#else
H5Oclose(obj);
-#endif
}
DEFUN_DLD (h5readatt_octave, args, nargout, "h5readatt_octave(<File_Name>,<DataSet_Name>,<Attribute_Name>)")
@@ -30,7 +24,7 @@ DEFUN_DLD (h5readatt_octave, args, nargout, "h5readatt_octave(<File_Name>,<DataS
}
//suppress hdf5 error output
- H5Eset_auto1(NULL, NULL);
+ H5Eset_auto2(H5E_DEFAULT, NULL, NULL);
hid_t file = H5Fopen( args(0).string_value().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT );
if (file==-1)
@@ -39,17 +33,7 @@ DEFUN_DLD (h5readatt_octave, args, nargout, "h5readatt_octave(<File_Name>,<DataS
return retval;
}
-#if ((H5_VERS_MAJOR == 1) && (H5_VERS_MINOR == 6))
- // this special treatment is necessary because Win32-Octave ships with a very old hdf5 version (1.6.10)
- hid_t obj = -1;
- //try opening the group
- obj = H5Gopen(file, args(1).string_value().c_str());
- //try opening the dataset if group failed
- if (obj==-1)
- obj = H5Dopen(file, args(1).string_value().c_str());
-#else
hid_t obj = H5Oopen(file, args(1).string_value().c_str(), H5P_DEFAULT);
-#endif
if (obj==-1)
{
@@ -59,7 +43,8 @@ DEFUN_DLD (h5readatt_octave, args, nargout, "h5readatt_octave(<File_Name>,<DataS
return retval;
}
- hid_t attr = H5Aopen_name(obj, args(2).string_value().c_str());
+ hid_t attr = H5Aopen_by_name(obj, ".", args(2).string_value().c_str(), H5P_DEFAULT, H5P_DEFAULT);
+
if (attr==-1)
{
CloseH5Object(obj);
diff --git a/openEMS/matlab/plotRefl.m b/openEMS/matlab/plotRefl.m
new file mode 100644
index 0000000..5f88631
--- /dev/null
+++ b/openEMS/matlab/plotRefl.m
@@ -0,0 +1,145 @@
+function h = plotRefl(port, varargin)
+
+% h = plotRefl(port,varargin)
+%
+% plot the reflection coefficient of a port into a Smith chart.
+% left and right facing triangles mark the lower and upper cutoff
+% frequency of the pass bands. An asterisk marks the frequnecy with
+% the lowest reflection.
+%
+% input:
+% port: port data structure. Call calcPort with an appropriate
+% frequency vector before calling this routine
+%
+% output: graphics handle for further modification of the plot.
+%
+% variable input:
+% 'precision': - number of decimal places (floating point precision)
+% for the frequency (always in MHz), default is 0
+% 'threshold': - Threshold value (in dB) for the upper and lower
+% cutoff frequency, default is -3
+% example:
+% myport = calcPort(myport, Sim_Path, linspace(f_0-f_c, f_0+f_c, 200));
+% plotRefl(myport);
+%
+% See also calcPort
+%
+% openEMS matlab interface
+% -----------------------
+% author: Georg Michel
+
+%defaults
+precision = 0;
+threshold = -3;
+
+
+for n=1:2:numel(varargin)
+ if (strcmp(varargin{n},'precision')==1);
+ precision = varargin{n+1};
+ elseif (strcmp(varargin{n},'threshold')==1);
+ threshold = varargin{n+1};
+ else
+ warning('openEMS:polarFF',['unknown argument key: ''' varargin{n} '''']);
+ end
+end
+
+
+if ~isfield(port, 'uf')
+ error('Cannot plot the reflection coefficient. Please call calcPort first.');
+end
+
+s11 = port.uf.ref ./ port.uf.inc;
+ffmt = ['%.', num2str(precision), 'f'];
+
+
+
+figure; %new figure
+
+plot([-1, 1], [0, 0], 'k');
+
+
+axis ([-1.15, 1.15, -1.15, 1.15], "square");
+axis off;
+hold on
+
+ReZ = [.2; .5; 1; 2];
+ImZ = 1i * [1 2 5 2];
+Z = bsxfun(@plus, ReZ, linspace(-ImZ, ImZ, 256));
+Gamma = (Z-1)./(Z+1);
+plot(Gamma.', 'k');
+
+ReZ = [.5 .5 1 1 2 2 5 5 10 10];
+ImZ = 1i * [-.2; .2; -.5; .5; -1; 1; -2; 2; -5; 5];
+Z = bsxfun(@plus, linspace(0, ReZ, 256), ImZ);
+Gamma = (Z-1)./(Z+1);
+plot(Gamma.', 'k');
+
+
+angle = linspace (0, 2 * pi, 256); ReZ = [0 5 10];
+center = ReZ ./ (ReZ + 1);
+radius = 1 ./ (ReZ + 1);
+plot(bsxfun(@plus, bsxfun(@times, radius, cos(angle.')), center), bsxfun(@times, radius, sin(angle.')), 'k');
+
+
+% resistance
+ReZ = [0.2 0.5 1 2 5 10]; ImZ = zeros (1, length (ReZ));
+rho = (ReZ.^2 + ImZ.^2 - 1 + 2i * ImZ) ./ ((ReZ + 1).^2 + ImZ.^2);
+
+xoffset = [0.1 0.1 0.05 0.05 0.05 0.075];
+yoffset = -0.03;
+
+for idx = 1:length (ReZ)
+ text (real (rho(idx)) - xoffset(idx), ...
+ imag (rho(idx)) - yoffset, num2str (ReZ(idx)));
+end
+
+% reactance
+ReZ = [-0.06 -0.06 -0.06 -0.12 -0.5];
+ImZ = [0.2 0.5 1 2 5];
+
+
+rho = (ReZ.^2 + ImZ.^2 - 1 + 2i * ImZ) ./ ((ReZ + 1).^2 + ImZ.^2);
+
+for idx = 1:length (ImZ)
+ text (real (rho(idx)), imag (rho(idx)), [num2str(ImZ(idx)), "j"]);
+ text (real (rho(idx)), -imag (rho(idx)), [num2str(-ImZ(idx)), "j"]); end
+
+% zero
+rho = (-0.05.^2 + 0.^2 - 1) ./ ((-0.05 + 1).^2 + 0.^2);
+
+text (real (rho), imag (rho), '0');
+
+s11dB = 20*log10(abs(s11));
+
+upperind = s11dB(1:end-1) < threshold & s11dB(2:end) > threshold;
+lowerind = s11dB(1:end-1) > threshold & s11dB(2:end) < threshold;
+minind = nthargout(2, @min, s11dB);
+handle1 = plot(s11(lowerind),['<','b']);
+handle2 = plot(s11(upperind),['>','b']);
+handle3 = plot(s11(minind),['*', 'b']);
+llegend = num2str(port.f(lowerind)(1)/1e6, ffmt);
+ulegend = num2str(port.f(upperind)(1)/1e6, ffmt);
+
+if nnz(lowerind) > 1
+ for i= 2:nnz(lowerind)
+ llegend = strjoin({llegend, num2str(port.f(lowerind)(i)/1e6, ffmt)}, ', ');
+ end
+end
+
+if nnz(upperind) > 1
+ for i= 2:nnz(upperind)
+ ulegend = strjoin({ulegend, num2str(port.f(upperind)(i)/1e6, ffmt)}, ', ');
+ end
+end
+
+legend([handle1, handle2, handle3], {[llegend, " MHz"], ...
+ [ulegend, " MHz"], ...
+ [num2str(20*log10(abs(s11(minind))), "%4.0f"), ...
+ "dB @ ", num2str(port.f(minind)/1e6, ffmt), " MHz"]});
+h = plot(s11);
+
+if (nargout == 0)
+ clear h;
+end
+
+end \ No newline at end of file
diff --git a/openEMS/matlab/setup.m b/openEMS/matlab/setup.m
index df56b25..e208ad4 100644
--- a/openEMS/matlab/setup.m
+++ b/openEMS/matlab/setup.m
@@ -5,7 +5,7 @@ function setup()
%
% openEMS matlab/octave interface
% -----------------------
-% author: Thorsten Liebig (2011)
+% author: Thorsten Liebig (2011-2017)
disp('setting up openEMS matlab/octave interface')
@@ -16,18 +16,22 @@ cd(dir);
if isOctave()
disp('compiling oct files')
- fflush(stdout)
+ fflush(stdout);
if isunix
- [res, fn] = unix('find /usr/lib -name libhdf5.so');
- if length(fn)>0
- [hdf5lib_dir, hdf5lib_fn] = fileparts(fn);
+ [res, fn_so] = unix('find /usr/lib -name libhdf5.so');
+ [res, fn_h] = unix('find /usr/include -name hdf5.h');
+ if length(fn_so)>0 && length(fn_h)>0
+ [hdf5lib_dir, hdf5lib_fn] = fileparts(fn_so);
disp(["HDF5 library path found at: " hdf5lib_dir])
- mkoctfile(["-L" hdf5lib_dir ],"-lhdf5 -DH5_USE_16_API", "h5readatt_octave.cc")
+
+ [hdf5inc_dir, hdf5inc_fn] = fileparts(fn_h);
+ disp(["HDF5 include path found at: " hdf5inc_dir])
+ mkoctfile(["-L" hdf5lib_dir " -I" hdf5inc_dir],"-lhdf5", "h5readatt_octave.cc")
else
- mkoctfile -lhdf5 -DH5_USE_16_API h5readatt_octave.cc
+ mkoctfile -lhdf5 h5readatt_octave.cc
end
else
- mkoctfile -lhdf5 -DH5_USE_16_API h5readatt_octave.cc
+ mkoctfile -lhdf5 h5readatt_octave.cc
end
else
disp('Matlab does not need this function. It is Octave only.')
diff --git a/openEMS/nf2ff/CMakeLists.txt b/openEMS/nf2ff/CMakeLists.txt
index d43d519..86811b4 100644
--- a/openEMS/nf2ff/CMakeLists.txt
+++ b/openEMS/nf2ff/CMakeLists.txt
@@ -38,6 +38,7 @@ TARGET_LINK_LIBRARIES( nf2ff
tinyxml
${HDF5_LIBRARIES}
${Boost_LIBRARIES}
+ ${MPI_LIBRARIES}
)
ADD_EXECUTABLE( nf2ff_bin main.cpp )
@@ -47,5 +48,7 @@ TARGET_LINK_LIBRARIES(nf2ff_bin nf2ff)
INSTALL(TARGETS nf2ff_bin DESTINATION bin)
INSTALL(TARGETS nf2ff DESTINATION lib${LIB_SUFFIX})
+INSTALL(FILES ${HEADERS} DESTINATION include/openEMS)
+
#TODO tarball, debug, release
diff --git a/openEMS/openEMS_MPI.pro b/openEMS/openEMS_MPI.pro
deleted file mode 100644
index d564c2b..0000000
--- a/openEMS/openEMS_MPI.pro
+++ /dev/null
@@ -1,5 +0,0 @@
-# enable MPI support
-!win32:CONFIG += MPI_SUPPORT
-
-include(openEMS.pro)
-
diff --git a/openEMS/openems.cpp b/openEMS/openems.cpp
index 72476d3..1d08ef7 100644
--- a/openEMS/openems.cpp
+++ b/openEMS/openems.cpp
@@ -418,11 +418,13 @@ bool openEMS::SetupProcessing()
}
if (CylinderCoords)
proc->SetMeshType(Processing::CYLINDRICAL_MESH);
- if ((pb->GetProbeType()==1) || (pb->GetProbeType()==3) || (pb->GetProbeType()==11))
+ if ((pb->GetProbeType()==1) || (pb->GetProbeType()==3))
{
proc->SetDualTime(true);
proc->SetDualMesh(true);
}
+ if (pb->GetProbeType()==11)
+ proc->SetDualTime(true);
proc->SetProcessInterval(Nyquist/m_OverSampling);
if (pb->GetStartTime()>0 || pb->GetStopTime()>0)
proc->SetProcessStartStopTime(pb->GetStartTime(), pb->GetStopTime());
@@ -715,9 +717,7 @@ bool openEMS::Parse_XML_FDTDSetup(TiXmlElement* FDTD_Opts)
ihelp = 0;
FDTD_Opts->QueryIntAttribute("OverSampling",&ihelp);
- if (ihelp<2)
- this->SetOverSampling(2);
- else
+ if (ihelp>1)
this->SetOverSampling(ihelp);
// check for cell constant material averaging
@@ -821,6 +821,7 @@ bool openEMS::Parse_XML_FDTDSetup(TiXmlElement* FDTD_Opts)
this->SetTimeStep(dhelp);
if (FDTD_Opts->QueryDoubleAttribute("TimeStepFactor",&dhelp)==TIXML_SUCCESS)
this->SetTimeStepFactor(dhelp);
+ return true;
}
void openEMS::SetGaussExcite(double f0, double fc)
diff --git a/openEMS/tools/hdf5_file_reader.cpp b/openEMS/tools/hdf5_file_reader.cpp
index 5e81614..a6e1359 100644
--- a/openEMS/tools/hdf5_file_reader.cpp
+++ b/openEMS/tools/hdf5_file_reader.cpp
@@ -68,7 +68,7 @@ bool HDF5_File_Reader::OpenGroup(hid_t &file, hid_t &group, string groupName)
return 0;
}
- group = H5Gopen(file, groupName.c_str() );
+ group = H5Gopen(file, groupName.c_str(), H5P_DEFAULT );
if (group<0)
{
cerr << "HDF5_File_Reader::OpenGroup: can't open group """ << groupName << """" << endl;
@@ -198,7 +198,7 @@ bool HDF5_File_Reader::ReadDataSet(string ds_name, hsize_t &nDim, hsize_t* &dims
return false;
}
- hid_t dataset = H5Dopen(hdf5_file, ds_name.c_str() );
+ hid_t dataset = H5Dopen(hdf5_file, ds_name.c_str(), H5P_DEFAULT );
if (dataset<0)
{
cerr << "HDF5_File_Reader::ReadDataSet: dataset not found" << endl;
diff --git a/openEMS/tools/hdf5_file_writer.cpp b/openEMS/tools/hdf5_file_writer.cpp
index 81c5692..e3a2113 100644
--- a/openEMS/tools/hdf5_file_writer.cpp
+++ b/openEMS/tools/hdf5_file_writer.cpp
@@ -52,7 +52,7 @@ hid_t HDF5_File_Writer::OpenGroup(hid_t hdf5_file, string group)
vector<string> results;
boost::split(results, group, boost::is_any_of("/"));
- hid_t grp=H5Gopen(hdf5_file,"/");
+ hid_t grp=H5Gopen(hdf5_file,"/", H5P_DEFAULT);
if (grp<0)
{
cerr << "HDF5_File_Writer::OpenGroup: Error, opening root group " << endl;
@@ -66,7 +66,7 @@ hid_t HDF5_File_Writer::OpenGroup(hid_t hdf5_file, string group)
{
if (H5Lexists(grp, results.at(n).c_str(), H5P_DEFAULT))
{
- grp = H5Gopen(grp, results.at(n).c_str());
+ grp = H5Gopen(grp, results.at(n).c_str(), H5P_DEFAULT);
H5Gclose(old_grp);
if (grp<0)
{
@@ -76,7 +76,7 @@ hid_t HDF5_File_Writer::OpenGroup(hid_t hdf5_file, string group)
}
else
{
- grp = H5Gcreate(grp,results.at(n).c_str(),0);
+ grp = H5Gcreate(grp,results.at(n).c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
H5Gclose(old_grp);
if (grp<0)
{
@@ -140,7 +140,7 @@ bool HDF5_File_Writer::WriteRectMesh(unsigned int const* numLines, float const*
return false;
}
- hid_t mesh_grp = H5Gcreate(hdf5_file,"/Mesh",0);
+ hid_t mesh_grp = H5Gcreate(hdf5_file,"/Mesh", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
if (mesh_grp<0)
{
cerr << "HDF5_File_Writer::WriteRectMesh: Error, creating group ""/Mesh"" failed" << endl;
@@ -165,7 +165,7 @@ bool HDF5_File_Writer::WriteRectMesh(unsigned int const* numLines, float const*
{
hsize_t dims[1]={numLines[n]};
hid_t space = H5Screate_simple(1, dims, NULL);
- hid_t dataset = H5Dcreate(mesh_grp, names[n].c_str(), H5T_NATIVE_FLOAT, space, H5P_DEFAULT);
+ hid_t dataset = H5Dcreate(mesh_grp, names[n].c_str(), H5T_NATIVE_FLOAT, space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
float* array = new float[numLines[n]];
for (unsigned int i=0; i<numLines[n]; ++i)
{
@@ -408,7 +408,7 @@ bool HDF5_File_Writer::WriteData(std::string dataSetName, hid_t mem_type, void
for (size_t n=0;n<dim;++n)
dims[n]=datasize[n];
hid_t space = H5Screate_simple(dim, dims, NULL);
- hid_t dataset = H5Dcreate(group, dataSetName.c_str(), mem_type, space, H5P_DEFAULT);
+ hid_t dataset = H5Dcreate(group, dataSetName.c_str(), mem_type, space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
if (H5Dwrite(dataset, mem_type, space, H5P_DEFAULT, H5P_DEFAULT, field_buf))
{
cerr << "HDF5_File_Writer::WriteData: Error, writing to dataset failed" << endl;
@@ -451,7 +451,7 @@ bool HDF5_File_Writer::WriteAtrribute(std::string locName, std::string attr_name
hid_t dataspace_id = H5Screate_simple(1, &size, NULL);
/* Create a dataset attribute. */
- hid_t attribute_id = H5Acreate(loc, attr_name.c_str(), mem_type, dataspace_id,H5P_DEFAULT);
+ hid_t attribute_id = H5Acreate(loc, attr_name.c_str(), mem_type, dataspace_id, H5P_DEFAULT, H5P_DEFAULT);
if (attribute_id<0)
{
cerr << "HDF5_File_Writer::WriteAtrribute: Error, failed to create the attrbute" << endl;