/*! \file Configuration.cpp
*/

#include <cmath>
#include <cstdio>
#include <fstream>
#include <string>
#include "include/List.h"
#include "include/Parsers.h"
#include "include/Configuration.h"

using namespace std;

GeometryConfiguration::GeometryConfiguration(
		int nsph, int lm, int _in_pol, int _npnt, int _npntts, int isam,
		double *x, double *y, double *z,
		double in_th_start, double in_th_step, double in_th_end,
		double sc_th_start, double sc_th_step, double sc_th_end,
		double in_ph_start, double in_ph_step, double in_ph_end,
		double sc_ph_start, double sc_ph_step, double sc_ph_end,
		int _jwtm
	) {
	number_of_spheres = nsph;
	l_max = lm;
	in_pol = _in_pol;
	npnt = _npnt;
	npntts = _npntts;
	meridional_type = isam;
	in_theta_start = in_th_start;
	in_theta_step = in_th_step;
	in_theta_end = in_th_end;
	in_phi_start = in_ph_start;
	in_phi_step = in_ph_step;
	in_phi_end = in_ph_end;
	sc_theta_start = sc_th_start;
	sc_theta_step = sc_th_step;
	sc_theta_end = sc_th_end;
	sc_phi_start = sc_ph_start;
	sc_phi_step = sc_ph_step;
	sc_phi_end = sc_ph_end;
	jwtm = _jwtm;
	sph_x = x;
	sph_y = y;
	sph_z = z;
}

GeometryConfiguration::~GeometryConfiguration() {
	delete[] sph_x;
	delete[] sph_y;
	delete[] sph_z;
}

GeometryConfiguration* GeometryConfiguration::from_legacy(string file_name) {
	int num_lines = 0;
	int last_read_line = 0;
	string *file_lines;
	try {
		file_lines = load_file(file_name, &num_lines);
	} catch (exception &ex) {
		throw OpenConfigurationFileException(file_name);
	}
	int nsph, lm, _in_pol, _npnt, _npntts, isam;
	sscanf(
			file_lines[last_read_line++].c_str(),
			" %d %d %d %d %d %d",
			&nsph, &lm, &_in_pol, &_npnt, &_npntts, &isam
	);
	double *x, *y, *z;
	x = new double[nsph];
	y = new double[nsph];
	z = new double[nsph];
	if (nsph == 1) {
		x[0] = 0.0;
		y[0] = 0.0;
		z[0] = 0.0;
	} else {
		for (int i = 0; i < nsph; i++) {
			double sph_x, sph_y, sph_z;
			int sph_x_exp, sph_y_exp, sph_z_exp;
			sscanf(
					file_lines[last_read_line++].c_str(),
					" %lf D%d %lf D%d %lf D%d",
					&sph_x, &sph_x_exp, &sph_y, &sph_y_exp, &sph_z, &sph_z_exp
			);
			x[i] = sph_x * pow(10.0, 1.0 * sph_x_exp);
			y[i] = sph_y * pow(10.0, 1.0 * sph_y_exp);
			z[i] = sph_z * pow(10.0, 1.0 * sph_z_exp);
		}
	}
	double in_th_start, in_th_end, in_th_step, sc_th_start, sc_th_end, sc_th_step;
	int in_th_start_exp, in_th_end_exp, in_th_step_exp, sc_th_start_exp, sc_th_end_exp, sc_th_step_exp;
	sscanf(
			file_lines[last_read_line++].c_str(),
			" %lf D%d %lf D%d %lf D%d %lf D%d %lf D%d %lf D%d",
			&in_th_start, &in_th_start_exp,
			&in_th_step, &in_th_step_exp,
			&in_th_end, &in_th_end_exp,
			&sc_th_start, &sc_th_start_exp,
			&sc_th_step, &sc_th_step_exp,
			&sc_th_end, &sc_th_end_exp
	);
	in_th_start *= pow(10.0, 1.0 * in_th_start_exp);
	in_th_step *= pow(10.0, 1.0 * in_th_step_exp);
	in_th_end *= pow(10.0, 1.0 * in_th_end_exp);
	sc_th_start *= pow(10.0, 1.0 * sc_th_start_exp);
	sc_th_step *= pow(10.0, 1.0 * sc_th_step_exp);
	sc_th_end *= pow(10.0, 1.0 * sc_th_end_exp);
	double in_ph_start, in_ph_end, in_ph_step, sc_ph_start, sc_ph_end, sc_ph_step;
	int in_ph_start_exp, in_ph_end_exp, in_ph_step_exp, sc_ph_start_exp, sc_ph_end_exp, sc_ph_step_exp;
	sscanf(
			file_lines[last_read_line++].c_str(),
			" %lf D%d %lf D%d %lf D%d %lf D%d %lf D%d %lf D%d",
			&in_ph_start, &in_ph_start_exp,
			&in_ph_step, &in_ph_step_exp,
			&in_ph_end, &in_ph_end_exp,
			&sc_ph_start, &sc_ph_start_exp,
			&sc_ph_step, &sc_ph_step_exp,
			&sc_ph_end, &sc_ph_end_exp
	);
	in_ph_start *= pow(10.0, 1.0 * in_ph_start_exp);
	in_ph_step *= pow(10.0, 1.0 * in_ph_step_exp);
	in_ph_end *= pow(10.0, 1.0 * in_ph_end_exp);
	sc_ph_start *= pow(10.0, 1.0 * sc_ph_start_exp);
	sc_ph_step *= pow(10.0, 1.0 * sc_ph_step_exp);
	sc_ph_end *= pow(10.0, 1.0 * sc_ph_end_exp);
	int _jwtm;
	sscanf(file_lines[last_read_line++].c_str(), " %d", &_jwtm);
	GeometryConfiguration *conf = new GeometryConfiguration(
			nsph, lm, _in_pol, _npnt, _npntts, isam,
			x, y, z,
			in_th_start, in_th_step, in_th_end,
			sc_th_start, sc_th_step, sc_th_end,
			in_ph_start, in_ph_step, in_ph_end,
			sc_ph_start, sc_ph_step, sc_ph_end,
			_jwtm
	);
	return conf;
}

