update parallel class

superdirectivity.cpp

@@ -262,17 +262,19 @@ Real sphere_ml_parallel::get_directivity_edipole_x(Real kz, const std::vector<Re
 	Complex ampl_in[3*_degree], ampl_out[3*_degree];
 	
 	if (kz <= kR[0]) {
-		// s-matrix of a the multilayer
+			// s-matrix of a the multilayer
 		sphere_sm(kR[0], eps[0], eps[1], rtml_out);
+		//cout << "sphere coeff:\n";
+		//for (int i=0; i<deg_max; ++i) {cout << i << " " << rtml_out[8*i+2] << " " << rtml_out[8*i+3] << endl;}
 		for (int k=1; k<ni; ++k) {
 			sphere_sm(kR[k], eps[k], eps[k+1], rt);
 			compose_sm(rtml_out, rt, rtml_out);
 		}
-		// dipole amplitudes
+			// dipole amplitudes
 		kzn = kz*sqrt(eps[0]);
 		if (arg(kzn) < -1e-14) kzn = -kzn;
 		get_edipole_amp_x(kzn, ampl_out, 1);
-		// output amplitudes
+			// output amplitudes
 		for (int n=0; n<deg_max; ++n) {
 			ampl_out[3*n] *= rtml_out[8*n+2]; // e, m = +-1
 			ampl_out[3*n+1] *= rtml_out[8*n+3]; // h, m = +-1
@@ -280,10 +282,10 @@ Real sphere_ml_parallel::get_directivity_edipole_x(Real kz, const std::vector<Re
 		}
 	}
 	else if (kz <= kR[kR.size()-1]) { // the dipole is located in a layer
-																		// determine a dipole layer
+			// determine a dipole layer
 		int kd = 0;
 		while (kz > kR[++kd]); // determine a first interface with radius larger than the dipole position
-													 // calculate s-matrices of multilayers inside and outside the dipole position:
+			// calculate s-matrices of multilayers inside and outside the dipole position:
 		sphere_sm(kR[0], eps[0], eps[1], rtml_in);
 		for (int k=1; k<kd; ++k) {
 			sphere_sm(kR[k], eps[k], eps[k+1], rt);
@@ -294,34 +296,34 @@ Real sphere_ml_parallel::get_directivity_edipole_x(Real kz, const std::vector<Re
 			sphere_sm(kR[k], eps[k], eps[k+1], rt);
 			compose_sm(rtml_out, rt, rtml_out);
 		}
-		// dipole amplitudes
+			// dipole amplitudes
 		kzn = kz*sqrt(eps[kd]);
 		if (arg(kzn) < -1e-14) kzn = -kzn;
 		get_edipole_amp_x(kzn, ampl_in, 0);
 		get_edipole_amp_x(kzn, ampl_out, 1);
-		// output amplitudes
+			// output self-consistent amplitudes
 		for (int n=0; n<deg_max; ++n) {
 			ampl_out[3*n] = ( ampl_in[3*n]*rtml_in[8*n+6] + ampl_out[3*n] ) * rtml_out[8*n+2]
-				/ (Real(1) - rtml_in[8*n+6]*rtml_out[8*n]); // e, m = 1
+										/ (Real(1) - rtml_in[8*n+6]*rtml_out[8*n]); // e, m = 1
 			ampl_out[3*n+1] = ( ampl_in[3*n+1]*rtml_in[8*n+7] + ampl_out[3*n+1] ) * rtml_out[8*n+3]
-				/ (Real(1) - rtml_in[8*n+7]*rtml_out[8*n+1]); // h, m = 1
+										/ (Real(1) - rtml_in[8*n+7]*rtml_out[8*n+1]); // h, m = 1
 			ampl_out[3*n+2] = ( ampl_in[3*n+2]*rtml_in[8*n+7] + ampl_out[3*n+2] ) * rtml_out[8*n+3]
-				/ (Real(1) - rtml_in[8*n+7]*rtml_out[8*n+1]); // h, m = 0
+										/ (Real(1) - rtml_in[8*n+7]*rtml_out[8*n+1]); // h, m = 0
 		}
 	}
 	else { // the dipole is located outside the multilayer
-				 // s-matrix of a the multilayer
+			// s-matrix of a the multilayer
 		sphere_sm(kR[0], eps[0], eps[1], rtml_out);
 		for (int k=1; k<ni; ++k) {
 			sphere_sm(kR[k], eps[k], eps[k+1], rt);
 			compose_sm(rtml_out, rt, rtml_out);
 		}
-		// dipole amplitudes
+			// dipole amplitudes
 		kzn = kz*sqrt(eps[ni]);
 		if (arg(kzn) < -1e-14) kzn = -kzn;
 		get_edipole_amp_x(kzn, ampl_in, 0);
 		get_edipole_amp_x(kzn, ampl_out, 1);
-		// output amplitudes
+			// output amplitudes
 		for (int n=0; n<deg_max; ++n) {
 			ampl_out[3*n] += ampl_in[3*n] * rtml_out[8*n+6]; // e, m = +-1
 			ampl_out[3*n+1] += ampl_in[3*n+1] * rtml_out[8*n+7]; // h, m = +-1

superdirectivity.h

@@ -14,7 +14,7 @@
 Real get_directivity_sphml_dZx_epsRe(std::vector<Real> &param, void *cls);
 Real get_directivity_sphml_dZx_epsRe_parallel(std::vector<Real> &param, void *cls);
 
-const int _degree = 15;
+const int _degree = 30;
 
 class sphere_ml {
 public: