meep Namespace Reference

Classes
class	Boundary_Region
struct	fieldop_data
struct	Dft_Chunk_Data
class	Fields
class	Fields_Chunk
class	Grace_Point
struct	h5_output_data
struct	rintegrand_data
class	H5file
	H5file.cpp: HDF5 file I/O. More...
struct	integrate_data
struct	rfun_wrap_data
class	Material_Function
class	Simple_Material_Function
class	Dft_Chunk
class	Dft_Flux
class	Flux_Vol
class	Polarizability_Identifier
class	Grace
class	Polarizability
class	Polarization
class	Bandsdata
	Yet another all public class. More...
class	initialize
class	Monitor_Point
struct	Src_Vol_Chunkloop_Data
class	Src_Time
	Time-dependence of a current source, intended to be overridden by subclasses. More...
class	Gaussian_Src_Time
	Gaussian-envelope source with given frequency, width, peak-time, cutoff I would guess that this provides a dipole source, with the following behavior: tt = time - peak_time. More...
class	Continuous_Src_Time
class	Custom_Src_Time
	user-specified source function with start and end times More...
class	Src_Vol
class	Structure_Chunk
class	Structure
struct	signed_direction
class	vec
	A general use class for storing and manipulating low rank vectors. More...
class	ivec
class	Geometric_Volume
class	Volume
class	Symmetry
class	Geometric_Volume_List
Typedefs
typedef double *	pdouble
typedef void(*	bicgstab_op )(const double *x, double *y, void *data)
typedef double(*	fx_func )(const vec &)
typedef void(*	field_chunkloop )(Fields_Chunk *fc, int ichunk, component cgrid, ivec is, ivec ie, vec s0, vec s1, vec e0, vec e1, double dV0, double dV1, ivec shift, complex< double > shift_phase, const Symmetry &S, int sn, void *chunkloop_data)
	field_chunkloop function pointer definition.
typedef complex< double >(*	field_function )(const complex< double > *Fields, const vec &loc, void *integrand_data_)
typedef double(*	field_rfunction )(const complex< double > *Fields, const vec &loc, void *integrand_data_)
typedef Structure_Chunk *	Structure_Chunk_ptr
Enumerations
enum	in_or_out { e_incoming = 0, e_outgoing }
enum	connect_phase { e_CONNECT_PHASE = 0, e_CONNECT_NEGATE = 1, e_CONNECT_COPY = 2 }
enum	boundary_condition { e_periodic = 0, e_metallic, e_magnetic, e_none }
enum	time_sink {   e_connecting, e_stepping, e_boundaries, e_mpiTime,   e_fieldOutput, e_fourierTransforming, e_other }
enum	Grace_type { e_XY, e_ERROR_BARS }
enum	component {   Ex = 0, Ey, Er, Ep,   Ez, Hx, Hy, Hr,   Hp, Hz, Dx, Dy,   Dr, Dp, Dz, Dielectric }
enum	derived_component {   Sx = 100, Sy, Sr, Sp,   Sz, EnergyDensity, D_EnergyDensity, H_EnergyDensity }
enum	ndim { D1 = 0, D2, D3, Dcyl }
enum	field_type { E_stuff = 0, H_stuff = 1, D_stuff = 2, P_stuff = 3 }
enum	boundary_side { High = 0, Low }
enum	direction {   X = 0, Y, Z, R,   P, NO_DIRECTION }
Functions
static vec	sphere_pt (const vec &cent, double R, int n, double &weight)
static void	anisoaverage (Material_Function &epsilon, const Geometric_Volume dV, component ec, double invepsrow[3], double tol, int maxeval)
int	do_harminv (complex< double > *data, int n, double dt, double fmin, double fmax, int maxbands, complex< double > *amps, double *freq_re, double *freq_im, double *errors, double spectral_density, double Q_thresh, double rel_err_thresh, double err_thresh, double rel_amp_thresh, double amp_thresh)
	Performs harmonic inversion of a given signal.
static double	dot (int n, const double *x, const double *y)
static double	norm2 (int n, const double *x)
static void	xpay (int n, double *x, double a, const double *y)
int	bicgstabL (const int L, const int n, double *x, bicgstab_op A, void *Adata, const double *b, const double tol, int *iters, double *work, const bool quiet)
Boundary_Region	pml (double thickness, direction d, boundary_side side)
Boundary_Region	pml (double thickness, direction d)
Boundary_Region	pml (double thickness)
static void	handle_control_c (int i)
void	deal_with_ctrl_c (int stop_now=2)
	Allows you to hit ctrl-C to tell your calculation to stop and clean up.
static void	Fields_to_array (const Fields &f, complex< double > *x)
static void	array_to_Fields (const complex< double > *x, Fields &f)
static void	fieldop (const double *xr, double *yr, void *data_)
static void	add_Dft_Chunkloop (meep::Fields_Chunk *fc, int ichunk, component cgrid, ivec is, ivec ie, vec s0, vec s1, vec e0, vec e1, double dV0, double dV1, ivec shift, complex< double > shift_phase, const Symmetry &S, int sn, void *chunkloop_data)
static int	Dft_Chunks_Ntotal (Dft_Chunk *Dft_Chunks, int *my_start)
void	save_dft_hdf5 (Dft_Chunk *Dft_Chunks, component c, H5file *file, const char *dprefix)
void	load_dft_hdf5 (Dft_Chunk *Dft_Chunks, component c, H5file *file, const char *dprefix)
static complex< double >	dot_integrand (const complex< double > *Fields, const vec &loc, void *data_)
static void	thermo_chunkloop (Fields_Chunk *fc, int ichunk, component cgrid, ivec is, ivec ie, vec s0, vec s1, vec e0, vec e1, double dV0, double dV1, ivec shift, complex< double > shift_phase, const Symmetry &S, int sn, void *sum_)
static complex< double >	dot3_max_integrand (const complex< double > *Fields, const vec &loc, void *data_)
static complex< double >	dot_fx_integrand (const complex< double > *Fields, const vec &loc, void *data_)
complex< double >	one (const vec &v)
field_rfunction	derived_component_func (derived_component c, const Volume &v, int &nFields, component cs[6])
static bool	cross_negative (direction a, direction b)
static direction	cross (direction a, direction b)
static bool	is_tm (component c)
static bool	is_like (ndim d, component c1, component c2)
static void	h5_findsize_chunkloop (Fields_Chunk *fc, int ichnk, component cgrid, ivec is, ivec ie, vec s0, vec s1, vec e0, vec e1, double dV0, double dV1, ivec shift, complex< double > shift_phase, const Symmetry &S, int sn, void *data_)
static void	h5_output_chunkloop (Fields_Chunk *fc, int ichnk, component cgrid, ivec is, ivec ie, vec s0, vec s1, vec e0, vec e1, double dV0, double dV1, ivec shift, complex< double > shift_phase, const Symmetry &S, int sn, void *data_)
static complex< double >	rintegrand_fun (const complex< double > *Fields, const vec &loc, void *data_)
static complex< double >	component_fun (const complex< double > *Fields, const vec &loc, void *data_)
static double	J (int m, double kr)
static double	Jprime (int m, double kr)
static double	Jroot (int m, int n)
static double	Jmax (int m, int n)
static complex< double >	JJ (const vec &v)
static complex< double >	JP (const vec &v)
static complex< double >	stupid_A (const vec &p)
static complex< double >	stupider_A (const vec &p)
static void	integrate_chunkloop (Fields_Chunk *fc, int ichunk, component cgrid, ivec is, ivec ie, vec s0, vec s1, vec e0, vec e1, double dV0, double dV1, ivec shift, complex< double > shift_phase, const Symmetry &S, int sn, void *data_)
static complex< double >	rfun_wrap (const complex< double > *Fields, const vec &loc, void *data_)
static complex< double >	return_the_field (const complex< double > *Fields, const vec &loc, void *integrand_data_)
static ivec	vec2diel_floor (const vec &v, double a, const ivec &equal_shift)
static ivec	vec2diel_ceil (const vec &v, double a, const ivec &equal_shift)
static int	iabs (int i)
const char *	make_output_directory (const char *exename, const char *jobname=NULL)
	Utility function to parse the executable name use it to come up with a directory name, avoiding overwriting any existing directory, unless the source file hasn't changed.
void	trash_output_directory (const char *dirname)
FILE *	create_output_file (const char *dirname, const char *fname)
double	max (double a, double b)
double	min (double a, double b)
int	max (int a, int b)
int	min (int a, int b)
static int	abs (int a)
static double	abs (double a)
static int	my_round (double x)
int	small_r_metal (int m)
int	rmin_bulk (int m)
Symmetry	r_to_minus_r_Symmetry (int m)
double	wall_time (void)
void	abort (const char *fmt,...)
void	send (int from, int to, double *data, int size)
void	broadcast (int from, double *data, int size)
void	broadcast (int from, char *data, int size)
void	broadcast (int from, complex< double > *data, int size)
void	broadcast (int from, int *data, int size)
complex< double >	broadcast (int from, complex< double > data)
double	broadcast (int from, double data)
int	broadcast (int from, int data)
bool	broadcast (int from, bool b)
double	max_to_master (double in)
double	max_to_all (double in)
int	max_to_all (int in)
ivec	max_to_all (const ivec &v)
double	sum_to_master (double in)
double	sum_to_all (double in)
void	sum_to_all (const double *in, double *out, int size)
long double	sum_to_all (long double in)
int	sum_to_all (int in)
int	partial_sum_to_all (int in)
complex< double >	sum_to_all (complex< double > in)
complex< long double >	sum_to_all (complex< long double > in)
bool	or_to_all (bool in)
bool	and_to_all (bool in)
void	all_wait ()
int	my_rank ()
int	count_processors ()
void	master_printf (const char *fmt,...)
void	debug_printf (const char *fmt,...)
void	master_fprintf (FILE *f, const char *fmt,...)
FILE *	master_fopen (const char *name, const char *mode)
void	master_fclose (FILE *f)
void	begin_critical_section (int tag)
void	end_critical_section (int tag)
int	am_master ()
static void	cp (const char *a, const char *b)
static bool	is_ok_dir (const char *dirname)
double	expi (int cmp, double x)
bool	Src_Times_equal (const Src_Time &t1, const Src_Time &t2)
	this function is necessary to make equality commutative .
void	Src_Vol_chunkloop (Fields_Chunk *fc, int ichunk, component c, ivec is, ivec ie, vec s0, vec s1, vec e0, vec e1, double dV0, double dV1, ivec shift, complex< double > shift_phase, const Symmetry &S, int sn, void *data_)
	Adding source volumes can be treated as a kind of "integration" problem, since we need to loop over all the chunks that intersect the source volume, with appropriate interpolation weights at the boundaries so that the integral of the current is fixed regardless of resolution.
