Complete migration of edfb.f to edfb.cpp (b525842a) · Commits · Giacomo Mulas / NP_TMcode

src/sphere/edfb.cpp

+273 −13

Original line number	Diff line number	Diff line
		@@ -50,19 +50,19 @@ int main(int argc, char **argv) {

		// Input file reading section
		int num_lines = 0;
		int last_read_line; // Keep track of where INXI left the input stream
		int last_read_line = 0; // Keep track of where the input stream was left
		string *file_lines = load_file("../../test_data/sphere/DEDFB", &num_lines);

		// Configuration code
		int nsph, ies;
		sscanf(file_lines[0].c_str(), " %d %d", &nsph, &ies);
		sscanf(file_lines[last_read_line].c_str(), " %d %d", &nsph, &ies);
		if (ies != 0) ies = 1;
		double exdc, wp, xip;
		int exdc_exp, wp_exp, xip_exp;
		int idfc, nxi, instpc, insn;
		int nsh;
		sscanf(
		file_lines[1].c_str(),
		file_lines[++last_read_line].c_str(),
		" %9lf D%d %9lf D%d %8lf D%d %d %d %d %d",
		&exdc, &exdc_exp,
		&wp, &wp_exp,
		@@ -78,28 +78,274 @@ int main(int argc, char **argv) {
		const double pigt = acos(0.0) * 4.0;
		const double evc = 6.5821188e-16;
		if (idfc >= 0) {
		printf("Not walked by input data.\n");
		// Not walked by default input data
		// This part of the code in not tested
		vss = new double[nxi];
		xiv = new double[nxi];
		pus = new double[nxi];
		evs = new double[nxi];
		wns = new double[nxi];
		wls = new double[nxi];
		if (instpc == 0) { // The variable vector is explicitly defined
		double vs;
		int vs_exp;
		for (int jxi_r = 0; jxi_r < nxi; jxi_r++) {
		sscanf(file_lines[++last_read_line].c_str(), " %lf D%d", &vs, &vs_exp);
		vs *= pow(10.0, vs_exp);
		vss[jxi_r] = vs;
		}
		switch (insn) {
		case 1: //xi vector definition
		vns[insn - 1] = "XIV";
		fprintf(output, " JXI XIV WNS WLS PUS EVS\n");
		for (int jxi210w = 0; jxi210w < nxi; jxi210w++) {
		xiv[jxi210w] = vss[jxi210w];
		pus[jxi210w] = xiv[jxi210w] * wp;
		evs[jxi210w] = pus[jxi210w] * evc;
		wns[jxi210w] = pus[jxi210w] / 3.0e8;
		wls[jxi210w] = pigt / wns[jxi210w];
		fprintf(
		output,
		"%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n",
		(jxi210w + 1),
		xiv[jxi210w],
		wns[jxi210w],
		wls[jxi210w],
		pus[jxi210w],
		evs[jxi210w]
		);
		}
		break;
		case 2: //wave number vector definition
		vns[insn - 1] = "WNS";
		fprintf(output, " JXI WNS WLS PUS EVS XIV\n");
		for (int jxi230w = 0; jxi230w < nxi; jxi230w++) {
		wns[jxi230w] = vss[jxi230w];
		wls[jxi230w] = pigt / wns[jxi230w];
		xiv[jxi230w] = 3.0e8 * wns[jxi230w] / wp;
		pus[jxi230w] = xiv[jxi230w] * wp;
		evs[jxi230w] = pus[jxi230w] * evc;
		fprintf(
		output,
		"%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n",
		(jxi230w + 1),
		wns[jxi230w],
		wls[jxi230w],
		pus[jxi230w],
		evs[jxi230w],
		xiv[jxi230w]
		);
		}
		break;
		case 3: //wavelength vector definition
		vns[insn - 1] = "WLS";
		fprintf(output, " JXI WLS WNS PUS EVS XIV\n");
		for (int jxi250w = 0; jxi250w < nxi; jxi250w++) {
		wls[jxi250w] = vss[jxi250w];
		wns[jxi250w] = pigt / wls[jxi250w];
		xiv[jxi250w] = 3.0e8 * wns[jxi250w] / wp;
		pus[jxi250w] = xiv[jxi250w] * wp;
		evs[jxi250w] = pus[jxi250w] * evc;
		fprintf(
		output,
		"%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n",
		(jxi250w + 1),
		wls[jxi250w],
		wns[jxi250w],
		pus[jxi250w],
		evs[jxi250w],
		xiv[jxi250w]
		);
		}
		break;
		case 4: //pu vector definition
		vns[insn - 1] = "PUS";
		fprintf(output, " JXI PUS WNS WLS EVS XIV\n");
		for (int jxi270w = 0; jxi270w < nxi; jxi270w++) {
		pus[jxi270w] = vss[jxi270w];
		xiv[jxi270w] = pus[jxi270w] / wp;
		wns[jxi270w] = pus[jxi270w] / 3.0e8;
		wls[jxi270w] = pigt / wns[jxi270w];
		evs[jxi270w] = pus[jxi270w] * evc;
		fprintf(
		output,
		"%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n",
		(jxi270w + 1),
		pus[jxi270w],
		wns[jxi270w],
		wls[jxi270w],
		evs[jxi270w],
		xiv[jxi270w]
		);
		}
		break;
		case 5: //eV vector definition
		vns[insn - 1] = "EVS";
		fprintf(output, " JXI EVS WNS WLS PUS XIV\n");
		for (int jxi290w = 0; jxi290w < nxi; jxi290w++) {
		evs[jxi290w] = vss[jxi290w];
		pus[jxi290w] = evs[jxi290w] / evc;
		xiv[jxi290w] = pus[jxi290w] / wp;
		wns[jxi290w] = pus[jxi290w] / 3.0e8;
		wls[jxi290w] = pigt / wns[jxi290w];
		fprintf(
		output,
		"%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n",
		(jxi290w + 1),
		evs[jxi290w],
		wns[jxi290w],
		wls[jxi290w],
		pus[jxi290w],
		xiv[jxi290w]
		);
		}
		break;
		}
		} else { // The variable vector needs to be computed in steps
		double vs, vs_step;
		int vs_exp, vs_step_exp;
		sscanf(file_lines[++last_read_line].c_str(), " %lf D%d %lf D%d", &vs, &vs_exp, &vs_step, &vs_step_exp);
		vs *= pow(10.0, vs_exp);
		vs_step *= pow(10.0, vs_step_exp);
		switch (insn) {
		case 1: //xi vector definition
		vns[insn - 1] = "XIV";
		fprintf(output, " JXI XIV WNS WLS PUS EVS\n");
		for (int jxi110w = 0; jxi110w < nxi; jxi110w++) {
		vss[jxi110w] = vs;
		xiv[jxi110w] = vss[jxi110w];
		pus[jxi110w] = xiv[jxi110w] * wp;
		wns[jxi110w] = pus[jxi110w] / 3.0e8;
		evs[jxi110w] = pus[jxi110w] * evc;
		wls[jxi110w] = pigt / wns[jxi110w];
		fprintf(
		output,
		"%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n",
		(jxi110w + 1),
		xiv[jxi110w],
		wns[jxi110w],
		wls[jxi110w],
		pus[jxi110w],
		evs[jxi110w]
		);
		vs += vs_step;
		}
		break;
		case 2: //wave number vector definition
		vns[insn - 1] = "WNS";
		fprintf(output, " JXI WNS WLS PUS EVS XIV\n");
		for (int jxi130w = 0; jxi130w < nxi; jxi130w++) {
		vss[jxi130w] = vs;
		wns[jxi130w] = vss[jxi130w];
		xiv[jxi130w] = 3.0e8 * wns[jxi130w] / wp;
		pus[jxi130w] = xiv[jxi130w] * wp;
		wls[jxi130w] = pigt / wns[jxi130w];
		evs[jxi130w] = pus[jxi130w] * evc;
		fprintf(
		output,
		"%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n",
		(jxi130w + 1),
		wns[jxi130w],
		wls[jxi130w],
		pus[jxi130w],
		evs[jxi130w],
		xiv[jxi130w]
		);
		vs += vs_step;
		}
		break;
		case 3: //wavelength vector definition
		vns[insn - 1] = "WLS";
		fprintf(output, " JXI WLS WNS PUS EVS XIV\n");
		for (int jxi150w = 0; jxi150w < nxi; jxi150w++) {
		vss[jxi150w] = vs;
		wls[jxi150w] = vss[jxi150w];
		wns[jxi150w] = pigt / wls[jxi150w];
		xiv[jxi150w] = 3.0e8 * wns[jxi150w] / wp;
		pus[jxi150w] = xiv[jxi150w] * wp;
		evs[jxi150w] = pus[jxi150w] * evc;
		fprintf(
		output,
		"%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n",
		(jxi150w + 1),
		wls[jxi150w],
		wns[jxi150w],
		pus[jxi150w],
		evs[jxi150w],
		xiv[jxi150w]
		);
		vs += vs_step;
		}
		break;
		case 4: //pu vector definition
		vns[insn - 1] = "PUS";
		fprintf(output, " JXI PUS WNS WLS EVS XIV\n");
		for (int jxi170w = 0; jxi170w < nxi; jxi170w++) {
		vss[jxi170w] = vs;
		pus[jxi170w] = vss[jxi170w];
		xiv[jxi170w] = pus[jxi170w] / wp;
		wns[jxi170w] = pus[jxi170w] / 3.0e8;
		wls[jxi170w] = pigt / wns[jxi170w];
		evs[jxi170w] = pus[jxi170w] * evc;
		fprintf(
		output,
		"%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n",
		(jxi170w + 1),
		pus[jxi170w],
		wns[jxi170w],
		wls[jxi170w],
		evs[jxi170w],
		xiv[jxi170w]
		);
		vs += vs_step;
		}
		break;
		case 5: //eV vector definition
		vns[insn - 1] = "EVS";
		fprintf(output, " JXI EVS WNS WLS PUS XIV\n");
		for (int jxi190w = 0; jxi190w < nxi; jxi190w++) {
		vss[jxi190w] = vs;
		evs[jxi190w] = vss[jxi190w];
		pus[jxi190w] = evs[jxi190w] / evc;
		xiv[jxi190w] = pus[jxi190w] / wp;
		wns[jxi190w] = pus[jxi190w] / 3.0e8;
		wls[jxi190w] = pigt / wns[jxi190w];
		fprintf(
		output,
		"%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n",
		(jxi190w + 1),
		evs[jxi190w],
		wns[jxi190w],
		wls[jxi190w],
		pus[jxi190w],
		xiv[jxi190w]
		);
		vs += vs_step;
		}
		break;
		}
		}
		// End of the untested code section.
		} else {
		if (instpc < 1) {
		// In this case the XI vector is explicitly defined.
		// Test input comes this way.
		vns[insn] = "XIV";
		List<double> xi_vector;
		double xi;
		double xi, pu, wn;
		int xi_exp;
		sscanf(file_lines[2].c_str(), " %9lE D%d", &xi, &xi_exp);
		vns[insn - 1] = "XIV";
		List<double> xi_vector;
		sscanf(file_lines[++last_read_line].c_str(), " %9lE D%d", &xi, &xi_exp);
		xi *= pow(10.0, xi_exp);
		xi_vector.set(0, xi);
		for (int jxi310 = 1; jxi310 < nxi; jxi310++) {
		sscanf(file_lines[2 + jxi310].c_str(), " %9lE D%d", &xi, &xi_exp);
		sscanf(file_lines[++last_read_line].c_str(), " %9lE D%d", &xi, &xi_exp);
		xi *= pow(10.0, xi_exp);
		xi_vector.append(xi);
		last_read_line = 2 + jxi310;
		}
		vss = xi_vector.to_array();
		xiv = xi_vector.to_array();
		double pu = xip + wp;
		double wn = pu / 3.0e8;
		pu = xip + wp;
		wn = pu / 3.0e8;
		fprintf(output, " XIP WN WL PU EV\n");
		fprintf(output, " %13.4lE", xip);
		fprintf(output, "%13.4lE", wn);
		@@ -198,7 +444,7 @@ int main(int argc, char **argv) {
		if (iog[i473 - 1] != i473) continue;
		ici = (nshl[i473 - 1] + 1) / 2;
		if (i473 == 1) ici += ies;
		fprintf(output, " SPHERE N. %d\n", i473);
		fprintf(output, " SPHERE N. %4d\n", i473);
		for (int ic472 = 0; ic472 < ici; ic472++) {
		double dc0_real = dc0m[ic472][i473 - 1][0].real(), dc0_img = dc0m[ic472][i473 - 1][0].imag();
		fprintf(output, "%5d %12.4lE%12.4lE\n", (ic472 + 1), dc0_real, dc0_img);
		@@ -206,6 +452,20 @@ int main(int argc, char **argv) {
		}
		} else {
		fprintf(output, " DIELECTRIC FUNCTIONS\n");
		for (int i478 = 1; i478 <= nsph; i478++) {
		if (iog[i478 - 1] != i478) continue;
		ici = (nshl[i478 - 1] + 1) / 2;
		if (i478 == 1) ici += ies;
		fprintf(output, " SPHERE N. %4d\n", i478);
		for (int ic477 = 1; ic477 <= ici; ic477++) {
		fprintf(output, " NONTRANSITION LAYER N. %2d , SCALE = %3c\n", ic477, vns[insn - 1].c_str());
		for (int jxi476 = 0; jxi476 < nxi; jxi476++) {
		double dc0_real = dc0m[ic477 - 1][i478 - 1][jxi476].real();
		double dc0_img = dc0m[ic477 - 1][i478 - 1][jxi476].imag();
		fprintf(output, "%5d (%12.4lE,%12.4lE)\n", (jxi476 + 1), dc0_real, dc0_img);
		}
		}
		}
		}
		fclose(output);
		return 0;