ScattererConfiguration::ScattererConfiguration(
		int nsph,
		double *scale_vector,
		int nxi,
		string variable_name,
		int *iog_vector,
		double *ros_vector,
		int *nshl_vector,
		double **rcf_vector,
		int dielectric_func_type,
		complex<double> ***dc_matrix,
		bool is_external,
		double ex,
		double w,
		double x
		) {
	number_of_spheres = nsph;
	scale_vec = scale_vector;
	number_of_scales = nxi;
	reference_variable_name = variable_name;
	iog_vec = iog_vector;
	radii_of_spheres = ros_vector;
	nshl_vec = nshl_vector;
	rcf = rcf_vector;
	idfc = dielectric_func_type,
	dc0_matrix = dc_matrix;
	use_external_sphere = is_external;
	exdc = ex;
	wp = w;
	xip = x;
}

ScattererConfiguration::~ScattererConfiguration() {
	int max_ici = 0;
	for (int i = 1; i <= number_of_spheres; i++) {
		if (iog_vec[i - 1] < i) continue;
		int nsh = nshl_vec[i - 1];
		if (max_ici < (nsh + 1) / 2) max_ici = (nsh + 1) / 2;
	}
	for (int i = 0; i < max_ici; i++) {
		for (int j = 0; j < number_of_spheres; j++) {
			delete[] dc0_matrix[i][j];
		}
	}
	for (int i = 0; i < number_of_spheres; i++) {
		delete[] rcf[i];
	}
	delete[] nshl_vec;
	delete[] radii_of_spheres;
	delete[] iog_vec;
	delete[] scale_vec;
}

