nf2ff: allow Radius, Eps_r and Mue_r definition

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

matlab/CalcNF2FF.m

@@ -25,6 +25,9 @@ function nf2ff = CalcNF2FF(nf2ff, Sim_Path, freq, theta, phi, varargin)
 % 'Outfile': alternative nf2ff result hdf5 file name
 %            default is: <nf2ff.name>.h5
 % 'Verbose': set verbose level for the nf2ff calculation 0-2 supported
+% 'Radius':  specify the radius for the nf2ff
+% 'Eps_r':   specify the relative electric permittivity for the nf2ff
+% 'Mue_r':   specify the relative magnetic permeability for the nf2ff
 %
 % See also: CreateNF2FFBox, ReadNF2FF
 %

matlab/private/ReadNF2FF.m

@@ -23,6 +23,17 @@ nf2ff.freq = ReadHDF5Attribute(file,'/nf2ff','Frequency');
 nf2ff.Prad = ReadHDF5Attribute(file,'/nf2ff','Prad');
 nf2ff.Dmax = ReadHDF5Attribute(file,'/nf2ff','Dmax');
 
+try
+    nf2ff.Eps_r = ReadHDF5Attribute(file,'/nf2ff','Eps_r');
+catch
+    nf2ff.Eps_r = ones(size(nf2ff.freq));
+end
+try
+    nf2ff.Mue_r = ReadHDF5Attribute(file,'/nf2ff','Mue_r');
+catch
+    nf2ff.Mue_r = ones(size(nf2ff.freq));
+end
+
 if isOctave
     hdf = load( '-hdf5', file );
     for n=1:numel(nf2ff.freq)

nf2ff/nf2ff.cpp

@@ -72,6 +72,44 @@ nf2ff::~nf2ff()
 	m_theta = NULL;
 }
 
+void nf2ff::SetRadius(float radius)
+{
+	m_radius = radius;
+	for (size_t fn=0;fn<m_nf2ff.size();++fn)
+		m_nf2ff.at(fn)->SetRadius(radius);
+}
+
+void nf2ff::SetPermittivity(vector<float> permittivity)
+{
+	if (permittivity.size()==0)
+		return;
+
+	if (permittivity.size()!=m_freq.size())
+	{
+		cerr << __func__ << ": Error, permittivity vector size must match number of set frequencies! skipping!" << endl;
+		return;
+	}
+	m_permittivity = permittivity;
+	for (size_t fn=0;fn<m_nf2ff.size();++fn)
+		m_nf2ff.at(fn)->SetPermittivity(permittivity.at(fn));
+}
+
+void nf2ff::SetPermeability(vector<float> permeability)
+{
+	if (permeability.size()==0)
+		return;
+
+	if (permeability.size()!=m_freq.size())
+	{
+		cerr << __func__ << ": Error, permeability vector size must match number of set frequencies! skipping!" << endl;
+		return;
+	}
+	m_permeability = permeability;
+	for (size_t fn=0;fn<m_nf2ff.size();++fn)
+		m_nf2ff.at(fn)->SetPermeability(permeability.at(fn));
+
+}
+
 bool nf2ff::AnalyseXMLNode(TiXmlElement* ti_nf2ff)
 {
 	if (ti_nf2ff==NULL)
@@ -160,6 +198,18 @@ bool nf2ff::AnalyseXMLNode(TiXmlElement* ti_nf2ff)
 	nf2ff* l_nf2ff = new nf2ff(freq,theta,phi,center,numThreads);
 	l_nf2ff->SetVerboseLevel(Verbose);
 
+	attr = ti_nf2ff->Attribute("Eps_r");
+	if (attr!=NULL)
+		l_nf2ff->SetPermittivity(SplitString2Float(attr));
+
+	attr = ti_nf2ff->Attribute("Mue_r");
+	if (attr!=NULL)
+		l_nf2ff->SetPermeability(SplitString2Float(attr));
+
+	float radius = 1;
+	if (ti_nf2ff->QueryFloatAttribute("Radius",&radius) ==  TIXML_SUCCESS)
+		l_nf2ff->SetRadius(radius);
+
 	TiXmlElement* ti_Planes = ti_nf2ff->FirstChildElement();
 	string E_name;
 	string H_name;
@@ -538,12 +588,32 @@ bool nf2ff::Write2HDF5(string filename)
 	for (size_t fn=0;fn<m_freq.size();++fn)
 		buffer[fn] = GetTotalRadPower(fn);
 	hdf_file.WriteAtrribute("/nf2ff", "Prad",buffer,m_freq.size());
+	delete[] buffer;
 
 	//write max directivity attribute
+	buffer = new double[m_freq.size()];
 	for (size_t fn=0;fn<m_freq.size();++fn)
 		buffer[fn] = GetMaxDirectivity(fn);
 	hdf_file.WriteAtrribute("/nf2ff", "Dmax",buffer,m_freq.size());
-
 	delete[] buffer;
+
+	if (m_permittivity.size()>0)
+	{
+		buffer = new double[m_permittivity.size()];
+		for (size_t n=0;n<m_permittivity.size();++n)
+			buffer[n] = m_permittivity.at(n);
+		hdf_file.WriteAtrribute("/nf2ff", "Eps_r",buffer,m_permittivity.size());
+		delete[] buffer;
+	}
+
+	if (m_permeability.size()>0)
+	{
+		buffer = new double[m_permeability.size()];
+		for (size_t n=0;n<m_permeability.size();++n)
+			buffer[n] = m_permeability.at(n);
+		hdf_file.WriteAtrribute("/nf2ff", "Mue_r",buffer,m_permeability.size());
+		delete[] buffer;
+	}
+
 	return true;
 }

