nf2ff: new support for cylindrical mesh & added nf2ff phase center definition

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

matlab/CalcNF2FF.m

@@ -4,6 +4,10 @@ function nf2ff = CalcNF2FF(nf2ff, Sim_Path, freq, theta, phi, varargin)
 % Calculate the near-field to far-field transformation created by
 % CreateNF2FFBox
 %
+% IMPORTANT:
+% Make sure to define the correct nf2ff phase center, aka. central antenna
+% position! See optional parameter below!! Default is [0 0 0]
+%
 % parameter:
 % nf2ff:    data structure created by CreateNF2FFBox
 % Sim_Path: path to simulation data
@@ -11,6 +15,10 @@ function nf2ff = CalcNF2FF(nf2ff, Sim_Path, freq, theta, phi, varargin)
 % theta,phi: spherical coordinates to evaluate the far-field on
 %
 % optional paramater:
+% 'Center': nf2ff phase center, default is [0 0 0]
+%           !! Make sure the center is never outside of your nf2ff box!!
+%           Definition is the correct coordinate system necessary
+%           --> either Cartesian or cylindrical coordinates
 % 'Mode':   'Mode', 0 -> read only, if data already exist (default)
 %           'Mode', 1 -> calculate anyway, overwrite existing
 %           'Mode', 2 -> read only, fail if not existing

matlab/CreateNF2FFBox.m

@@ -27,6 +27,7 @@ end
 directions = ones(6,1);
 add_args = {};
 dump_type = 0;
+dump_mode = 1;
 
 for n=1:numel(varargin)/2
     if strcmp(varargin{2*n-1},'Frequency')