ScattererConfiguration* ScattererConfiguration::from_binary(string file_name, string mode) {
	int nsph;
	int *iog;
	double _exdc, _wp, _xip;
	int _idfc, nxi;
	int *nshl_vector;
	double *xi_vec;
	double *ros_vector;
	double **rcf_vector;
	complex<double> ***dc0m;
	if (mode.compare("LEGACY") == 0) { // Legacy mode was chosen.
		int max_ici = 0;
		fstream input;
		input.open(file_name.c_str(), ios::in | ios::binary);
		if (input.is_open()) {
			input.read(reinterpret_cast<char *>(&nsph), sizeof(int));
			iog = new int[nsph];
			for (int i = 0; i < nsph; i++) {
				input.read(reinterpret_cast<char *>(&(iog[i])), sizeof(int));
			}
			input.read(reinterpret_cast<char *>(&_exdc), sizeof(double));
			input.read(reinterpret_cast<char *>(&_wp), sizeof(double));
			input.read(reinterpret_cast<char *>(&_xip), sizeof(double));
			input.read(reinterpret_cast<char *>(&_idfc), sizeof(int));
			input.read(reinterpret_cast<char *>(&nxi), sizeof(int));
			try {
				xi_vec = new double[nxi];
			} catch (bad_alloc &ex) {
				throw UnrecognizedConfigurationException("Wrong parameter set: invalid number of scales " + nxi);
			}
			for (int i = 0; i < nxi; i++) {
				input.read(reinterpret_cast<char *>(&(xi_vec[i])), sizeof(double));
			}
			nshl_vector = new int[nsph];
			ros_vector = new double[nsph];
			rcf_vector = new double*[nsph];
			for (int i115 = 1; i115 <= nsph; i115++) {
				if (iog[i115 - 1] < i115) continue;
				input.read(reinterpret_cast<char *>(&(nshl_vector[i115 - 1])), sizeof(int));
				input.read(reinterpret_cast<char *>(&(ros_vector[i115 - 1])), sizeof(double));
				int nsh = nshl_vector[i115 - 1];
				if (max_ici < (nsh + 1) / 2) max_ici = (nsh + 1) / 2;
				try {
					rcf_vector[i115 - 1] = new double[nsh];
				} catch (bad_alloc &ex) {
					throw UnrecognizedConfigurationException("Wrong parameter set: invalid number of layers " + nsh);
				}
				for (int nsi = 0; nsi < nsh; nsi++) {
					input.read(reinterpret_cast<char *>(&(rcf_vector[i115 - 1][nsi])), sizeof(double));
				}
			}
			dc0m = new complex<double>**[max_ici];
			for (int dim1 = 0; dim1 < max_ici; dim1++) {
				dc0m[dim1] = new complex<double>*[nsph];
			    for (int dim2 = 0; dim2 < nsph; dim2++) {
			    	dc0m[dim1][dim2] = new complex<double>[nxi];
			    }
			}
			for (int jxi468 = 1; jxi468 <= nxi; jxi468++) {
				if (_idfc != 0 && jxi468 > 1) continue;
			    for (int i162 = 1; i162 <= nsph; i162++) {
			    	if (iog[i162 - 1] < i162) continue;
			    	int nsh = nshl_vector[i162 - 1];
			    	int ici = (nsh + 1) / 2; // QUESTION: is integer division really intended here?
			    	for (int i157 = 0; i157 < ici; i157++) {
			    		double dc0_real, dc0_img;
			    		input.read(reinterpret_cast<char *>(&dc0_real), sizeof(double));
			    		input.read(reinterpret_cast<char *>(&dc0_img), sizeof(double));
			    		dc0m[i157][i162 - 1][jxi468 - 1] = dc0_real + 1i * dc0_img;
			    	}
			    }
			}
			input.close();
		} else { // Opening of the input file did not succeed.
			throw OpenConfigurationFileException(file_name);
		}
	} else { // A different binary format was chosen.
		//TODO: this part is not yet implemented.
		//      Functions to write optimized file formats may be invoked here.
	}
	ScattererConfiguration *conf = new ScattererConfiguration(
			nsph,
			xi_vec,
			nxi,
			"XIV",
			iog,
			ros_vector,
			nshl_vector,
			rcf_vector,
			_idfc,
			dc0m,
			false,
			_exdc,
			_wp,
			_xip
	);
	return conf;
}

ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_name) {
	int num_lines = 0;
	int last_read_line = 0;
	string *file_lines;
	try {
		file_lines = load_file(dedfb_file_name, &num_lines);
	} catch (exception &ex) {
		throw OpenConfigurationFileException(dedfb_file_name);
	}
	int nsph, ies;
	int max_ici = 0;
	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;
	sscanf(
			file_lines[++last_read_line].c_str(),
			" %lf D%d %lf D%d %lf D%d %d %d %d %d",
			&_exdc, &exdc_exp,
			&_wp, &wp_exp,
			&_xip, &xip_exp,
			&_idfc, &nxi, &instpc, &insn
	);
	_exdc *= pow(10.0, 1.0 * 1.0 * exdc_exp);
	_wp *= pow(10.0, 1.0 * wp_exp);
	_xip *= pow(10.0, 1.0 * xip_exp);

	double *variable_vector;
	string variable_name;

	if (_idfc < 0) { // Diel. functions at XIP value and XI is scale factor
		variable_name = "XIV";
		if (instpc < 1) { // The variable vector is explicitly defined.
			double xi;
			int xi_exp;
			List<double> xi_vector;
			sscanf(file_lines[++last_read_line].c_str(), " %9lE D%d", &xi, &xi_exp);
			xi *= pow(10.0, 1.0 * xi_exp);
			xi_vector.set(0, xi);
			for (int jxi310 = 1; jxi310 < nxi; jxi310++) {
				sscanf(file_lines[++last_read_line].c_str(), " %9lE D%d", &xi, &xi_exp);
				xi *= pow(10.0, 1.0 * xi_exp);
				xi_vector.append(xi);
			}
			variable_vector = xi_vector.to_array();
		} else { // instpc >= 1: the variable vector is defined in steps
			double xi, xi_step;
			int xi_exp, xi_step_exp;
			sscanf(file_lines[++last_read_line].c_str(), " %9lE D%d %9lE D%d", &xi, &xi_exp, &xi_step, &xi_step_exp);
			xi *= pow(10.0, 1.0 * xi_exp);
			xi_step *= pow(10.0, 1.0 * xi_step_exp);
			variable_vector = new double[nxi];
			for (int jxi320 = 0; jxi320 < nxi; jxi320++) {
				variable_vector[jxi320] = xi;
				xi += xi_step;
			}
		}
	} else { // idfc >= 0
		variable_vector = 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, 1.0 * vs_exp);
				variable_vector[jxi_r] = vs;
			}
			switch (insn) {
			case 1: //xi vector definition
				variable_name = "XIV";
				break;
			case 2: //wave number vector definition
				variable_name = "WNS";
		    	break;
			case 3: //wavelength vector definition
				variable_name = "WLS";
				break;
			case 4: //pu vector definition
				variable_name = "PUS";
				break;
			case 5: //eV vector definition
				variable_name = "EVS";
				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, 1.0 * vs_exp);
			vs_step *= pow(10.0, 1.0 * vs_step_exp);
			for (int jxi110w = 0; jxi110w < nxi; jxi110w++) {
				variable_vector[jxi110w] = vs;
				vs += vs_step;
			}
			switch (insn) {
			case 1: //xi vector definition
				variable_name = "XIV";
				break;
			case 2: //wave number vector definition
				variable_name = "WNS";
				break;
			case 3: //wavelength vector definition
				variable_name = "WLS";
				break;
			case 4: //pu vector definition
				variable_name = "PUS";
				break;
			case 5: //eV vector definition
				variable_name = "EVS";
				break;
			}
		}
	}
	last_read_line++;
	int *iog_vector = new int[nsph];
	double *ros_vector = new double[nsph];
	double **rcf_vector = new double*[nsph];
	int *nshl_vector = new int[nsph];
	for (int i = 0; i < nsph; i++) {
		string read_format = "";
		for (int j = 0; j < i; j++) read_format += " %*d";
		read_format += " %d";
		sscanf(file_lines[last_read_line].c_str(), read_format.c_str(), (iog_vector + i));
	}
	for (int i113 = 1; i113 <= nsph; i113++) {
		int i_val, nsh;
	    double ros_val;
	    int ros_val_exp;
	    if (iog_vector[i113 - 1] < i113) continue;
	    sscanf(file_lines[++last_read_line].c_str(), " %d %lf D%d", &i_val, &ros_val, &ros_val_exp);
	    nshl_vector[i113 - 1] = i_val;
	    if (max_ici < (i_val + 1) / 2) max_ici = (i_val + 1) / 2;
	    ros_vector[i113 - 1] = ros_val * pow(10.0, 1.0 * ros_val_exp);
	    nsh = nshl_vector[i113 - 1];
	    if (i113 == 1) nsh += ies;
	    rcf_vector[i113 - 1] = new double[nsh];
	    for (int ns = 0; ns < nsh; ns++) {
	    	double ns_rcf;
	    	int ns_rcf_exp;
	    	sscanf(file_lines[++last_read_line].c_str(), " %lf D%d", &ns_rcf, &ns_rcf_exp);
	    	rcf_vector[i113 -1][ns] = ns_rcf * pow(10.0, 1.0 * ns_rcf_exp);
	    }
	}
	complex<double> ***dc0m = new complex<double>**[max_ici];
	for (int dim1 = 0; dim1 < max_ici; dim1++) {
		dc0m[dim1] = new complex<double>*[nsph];
	    for (int dim2 = 0; dim2 < nsph; dim2++) {
	    	dc0m[dim1][dim2] = new complex<double>[nxi];
	    }
	}
	for (int jxi468 = 1; jxi468 <= nxi; jxi468++) {
		if (_idfc != 0 && jxi468 > 1) continue;
	    for (int i162 = 1; i162 <= nsph; i162++) {
	    	if (iog_vector[i162 - 1] < i162) continue;
	    	int nsh = nshl_vector[i162 - 1];
	    	int ici = (nsh + 1) / 2; // QUESTION: is integer division really intended here?
	    	if (i162 == 1) ici = ici + ies;
	    	for (int i157 = 0; i157 < ici; i157++) {
	    		double dc0_real, dc0_img;
	    		int dc0_real_exp, dc0_img_exp;
	    		sscanf(file_lines[++last_read_line].c_str(), " (%lf D%d, %lf D%d)", &dc0_real, &dc0_real_exp, &dc0_img, &dc0_img_exp);
	    		dc0_real *= pow(10.0, 1.0 * dc0_real_exp);
	    		dc0_img *= pow(10.0, 1.0 * dc0_img_exp);
	    		dc0m[i157][i162 - 1][jxi468 - 1] = dc0_real + 1i * dc0_img;
	    	}
	    }
	}

	ScattererConfiguration *config = new ScattererConfiguration(
			nsph,
			variable_vector,
			nxi,
			variable_name,
			iog_vector,
			ros_vector,
			nshl_vector,
			rcf_vector,
			_idfc,
			dc0m,
			(ies > 0),
			_exdc,
			_wp,
			_xip
	);
	return config;
}