static double	it (int cmp, double *(f[2]), int ind)
int	rstart_0 (const Volume &v, double m)
int	rstart_1 (const Volume &v, double m)
static const char *	ts2n (time_sink s)
static void	pt (double ts[], time_sink s)
double	calc_nonlinear_inveps (const double Dsqr, const double Di, const double inveps, const double chi2, const double chi3)
const char *	dimension_name (ndim dim)
const char *	direction_name (direction d)
const char *	component_name (component c)
const char *	component_name (derived_component c)
const char *	component_name (int c)
vec	min (const vec &vec1, const vec &vec2)
vec	max (const vec &vec1, const vec &vec2)
ivec	min (const ivec &ivec1, const ivec &ivec2)
ivec	max (const ivec &ivec1, const ivec &ivec2)
static direction	yucky_dir (ndim dim, int n)
static void	stupidsort (int *ind, double *w, int l)
static void	stupidsort (ivec *locs, double *w, int l)
Geometric_Volume	empty_volume (ndim dim)
Volume	volone (double zsize, double a)
Volume	voltwo (double xsize, double ysize, double a)
Volume	vol1d (double zsize, double a)
Volume	vol2d (double xsize, double ysize, double a)
Volume	vol3d (double xsize, double ysize, double zsize, double a)
Volume	volcyl (double rsize, double zsize, double a)
Symmetry	rotate4 (direction axis, const Volume &v)
Symmetry	rotate2 (direction axis, const Volume &v)
Symmetry	mirror (direction axis, const Volume &v)
Symmetry	r_to_minus_r_symmetry (double m)
Symmetry	identity ()
static double	poynting_fun (const complex< double > *Fields, const vec &loc, void *data_)
static double	energy_fun (const complex< double > *Fields, const vec &loc, void *data_)
int	number_of_directions (ndim dim)
direction	start_at_direction (ndim dim)
direction	stop_at_direction (ndim dim)
signed_direction	flip (signed_direction d)
bool	has_direction (ndim dim, direction d)
bool	coordinate_mismatch (ndim dim, direction d)
bool	is_electric (component c)
bool	is_magnetic (component c)
bool	is_D (component c)
bool	is_derived (int c)
bool	is_poynting (derived_component c)
bool	is_energydensity (derived_component c)
field_type	type (component c)
direction	component_direction (int c)
int	direction_component (int c, direction d)
direction	component_direction (component c)
direction	component_direction (derived_component c)
component	direction_component (component c, direction d)
derived_component	direction_component (derived_component c, direction d)
bool	coordinate_mismatch (ndim dim, component c)
bool	coordinate_mismatch (ndim dim, derived_component c)
vec	veccyl (double rr, double zz)
double	abs (const vec &v)
vec	zero_vec (ndim di)
vec	unit_vec (ndim di, direction d)
vec	clean_vec (const vec &v, double val_unused=0.0)
ivec	iveccyl (int xx, int yy)
ivec	zero_ivec (ndim di)
ivec	one_ivec (ndim di)
ivec	unit_ivec (ndim di, direction d)
Symmetry	r_to_minus_r_Symmetry (double m)
Variables
static const int	num_sphere_quad [3] = { 2, 12, 50 }
static const double	sphere_quad [3][50][4]
int	interrupt = 0
	Incremented when ctrl-C is called.(starts with a zero).
static int	kill_time = 2
const double	infinity = 1e20
	this should be big enough for us
const double	nan = -1.23454321e123
	ideally, a value never encountered in practice
const char	Grace_header []
static double	ktrans
static double	kax
static int	m_for_J
static complex< double >(*	A_static )(component, const vec &)
static component	c_static
static Fields *	f_static
static complex< double >	prefactor_static
bool	quiet = false
	if true, suppress all non-error messages from Meep
const double	pi = 3.141592653589793238462643383276
const int	num_bandpts = 32
static FILE *	debf = NULL
const char	symlink_name [] = "latest_output"
const double	Cmax = 0.5
const int	NUM_FIELD_COMPONENTS = 15
const int	NUM_FIELD_TYPES = 4
vec	zero_vec (ndim)
ivec	zero_ivec (ndim)
ivec	one_ivec (ndim)

---

Typedef Documentation

typedef void(* meep::bicgstab_op)(const double *x, double *y, void *data)

typedef void(* meep::field_chunkloop)(Fields_Chunk *fc, int ichunk, component cgrid, ivec is, ivec ie, vec s0, vec s1, vec e0, vec e1, double dV0, double dV1, ivec shift, complex< double > shift_phase, const Symmetry &S, int sn, void *chunkloop_data)

field_chunkloop function pointer definition.

Parameters:

	fc	pointer to the Fields_Chunk to be operated on.
	ichun	chunk index?
	cgrid	?
	is	?
	ie	?
	e0	?
	e1	?
	dV0	?
	dV1	?
	shift	?
	shift_phase	?
	S	a symmetry.
	sn	?
	chunkloop_data	Data needed for the operation presumably.

typedef complex<double>(* meep::field_function)(const complex< double > *Fields, const vec &loc, void *integrand_data_)

typedef double(* meep::field_rfunction)(const complex< double > *Fields, const vec &loc, void *integrand_data_)

typedef double(* meep::fx_func)(const vec &)

typedef double* meep::pdouble

typedef Structure_Chunk* meep::Structure_Chunk_ptr

---

Enumeration Type Documentation

enum meep::boundary_condition

Enumerator:

e_periodic
e_metallic
e_magnetic
e_none

{ e_periodic=0, e_metallic, e_magnetic, e_none };

enum meep::boundary_side

Enumerator:

High
Low

{ High=0, Low };

enum meep::component

Enumerator:

Ex
Ey
Er
Ep
Ez
Hx
Hy
Hr
Hp
Hz
Dx
Dy
Dr
Dp
Dz
Dielectric

               { Ex=0, Ey, Er, Ep, Ez, Hx, Hy, Hr, Hp, Hz,
                 Dx, Dy, Dr, Dp, Dz, Dielectric };

enum meep::connect_phase

Enumerator:

e_CONNECT_PHASE
e_CONNECT_NEGATE
e_CONNECT_COPY

{ e_CONNECT_PHASE = 0, e_CONNECT_NEGATE=1, e_CONNECT_COPY=2 };

enum meep::derived_component

Enumerator:

Sx
Sy
Sr
Sp
Sz
EnergyDensity
D_EnergyDensity
H_EnergyDensity

                       { Sx=100, Sy, Sr, Sp, Sz, EnergyDensity,
                         D_EnergyDensity, H_EnergyDensity };

enum meep::direction

Enumerator:

X
Y
Z
R
P
NO_DIRECTION

{ X=0,Y,Z,R,P, NO_DIRECTION };

enum meep::field_type

Enumerator:

E_stuff
H_stuff
D_stuff
P_stuff

{ E_stuff=0, H_stuff=1, D_stuff=2, P_stuff=3 };

enum meep::Grace_type

Enumerator:

e_XY
e_ERROR_BARS

{ e_XY, e_ERROR_BARS };

enum meep::in_or_out

Enumerator:

e_incoming
e_outgoing

{ e_incoming=0, e_outgoing };

enum meep::ndim

Enumerator:

D1
D2
D3
Dcyl

{ D1=0, D2, D3, Dcyl };

enum meep::time_sink

Enumerator:

e_connecting
e_stepping
e_boundaries
e_mpiTime
e_fieldOutput
e_fourierTransforming
e_other

                 { e_connecting, e_stepping, e_boundaries, e_mpiTime,
                   e_fieldOutput, e_fourierTransforming, e_other };

---

Function Documentation

void meep::abort	(	const char *	fmt,
			...
	)

                                 {
  va_list ap;
  va_start(ap, fmt);
  fprintf(stderr, "meep: ");
  vfprintf(stderr, fmt, ap);
  va_end(ap);
  if (fmt[strlen(fmt) - 1] != '\n') fputc('\n', stderr); // force newline
#ifdef HAVE_MPI
  MPI_Abort(MPI_COMM_WORLD, 1);
#endif
  exit(1);
}

double meep::abs

(

const vec &

v

)

[inline]

{ return sqrt(v & v); }

static double meep::abs

(

double

a

)

[inline, static]

{ return fabs(a); }

static int meep::abs

(

int

a

)

[inline, static]

{ return a < 0 ? -a : a; }

static void meep::add_Dft_Chunkloop	(	meep::Fields_Chunk *	fc,
		int	ichunk,
		component	cgrid,
		ivec	is,
		ivec	ie,
		vec	s0,
		vec	s1,
		vec	e0,
		vec	e1,
		double	dV0,
		double	dV1,
		ivec	shift,
		complex< double >	shift_phase,
		const Symmetry &	S,
		int	sn,
		void *	chunkloop_data
	)			[static]

  {
        Dft_Chunk_Data *data = (Dft_Chunk_Data *) chunkloop_data;
        (void) shift; (void) ichunk; // unused

        if (cgrid != Dielectric) abort("dft chunks should use the Dielectric grid");
  
        component c = S.transform(data->c, -sn);
        if (c >= NUM_FIELD_COMPONENTS || !fc->m_f[c][0])
          return; // this chunk doesn't have component c

        data->Dft_Chunks = new Dft_Chunk(fc,is,ie,s0,s1,e0,e1,dV0,dV1,
                                                                         shift_phase * S.phase_shift(c, sn),
                                                                         c, chunkloop_data);
  }

void void meep::all_wait

(

)

                {
#ifdef HAVE_MPI
  MPI_Barrier(MPI_COMM_WORLD);
#endif
}

int meep::am_master

(

)

[inline]

{ return my_rank() == 0; }

bool meep::and_to_all

(

bool

in

)

                         {
  int out = in;
#ifdef HAVE_MPI
  MPI_Allreduce(&in,&out,1,MPI_INT,MPI_LAND,MPI_COMM_WORLD);
#endif
  return (bool) out;
}

static void meep::anisoaverage	(	Material_Function &	epsilon,
		const Geometric_Volume	dV,
		component	ec,
		double	invepsrow[3],
		double	tol,
		int	maxeval
	)			[static]

                                 {
          IVEC_LOOP_LOC(m_v, here);
          m_inveps[c][c_d][i] = 1/epsilon.eps(here);
        }