nf2ff/nf2ff.h

@@ -37,6 +37,10 @@ public:
 
 	bool AnalyseFile(string E_Field_file, string H_Field_file);
 
+	void SetRadius(float radius);
+	void SetPermittivity(vector<float> permittivity);
+	void SetPermeability(vector<float> permeability);
+
 	double GetTotalRadPower(size_t f_idx) const {return m_nf2ff.at(f_idx)->GetTotalRadPower();}
 	double GetMaxDirectivity(size_t f_idx) const {return m_nf2ff.at(f_idx)->GetMaxDirectivity();}
 
@@ -54,6 +58,8 @@ public:
 
 protected:
 	vector<float> m_freq;
+	vector<float> m_permittivity;
+	vector<float> m_permeability;
 	unsigned int m_numTheta;
 	unsigned int m_numPhi;
 	float* m_theta;

nf2ff/nf2ff_calc.cpp

@@ -113,7 +113,7 @@ void nf2ff_calc_thread::operator()()
 	float sinT,sinP;
 	float cosP,cosT;
 	float r_cos_psi;
-	float k = 2*M_PI*m_nf_calc->m_freq/__C0__;
+	float k = 2*M_PI*m_nf_calc->m_freq/__C0__*sqrt(m_nf_calc->m_permittivity*m_nf_calc->m_permeability);
 	complex<float> exp_jkr;
 	complex<float> _I_(0,1);
 	for (unsigned int tn=0;tn<m_nf_calc->m_numTheta;++tn)
@@ -162,6 +162,8 @@ void nf2ff_calc_thread::operator()()
 nf2ff_calc::nf2ff_calc(float freq, vector<float> theta, vector<float> phi, vector<float> center)
 {
 	m_freq = freq;
+	m_permittivity = 1;
+	m_permeability = 1;
 
 	m_numTheta = theta.size();
 	m_theta = new float[m_numTheta];
@@ -375,11 +377,12 @@ bool nf2ff_calc::AddPlane(float **lines, unsigned int* numLines, complex<float>*
 	Delete_N_3DArray(Ms,numLines);
 
 	// calc equations 8.23a/b and 8.24a/b
-	float k = 2*M_PI*m_freq/__C0__;
+	float k = 2*M_PI*m_freq/__C0__*sqrt(m_permittivity*m_permeability);
 	complex<float> factor(0,k/4.0/M_PI/m_radius);
 	complex<float> f_exp(0,-1*k*m_radius);
 	factor *= exp(f_exp);
-	complex<float> Z0 = __Z0__;
+	float fZ0 = __Z0__ * sqrt(m_permeability/m_permittivity);
+	complex<float> Z0 = fZ0;
 	float P_max = 0;
 	for (unsigned int tn=0;tn<m_numTheta;++tn)
 		for (unsigned int pn=0;pn<m_numPhi;++pn)
@@ -390,7 +393,7 @@ bool nf2ff_calc::AddPlane(float **lines, unsigned int* numLines, complex<float>*
 			m_H_theta[tn][pn] += factor*(Np[tn][pn] - Lt[tn][pn]/Z0);
 			m_H_phi[tn][pn] -= factor*(Nt[tn][pn] + Lp[tn][pn]/Z0);
 
-			m_P_rad[tn][pn] = m_radius*m_radius/(2*__Z0__) * abs((m_E_theta[tn][pn]*conj(m_E_theta[tn][pn])+m_E_phi[tn][pn]*conj(m_E_phi[tn][pn])));
+			m_P_rad[tn][pn] = m_radius*m_radius/(2*fZ0) * abs((m_E_theta[tn][pn]*conj(m_E_theta[tn][pn])+m_E_phi[tn][pn]*conj(m_E_phi[tn][pn])));
 			if (m_P_rad[tn][pn]>P_max)
 				P_max = m_P_rad[tn][pn];
 		}

nf2ff/nf2ff_calc.h

@@ -75,6 +75,10 @@ public:
 	nf2ff_calc(float freq, vector<float> theta, vector<float> phi, vector<float> center);
 	~nf2ff_calc();
 
+	void SetRadius(float radius) {m_radius=radius;}
+	void SetPermittivity(float permittivity) {m_permittivity=permittivity;}
+	void SetPermeability(float permeability) {m_permeability=permeability;}
+
 	double GetTotalRadPower() const {return m_radPower;}
 	double GetMaxDirectivity() const {return m_maxDir;}
 
@@ -91,6 +95,9 @@ protected:
 	float m_freq;
 	float m_radius;
 
+	float m_permittivity; //relative electric permittivity
+	float m_permeability; //relative magnetic permeability
+
 	double m_radPower;
 	double m_maxDir;