maxwell_eps.c

MIT开发出来的计算光子晶体的软件

/* Copyright (C) 1999, 2000, 2001, 2002, Massachusetts Institute of Technology.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include "../config.h"
#include <check.h>
#include <mpiglue.h>
#include <mpi_utils.h>

#include "maxwell.h"

/**************************************************************************/

/* Lapack eigenvalue functions */
#ifdef SCALAR_SINGLE_PREC
#  define HEEV F77_FUNC(cheev,CHEEV)
#  define SYEV F77_FUNC(ssyev,SSYEV)
#else
#  define HEEV F77_FUNC(zheev,ZHEEV)
#  define SYEV F77_FUNC(dsyev,DSYEV)
#endif
extern void HEEV(char *, char *, int *, scalar *, int *, real *,
		 scalar *, int *, real *, int *);
extern void SYEV(char *, char *, int *, real *, int *, real *,
		 real *, int *, int *);

/* compute the 3 real eigenvalues of the matrix V */
void maxwell_sym_matrix_eigs(real eigs[3], const symmetric_matrix *V)
{
     int n = 3, nw = 9, info;
#if defined(WITH_HERMITIAN_EPSILON)
     scalar_complex Vm[3][3], W[9], W2[9];
     CASSIGN_SCALAR(Vm[0][0], V->m00, 0);
     CASSIGN_SCALAR(Vm[1][1], V->m11, 0);
     CASSIGN_SCALAR(Vm[2][2], V->m22, 0);
     Vm[0][1] = V->m01; CASSIGN_CONJ(Vm[1][0], V->m01);
     Vm[0][2] = V->m02; CASSIGN_CONJ(Vm[2][0], V->m02);
     Vm[1][2] = V->m12; CASSIGN_CONJ(Vm[2][1], V->m12);
     HEEV("V", "U", &n, &Vm[0][0], &n, eigs, W, &nw, W2, &info);
#else
     real Vm[3][3], W[9];
     Vm[0][0] = V->m00;
     Vm[1][1] = V->m11;
     Vm[2][2] = V->m22;
     Vm[0][1] = Vm[1][0] = V->m01;
     Vm[0][2] = Vm[2][0] = V->m02;
     Vm[1][2] = Vm[2][1] = V->m12;
     SYEV("V", "U", &n, &Vm[0][0], &n, eigs, W, &nw, &info);
#endif
     CHECK(info >= 0, "invalid argument in heev");
     CHECK(info <= 0, "failure to converge in heev");
}

/* Set Vinv = inverse of V, where both V and Vinv are real-symmetric
   (or possibly complex-Hermitian) matrices. */
