new nf2ff application, replacing the matlab version
Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>
This commit is contained in:
Thorsten Liebig
parent
5d0f08aaec
commit
4c249bd4a3
44
nf2ff/main.cpp
Normal file
44
nf2ff/main.cpp
Normal file
@ -0,0 +1,44 @@
|
||||
/*
|
||||
* Copyright (C) 2012 Thorsten Liebig (Thorsten.Liebig@gmx.de)
|
||||
*
|
||||
* This program is free software: you can redistribute it and/or modify
|
||||
* it under the terms of the GNU General Public License as published by
|
||||
* the Free Software Foundation, either version 3 of the License, or
|
||||
* (at your option) any later version.
|
||||
*
|
||||
* This program is distributed in the hope that it will be useful,
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
* GNU General Public License for more details.
|
||||
*
|
||||
* You should have received a copy of the GNU General Public License
|
||||
* along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include <iostream>
|
||||
|
||||
#include "nf2ff.h"
|
||||
|
||||
using namespace std;
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
cout << " ---------------------------------------------------------------------- " << endl;
|
||||
cout << " | nf2ff, near-field to far-field transformation for openEMS " << endl;
|
||||
cout << " | (C) 2012 Thorsten Liebig <thorsten.liebig@gmx.de> GPL license" << endl;
|
||||
cout << " ---------------------------------------------------------------------- " << endl;
|
||||
|
||||
if (argc<=1)
|
||||
{
|
||||
cout << " Usage: nf2ff <nf2ff-xml-file>" << endl << endl;
|
||||
cout << endl;
|
||||
exit(-1);
|
||||
}
|
||||
|
||||
if (argc>=2)
|
||||
{
|
||||
return !nf2ff::AnalyseXMLFile(argv[argc-1]);
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
365
nf2ff/nf2ff.cpp
Normal file
365
nf2ff/nf2ff.cpp
Normal file
@ -0,0 +1,365 @@
|
||||
/*
|
||||
* Copyright (C) 2012 Thorsten Liebig (Thorsten.Liebig@gmx.de)
|
||||
*
|
||||
* This program is free software: you can redistribute it and/or modify
|
||||
* it under the terms of the GNU General Public License as published by
|
||||
* the Free Software Foundation, either version 3 of the License, or
|
||||
* (at your option) any later version.
|
||||
*
|
||||
* This program is distributed in the hope that it will be useful,
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
* GNU General Public License for more details.
|
||||
*
|
||||
* You should have received a copy of the GNU General Public License
|
||||
* along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "nf2ff.h"
|
||||
#include "nf2ff_calc.h"
|
||||
#include "../tools/array_ops.h"
|
||||
#include "../tools/useful.h"
|
||||
#include "../tools/hdf5_file_reader.h"
|
||||
#include "../tools/hdf5_file_writer.h"
|
||||
#include <hdf5.h>
|
||||
#include <boost/algorithm/string.hpp>
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <vector>
|
||||
#include <cmath>
|
||||
#include <complex>
|
||||
#include <iostream>
|
||||
#include <sstream>
|
||||
|
||||
//external libs
|
||||
#include "tinyxml.h"
|
||||
|
||||
nf2ff::nf2ff(vector<float> freq, vector<float> theta, vector<float> phi, unsigned int numThreads)
|
||||
{
|
||||
m_freq = freq;
|
||||
|
||||
m_numTheta = theta.size();
|
||||
m_theta = new float[m_numTheta];
|
||||
for (size_t n=0;n<m_numTheta;++n)
|
||||
m_theta[n]=theta.at(n);
|
||||
|
||||
m_numPhi = phi.size();
|
||||
m_phi = new float[m_numPhi];
|
||||
for (size_t n=0;n<m_numPhi;++n)
|
||||
m_phi[n]=phi.at(n);
|
||||
|
||||
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);
|
||||
if (numThreads)
|
||||
m_nf2ff.at(fn)->SetNumThreads(numThreads);
|
||||
}
|
||||
|
||||
m_radius = 1;
|
||||
m_Verbose = 0;
|
||||
}
|
||||
|
||||
nf2ff::~nf2ff()
|
||||
{
|
||||
m_freq.clear();
|
||||
for (size_t fn=0;fn<m_nf2ff.size();++fn)
|
||||
delete m_nf2ff.at(fn);
|
||||
m_nf2ff.clear();
|
||||
|
||||
delete[] m_phi;
|
||||
m_phi = NULL;
|
||||
delete[] m_theta;
|
||||
m_theta = NULL;
|
||||
}
|
||||
|
||||
bool nf2ff::AnalyseXMLNode(TiXmlElement* ti_nf2ff)
|
||||
{
|
||||
if (ti_nf2ff==NULL)
|
||||
return false;
|
||||
|
||||
unsigned int numThreads=0;
|
||||
int ihelp=0;
|
||||
if (ti_nf2ff->QueryIntAttribute("NumThreads",&ihelp) == TIXML_SUCCESS)
|
||||
{
|
||||
numThreads = ihelp;
|
||||
cerr << "nf2ff: Set number of threads to: " << numThreads << endl;
|
||||
}
|
||||
int Verbose=0;
|
||||
if (ti_nf2ff->QueryIntAttribute("Verbose",&Verbose) == TIXML_SUCCESS)
|
||||
cerr << "nf2ff: Set verbose level to " << Verbose << endl;
|
||||
else
|
||||
Verbose = 0;
|
||||
|
||||
const char* attr = NULL;
|
||||
attr = ti_nf2ff->Attribute("freq");
|
||||
if (attr==NULL)
|
||||
{
|
||||
cerr << "nf2ff::AnalyseXMLNode: Can't read frequency inforamtions ... " << endl;
|
||||
return false;
|
||||
}
|
||||
vector<float> freq = SplitString2Float(attr);
|
||||
attr = ti_nf2ff->Attribute("Outfile");
|
||||
if (attr==NULL)
|
||||
{
|
||||
cerr << "nf2ff::AnalyseXMLNode: Can't read frequency inforamtions ... " << endl;
|
||||
return false;
|
||||
}
|
||||
string outfile = string(attr);
|
||||
if (outfile.empty())
|
||||
{
|
||||
cerr << "nf2ff::AnalyseXMLNode: outfile is empty, skipping nf2ff... " << endl;
|
||||
return false;
|
||||
}
|
||||
|
||||
TiXmlElement* ti_theta = ti_nf2ff->FirstChildElement("theta");
|
||||
if (ti_theta==NULL)
|
||||
{
|
||||
cerr << "nf2ff::AnalyseXMLNode: Can't read theta values ... " << endl;
|
||||
return false;
|
||||
}
|
||||
TiXmlNode* ti_theta_node = ti_theta->FirstChild();
|
||||
if (ti_theta_node==NULL)
|
||||
{
|
||||
cerr << "nf2ff::AnalyseXMLNode: Can't read theta text child ... " << endl;
|
||||
return false;
|
||||
}
|
||||
TiXmlText* ti_theta_text = ti_theta_node->ToText();
|
||||
if (ti_theta_text==NULL)
|
||||
{
|
||||
cerr << "nf2ff::AnalyseXMLNode: Can't read theta text values ... " << endl;
|
||||
return false;
|
||||
}
|
||||
vector<float> theta = SplitString2Float(ti_theta_text->Value());
|
||||
|
||||
TiXmlElement* ti_phi = ti_nf2ff->FirstChildElement("phi");
|
||||
if (ti_phi==NULL)
|
||||
{
|
||||
cerr << "nf2ff::AnalyseXMLNode: Can't read phi values ... " << endl;
|
||||
return false;
|
||||
}
|
||||
TiXmlNode* ti_phi_node = ti_phi->FirstChild();
|
||||
if (ti_phi_node==NULL)
|
||||
{
|
||||
cerr << "nf2ff::AnalyseXMLNode: Can't read phi text child ... " << endl;
|
||||
return false;
|
||||
}
|
||||
TiXmlText* ti_phi_text = ti_phi_node->ToText();
|
||||
if (ti_phi_text==NULL)
|
||||
{
|
||||
cerr << "nf2ff::AnalyseXMLNode: Can't read phi text values ... " << endl;
|
||||
return false;
|
||||
}
|
||||
vector<float> phi = SplitString2Float(ti_phi_text->Value());
|
||||
|
||||
nf2ff* l_nf2ff = new nf2ff(freq,theta,phi,numThreads);
|
||||
l_nf2ff->SetVerboseLevel(Verbose);
|
||||
|
||||
TiXmlElement* ti_Planes = ti_nf2ff->FirstChildElement();
|
||||
string E_name;
|
||||
string H_name;
|
||||
while (ti_Planes!=NULL)
|
||||
{
|
||||
E_name = string(ti_Planes->Attribute("E_Field"));
|
||||
H_name = string(ti_Planes->Attribute("H_Field"));
|
||||
if ((!E_name.empty()) && (!H_name.empty()))
|
||||
{
|
||||
if (l_nf2ff->AnalyseFile(E_name,H_name)==false)
|
||||
{
|
||||
cerr << "nf2ff::AnalyseXMLNode: Error, analysing Plane ... " << endl;
|
||||
return false;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
cerr << "nf2ff::AnalyseXMLNode: Error, invalid plane entry ... " << endl;
|
||||
return false;
|
||||
}
|
||||
ti_Planes = ti_Planes->NextSiblingElement("Planes");
|
||||
}
|
||||
l_nf2ff->Write2HDF5(outfile);
|
||||
delete l_nf2ff;
|
||||
return true;
|
||||
}
|
||||
|
||||
bool nf2ff::AnalyseXMLFile(string filename)
|
||||
{
|
||||
TiXmlDocument doc(filename.c_str());
|
||||
if (!doc.LoadFile())
|
||||
{
|
||||
cerr << "nf2ff::AnalyseXMLFile: Error loading xml-file failed!!! File: " << filename << endl;
|
||||
return false;
|
||||
}
|
||||
TiXmlElement* ti_nf2ff = doc.FirstChildElement("nf2ff");
|
||||
if (ti_nf2ff==NULL)
|
||||
{
|
||||
cerr << "nf2ff::AnalyseXMLFile: Can't read nf2ff ... " << endl;
|
||||
return false;
|
||||
}
|
||||
|
||||
return AnalyseXMLNode(ti_nf2ff);
|
||||
}
|
||||
|
||||
bool nf2ff::AnalyseFile(string E_Field_file, string H_Field_file)
|
||||
{
|
||||
HDF5_File_Reader E_file(E_Field_file);
|
||||
HDF5_File_Reader H_file(H_Field_file);
|
||||
|
||||
if (m_Verbose>0)
|
||||
cerr << "nf2ff: Reading planes: " << E_Field_file << " & " << E_Field_file << endl;
|
||||
|
||||
// read E-mesh
|
||||
float* E_lines[3]={NULL,NULL,NULL};
|
||||
unsigned int E_numLines[3];
|
||||
int E_meshType;
|
||||
if (E_file.ReadMesh(E_lines, E_numLines, E_meshType) == false)
|
||||
{
|
||||
cerr << "nf2ff::AnalyseFile: Error reading E-field mesh..." << endl;
|
||||
return false;
|
||||
}
|
||||
|
||||
// read H-mesh
|
||||
float* H_lines[3]={NULL,NULL,NULL};
|
||||
unsigned int H_numLines[3];
|
||||
int H_meshType;
|
||||
if (H_file.ReadMesh(H_lines, H_numLines, H_meshType) == false)
|
||||
{
|
||||
cerr << "nf2ff::AnalyseFile: Error reading H-Field mesh..." << endl;
|
||||
return false;
|
||||
}
|
||||
|
||||
// compare E/H meshs
|
||||
if (E_meshType!=H_meshType)
|
||||
{
|
||||
cerr << "nf2ff::AnalyseFile: Error mesh types don't agree" << endl;
|
||||
return false;
|
||||
}
|
||||
if ((E_numLines[0]!=H_numLines[0]) || (E_numLines[1]!=H_numLines[1]) || (E_numLines[2]!=H_numLines[2]))
|
||||
{
|
||||
cerr << "nf2ff::AnalyseFile: Error mesh dimensions don't agree" << endl;
|
||||
return false;
|
||||
}
|
||||
for (int n=0;n<3;++n)
|
||||
for (unsigned int m=0;m<E_numLines[n];++m)
|
||||
if (E_lines[n][m]!=H_lines[n][m])
|
||||
{
|
||||
cerr << "nf2ff::AnalyseFile: Error mesh lines don't agree" << endl;
|
||||
return false;
|
||||
}
|
||||
|
||||
if (m_Verbose>1)
|
||||
cerr << "nf2ff: Data-Size: " << E_numLines[0] << "x" << E_numLines[1] << "x" << E_numLines[2] << endl;
|
||||
if (m_Verbose>1)
|
||||
cerr << "nf2ff: calculate dft..." << endl;
|
||||
|
||||
unsigned int data_size[4];
|
||||
vector<complex<float>****> E_fd_data;
|
||||
E_file.CalcFDVectorData(m_freq,E_fd_data,data_size);
|
||||
|
||||
vector<complex<float>****> H_fd_data;
|
||||
H_file.CalcFDVectorData(m_freq,H_fd_data,data_size);
|
||||
|
||||
if (m_Verbose>0)
|
||||
cerr << "nf2ff: Analysing far-field for " << m_nf2ff.size() << " frequencies. " << endl;
|
||||
|
||||
for (size_t fn=0;fn<m_nf2ff.size();++fn)
|
||||
{
|
||||
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));
|
||||
if (m_Verbose>1)
|
||||
cerr << " done." << endl;
|
||||
}
|
||||
for (int n=0;n<3;++n)
|
||||
{
|
||||
delete[] H_lines[n];
|
||||
delete[] E_lines[n];
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
bool nf2ff::Write2HDF5(string filename)
|
||||
{
|
||||
HDF5_File_Writer hdf_file(filename);
|
||||
|
||||
//write mesh information
|
||||
hdf_file.SetCurrentGroup("/Mesh");
|
||||
size_t meshsize[1]={m_numTheta};
|
||||
if (hdf_file.WriteData(string("theta"),m_theta,1,meshsize)==false)
|
||||
return false;
|
||||
meshsize[0]=m_numPhi;
|
||||
if (hdf_file.WriteData(string("phi"),m_phi,1,meshsize)==false)
|
||||
return false;
|
||||
meshsize[0]=1;
|
||||
float rad[1]={m_radius};
|
||||
if (hdf_file.WriteData(string("r"),rad,1,meshsize)==false)
|
||||
return false;
|
||||
|
||||
float attr_value = 2;
|
||||
hdf_file.WriteAtrribute("/Mesh", "MeshType", &attr_value, 1);
|
||||
|
||||
//write field data
|
||||
size_t dim = 2;
|
||||
size_t pos = 0;
|
||||
size_t datasize[2]={m_numPhi,m_numTheta};
|
||||
size_t size = datasize[0]*datasize[1];
|
||||
float* buffer = new float[size];
|
||||
complex<float>** field_data;
|
||||
string field_names[2]={"E_theta", "E_phi"};
|
||||
for (int n=0;n<2;++n)
|
||||
{
|
||||
hdf_file.SetCurrentGroup("/nf2ff/" + field_names[n] + "/FD");
|
||||
for (size_t fn=0;fn<m_freq.size();++fn)
|
||||
{
|
||||
stringstream ss;
|
||||
ss << "f" << fn;
|
||||
pos = 0;
|
||||
if (n==0)
|
||||
field_data = GetETheta(fn);
|
||||
else
|
||||
field_data = GetEPhi(fn);
|
||||
for (size_t j=0;j<m_numPhi;++j)
|
||||
for (size_t i=0;i<m_numTheta;++i)
|
||||
{
|
||||
buffer[pos++]=real(field_data[i][j]);
|
||||
}
|
||||
if (hdf_file.WriteData(ss.str() + "_real",buffer,dim,datasize)==false)
|
||||
{
|
||||
delete[] buffer;
|
||||
cerr << "nf2ff::Write2HDF5: Error writing field data" << endl;
|
||||
return false;
|
||||
}
|
||||
|
||||
pos = 0;
|
||||
for (size_t j=0;j<m_numPhi;++j)
|
||||
for (size_t i=0;i<m_numTheta;++i)
|
||||
{
|
||||
buffer[pos++]=imag(field_data[i][j]);
|
||||
}
|
||||
if (hdf_file.WriteData(ss.str() + "_imag",buffer,dim,datasize)==false)
|
||||
{
|
||||
delete[] buffer;
|
||||
cerr << "nf2ff::Write2HDF5: Error writing field data" << endl;
|
||||
return false;
|
||||
}
|
||||
}
|
||||
}
|
||||
delete[] buffer;
|
||||
|
||||
//write frequency attribute
|
||||
hdf_file.WriteAtrribute("/nf2ff", "Frequency",m_freq);
|
||||
|
||||
buffer = new float[m_freq.size()];
|
||||
//write radiated power attribute
|
||||
for (size_t fn=0;fn<m_freq.size();++fn)
|
||||
buffer[fn] = GetRadPower(fn);
|
||||
hdf_file.WriteAtrribute("/nf2ff", "Prad",buffer,m_freq.size());
|
||||
|
||||
//write max directivity attribute
|
||||
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;
|
||||
return true;
|
||||
}
|
||||
65
nf2ff/nf2ff.h
Normal file
65
nf2ff/nf2ff.h
Normal file
@ -0,0 +1,65 @@
|
||||
/*
|
||||
* Copyright (C) 2012 Thorsten Liebig (Thorsten.Liebig@gmx.de)
|
||||
*
|
||||
* This program is free software: you can redistribute it and/or modify
|
||||
* it under the terms of the GNU General Public License as published by
|
||||
* the Free Software Foundation, either version 3 of the License, or
|
||||
* (at your option) any later version.
|
||||
*
|
||||
* This program is distributed in the hope that it will be useful,
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
* GNU General Public License for more details.
|
||||
*
|
||||
* You should have received a copy of the GNU General Public License
|
||||
* along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef NF2FF_H
|
||||
#define NF2FF_H
|
||||
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <vector>
|
||||
#include <cmath>
|
||||
#include <complex>
|
||||
#include "nf2ff_calc.h"
|
||||
|
||||
using namespace std;
|
||||
|
||||
class TiXmlElement;
|
||||
|
||||
class nf2ff
|
||||
{
|
||||
public:
|
||||
nf2ff(vector<float> freq, vector<float> theta, vector<float> phi, unsigned int numThreads=0);
|
||||
~nf2ff();
|
||||
|
||||
bool AnalyseFile(string E_Field_file, string H_Field_file);
|
||||
|
||||
float GetRadPower(size_t f_idx) const {return m_nf2ff.at(f_idx)->GetRadPower();}
|
||||
float GetMaxDirectivity(size_t f_idx) const {return m_nf2ff.at(f_idx)->GetMaxDirectivity();}
|
||||
|
||||
complex<float>** GetETheta(size_t f_idx) const {return m_nf2ff.at(f_idx)->GetETheta();}
|
||||
complex<float>** GetEPhi(size_t f_idx) const {return m_nf2ff.at(f_idx)->GetEPhi();}
|
||||
|
||||
//! Write results to a hdf5 file
|
||||
bool Write2HDF5(string filename);
|
||||
|
||||
void SetVerboseLevel(int level) {m_Verbose=level;}
|
||||
|
||||
static bool AnalyseXMLNode(TiXmlElement* ti_nf2ff);
|
||||
static bool AnalyseXMLFile(string filename);
|
||||
|
||||
protected:
|
||||
vector<float> m_freq;
|
||||
unsigned int m_numTheta;
|
||||
unsigned int m_numPhi;
|
||||
float* m_theta;
|
||||
float* m_phi;
|
||||
float m_radius;
|
||||
int m_Verbose;
|
||||
vector<nf2ff_calc*> m_nf2ff;
|
||||
};
|
||||
|
||||
#endif // NF2FF_H
|
||||
60
nf2ff/nf2ff.pro
Normal file
60
nf2ff/nf2ff.pro
Normal file
@ -0,0 +1,60 @@
|
||||
TARGET = nf2ff
|
||||
CONFIG += console
|
||||
CONFIG -= app_bundle qt
|
||||
TEMPLATE = app
|
||||
OBJECTS_DIR = obj
|
||||
INCLUDEPATH += .
|
||||
INCLUDEPATH += ../../tinyxml
|
||||
CONFIG += debug_and_release
|
||||
|
||||
win32 {
|
||||
INCLUDEPATH += ../../hdf5/include ../../hdf5/include/cpp ../../boost/include/boost-1_42
|
||||
LIBS += ../../hdf5/lib/hdf5.lib
|
||||
LIBS += ../../boost/lib/libboost_thread-mgw44-mt.lib
|
||||
LIBS += ../../tinyxml/release/libtinyxml2.a
|
||||
}
|
||||
!win32 {
|
||||
LIBS += -lboost_thread
|
||||
LIBS += -lhdf5
|
||||
LIBS += ../../tinyxml/libtinyxml.so
|
||||
}
|
||||
QMAKE_LFLAGS += \'-Wl,-rpath,\$$ORIGIN/../../tinyxml\'
|
||||
|
||||
TOOLSPATH = ../tools
|
||||
|
||||
#### SOURCES ################################################################
|
||||
SOURCES += main.cpp \
|
||||
nf2ff.cpp \
|
||||
nf2ff_calc.cpp
|
||||
|
||||
# tools
|
||||
SOURCES += $$TOOLSPATH/global.cpp \
|
||||
$$TOOLSPATH/useful.cpp \
|
||||
$$TOOLSPATH/array_ops.cpp \
|
||||
$$TOOLSPATH/hdf5_file_reader.cpp \
|
||||
$$TOOLSPATH/hdf5_file_writer.cpp
|
||||
|
||||
#### HEADERS ################################################################
|
||||
HEADERS += nf2ff.h \
|
||||
nf2ff_calc.h
|
||||
|
||||
# tools
|
||||
HEADERS += $$TOOLSPATH/constants.h \
|
||||
$$TOOLSPATH/array_ops.h \
|
||||
$$TOOLSPATH/global.h \
|
||||
$$TOOLSPATH/useful.h \
|
||||
$$TOOLSPATH/aligned_allocator.h \
|
||||
$$TOOLSPATH/hdf5_file_reader.h \
|
||||
$$TOOLSPATH/hdf5_file_writer.h
|
||||
|
||||
QMAKE_CXXFLAGS_RELEASE = -O3 \
|
||||
-g \
|
||||
-march=native
|
||||
QMAKE_CXXFLAGS_DEBUG = -O0 \
|
||||
-g \
|
||||
-march=native
|
||||
|
||||
# add git revision
|
||||
# QMAKE_CXXFLAGS += -DGIT_VERSION=\\\"`git describe --tags`\\\"
|
||||
|
||||
|
||||
353
nf2ff/nf2ff_calc.cpp
Normal file
353
nf2ff/nf2ff_calc.cpp
Normal file
@ -0,0 +1,353 @@
|
||||
/*
|
||||
* Copyright (C) 2012 Thorsten Liebig (Thorsten.Liebig@gmx.de)
|
||||
*
|
||||
* This program is free software: you can redistribute it and/or modify
|
||||
* it under the terms of the GNU General Public License as published by
|
||||
* the Free Software Foundation, either version 3 of the License, or
|
||||
* (at your option) any later version.
|
||||
*
|
||||
* This program is distributed in the hope that it will be useful,
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
* GNU General Public License for more details.
|
||||
*
|
||||
* You should have received a copy of the GNU General Public License
|
||||
* along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "nf2ff_calc.h"
|
||||
#include "../tools/array_ops.h"
|
||||
#include "../tools/useful.h"
|
||||
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <vector>
|
||||
#include <cmath>
|
||||
#include <complex>
|
||||
#include <iostream>
|
||||
#include <sstream>
|
||||
|
||||
nf2ff_calc_thread::nf2ff_calc_thread(nf2ff_calc* nfc, unsigned int start, unsigned int stop, unsigned int threadID, nf2ff_data &data)
|
||||
{
|
||||
m_nf_calc = nfc;
|
||||
m_start = start;
|
||||
m_stop = stop;
|
||||
m_threadID = threadID;
|
||||
m_data = data;
|
||||
}
|
||||
|
||||
void nf2ff_calc_thread::operator()()
|
||||
{
|
||||
m_nf_calc->m_Barrier->wait(); // start
|
||||
|
||||
int ny = m_data.ny;
|
||||
int nP = (ny+1)%3;
|
||||
int nPP = (ny+2)%3;
|
||||
|
||||
unsigned int* numLines = m_data.numLines;
|
||||
float* normDir = m_data.normDir;
|
||||
float **lines = m_data.lines;
|
||||
float* edge_length_P = m_data.edge_length_P;
|
||||
float* edge_length_PP = m_data.edge_length_PP;
|
||||
|
||||
unsigned int pos[3];
|
||||
unsigned int pos_t=0;
|
||||
unsigned int num_t=m_stop-m_start+1;
|
||||
|
||||
|
||||
complex<float>**** Js=m_data.Js;
|
||||
complex<float>**** Ms=m_data.Ms;
|
||||
complex<float>**** E_field=m_data.E_field;
|
||||
complex<float>**** H_field=m_data.H_field;
|
||||
|
||||
// calc Js and Ms (eq. 8.15a/b)
|
||||
pos[ny]=0;
|
||||
for (pos_t=0; pos_t<num_t; ++pos_t)
|
||||
{
|
||||
pos[nP] = m_start+pos_t;
|
||||
for (pos[nPP]=0; pos[nPP]<numLines[nPP]; ++pos[nPP])
|
||||
{
|
||||
// Js = n x H
|
||||
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]];
|
||||
Js[1][pos[0]][pos[1]][pos[2]] = normDir[2]*H_field[0][pos[0]][pos[1]][pos[2]] - normDir[0]*H_field[2][pos[0]][pos[1]][pos[2]];
|
||||
Js[2][pos[0]][pos[1]][pos[2]] = normDir[0]*H_field[1][pos[0]][pos[1]][pos[2]] - normDir[1]*H_field[0][pos[0]][pos[1]][pos[2]];
|
||||
|
||||
// Ms = -n x E
|
||||
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]];
|
||||
}
|
||||
}
|
||||
|
||||
complex<float>** m_Nt=m_data.m_Nt;
|
||||
complex<float>** m_Np=m_data.m_Np;
|
||||
complex<float>** m_Lt=m_data.m_Lt;
|
||||
complex<float>** m_Lp=m_data.m_Lp;
|
||||
|
||||
// calc local Nt,Np,Lt and Lp
|
||||
float area;
|
||||
float cosT_cosP,cosP_sinT;
|
||||
float cosT_sinP,sinT_sinP;
|
||||
float sinT,sinP;
|
||||
float cosP,cosT;
|
||||
float r_cos_psi;
|
||||
float k = 2*M_PI*m_nf_calc->m_freq/__C0__;
|
||||
complex<float> exp_jkr;
|
||||
complex<float> _I_(0,1);
|
||||
for (unsigned int tn=0;tn<m_nf_calc->m_numTheta;++tn)
|
||||
for (unsigned int pn=0;pn<m_nf_calc->m_numPhi;++pn)
|
||||
{
|
||||
sinT = sin(m_nf_calc->m_theta[tn]);
|
||||
sinP = sin(m_nf_calc->m_phi[pn]);
|
||||
cosT = cos(m_nf_calc->m_theta[tn]);
|
||||
cosP = cos(m_nf_calc->m_phi[pn]);
|
||||
cosT_cosP = cosT*cosP;
|
||||
cosT_sinP = cosT*sinP;
|
||||
cosP_sinT = cosP*sinT;
|
||||
sinT_sinP = sinP*sinT;
|
||||
|
||||
for (pos_t=0; pos_t<num_t; ++pos_t)
|
||||
{
|
||||
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;
|
||||
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 \
|
||||
- Js[2][pos[0]][pos[1]][pos[2]]*sinT);
|
||||
m_Np[tn][pn] += area*exp_jkr*(Js[1][pos[0]][pos[1]][pos[2]]*cosP - Js[0][pos[0]][pos[1]][pos[2]]*sinP);
|
||||
|
||||
m_Lt[tn][pn] += area*exp_jkr*(Ms[0][pos[0]][pos[1]][pos[2]]*cosT_cosP + Ms[1][pos[0]][pos[1]][pos[2]]*cosT_sinP \
|
||||
- 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);
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
m_nf_calc->m_Barrier->wait(); //combine all thread local Nt,Np,Lt and Lp
|
||||
|
||||
m_nf_calc->m_Barrier->wait(); //wait for termination
|
||||
}
|
||||
|
||||
|
||||
/***********************************************************************/
|
||||
|
||||
|
||||
nf2ff_calc::nf2ff_calc(float freq, vector<float> theta, vector<float> phi)
|
||||
{
|
||||
m_freq = freq;
|
||||
|
||||
m_numTheta = theta.size();
|
||||
m_theta = new float[m_numTheta];
|
||||
for (size_t n=0;n<m_numTheta;++n)
|
||||
m_theta[n]=theta.at(n);
|
||||
|
||||
m_numPhi = phi.size();
|
||||
m_phi = new float[m_numPhi];
|
||||
for (size_t n=0;n<m_numPhi;++n)
|
||||
m_phi[n]=phi.at(n);
|
||||
|
||||
unsigned int numLines[2] = {m_numTheta, m_numPhi};
|
||||
m_E_theta = Create2DArray<std::complex<float> >(numLines);
|
||||
m_E_phi = Create2DArray<std::complex<float> >(numLines);
|
||||
m_H_theta = Create2DArray<std::complex<float> >(numLines);
|
||||
m_H_phi = Create2DArray<std::complex<float> >(numLines);
|
||||
m_P_rad = Create2DArray<float>(numLines);
|
||||
|
||||
m_centerCoord[0]=m_centerCoord[1]=m_centerCoord[2]=0;
|
||||
|
||||
m_radPower = 0;
|
||||
m_maxDir = 0;
|
||||
m_radius = 1;
|
||||
|
||||
m_Barrier = NULL;
|
||||
m_numThreads = boost::thread::hardware_concurrency();
|
||||
}
|
||||
|
||||
nf2ff_calc::~nf2ff_calc()
|
||||
{
|
||||
delete[] m_phi;
|
||||
m_phi = NULL;
|
||||
delete[] m_theta;
|
||||
m_theta = NULL;
|
||||
|
||||
unsigned int numLines[2] = {m_numTheta, m_numPhi};
|
||||
Delete2DArray(m_E_theta,numLines);
|
||||
m_E_theta = NULL;
|
||||
Delete2DArray(m_E_phi,numLines);
|
||||
m_E_phi = NULL;
|
||||
Delete2DArray(m_H_theta,numLines);
|
||||
m_H_theta = NULL;
|
||||
Delete2DArray(m_H_phi,numLines);
|
||||
m_H_phi = NULL;
|
||||
Delete2DArray(m_P_rad,numLines);
|
||||
m_P_rad = NULL;
|
||||
|
||||
delete m_Barrier;
|
||||
m_Barrier = NULL;
|
||||
}
|
||||
|
||||
bool nf2ff_calc::AddPlane(float **lines, unsigned int* numLines, complex<float>**** E_field, complex<float>**** H_field)
|
||||
{
|
||||
//find normal direction
|
||||
int ny = -1;
|
||||
int nP,nPP;
|
||||
for (int n=0;n<3;++n)
|
||||
{
|
||||
nP = (n+1)%3;
|
||||
nPP = (n+2)%3;
|
||||
if ((numLines[n]==1) && (numLines[nP]>2) && (numLines[nPP]>2))
|
||||
ny=n;
|
||||
}
|
||||
nP = (ny+1)%3;
|
||||
nPP = (ny+2)%3;
|
||||
if (ny<0)
|
||||
{
|
||||
cerr << "nf2ff_calc::AddPlane: Error can't determine normal direction..." << endl;
|
||||
return false;
|
||||
}
|
||||
|
||||
complex<float>**** Js = Create_N_3DArray<complex<float> >(numLines);
|
||||
complex<float>**** Ms = Create_N_3DArray<complex<float> >(numLines);
|
||||
|
||||
float normDir[3]= {0,0,0};
|
||||
if (lines[ny][0]>=m_centerCoord[ny])
|
||||
normDir[ny]=1;
|
||||
else
|
||||
normDir[ny]=-1;
|
||||
unsigned int pos[3];
|
||||
|
||||
float edge_length_P[numLines[nP]];
|
||||
for (unsigned int n=1;n<numLines[nP]-1;++n)
|
||||
edge_length_P[n]=0.5*(lines[nP][n+1]-lines[nP][n-1]);
|
||||
edge_length_P[0]=0.5*(lines[nP][1]-lines[nP][0]);
|
||||
edge_length_P[numLines[nP]-1]=0.5*(lines[nP][numLines[nP]-1]-lines[nP][numLines[nP]-2]);
|
||||
|
||||
float edge_length_PP[numLines[nPP]];
|
||||
for (unsigned int n=1;n<numLines[nPP]-1;++n)
|
||||
edge_length_PP[n]=0.5*(lines[nPP][n+1]-lines[nPP][n-1]);
|
||||
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]);
|
||||
|
||||
complex<float> power = 0;
|
||||
float area;
|
||||
for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
|
||||
for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
|
||||
for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
|
||||
{
|
||||
area = edge_length_P[pos[nP]]*edge_length_PP[pos[nPP]];
|
||||
power = (E_field[nP][pos[0]][pos[1]][pos[2]]*conj(H_field[nPP][pos[0]][pos[1]][pos[2]]) \
|
||||
- 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
|
||||
vector<unsigned int> jpt = AssignJobs2Threads(numLines[nP], m_numThreads, true);
|
||||
m_numThreads = jpt.size();
|
||||
nf2ff_data thread_data[m_numThreads];
|
||||
m_Barrier = new boost::barrier(m_numThreads+1); // numThread workers + 1 controller
|
||||
unsigned int start=0;
|
||||
unsigned int stop=jpt.at(0)-1;
|
||||
for (unsigned int n=0; n<m_numThreads; n++)
|
||||
{
|
||||
thread_data[n].ny=ny;
|
||||
thread_data[n].normDir=normDir;
|
||||
thread_data[n].numLines=numLines;
|
||||
thread_data[n].lines=lines;
|
||||
thread_data[n].edge_length_P=edge_length_P;
|
||||
thread_data[n].edge_length_PP=edge_length_PP;
|
||||
thread_data[n].E_field=E_field;
|
||||
thread_data[n].H_field=H_field;
|
||||
thread_data[n].Js=Js;
|
||||
thread_data[n].Ms=Ms;
|
||||
thread_data[n].m_Nt=Create2DArray<complex<float> >(numAngles);
|
||||
thread_data[n].m_Np=Create2DArray<complex<float> >(numAngles);
|
||||
thread_data[n].m_Lt=Create2DArray<complex<float> >(numAngles);
|
||||
thread_data[n].m_Lp=Create2DArray<complex<float> >(numAngles);
|
||||
|
||||
boost::thread *t = new boost::thread( nf2ff_calc_thread(this,start,stop,n,thread_data[n]) );
|
||||
|
||||
m_thread_group.add_thread( t );
|
||||
|
||||
start = stop+1;
|
||||
if (n<m_numThreads-1)
|
||||
stop = start + jpt.at(n+1)-1;
|
||||
}
|
||||
//all threads a running and waiting for the barrier
|
||||
|
||||
m_Barrier->wait(); //start
|
||||
|
||||
// threads: calc Js and Ms (eq. 8.15a/b)
|
||||
// threads calc their local Nt,Np,Lt and Lp
|
||||
|
||||
m_Barrier->wait(); //combine all thread local Nt,Np,Lt and Lp
|
||||
|
||||
//cleanup E- & H-Fields
|
||||
Delete_N_3DArray(E_field,numLines);
|
||||
Delete_N_3DArray(H_field,numLines);
|
||||
|
||||
complex<float>** Nt = Create2DArray<complex<float> >(numAngles);
|
||||
complex<float>** Np = Create2DArray<complex<float> >(numAngles);
|
||||
complex<float>** Lt = Create2DArray<complex<float> >(numAngles);
|
||||
complex<float>** Lp = Create2DArray<complex<float> >(numAngles);
|
||||
|
||||
for (unsigned int n=0; n<m_numThreads; n++)
|
||||
{
|
||||
for (unsigned int tn=0;tn<m_numTheta;++tn)
|
||||
for (unsigned int pn=0;pn<m_numPhi;++pn)
|
||||
{
|
||||
Nt[tn][pn] += thread_data[n].m_Nt[tn][pn];
|
||||
Np[tn][pn] += thread_data[n].m_Np[tn][pn];
|
||||
Lt[tn][pn] += thread_data[n].m_Lt[tn][pn];
|
||||
Lp[tn][pn] += thread_data[n].m_Lp[tn][pn];
|
||||
}
|
||||
Delete2DArray(thread_data[n].m_Nt,numAngles);
|
||||
Delete2DArray(thread_data[n].m_Np,numAngles);
|
||||
Delete2DArray(thread_data[n].m_Lt,numAngles);
|
||||
Delete2DArray(thread_data[n].m_Lp,numAngles);
|
||||
}
|
||||
|
||||
m_Barrier->wait(); //wait for termination
|
||||
m_thread_group.join_all(); // wait for termination
|
||||
delete m_Barrier;
|
||||
m_Barrier = NULL;
|
||||
|
||||
//cleanup Js & Ms
|
||||
Delete_N_3DArray(Js,numLines);
|
||||
Delete_N_3DArray(Ms,numLines);
|
||||
|
||||
// calc equations 8.23a/b and 8.24a/b
|
||||
float k = 2*M_PI*m_freq/__C0__;
|
||||
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 P_max = 0;
|
||||
for (unsigned int tn=0;tn<m_numTheta;++tn)
|
||||
for (unsigned int pn=0;pn<m_numPhi;++pn)
|
||||
{
|
||||
m_E_theta[tn][pn] -= factor*(Lp[tn][pn] + Z0*Nt[tn][pn]);
|
||||
m_E_phi[tn][pn] += factor*(Lt[tn][pn] - Z0*Np[tn][pn]);
|
||||
|
||||
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])));
|
||||
if (m_P_rad[tn][pn]>P_max)
|
||||
P_max = m_P_rad[tn][pn];
|
||||
}
|
||||
|
||||
//cleanup Nx and Lx
|
||||
Delete2DArray(Nt,numAngles);
|
||||
Delete2DArray(Np,numAngles);
|
||||
Delete2DArray(Lt,numAngles);
|
||||
Delete2DArray(Lp,numAngles);
|
||||
|
||||
m_maxDir = 4*M_PI*P_max / m_radPower;
|
||||
|
||||
return true;
|
||||
}
|
||||
114
nf2ff/nf2ff_calc.h
Normal file
114
nf2ff/nf2ff_calc.h
Normal file
@ -0,0 +1,114 @@
|
||||
/*
|
||||
* Copyright (C) 2012 Thorsten Liebig (Thorsten.Liebig@gmx.de)
|
||||
*
|
||||
* This program is free software: you can redistribute it and/or modify
|
||||
* it under the terms of the GNU General Public License as published by
|
||||
* the Free Software Foundation, either version 3 of the License, or
|
||||
* (at your option) any later version.
|
||||
*
|
||||
* This program is distributed in the hope that it will be useful,
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
* GNU General Public License for more details.
|
||||
*
|
||||
* You should have received a copy of the GNU General Public License
|
||||
* along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef NF2FF_CALC_H
|
||||
#define NF2FF_CALC_H
|
||||
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <vector>
|
||||
#include <cmath>
|
||||
#include <complex>
|
||||
#include <boost/thread.hpp>
|
||||
|
||||
using namespace std;
|
||||
|
||||
class nf2ff_calc;
|
||||
|
||||
// data structure to exchange data between thread-controller and worker-threads
|
||||
typedef struct
|
||||
{
|
||||
//local working data IN
|
||||
int ny;
|
||||
float* normDir;
|
||||
unsigned int* numLines;
|
||||
float **lines;
|
||||
float* edge_length_P;
|
||||
float* edge_length_PP;
|
||||
|
||||
complex<float>**** E_field;
|
||||
complex<float>**** H_field;
|
||||
complex<float>**** Js;
|
||||
complex<float>**** Ms;
|
||||
|
||||
//local working data OUT
|
||||
complex<float>** m_Nt;
|
||||
complex<float>** m_Np;
|
||||
complex<float>** m_Lt;
|
||||
complex<float>** m_Lp;
|
||||
|
||||
} nf2ff_data;
|
||||
|
||||
class nf2ff_calc_thread
|
||||
{
|
||||
public:
|
||||
nf2ff_calc_thread(nf2ff_calc* nfc, unsigned int start, unsigned int stop, unsigned int threadID, nf2ff_data &data);
|
||||
void operator()();
|
||||
|
||||
protected:
|
||||
unsigned int m_start, m_stop, m_threadID;
|
||||
nf2ff_calc *m_nf_calc;
|
||||
|
||||
nf2ff_data m_data;
|
||||
};
|
||||
|
||||
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 GetRadPower() const {return m_radPower;}
|
||||
float GetMaxDirectivity() const {return m_maxDir;}
|
||||
|
||||
complex<float>** GetETheta() const {return m_E_theta;}
|
||||
complex<float>** GetEPhi() const {return m_E_phi;}
|
||||
|
||||
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);
|
||||
|
||||
protected:
|
||||
float m_freq;
|
||||
float m_radius;
|
||||
|
||||
float m_radPower;
|
||||
float m_maxDir;
|
||||
|
||||
complex<float>** m_E_theta;
|
||||
complex<float>** m_E_phi;
|
||||
complex<float>** m_H_theta;
|
||||
complex<float>** m_H_phi;
|
||||
float** m_P_rad;
|
||||
|
||||
float m_centerCoord[3];
|
||||
unsigned int m_numTheta;
|
||||
unsigned int m_numPhi;
|
||||
float* m_theta;
|
||||
float* m_phi;
|
||||
|
||||
//boost multi-threading
|
||||
unsigned int m_numThreads;
|
||||
boost::thread_group m_thread_group;
|
||||
boost::barrier *m_Barrier;
|
||||
};
|
||||
|
||||
|
||||
#endif // NF2FF_CALC_H
|
||||
Loading…
Reference in New Issue
Block a user