Use classes for dynamic management of common blocks

/*! \file Commons.cpp

 *

 */

#include "include/Commons.h"

using namespace std;

C1::C1(int ns, int l_max) {

	nlmmt = 2 * (l_max * (l_max + 2));

	nsph = ns;

	lm = l_max;

	rmi = new complex<double>*[lm];

	rei = new complex<double>*[lm];

	for (int ri = 0; ri < lm; ri++) {

		rmi[ri] = new complex<double>[nsph];

		rei[ri] = new complex<double>[nsph];

	}

	w = new complex<double>*[nlmmt];

	for (int wi = 0; wi < nlmmt; wi++) w[wi] = new complex<double>[4];

	vints = new complex<double>*[nsph];

	rc = new double*[nsph];

	for (int vi = 0; vi < nsph; vi++) {

		rc[vi] = new double[lm];

		vints[vi] = new complex<double>[16];

	}

	fsas = new complex<double>[nsph];

	sscs = new double[nsph];

	sexs = new double[nsph];

	sabs = new double[nsph];

	sqscs = new double[nsph];

	sqexs = new double[nsph];

	sqabs = new double[nsph];

	gcsv = new double[nsph];;

	rxx = new double[nsph];

	ryy = new double[nsph];

	rzz = new double[nsph];

	ros = new double[nsph];

	iog = new int[nsph];

	nshl = new int[nsph];

	sas = new complex<double>**[nsph];

	for (int si = 0; si < nsph; si++) {

		sas[si] = new complex<double>*[2];

		sas[si][0] = new complex<double>[2];

		sas[si][1] = new complex<double>[2];

	}

}

C1::~C1() {

	delete[] rmi;

	delete[] rei;

	for (int wi = 1; wi <= nlmmt; wi++) delete[] w[wi];

	for (int vi = 1; vi <= nsph; vi++) {

		delete[] rc[nsph - vi];

		delete[] vints[nsph - vi];

	}

	for (int si = 1; si <= nsph; si++) {

		delete[] sas[nsph - si][1];

		delete[] sas[nsph - si][0];

	}

	delete[] fsas;

	delete[] sscs;

	delete[] sexs;

	delete[] sabs;

	delete[] sqscs;

	delete[] sqexs;

	delete[] sqabs;

	delete[] gcsv;

	delete[] rxx;

	delete[] ryy;

	delete[] rzz;

	delete[] ros;

	delete[] iog;

	delete[] nshl;

}

C2::C2(int ns, int nl, int npnt, int npntts) {

	nsph = ns;

	int max_n = (npnt > npntts) ? npnt : npntts;

	nhspo = 2 * max_n - 1;

	ris = new complex<double>[nhspo];

	dlri = new complex<double>[nhspo];

	vkt = new complex<double>[nsph];

	dc0 = new complex<double>[nl];

	vsz = new double[nsph];

}

C2::~C2() {

	delete[] ris;

	delete[] dlri;

	delete[] vkt;

	delete[] dc0;

	delete[] vsz;

}

$(SUBDIRS):

	$(MAKE) -C $@

np_sphere: $(BUILLDDIR)/sphere/np_sphere.o $(BUILLDDIR)/sphere/Configuration.o $(BUILLDDIR)/sphere/Parsers.o $(BUILLDDIR)/sphere/sphere.o

	$(CC) -o $(BUILDDIR)/sphere/np_sphere $(BUILDDIR)/sphere/np_sphere.o $(BUILDDIR)/sphere/Configuration.o $(BUILDDIR)/sphere/Parsers.o $(BUILDDIR)/sphere/sphere.o

np_sphere: $(BUILDDIR)/sphere/np_sphere.o $(BUILDDIR)/sphere/Commons.o $(BUILDDIR)/sphere/Configuration.o $(BUILDDIR)/sphere/Parsers.o $(BUILDDIR)/sphere/sphere.o

	$(CC) -o $(BUILDDIR)/sphere/np_sphere $(BUILDDIR)/sphere/np_sphere.o $(BUILDDIR)/sphere/Commons.o $(BUILDDIR)/sphere/Configuration.o $(BUILDDIR)/sphere/Parsers.o $(BUILDDIR)/sphere/sphere.o

clean:

	rm -f $(BUILDDIR)/cluster/*.o

.PHONY: all $(SUBDIRS)

$(BUILLDDIR)/sphere/np_sphere.o:

$(BUILDDIR)/sphere/np_sphere.o:

	$(CC) -c np_sphere.cpp -o $(BUILDDIR)/sphere/np_sphere.o

$(BUILLDDIR)/sphere/Configuration.o:

$(BUILDDIR)/sphere/Commons.o:

	$(CC) -c Commons.cpp -o $(BUILDDIR)/sphere/Commons.o

$(BUILDDIR)/sphere/Configuration.o:

	$(CC) -c Configuration.cpp -o $(BUILDDIR)/sphere/Configuration.o

$(BUILLDDIR)/sphere/Parsers.o:

$(BUILDDIR)/sphere/Parsers.o:

	$(CC) -c Parsers.cpp -o $(BUILDDIR)/sphere/Parsers.o

$(BUILLDDIR)/sphere/sphere.o:

$(BUILDDIR)/sphere/sphere.o:

	$(CC) -c sphere/sphere.cpp -o $(BUILDDIR)/sphere/sphere.o

/*! \file Commons.h

 *

 * \brief C++ porting of common data structures.

 *

 * Many functions of the original FORTRAN code share complex data blocks in

 * form of COMMON blocks. This poses the limit of freezing the structure of

 * the data blocks in the code, therefore implying the necessity to modify

 * the code to adapt it to the input and to recompile before running. C++,

 * on the contrary, offers the possibility to represent the necessary data

 * structures as classes that can dinamically instantiate the shared information

 * in the most convenient format for the current configuration. This approach

 * adds an abstraction layer that lifts the need to modify and recompile the

 * code depending on the input.

 *

 */

