Use abs() function to calculate s0mag

src/cluster/cluster.cpp

@@ -368,7 +368,7 @@ void cluster() {
 						s0 = c1->fsas[i] * exri;
 						qschu = s0.imag() * csch;
 						pschu = s0.real() * csch;
-						s0mag = sqrt((s0.real() + s0.imag()) * (s0.real() - s0.imag())) * cs0;
+						s0mag = abs(s0) * cs0;
 						fprintf(output, "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschu, pschu, s0mag);
 						double rapr = c1->sexs[i] - gaps[i];
 						double cosav = gaps[i] / c1->sscs[i];
@@ -382,7 +382,7 @@ void cluster() {
 				s0 = c3->tfsas * exri;
 				qschu = s0.imag() * csch;
 				pschu = s0.real() * csch;
-				s0mag = sqrt((s0.real() + s0.imag()) * (s0.real() - s0.imag())) * cs0;
+				s0mag = abs(s0) * cs0;
 				fprintf(output, "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschu, pschu, s0mag);
 				tppoan.write(reinterpret_cast<char *>(&vk), sizeof(double));
 				pcrsm0(vk, exri, inpol, c1, c1ao, c4);
@@ -731,7 +731,7 @@ void cluster() {
 									s0 = c1ao->fsac[ilr290 - 1][ilr290 - 1] * exri;
 									double qschu = s0.imag() * csch;
 									double pschu = s0.real() * csch;
-									s0mag = sqrt((s0.real() + s0.imag()) * (s0.real() - s0.imag())) * cs0;
+									s0mag = abs(s0) * cs0;
 									double refinr = c1ao->fsac[ilr290 - 1][ilr290 - 1].real() / c3->tfsas.real();
 									double extcor = c1ao->fsac[ilr290 - 1][ilr290 - 1].imag() / c3->tfsas.imag();
 									if (inpol == 0) {

src/sphere/sphere.cpp

@@ -277,7 +277,6 @@ void sphere() {
 					}
 				}
 				double cs0 = 0.25 * vk * vk * vk / half_pi;
-				//printf("DEBUG: cs0 = %lE\n", cs0);
 				sscr0(tfsas, nsph, gconf->l_max, vk, exri, c1);
 				printf("DEBUG: TFSAS = (%lE,%lE)\n", tfsas.real(), tfsas.imag());
 				double sqk = vk * vk * sconf->exdc;
@@ -317,15 +316,10 @@ void sphere() {
 						);
 						fprintf(output, "  FSAS=%15.7lE%15.7lE\n", c1->fsas[i170 - 1].real(), c1->fsas[i170 - 1].imag());
 						double csch = 2.0 * vk * sqsfi / c1->gcsv[i170 -1];
-						//printf("DEBUG: csch = %lE\n", csch);
 						s0 = c1->fsas[i170 - 1] * exri;
-						//printf("DEBUG: s0 = (%lE,%lE)\n", s0.real(), s0.imag());
 						double qschu = csch * s0.imag();
-						//printf("DEBUG: qschu = %lE\n", qschu);
 						double pschu = csch * s0.real();
-						//printf("DEBUG: pschu = %lE\n", pschu);
-						double s0mag = cs0 * sqrt((s0.real() + s0.imag()) * (s0.real() - s0.imag()));
-						//printf("DEBUG: s0mag = %lE\n", s0mag);
+						double s0mag = cs0 * abs(s0);
 						fprintf(
 								output,
 								"  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
@@ -364,7 +358,7 @@ void sphere() {
 					s0 = tfsas * exri;
 					double qschu = csch * s0.imag();
 					double pschu = csch * s0.real();
-					double s0mag = cs0 * sqrt((s0.real() + s0.imag()) * (s0.real() - s0.imag()));
+					double s0mag = cs0 * abs(s0);
 					fprintf(
 							output,
 							"  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",