@@ -43,14 +44,20 @@ nf2ff.filenames_E = {[name '_E_xn'],[name '_E_xp'],[name '_E_yn'],[name '_E_yp']
 nf2ff.filenames_H = {[name '_H_xn'],[name '_H_xp'],[name '_H_yn'],[name '_H_yp'],[name '_H_zn'],[name '_H_zp']};
 nf2ff.directions = directions;
 
+if (isfield(CSX,'ATTRIBUTE'))
+    if (isfield(CSX.ATTRIBUTE,'CoordSystem'))
+        nf2ff.CoordSystem = CSX.ATTRIBUTE.CoordSystem;
+    end
+end
+
 for nd = 1:3
     pos = 2*nd-1;
     if (directions(pos))
         l_start = start;
         l_stop = stop;
         l_stop(nd) = start(nd);
-        CSX = AddBox( AddDump(CSX,nf2ff.filenames_E{pos},'DumpType',dump_type,'DumpMode',2,'FileType',1,add_args{:}), nf2ff.filenames_E{pos}, 0, l_start, l_stop );
+        CSX = AddBox( AddDump(CSX,nf2ff.filenames_E{pos},'DumpType',dump_type,'DumpMode',dump_mode,'FileType',1,add_args{:}), nf2ff.filenames_E{pos}, 0, l_start, l_stop );
-        CSX = AddBox( AddDump(CSX,nf2ff.filenames_H{pos},'DumpType',dump_type+1,'DumpMode',2,'FileType',1,add_args{:}), nf2ff.filenames_H{pos}, 0, l_start, l_stop );
+        CSX = AddBox( AddDump(CSX,nf2ff.filenames_H{pos},'DumpType',dump_type+1,'DumpMode',dump_mode,'FileType',1,add_args{:}), nf2ff.filenames_H{pos}, 0, l_start, l_stop );
     else
         nf2ff.filenames_E{pos}='';
         nf2ff.filenames_H{pos}='';
@@ -60,8 +67,8 @@ for nd = 1:3
         l_start = start;
         l_stop = stop;
         l_start(nd) = stop(nd);
-        CSX = AddBox( AddDump(CSX,nf2ff.filenames_E{pos},'DumpType',dump_type,'DumpMode',2,'FileType',1,add_args{:}), nf2ff.filenames_E{pos}, 0, l_start, l_stop );
+        CSX = AddBox( AddDump(CSX,nf2ff.filenames_E{pos},'DumpType',dump_type,'DumpMode',dump_mode,'FileType',1,add_args{:}), nf2ff.filenames_E{pos}, 0, l_start, l_stop );
-        CSX = AddBox( AddDump(CSX,nf2ff.filenames_H{pos},'DumpType',dump_type+1,'DumpMode',2,'FileType',1,add_args{:}), nf2ff.filenames_H{pos}, 0, l_start, l_stop );
+        CSX = AddBox( AddDump(CSX,nf2ff.filenames_H{pos},'DumpType',dump_type+1,'DumpMode',dump_mode,'FileType',1,add_args{:}), nf2ff.filenames_H{pos}, 0, l_start, l_stop );
     else
         nf2ff.filenames_E{pos}='';
         nf2ff.filenames_H{pos}='';

nf2ff/nf2ff.cpp

@@ -34,7 +34,7 @@
 //external libs
 #include "tinyxml.h"
 
-nf2ff::nf2ff(vector<float> freq, vector<float> theta, vector<float> phi, unsigned int numThreads)
+nf2ff::nf2ff(vector<float> freq, vector<float> theta, vector<float> phi, vector<float> center, unsigned int numThreads)
 {
 	m_freq = freq;
 
@@ -51,11 +51,10 @@ nf2ff::nf2ff(vector<float> freq, vector<float> theta, vector<float> phi, unsigne
 	m_nf2ff.resize(freq.size(),NULL);
 	for (size_t fn=0;fn<freq.size();++fn)
 	{
-		m_nf2ff.at(fn) = new nf2ff_calc(freq.at(fn),theta, phi);
+		m_nf2ff.at(fn) = new nf2ff_calc(freq.at(fn),theta, phi, center);
 		if (numThreads)
 			m_nf2ff.at(fn)->SetNumThreads(numThreads);
 	}
-
 	m_radius = 1;
 	m_Verbose = 0;
 }
@@ -99,6 +98,12 @@ bool nf2ff::AnalyseXMLNode(TiXmlElement* ti_nf2ff)
 		return false;
 	}
 	vector<float> freq = SplitString2Float(attr);
+
+	vector<float> center;
+	attr = ti_nf2ff->Attribute("Center");
+	if (attr!=NULL)
+		center = SplitString2Float(attr);
+
 	attr = ti_nf2ff->Attribute("Outfile");
 	if (attr==NULL)
 	{
@@ -152,7 +157,7 @@ bool nf2ff::AnalyseXMLNode(TiXmlElement* ti_nf2ff)
 	}
 	vector<float> phi = SplitString2Float(ti_phi_text->Value());
 
-	nf2ff* l_nf2ff = new nf2ff(freq,theta,phi,numThreads);
+	nf2ff* l_nf2ff = new nf2ff(freq,theta,phi,center,numThreads);
 	l_nf2ff->SetVerboseLevel(Verbose);
 
 	TiXmlElement* ti_Planes = ti_nf2ff->FirstChildElement();
@@ -266,7 +271,7 @@ bool nf2ff::AnalyseFile(string E_Field_file, string H_Field_file)
 	{
 		if (m_Verbose>1)
 			cerr << "nf2ff: f = " << m_freq.at(fn) << "Hz (" << fn+1 << "/" << m_nf2ff.size() << ") ...";
-		m_nf2ff.at(fn)->AddPlane(E_lines, E_numLines, E_fd_data.at(fn), H_fd_data.at(fn));
+		m_nf2ff.at(fn)->AddPlane(E_lines, E_numLines, E_fd_data.at(fn), H_fd_data.at(fn),E_meshType);
 		if (m_Verbose>1)
 			cerr << " done." << endl;
 	}

nf2ff/nf2ff.h

@@ -32,7 +32,7 @@ class TiXmlElement;
 class nf2ff
 {
 public:
-	nf2ff(vector<float> freq, vector<float> theta, vector<float> phi, unsigned int numThreads=0);
+	nf2ff(vector<float> freq, vector<float> theta, vector<float> phi, vector<float> center, unsigned int numThreads=0);
 	~nf2ff();
 
 	bool AnalyseFile(string E_Field_file, string H_Field_file);

nf2ff/nf2ff_calc.cpp

@@ -60,6 +60,8 @@ void nf2ff_calc_thread::operator()()
 	complex<float>**** E_field=m_data.E_field;
 	complex<float>**** H_field=m_data.H_field;
 
+	int mesh_type = m_data.mesh_type;
+
 	// calc Js and Ms (eq. 8.15a/b)
 	pos[ny]=0;
 	for (pos_t=0; pos_t<num_t; ++pos_t)
@@ -76,6 +78,20 @@ void nf2ff_calc_thread::operator()()
 			Ms[0][pos[0]][pos[1]][pos[2]] = normDir[2]*E_field[1][pos[0]][pos[1]][pos[2]] - normDir[1]*E_field[2][pos[0]][pos[1]][pos[2]];
 			Ms[1][pos[0]][pos[1]][pos[2]] = normDir[0]*E_field[2][pos[0]][pos[1]][pos[2]] - normDir[2]*E_field[0][pos[0]][pos[1]][pos[2]];
 			Ms[2][pos[0]][pos[1]][pos[2]] = normDir[1]*E_field[0][pos[0]][pos[1]][pos[2]] - normDir[0]*E_field[1][pos[0]][pos[1]][pos[2]];
+
+			//transform to cartesian coordinates
+			if (mesh_type==1)
+			{
+				Js[0][pos[0]][pos[1]][pos[2]] = (normDir[1]*H_field[2][pos[0]][pos[1]][pos[2]] - normDir[2]*H_field[1][pos[0]][pos[1]][pos[2]])*cos(lines[1][pos[1]]) \
+						- (normDir[2]*H_field[0][pos[0]][pos[1]][pos[2]] - normDir[0]*H_field[2][pos[0]][pos[1]][pos[2]])*sin(lines[1][pos[1]]);
+				Js[1][pos[0]][pos[1]][pos[2]] = (normDir[1]*H_field[2][pos[0]][pos[1]][pos[2]] - normDir[2]*H_field[1][pos[0]][pos[1]][pos[2]])*sin(lines[1][pos[1]]) \
+						+ (normDir[2]*H_field[0][pos[0]][pos[1]][pos[2]] - normDir[0]*H_field[2][pos[0]][pos[1]][pos[2]])*cos(lines[1][pos[1]]);
+
+				Ms[0][pos[0]][pos[1]][pos[2]] = (normDir[2]*E_field[1][pos[0]][pos[1]][pos[2]] - normDir[1]*E_field[2][pos[0]][pos[1]][pos[2]])*cos(lines[1][pos[1]]) \
+						- (normDir[0]*E_field[2][pos[0]][pos[1]][pos[2]] - normDir[2]*E_field[0][pos[0]][pos[1]][pos[2]])*sin(lines[1][pos[1]]);
+				Ms[1][pos[0]][pos[1]][pos[2]] = (normDir[2]*E_field[1][pos[0]][pos[1]][pos[2]] - normDir[1]*E_field[2][pos[0]][pos[1]][pos[2]])*sin(lines[1][pos[1]]) \
+						+ (normDir[0]*E_field[2][pos[0]][pos[1]][pos[2]] - normDir[2]*E_field[0][pos[0]][pos[1]][pos[2]])*cos(lines[1][pos[1]]);
+			}
 		}
 	}
 