void maxwell_sym_matrix_invert(symmetric_matrix *Vinv, 
			       const symmetric_matrix *V)
{
     real m00 = V->m00, m11 = V->m11, m22 = V->m22;

#if defined(WITH_HERMITIAN_EPSILON)
     scalar_complex m01 = V->m01, m02 = V->m02, m12 = V->m12;

     if (m01.re == 0.0 && m01.im == 0.0 &&
	 m02.re == 0.0 && m02.im == 0.0 &&
	 m12.re == 0.0 && m12.im == 0.0) {
	  /* optimize common case of a diagonal matrix: */
	  Vinv->m00 = 1.0 / m00;
	  Vinv->m11 = 1.0 / m11;
	  Vinv->m22 = 1.0 / m22;
	  CASSIGN_ZERO(Vinv->m01);
	  CASSIGN_ZERO(Vinv->m02);
	  CASSIGN_ZERO(Vinv->m12);
     }
     else {
	  double detinv;
	  
	  /* compute the determinant: */
	  detinv = m00*m11*m22 - m11*CSCALAR_NORMSQR(m02) -
	       CSCALAR_NORMSQR(m01)*m22 - CSCALAR_NORMSQR(m12)*m00 +
	       2.0 * ((m01.re * m12.re - m01.im * m12.im) * m02.re +
		      (m01.re * m12.im + m01.im * m12.re) * m02.im);
	  
	  CHECK(detinv != 0.0, "singular 3x3 matrix");
	  
	  detinv = 1.0/detinv;
	  
	  Vinv->m00 = detinv * (m11*m22 - CSCALAR_NORMSQR(m12));
	  Vinv->m11 = detinv * (m00*m22 - CSCALAR_NORMSQR(m02));
	  Vinv->m22 = detinv * (m11*m00 - CSCALAR_NORMSQR(m01));
	  
	  CASSIGN_SCALAR(Vinv->m02,
			 detinv * (m01.re*m12.re-m01.im*m12.im - m11*m02.re),
			 -detinv*(-m01.re*m12.im-m01.im*m12.re + m11*m02.im));

	  CASSIGN_SCALAR(Vinv->m01,
			 detinv * (m12.re*m02.re+m12.im*m02.im - m22*m01.re),
			 -detinv * (m12.im*m02.re-m12.re*m02.im + m22*m01.im));

	  CASSIGN_SCALAR(Vinv->m12,
			 detinv * (m01.re*m02.re+m01.im*m02.im - m00*m12.re),
			 -detinv * (m01.im*m02.re-m01.re*m02.im + m00*m12.im));
     }
     
#else /* real matrix */
     real m01 = V->m01, m02 = V->m02, m12 = V->m12;

     if (m01 == 0.0 && m02 == 0.0 && m12 == 0.0) {
	  /* optimize common case of a diagonal matrix: */
	  Vinv->m00 = 1.0 / m00;
	  Vinv->m11 = 1.0 / m11;
	  Vinv->m22 = 1.0 / m22;
	  Vinv->m01 = Vinv->m02 = Vinv->m12 = 0.0;
     }
     else {
	  double detinv;
	  
	  /* compute the determinant: */
	  detinv = m00*m11*m22 - m02*m11*m02 + 2.0 * m01*m12*m02 -
	       m01*m01*m22 - m12*m12*m00;
	  
	  CHECK(detinv != 0.0, "singular 3x3 matrix");
	  
	  detinv = 1.0/detinv;
	  
	  Vinv->m00 = detinv * (m11*m22 - m12*m12);
	  Vinv->m11 = detinv * (m00*m22 - m02*m02);
	  Vinv->m22 = detinv * (m11*m00 - m01*m01);
	  
	  Vinv->m02 = detinv * (m01*m12 - m11*m02);
	  Vinv->m01 = detinv * (m12*m02 - m01*m22);
	  Vinv->m12 = detinv * (m01*m02 - m00*m12);
     }
#endif /* real matrix */
}

/* Returns whether or not V is positive-definite. */
static int sym_matrix_positive_definite(symmetric_matrix *V)
{
     real det2, det3;
     real m00 = V->m00, m11 = V->m11, m22 = V->m22;

#if defined(WITH_HERMITIAN_EPSILON)
     scalar_complex m01 = V->m01, m02 = V->m02, m12 = V->m12;

     det2 = m00*m11 - CSCALAR_NORMSQR(m01);
     det3 = det2*m22 - m11*CSCALAR_NORMSQR(m02) - CSCALAR_NORMSQR(m12)*m00 +
	  2.0 * ((m01.re * m12.re - m01.im * m12.im) * m02.re +
		 (m01.re * m12.im + m01.im * m12.re) * m02.im);
#else /* real matrix */
     real m01 = V->m01, m02 = V->m02, m12 = V->m12;

     det2 = m00*m11 - m01*m01;
     det3 = det2*m22 - m02*m11*m02 + 2.0 * m01*m12*m02 - m12*m12*m00;
#endif /* real matrix */
     
     return (m00 > 0.0 && det2 > 0.0 && det3 > 0.0);
}

/**************************************************************************/

int check_maxwell_dielectric(maxwell_data *d,
			     int negative_epsilon_okp)
{
     int i, require_2d;

     require_2d = d->nz == 1 && (d->parity & (EVEN_Z_PARITY | ODD_Z_PARITY));

     for (i = 0; i < d->fft_output_size; ++i) {
	  if (!negative_epsilon_okp &&
	      !sym_matrix_positive_definite(d->eps_inv + i))
	       return 1;
	  if (require_2d) {
#if defined(WITH_HERMITIAN_EPSILON)
	       if (d->eps_inv[i].m02.re != 0.0 ||
		   d->eps_inv[i].m02.im != 0.0 ||
		   d->eps_inv[i].m12.re != 0.0 ||
		   d->eps_inv[i].m12.im != 0.0)
		    return 2;
#else /* real matrix */
	       if (d->eps_inv[i].m02 != 0.0 || d->eps_inv[i].m12 != 0.0)
		    return 2;
#endif /* real matrix */
	  }
     }
     return 0;
}

