diff --git a/src/sphere/edfb.cpp b/src/sphere/edfb.cpp index a40dea360e6cd01e076d8792493d47f9f64638f5..574c089e6245b7e82eb7b0c24fbaf3875e8d9d4b 100644 --- a/src/sphere/edfb.cpp +++ b/src/sphere/edfb.cpp @@ -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;