void ScattererConfiguration::print() {
	int ies = (use_external_sphere)? 1 : 0;
	int max_ici = 0;
	printf("### CONFIGURATION DATA ###\n");
	printf("NSPH  = %d\n", number_of_spheres);
	printf("ROS   = [");
	for (int i = 0; i < number_of_spheres; i++) printf("\t%lg", radii_of_spheres[i]);
	printf("\t]\n");
	printf("IOG   = [");
	for (int i = 0; i < number_of_spheres; i++) printf("\t%d", iog_vec[i]);
	printf("\t]\n");
	printf("NSHL  = [");
	for (int i = 0; i < number_of_spheres; i++) printf("\t%d", nshl_vec[i]);
	printf("\t]\n");
	printf("RCF   = [\n");
	for (int i = 1; i <= number_of_spheres; i++) {
		int nsh = nshl_vec[i - 1];
	    if (i == 1) nsh += ies;
	    if (max_ici < (nsh + 1) / 2) max_ici = (nsh + 1) / 2;
		printf("         [");
	    for (int ns = 0; ns < nsh; ns++) {
	    	printf("\t%lg", rcf[i - 1][ns]);
	    }
	    printf("\t]\n");
	}
	printf("        ]\n");
	printf("SCALE = %s\n", reference_variable_name.c_str());
	printf("NXI   = %d\n", number_of_scales);
	printf("VEC   = [");
	for (int i = 0; i < number_of_scales; i++) printf("\t%lg", scale_vec[i]);
	printf("\t]\n");
	printf("DC0M  = [\n");
	for (int i = 0; i < max_ici; i++) {
		printf("         [\n");
	    for (int j = 0; j < number_of_spheres; j++) {
			printf("          [");
			for (int k = 0; k < number_of_scales; k++) {
				if (idfc != 0 and k > 0) continue;
		    	printf("\t%lg + i(%lg)", dc0_matrix[i][j][k].real(), dc0_matrix[i][j][k].imag());
			}
			printf("\t]\n");
	    }
	    printf("         ]\n");
	}
	printf("        ]\n");
}