@@ -84,6 +100,12 @@ void nf2ff_calc_thread::operator()()
 	complex<float>** m_Lt=m_data.m_Lt;
 	complex<float>** m_Lp=m_data.m_Lp;
 
+	float center[3] = {m_nf_calc->m_centerCoord[0],m_nf_calc->m_centerCoord[1],m_nf_calc->m_centerCoord[2]};
+	if (mesh_type==1)
+	{
+		center[0] = m_nf_calc->m_centerCoord[0]*cos(m_nf_calc->m_centerCoord[1]);
+		center[1] = m_nf_calc->m_centerCoord[0]*sin(m_nf_calc->m_centerCoord[1]);
+	}
 	// calc local Nt,Np,Lt and Lp
 	float area;
 	float cosT_cosP,cosP_sinT;
@@ -111,7 +133,10 @@ void nf2ff_calc_thread::operator()()
 				pos[nP] = m_start+pos_t;
 				for (pos[nPP]=0; pos[nPP]<numLines[nPP]; ++pos[nPP])
 				{
-						r_cos_psi = lines[0][pos[0]]*cosP_sinT + lines[1][pos[1]]*sinT_sinP + lines[2][pos[2]]*cosT;
+					if (mesh_type==0)
+						r_cos_psi = (lines[0][pos[0]]-center[0])*cosP_sinT + (lines[1][pos[1]]-center[1])*sinT_sinP + (lines[2][pos[2]]-center[2])*cosT;
+					else
+						r_cos_psi = ((lines[0][pos[0]]*cos(lines[1][pos[1]]))-center[0])*cosP_sinT + ((lines[0][pos[0]]*sin(lines[1][pos[1]]))-center[1])*sinT_sinP + (lines[2][pos[2]]-center[2])*cosT;
 					exp_jkr = exp(_I_*k*r_cos_psi);
 					area = edge_length_P[pos[nP]]*edge_length_PP[pos[nPP]];
 					m_Nt[tn][pn] += area*exp_jkr*(Js[0][pos[0]][pos[1]][pos[2]]*cosT_cosP + Js[1][pos[0]][pos[1]][pos[2]]*cosT_sinP \
@@ -122,7 +147,6 @@ void nf2ff_calc_thread::operator()()
 												  - Ms[2][pos[0]][pos[1]][pos[2]]*sinT);
 					m_Lp[tn][pn] += area*exp_jkr*(Ms[1][pos[0]][pos[1]][pos[2]]*cosP - Ms[0][pos[0]][pos[1]][pos[2]]*sinP);
 				}