#endif
      }
  } else {
    const double smoothing_diameter = 1.0; // FIXME: make this user-changable?

    // may take a long time in 3d, so prepare to print status messages
    int npixels = 0, ipixel = 0;
    FOR_ELECTRIC_COMPONENTS(c) if (m_v.has_field(c)) {
      int loop_npixels = 0;
      LOOP_OVER_VOL(m_v, c, i) {
        loop_npixels = loop_n1 * loop_n2 * loop_n3;
        goto breakout;
      }
    breakout: npixels += loop_npixels;
    }
    double last_output_time = wall_time();

    FOR_ELECTRIC_COMPONENTS(c)
      if (m_v.has_field(c)) {
        FOR_ELECTRIC_COMPONENTS(c2) if (m_v.has_field(c2)) {
          direction d = component_direction(c2);
          if (!m_inveps[c][d]) m_inveps[c][d] = new double[m_v.ntot()];
          if (!m_inveps[c][d]) abort("Memory allocation error.\n");
        }

static void meep::array_to_Fields	(	const complex< double > *	x,
		Fields &	f
	)			[static]

{
  int ix = 0;
  for (int i=0;i<f.m_num_chunks;i++)
    if (f.m_chunks[i]->is_mine())
      FOR_COMPONENTS(c)
        if (f.m_chunks[i]->m_f[c][0] && (is_D(c) || is_magnetic(c))) {
          double *fs[3][2];
          DOCMP2 { fs[0][cmp] = f.m_chunks[i]->m_f[c][cmp];
                   fs[1][cmp] = f.m_chunks[i]->m_f_p_pml[c][cmp];
                   fs[2][cmp] = f.m_chunks[i]->m_f_m_pml[c][cmp]; }
          for (int k = 0; k < 3; ++k)
            if (fs[k][0]) {
              LOOP_OVER_VOL_OWNED(f.m_chunks[i]->m_vol, c, idx) {
                fs[k][0][idx] = real(x[ix]);
                fs[k][1][idx] = imag(x[ix]);
                ++ix;
              }
            }
        }
  f.force_consistency(H_stuff);
  f.force_consistency(D_stuff);
}

void meep::begin_critical_section

(

int

tag

)

{
#ifdef HAVE_MPI
     int process_rank;
     MPI_Comm_rank(MPI_COMM_WORLD, &process_rank);
     if (process_rank > 0) { /* wait for a message before continuing */
          MPI_Status status;
          int recv_tag = tag - 1; /* initialize to wrong value */
          MPI_Recv(&recv_tag, 1, MPI_INT, process_rank - 1, tag, 
                   MPI_COMM_WORLD, &status);
          if (recv_tag != tag) abort("invalid tag received in begin_critical_section");
     }
#else
     UNUSED(tag);
#endif
}

int meep::bicgstabL	(	const int	L,
		const int	n,
		double *	x,
		bicgstab_op	A,
		void *	Adata,
		const double *	b,
		const double	tol,
		int *	iters,
		double *	work,
		const bool	quiet
	)

{
  if (!work) return (2*L+3)*n; // required workspace

  pdouble *r = new pdouble[L+1];
  pdouble *u = new pdouble[L+1];
  for (int i = 0; i <= L; ++i) {
    r[i] = work + i * n;
    u[i] = work + (L+1 + i) * n;
  }

  double bnrm = norm2(n, b);
  if (bnrm == 0.0) bnrm = 1.0;
  
  int iter = 0;
  double last_output_wall_time = wall_time();

  double *gamma = new double[L + 1];
  double *gamma_p = new double[L + 1];
  double *gamma_pp = new double[L + 1];

  double *tau = new double[L * L];
  double *sigma = new double[L + 1];

  int ierr = 0; // error code to return, if any
  const double breaktol = 1e-30;

  /**** FIXME: check for breakdown conditions(?) during iteration  ****/

  // rtilde = r[0] = b - Ax
  double *rtilde = work + (2*L+2) * n;
  A(x, r[0], Adata);
  for (int m = 0; m < n; ++m) rtilde[m] = r[0][m] = b[m] - r[0][m];

  { /* Sleipjen normalizes rtilde in his code; it seems to help slightly */
    double s = 1.0 / norm2(n, rtilde);
    for (int m = 0; m < n; ++m) rtilde[m] *= s;
  }

  memset(u[0], 0, sizeof(double) * n); // u[0] = 0

  double rho = 1.0, alpha = 0, omega = 1;

  double resid;
  while ((resid = norm2(n, r[0])) > tol * bnrm) {
    ++iter;
    if (!quiet && wall_time() > last_output_wall_time + MIN_OUTPUT_TIME) {
      master_printf("residual[%d] = %g\n", iter, resid / bnrm);
      last_output_wall_time = wall_time();
    }

    rho = -omega * rho;
    for (int j = 0; j < L; ++j) {
      if (fabs(rho) < breaktol) { ierr = -1; goto finish; }
      double rho1 = dot(n, r[j], rtilde);
      double beta = alpha * rho1 / rho;
      rho = rho1;
      for (int i = 0; i <= j; ++i)
        for (int m = 0; m < n; ++m) u[i][m] = r[i][m] - beta * u[i][m];
      A(u[j], u[j+1], Adata);
      alpha = rho / dot(n, u[j+1], rtilde);
      for (int i = 0; i <= j; ++i)
        xpay(n, r[i], -alpha, u[i+1]);
      A(r[j], r[j+1], Adata);
      xpay(n, x, alpha, u[0]);
    }

    for (int j = 1; j <= L; ++j) {
      for (int i = 1; i < j; ++i) {
        int ij = (j-1)*L + (i-1);
        tau[ij] = dot(n, r[j], r[i]) / sigma[i];
        xpay(n, r[j], -tau[ij], r[i]);
      }
      sigma[j] = dot(n, r[j],r[j]);
      gamma_p[j] = dot(n, r[0], r[j]) / sigma[j];
    }

    omega = gamma[L] = gamma_p[L];
    for (int j = L-1; j >= 1; --j) {
      gamma[j] = gamma_p[j];
      for (int i = j+1; i <= L; ++i)
        gamma[j] -= tau[(i-1)*L + (j-1)] * gamma[i];
    }
    for (int j = 1; j < L; ++j) {
      gamma_pp[j] = gamma[j+1];
      for (int i = j+1; i < L; ++i)
        gamma_pp[j] += tau[(i-1)*L + (j-1)] * gamma[i+1];
    }
    
    xpay(n, x, gamma[1], r[0]);
    xpay(n, r[0], -gamma_p[L], r[L]);
    xpay(n, u[0], -gamma[L], u[L]);
    for (int j = 1; j < L; ++j) { /* TODO: use blas DGEMV (for L > 2) */
      xpay(n, x, gamma_pp[j], r[j]);
      xpay(n, r[0], -gamma_p[j], r[j]);
      xpay(n, u[0], -gamma[j], u[j]);
    }

    if (iter == *iters) { ierr = 1; break; }
  }

  if (!quiet) master_printf("final residual = %g\n", norm2(n, r[0]) / bnrm);

 finish:
  delete[] sigma;
  delete[] tau;
  delete[] gamma_pp;
  delete[] gamma_p;
  delete[] gamma;
  delete[] u;
  delete[] r;

  *iters = iter;
  return ierr;
}

bool meep::broadcast	(	int	from,
		bool	b
	)

                                 {
  return broadcast(from, (int) b);
}

int meep::broadcast	(	int	from,
		int	data
	)

                                  {
#ifdef HAVE_MPI
  MPI_Bcast(&data, 1, MPI_INT, from, MPI_COMM_WORLD);
#else
  UNUSED(from);
#endif
  return data;
}

double meep::broadcast	(	int	from,
		double	data
	)

                                        {
#ifdef HAVE_MPI
  MPI_Bcast(&data, 1, MPI_DOUBLE, from, MPI_COMM_WORLD);
#else
  UNUSED(from);
#endif
  return data;
}

complex< double > meep::broadcast	(	int	from,
		complex< double >	data
	)

                                                          {
#ifdef HAVE_MPI
  MPI_Bcast(&data, 2, MPI_DOUBLE, from, MPI_COMM_WORLD);
#else
  UNUSED(from);
#endif
  return data;
}

void meep::broadcast	(	int	from,
		int *	data,
		int	size
	)

                                              {
#ifdef HAVE_MPI
  if (size == 0) return;
  MPI_Bcast(data, size, MPI_INT, from, MPI_COMM_WORLD);
#else
  UNUSED(from);
  UNUSED(data);
  UNUSED(size);
#endif
}

void meep::broadcast	(	int	from,
		complex< double > *	data,
		int	size
	)

                                                          {
#ifdef HAVE_MPI
  if (size == 0) return;
  MPI_Bcast(data, 2*size, MPI_DOUBLE, from, MPI_COMM_WORLD);
#else
  UNUSED(from);
  UNUSED(data);
  UNUSED(size);
#endif
}

void meep::broadcast	(	int	from,
		char *	data,
		int	size
	)

                                               {
#ifdef HAVE_MPI
  if (size == 0) return;
  MPI_Bcast(data, size, MPI_CHAR, from, MPI_COMM_WORLD);
#else
  UNUSED(from);
  UNUSED(data);
  UNUSED(size);
#endif
}

void meep::broadcast	(	int	from,
		double *	data,
		int	size
	)

                                                 {
#ifdef HAVE_MPI
  if (size == 0) return;
  MPI_Bcast(data, size, MPI_DOUBLE, from, MPI_COMM_WORLD);
#else
  UNUSED(from);
  UNUSED(data);
  UNUSED(size);
#endif
}

double meep::calc_nonlinear_inveps	(	const double	Dsqr,
		const double	Di,
		const double	inveps,
		const double	chi2,
		const double	chi3
	)			[inline]

Given Dsqr = |D|^2 and Di = component of D, compute the factor f so that Ei = inveps * f * Di. In principle, this would involve solving a cubic equation, but instead we use a Pade approximant that is accurate to several orders. This is inaccurate if the nonlinear index change is large, of course, but in that case the chi2/chi3 power-series expansion isn't accurate anyway, so the cubic isn't physical there either.

                                                                          {
  double c2 = Di*chi2*(inveps*inveps);
  double c3 = Dsqr*chi3*(inveps*(inveps*inveps));
  return (1 + c2 + 2*c3)/(1 + 2*c2 + 3*c3);
}

vec meep::clean_vec	(	const vec &	v,
		double	val_unused = 0.0
	)			[inline]

                                                            {
  vec vc(v.dim, val_unused);
  LOOP_OVER_DIRECTIONS(v.dim, d) vc.set_direction(d, v.in_direction(d));
  return vc;
}

direction meep::component_direction

(

derived_component

c

)

[inline]

                                                          {
  switch (c) {
  case Sx: return X;
  case Sy: return Y;
  case Sz: return Z;
  case Sr: return R;
  case Sp: return P;
  case EnergyDensity: case D_EnergyDensity: case H_EnergyDensity:
    return NO_DIRECTION;
  }
  return X; // This code is never reached...
}

direction meep::component_direction

(

component

c

)

[inline]

                                                  {
  switch (c) {
  case Ex: case Hx: case Dx: return X;
  case Ey: case Hy: case Dy: return Y;
  case Ez: case Hz: case Dz: return Z;
  case Er: case Hr: case Dr: return R;
  case Ep: case Hp: case Dp: return P;
  case Dielectric: return NO_DIRECTION;
  }
  return X; // This code is never reached...
}

direction meep::component_direction

(

int

c

)

[inline]

                                            {
  if (is_derived(c))
    return component_direction(derived_component(c));
  else
    return component_direction(component(c));
}

static complex<double> meep::component_fun	(	const complex< double > *	Fields,
		const vec &	loc,
		void *	data_
	)			[static]

{
     (void) loc; // unused
     (void) data_; // unused
     return Fields[0];
}

const char * meep::component_name

(

int

c

)

                                  {
  return (is_derived(c) ? component_name(derived_component(c))
          : component_name(component(c)));
}

const char * meep::component_name

(

derived_component

c

)

                                                {
  if (!is_derived(int(c))) return component_name(component(c));
  switch (c) {
  case Sr: return "sr";
  case Sp: return "sp";
  case Sz: return "sz";
  case Sx: return "sx";
  case Sy: return "sy";
  case EnergyDensity: return "energy";
  case D_EnergyDensity: return "denergy";
  case H_EnergyDensity: return "henergy";
  }
  return "Error in component_name";
}

const char * meep::component_name

(

component

c

)

                                        {
  if (is_derived(int(c))) return component_name(derived_component(c));
  switch (c) {
  case Er: return "er";
  case Ep: return "ep";
  case Ez: return "ez";
  case Hr: return "hr";
  case Hp: return "hp";
  case Hz: return "hz";
  case Ex: return "ex";
  case Ey: return "ey";
  case Hx: return "hx";
  case Hy: return "hy";
  case Dx: return "dx";
  case Dy: return "dy";
  case Dz: return "dz";
  case Dr: return "dr";
  case Dp: return "dp";
  case Dielectric: return "eps";
  }
  return "Error in component_name";
}

bool meep::coordinate_mismatch	(	ndim	dim,
		derived_component	c
	)			[inline]

                                                               {
  return coordinate_mismatch(dim, component_direction(c));
}

bool meep::coordinate_mismatch	(	ndim	dim,
		component	c
	)			[inline]

                                                       {
  return coordinate_mismatch(dim, component_direction(c));
}

bool meep::coordinate_mismatch	(	ndim	dim,
		direction	d
	)			[inline]

                                                       {
  return (d != NO_DIRECTION &&
          ((dim >= D1 && dim <= D3 && d != X && d != Y && d != Z) ||
           (dim == Dcyl && d != R && d != P && d != Z)));
}

int meep::count_processors

(

)

                       {
#ifdef HAVE_MPI
  int n;
  MPI_Comm_size(MPI_COMM_WORLD, &n);
  return n;
#else
  return 1;
#endif
}

static void meep::cp	(	const char *	a,
		const char *	b
	)			[static]

                                             {
  FILE *fa = fopen(a,"r");
  FILE *fb = fopen(b,"w");
  if (!fa || !fb) return;
  int ca;
  while (1) {
    ca = getc(fa);
    if (ca == EOF) break;
    putc(ca,fb);
  }
  fclose(fa);
  fclose(fb);
}

FILE * meep::create_output_file	(	const char *	dirname,
		const char *	fname
	)

                                                                 {
  const int buflen = 300;
  char n[buflen];
  snprintf(n, buflen, "%s/%s", dirname, fname);
  FILE *o = master_fopen(n, "w");
  if (!o) abort("Unable to create file %s!\n", n);
  return o;
}

static direction meep::cross	(	direction	a,
		direction	b
	)			[inline, static]

                                                          {
    if (a < R && b < R) return (direction)((3+2*a-b)%3);
    return (direction) (2 + (3+2*(a-2)-(b-2))%3);
  }

static bool meep::cross_negative	(	direction	a,
		direction	b
	)			[inline, static]

                                                              {
    return ((3+b-a)%3) == 2;
  }

void meep::deal_with_ctrl_c

(

int

stop_now

)

Allows you to hit ctrl-C to tell your calculation to stop and clean up.

                                    {
  kill_time = stop_now;
  if (signal(SIGINT, handle_control_c) == SIG_IGN)
    signal(SIGINT, SIG_IGN); // ignore if parent process was ignoring
}

void void meep::debug_printf	(	const char *	fmt,
			...
	)

                                        {
  va_list ap;
  va_start(ap, fmt);
  if (debf == NULL) {
    char temp[50];
    snprintf(temp, 50, "debug_out_%d", my_rank());
    debf = fopen(temp,"w");
    if (!debf) abort("Unable to open debug output %s\n", temp);
  }
  vfprintf(debf, fmt, ap);
  fflush(debf);
  va_end(ap);
}

field_rfunction meep::derived_component_func	(	derived_component	c,
		const Volume &	v,
		int &	nFields,
		component	cs[6]
	)

                                                                      {
  switch (c) {
  case Sx: case Sy: case Sz: case Sr: case Sp:
    switch (c) {
    case Sx: cs[0] = Ey; cs[1] = Hz; break;
    case Sy: cs[0] = Ez; cs[1] = Hx; break;
    case Sz: cs[0] = Ex; cs[1] = Hy; break;
    case Sr: cs[0] = Ep; cs[1] = Hz; break;
    case Sp: cs[0] = Ez; cs[1] = Hr; break;
    default: break; // never reached
    }
    nFields = 4;
    cs[2] = direction_component(Ex, component_direction(cs[1]));
    cs[3] = direction_component(Hx, component_direction(cs[0]));
    return poynting_fun;

  case EnergyDensity: case D_EnergyDensity: case H_EnergyDensity:
    nFields = 0;
    if (c != H_EnergyDensity)
      FOR_ELECTRIC_COMPONENTS(c0) if (v.has_field(c0)) {
        cs[nFields++] = c0;
        cs[nFields++] = direction_component(Dx, component_direction(c0));
      }
    if (c != D_EnergyDensity)
      FOR_MAGNETIC_COMPONENTS(c0) if (v.has_field(c0)) {
        cs[nFields++] = c0;
        cs[nFields++] = c0;
      }
    if (nFields > 6) abort("too many field components");
    return energy_fun;

  default: 
    abort("unknown derived_component in derived_component_func");
  }
  return 0;
}

static int meep::Dft_Chunks_Ntotal	(	Dft_Chunk *	Dft_Chunks,
		int *	my_start
	)			[static]

                                                                     {
        int n = 0;
        for (Dft_Chunk *cur = Dft_Chunks; cur; cur = cur->m_next_in_dft)
          n += cur->m_N * cur->m_Nomega * 2;
        *my_start = partial_sum_to_all(n) - n; // sum(n) for processes before this
        return sum_to_all(n);
  }

const char * meep::dimension_name

(

ndim

dim

)

                                     {
  switch (dim) {
  case D1: return "1D";
  case D2: return "2D";
  case D3: return "3D";
  case Dcyl: return "Cylindrical";
  }
  return "Error in dimension_name";
}

derived_component meep::direction_component	(	derived_component	c,
		direction	d
	)			[inline]

                                                                               {
  derived_component start_point;
  if (is_poynting(c)) start_point = Sx;
  else if (is_energydensity(c) && d == NO_DIRECTION) return c;
  else abort("unknown field component %d", c);
  switch (d) {
  case X: return start_point;
  case Y: return (derived_component) (start_point + 1);
  case Z: return (derived_component) (start_point + 4);
  case R: return (derived_component) (start_point + 2);
  case P: return (derived_component) (start_point + 3);
  case NO_DIRECTION: abort("vector %d derived_component in NO_DIRECTION", c);
  }
  return Sx; // This is never reached.
}

component meep::direction_component	(	component	c,
		direction	d
	)			[inline]

                                                               {
  component start_point;
  if (is_electric(c)) start_point = Ex;
  else if (is_magnetic(c)) start_point = Hx;
  else if (is_D(c)) start_point = Dx;
  else if (c == Dielectric && d == NO_DIRECTION) return Dielectric;
  else abort("unknown field component %d", c);
  switch (d) {
  case X: return start_point;
  case Y: return (component) (start_point + 1);
  case Z: return (component) (start_point + 4);
  case R: return (component) (start_point + 2);
  case P: return (component) (start_point + 3);
  case NO_DIRECTION: abort("vector %d component in NO_DIRECTION", c);
  }
  return Ex; // This is never reached.
}

int meep::direction_component	(	int	c,
		direction	d
	)			[inline]

                                                   {
  if (is_derived(c))
    return int(direction_component(derived_component(c), d));
  else
    return int(direction_component(component(c), d));
}

const char * meep::direction_name

(

direction

d

)

                                        {
  switch (d) {
  case X: return "x";
  case Y: return "y";
  case Z: return "z";
  case R: return "r";
  case P: return "phi";
  case NO_DIRECTION: return "no_direction";
  }
  return "Error in direction_name";
}

int meep::do_harminv	(	complex< double > *	data,
		int	n,
		double	dt,
		double	fmin,
		double	fmax,
		int	maxbands,
		complex< double > *	amps,
		double *	freq_re,
		double *	freq_im,
		double *	errors,
		double	spectral_density,
		double	Q_thresh,
		double	rel_err_thresh,
		double	err_thresh,
		double	rel_amp_thresh,
		double	amp_thresh
	)

Performs harmonic inversion of a given signal.

Perform harmonic inversion.

See harminv documentation for detail.

Parameters:

	data	Complex data for harmonic inversion.
	n	The size of the data for inversion.
	dt	Size of the timestep-> gets you units right.
	fmin	Minimum frequency.
	fmax	Maximum frequency.
	maxbands	Maximum number of modes to find.
	amps	Complex amplitudes.
	freq_re	Pointer to real frequency of modes.
	freq_im	Pointer to imaginary frequencies.
	errors	Pointer to errors array.
	spectral_density	No idea.
	Q_thresh	No idea.
	rel_err_thress	No idea.
	err_thresh	No idea.
	rel_amp_thresh	No idea.
	amp_thres	No idea.

  {
#ifndef HAVE_HARMINV
    abort("compiled without Harminv library, required for do_harminv");
    return 0;
#else
    int numfreqs = int(fabs(fmax-fmin)*dt*n*spectral_density); // c.f. harminv
    if (numfreqs > 150) numfreqs = 150; // prevent matrices from getting too big
    if (numfreqs < 2) numfreqs = 2;

    if (maxbands > numfreqs)
      numfreqs = maxbands;

    // check for all zeros in input
    // data is a size n array.
    {
      int i;
      for (i = 0; i < n && data[i] == 0.0; i++)
        ;
      if (i == n)
        return 0;
    }

#if 0
    // debugging: save data file and arguments for standalone harminv program
    {
      FILE *f = fopen("harminv.dat", "w");
      fprintf(f, "# -f %d -t %g %g-%g -Q %e -e %e -E %e -a %e -A %e -F\n",
              numfreqs, dt, fmin, fmax,
              Q_thresh, rel_err_thresh, err_thresh, rel_amp_thresh, amp_thresh);
      for (int i = 0; i < n; ++i)
        fprintf(f, "%g%+gi\n", real(data[i]), imag(data[i]));
      fclose(f);
    }
#endif

    harminv_data hd = harminv_data_create(n, data, fmin*dt, fmax*dt, numfreqs);
    harminv_solve(hd);

    int nf = harminv_get_num_freqs(hd);
    if (nf == 0) return 0;
    int *fsort = new int[nf]; // indices of frequencies, sorted as needed
  
    for (int i = 0; i < nf; ++i)
      fsort[i] = i;
    for (int i = 0; i < nf; ++i) // sort in increasing order of error
      for (int j = i + 1; j < nf; ++j) 
        if (harminv_get_freq_error(hd, fsort[i]) >
            harminv_get_freq_error(hd, fsort[j])) {
          int k = fsort[i];
          fsort[i] = fsort[j];
          fsort[j] = k;
        }
  
    double min_err = harminv_get_freq_error(hd, fsort[0]);
    double max_amp = abs(harminv_get_amplitude(hd, 0));
    for (int i = 1; i < nf; ++i) {
      double amp = abs(harminv_get_amplitude(hd, i));
      if (max_amp < amp)
        max_amp = amp;
    }
    { // eliminate modes that fall outside the various thresholds:
      int j = 0;
      for (int i = 0; i < nf; ++i) {
        double f = abs(harminv_get_freq(hd, fsort[i]) / dt);
        double err = harminv_get_freq_error(hd, fsort[i]);
        double amp = abs(harminv_get_amplitude(hd, fsort[i]));
        if (f >= fmin && f <= fmax
            && abs(harminv_get_Q(hd, fsort[i])) > Q_thresh
            && err < err_thresh
            && err < rel_err_thresh * min_err
            && amp > amp_thresh
            && amp > rel_amp_thresh * max_amp) {
          fsort[j++] = fsort[i];
        }
      }
      nf = j;
    }
    { // eliminate positive/negative frequency pairs
      // set indices to -1 for frequencies to be eliminated
      for (int i = 0; i < nf; ++i) 
        if (fsort[i] != -1) { // i hasn't been eliminated yet
          double f = harminv_get_freq(hd, fsort[i]);
          if (f < 0.0) {
            double kdiff = -2 * f;
            int kpos = i;
            for (int k = 0; k < nf; ++k) // search for closest positive freq.
              if (fsort[k] != -1) { // k hasn't been eliminated yet
                double fdiff = abs(harminv_get_freq(hd, fsort[k]) + f);
                if (fdiff < kdiff) {
                  kpos = k;
                  kdiff = fdiff;
                }
              }
            if (kpos != i && kdiff < 2.0 / n) { // consider them the same
              // pick the one with the smaller error
              if (harminv_get_freq_error(hd, fsort[i]) <
                  harminv_get_freq_error(hd, fsort[kpos]))
                fsort[kpos] = -1;
              else
                fsort[i] = -1;
            }
          }
        }
      int j = 0;
      for (int i = 0; i < nf; ++i) // remove the eliminated indices
        if (fsort[i] != -1)
          fsort[j++] = fsort[i];
      nf = j;
    }
  
    if (nf > maxbands)
      nf = maxbands;
  
    // sort again, this time in increasing order of freq:
    for (int i = 0; i < nf; ++i) // simple O(nf^2) sort
      for (int j = i + 1; j < nf; ++j) 
        if (abs(harminv_get_freq(hd, fsort[i])) >
            abs(harminv_get_freq(hd, fsort[j]))) {
          int k = fsort[i];
          fsort[i] = fsort[j];
          fsort[j] = k;
        }
  
    for (int i = 0; i < nf; ++i) {
      complex<double> freq = harminv_get_omega(hd, fsort[i]) / (2*pi*dt);
      freq_re[i] = abs(real(freq));
      freq_im[i] = imag(freq);
      amps[i] = harminv_get_amplitude(hd, fsort[i]);
      if (errors)
        errors[i] = harminv_get_freq_error(hd, fsort[i]);
    }

    delete[] fsort;
    harminv_data_destroy(hd);
    return nf;
#endif
  }

static double meep::dot	(	int	n,
		const double *	x,
		const double *	y
	)			[static]

{
  double sum = 0;
  for (int i = 0; i < n; ++i) sum += x[i] * y[i];
  return sum_to_all(sum);
}

static complex<double> meep::dot3_max_integrand	(	const complex< double > *	Fields,
		const vec &	loc,
		void *	data_
	)			[static]

{
  (void) loc; (void) data_; // unused;
  return (real(conj(Fields[0]) * Fields[3]) +
          real(conj(Fields[1]) * Fields[4]) +
          real(conj(Fields[2]) * Fields[5]));
}

static complex<double> meep::dot_fx_integrand	(	const complex< double > *	Fields,
		const vec &	loc,
		void *	data_
	)			[static]

                                                                     {
  fx_func fx = (fx_func) data_;
  return (real(conj(Fields[0]) * Fields[1]) * fx(loc));
}

static complex<double> meep::dot_integrand	(	const complex< double > *	Fields,
		const vec &	loc,
		void *	data_
	)			[static]

{
  (void) loc; (void) data_; // unused;
  return real(conj(Fields[0]) * Fields[1]);
}

Geometric_Volume meep::empty_volume

(

ndim

dim

)

                                        {
  Geometric_Volume out(dim);
  LOOP_OVER_DIRECTIONS(dim,d) {
    out.set_direction_max(d,0.0);
    out.set_direction_min(d,0.0);
  }
  return out;
}

void meep::end_critical_section

(

int

tag

)

{
#ifdef HAVE_MPI
     int process_rank, num_procs;
     MPI_Comm_rank(MPI_COMM_WORLD, &process_rank);
     MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
     if (process_rank != num_procs - 1) { /* send a message to next process */
          MPI_Send(&tag, 1, MPI_INT, process_rank + 1, tag, 
                   MPI_COMM_WORLD);
     }
#else
     UNUSED(tag);
#endif
}

static double meep::energy_fun	(	const complex< double > *	Fields,
		const vec &	loc,
		void *	data_
	)			[static]

{
     (void) loc; // unused
     int nFields = *((int *) data_) / 2;
     double sum = 0;
     for (int k = 0; k < nFields; ++k)
       sum += real(conj(Fields[2*k]) * Fields[2*k+1]);
     return sum * 0.5;
}

double meep::expi	(	int	cmp,
		double	x
	)			[inline]

                                        {
         return (cmp) ? cos(x) : sin(x);
  }

static void meep::fieldop	(	const double *	xr,
		double *	yr,
		void *	data_
	)			[static]

{
  const complex<double> *x = reinterpret_cast<const complex<double>*>(xr);
  complex<double> *y = reinterpret_cast<complex<double>*>(yr);
  fieldop_data *data = (fieldop_data *) data_;
  array_to_Fields(x, *data->f);
  data->f->step();
  Fields_to_array(*data->f, y);
  int n = data->n;
  double dt_inv = 1.0 / data->f->m_dt;
  complex<double> iomega = data->iomega;
  for (int i = 0; i < n; ++i) y[i] = (y[i] - x[i]) * dt_inv + iomega * x[i];
  data->iters++;
}

static void meep::Fields_to_array	(	const Fields &	f,
		complex< double > *	x
	)			[static]

{
  int ix = 0;
  for (int i=0;i<f.m_num_chunks;i++)
    if (f.m_chunks[i]->is_mine())
      FOR_COMPONENTS(c)
        if (f.m_chunks[i]->m_f[c][0] && (is_D(c) || is_magnetic(c))) {
          double *fs[3][2];
          // TODO: this is somewhat wasteful, because at least
          // one of f/m_f_p_pml/m_f_m_pml is redundant.
          DOCMP2 { fs[0][cmp] = f.m_chunks[i]->m_f[c][cmp];
                   fs[1][cmp] = f.m_chunks[i]->m_f_p_pml[c][cmp];
                   fs[2][cmp] = f.m_chunks[i]->m_f_m_pml[c][cmp]; }
          for (int k = 0; k < 3; ++k)
            if (fs[k][0]) {
              LOOP_OVER_VOL_OWNED(f.m_chunks[i]->m_vol, c, idx)
                x[ix++] = complex<double>(fs[k][0][idx], fs[k][1][idx]);
            }
        }
}

signed_direction meep::flip

(

signed_direction

d

)

[inline]

                                                 {
  signed_direction d2 = d;
  d2.flipped = !d.flipped;
  return d2;
}

static void meep::h5_findsize_chunkloop	(	Fields_Chunk *	fc,
		int	ichnk,
		component	cgrid,
		ivec	is,
		ivec	ie,
		vec	s0,
		vec	s1,
		vec	e0,
		vec	e1,
		double	dV0,
		double	dV1,
		ivec	shift,
		complex< double >	shift_phase,
		const Symmetry &	S,
		int	sn,
		void *	data_
	)			[static]

{
  h5_output_data *data = (h5_output_data *) data_;
  ivec isS = S.transform(is, sn) + shift;
  ivec ieS = S.transform(ie, sn) + shift;
  data->min_corner = min(data->min_corner, min(isS, ieS));
  data->max_corner = max(data->max_corner, max(isS, ieS));
  data->num_chunks++;
  int bufsz = 1;
  LOOP_OVER_DIRECTIONS(fc->m_vol.m_dim, d)
    bufsz *= (ie.in_direction(d) - is.in_direction(d)) / 2 + 1;
  data->bufsz = max(data->bufsz, bufsz);
}

static void meep::h5_output_chunkloop	(	Fields_Chunk *	fc,
		int	ichnk,
		component	cgrid,
		ivec	is,
		ivec	ie,
		vec	s0,
		vec	s1,
		vec	e0,
		vec	e1,
		double	dV0,
		double	dV1,
		ivec	shift,
		complex< double >	shift_phase,
		const Symmetry &	S,
		int	sn,
		void *	data_
	)			[static]

{
  h5_output_data *data = (h5_output_data *) data_;

  //-----------------------------------------------------------------------//
  // Find output chunk dimensions and strides, etc.

  int start[3]={0,0,0}, count[3]={1,1,1};
  int offset[3]={0,0,0}, stride[3]={1,1,1};

  ivec isS = S.transform(is, sn) + shift;
  ivec ieS = S.transform(ie, sn) + shift;
  
  // figure out what yucky_directions (in LOOP_OVER_IVECS)
  // correspond to what directions in the transformed vectors (in output).
  ivec permute(zero_ivec(fc->m_vol.m_dim));
  for (int i = 0; i < 3; ++i) 
    permute.set_direction(fc->m_vol.yucky_direction(i), i);
  permute = S.transform_unshifted(permute, sn);
  LOOP_OVER_DIRECTIONS(permute.dim, d)
    permute.set_direction(d, abs(permute.in_direction(d)));
  
  // compute the size of the chunk to output, and its strides etc.
  for (int i = 0; i < data->rank; ++i) {
    direction d = data->ds[i];
    int isd = isS.in_direction(d), ied = ieS.in_direction(d);
    start[i] = (min(isd, ied) - data->min_corner.in_direction(d)) / 2;
    count[i] = abs(ied - isd) / 2 + 1;
    int j = permute.in_direction(d);
    if (ied < isd) offset[permute.in_direction(d)] = count[i] - 1;
  }
  for (int i = 0; i < data->rank; ++i) {
    direction d = data->ds[i];
    int j = permute.in_direction(d);
    for (int k = i + 1; k < data->rank; ++k) stride[j] *= count[k];
    offset[j] *= stride[j];
    if (offset[j]) stride[j] *= -1;
  }
  
  //-----------------------------------------------------------------------//
  // Compute the function to output, exactly as in Fields::integrate,
  // except that here we store its values in a buffer instead of integrating.

  int *off = data->offsets;
  component *cS = data->cS;
  complex<double> *Fields = data->Fields, *ph = data->ph;
  complex<long double> sum = 0.0;
  double maxabs = 0;
  const component *iecs = data->inveps_cs;
  const direction *ieds = data->inveps_ds;
  int *ieos = data->inveps_offsets;

  for (int i = 0; i < data->num_Fields; ++i) {
    cS[i] = S.transform(data->components[i], -sn);
    if (cS[i] == Dielectric)
      ph[i] = 1.0;
    else {
      fc->m_vol.yee2diel_offsets(cS[i], off[2*i], off[2*i+1]);
      ph[i] = shift_phase * S.phase_shift(cS[i], sn);
    }
  }
  for (int k = 0; k < data->ninveps; ++k)
    fc->m_vol.yee2diel_offsets(iecs[k], ieos[2*k], ieos[2*k+1]);

  vec rshift(shift * (0.5*fc->m_vol.m_inva));
  LOOP_OVER_IVECS(fc->m_vol, is, ie, idx) {
    IVEC_LOOP_LOC(fc->m_vol, loc);
    loc = S.transform(loc, sn) + rshift;

    for (int i = 0; i < data->num_Fields; ++i) {
      if (cS[i] == Dielectric) {
        double tr = 0.0;
        for (int k = 0; k < data->ninveps; ++k)
          tr += (fc->m_structChunk->m_inveps[iecs[k]][ieds[k]][idx]
                 + fc->m_structChunk->m_inveps[iecs[k]][ieds[k]][idx+ieos[2*k]]
                 + fc->m_structChunk->m_inveps[iecs[k]][ieds[k]][idx+ieos[1+2*k]]
                 + fc->m_structChunk->m_inveps[iecs[k]][ieds[k]][idx+ieos[2*k]+ieos[1+2*k]]);
        Fields[i] = (4 * data->ninveps) / tr;
      }
      else {
        double f[2];
        for (int k = 0; k < 2; ++k)
          if (fc->m_f[cS[i]][k])
            f[k] = 0.25 * (fc->m_f[cS[i]][k][idx]
                           + fc->m_f[cS[i]][k][idx+off[2*i]]
                           + fc->m_f[cS[i]][k][idx+off[2*i+1]]
                           + fc->m_f[cS[i]][k][idx+off[2*i]+off[2*i+1]]);
          else
            f[k] = 0;
        Fields[i] = complex<double>(f[0], f[1]) * ph[i];
      }
    }

    complex<double> fun = data->fun(Fields, loc, data->fun_data_);
    int idx2 = ((((offset[0] + offset[1] + offset[2])
                  + loop_i1 * stride[0]) 
                 + loop_i2 * stride[1]) + loop_i3 * stride[2]);
    data->buf[idx2] = data->reim ? imag(fun) : real(fun);
  }

  //-----------------------------------------------------------------------//

  data->file->write_chunk(data->rank, start, count, data->buf);
}

static void meep::handle_control_c

(

int

i

)

[static]

                                    {
  (void) i; // unused: should equal SIGINT
  interrupt++;
  if (interrupt >= kill_time) {
    abort("interrupted");
  } else if (interrupt + 1 == kill_time) {
    printf("Be patient... hit ctrl-C one more time to kill me.\n");
  } else {
    printf("Be patient... hit ctrl-C %d more times to kill me.\n", kill_time - interrupt);
  }
}

bool meep::has_direction	(	ndim	dim,
		direction	d
	)			[inline]

                                                 {
  LOOP_OVER_DIRECTIONS(dim, dd) if (dd == d) return true;
  return false;
}

static int meep::iabs

(

int

i

)

[inline, static]

{ return (i < 0 ? -i : i); }

Symmetry meep::identity

(

)

                    {
  return Symmetry();
}

static void meep::integrate_chunkloop	(	Fields_Chunk *	fc,
		int	ichunk,
		component	cgrid,
		ivec	is,
		ivec	ie,
		vec	s0,
		vec	s1,
		vec	e0,
		vec	e1,
		double	dV0,
		double	dV1,
		ivec	shift,
		complex< double >	shift_phase,
		const Symmetry &	S,
		int	sn,
		void *	data_
	)			[static]

{
  integrate_data *data = (integrate_data *) data_;
  int *off = data->offsets;
  component *cS = data->cS;
  complex<double> *Fields = data->Fields, *ph = data->ph;
  complex<long double> sum = 0.0;
  double maxabs = 0;
  const component *iecs = data->inveps_cs;
  const direction *ieds = data->inveps_ds;
  int *ieos = data->inveps_offsets;

  for (int i = 0; i < data->num_Fields; ++i) {
    cS[i] = S.transform(data->components[i], -sn);
    if (cS[i] == Dielectric)
      ph[i] = 1.0;
    else {
      if (cgrid == Dielectric)
                fc->m_vol.yee2diel_offsets(cS[i], off[2*i], off[2*i+1]);
      ph[i] = shift_phase * S.phase_shift(cS[i], sn);
    }
  }
  for (int k = 0; k < data->ninveps; ++k)
    fc->m_vol.yee2diel_offsets(iecs[k], ieos[2*k], ieos[2*k+1]);

  vec rshift(shift * (0.5*fc->m_vol.m_inva));
  LOOP_OVER_IVECS(fc->m_vol, is, ie, idx) {
    IVEC_LOOP_LOC(fc->m_vol, loc);
    loc = S.transform(loc, sn) + rshift;

    for (int i = 0; i < data->num_Fields; ++i) {
      if (cS[i] == Dielectric) {
        double tr = 0.0;
        for (int k = 0; k < data->ninveps; ++k)
          tr += (fc->m_structChunk->m_inveps[iecs[k]][ieds[k]][idx]
                 + fc->m_structChunk->m_inveps[iecs[k]][ieds[k]][idx+ieos[2*k]]
                 + fc->m_structChunk->m_inveps[iecs[k]][ieds[k]][idx+ieos[1+2*k]]
                 + fc->m_structChunk->m_inveps[iecs[k]][ieds[k]][idx+ieos[2*k]+ieos[1+2*k]]);
        Fields[i] = (4 * data->ninveps) / tr;
      }
      else {
        double f[2];
        for (int k = 0; k < 2; ++k)
          if (fc->m_f[cS[i]][k])
            f[k] = 0.25 * (fc->m_f[cS[i]][k][idx]
                           + fc->m_f[cS[i]][k][idx+off[2*i]]
                           + fc->m_f[cS[i]][k][idx+off[2*i+1]]
                           + fc->m_f[cS[i]][k][idx+off[2*i]+off[2*i+1]]);
          else
            f[k] = 0;
        Fields[i] = complex<double>(f[0], f[1]) * ph[i];
      }
    }

    complex<double> integrand = 
      data->integrand(Fields, loc, data->integrand_data_);
    maxabs = max(maxabs, abs(integrand));
    sum += integrand * IVEC_LOOP_WEIGHT(s0, s1, e0, e1, dV0 + dV1 * loop_i2);
  }

  data->maxabs = max(data->maxabs, maxabs);
  data->sum += sum;
}

bool meep::is_D

(

component

c

)

[inline]

{ return c >= Dx && c < Dielectric; }

bool meep::is_derived

(

int

c

)

[inline]

{ return c >= Sx; }

bool meep::is_electric

(

component

c

)

[inline]

{ return c < Hx; }

bool meep::is_energydensity

(

derived_component

c

)

[inline]

{ return c>=EnergyDensity; }

static bool meep::is_like	(	ndim	d,
		component	c1,
		component	c2
	)			[static]

                                                          {
    if (d != D2) return true;
    return !(is_tm(c1) ^ is_tm(c2));
  }

bool meep::is_magnetic

(

component

c

)

[inline]

{ return c >= Hx && c < Dx; }

static bool meep::is_ok_dir

(

const char *

dirname

)

[static]

                                           {
  DIR *dir;
  bool direxists;
  if (am_master()) {
    direxists = (dir = opendir(dirname)) != NULL;
    if (direxists) closedir(dir);
    else mkdir(dirname, 00777);
  }
  direxists = broadcast(0, direxists);
  return !direxists;
}

bool meep::is_poynting

(

derived_component

c

)

[inline]

{ return c < EnergyDensity; }

static bool meep::is_tm

(

component

c

)

[static]

                                 {
    switch (c) {
    case Hx: case Hy: case Ez: case Dz: return true;
    default: return false;
    }
    return false;
  }

static double meep::it	(	int	cmp,
		double *	f[2],
		int	ind
	)			[inline, static]

                                                          {
  return (f[1-cmp]) ? (1-2*cmp)*f[1-cmp][ind] : 0;
}

ivec meep::iveccyl	(	int	xx,
		int	yy
	)			[inline]

                                    {
  ivec v(Dcyl); v.t[R] = rr; v.t[Z] = zz; return v;
}

static double meep::J	(	int	m,
		double	kr
	)			[static]

                                  {
#ifdef HAVE_LIBGSL
  return gsl_sf_bessel_Jn(m, kr);
#else
  abort("not compiled with GSL, required for Bessel functions");
  return 0;
#endif
}

static complex<double> meep::JJ

(

const vec &

v

)

[static]

                                        {
  return polar(J(m_for_J, ktrans*v.r()),kax*v.r());
}

static double meep::Jmax	(	int	m,
		int	n
	)			[static]

                                 {
  double rlow, rhigh = Jroot(m,n), rtry;
  if (n == 0) rlow = 0;
  else rlow = Jroot(m, n-1);
  double jplow = Jprime(m,rlow), jptry;
  do {
    rtry = rlow + (rhigh - rlow)*0.5;
    jptry = Jprime(m,rtry);
    if (jplow*jptry < 0) rhigh = rtry;
    else rlow = rtry;
  } while (rhigh - rlow > rhigh*1e-15);
  return rtry;
}

static complex<double> meep::JP

(

const vec &

v

)

[static]

                                        {
  return polar(Jprime(m_for_J, ktrans*v.r()),kax*v.r());
}

static double meep::Jprime	(	int	m,
		double	kr
	)			[static]

                                       { 
  if (m) return 0.5*(J(m-1,kr)-J(m+1,kr));
  else return -J(1,kr);
}

static double meep::Jroot	(	int	m,
		int	n
	)			[static]

                                  {
#ifdef HAVE_LIBGSL
  return gsl_sf_bessel_zero_Jnu(m, n+1);
#else
  abort("not compiled with GSL, required for Bessel functions");
  return 0;
#endif
}

void meep::load_dft_hdf5	(	Dft_Chunk *	Dft_Chunks,
		component	c,
		H5file *	file,
		const char *	dprefix
	)

                                                              {
        int istart;
        int n = Dft_Chunks_Ntotal(Dft_Chunks, &istart);

        char dataname[1024];
        snprintf(dataname, 1024, "%s%s" "%s_dft", 
                         dprefix ? dprefix : "", dprefix && dprefix[0] ? "_" : "",
                         component_name(c));
        int file_rank, file_dims;
        file->read_size(dataname, &file_rank, &file_dims, 1);
        if (file_rank != 1 || file_dims != n)
          abort("incorrect dataset size (%d vs. %d) in load_dft_hdf5 %s:%s", file_dims, n, file->file_name(), dataname);
  
        for (Dft_Chunk *cur = Dft_Chunks; cur; cur = cur->m_next_in_dft) {
          int Nchunk = cur->m_N * cur->m_Nomega * 2;
          file->read_chunk(1, &istart, &Nchunk, (double *) cur->m_dft);
          istart += Nchunk;
        }
  }

const char * meep::make_output_directory	(	const char *	exename,
		const char *	jobname = NULL
	)

Utility function to parse the executable name use it to come up with a directory name, avoiding overwriting any existing directory, unless the source file hasn't changed.

                                                                            {
  const int buflen = 300;
  char basename[buflen];
  const char * const evil_suffs[] = { ".dac", ".cpp", ".cc", ".cxx", ".C" };
  char stripped_name[buflen];
  const char *bnp = exename; // stripped_name holds the actual name of the executable (dirs removed).
  const char *t;
  for (t=exename;*t;t++) {
    if (*t == '/') bnp = t+1;
  }
  
  snprintf(stripped_name, buflen, "%s", bnp);
  for (int i = 0; i < (int)(sizeof(evil_suffs) / sizeof(evil_suffs[0])); ++i) {
    int sufflen = strlen(evil_suffs[i]);
    if (strcmp(stripped_name + strlen(stripped_name) - sufflen, 
               evil_suffs[i]) == 0
        && strlen(stripped_name) > size_t(sufflen)) {
      stripped_name[strlen(stripped_name) - sufflen] = (char)0;
      break;
    }
  }

  char sourcename[buflen]; // Holds the "example.cpp" filename.
  snprintf(sourcename, buflen, "%s.cpp", stripped_name);

  if (jobname != NULL) {
    snprintf(basename, buflen, "%s", jobname);
  }
  else {
    snprintf(basename, buflen, "%s", stripped_name);
  }

  static char outdirname[buflen];
  snprintf(outdirname, buflen, "%s-out", basename);
  {
    int i = 0;
    while (!is_ok_dir(outdirname)) {
      if (!quiet)
        master_printf("Output directory %s already exists!\n", outdirname);
      snprintf(outdirname, buflen, "%s-out-%d", basename, i++);
    }
  }
  char outsrcname[buflen];
  snprintf(outsrcname, buflen, "%s/%s", outdirname, sourcename);
  cp(sourcename, outsrcname);

  return outdirname;
}

void meep::master_fclose

(

FILE *

f

)

                            {
  if (am_master()) fclose(f);
}

void void void FILE * meep::master_fopen	(	const char *	name,
		const char *	mode
	)

                                                       {
  FILE *f = am_master() ? fopen(name, mode) : 0;

  /* other processes need to know if fopen returned zero, in order
     to abort if fopen failed.  If fopen was successfully, just return
     a random non-zero pointer (which is never used except to compare to zero)
     on non-master processes */
  if (broadcast(0, bool(f != 0)) && !am_master())
    f = (FILE *) name;
  return f;
}

void void void meep::master_fprintf	(	FILE *	f,
		const char *	fmt,
			...
	)

                                                   {
  va_list ap;
  va_start(ap, fmt);
  if (am_master()) { vfprintf(f, fmt, ap); fflush(f); }
  va_end(ap);
}

void meep::master_printf	(	const char *	fmt,
			...
	)

                                         {
  va_list ap;
  va_start(ap, fmt);
  if (am_master()) { vprintf(fmt, ap); fflush(stdout); }
  va_end(ap);
}

ivec meep::max	(	const ivec &	ivec1,
		const ivec &	ivec2
	)

                                               {
  ivec m(ivec1.dim);
  LOOP_OVER_DIRECTIONS(ivec1.dim, d)
    m.set_direction(d, max(ivec1.in_direction(d), ivec2.in_direction(d)));
  return m;
}

vec meep::max	(	const vec &	vec1,
		const vec &	vec2
	)

                                          {
  vec m(vec1.dim);
  LOOP_OVER_DIRECTIONS(vec1.dim, d)
    m.set_direction(d, max(vec1.in_direction(d), vec2.in_direction(d)));
  return m;
}

int meep::max	(	int	a,
		int	b
	)			[inline]

{ return (a > b) ? a : b; }

double meep::max	(	double	a,
		double	b
	)			[inline]

{ return (a > b) ? a : b; }

ivec meep::max_to_all

(

const ivec &

v

)

                               {
  int in[5], out[5];
  for (int i=0; i<5; ++i) in[i] = out[i] = v.in_direction(direction(i));
#ifdef HAVE_MPI
  MPI_Allreduce(&in,&out,5,MPI_INT,MPI_MAX,MPI_COMM_WORLD);
#endif
  ivec vout(v.dim);
  for (int i=0; i<5; ++i) vout.set_direction(direction(i), out[i]);
  return vout;
}

int meep::max_to_all

(

int

in

)

                       {
  int out = in;
#ifdef HAVE_MPI
  MPI_Allreduce(&in,&out,1,MPI_INT,MPI_MAX,MPI_COMM_WORLD);
#endif
  return out;
}

double meep::max_to_all

(

double

in

)

                             {
  double out = in;
#ifdef HAVE_MPI
  MPI_Allreduce(&in,&out,1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD);
#endif
  return out;
}

double meep::max_to_master

(

double

in

)

                                {
  double out = in;
#ifdef HAVE_MPI
  MPI_Reduce(&in,&out,1,MPI_DOUBLE,MPI_MAX,0,MPI_COMM_WORLD);
#endif
  return out;
}

ivec meep::min	(	const ivec &	ivec1,
		const ivec &	ivec2
	)

                                               {
  ivec m(ivec1.dim);
  LOOP_OVER_DIRECTIONS(ivec1.dim, d)
    m.set_direction(d, min(ivec1.in_direction(d), ivec2.in_direction(d)));
  return m;
}

vec meep::min	(	const vec &	vec1,
		const vec &	vec2
	)

                                          {
  vec m(vec1.dim);
  LOOP_OVER_DIRECTIONS(vec1.dim, d)
    m.set_direction(d, min(vec1.in_direction(d), vec2.in_direction(d)));
  return m;
}

int meep::min	(	int	a,
		int	b
	)			[inline]

{ return (a < b) ? a : b; }

double meep::min	(	double	a,
		double	b
	)			[inline]

{ return (a < b) ? a : b; }

Symmetry meep::mirror	(	direction	axis,
		const Volume &	v
	)

                                                 {
  Symmetry s = identity();
  s.m_g = 2;
  s.m_S[axis].flipped = true;
  s.m_symmetry_point = v.center();
  s.m_i_symmetry_point = v.icenter();
  return s;
}

int meep::my_rank

(

)

              {
#ifdef HAVE_MPI
  int rank;
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  return rank;
#else
  return 0;
#endif
}

static int meep::my_round

(

double

x

)

[inline, static]

                                     {
  return int(floor(fabs(x) + 0.5) * (x < 0 ? -1 : 1));
}

static double meep::norm2	(	int	n,
		const double *	x
	)			[static]

{ return sqrt(dot(n, x, x)); }

int meep::number_of_directions

(

ndim

dim

)

[inline]

                                          {
  return (int) (dim + 1 - 2 * (dim == Dcyl));
}

complex<double> meep::one

(

const vec &

v

)

Warning:

Badly named?

{(void) v; return 1.0;}

ivec meep::one_ivec

(

ndim

di

)

[inline]

                              {
  ivec v; v.dim = di; LOOP_OVER_DIRECTIONS(di, d) v.set_direction(d, 1);
  return v;
}

bool meep::or_to_all

(

bool

in

)

                        {
  int out = in;
#ifdef HAVE_MPI
  MPI_Allreduce(&in,&out,1,MPI_INT,MPI_LOR,MPI_COMM_WORLD);
#endif
  return (bool) out;
}

int meep::partial_sum_to_all

(

int

in

)

                               {
  int out = in;
#ifdef HAVE_MPI
  MPI_Scan(&in,&out,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
#endif
  return out;
}

Boundary_Region meep::pml

(

double

thickness

)

                                      {
  Boundary_Region r;
  for (int id = 0; id < 5; ++id)
    r = r + pml(thickness, (direction) id);
  return r;
}

Boundary_Region meep::pml	(	double	thickness,
		direction	d
	)

                                                   {
  return (pml(thickness, d, Low) + pml(thickness, d, High));
}

Boundary_Region meep::pml	(	double	thickness,
		direction	d,
		boundary_side	side
	)

                                                                       {
  return Boundary_Region(Boundary_Region::PML, thickness, 1.0, d, side, NULL);
}

static double meep::poynting_fun	(	const complex< double > *	Fields,
		const vec &	loc,
		void *	data_
	)			[static]

{
     (void) loc; // unused
     (void) data_; // unused
     return (real(conj(Fields[0]) * Fields[1])
             - real(conj(Fields[2])*Fields[3]));
}

static void meep::pt	(	double	ts[],
		time_sink	s
	)			[static]

                                         {
  if (ts[s]) master_printf("    %18s: %g s\n", ts2n(s), ts[s]);
}

Symmetry meep::r_to_minus_r_Symmetry

(

double

m

)

Symmetry meep::r_to_minus_r_symmetry

(

double

m

)

                                         {
  Symmetry s = identity();
  s.m_g = 2;
  s.m_S[R].flipped = true;
  s.m_S[P].flipped = true;
  s.m_symmetry_point = zero_vec(Dcyl);
  s.m_i_symmetry_point = zero_ivec(Dcyl);
  if (m == int(m)) // phase is purely real (+/- 1) when m an integer
    s.m_ph = (int(m) & 1) ? -1.0 : 1.0;
  else
    s.m_ph = polar(1.0, m * pi); // general case
  return s;
}

Symmetry meep::r_to_minus_r_Symmetry

(

int

m

)

static complex<double> meep::return_the_field	(	const complex< double > *	Fields,
		const vec &	loc,
		void *	integrand_data_
	)			[static]

{
  (void) integrand_data_; (void) loc; // unused
  return Fields[0];
}

static complex<double> meep::rfun_wrap	(	const complex< double > *	Fields,
		const vec &	loc,
		void *	data_
	)			[static]

                                                              {
  rfun_wrap_data *data = (rfun_wrap_data *) data_;
  return data->integrand(Fields, loc, data->integrand_data);
}

static complex<double> meep::rintegrand_fun	(	const complex< double > *	Fields,
		const vec &	loc,
		void *	data_
	)			[static]

{
  rintegrand_data *data = (rintegrand_data *) data_;
  return data->fun(Fields, loc, data->fun_data_);
}

int meep::rmin_bulk

(

int

m

)

[inline]

                            {
  int r = 1 + small_r_metal(m);
  if (r < 1) r = 1;
  return r;
}

Symmetry meep::rotate2	(	direction	axis,
		const Volume &	v
	)

                                                  {
  Symmetry s = identity();
  if (axis > 2) abort("Can only rotate2 in 2D or 3D.\n");
  s.m_g = 2;
  s.m_S[(axis+1)%3].flipped = true;
  s.m_S[(axis+2)%3].flipped = true;
  s.m_symmetry_point = v.center();
  s.m_i_symmetry_point = v.icenter();
  return s;
}

Symmetry meep::rotate4	(	direction	axis,
		const Volume &	v
	)

                                                  {
  Symmetry s = identity();
  if (axis > 2) abort("Can only rotate4 in 2D or 3D.\n");
  s.m_g = 4;
  FOR_DIRECTIONS(d) {
    s.m_S[d].d = d;
    s.m_S[d].flipped = false;
  }
  s.m_S[(axis+1)%3].d = (direction)((axis+2)%3);
  s.m_S[(axis+1)%3].flipped = true;
  s.m_S[(axis+2)%3].d = (direction)((axis+1)%3);
  s.m_symmetry_point = v.center();
  s.m_i_symmetry_point = v.icenter();
  return s;
}

int meep::rstart_0	(	const Volume &	v,
		double	m
	)			[inline]

                                               {
  return (int) max(0.0, m - (int)(v.origin_r()*v.m_a+0.5) - 1.0);
}

int meep::rstart_1	(	const Volume &	v,
		double	m
	)			[inline]

                                               {
  return (int) max(1.0, m - (int)(v.origin_r()*v.m_a+0.5));
}

void meep::save_dft_hdf5	(	Dft_Chunk *	Dft_Chunks,
		component	c,
		H5file *	file,
		const char *	dprefix
	)

                                                              {
        int istart;
        int n = Dft_Chunks_Ntotal(Dft_Chunks, &istart);

        char dataname[1024];
        snprintf(dataname, 1024, "%s%s" "%s_dft", 
                         dprefix ? dprefix : "", dprefix && dprefix[0] ? "_" : "", 
                         component_name(c));
        file->create_data(dataname, 1, &n);

        for (Dft_Chunk *cur = Dft_Chunks; cur; cur = cur->m_next_in_dft) {
          int Nchunk = cur->m_N * cur->m_Nomega * 2;
          file->write_chunk(1, &istart, &Nchunk, (double *) cur->m_dft);
          istart += Nchunk;
        }
        file->done_writing_chunks();
  }

void meep::send	(	int	from,
		int	to,
		double *	data,
		int	size
	)

                                                    {
#ifdef HAVE_MPI
  if (from == to) return;
  if (size == 0) return;
  const int me = my_rank();
  if (from == me) MPI_Send(data, size, MPI_DOUBLE, to, 1, MPI_COMM_WORLD);
  MPI_Status stat;
  if (to == me) MPI_Recv(data, size, MPI_DOUBLE, from, 1, MPI_COMM_WORLD, &stat);
#else
  UNUSED(from);
  UNUSED(to);
  UNUSED(data);
  UNUSED(size);
#endif
}

int meep::small_r_metal

(

int

m

)

[inline]

                                {
  return m-1;
}

static vec meep::sphere_pt	(	const vec &	cent,
		double	R,
		int	n,
		double &	weight
	)			[static]

00095          :
00096     while ((fabs(meps-old_meps) > tol*old_meps) && (fabs(minveps-old_minveps) > tol*old_minveps)) {
00097       old_meps=meps; old_minveps=minveps;
00098       meps = minveps = 0;
00099       for (int j=0; j < ms; j++)
00100         for (int i=0; i < ms; i++) {
00101           double ep = eps(gv.get_min_corner() + vec(i*d.x()/ms, j*d.y()/ms));
00102           meps += ep; minveps += 1/ep;
00103         }
00104       meps /= ms*ms;
00105       minveps /= ms*ms;
00106       ms *= 2; 
00107       if (maxeval && (iter += ms*ms) >= maxeval) return;
00108     }
00109     break;
00110   case Dcyl:
00111     while ((fabs(meps-old_meps) > tol*old_meps) && (fabs(minveps-old_minveps) > tol*old_minveps)) {
00112       old_meps=meps; old_minveps=minveps;
00113       meps = minveps = 0;
00114       double sumvol = 0;
00115       for (int j=0; j < ms; j++)
00116         for (int i=0; i < ms; i++) {
00117           double r = gv.get_min_corner().r() + i*d.r()/ms;
00118           double ep = eps(gv.get_min_corner() + veccyl(i*d.r()/ms, j*d.z()/ms));
00119           sumvol += r;
00120           meps += ep * r; minveps += r/ep;
00121         }
00122       meps /= sumvol;
00123       minveps /= sumvol;
      ms *= 2; 

bool meep::Src_Times_equal	(	const Src_Time &	t1,
		const Src_Time &	t2
	)

this function is necessary to make equality commutative .

.. ugh warning Grok Me!

  {
        return t1.is_equal(t2) && t2.is_equal(t1);
  }

void meep::Src_Vol_chunkloop	(	Fields_Chunk *	fc,
		int	ichunk,
		component	c,
		ivec	is,
		ivec	ie,
		vec	s0,
		vec	s1,
		vec	e0,
		vec	e1,
		double	dV0,
		double	dV1,
		ivec	shift,
		complex< double >	shift_phase,
		const Symmetry &	S,
		int	sn,
		void *	data_
	)

Adding source volumes can be treated as a kind of "integration" problem, since we need to loop over all the chunks that intersect the source volume, with appropriate interpolation weights at the boundaries so that the integral of the current is fixed regardless of resolution.

Unlike most uses of Fields::loop_in_chunks, however, we set use_Symmetry=false: we only find the intersection of the volume with the untransformed chunks (since the transformed versions are implicit).

I put this in the header becuase it's declared static, and called in fields, and this was a shortcut. This, and other static functions, need to check if it really needs to be static. If so, then do it properly...

Parameters:

	fc	pointer to Fields_Chunk.
	ichunk	dunno.
	c	A component, why this needs to be spec'd I don't know.
	is	dunno.
	ie	dunno.
	s0	dunno.
	s1	dunno.
	e0	dunno.
	e1	dunno.
	dV0	dunno.
	dV1	dunno.
	shift	dunno.
	shift_phase	dunno.
	S	reference to a Symmetry.
	sn	dunno.
	data	Okay, declared with a void* type, but there is a

Warning:

I removed the static qualifier here becuase I saw no need for it. need to check performance, also compare parallel code, make sure nothing is broken.

  {
        Src_Vol_Chunkloop_Data* data = (Src_Vol_Chunkloop_Data *) data_;
        (void) S; (void) sn; // these should be the identity
        (void) dV0; (void) dV1; // volume weighting is included in data->amp
        (void) ichunk;

        int npts = 1;
        LOOP_OVER_DIRECTIONS(is.dim, d)
          npts *= (ie.in_direction(d) - is.in_direction(d)) / 2 + 1;
        int *index_array = new int[npts];
        complex<double> *amps_array = new complex<double>[npts];

        complex<double> amp = data->amp * conj(shift_phase);

        direction cd = component_direction(c);
        double inva = fc->m_vol.m_inva;
        int idx_vol = 0;
        LOOP_OVER_IVECS(fc->m_vol, is, ie, idx) {
          IVEC_LOOP_LOC(fc->m_vol, loc);
          loc += shift * (0.5*inva) - data->center;

          amps_array[idx_vol] = IVEC_LOOP_WEIGHT(s0,s1,e0,e1,1) * amp * data->A(loc);

          /* for "D" sources, multiply by epsilon.  FIXME: this is not quite
                 right because it doesn't handle non-diagonal inveps! */
          if (is_D(c) && fc->m_structChunk->m_inveps[c][cd]) 
                amps_array[idx_vol] /= fc->m_structChunk->m_inveps[c][cd][idx];

          index_array[idx_vol++] = idx;
        }

        if (idx_vol != npts)
          abort("add_volume_source: computed wrong npts (%d vs. %d)", npts, idx_vol);

        Src_Vol *tmp = new Src_Vol(c, data->src, npts, index_array, amps_array);
        if (is_magnetic(c))
          fc->m_h_sources = tmp->add_to(fc->m_h_sources);
        else
          fc->m_e_sources = tmp->add_to(fc->m_e_sources);
  }

direction meep::start_at_direction

(

ndim

dim

)

[inline]

                                              {
  return (direction) (((dim == D1) || (dim == Dcyl)) ? 2 : 0);
}

direction meep::stop_at_direction

(

ndim

dim

)

[inline]

                                             {
  return (direction) (dim + 1 + 2 * (dim == D1));
}

static complex<double> meep::stupid_A

(

const vec &

p

)

[static]

                                              {
  return A_static(c_static, p);
}

static complex<double> meep::stupider_A

(

const vec &

p

)

[static]

                                                {
  return prefactor_static*f_static->get_field(c_static, p);
}

static void meep::stupidsort	(	ivec *	locs,
		double *	w,
		int	l
	)			[inline, static]

                                                            {
  while (l) {
    if (fabs(w[0]) < 2e-15) {
      w[0] = w[l-1];
      locs[0] = locs[l-1];
      w[l-1] = 0.0;
      locs[l-1] = 0;
    } else {
      w += 1;
      locs += 1;
    }
    l -= 1;
  }
}

static void meep::stupidsort	(	int *	ind,
		double *	w,
		int	l
	)			[inline, static]

                                                          {
  while (l) {
    if (fabs(w[0]) < 2e-15) {
      w[0] = w[l-1];
      ind[0] = ind[l-1];
      w[l-1] = 0.0;
      ind[l-1] = 0;
    } else {
      w += 1;
      ind += 1;
    }
    l -= 1;
  }
}

complex< long double > meep::sum_to_all

(

complex< long double >

in

)

                                                         {
  complex<long double> out = in;
#ifdef HAVE_MPI
  MPI_Allreduce(&in,&out,2,MPI_LONG_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
#endif
  return out;
}

complex< double > meep::sum_to_all

(

complex< double >

in

)

                                               {
  complex<double> out = in;
#ifdef HAVE_MPI
  MPI_Allreduce(&in,&out,2,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
#endif
  return out;
}

int meep::sum_to_all

(

int

in

)

                       {
  int out = in;
#ifdef HAVE_MPI
  MPI_Allreduce(&in,&out,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
#endif
  return out;
}

long double meep::sum_to_all

(

long double

in

)

                                       {
  long double out = in;
#ifdef HAVE_MPI
  MPI_Allreduce(&in,&out,1,MPI_LONG_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
#endif
  return out;
}

void meep::sum_to_all	(	const double *	in,
		double *	out,
		int	size
	)

                                                         {
  memcpy(out, in, sizeof(double) * size);
#ifdef HAVE_MPI
  MPI_Allreduce((void*) in, out, size, MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
#endif
}

double meep::sum_to_all

(

double

in

)

                             {
  double out = in;
#ifdef HAVE_MPI
  MPI_Allreduce(&in,&out,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
#endif
  return out;
}

double meep::sum_to_master

(

double

in

)

                                {
  double out = in;
#ifdef HAVE_MPI
  MPI_Reduce(&in,&out,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
#endif
  return out;
}

static void meep::thermo_chunkloop	(	Fields_Chunk *	fc,
		int	ichunk,
		component	cgrid,
		ivec	is,
		ivec	ie,
		vec	s0,
		vec	s1,
		vec	e0,
		vec	e1,
		double	dV0,
		double	dV1,
		ivec	shift,
		complex< double >	shift_phase,
		const Symmetry &	S,
		int	sn,
		void *	sum_
	)			[static]

                                         {
  long double *sum = (long double *) sum_;
  (void) shift; (void) shift_phase; (void) S; (void) sn; // unused
  if (fc->m_pol && fc->m_pol->m_energy[cgrid])
    LOOP_OVER_IVECS(fc->m_vol, is, ie, idx)
      *sum += IVEC_LOOP_WEIGHT(s0, s1, e0, e1, dV0 + dV1 * loop_i2)
          * fc->m_pol->m_energy[cgrid][idx];
}

void meep::trash_output_directory

(

const char *

dirname

)

                                                 {
  if (am_master()) mkdir(dirname, 00777);
}

static const char* meep::ts2n

(

time_sink

s

)

[static]

                                     {
  switch (s) {
  case e_stepping: return "time stepping";
  case e_connecting: return "connnecting chunks";
  case e_boundaries: return "copying borders";
  case e_mpiTime: return "communicating";
  case e_fieldOutput: return "outputting Fields";
  case e_fourierTransforming: return "Fourier transforming";
  case e_other: break;
  }
  return "everything else";
}

field_type meep::type

(

component

c

)

[inline]

                                    {
  if (is_electric(c)) return E_stuff;
  else if (is_magnetic(c)) return H_stuff;
  else if (is_D(c)) return D_stuff;
  abort("Invalid field in type.\n");
  return E_stuff; // This is never reached.
}

ivec meep::unit_ivec	(	ndim	di,
		direction	d
	)			[inline]

                                            {
  ivec v(zero_ivec(di));
  v.set_direction(d, 1);
  return v;
}

vec meep::unit_vec	(	ndim	di,
		direction	d
	)			[inline]

                                          {
  vec v(zero_vec(di));
  v.set_direction(d, 1.0);
  return v;
}

static ivec meep::vec2diel_ceil	(	const vec &	v,
		double	a,
		const ivec &	equal_shift
	)			[static]

                                                                           {
  ivec iv(v.dim);
  LOOP_OVER_DIRECTIONS(v.dim, d) {
    iv.set_direction(d, 1+2*int(ceil(v.in_direction(d)*a-.5)));
    if (iv.in_direction(d) == v.in_direction(d))
      iv.set_direction(d, iv.in_direction(d) + equal_shift.in_direction(d));
  }
  return iv;
}

static ivec meep::vec2diel_floor	(	const vec &	v,
		double	a,
		const ivec &	equal_shift
	)			[static]

                                                                            {
  ivec iv(v.dim);
  LOOP_OVER_DIRECTIONS(v.dim, d) {
    iv.set_direction(d, 1+2*int(floor(v.in_direction(d)*a-.5)));
    if (iv.in_direction(d) == v.in_direction(d))
      iv.set_direction(d, iv.in_direction(d) + equal_shift.in_direction(d));
  }
  return iv;
}

vec meep::veccyl	(	double	rr,
		double	zz
	)			[inline]

                                        {
  vec v(Dcyl); v.t[R] = rr; v.t[Z] = zz; return v;
}

Volume meep::vol1d	(	double	zsize,
		double	a
	)

                                     {
  return volone(zsize, a);
}

Volume meep::vol2d	(	double	xsize,
		double	ysize,
		double	a
	)

                                                   {
  return voltwo(xsize, ysize, a);
}

Volume meep::vol3d	(	double	xsize,
		double	ysize,
		double	zsize,
		double	a
	)

                                                                 {
  return Volume(D3, a,(xsize==0)?1:(int) (xsize*a + 0.5),
                      (ysize==0)?1:(int) (ysize*a + 0.5),
                      (zsize==0)?1:(int) (zsize*a + 0.5));
}

Volume meep::volcyl	(	double	rsize,
		double	zsize,
		double	a
	)

                                                    {
  if (zsize == 0.0) return Volume(Dcyl, a, (int) (rsize*a + 0.5), 0, 1);
  else return Volume(Dcyl, a, (int) (rsize*a + 0.5), 0, (int) (zsize*a + 0.5));
}

Volume meep::volone	(	double	zsize,
		double	a
	)

                                      {
  return Volume(D1, a, 0, 0, (int) (zsize*a + 0.5));
}

Volume meep::voltwo	(	double	xsize,
		double	ysize,
		double	a
	)

                                                    {
  return Volume(D2, a, (xsize==0)?1:(int) (xsize*a + 0.5),
                       (ysize==0)?1:(int) (ysize*a + 0.5),0);
}

double meep::wall_time

(

void

)

                       {
#ifdef HAVE_MPI
  return MPI_Wtime();
#elif HAVE_GETTIMEOFDAY
  struct timeval tv;
  gettimeofday(&tv, 0);
  return(tv.tv_sec + tv.tv_usec * 1e-6);
#else
  return (clock() * 1.0 / CLOCKS_PER_SECOND);
#endif
}

static void meep::xpay	(	int	n,
		double *	x,
		double	a,
		const double *	y
	)			[static]

                                                              {
  for (int m = 0; m < n; ++m) x[m] += a * y[m];
}

static direction meep::yucky_dir	(	ndim	dim,
		int	n
	)			[static]

                                            {
  if (dim == Dcyl)
    switch (n) {
    case 0: return P;
    case 1: return R;
    case 2: return Z;
    }
  else if (dim == D2)
    return (direction) ((n + 2) % 3); /* n = 0,1,2 gives Z, X, Y */
  return (direction) n ;
}

ivec meep::zero_ivec

(

ndim

di

)

[inline]

                               {
  ivec v; v.dim = di; LOOP_OVER_DIRECTIONS(di, d) v.set_direction(d, 0);
  return v;
}

vec meep::zero_vec

(

ndim

di

)

[inline]

                             {
  vec v(di); LOOP_OVER_DIRECTIONS(di, d) v.set_direction(d, 0.0);
  return v;
}

---

Variable Documentation

complex<double>(* meep::A_static)(component, const vec &) [static]

component meep::c_static [static]

const double meep::Cmax = 0.5

FILE* meep::debf = NULL [static]

Fields* meep::f_static [static]

const char meep::Grace_header[]

Initial value:

 "# Grace project file\
#\
@page size 792, 612\
@default symbol size 0.330000\
@g0 on\
@with g0\
"

const double meep::infinity = 1e20

this should be big enough for us

int meep::interrupt = 0

Incremented when ctrl-C is called.(starts with a zero).

double meep::kax [static]

int meep::kill_time = 2 [static]

double meep::ktrans [static]

int meep::m_for_J [static]

const double meep::nan = -1.23454321e123

ideally, a value never encountered in practice

const int meep::num_bandpts = 32

const int meep::NUM_FIELD_COMPONENTS = 15

const int meep::NUM_FIELD_TYPES = 4

const int meep::num_sphere_quad[3] = { 2, 12, 50 } [static]

ivec meep::one_ivec(ndim)

const double meep::pi = 3.141592653589793238462643383276

complex<double> meep::prefactor_static [static]

bool meep::quiet = false

if true, suppress all non-error messages from Meep

const double meep::sphere_quad[3][50][4] [static]

const char meep::symlink_name[] = "latest_output"

ivec meep::zero_ivec(ndim)

vec meep::zero_vec(ndim)