void ScattererConfiguration::write_binary(string file_name, string mode) {
	const double two_pi = acos(0.0) * 4.0;
	const double evc = 6.5821188e-16;
	int max_ici = 0;
	if (mode.compare("LEGACY") == 0) { // Legacy mode was chosen.
		fstream output;
		int ies = (use_external_sphere)? 1 : 0;
		double *xi_vec;
		if (reference_variable_name.compare("XIV") == 0) xi_vec = scale_vec;
		else {
			xi_vec = new double[number_of_scales];
			if (reference_variable_name.compare("WNS") == 0) {
				for (int i = 0; i < number_of_scales; i++)
					xi_vec[i] = 3.0e8 * scale_vec[i] / wp;
			} else if (reference_variable_name.compare("WLS") == 0) {
				for (int i = 0; i < number_of_scales; i++) {
					double wn = two_pi / scale_vec[i];
					xi_vec[i] = 3.0e8 * wn / wp;
				}
			} else if (reference_variable_name.compare("PUS") == 0) {
				for (int i = 0; i < number_of_scales; i++)
					xi_vec[i] = scale_vec[i] / wp;
			} else if (reference_variable_name.compare("EVS") == 0) {
				for (int i = 0; i < number_of_scales; i++) {
					double pu = scale_vec[i] / evc;
					xi_vec[i] = pu / wp;
				}
			} else {
				throw UnrecognizedConfigurationException(
						"Wrong parameter set: unrecognized scale type "
						+ reference_variable_name
				);
			}
		}
		output.open(file_name.c_str(), ios::out | ios::binary);
		output.write(reinterpret_cast<char *>(&number_of_spheres), sizeof(int));
		for (int i = 0; i < number_of_spheres; i++)
			output.write(reinterpret_cast<char *>(&(iog_vec[i])), sizeof(int));
		output.write(reinterpret_cast<char *>(&exdc), sizeof(double));
		output.write(reinterpret_cast<char *>(&wp), sizeof(double));
		output.write(reinterpret_cast<char *>(&xip), sizeof(double));
		output.write(reinterpret_cast<char *>(&idfc), sizeof(int));
		output.write(reinterpret_cast<char *>(&number_of_scales), sizeof(int));
		for (int i = 0; i < number_of_scales; i++)
			output.write(reinterpret_cast<char *>(&(xi_vec[i])), sizeof(double));
		for (int i115 = 1; i115 <= number_of_spheres; i115++) {
			if (iog_vec[i115 - 1] < i115) continue;
			output.write(reinterpret_cast<char *>(&(nshl_vec[i115 - 1])), sizeof(int));
			output.write(reinterpret_cast<char *>(&(radii_of_spheres[i115 - 1])), sizeof(double));
			int nsh = nshl_vec[i115 - 1];
			if (i115 == 1) nsh += ies;
			if (max_ici < (nsh + 1) / 2) max_ici = (nsh + 1) / 2;
			for (int nsi = 0; nsi < nsh; nsi++)
				output.write(reinterpret_cast<char *>(&(rcf[i115 - 1][nsi])), sizeof(double));
		}
		for (int jxi468 = 1; jxi468 <= number_of_scales; jxi468++) {
			if (idfc != 0 && jxi468 > 1) continue;
		    for (int i162 = 1; i162 <= number_of_spheres; i162++) {
		    	if (iog_vec[i162 - 1] < i162) continue;
		    	int nsh = nshl_vec[i162 - 1];
		    	int ici = (nsh + 1) / 2; // QUESTION: is integer division really intended here?
		    	if (i162 == 1) ici = ici + ies;
		    	for (int i157 = 0; i157 < ici; i157++) {
		    		double dc0_real, dc0_img;
		    		dc0_real = dc0_matrix[i157][i162 - 1][jxi468 - 1].real();
		    		dc0_img = dc0_matrix[i157][i162 - 1][jxi468 - 1].imag();
		    		// The FORTRAN code writes the complex numbers as a 16-byte long binary stream.
		    		// Here we assume that the 16 bytes are equally split in 8 bytes to represent the
		    		// real part and 8 bytes to represent the imaginary one.
		    		output.write(reinterpret_cast<char *>(&dc0_real), sizeof(double));
		    		output.write(reinterpret_cast<char *>(&dc0_img), sizeof(double));
		    	}
		    }
		}
		output.close();
	}
}