-
 			}
 		}
 
@@ -135,7 +159,7 @@ void nf2ff_calc_thread::operator()()
 /***********************************************************************/
 
 
-nf2ff_calc::nf2ff_calc(float freq, vector<float> theta, vector<float> phi)
+nf2ff_calc::nf2ff_calc(float freq, vector<float> theta, vector<float> phi, vector<float> center)
 {
 	m_freq = freq;
 
@@ -156,7 +180,19 @@ nf2ff_calc::nf2ff_calc(float freq, vector<float> theta, vector<float> phi)
 	m_H_phi = Create2DArray<std::complex<float> >(numLines);
 	m_P_rad = Create2DArray<float>(numLines);
 
-	m_centerCoord[0]=m_centerCoord[1]=m_centerCoord[2]=0;
+	if (center.size()==3)
+	{
+		m_centerCoord[0]=center.at(0);
+		m_centerCoord[1]=center.at(1);
+		m_centerCoord[2]=center.at(2);
+	}
+	else if (center.size()>0)
+	{
+		cerr << "nf2ff_calc::nf2ff_calc: Warning: Center coordinates error, ignoring!" << endl;
+		m_centerCoord[0]=m_centerCoord[1]=m_centerCoord[2]=0.0;
+	}
+	else
+		m_centerCoord[0]=m_centerCoord[1]=m_centerCoord[2]=0.0;
 
 	m_radPower = 0;
 	m_maxDir = 0;
@@ -189,7 +225,7 @@ nf2ff_calc::~nf2ff_calc()
 	m_Barrier = NULL;
 }
 
