|
libRWG and libTRDT are C++ libraries for
use in developing codes for boundary-element analysis of
3D (libRWG) and 2D (libTDRT) electromagnetic
scattering problems. These libraries implement algorithms for
(1) reading in geometry files describing mesh discretizations
of 3D and 2D objects (scatterers), (2) numerically evaluating
the integrals needed to compute the elements of the matrices and
vectors involved in standard BEM formulations of EM scattering
problems (the EFIE and PMCHW formulations), and (3) performing
a variety of postprocessing and visualization techniques, including
computing the fields scattered from objects illuminated by given
sources.
libRWG and libTDRT include support
for both real frequencies (oscillatory fields, i.e. radar problems,
antenna problems, and other ordinary electromagnetic scattering
problems) and imaginary frequencies (decaying fields, as used
e.g. in numerical computations of Casimir interactions).
libRWG and libTDRT are named after
the basis functions they support. libRWG works with
the "RWG" basis functions introduced by Rao,
Wilton, and Glisson
(IEEE Transactions on Antennas and Propagation AP-30 409 (1982)). libTDRT works with Two
Dimensional Roof-Top functions. For a
detailed description of these basis functions and their use
in solving electromagnetic scattering problems, see the
overview section below.
Although libRWG and libTDRT use
separate code bases, their structures and APIs are sufficiently
similar that we can describe them both together in this document. Most
users will only need to use one or the other depending on whether
the scattering problem in question is of fundamentally 3D or 2D
character.
libRWG and libTDRT work with
objects and geometries. An object is a discretized
representation of a single, isolated, compact scatterer.
(libRWG and libTDRT don't do the
discretization for you -- you do that on your own, using external
meshing software, and supply the resulting meshfiles to
libRWG/libTDRT). A geometry is simply a
collection of objects, which may be displaced or rotated relative
to one another using methods supplied by libRWG/libTDRT.
The algorithms implemented by libRWG/libTDRT
can then be divided into two granularity levels, described
by two portions of the API:
-
a high-level or geometry-level API,
containing algorithms that operate on a geometry
as a whole, and
-
a low-level API, containing algorithms that
operate on individual objects in the geometry and the
subelements of which those objects are comprised.
Roughly speaking, if you just want to use libRWG/libTDRT
to solve scattering problems and compute the scattered fields, you
will only need to use the high-level API. If you want to use these
libraries to develop your own new methods for solving scattering
problems, you will want to use the low-level API.
Table Of Contents
1. An Overview of RWG and TDRT Basis Functions
and Their Use in Solving Electromagnetic
Scattering Problems
This discussion is sufficiently involved that I think I will put it
into a PDF document instead of trying to do in on a webpage:
2. High-Level API
C++ code files that access libRWG and/or
libTDRT must #include the appropriate
header files.
#include <libRWG.h>
#include <libTDRT.h>
2.1 Creating a Geometry
The primary object exported by libRWG is a C++
class called RWGGeometry (for libTDRT,
it's TDRTGeometry). Before doing anything else,
you should create an instance of one of these classes. The
class constructor takes a single string argument, the name
of a geometry file describing the geometry.
RWGGeometry *R;
TDRTGeometry *T;
R=new RWGGeometry("Spheres.rwggeo");
T=new TDRTGeometry("Cylinder.tdgeo");
For details on the format of the .rwggeo and
.tdgeo files used to specify 3D and 2D geometries,
see
GeometryFiles.
2.2 Creating the BEM matrix for scattering problems
The first step in solving a scattering problem is to form the BEM matrix
at a given frequency (in the 3D case) or at a given frequency and
z-wavevector (in the 2D case). (For brevity in what follows we
will abbreviate this full phrase to "at a given frequency.")
The RWGGeometry and TDRTGeometry classes
offer class methods for assembling the BEM matrices, at real or
imaginary frequencies, in either the EFIE or PMCHW formulations.
The first step is to call AllocateBEMMatrix() to
allocate a matrix of the size appropriate for the BEM matrix. (The
return value of this routine is a pointer to a type called
HMatrix, which is a very simple little wrapper class
I put together to abstract away some of the details of how matrices
are stored and manipulated.) After you have allocated a matrix, you
call AssembleBEMMatrix() to fill it in.
Class Methods for Allocating and
Assembling BEM Matrices
| Dimension | Class Methods |
| 3D |
HMatrix *RWGGeometry::AllocateBEMMatrix(int RealFreq);
void RWGGeometry::AssembleBEMMatrix(double Frequency, int RealFreq,
int nThread, HMatrix *M);
|
| 2D |
HMatrix *TDRTGeometry::AllocateBEMMatrix(int RealFreq);
void TDRTGeometry::AssembleBEMMatrix(double Frequency, int RealFreq,
double q, int nThread, HMatrix *M);
|
The various parameters in these routines are summarized in the
following table.
Parameters in BEM Allocation / Assembly Routines
| Parameter | Type | Description |
Frequency |
double |
Real or imaginary frequency (Omega or Beta) depending on
the value of the RealFreq flag.
|
RealFreq |
int |
Set RealFreq=1 for real-frequency problems,
RealFreq=0 for imaginary-frequency problems.
(libRWG.h and libTDRT.h define
constants REAL_FREQ and IMAG_FREQ,
with values of 1 and 0 respectively, for use as values for
the RealFreq parameter to library routines.)
|
q |
double |
z wavevector for 2D routines. |
nThread |
int |
The number of threads spawned to assemble the matrix.
Should be set less than or equal to the number of
available cores on your machine.
|
M |
HMatrix * |
Pointer to an instance of the HMatrix class.
This is storage for the BEM matrix, with details
about the size of the matrix, and whether its entries
are real or complex, abstracted away. The matrix must be
created by a call to AllocateBEMMatrix
before it can be filled in by AssembleBEMMatrix.
|
2.3 Computing the RHS vector for scattering problems
The next step in solving a scattering problem is to construct
the RHS vector (the vector we called B in the
memo). This requires the user to supply a routine for computing
the incident field at an arbitrary point. This routine must have
the following prototype:
void EHFunc(double *R, void *UserData, cdouble *EH);
where the parameters are described in the following table:
Parameters to EHFunc Routine
| Parameter |
Input/Output |
Description |
|---|
R |
Input |
Cartesian coordinates of point at which
incident fields are to be evaluated.
R[0..2]=x,y,z.
|
UserData |
Input |
May be used to pass any additional
parameters your incident-field routine
may need. |
EH |
Output |
On return, EH[0..2] must
contain the 3 Cartesian components of
the incident electric field at R.
For PMCHW routines, EH[3..5] must
also contain the 3 Cartesian components of the
incident magnetic field at R.
|
Two particularly common types of incident field are (a) plane
waves and (b) the field of a point source (electric or
magnetic dipole). My little auxiliary library
libIncField
contains prefabricated EHFunc functions for each
of these situations, so most users will never have to write their
own EFFunc routine.
(For those who do, we pause to note the following with regard to the
user-supplied EHFunc routine:
- There is no explicit frequency parameter to
EHFunc.
If your incident field depends on frequency, you will need to
pass that information using the UserData mechanism.
(And nobody will check to make sure that the frequency you use
when assembling the RHS is the same as the frequency you used
to assemble the BEM matrix, so be careful.)
- If all the objects in your geometry are perfect metals, the
BEM formulation of the scattering problem makes no reference
to the incident magnetic field, and hence your incident field
routine may leave
EH[3..5] uninitialized. (Otherwise,
your incident field routine must fill in all 6 slots in the
EH vector.)
- For 3D imaginary-frequency problems, the incident fields are
real-valued. In this case the imaginary components of
EH are ignored.)
To simplify function prototypes, libRWG.h
and libTDRT.h define a type for pointers to
EHFunc functions:
typedef void (*EHFuncType)(double *R, void *UserData, cdouble *EH);
Once you have created your EHFunc routine, you can
allocate space for, and initialize, an RHS vector using routines
similar to those you used to allocate and assemble the BEM
matrix. RHS vectors are stored using a data type named
HVector, designed for interoperation with
HMatrix.
Class Methods for Allocating and Assembling RHS Vectors
| Dimension | Class Method |
| 3D |
HVector *RWGGeometry::AllocateRHSVector(int RealFreq);
void RWGGeometry::AssembleRHSVector( EHFuncType EHFunc, void *UserData, int RealFreq,
int nThread, HVector *B);
|
| 3D |
HVector *TDRTGeometry::AllocateRHSVector(int RealFreq);
void TDRTGeometry::AssembleRHSVector( EHFuncType EHFunc, void *UserData, int RealFreq,
int nThread, HVector *B);
|
2.4 Solving scattering problems
Once you have assembled the BEM matrix and the RHS vector for a given
incident field, you must solve the linear system to compute the induced
current distribution on the object surfaces. libRWG and
libTDRT don't care how you do this, but the most direct way
is to use the following class methods provided by the
HMatrix class:
void HMatrix::LUFactorize();
int HMatrix::LUSolve(HVector *X);
The first of these replaces an HMatrix with its
LU-factorization as computed by LAPACK. The second uses the
LU-factorization to solve a linear system with RHS vector
X. (The solution vector is overwritten on the
input vector, so that when the routine returns X
is the solution to the system).
It is important to note that you only need to form and factorize
the BEM matrix once at a given frequency. Once you have done this,
you can solve the scattering problem for any number of incident fields
without recomputing the matrix (indeed, this fact is one of the key
virtues of the BEM in the first place). You will only need to recompute and
refactorize the BEM matrix if you change the frequency or apply geometrical
transformations to the geometry (see below).
Here is a code snippet illustrating how you might solve the
scattering problem for several different incident fields. This
is for a 3D imaginary-frequency problem.
RWGGeometry *G;
HMatrix *M;
HVector *X;
int nThread;
double Beta;
Beta=0.1;
nThread=GetNumProcs();
/* create the geometry */
G=new RWGGeometry("Spheres.rwggeo");
/* allocate, assemble, factorize BEM matrix */
M=G->AllocateBEMMatrix(IMAG_FREQ);
G->AssembleBEMMatrix(Beta, IMAG_FREQ, nThread, M);
M->LUFactorize();
/* form RHS for first incident field and solve scattering problem */
X=G->AllocateRHSVector(IMAG_FREQ);
G->AssembleRHSVector(MyEHFunc1, MyUserData1, IMAG_FREQ, nThread, X);
M->LUSolve(X);
/* do various things with the solution to that problem ... */
/* form RHS for second incident field and solve scattering problem */
G->AssembleRHSVector(MyEHFunc2, MyUserData2, IMAG_FREQ, nThread, X);
M->LUSolve(X);
/* now do various things with the solution to _that_ problem ... */
2.5 Computing scattered fields
After you have solved a scattering problem to obtain the solution vector
K, you can evaluate the scattered fields anywhere in space
using libRWG/libTDRT class methods.
Class Methods for Computing Scattered Fields
| Dimension | Class Method |
| 3D |
void RWGGeometry::GetFields( double *X, int WhichObject, double Frequency, int RealFreq,
HVector *X, int nThread, cdouble *EH);
void RWGGeometry::GetFields( double *X, double Frequency, int RealFreq,
HVector *X, int nThread, cdouble *EH);
|
| 2D |
void TDRTGeometry::GetFields( double *X, int WhichObject, double Frequency, int RealFreq, double q,
HVector *X, int nThread, cdouble *EH);
TDRTGeometry::GetFields( double *X, double Frequency, int RealFreq, double q,
HVector *X, int nThread, cdouble *EH);
|
Parameters in Scattered Field Routines
| Parameter | Type | Description |
X |
double[3] |
Cartesian coordinates of evaluation point. X[0..2]=x,y,z. |
WhichObject |
int |
If the evaluation point X lies in the
interior of one of the scattering objects in the
geometry, you must set WhichObject equal to
the index of that object (0 if the object is
the first object listed in the .rwggeo file,
1 if the second object, etc.). If the evaluation
point does not lie in the interior of a scattering
object, set WhichObject=-1, or use the
second of the two function prototypes above (in which the
WhichObject parameter does not appear).
|
Frequency |
double |
Real or imaginary frequency (Omega or Beta) depending on
the value of the RealFreq flag.
|
RealFreq |
int |
Set to 1 or 0 to select real or imaginary frequency. |
q |
double |
z wavevector for 2D routines. |
X |
double *, cdouble * |
Linear BEM system solution vector.
|
nThread |
int |
The number of threads spawned to compute the fields.
Should be set less than or equal to the number of
available cores on your machine.
|
EH |
double[6], cdouble[6] |
On entry, EH must point to a buffer with space
for 6 double (3D imaginary frequency) or
cdouble (all others) values.
On return, EH[0..2] are the three Cartesian
components of the scattered electric field at the evaluation
point X, and EH[3..5] are the
components of the scattered magnetic field.
|
2.6 Computing Dyadic Green's Functions and Stress Tensor VEVs
The electric and magnetic dyadic Green's functions are defined to be
the electric and magnetic fields arising from unit-strength electric
and magnetic sources. It is easy enough to compute these by doing six
separate libRWG scattering calculations. For the special
case in which the source and destination points in the Green's function
coincide, the six-scattering-calculation procedure is encoded in a
libRWG routine. (Needless to say, we are talking here about
the dyadic Green's functions with the vacuum component subtracted off,
as is natural for boundary-element analysis, so that G(x,xp) is
finite as x → xp.)
The prototype is
void RWGGeometry::GetGij(double *R, HMatrix *M, HVector *KN,
double Frequency, int RealFreq,
int nThread,
cdouble GE[3][3], cdouble GM[3][3]);
where
-
R[0..2] are the cartesian components of the
evaluation point
-
M is the LU-factorized BEM matrix
-
KN is a scratch vector: you have to
allocate KN via a call to AllocateRHSVector(),
but you don't have to do anything else to it
-
Frequency, RealFreq, nThread
are as in other routines described above
-
GE and GM must be user-allocated
3x3 arrays of cdoubles. On return, GE[i][j]
and GM[i][j] are the (i,j) components
of the electric and magnetic dyadic Green's functions at spatial
point R (i.e. both spatial arguments of the Green's
function are set equal to R.) Note that GE[i][j]
is the i component of the scattered electric field at
R due to a point electric current source at R
pointing in the j direction.
The next obvious requirement for Casimir applications is a routine
that evaluates the (vacuum expectation value of the) Maxwell (van der Waals)
stress tensor at a point in space. (This is done by evaluating the
electric and magnetic dyadic Green's functions and combining them
appropriately.)
The libRWG routine for doing this is
void GetTijNj(double R[3], double nHat[3],
HMatrix *M, HVector *KN,
double Frequency, int RealFreq,
int nThread,
double TE[3], double TM[3]);
where
-
R[0..2] are the cartesian components of the
evaluation point
-
nHat[0..2] are the cartesian components of the
vector normal to the bounding surface over which we are integrating
at R.
-
M is the LU-factorized BEM matrix
-
KN is a scratch vector allocated
by a call to AllocateRHSVector()
-
Frequency, RealFreq, nThread
are as in other routines described above
-
TE and TM are the cartesian components
of the force-per-unit-area on an infinitesimal patch of surface
at R with surface normal nHat.
2.7 Transforming a Geometry
2.8 Exporting to MATLAB
For some purposes it is useful to be able to play with the BEM system
matrices and vectors in MATLAB. For this purpose you can use the
ToMatlab() class method provided by both HMatrix
and HVector. This class method takes a string argument
and produces a .m file that you can execute at the
MATLAB command line to import the matrix or vector in
question:
HMatrix *MyMatrix;
HVector *MyVector;
...
MyMatrix->ToMatlab("M");
MyVector->ToMatlab("V");
This will produces files named Load_M.hd5, Load_V.hd5
and Load_M.m, Load_V.m in your current
working directory. The former are binary files containing the
actual data, while the latter are M-files you can execute
at the MATLAB prompt:
>> Load_M;
>> Load_V;
Now your MATLAB session contains a matrix named M
and a vector named V.
Note that the ToMatlab class method supports
printf-style variable argument lists:
int N=10;
int Omega=2.0;
MyMatrix->ToMatlab("M_%i_%g",N,Omega);
This will generate files M_10_2.0.hd5 and
M_10_2.0.m.
2.9 Visualization
2.9.1 Visualizing Geometries
2.9.2 Visualizing Induced Current Distributions
RWGGeometry offers a class method named
PlotSurfaceCurrents that will produce GMSH
files suitable for visualizing the surface currents induced by
the incident fields on the scatterers in your geometry.
Class Methods for Visualization
| Dimension | Class Method |
| 3D |
void *RWGGeometry::PlotSurfaceCurrents(HVector *KN, double Frequency, int RealFreq,
const char *format, ...);
|
Parameters in Visualization Routines
| Parameter | Type | Description |
KN |
HVector * |
The vector of expansion coefficients of the surface
current distribution to plot.
|
Frequency |
double |
Real or imaginary frequency depending on the
RealFreq parameter.
|
RealFreq |
int |
Set RealFreq=1 for real-frequency problems,
RealFreq=0 for imaginary-frequency problems.
|
format |
char * |
A printf-style format string specifying the
filename to which the visualization data will be written.
(The file may be subsequently opened in GMSH
for viewing). Any additional arguments needed by the
format string are passed immediately
after format.
|
2.10 Computing Induced Dipole Moments
RWGGeometry offers a class method named
GetDipoleMoments for computing the dipole
moments (electric and magnetic) induced by the incident field
on the objects in your scattering geometry.
void RWGGeometry::GetDipoleMoments(double Frequency, int RealFreq, HVector *KN, int nThread, cdouble (*PM)[6]);
The first four parameters are inputs. For the
PM parameter, you pass in a
two-dimensional array of dimension
(number of objects in geometry) x (6).
In other words, declare this array like this:
cdouble PM[G->NumObjects][6];
On return, PM is filled in as follows:
-
PM[0][0..2] are the x, y, z
components of the electric dipole moment induced
on the first scattering object in your geometry
(the first OBJECT in the .rwggeo
file.)
-
PM[0][3..5] are the x, y, z
components of the magnetic dipole moment induced on
the first scattering object in your geometry.
-
PM[1][0..2] are the x, y, z
components of the electric dipole moment induced
on the second scattering object in your geometry ,
- etc.
2.11 A complete real-frequency example: Radar cross-sections of
arbitrary objects
This program takes two command-line arguments: an
.rwggeo file and a (real) frequency. The program
illuminates the scattering geometry with plane waves traveling
in the positive z direction and (a) linearly polarized
with E-field in the x direction,
(b) linearly polarized with E-field in the y
direction, and (c) left-circularly polarized. For each
polarization, the program solves the BEM system, generates
a .pp file that may be opened in GMSH to
view the induced current distribution, and generates a
.out file containing RCS data that may be plotted
in gnuplot.
/*
* PlotRCS.cc -- an example program to demonstrate how to use the
* -- real-frequency algorithms in libRWG to compute the
* -- radar cross section for scattering of a plane wave
* -- from an arbitrary scattering geometry
*
* homer reid -- 12/2009
*/
#include
#include
#include
#include
#include
#include
#include
#include
#include
/***************************************************************/
/* constants ***************************************************/
/***************************************************************/
#define DTHETA 0.1 // angular increment at which to plot RCS
#define R_RCS 25.0 // radius at which to plot RCS
#define R_RCS2 (R_RCS*R_RCS)
/***************************************************************/
/* main function *********************************************/
/***************************************************************/
int main(int argc, char *argv[])
{
RWGGeometry *G;
HMatrix *M;
HVector *KN_XPol, *KN_YPol, *KN_CPol;
PlaneWaveData PWD;
int nThread;
double Omega, Theta, R[3], PV[3];
double RCS_XPol_Par, RCS_XPol_Perp;
double RCS_YPol_Par, RCS_YPol_Perp;
double RCS_CPol_Par, RCS_CPol_Perp;
FILE *f;
cdouble EH[6], *E, *H;
cdouble PM[1][6]; // if geometry had 3 objects, we would say PM[3][6] instead
/***************************************************************/
/* 1. create the geometry */
/***************************************************************/
if ( argc!=3 )
{ fprintf(stderr,"usage: %s MyGeometry.rwggeo Omega\n",argv[0]);
exit(1);
};
SetLogFileName("PlotRCS.log");
G=new RWGGeometry(argv[1]);
sscanf(argv[2],"%le",&Omega);
printf("Computing RCS for geometry %s at Omega=%g\n",G->GeoFileName,Omega);
nThread=GetNumProcs(); /* GetNumProcs is a routine in libhrutil*/
/***************************************************************/
/* 2. Initialize acceleration tables to speed up computation */
/* of matrix elements */
/***************************************************************/
G->PreCompute(nThread);
/***************************************************************/
/* 2. assemble and factorize the BEM matrix */
/* (and incidentally export it in MATLAB-readable format) */
/***************************************************************/
M=G->AllocateBEMMatrix(REAL_FREQ);
G->AssembleBEMMatrix(Omega, REAL_FREQ, nThread, M);
// M->ToMatlab("%s_%g",GetFileBase(argv[1]),Omega); /* export to matlab*/
M->LUFactorize();
/***************************************************************/
/* 3. assemble RHS vectors corresponding to plane waves */
/* traveling in the positive Z direction and */
/* (1) polarized with E-field parallel to the X axis */
/* (2) polarized with E-field parallel to the Y axis */
/* (3) circularly polarized (E \propto X + iY ) */
/***************************************************************/
PWD.nHat[0]=0.0; PWD.nHat[1]=0.0; PWD.nHat[2]=1.0;
PWD.Frequency=Omega;
PWD.RealFreq=REAL_FREQ;
PWD.Eps=PWD.Mu=1.0;
PWD.E0[0]=1.0; PWD.E0[1]=0.0; PWD.E0[2]=0.0;
KN_XPol=G->AllocateRHSVector(REAL_FREQ);
G->AssembleRHSVector(EHPlaneWave, (void *)&PWD, REAL_FREQ, nThread, KN_XPol);
PWD.E0[0]=0.0; PWD.E0[1]=1.0; PWD.E0[2]=0.0;
KN_YPol=G->AllocateRHSVector(REAL_FREQ);
G->AssembleRHSVector(EHPlaneWave, (void *)&PWD, REAL_FREQ, nThread, KN_YPol);
PWD.E0[0]=M_SQRT1_2; PWD.E0[1]=1.0fi*M_SQRT1_2; PWD.E0[2]=0.0;
KN_CPol=G->AllocateRHSVector(REAL_FREQ);
G->AssembleRHSVector(EHPlaneWave, (void *)&PWD, REAL_FREQ, nThread, KN_CPol);
/***************************************************************/
/* 4. solve the BEM systems for the three RHS vectors */
/***************************************************************/
M->LUSolve(KN_XPol);
M->LUSolve(KN_YPol);
M->LUSolve(KN_CPol);
/***************************************************************/
/* generate plots of induced surface currents and charges */
/***************************************************************/
G->PlotSurfaceCurrents(KN_XPol, Omega, REAL_FREQ, "%s_%g_K_XPol.pp",
GetFileBase(G->GeoFileName),Omega);
G->PlotSurfaceCurrents(KN_YPol, Omega, REAL_FREQ, "%s_%g_K_YPol.pp",
GetFileBase(G->GeoFileName),Omega);
G->PlotSurfaceCurrents(KN_CPol, Omega, REAL_FREQ, "%s_%g_K_CPol.pp",
GetFileBase(G->GeoFileName),Omega);
/***************************************************************/
/* print induced dipole moments *****************************/
/***************************************************************/
G->GetDipoleMoments(Omega, 1, KN_XPol, nThread, PM);
printf("XPol electric dipole moment: (%+7.3e,%7.3e)i + (%+7.3e,%7.3e)j + (%+7.3e,%7.3e)k\n",
__real__ PM[0][0], __imag__ PM[0][0],
__real__ PM[0][1], __imag__ PM[0][1],
__real__ PM[0][2], __imag__ PM[0][2]);
printf("XPol magnetic dipole moment: (%+7.3e,%7.3e)i + (%+7.3e,%7.3e)j + (%+7.3e,%7.3e)k\n",
__real__ PM[0][3], __imag__ PM[0][3],
__real__ PM[0][4], __imag__ PM[0][4],
__real__ PM[0][5], __imag__ PM[0][5]);
G->GetDipoleMoments(Omega, 1, KN_YPol, nThread, PM);
printf("YPol electric dipole moment: (%+7.3e,%7.3e)i + (%+7.3e,%7.3e)j + (%+7.3e,%7.3e)k\n",
__real__ PM[0][0], __imag__ PM[0][0],
__real__ PM[0][1], __imag__ PM[0][1],
__real__ PM[0][2], __imag__ PM[0][2]);
printf("YPol magnetic dipole moment: (%+7.3e,%7.3e)i + (%+7.3e,%7.3e)j + (%+7.3e,%7.3e)k\n",
__real__ PM[0][3], __imag__ PM[0][3],
__real__ PM[0][4], __imag__ PM[0][4],
__real__ PM[0][5], __imag__ PM[0][5]);
G->GetDipoleMoments(Omega, 1, KN_CPol, nThread, PM);
printf("CPol electric dipole moment: (%+7.3e,%7.3e)i + (%+7.3e,%7.3e)j + (%+7.3e,%7.3e)k\n",
__real__ PM[0][0], __imag__ PM[0][0],
__real__ PM[0][1], __imag__ PM[0][1],
__real__ PM[0][2], __imag__ PM[0][2]);
printf("CPol magnetic dipole moment: (%+7.3e,%7.3e)i + (%+7.3e,%7.3e)j + (%+7.3e,%7.3e)k\n",
__real__ PM[0][3], __imag__ PM[0][3],
__real__ PM[0][4], __imag__ PM[0][4],
__real__ PM[0][5], __imag__ PM[0][5]);
/***************************************************************/
/* 5. plot the RCS on a half-circle of radius R_RCS (=100) */
/* centered at the origin and lying in the XZ plane */
/***************************************************************/
f=vfopen("%s_%g.rcs","w",GetFileBase(argv[1]),Omega);
for(Theta=0.0; Theta<=M_PI; Theta+=DTHETA)
{
R[0]=R_RCS*sin(Theta);
R[1]=0.0;
R[2]=R_RCS*cos(Theta);
fprintf(f,"%e ",Theta);
G->GetFields(R, Omega, REAL_FREQ, KN_XPol, nThread, EH);
E=EH; H=EH+3;
RCS_XPol_Par = R_RCS2*( __real__ ( E[0]* ~E[0] ));
RCS_XPol_Perp= R_RCS2*( __real__ ( E[1]* ~E[1] ));
fprintf(f,"%e %e %e %e %e %e ",
CMag2(EH[0]),CMag2(EH[1]),CMag2(EH[2]),
CMag2(EH[3]),CMag2(EH[4]),CMag2(EH[5]));
G->GetFields(R, Omega, REAL_FREQ, KN_YPol, nThread, EH);
E=EH; H=EH+3;
RCS_YPol_Par = R_RCS2*( __real__ ( E[0]* ~E[0] ));
RCS_YPol_Perp= R_RCS2*( __real__ ( E[1]* ~E[1] ));
fprintf(f,"%e %e %e %e %e %e ",
CMag2(EH[0]),CMag2(EH[1]),CMag2(EH[2]),
CMag2(EH[3]),CMag2(EH[4]),CMag2(EH[5]));
G->GetFields(R, Omega, REAL_FREQ, KN_CPol, nThread, EH);
E=EH; H=EH+3;
RCS_CPol_Par = R_RCS2*( __real__ ( E[0]* ~E[0] ));
RCS_CPol_Perp= R_RCS2*( __real__ ( E[1]* ~E[1] ));
fprintf(f,"%e %e %e %e %e %e ",
CMag2(EH[0]),CMag2(EH[1]),CMag2(EH[2]),
CMag2(EH[3]),CMag2(EH[4]),CMag2(EH[5]));
fprintf(f,"\n");
};
fclose(f);
printf("\nThank you for your support.\n");
}
2.12 A complete imaginary-frequency example: Computing the
vacuum expectation value of the electromagnetic stress-energy
tensor
/*
* StressTensor.cc -- an example program to demonstrate how to use the
* -- imaginary-frequency algorithms of libRWG to
* -- compute the vacuum expectation value of the
* -- electromagnetic stress-energy tensor
*
* homer reid -- 11/2009
*/
#include
#include
#include
#include
#include
#include "libRWG.h"
#include "libIncField.h"
#include "libhrutil.h"
/***************************************************************/
/* prototypes for lapack routines ******************************/
/***************************************************************/
extern "C" {
int dpotrf_(const char *uplo, int *n, double *a, int *lda, int *info);
int dpotrs_(const char *uplo, int *n, int *nrhs, double *a, int *lda,
double *b, int *ldb, int *info);
}
/***************************************************************/
/* routine that evaluates the electromagnetic stress-energy */
/* tensor. */
/* inputs: */
/* G: pointer to RWGGeometry */
/* Beta: imaginary frequency */
/* M: BEM matrix as factorized by dpotrf */
/* R: evaluation point */
/* nThread: number of threads to spawn when assembling RHS */
/* outputs: */
/* TE[3][3]: electric stress-energy tensor (TE[i][j]=TE_{ij})*/
/* TM[3][3]: matgnetic stress-energy tensor (TM[i][j]=TM_{ij})*/
/***************************************************************/
void GetStressTensor(RWGGeometry *G, double Beta, double *M, double R[3],
int nThread, double TE[3][3], double TM[3][3])
{
static int NBF;
static double *X=0;
int iOne=1, info, i, j;
struct PointSourceData MyPSD, *PSD=&MyPSD;
double EH[6];
/* on first invocation only, allocate memory to store the RHS vector */
if (X==0)
{ NBF=G->GetDimension();
X=(double *)malloc(NBF*sizeof(double));
};
/* configure some properties of the point source */
memcpy(PSD->X0, R, 3*sizeof(double)); /* point source is located at R */
PSD->FreqType=LIF_IMAG_FREQ; /* point source has imag frequency */
PSD->Frequency=Beta; /* and that frequency is beta */
/*------------------------------------------------------------------*/
/*- step 1: compute all components of electric stress-energy tensor */
/*------------------------------------------------------------------*/
PSD->SourceType=LIF_TYPE_PSEC; /* for this, point source is a PSEC */
for(j=0; j<3; j++)
{
/* set point source to point in the j direction */
memset(PSD->nHat,0,3*sizeof(double));
PSD->nHat[j]=1.0;
/* form RHS vector using incident field of point source */
G->AssembleRHS_IF(EHPointSource, (void *)PSD, nThread, X);
/* solve the BEM system to compute current source BF weights */
dpotrs_("U", &NBF, &iOne, M, &NBF, X, &NBF, &info);
/* compute the scattered electric field back at X */
G->GetFields_IF(R, Beta, X, nThread, EH);
/* the three components of the electric field at X are the */
/* jth column of the TE matrix */
for(i=0; i<3; i++)
TE[i][j] = EH[i];
};
/*------------------------------------------------------------------*/
/* step 2: compute all components of magnetic stress-energy tensor */
/*------------------------------------------------------------------*/
PSD->SourceType=LIF_TYPE_PSMC; /* now the point source is a PSMC */
for(j=0; j<3; j++)
{
/* set point source to point in the j direction */
memset(PSD->nHat,0,3*sizeof(double));
PSD->nHat[j]=1.0;
/* form RHS vector using incident field of point source */
G->AssembleRHS_IF(EHPointSource, (void *)PSD, nThread, X);
/* solve the BEM system to compute current source BF weights */
dpotrs_("U", &NBF, &iOne, M, &NBF, X, &NBF, &info);
/* compute the scattered electric field back at X */
G->GetFields_IF(R, Beta, X, nThread, EH);
/* the three components of the magnetic field at X are the */
/* jth column of the TH matrix */
for(i=0; i<3; i++)
TM[i][j] = EH[i+3];
};
}
/***************************************************************/
/* main function *********************************************/
/***************************************************************/
int main()
{
RWGGeometry *G;
int i, NBF, nThread, info;
double Beta;
double *M, R[3], TE[3][3], TM[3][3];
/***************************************************************/
/* 1. create the geometry */
/***************************************************************/
G=new RWGGeometry("Plates.rwggeo");
Beta=0.1;
nThread=GetNumProcs(); /* GetNumProcs is a routine in libhrutil*/
/***************************************************************/
/* 2. assemble and factorize the BEM matrix */
/***************************************************************/
NBF=G->GetDimension();
M=(double *)malloc(NBF*NBF*sizeof(double));
G->AssembleM_IF(Beta,nThread,M);
dpotrf_("U", &NBF, M, &NBF, &info);
if (info) /* this should only happen in case of defective geometries */
ErrExit("BEM matrix was not positive definite");
/***************************************************************/
/* 3. get the stress tensors at R=(1,2,3) **********************/
/***************************************************************/
R[0]=1.0;
R[1]=2.0;
R[2]=3.0;
GetStressTensor(G, Beta, M, R, nThread, TE, TM);
/***************************************************************/
/* print results ***********************************************/
/***************************************************************/
printf("Electric stress-energy tensor: \n");
for(i=0; i<3; i++)
printf("%+7.3e %+7.3e %+7.3e\n",TE[i][0],TE[i][1],TE[i][2]);
printf("Magnetic stress-energy tensor: \n");
for(i=0; i<3; i++)
printf("%+7.3e %+7.3e %+7.3e\n",TM[i][0],TM[i][1],TM[i][2]);
printf("\nThank you for your support.\n");
}
3. Low-Level API
3.1 Creating an Object
3.2 Computing Matrix Elements
|
