Relax rules on needs_P for gyrotropic media (Bug#1117)

seewhydee · Apr 22, 2020 · 6f0928f · 6f0928f
1 parent ed3f3dd
commit 6f0928f
Show file tree

Hide file tree

Showing 2 changed files with 6 additions and 5 deletions.
diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp
@@ -1190,6 +1190,7 @@ static bool susceptibility_equiv(const susceptibility *o0, const susceptibility
   if (!vector3_equal(o0->bias, o->bias)) return 0;
   if (o0->frequency != o->frequency) return 0;
   if (o0->gamma != o->gamma) return 0;
+  if (o0->alpha != o->alpha) return 0;
   if (o0->noise_amp != o->noise_amp) return 0;
   if (o0->drude != o->drude) return 0;
   if (o0->saturated_gyrotropy != o->saturated_gyrotropy) return 0;

diff --git a/src/susceptibility.cpp b/src/susceptibility.cpp
@@ -434,7 +434,7 @@ void *gyrotropic_susceptibility::copy_internal_data(void *data) const {
 bool gyrotropic_susceptibility::needs_P(component c, int cmp, realnum *W[NUM_FIELD_COMPONENTS][2]) const {
   if (!is_electric(c) && !is_magnetic(c)) return false;
   direction d0 = component_direction(c);
-  return (d0 == X || d0 == Y || d0 == Z) && sigma[c][d0] && W[c][cmp];
+  return (d0 == X || d0 == Y || d0 == Z) && W[c][cmp];
 }
 
 // Similar to the OFFDIAG macro, but without averaging sigma.
@@ -470,8 +470,8 @@ void gyrotropic_susceptibility::update_P(realnum *W[NUM_FIELD_COMPONENTS][2],
                                  invdet * (gd * gd + gz * gz)}};
 
       FOR_COMPONENTS(c) DOCMP2 {
-        if (d->P[c][cmp][0]) {
-          const direction d0 = component_direction(c);
+	const direction d0 = component_direction(c);
+        if (d->P[c][cmp][0] && sigma[c][d0]) {
           const realnum *w0 = W[c][cmp], *s = sigma[c][d0];
 
           if (!w0 || !s || (d0 != X && d0 != Y && d0 != Z))
@@ -531,8 +531,8 @@ void gyrotropic_susceptibility::update_P(realnum *W[NUM_FIELD_COMPONENTS][2],
                                  invdet * (gd * gd + gz * gz)}};
 
       FOR_COMPONENTS(c) DOCMP2 {
-        if (d->P[c][cmp][0]) {
-          const direction d0 = component_direction(c);
+	const direction d0 = component_direction(c);
+        if (d->P[c][cmp][0] && sigma[c][d0]) {
           const realnum *w0 = W[c][cmp], *s = sigma[c][d0];
 
           if (!w0 || !s || (d0 != X && d0 != Y && d0 != Z))