-bool nf2ff_calc::AddPlane(float **lines, unsigned int* numLines, complex<float>**** E_field, complex<float>**** H_field)
+bool nf2ff_calc::AddPlane(float **lines, unsigned int* numLines, complex<float>**** E_field, complex<float>**** H_field, int MeshType)
 {
 	//find normal direction
 	int ny = -1;
@@ -231,6 +267,24 @@ bool nf2ff_calc::AddPlane(float **lines, unsigned int* numLines, complex<float>*
 	edge_length_PP[0]=0.5*(lines[nPP][1]-lines[nPP][0]);
 	edge_length_PP[numLines[nPP]-1]=0.5*(lines[nPP][numLines[nPP]-1]-lines[nPP][numLines[nPP]-2]);
 
+	//check for cylindrical mesh
+	if (MeshType==1)
+	{
+		if (ny==0) //surface a-z
+		{
+			for (unsigned int n=0;n<numLines[nP];++n)
+				edge_length_P[n]*=lines[0][0]; //angle-width * radius
+		}
+		else if (ny==2) //surface r-a
+		{
+			//calculate: area = delta_angle * delta_radius * center_radius
+			for (unsigned int n=1;n<numLines[nP]-1;++n)
+				edge_length_P[n]*=lines[nP][n];  //radius-width * center-radius
+			edge_length_P[0]*=(lines[nP][0]+0.5*edge_length_P[0]);
+			edge_length_P[numLines[nP]-1]*=(lines[nP][numLines[nP]-1]-0.5*edge_length_P[numLines[nP]-1]);
+		}
+	}
+
 	complex<float> power = 0;
 	float area;
 	for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
@@ -242,7 +296,6 @@ bool nf2ff_calc::AddPlane(float **lines, unsigned int* numLines, complex<float>*
 						 - E_field[nPP][pos[0]][pos[1]][pos[2]]*conj(H_field[nP][pos[0]][pos[1]][pos[2]]));
 				m_radPower += 0.5*area*real(power)*normDir[ny];
 			}
-
 	unsigned int numAngles[2] = {m_numTheta, m_numPhi};
 
 	// setup multi-threading jobs
@@ -255,6 +308,7 @@ bool nf2ff_calc::AddPlane(float **lines, unsigned int* numLines, complex<float>*
 	for (unsigned int n=0; n<m_numThreads; n++)
 	{
 		thread_data[n].ny=ny;
+		thread_data[n].mesh_type = MeshType;
 		thread_data[n].normDir=normDir;
 		thread_data[n].numLines=numLines;
 		thread_data[n].lines=lines;

nf2ff/nf2ff_calc.h

@@ -34,6 +34,7 @@ typedef struct
 {
 	//local working data IN
 	int ny;
+	int mesh_type;
 	float* normDir;
 	unsigned int* numLines;
 	float **lines;
@@ -71,7 +72,7 @@ class nf2ff_calc
 	// allow full data access to nf2ff_calc_thread class
 	friend class nf2ff_calc_thread;
 public:
-	nf2ff_calc(float freq, vector<float> theta, vector<float> phi);
+	nf2ff_calc(float freq, vector<float> theta, vector<float> phi, vector<float> center);
 	~nf2ff_calc();
 
 	float GetRadPower() const {return m_radPower;}
@@ -83,7 +84,7 @@ public:
 	unsigned int GetNumThreads() const {return m_numThreads;}
 	void SetNumThreads(unsigned int n) {m_numThreads=n;}
 
-	bool AddPlane(float **lines, unsigned int* numLines, complex<float>**** E_field, complex<float>**** H_field);
+	bool AddPlane(float **lines, unsigned int* numLines, complex<float>**** E_field, complex<float>**** H_field, int MeshType=0);
 
 protected:
 	float m_freq;