Use the configuration counter to handle sphere structural parameters (fa46f1a3) · Commits · Giacomo Mulas / NP_TMcode

src/cluster/cluster.cpp

+5 −2

Original line number	Original line	Diff line number	Diff line
	@@ -681,6 +681,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
	int jwtm = gconf->jwtm;		int jwtm = gconf->jwtm;
	np_int ndit = 2 * nsph * cid->c4->nlim;		np_int ndit = 2 * nsph * cid->c4->nlim;
	int isq, ibf;		int isq, ibf;
			int last_configuration;

	#ifdef USE_NVTX		#ifdef USE_NVTX
	nvtxRangePush("Prepare matrix calculation");		nvtxRangePush("Prepare matrix calculation");
	@@ -710,6 +711,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
	return jer;		return jer;
	// break; // rewrite this to go to the end of the function, to free locally allocated variables and return jer		// break; // rewrite this to go to the end of the function, to free locally allocated variables and return jer
	}		}
			last_configuration = 0;
	for (int i132 = 1; i132 <= nsph; i132++) {		for (int i132 = 1; i132 <= nsph; i132++) {
	int iogi = cid->c1->iog[i132 - 1];		int iogi = cid->c1->iog[i132 - 1];
	if (iogi != i132) {		if (iogi != i132) {
	@@ -718,7 +720,8 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
	cid->c1->rei[l123 - 1][i132 - 1] = cid->c1->rei[l123 - 1][iogi - 1];		cid->c1->rei[l123 - 1][i132 - 1] = cid->c1->rei[l123 - 1][iogi - 1];
	} // l123 loop		} // l123 loop
	} else {		} else {
	int nsh = cid->c1->nshl[i132 - 1];		last_configuration++;
			int nsh = cid->c1->nshl[last_configuration - 1];
	int ici = (nsh + 1) / 2;		int ici = (nsh + 1) / 2;
	if (idfc == 0) {		if (idfc == 0) {
	for (int ic = 0; ic < ici; ic++)		for (int ic = 0; ic < ici; ic++)
	@@ -731,7 +734,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
	}		}
	if (nsh % 2 == 0) cid->c2->dc0[ici] = exdc;		if (nsh % 2 == 0) cid->c2->dc0[ici] = exdc;
	dme(		dme(
	cid->c4->li, i132, npnt, npntts, vkarg, exdc, exri,		cid->c4->li, last_configuration, npnt, npntts, vkarg, exdc, exri,
	cid->c1, cid->c2, jer, lcalc, cid->arg		cid->c1, cid->c2, jer, lcalc, cid->arg
	);		);
	if (jer != 0) {		if (jer != 0) {

src/libnptm/Configuration.cpp

+74 −67

Original line number	Original line	Diff line number	Diff line
	@@ -429,7 +429,7 @@ ScattererConfiguration::ScattererConfiguration(const mixMPI *mpidata)
	MPI_Bcast(&_xip, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);		MPI_Bcast(&_xip, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
	_iog_vec = new int[_number_of_spheres]();		_iog_vec = new int[_number_of_spheres]();
	MPI_Bcast(_iog_vec, _number_of_spheres, MPI_INT, 0, MPI_COMM_WORLD);		MPI_Bcast(_iog_vec, _number_of_spheres, MPI_INT, 0, MPI_COMM_WORLD);
	_radii_of_spheres = new double[_number_of_spheres]();		_radii_of_spheres = new double[_configurations]();
	MPI_Bcast(_radii_of_spheres, _configurations, MPI_DOUBLE, 0, MPI_COMM_WORLD);		MPI_Bcast(_radii_of_spheres, _configurations, MPI_DOUBLE, 0, MPI_COMM_WORLD);
	_nshl_vec = new int[_configurations]();		_nshl_vec = new int[_configurations]();
	MPI_Bcast(_nshl_vec, _configurations, MPI_INT, 0, MPI_COMM_WORLD);		MPI_Bcast(_nshl_vec, _configurations, MPI_INT, 0, MPI_COMM_WORLD);
	@@ -513,6 +513,7 @@ ScattererConfiguration* ScattererConfiguration::from_binary(const std::string& f
	ScattererConfiguration* ScattererConfiguration::from_dedfb(const std::string& dedfb_file_name) {		ScattererConfiguration* ScattererConfiguration::from_dedfb(const std::string& dedfb_file_name) {
	int num_lines = 0;		int num_lines = 0;
	int last_read_line = 0;		int last_read_line = 0;
			int last_configuration;
	string *file_lines;		string *file_lines;
	regex re;		regex re;
	smatch m;		smatch m;
	@@ -687,32 +688,34 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(const std::string& de
	}		}
	}		}
	int index = 0;		int index = 0;
	double *ros_vector = new double[fnsph]();		double *ros_vector = new double[configurations]();
	int *nshl_vector = new int[fnsph]();		int *nshl_vector = new int[configurations]();
	double *rcf_vector = new double[fnsph];		double *rcf_vector = new double[configurations];
			last_configuration = 0;
	for (int i113 = 1; i113 <= fnsph; i113++) {		for (int i113 = 1; i113 <= fnsph; i113++) {
	if (iog_vector[i113 - 1] < i113) continue;		if (iog_vector[i113 - 1] < i113) continue;
			else last_configuration++;
	str_target = file_lines[++last_read_line];		str_target = file_lines[++last_read_line];
	re = regex("[0-9]+");		re = regex("[0-9]+");
	regex_search(str_target, m, re);		regex_search(str_target, m, re);
	nshl_vector[i113 - 1] = stoi(m.str());		nshl_vector[last_configuration - 1] = stoi(m.str());
	str_target = m.suffix().str();		str_target = m.suffix().str();
	re = regex("-?[0-9]+\\.[0-9]+([eEdD][-+]?[0-9]+)?");		re = regex("-?[0-9]+\\.[0-9]+([eEdD][-+]?[0-9]+)?");
	regex_search(str_target, m, re);		regex_search(str_target, m, re);
	string str_number = m.str();		string str_number = m.str();
	str_number = regex_replace(str_number, regex("D"), "e");		str_number = regex_replace(str_number, regex("D"), "e");
	str_number = regex_replace(str_number, regex("d"), "e");		str_number = regex_replace(str_number, regex("d"), "e");
	ros_vector[i113 - 1] = stod(str_number);		ros_vector[last_configuration - 1] = stod(str_number);
	int nsh = nshl_vector[i113 - 1];		int nsh = nshl_vector[last_configuration - 1];
	if (i113 == 1) nsh += fies;		if (i113 == 1) nsh += fies;
	rcf_vector[i113 - 1] = new double[nsh]();		rcf_vector[last_configuration - 1] = new double[nsh]();
	for (int ns112 = 0; ns112 < nsh; ns112++) {		for (int ns112 = 0; ns112 < nsh; ns112++) {
	str_target = file_lines[++last_read_line];		str_target = file_lines[++last_read_line];
	regex_search(str_target, m, re);		regex_search(str_target, m, re);
	str_number = m.str();		str_number = m.str();
	str_number = regex_replace(str_number, regex("D"), "e");		str_number = regex_replace(str_number, regex("D"), "e");
	str_number = regex_replace(str_number, regex("d"), "e");		str_number = regex_replace(str_number, regex("d"), "e");
	rcf_vector[i113 - 1][ns112] = stod(str_number);		rcf_vector[last_configuration - 1][ns112] = stod(str_number);
	} // ns112 loop		} // ns112 loop
	} // i113 loop		} // i113 loop

	@@ -725,9 +728,11 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(const std::string& de
	} // di loop		} // di loop
	for (int jxi468 = 1; jxi468 <= fnxi; jxi468++) {		for (int jxi468 = 1; jxi468 <= fnxi; jxi468++) {
	if (fidfc != 0 && jxi468 > 1) continue; // jxi468 loop		if (fidfc != 0 && jxi468 > 1) continue; // jxi468 loop
			last_configuration = 0;
	for (int i162 = 1; i162 <= fnsph; i162++) {		for (int i162 = 1; i162 <= fnsph; i162++) {
	if (iog_vector[i162 - 1] >= i162) {		if (iog_vector[i162 - 1] >= i162) {
	int nsh = nshl_vector[i162 - 1];		last_configuration++;
			int nsh = nshl_vector[last_configuration - 1];
	int ici = (nsh + 1) / 2;		int ici = (nsh + 1) / 2;
	if (i162 == 1) ici += fies;		if (i162 == 1) ici += fies;
	for (int ic157 = 0; ic157 < ici; ic157++) {		for (int ic157 = 0; ic157 < ici; ic157++) {
	@@ -781,6 +786,7 @@ ScattererConfiguration* ScattererConfiguration::from_hdf5(const std::string& fil
	string str_name, str_type;		string str_name, str_type;
	if (status == 0) {		if (status == 0) {
	int nsph, ies;		int nsph, ies;
			int last_configuration;
	int *iog;		int *iog;
	double _exdc, _wp, _xip;		double _exdc, _wp, _xip;
	int _idfc, nxi;		int _idfc, nxi;
	@@ -804,22 +810,24 @@ ScattererConfiguration* ScattererConfiguration::from_hdf5(const std::string& fil
	if (iog[ci] < ci + 1) continue;		if (iog[ci] < ci + 1) continue;
	configuration_count++;		configuration_count++;
	}		}
	nshl_vector = new int[nsph]();		nshl_vector = new int[configuration_count]();
	ros_vector = new double[nsph]();		ros_vector = new double[configuration_count]();
	rcf_vector = new double*[nsph];		rcf_vector = new double*[configuration_count];
	for (int i115 = 1; i115 <= nsph; i115++) {		last_configuration = 0;
			for (int i115 = 1; i115 <= configuration_count; i115++) {
	if (iog[i115 - 1] < i115) continue;		if (iog[i115 - 1] < i115) continue;
	str_name = "NSHL_" + to_string(iog[i115 - 1]);		else last_configuration++;
			str_name = "NSHL_" + to_string(last_configuration);
	str_type = "INT32_(1)";		str_type = "INT32_(1)";
	status = hdf_file->read(str_name, str_type, (nshl_vector + i115 - 1));		status = hdf_file->read(str_name, str_type, (nshl_vector + last_configuration - 1));
	str_name = "ROS_" + to_string(iog[i115 - 1]);		str_name = "ROS_" + to_string(last_configuration);
	str_type = "FLOAT64_(1)";		str_type = "FLOAT64_(1)";
	status = hdf_file->read(str_name, str_type, (ros_vector + i115 - 1));		status = hdf_file->read(str_name, str_type, (ros_vector + last_configuration - 1));
	int nsh = nshl_vector[i115 - 1];		int nsh = nshl_vector[last_configuration - 1];
	rcf_vector[i115 - 1] = new double[nsh]();		rcf_vector[last_configuration - 1] = new double[nsh]();
	str_name = "RCF_" + to_string(iog[i115 - 1]);		str_name = "RCF_" + to_string(last_configuration);
	str_type = "FLOAT64_(" + to_string(nsh) + ")";		str_type = "FLOAT64_(" + to_string(nsh) + ")";
	status = hdf_file->read(str_name, str_type, (rcf_vector[i115 - 1]));		status = hdf_file->read(str_name, str_type, (rcf_vector[last_configuration - 1]));
	}		}
	xi_vec = new double[nxi]();		xi_vec = new double[nxi]();
	str_name = "XIVEC";		str_name = "XIVEC";
	@@ -874,6 +882,7 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(const std::string& f
	int _nsph;		int _nsph;
	int *_iog_vec;		int *_iog_vec;
	int _ies;		int _ies;
			int last_configuration;
	double _exdc, _wp, _xip;		double _exdc, _wp, _xip;
	int _idfc, _nxi;		int _idfc, _nxi;
	int *_nshl_vec;		int *_nshl_vec;
	@@ -901,18 +910,20 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(const std::string& f
	for (int i = 1; i <= _nsph; i++) {		for (int i = 1; i <= _nsph; i++) {
	if (_iog_vec[i - 1] >= i) configurations++;		if (_iog_vec[i - 1] >= i) configurations++;
	}		}
	_nshl_vec = new int[_nsph]();		_nshl_vec = new int[configurations]();
	_ros_vec = new double[_nsph]();		_ros_vec = new double[configurations]();
	_rcf_vec = new double*[_nsph];		_rcf_vec = new double*[configurations];
			last_configuration = 0;
	for (int i115 = 1; i115 <= _nsph; i115++) {		for (int i115 = 1; i115 <= _nsph; i115++) {
	if (_iog_vec[i115 - 1] < i115) continue;		if (_iog_vec[i115 - 1] < i115) continue;
	input.read(reinterpret_cast<char *>(&(_nshl_vec[i115 - 1])), sizeof(int));		else last_configuration++;
	input.read(reinterpret_cast<char *>(&(_ros_vec[i115 - 1])), sizeof(double));		input.read(reinterpret_cast<char *>(&(_nshl_vec[last_configuration - 1])), sizeof(int));
	int nsh = _nshl_vec[i115 - 1];		input.read(reinterpret_cast<char *>(&(_ros_vec[last_configuration - 1])), sizeof(double));
			int nsh = _nshl_vec[last_configuration - 1];
	if (i115 == 1) nsh += _ies;		if (i115 == 1) nsh += _ies;
	_rcf_vec[i115 - 1] = new double[nsh]();		_rcf_vec[last_configuration - 1] = new double[nsh]();
	for (int nsi = 0; nsi < nsh; nsi++)		for (int nsi = 0; nsi < nsh; nsi++)
	input.read(reinterpret_cast<char *>(&(_rcf_vec[i115 - 1][nsi])), sizeof(double));		input.read(reinterpret_cast<char *>(&(_rcf_vec[last_configuration - 1][nsi])), sizeof(double));
	}		}
	_dc0m = new dcomplex**[configurations];		_dc0m = new dcomplex**[configurations];
	int dim3 = (_idfc == 0) ? _nxi : 1;		int dim3 = (_idfc == 0) ? _nxi : 1;
	@@ -922,9 +933,11 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(const std::string& f
	} // di loop		} // di loop
	for (int jxi468 = 1; jxi468 <= _nxi; jxi468++) {		for (int jxi468 = 1; jxi468 <= _nxi; jxi468++) {
	if (_idfc != 0 && jxi468 > 1) continue;		if (_idfc != 0 && jxi468 > 1) continue;
			last_configuration = 0;
	for (int i162 = 1; i162 <= _nsph; i162++) {		for (int i162 = 1; i162 <= _nsph; i162++) {
	if (_iog_vec[i162 - 1] < i162) continue;		if (_iog_vec[i162 - 1] < i162) continue;
	int nsh = _nshl_vec[i162 - 1];		else last_configuration++;
			int nsh = _nshl_vec[last_configuration];
	int ici = (nsh + 1) / 2;		int ici = (nsh + 1) / 2;
	if (i162 == 1) ici = ici + _ies;		if (i162 == 1) ici = ici + _ies;
	for (int i157 = 0; i157 < ici; i157++) {		for (int i157 = 0; i157 < ici; i157++) {
	@@ -957,26 +970,6 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(const std::string& f
	return conf;		return conf;
	}		}

	/*
	double ScattererConfiguration::get_param(const std::string& param_name) {
	double value;
	if (param_name.compare("number_of_spheres") == 0) value = (double)number_of_spheres;
	else if (param_name.compare("nsph") == 0) value = (double)number_of_spheres;
	else if (param_name.compare("configurations") == 0) value = (double)configurations;
	else if (param_name.compare("number_of_scales") == 0) value = (double)number_of_scales;
	else if (param_name.compare("nxi") == 0) value = (double)number_of_scales;
	else if (param_name.compare("idfc") == 0) value = (double)idfc;
	else if (param_name.compare("exdc") == 0) value = exdc;
	else if (param_name.compare("wp") == 0) value = wp;
	else if (param_name.compare("xip") == 0) value = xip;
	else {
	string message = "unrecognized parameter \"" + param_name + "\"";
	throw UnrecognizedParameterException(message);
	}
	return value;
	}
	*/

	void ScattererConfiguration::print() {		void ScattererConfiguration::print() {
	int ies = (_use_external_sphere)? 1 : 0;		int ies = (_use_external_sphere)? 1 : 0;
	printf("### CONFIGURATION DATA ###\n");		printf("### CONFIGURATION DATA ###\n");
	@@ -1035,6 +1028,7 @@ void ScattererConfiguration::write_binary(const std::string& file_name, const st

	void ScattererConfiguration::write_hdf5(const std::string& file_name) {		void ScattererConfiguration::write_hdf5(const std::string& file_name) {
	int ies = (use_external_sphere)? 1 : 0;		int ies = (use_external_sphere)? 1 : 0;
			int last_configuration;
	List<string> *rec_name_list = new List<string>(1);		List<string> *rec_name_list = new List<string>(1);
	List<string> *rec_type_list = new List<string>(1);		List<string> *rec_type_list = new List<string>(1);
	List<void > rec_ptr_list = new List<void *>(1);		List<void > rec_ptr_list = new List<void *>(1);
	@@ -1068,23 +1062,25 @@ void ScattererConfiguration::write_hdf5(const std::string& file_name) {
	str_type = "FLOAT64_(" + to_string(number_of_scales) + ")";		str_type = "FLOAT64_(" + to_string(number_of_scales) + ")";
	rec_type_list->append(str_type);		rec_type_list->append(str_type);
	rec_ptr_list->append(_scale_vec);		rec_ptr_list->append(_scale_vec);
			last_configuration = 0;
	for (int i115 = 1; i115 <= _number_of_spheres; i115++) {		for (int i115 = 1; i115 <= _number_of_spheres; i115++) {
	if (_iog_vec[i115 - 1] < i115) continue;		if (_iog_vec[i115 - 1] < i115) continue;
	str_name = "NSHL_" + to_string(i115);		else last_configuration++;
			str_name = "NSHL_" + to_string(last_configuration);
	rec_name_list->append(str_name);		rec_name_list->append(str_name);
	rec_type_list->append("INT32_(1)");		rec_type_list->append("INT32_(1)");
	rec_ptr_list->append(&(_nshl_vec[i115 - 1])); // was not from IOG		rec_ptr_list->append(&(_nshl_vec[last_configuration - 1])); // was not from IOG
	str_name = "ROS_" + to_string(i115);		str_name = "ROS_" + to_string(last_configuration);
	rec_name_list->append(str_name);		rec_name_list->append(str_name);
	rec_type_list->append("FLOAT64_(1)");		rec_type_list->append("FLOAT64_(1)");
	rec_ptr_list->append(&(_radii_of_spheres[i115 - 1])); // was not from IOG		rec_ptr_list->append(&(_radii_of_spheres[last_configuration - 1])); // was not from IOG
	int nsh = _nshl_vec[i115 - 1]; // was not from IOG		int nsh = _nshl_vec[last_configuration - 1]; // was not from IOG
	if (i115 == 1) nsh += ies;		if (i115 == 1) nsh += ies;
	str_name = "RCF_" + to_string(i115); // was not from IOG		str_name = "RCF_" + to_string(last_configuration); // was not from IOG
	str_type = "FLOAT64_(" + to_string(nsh) + ")";		str_type = "FLOAT64_(" + to_string(nsh) + ")";
	rec_name_list->append(str_name);		rec_name_list->append(str_name);
	rec_type_list->append(str_type);		rec_type_list->append(str_type);
	rec_ptr_list->append(&(_rcf[i115 - 1][0])); // was not from IOG		rec_ptr_list->append(&(_rcf[last_configuration - 1][0])); // was not from IOG
	}		}

	int dim3 = (idfc == 0) ? _number_of_scales : 1;		int dim3 = (idfc == 0) ? _number_of_scales : 1;
	@@ -1093,9 +1089,11 @@ void ScattererConfiguration::write_hdf5(const std::string& file_name) {
	int dc0_index = 0;		int dc0_index = 0;
	for (int jxi468 = 1; jxi468 <= _number_of_scales; jxi468++) {		for (int jxi468 = 1; jxi468 <= _number_of_scales; jxi468++) {
	if (_idfc != 0 && jxi468 > 1) continue;		if (_idfc != 0 && jxi468 > 1) continue;
			last_configuration = 0;
	for (int i162 = 1; i162 <= _number_of_spheres; i162++) {		for (int i162 = 1; i162 <= _number_of_spheres; i162++) {
	if (_iog_vec[i162 - 1] < i162) continue;		if (_iog_vec[i162 - 1] < i162) continue;
	int nsh = _nshl_vec[i162 - 1]; // was not from IOG		else last_configuration++;
			int nsh = _nshl_vec[last_configuration - 1]; // was not from IOG
	int ici = (nsh + 1) / 2;		int ici = (nsh + 1) / 2;
	if (i162 == 1) ici = ici + ies;		if (i162 == 1) ici = ici + ies;
	for (int i157 = 0; i157 < ici; i157++) {		for (int i157 = 0; i157 < ici; i157++) {
	@@ -1138,6 +1136,7 @@ void ScattererConfiguration::write_hdf5(const std::string& file_name) {
	void ScattererConfiguration::write_legacy(const std::string& file_name) {		void ScattererConfiguration::write_legacy(const std::string& file_name) {
	fstream output;		fstream output;
	int ies = (_use_external_sphere)? 1 : 0;		int ies = (_use_external_sphere)? 1 : 0;
			int last_configuration;
	output.open(file_name.c_str(), ios::out \| ios::binary);		output.open(file_name.c_str(), ios::out \| ios::binary);
	output.write(reinterpret_cast<char *>(&_number_of_spheres), sizeof(int));		output.write(reinterpret_cast<char *>(&_number_of_spheres), sizeof(int));
	output.write(reinterpret_cast<char *>(&ies), sizeof(int));		output.write(reinterpret_cast<char *>(&ies), sizeof(int));
	@@ -1150,20 +1149,24 @@ void ScattererConfiguration::write_legacy(const std::string& file_name) {
	output.write(reinterpret_cast<char *>(&_number_of_scales), sizeof(int));		output.write(reinterpret_cast<char *>(&_number_of_scales), sizeof(int));
	for (int i = 0; i < _number_of_scales; i++)		for (int i = 0; i < _number_of_scales; i++)
	output.write(reinterpret_cast<char *>(&(_scale_vec[i])), sizeof(double));		output.write(reinterpret_cast<char *>(&(_scale_vec[i])), sizeof(double));
			last_configuration = 0;
	for (int i115 = 1; i115 <= _number_of_spheres; i115++) {		for (int i115 = 1; i115 <= _number_of_spheres; i115++) {
	if (_iog_vec[i115 - 1] < i115) continue;		if (_iog_vec[i115 - 1] < i115) continue;
	output.write(reinterpret_cast<char *>(&(_nshl_vec[i115 - 1])), sizeof(int)); // was not from IOG		else last_configuration++;
	output.write(reinterpret_cast<char *>(&(_radii_of_spheres[i115 - 1])), sizeof(double)); // was not from IOG		output.write(reinterpret_cast<char *>(&(_nshl_vec[last_configuration - 1])), sizeof(int)); // was not from IOG
	int nsh = _nshl_vec[i115 - 1]; // was not from IOG		output.write(reinterpret_cast<char *>(&(_radii_of_spheres[last_configuration - 1])), sizeof(double)); // was not from IOG
			int nsh = _nshl_vec[last_configuration - 1]; // was not from IOG
	if (i115 == 1) nsh += ies;		if (i115 == 1) nsh += ies;
	for (int nsi = 0; nsi < nsh; nsi++)		for (int nsi = 0; nsi < nsh; nsi++)
	output.write(reinterpret_cast<char *>(&(_rcf[i115 - 1][nsi])), sizeof(double)); // was not from IOG		output.write(reinterpret_cast<char *>(&(_rcf[last_configuration - 1][nsi])), sizeof(double)); // was not from IOG
	}		}
	for (int jxi468 = 1; jxi468 <= _number_of_scales; jxi468++) {		for (int jxi468 = 1; jxi468 <= _number_of_scales; jxi468++) {
	if (idfc != 0 && jxi468 > 1) continue;		if (idfc != 0 && jxi468 > 1) continue;
			last_configuration = 0;
	for (int i162 = 1; i162 <= _number_of_spheres; i162++) {		for (int i162 = 1; i162 <= _number_of_spheres; i162++) {
	if (_iog_vec[i162 - 1] < i162) continue;		if (_iog_vec[i162 - 1] < i162) continue;
	int nsh = _nshl_vec[i162 - 1]; // was not from IOG		else last_configuration++;
			int nsh = _nshl_vec[last_configuration - 1]; // was not from IOG
	int ici = (nsh + 1) / 2; // QUESTION: is integer division really intended here?		int ici = (nsh + 1) / 2; // QUESTION: is integer division really intended here?
	if (i162 == 1) ici = ici + ies;		if (i162 == 1) ici = ici + ies;
	for (int i157 = 0; i157 < ici; i157++) {		for (int i157 = 0; i157 < ici; i157++) {
	@@ -1186,6 +1189,7 @@ void ScattererConfiguration::write_formatted(const std::string& file_name) {
	const double two_pi = acos(0.0) * 4.0;		const double two_pi = acos(0.0) * 4.0;
	double *xi_vec;		double *xi_vec;
	int ici;		int ici;
			int last_configuration;
	int ies = (use_external_sphere)? 1: 0;		int ies = (use_external_sphere)? 1: 0;
	FILE *output = fopen(file_name.c_str(), "w");		FILE *output = fopen(file_name.c_str(), "w");
	int scale_type = -1;		int scale_type = -1;
	@@ -1332,9 +1336,11 @@ void ScattererConfiguration::write_formatted(const std::string& file_name) {
	}		}
	if (_idfc != 0) {		if (_idfc != 0) {
	fprintf(output, " DIELECTRIC CONSTANTS\n");		fprintf(output, " DIELECTRIC CONSTANTS\n");
			last_configuration = 0;
	for (int i473 = 1; i473 <= _number_of_spheres; i473++) {		for (int i473 = 1; i473 <= _number_of_spheres; i473++) {
	if (_iog_vec[i473 - 1] != i473) continue;		if (_iog_vec[i473 - 1] != i473) continue;
	ici = (_nshl_vec[i473 - 1] + 1) / 2;		else last_configuration++;
			ici = (_nshl_vec[last_configuration - 1] + 1) / 2;
	if (i473 == 1) ici += ies;		if (i473 == 1) ici += ies;
	fprintf(output, " SPHERE N.%4d\n", i473);		fprintf(output, " SPHERE N.%4d\n", i473);
	for (int ic472 = 0; ic472 < ici; ic472++) {		for (int ic472 = 0; ic472 < ici; ic472++) {
	@@ -1345,9 +1351,11 @@ void ScattererConfiguration::write_formatted(const std::string& file_name) {
	}		}
	} else {		} else {
	fprintf(output, " DIELECTRIC FUNCTIONS\n");		fprintf(output, " DIELECTRIC FUNCTIONS\n");
			last_configuration = 0;
	for (int i478 = 1; i478 <= _number_of_spheres; i478++) {		for (int i478 = 1; i478 <= _number_of_spheres; i478++) {
	if (_iog_vec[i478 - 1] != i478) continue;		if (_iog_vec[i478 - 1] != i478) continue;
	ici = (_nshl_vec[i478 - 1] + 1) / 2;		last_configuration++;
			ici = (_nshl_vec[last_configuration - 1] + 1) / 2;
	if (i478 == 1) ici += ies;		if (i478 == 1) ici += ies;
	fprintf(output, " SPHERE N.%4d\n", i478);		fprintf(output, " SPHERE N.%4d\n", i478);
	for (int ic477 = 1; ic477 <= ici; ic477++) {		for (int ic477 = 1; ic477 <= ici; ic477++) {
	@@ -1418,4 +1426,3 @@ bool ScattererConfiguration::operator ==(const ScattererConfiguration &other) {
	} // ri loop		} // ri loop
	return true;		return true;
	}		}