add support for multilevel-atomic susceptibility

libctl/meep.scm.in

@@ -119,19 +119,20 @@
 ; ****************************************************************
 ; Multilevel-atom nonlinear susceptibilities
 
-;(define-class transition no-parent
-;  (define-property from-level no-default 'integer non-negative?)
-;  (define-property to-level no-default 'integer non-negative?)
-;  (define-property transition-rate 0 'number) ; nonradiative rate (0 if none)
-;  (define-property frequency 0 'number) ; radiative frequency (0 if none)
-;  (define-property sigma-diag (vector3 1 1 1) 'vector3) ; per-transition sigma
-;  (define-property gamma 0 'number)) ; optical damping rate
-
-;(define (transition-time t) (transition-rate (/ t)))
-
-;(define-class multilevel-atom susceptibility
-;  (define-property initial-populations '() (make-list-type 'number))
-;  (define-property transitions '() (make-list-type 'transition)))
+(define-class transition no-parent
+ (define-property from-level no-default 'integer non-negative?)
+ (define-property to-level no-default 'integer non-negative?)
+ (define-property transition-rate 0 'number) ; nonradiative rate (0 if none)
+ (define-property frequency 0 'number) ; radiative frequency (0 if none)
+ (define-property sigma-diag (vector3 1 1 1) 'vector3) ; per-transition sigma
+ (define-property gamma 0 'number) ; optical damping rate
+ (define-property pumping-rate 0 'number)) ; pumping rate (0 if none)
+
+(define (transition-time t) (transition-rate (/ t)))
+
+(define-class multilevel-atom susceptibility
+ (define-property initial-populations '() (make-list-type 'number))
+ (define-property transitions '() (make-list-type 'transition)))
 
 ; ****************************************************************
 ; Add some predefined variables, for convenience:

libctl/structure.cpp

@@ -1199,11 +1199,9 @@ static bool susceptibility_equiv(const susceptibility *o0,
 				 const susceptibility *o)
 {
 if (o0->which_subclass != o->which_subclass) return 0;
-#if 0
 if (o0->which_subclass == susceptibility::MULTILEVEL_ATOM) {
 if (!multilevel_atom_equal(o0->subclass.multilevel_atom_data, o->subclass.multilevel_atom_data)) return 0;
 }
-#endif
 else if (o0->which_subclass == susceptibility::DRUDE_SUSCEPTIBILITY) {
 if (!drude_susceptibility_equal(o0->subclass.drude_susceptibility_data, o->subclass.drude_susceptibility_data)) return 0;
 }
@@ -1266,7 +1264,6 @@ void geom_epsilon::sigma_row(meep::component c, double sigrow[3],
     material_type_destroy(material);
 }
 
-#if 0
 /* make multilevel_susceptibility from scheme input data */
 static meep::susceptibility *make_multilevel_sus(const multilevel_atom *d) {
   if (!d || d->transitions.num_items == 0) return NULL;
@@ -1300,8 +1297,8 @@ static meep::susceptibility *make_multilevel_sus(const multilevel_atom *d) {
   for (int t = 0; t < d->transitions.num_items; ++t) {
     int i = d->transitions.items[t].from_level - minlev;
     int j = d->transitions.items[t].to_level - minlev;
-    Gamma[i*L+i] += d->transitions.items[t].transition_rate;
-    Gamma[j*L+i] -= d->transitions.items[t].transition_rate;
+    Gamma[i*L+i] += d->transitions.items[t].transition_rate + d->transitions.items[t].pumping_rate;
+    Gamma[j*L+i] -= d->transitions.items[t].transition_rate + d->transitions.items[t].pumping_rate;
   }
 
   // initial populations of each level
@@ -1348,7 +1345,6 @@ static meep::susceptibility *make_multilevel_sus(const multilevel_atom *d) {
 
   return s;
 }
-#endif
 
 // add a polarization to the list if it is not already there
 static pol *add_pol(pol *pols, const susceptibility *user_s)
@@ -1435,13 +1431,12 @@ void geom_epsilon::add_susceptibilities(meep::field_type ft,
 	}
 	break;
       }
-#if 0
       case susceptibility::MULTILEVEL_ATOM: {
 	multilevel_atom *d = p->user_s.subclass.multilevel_atom_data;
+	master_printf("multilevel susceptibility: number of items=%d\n",d->transitions.num_items);
 	sus = make_multilevel_sus(d);
 	break;
       }
-#endif
       default:
 	meep::abort("unknown susceptibility type");
     }

src/meep.hpp

@@ -220,7 +220,7 @@ class noisy_lorentzian_susceptibility : public lorentzian_susceptibility {
 
 class multilevel_susceptibility : public susceptibility {
 public:
-  multilevel_susceptibility() : L(0), T(0), Gamma(0), N0(0), alpha(0), omega(0), gamma(0) {}
+  multilevel_susceptibility() : L(0), T(0), Gamma(0), N0(0), alpha(0), omega(0), gamma(0), sigmat(0) {}
   multilevel_susceptibility(int L, int T,
 			    const realnum *Gamma,
 			    const realnum *N0,
@@ -266,7 +266,7 @@ class multilevel_susceptibility : public susceptibility {
 protected:
   int L; // number of atom levels
   int T; // number of optical transitions
-  realnum *Gamma; // LxL matrix of relaxation rates Gamma[i*L+j] from i -> j
+  realnum *Gamma; // LxL matrix of relaxation and pumping rates Gamma[i*L+j] from i -> j
   realnum *N0; // L initial populations
   realnum *alpha; // LxT matrix of transition coefficients 1/omega
   realnum *omega; // T transition frequencies

src/multilevel-atom.cpp

@@ -229,7 +229,7 @@ realnum *multilevel_susceptibility::cinternal_notowned_ptr(
 					int n, 
 					void *P_internal_data) const {
   multilevel_data *d = (multilevel_data *) P_internal_data;
-  if (!d->P[c][cmp] || inotowned < 0 || inotowned >= T) // never true
+  if (!d || !d->P[c][cmp] || inotowned < 0 || inotowned >= T)
     return NULL;
   return d->P[c][cmp][inotowned] + n;
 }
@@ -256,12 +256,12 @@ void multilevel_susceptibility::update_P
   realnum *Ntmp = d->Ntmp;
   LOOP_OVER_VOL_OWNED(gv, Centered, i) {
     realnum *N = d->N + i*L; // N at current point, to update
-    
+
     // Ntmp = (I - Gamma * dt/2) * N
     for (int l1 = 0; l1 < L; ++l1) {
-      Ntmp[l1] = (1.0 - Gamma[l1*L + l1]*dt2) * N[l1]; // diagonal term
-      for (int l2 = 0; l2 < l1; ++l2) Ntmp[l1] -= Gamma[l1*L+l2]*dt2 * N[l2];
-      for (int l2 = l1+1; l2 < L; ++l2) Ntmp[l1] -= Gamma[l1*L+l2]*dt2 * N[l2];
+      Ntmp[l1] = 0;
+      for (int l2 = 0; l2 < L; ++l2)
+	Ntmp[l1] += ((l1 == l2) - Gamma[l1*L+l2]*dt2) * N[l2];
     }
 
     // compute E*8 at point i
@@ -291,7 +291,7 @@ void multilevel_susceptibility::update_P
 	if (d->P[cdot[idot]][1]) {
 	  p = d->P[cdot[idot]][1][t]; pp = d->P_prev[cdot[idot]][1][t];
 	  dP = p[i]+p[i+o1[idot]]+p[i+o2[idot]]+p[i+o1[idot]+o2[idot]]
-	    + (pp[i]+pp[i+o1[idot]]+pp[i+o2[idot]]+pp[i+o1[idot]+o2[idot]]);
+	    - (pp[i]+pp[i+o1[idot]]+pp[i+o2[idot]]+pp[i+o1[idot]+o2[idot]]);
 	  EdP32 += dP * E8[idot][1];
 	}
       }
@@ -308,9 +308,9 @@ void multilevel_susceptibility::update_P
 
   // each P is updated as a damped harmonic oscillator
   for (int t = 0; t < T; ++t) {
-    const double omega2pi = 2*pi*omega[t], g2pi = gamma[t]*2*pi;
-    const double omega0dtsqr = omega2pi * omega2pi * dt * dt;
-    const double gamma1inv = 1 / (1 + g2pi*dt/2), gamma1 = (1 - g2pi*dt/2);
+    const double omega2pi = omega[t]*2*pi, g2pi = gamma[t]*2*pi;
+    const double omega0dtsqr = omega2pi*omega2pi*dt*dt;
+    const double gamma1inv = 1 / (1 + g2pi*dt2), gamma1 = (1 - g2pi*dt2);
 
     // figure out which levels this transition couples
     int lp = -1, lm = -1;
@@ -341,20 +341,19 @@ void multilevel_susceptibility::update_P
 	component c2 = direction_component(c, d2);
 	const realnum *w2 = W[c2][cmp];
 	const realnum *s2 = w2 ? sigma[c][d2] : NULL;
-	
-	if (s1 || s2) {
+
+	if (s1 || s2)
 	  abort("nondiagonal saturable gain is not yet supported");
-	}
 	else { // isotropic
 	  LOOP_OVER_VOL_OWNED(gv, c, i) {
 	    realnum pcur = p[i];
 	    const realnum *Ni = N + i*L;
 	    // dNi is population inversion for this transition
-	    double dNi = -0.25 * (Ni[lp]+Ni[lp+o1]+Ni[lp+o2]+Ni[lp+o1+o2]
-				  -Ni[lm]-Ni[lm+o1]-Ni[lm+o2]-Ni[lm+o1+o2]);
+	    double dNi = 0.25 * (Ni[lp]+Ni[lp+o1]+Ni[lp+o2]+Ni[lp+o1+o2]
+				 -Ni[lm]-Ni[lm+o1]-Ni[lm+o2]-Ni[lm+o1+o2]);
 	    p[i] = gamma1inv * (pcur * (2 - omega0dtsqr) 
-				- gamma1 * pp[i] 
-				+ omega0dtsqr * (st * s[i] * w[i])) * dNi;
+				- gamma1 * pp[i]
+				- omega0dtsqr * st * s[i] * w[i] * dNi);
 	    pp[i] = pcur;
 	  }
 	}