#ifndef SRC_INCLUDE_COMMONS_

#define SRC_INCLUDE_COMMONS_

#include <complex>

/*! \brief Representation of the FORTRAN C1 common blocks.

 *

 * C1 common blocks are used to store vector field expansions and geometric

 * properties, such as sphere sizes, positions and cross-sections. These are

 * used by functions such as `aps`, `diel`, `pwma`, `rabas`, `sscr0`, `sscr2`,

 * `wmamp`, `wmasp` and `dme`. QUESTION: correct?

 *

 * Due to the necessity to share the class contents with many functions that

 * need to read and update various fields, all shared members of C1 are declared

 * public (i.e., the class is just a structure with more advanced constructor

 * and destroyer). Further development may go in the direction of creating

 * a better encapsulation, either by unpacking the contents (recommended) or

 * by creating public methods to access protected fields (standard way, but not

 * recommended for HPC applications).

 */

class C1 {

protected:

	//! \brief Number of spheres.

	int nsph;

	//! \brief Maximum order of field expansion.

	int lm;

	//! \brief NLMMT. QUESTION: definition?

	int nlmmt;

public:

	std::complex<double> **rmi;

	std::complex<double> **rei;

	std::complex<double> **w;

	std::complex<double> *fsas;

	std::complex<double> **vints;

	double *sscs;

	double *sexs;

	double *sabs;

	double *sqscs;

	double *sqexs;

	double *sqabs;

	double *gcsv;

	double *rxx;

	double *ryy;

	double *rzz;

	double *ros;

	double **rc;

	int *iog;

	int *nshl;

	std::complex<double> ***sas;

	/*! \brief C1 instance constructor.

	 *

	 * \param ns: `int` Number of spheres.

	 * \param l_max: `int` Maximum order of field expansion.

	 */

	C1(int ns, int l_max);

	//! \brief C1 instance destroyer.

	~C1();

};

/*! \brief Representation of the FORTRAN C2 blocks.

 *

 */

class C2 {

	int nsph, nhspo;

public:

	std::complex<double> *ris;

	std::complex<double> *dlri;

	std::complex<double> *dc0;

	std::complex<double> *vkt;

	double *vsz;

	/*! \brief C2 instance constructor.

	 *

	 * \param ns: `int` Number of spheres.

	 * \param nl: `int`

	 * \param npnt: `int`

	 * \param npntts: `int`

	 */

	C2(int ns, int nl, int npnt, int npntts);

	//! \brief C2 instance destroyer.

	~C2();

};

#endif

#define SRC_INCLUDE_SPH_SUBS_H_

#include <complex>

//! \brief A structure representing the common block C1

struct C1 {

	std::complex<double> rmi[40][2];

	std::complex<double> rei[40][2];

	std::complex<double> w[3360][4];

	std::complex<double> fsas[2];

	std::complex<double> sas[2][2][2];

	std::complex<double> vints[2][16];

	double sscs[2];

	double sexs[2];

	double sabs[2];

	double sqscs[2];

	double sqexs[2];

	double sqabs[2];

	double gcsv[2];

	double rxx[1];

	double ryy[1];

	double rzz[1];

	double ros[2];

	double rc[2][8];

	int iog[2];

	int nshl[2];

};

//! \brief A structure representing the common block C1

struct C2 {

	std::complex<double> ris[999];

	std::complex<double> dlri[999];

	std::complex<double> dc0[5];

	std::complex<double> vkt[2];

	double vsz[2];

};

#include "Commons.h"

/*! \brief Conjugate of a double precision complex number

 *

		complex<double> *vint = new complex<double>[16];

		int jw;

		int nsph = gconf->number_of_spheres;

		C1 *c1 = new C1;

		C1 *c1 = new C1(nsph, gconf->l_max);

		for (int i = 0; i < nsph; i++) {

			c1->iog[i] = sconf->iog_vec[i];

			c1->nshl[i] = sconf->nshl_vec[i];

			c1->rei[i][0] = complex<double>(0.0, 0.0);

			c1->rei[i][1] = complex<double>(0.0, 0.0);

		}

		C2 *c2 = new C2;

		C2 *c2 = new C2(nsph, 5, gconf->npnt, gconf->npntts);

		argi = new double[1];

		args = new double[1];

		gaps = new double[2];