/**************************************************************************/

#define SHIFT3(x,y,z) {real SHIFT3_dummy = z; z = y; y = x; x = SHIFT3_dummy;}

/* Compute quadrature points and weights for integrating on the
   unit sphere.  x, y, z, and weight should be arrays of num_sq_pts
   values to hold the coordinates and weights of the quadrature points
   on output.   Currently, num_sq_pts = 12, 50, and 72 are supported. */
void spherical_quadrature_points(real *x, real *y, real *z,
				 real *weight, int num_sq_pts)
{
     int i,j,k,l, n = 0;
     real x0, y0, z0, w;

     if (num_sq_pts == 50) {
	  /* Computes quadrature points and weights for 50-point 11th degree
	     integration formula on a unit sphere.  This particular quadrature
	     formula has the advantage, for our purposes, of preserving the
	     symmetry group of an octahedraon (i.e. simple cubic symmetry, with
	     respect to the Cartesian xyz axes).  References:
	
	     A. D. McLaren, "Optimal Numerical Integration on a Sphere,"
	     Math. Comp. 17, pp. 361-383 (1963).
	     
	     Also in: Arthur H. Stroud, "Approximate Calculation of Multiple
	     Integrals" (Prentice Hall, 1971) (formula number U3:11-1). 

	     This code was written with the help of example code by
	     John Burkardt: 
	           http://www.psc.edu/~burkardt/src/stroud/stroud.html */
     
	  x0 = 1; y0 = z0 = 0;
	  w = 9216 / 725760.0;
	  for (i = 0; i < 2; ++i) {
	       x0 = -x0;
	       for (j = 0; j < 3; ++j) {
		    SHIFT3(x0,y0,z0);
		    x[n] = x0; y[n] = y0; z[n] = z0; weight[n++] = w;
	       }
	  }
	  
	  x0 = y0 = sqrt(0.5); z0 = 0;
	  w = 16384 / 725760.0;
	  for (i = 0; i < 2; ++i) {
	       x0 = -x0;
	       for (j = 0; j < 2; ++j) {
		    y0 = -y0;
		    for (k = 0; k < 3; ++k) {
			 SHIFT3(x0,y0,z0);
			 x[n] = x0; y[n] = y0; z[n] = z0; weight[n++] = w;
		    }
	       }
	  }
	  
	  x0 = y0 = z0 = sqrt(1.0 / 3.0);
	  w = 15309 / 725760.0;
	  for (i = 0; i < 2; ++i) {
	       x0 = -x0;
	       for (j = 0; j < 2; ++j) {
		    y0 = -y0;
		    for (k = 0; k < 2; ++k) {
			 z0 = -z0;
			 x[n] = x0; y[n] = y0; z[n] = z0; weight[n++] = w;
		    }
	       }
	  }
	  
	  x0 = y0 = sqrt(1.0 / 11.0); z0 = 3 * x0;
	  w = 14641 / 725760.0;
	  for (i = 0; i < 2; ++i) {
	       x0 = -x0;
	       for (j = 0; j < 2; ++j) {
		    y0 = -y0;
		    for (k = 0; k < 2; ++k) {
			 z0 = -z0;
			 for (l = 0; l < 3; ++l) {
			      SHIFT3(x0,y0,z0);
			      x[n] = x0; y[n] = y0; z[n] = z0; weight[n++] = w;
			 }
		    }
	       }
	  }
     }  
     else if (num_sq_pts == 72 || num_sq_pts == 12) {
	  /* As above (same references), but with a 72-point 14th degree
	     formula, this time with the symmetry group of an icosohedron.
	     (Stroud formula number U3:14-1.)  Alternatively, just use
	     the 12-point 5th degree formula consisting of the vertices
	     of a regular icosohedron. */