void ScattererConfiguration::write_formatted(string file_name) {
	const double evc = 6.5821188e-16;
	const double two_pi = acos(0.0) * 4.0;
	double *xi_vec;
	int ici;
	int ies = (use_external_sphere)? 1: 0;
	FILE *output = fopen(file_name.c_str(), "w");
	int scale_type = -1;
	if (reference_variable_name.compare("XIV") == 0) scale_type = 0;
	else if (reference_variable_name.compare("WNS") == 0) scale_type = 1;
	else if (reference_variable_name.compare("WLS") == 0) scale_type = 2;
	else if (reference_variable_name.compare("PUS") == 0) scale_type = 3;
	else if (reference_variable_name.compare("EVS") == 0) scale_type = 4;
	if (idfc >= 0) { // Dielectric functions are constant or they depend on XI
		double  *pu_vec, *ev_vec, *wn_vec, *wl_vec;
	    xi_vec = new double[number_of_scales];
	    pu_vec = new double[number_of_scales];
	    ev_vec = new double[number_of_scales];
	    wn_vec = new double[number_of_scales];
	    wl_vec = new double[number_of_scales];
		switch (scale_type) {
		case 0:
			fprintf(output, "  JXI     XIV          WNS          WLS          PUS          EVS\n");
			for (int i = 0; i < number_of_scales; i++) {
				xi_vec[i] = scale_vec[i];
				pu_vec[i] = xi_vec[i] * wp;
				ev_vec[i] = pu_vec[i] * evc;
				wn_vec[i] = pu_vec[i] / 3.0e8;
				wl_vec[i] = two_pi / wn_vec[i];
				fprintf(
						output,
						"%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n",
						(i + 1),
						xi_vec[i],
						wn_vec[i],
						wl_vec[i],
						pu_vec[i],
						ev_vec[i]
				);
			}
			break;
		case 1:
			fprintf(output, "  JXI     WNS          WLS          PUS          EVS          XIV\n");
			for (int i = 0; i < number_of_scales; i++) {
				wn_vec[i] = scale_vec[i];
				wl_vec[i] = two_pi / wn_vec[i];
				xi_vec[i] = 3.0e8 * wn_vec[i] / wp;
				pu_vec[i] = xi_vec[i] * wp;
				ev_vec[i] = pu_vec[i] * evc;
				fprintf(
						output,
						"%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n",
						(i + 1),
						wn_vec[i],
						wl_vec[i],
						pu_vec[i],
						ev_vec[i],
						xi_vec[i]
				);
			}
			break;
		case 2:
			fprintf(output, "  JXI     WLS          WNS          PUS          EVS          XIV\n");
			for (int i = 0; i < number_of_scales; i++) {
				wl_vec[i] = scale_vec[i];
				wn_vec[i] = two_pi / wl_vec[i];
				xi_vec[i] = 3.0e8 * wn_vec[i] / wp;
				pu_vec[i] = xi_vec[i] * wp;
				ev_vec[i] = pu_vec[i] * evc;
				fprintf(
						output,
						"%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n",
						(i + 1),
						wl_vec[i],
						wn_vec[i],
						pu_vec[i],
						ev_vec[i],
						xi_vec[i]
				);
			}
			break;
		case 3:
			fprintf(output, "  JXI     PUS          WNS          WLS          EVS          XIV\n");
			for (int i = 0; i < number_of_scales; i++) {
				pu_vec[i] = scale_vec[i];
				xi_vec[i] = pu_vec[i] / wp;
				wn_vec[i] = pu_vec[i] / 3.0e8;
				wl_vec[i] = two_pi / wn_vec[i];
				ev_vec[i] = pu_vec[i] * evc;
				fprintf(
						output,
						"%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n",
						(i + 1),
						pu_vec[i],
						wn_vec[i],
						wl_vec[i],
						ev_vec[i],
						xi_vec[i]
				);
			}
			break;
		case 4:
			fprintf(output, "  JXI     EVS          WNS          WLS          PUS          XIV\n");
			for (int i = 0; i < number_of_scales; i++) {
				ev_vec[i] = scale_vec[i];
				pu_vec[i] = ev_vec[i] / evc;
				xi_vec[i] = pu_vec[i] / wp;
				wn_vec[i] = pu_vec[i] / 3.0e8;
				wl_vec[i] = two_pi / wn_vec[i];
				fprintf(
						output,
						"%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n",
						(i + 1),
						ev_vec[i],
						wn_vec[i],
						wl_vec[i],
						pu_vec[i],
						xi_vec[i]
				);
			}
			break;
		default:
			throw UnrecognizedConfigurationException(
					"Wonrg parameter set: unrecognized scale variable type "
					+ reference_variable_name
			);
			break;
		}
	} else { // idfc < 0, Dielectric functions are at XIP and XI is scale for dimensions
		double pu, wn;
		xi_vec = scale_vec;
		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);
		fprintf(output, "%13.4lE", two_pi / wn);
		fprintf(output, "%13.4lE", pu);
		fprintf(output, "%13.4lE\n", pu * evc);
		fprintf(output, "  SCALE FACTORS XI\n");
		for (int i = 0; i < number_of_scales; i++)
			fprintf(output, "%5d%13.4lE\n", (i + 1), xi_vec[i]);
	}
	if (idfc != 0) {
		fprintf(output, "  DIELECTRIC CONSTANTS\n");
	    for (int i473 = 1; i473 <= number_of_spheres; i473++) {
	    	if (iog_vec[i473 - 1] != i473) continue;
	    	ici = (nshl_vec[i473 - 1] + 1) / 2;
	    	if (i473 == 1) ici += ies;
	    	fprintf(output, " SPHERE N. %4d\n", i473);
	    	for (int ic472 = 0; ic472 < ici; ic472++) {
	    		double dc0_real = dc0_matrix[ic472][i473 - 1][0].real(), dc0_img = dc0_matrix[ic472][i473 - 1][0].imag();
	    		fprintf(output, "%5d %12.4lE%12.4lE\n", (ic472 + 1), dc0_real, dc0_img);
	    	}
	    }
	  } else {
		  fprintf(output, "  DIELECTRIC FUNCTIONS\n");
		  for (int i478 = 1; i478 <= number_of_spheres; i478++) {
			  if (iog_vec[i478 - 1] != i478) continue;
			  ici = (nshl_vec[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 =  %3s\n", ic477, reference_variable_name.c_str());
				  for (int jxi476 = 0; jxi476 < number_of_scales; jxi476++) {
					  double dc0_real = dc0_matrix[ic477 - 1][i478 - 1][jxi476].real();
					  double dc0_img = dc0_matrix[ic477 - 1][i478 - 1][jxi476].imag();
					  fprintf(output, "%5d (%12.4lE,%12.4lE)\n", (jxi476 + 1), dc0_real, dc0_img);
				  }
			  }
		  }
	  }
	fclose(output);
}