From cb809b8225fc8090fd885722a4aca07e88376b7b Mon Sep 17 00:00:00 2001
From: "Mulas, Giacomo" <gmulas@oa-cagliari.inaf.it>
Date: Fri, 15 Mar 2024 10:01:11 +0100
Subject: [PATCH] - add constructors to the "Common" objects allowing them to
 be copied - turn C4 from a structure to a class (to do the above), move its
 initialisation to its constructor - eliminate "allocate_vectors()", redundant
 and including a (wrong) double allocation

---
 src/cluster/cluster.cpp |  15 +-
 src/include/Commons.h   | 519 ++++++++++++++++++++++------------------
 src/libnptm/Commons.cpp | 317 ++++++++++++++++++++++--
 3 files changed, 584 insertions(+), 267 deletions(-)

diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index 7d41b25e..ce5de2c0 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -107,20 +107,7 @@ void cluster(string config_file, string data_file, string output_path) {
       c1->rzz[c1i] = gconf->sph_z[c1i];
       c1->ros[c1i] = sconf->radii_of_spheres[c1i];
     }
-    C4 *c4 = new C4;
-    c4->li = gconf->li;
-    c4->le = gconf->le;
-    c4->lm = (c4->li > c4-> le) ? c4->li : c4->le;
-    c4->nv3j = (c4->lm * (c4->lm  + 1) * (2 * c4->lm + 7)) / 6;
-    c4->nsph = nsph;
-    // The following is needed to initialize C1_AddOns
-    c4->litpo = c4->li + c4->li + 1;
-    c4->litpos = c4->litpo * c4->litpo;
-    c4->lmtpo = c4->li + c4->le + 1;
-    c4->lmtpos = c4->lmtpo * c4->lmtpo;
-    c4->nlim = c4->li * (c4->li + 2);
-    c4->nlem = c4->le * (c4->le + 2);
-    c4->lmpo = c4->lm + 1;
+    C4 *c4 = new C4(gconf->li, gconf->le, nsph);
     C1_AddOns *c1ao = new C1_AddOns(c4);
     // End of add-ons initialization
     C6 *c6 = new C6(c4->lmtpo);
diff --git a/src/include/Commons.h b/src/include/Commons.h
index e9545444..77588510 100644
--- a/src/include/Commons.h
+++ b/src/include/Commons.h
@@ -36,153 +36,181 @@
  */
 class C1 {
 protected:
-	//! \brief Number of spheres.
-	int nsph;
-	//! \brief Maximum order of field expansion.
-	int lm;
-	//! \brief NLMMT. QUESTION: definition?
-	int nlmmt;
+  //! \brief Number of spheres.
+  int nsph;
+  //! \brief Maximum order of field expansion.
+  int lm;
+  //! \brief NLMMT. QUESTION: definition?
+  int nlmmt;
+  //! \brief Number of configurations
+  int configurations;
 public:
-	//! \brief QUESTION: definition?
-	dcomplex **rmi;
-	//! \brief QUESTION: definition?
-	dcomplex **rei;
-	//! \brief QUESTION: definition?
-	dcomplex **w;
-	//! \brief QUESTION: definition?
-	dcomplex *fsas;
-	//! \brief QUESTION: definition?
-	dcomplex **vints;
-	//! \brief QUESTION: definition?
-	double *sscs;
-	//! \brief QUESTION: definition?
-	double *sexs;
-	//! \brief QUESTION: definition?
-	double *sabs;
-	//! \brief QUESTION: definition?
-	double *sqscs;
-	//! \brief QUESTION: definition?
-	double *sqexs;
-	//! \brief QUESTION: definition?
-	double *sqabs;
-	//! \brief QUESTION: definition?
-	double *gcsv;
-	//! \brief Vector of sphere X coordinates.
-	double *rxx;
-	//! \brief Vector of sphere X coordinates.
-	double *ryy;
-	//! \brief Vector of sphere X coordinates.
-	double *rzz;
-	//! \brief Vector of sphere radii.
-	double *ros;
-	//! \brief Matrix of spherical layer break radii. QUESTION: correct?
-	double **rc;
-	//! \brief Vector of spherical component identifiers.
-	int *iog;
-	//! \brief Vector of numbers of spherical layers.
-	int *nshl;
-	//! \brief QUESTION: definition?
-	dcomplex ***sas;
+  //! \brief QUESTION: definition?
+  dcomplex **rmi;
+  //! \brief QUESTION: definition?
+  dcomplex **rei;
+  //! \brief QUESTION: definition?
+  dcomplex **w;
+  //! \brief QUESTION: definition?
+  dcomplex *fsas;
+  //! \brief QUESTION: definition?
+  dcomplex **vints;
+  //! \brief QUESTION: definition?
+  double *sscs;
+  //! \brief QUESTION: definition?
+  double *sexs;
+  //! \brief QUESTION: definition?
+  double *sabs;
+  //! \brief QUESTION: definition?
+  double *sqscs;
+  //! \brief QUESTION: definition?
+  double *sqexs;
+  //! \brief QUESTION: definition?
+  double *sqabs;
+  //! \brief QUESTION: definition?
+  double *gcsv;
+  //! \brief Vector of sphere X coordinates.
+  double *rxx;
+  //! \brief Vector of sphere X coordinates.
+  double *ryy;
+  //! \brief Vector of sphere X coordinates.
+  double *rzz;
+  //! \brief Vector of sphere radii.
+  double *ros;
+  //! \brief Matrix of spherical layer break radii. QUESTION: correct?
+  double **rc;
+  //! \brief Vector of spherical component identifiers.
+  int *iog;
+  //! \brief Vector of numbers of spherical layers.
+  int *nshl;
+  //! \brief QUESTION: definition?
+  dcomplex ***sas;
 
-	/*! \brief C1 instance constructor.
-	 *
-	 * \param ns: `int` Number of spheres.
-	 * \param l_max: `int` Maximum order of field expansion.
-	 * \param nshl: `int *` Array of number of layers in spheres.
-	 * \param iog: `int *` Vector of spherical units ID numbers.
-	 */
-	C1(int ns, int l_max, int *nshl, int *iog);
-
-	//! \brief C1 instance destroyer.
-	~C1();
+  /*! \brief C1 instance constructor.
+   *
+   * \param ns: `int` Number of spheres.
+   * \param l_max: `int` Maximum order of field expansion.
+   * \param nshl: `int *` Array of number of layers in spheres.
+   * \param iog: `int *` Vector of spherical units ID numbers.
+   */
+  C1(int ns, int l_max, int *nshl, int *iog);
+  /*! \brief C1 instance constructor copying all contents from a preexisting template
+   *
+   * \param rhs: `C1` preexisting template.
+   */
+  C1(const C1& rhs);
+  //! \brief C1 instance destroyer.
+  ~C1();
 };
 
 /*! \brief Representation of the FORTRAN C2 blocks.
  *
  */
 class C2 {
-	//! \brief Number of spheres.
-	int nsph;
-	//! \brief Number of required orders.
-	int nhspo;
+  //! \brief Number of spheres.
+  int nsph;
+  //! \brief Number of required orders.
+  int nhspo;
+  //! \brief QUESTION: what is nl?
+  int nl;
 public:
-	//! \brief QUESTION: definition?
-	dcomplex *ris;
-	//! \brief QUESTION: definition?
-	dcomplex *dlri;
-	//! \brief QUESTION: definition?
-	dcomplex *dc0;
-	//! \brief QUESTION: definition?
-	dcomplex *vkt;
-	//! Vector of scaled sizes. QUESTION: correct?
-	double *vsz;
+  //! \brief QUESTION: definition?
+  dcomplex *ris;
+  //! \brief QUESTION: definition?
+  dcomplex *dlri;
+  //! \brief QUESTION: definition?
+  dcomplex *dc0;
+  //! \brief QUESTION: definition?
+  dcomplex *vkt;
+  //! Vector of scaled sizes. QUESTION: correct?
+  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 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 constructor copying its contents from preexisting instance.
+   *
+   * \param rhs: `C2` object to copy contents from
+   */
+  C2(const C2& rhs);
 
-	//! \brief C2 instance destroyer.
-	~C2();
+  //! \brief C2 instance destroyer.
+  ~C2();
 };
 
 /*! \brief Representation of the FORTRAN C3 blocks.
  */
 class C3 {
 public:
-	//! \brief QUESTION: definition?
-	dcomplex tfsas;
-	//! \brief QUESTION: definition?
-	dcomplex **tsas;
-	//! \brief QUESTION: definition?
-	double gcs;
-	//! \brief QUESTION: definition?
-	double scs;
-	//! \brief QUESTION: definition?
-	double ecs;
-	//! \brief QUESTION: definition?
-	double acs;
+  //! \brief QUESTION: definition?
+  dcomplex tfsas;
+  //! \brief QUESTION: definition?
+  dcomplex **tsas;
+  //! \brief QUESTION: definition?
+  double gcs;
+  //! \brief QUESTION: definition?
+  double scs;
+  //! \brief QUESTION: definition?
+  double ecs;
+  //! \brief QUESTION: definition?
+  double acs;
 
-	/*! \brief C3 instance constructor.
-	 */
-	C3();
+  /*! \brief C3 instance constructor.
+   */
+  C3();
+  /*! \brief C3 instance constructor copying its contents from a preexisting object.
+   */
+  C3(const C3& rhs);
 
-	/*! \brief C3 instance destroyer.
-	 */
-	~C3();
+  /*! \brief C3 instance destroyer.
+   */
+  ~C3();
 };
 
 /*! \brief Representation of the FORTRAN C4 blocks.
  */
-struct C4 {
-	//! \brief QUESTION: definition?
-	int litpo;
-	//! \brief QUESTION: definition?
-	int litpos;
-	//! \brief Maximum field expansion order plus one. QUESTION: correct?
-	int lmpo;
-	//! \brief Twice maximum field expansion order plus one. QUESTION: correct?
-	int lmtpo;
-	//! \brief Square of `lmtpo`.
-	int lmtpos;
-	//! \brief QUESTION: definition?
-	int li;
-	//! \brief QUESTION: definition?
-	int nlim;
-	//! \brief QUESTION: definition?
-	int le;
-	//! \brief QUESTION: definition?
-	int nlem;
-	//! \brief Maximum field expansion order. QUESTION: correct?
-	int lm;
-	//! \brief Number of spheres.
-	int nsph;
-	//! \brief QUESTION: definition?
-	int nv3j;
+class C4 {
+public:
+  //! \brief QUESTION: definition?
+  int litpo;
+  //! \brief QUESTION: definition?
+  int litpos;
+  //! \brief Maximum field expansion order plus one. QUESTION: correct?
+  int lmpo;
+  //! \brief Twice maximum field expansion order plus one. QUESTION: correct?
+  int lmtpo;
+  //! \brief Square of `lmtpo`.
+  int lmtpos;
+  //! \brief QUESTION: definition?
+  int li;
+  //! \brief QUESTION: definition?
+  int nlim;
+  //! \brief QUESTION: definition?
+  int le;
+  //! \brief QUESTION: definition?
+  int nlem;
+  //! \brief Maximum field expansion order. QUESTION: correct?
+  int lm;
+  //! \brief Number of spheres.
+  int nsph;
+  //! \brief QUESTION: definition?
+  int nv3j;
+
+  /*! \brief C3 instance constructor.
+   */
+  C4(int li, int le, int nsph);
+  /*! \brief C3 instance constructor copying its contents from a preexisting object.
+   */
+  C4(const C4& rhs);
+
+  /*! \brief C3 instance destroyer.
+   */
+  ~C4();
 };
 
 /*! \brief Vectors and matrices that are specific to cluster C1 blocks.
@@ -190,131 +218,164 @@ struct C4 {
  */
 class C1_AddOns {
 protected:
-	//! \brief Number of spheres.
-	int nsph;
-	//! \brief QUESTION: definition?
-	int nlemt;
-	//! \brief Maximum expansion order plus one. QUESTION: correct?
-	int lmpo;
+  //! \brief Number of spheres.
+  int nsph;
+  //! \brief QUESTION: definition?
+  int nlemt;
+  //! \brief Maximum expansion order plus one. QUESTION: correct?
+  int lmpo;
+  //! \brief QUESTION: definition?
+  int litpo;
+  //! \brief QUESTION: definition?
+  int lmtpo;
+  //! \brief QUESTION: definition?
+  int litpos;
+  //! \brief QUESTION: definition?
+  int lmtpos;
+  //! \brief QUESTION: definition?
+  int nv3j;
+  //! \brief QUESTION: definition?
+  int lm;
 
-	/*! \brief Allocate the necessary common vectors depending on configuration.
-	 *
-	 * The size of the vectors and matrices defined in various common
-	 * blocks, and particularly in C1, depends on many settings of the
-	 * problem configuration, such as the number of spheres, the number
-	 * of layers the spheres are made of, the field expansion order and
-	 * others. This function collects the calculations needed to infer
-	 * the necessary amount of memory for these configurable elements,
-	 * thus making the class constructor more compact and easier to handle.
-	 *
-	 * \param c4: `C4 *` Pointer to a C4 structure.
-	 */
-	void allocate_vectors(C4 *c4);
+  /*! \brief Allocate the necessary common vectors depending on configuration.
+   *
+   * The size of the vectors and matrices defined in various common
+   * blocks, and particularly in C1, depends on many settings of the
+   * problem configuration, such as the number of spheres, the number
+   * of layers the spheres are made of, the field expansion order and
+   * others. This function collects the calculations needed to infer
+   * the necessary amount of memory for these configurable elements,
+   * thus making the class constructor more compact and easier to handle.
+   *
+   * \param c4: `C4 *` Pointer to a C4 structure.
+   */
+  //void allocate_vectors(C4 *c4);
 public:
-	//! \brief QUESTION: definition?
-	dcomplex *vh;
-	//! \brief QUESTION: definition?
-	dcomplex *vj0;
-	//! \brief QUESTION: definition?
-	dcomplex *vj;
-	//! \brief QUESTION: definition?
-	dcomplex *vyhj;
-	//! \brief QUESTION: definition?
-	dcomplex *vyj0;
-	//! \brief QUESTION: definition?
-	dcomplex **am0m;
-	//! \brief QUESTION: definition?
-	dcomplex *vint;
-	//! \brief QUESTION: definition?
-	dcomplex *vintm;
-	//! \brief QUESTION: definition?
-	dcomplex **vints;
-	//! \brief QUESTION: definition?
-	dcomplex *vintt;
-	//! \brief QUESTION: definition?
-	dcomplex **fsac;
-	//! \brief QUESTION: definition?
-	dcomplex **sac;
-	//! \brief QUESTION: definition?
-	dcomplex **fsacm;
-	//! \brief QUESTION: definition?
-	double *scsc;
-	//! \brief QUESTION: definition?
-	dcomplex *scscp;
-	//! \brief QUESTION: definition?
-	double *ecsc;
-	//! \brief QUESTION: definition?
-	double *ecscm;
-	//! \brief QUESTION: definition?
-	double *scscm;
-	//! \brief QUESTION: definition?
-	dcomplex *ecscp;
-	//! \brief QUESTION: definition?
-	dcomplex *scscpm;
-	//! \brief QUESTION: definition?
-	dcomplex *ecscpm;
-	//! \brief QUESTION: definition?
-	double *v3j0;
-	//! \brief QUESTION: definition?
-	double *sscs;
-	//! \brief QUESTION: definition?
-	int **ind3j;
+  //! \brief QUESTION: definition?
+  dcomplex *vh;
+  //! \brief QUESTION: definition?
+  dcomplex *vj0;
+  //! \brief QUESTION: definition?
+  dcomplex *vj;
+  //! \brief QUESTION: definition?
+  dcomplex *vyhj;
+  //! \brief QUESTION: definition?
+  dcomplex *vyj0;
+  //! \brief QUESTION: definition?
+  dcomplex **am0m;
+  //! \brief QUESTION: definition?
+  dcomplex *vint;
+  //! \brief QUESTION: definition?
+  dcomplex *vintm;
+  //! \brief QUESTION: definition?
+  dcomplex **vints;
+  //! \brief QUESTION: definition?
+  dcomplex *vintt;
+  //! \brief QUESTION: definition?
+  dcomplex **fsac;
+  //! \brief QUESTION: definition?
+  dcomplex **sac;
+  //! \brief QUESTION: definition?
+  dcomplex **fsacm;
+  //! \brief QUESTION: definition?
+  double *scsc;
+  //! \brief QUESTION: definition?
+  dcomplex *scscp;
+  //! \brief QUESTION: definition?
+  double *ecsc;
+  //! \brief QUESTION: definition?
+  double *ecscm;
+  //! \brief QUESTION: definition?
+  double *scscm;
+  //! \brief QUESTION: definition?
+  dcomplex *ecscp;
+  //! \brief QUESTION: definition?
+  dcomplex *scscpm;
+  //! \brief QUESTION: definition?
+  dcomplex *ecscpm;
+  //! \brief QUESTION: definition?
+  double *v3j0;
+  //! \brief QUESTION: definition?
+  double *sscs;
+  //! \brief QUESTION: definition?
+  int **ind3j;
 
-	/*! \brief C1_AddOns instance constructor.
-	 *
-	 * \param c4: `C4 *` Pointer to a C4 structure.
-	 */
-	C1_AddOns(C4 *c4);
+  /*! \brief C1_AddOns instance constructor.
+   *
+   * \param c4: `C4 *` Pointer to a C4 structure.
+   */
+  C1_AddOns(C4 *c4);
+  /*! \brief C1_AddOns instance constructor copying contents from a preexisting object.
+   *
+   * \param rhs: `C1_AddOns` preexisting object to copy from.
+   */
+  C1_AddOns(const C1_AddOns& rhs);
 
-	//! \brief C1_AddOns instance destroyer.
-	~C1_AddOns();
+  //! \brief C1_AddOns instance destroyer.
+  ~C1_AddOns();
 };
 
 /*! \brief Representation of the FORTRAN C6 blocks.
  */
 class C6 {
 public:
-	//! \brief QUESTION: definition?
-	double *rac3j;
+  //! \brief QUESTION: definition?
+  int lmtpo;
+  //! \brief QUESTION: definition?
+  double *rac3j;
 
-	/*! \brief C6 instance constructor.
-	 *
-	 * \param lmtpo: `int` QUESTION: definition?
-	 */
-	C6(int lmtpo);
+  /*! \brief C6 instance constructor.
+   *
+   * \param lmtpo: `int` QUESTION: definition?
+   */
+  C6(int lmtpo);
+  /*! \brief C6 instance constructor copying contents from preexisting object.
+   *
+   * \param lmtpo: `int` QUESTION: definition?
+   */
+  C6(const C6& rhs);
 
-	/*! \brief C6 instance destroyer.
-	 */
-	~C6();
+  /*! \brief C6 instance destroyer.
+   */
+  ~C6();
 };
 
 /*! \brief Representation of the FORTRAN C9 blocks.
  */
 class C9 {
 protected:
-	//! \brief Number of rows in the GIS and GLS matrices
-	int gis_size_0;
-	//! \brief Number of rows in the SAM matrix
-	int sam_size_0;
+  //! \brief Number of rows in the GIS and GLS matrices
+  int gis_size_0;
+  //! \brief Number of rows in the SAM matrix
+  int sam_size_0;
+  //! \brief QUESTION: definition?
+  int nlem;
+  //! \brief QUESTION: definition?
+  int nlemt;
 public:
-	//! \brief QUESTION: definition?
-	dcomplex **gis;
-	//! \brief QUESTION: definition?
-	dcomplex **gls;
-	//! \brief QUESTION: definition?
-	dcomplex **sam;
+  //! \brief QUESTION: definition?
+  dcomplex **gis;
+  //! \brief QUESTION: definition?
+  dcomplex **gls;
+  //! \brief QUESTION: definition?
+  dcomplex **sam;
 
-	/*! \brief C9 instance constructor.
-	 *
-	 * \param ndi: `int` QUESTION: definition?
-	 * \param nlem: `int` QUESTION: definition?
-	 * \param ndit: `int` QUESTION: definition?
-	 * \param nlemt: `int` QUESTION: definition?
-	 */
-	C9(int ndi, int nlem, int ndit, int nlemt);
+  /*! \brief C9 instance constructor.
+   *
+   * \param ndi: `int` QUESTION: definition?
+   * \param nlem: `int` QUESTION: definition?
+   * \param ndit: `int` QUESTION: definition?
+   * \param nlemt: `int` QUESTION: definition?
+   */
+  C9(int ndi, int nlem, int ndit, int nlemt);
+  /*! \brief C9 instance constructor copying contents from preexisting object.
+   *
+   * \param rhs: `C9` preexisting object to copy from
+   */
+  C9(const C9& rhs);
 
-	/*! \brief C9 instance destroyer.
-	 */
-	~C9();
+  /*! \brief C9 instance destroyer.
+   */
+  ~C9();
 };
 #endif
diff --git a/src/libnptm/Commons.cpp b/src/libnptm/Commons.cpp
index fc9498a8..24648680 100644
--- a/src/libnptm/Commons.cpp
+++ b/src/libnptm/Commons.cpp
@@ -72,6 +72,78 @@ C1::C1(int ns, int l_max, int *_nshl, int *_iog) {
   }
 }
 
+C1::C1(const C1& rhs) {
+  nlmmt = rhs.nlmmt;
+  nsph = rhs.nsph;
+  lm = rhs.lm;
+
+  rmi = new dcomplex*[lm];
+  rei = new dcomplex*[lm];
+  for (int ri = 0; ri < lm; ri++) {
+    rmi[ri] = new dcomplex[nsph]();
+    rei[ri] = new dcomplex[nsph]();
+    /*! Copy the contents from the template */
+    for (int rj=0; rj<nsph; rj++) {
+      rmi[ri][rj] = rhs.rmi[ri][rj];
+      rei[ri][rj] = rhs.rei[ri][rj];
+    }
+  }
+  w = new dcomplex*[nlmmt];
+  for (int wi = 0; wi < nlmmt; wi++) {
+    w[wi] = new dcomplex[4]();
+    for (int wj=0; wj<4; wj++) w[wi][wj] = rhs.w[wi][wj];
+  }
+  configurations = rhs.configurations;
+  vints = new dcomplex*[nsph];
+  rc = new double*[configurations];
+  nshl = new int[configurations]();
+  for (int ni=0; ni<configurations; ni++) nshl[ni] = rhs.nshl[ni];
+  iog = new int[nsph]();
+  for (int ni=0; ni<nsph; ni++) iog[ni] = rhs.iog[ni];
+  int conf_index = 0;
+  for (int vi = 0; vi < nsph; vi++) {
+    vints[vi] = new dcomplex[16]();
+    for (int vj=0; vj<16; vj++) vints[vi][vj] = rhs.vints[vi][vj];
+  }
+  for (int ri=0; ri<configurations; ri++) {
+    rc[ri] = new double[nshl[ri]]();
+    for (int rj=0; rj<nshl[ri]; rj++) rc[ri][rj] = rhs.rc[ri][rj];
+  }
+  fsas = new dcomplex[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]();
+
+  sas = new dcomplex**[nsph];
+  for (int si = 0; si < nsph; si++) {
+    sas[si] = new dcomplex*[2];
+    for (int sj=0; sj<2; sj++) {
+      sas[si][sj] = new dcomplex[2]();
+      for (int sk=0; sk<2; sk++) sas[si][sj][sk] = rhs.sas[si][sj][sk];
+    }
+    fsas[si] = rhs.fsas[si];
+    sscs[si] = rhs.sscs[si];
+    sexs[si] = rhs.sexs[si];
+    sabs[si] = rhs.sabs[si];
+    sqscs[si] = rhs.sqscs[si];
+    sqexs[si] = rhs.sqexs[si];
+    sqabs[si] = rhs.sqabs[si];
+    gcsv[si] = rhs.gcsv[si];
+    rxx[si] = rhs.rxx[si];
+    ryy[si] = rhs.ryy[si];
+    rzz[si] = rhs.rzz[si];
+    ros[si] = rhs.ros[si];    
+  }
+}
+
 C1::~C1() {
   delete[] rmi;
   delete[] rei;
@@ -114,11 +186,17 @@ C1_AddOns::C1_AddOns(C4 *c4) {
   nsph = c4->nsph;
   lmpo = c4->lmpo;
   nlemt = 2 * c4->nlem;
-  vh = new dcomplex[(nsph * nsph - 1) * c4->litpo]();
-  vj0 = new dcomplex[nsph * c4->lmtpo]();
+  litpo = c4->litpo;
+  lmtpo = c4->lmtpo;
+  litpos = c4->litpos;
+  lmtpos = c4->lmtpos;
+  nv3j = c4->nv3j;
+  lm = c4->lm;
+  vh = new dcomplex[(nsph * nsph - 1) * litpo]();
+  vj0 = new dcomplex[nsph * lmtpo]();
   vj = new dcomplex[1](); // QUESTION: is 1 really enough for a general case?
-  vyhj = new dcomplex[(nsph * nsph - 1) * c4->litpos]();
-  vyj0 = new dcomplex[nsph * c4->lmtpos]();
+  vyhj = new dcomplex[(nsph * nsph - 1) * litpos]();
+  vyj0 = new dcomplex[nsph * lmtpos]();
   am0m = new dcomplex*[nlemt];
   for (int ai = 0; ai < nlemt; ai++) {
     am0m[ai] = new dcomplex[nlemt]();
@@ -140,12 +218,103 @@ C1_AddOns::C1_AddOns(C4 *c4) {
   ecscp = new dcomplex[2]();
   scscpm = new dcomplex[2]();
   ecscpm = new dcomplex[2]();
-  allocate_vectors(c4);
+  v3j0 = new double[nv3j]();
+  ind3j = new int*[lmpo];
+  for (int ii = 0; ii < lmpo; ii++) ind3j[ii] = new int[lm]();
+  //allocate_vectors(c4);
+  sscs = new double[nsph]();
+  ecscm = new double[2]();
+  scscm = new double[2]();
+  scsc = new double[2]();
+  ecsc = new double[2]();
+}
+
+C1_AddOns::C1_AddOns(const C1_AddOns& rhs) {
+  nsph = rhs.nsph;
+  lmpo = rhs.lmpo;
+  nlemt = rhs.nlemt;
+  litpo = rhs.litpo;
+  lmtpo = rhs.lmtpo;
+  litpos = rhs.litpos;
+  lmtpos = rhs.lmtpos;
+  nv3j = rhs.nv3j;
+  lm = rhs.lm;
+  int vhsize=(nsph * nsph - 1) * litpo;
+  vh = new dcomplex[vhsize]();
+  for (int hi=0; hi<vhsize; hi++) vh[hi] = rhs.vh[hi];
+  int vj0size=nsph * lmtpo;
+  vj0 = new dcomplex[vj0size]();
+  for (int hi=0; hi<vj0size; hi++) vj0[hi] = rhs.vj0[hi];
+  vj = new dcomplex[1](); // QUESTION: is 1 really enough for a general case?
+  vj[0] = rhs.vj[0];
+  int vyhjsize=(nsph * nsph - 1) * litpos;
+  vyhj = new dcomplex[vyhjsize]();
+  for (int hi=0; hi<vyhjsize; hi++) vyhj[hi] = rhs.vyhj[hi];
+  int vyj0size = nsph * lmtpos;
+  vyj0 = new dcomplex[vyj0size]();
+  for (int hi=0; hi<vyj0size; hi++) vyj0[hi] = rhs.vyj0[hi];
+  am0m = new dcomplex*[nlemt];
+  for (int ai = 0; ai < nlemt; ai++) {
+    am0m[ai] = new dcomplex[nlemt]();
+    for (int aj = 0; aj < nlemt; aj++) am0m[ai][aj] = rhs.am0m[ai][aj];
+  }
+  vint = new dcomplex[16](); // QUESTION: is dimension 16 generally fixed?
+  vintm = new dcomplex[16]();
+  vintt = new dcomplex[16]();
+  for (int hi=0; hi<16; hi++) {
+    vint[hi] = rhs.vint[hi];
+    vintm[hi] = rhs.vintm[hi];
+    vintt[hi] = rhs.vintt[hi];
+  }
+  vints = new dcomplex*[nsph];
+  for (int vi = 0; vi < nsph; vi++) {
+    vints[vi] = new dcomplex[16]();
+    for (int hj=0; hj<16; hj++) vints[vi][hj] = rhs.vints[vi][hj];
+  }
+  fsac = new dcomplex*[2];
+  sac = new dcomplex*[2];
+  fsacm = new dcomplex*[2];
+  for (int fi = 0; fi < 2; fi++) {
+    fsac[fi] = new dcomplex[2]();
+    sac[fi] = new dcomplex[2]();
+    fsacm[fi] = new dcomplex[2]();
+    for (int fj = 0; fj < 2; fj++) {
+      fsac[fi][fj] = rhs.fsac[fi][fj];
+      sac[fi][fj] = rhs.sac[fi][fj];
+      fsacm[fi][fj] = rhs.fsacm[fi][fj];
+    }
+  }
+  scscp = new dcomplex[2]();
+  ecscp = new dcomplex[2]();
+  scscpm = new dcomplex[2]();
+  ecscpm = new dcomplex[2]();
+  v3j0 = new double[nv3j]();
+  for (int hi=0; hi<nv3j; hi++) v3j0[hi] = rhs.v3j0[hi];
+  ind3j = new int*[lmpo];
+  for (int ii = 0; ii < lmpo; ii++) {
+    ind3j[ii] = new int[lm]();
+    for (int ij=0; ij<lm; ij++) ind3j[ii][ij] = rhs.ind3j[ii][ij];
+  }
+  // vyhj = new dcomplex[int vyhjsize = (nsph * nsph - 1) * litpos]();
+  // for (int hi=0; hi<vyhjsize; hi++) vyhj[hi] = rhs.vyhj[hi]
+  // vyj0 = new dcomplex[int vyj0size=lmtpos * nsph]();
+  // for (int hi=0; hi<vyj0size; hi++) vyj0[hi] = rhs.vyj0[hi];
   sscs = new double[nsph]();
+  for (int hi=0; hi<nsph; hi++) sscs[hi] = rhs.sscs[hi];
   ecscm = new double[2]();
   scscm = new double[2]();
   scsc = new double[2]();
   ecsc = new double[2]();
+  for (int hi=0; hi<2; hi++) {
+    scscp[hi] = rhs.scscp[hi];
+    ecscp[hi] = rhs.ecscp[hi];
+    scscpm[hi] = rhs.scscpm[hi];
+    ecscpm[hi] = rhs.ecscpm[hi];
+    ecscm[hi] = rhs.ecscm[hi];
+    scscm[hi] = rhs.scscm[hi];
+    scsc[hi] = rhs.scsc[hi];
+    ecsc[hi] = rhs.ecsc[hi];
+  }
 }
 
 C1_AddOns::~C1_AddOns() {
@@ -184,25 +353,26 @@ C1_AddOns::~C1_AddOns() {
   delete[] ecsc;
 }
 
-void C1_AddOns::allocate_vectors(C4 *c4) {
-  // This calculates the size of v3j0
-  int lm = (c4->li > c4->le) ? c4->li : c4->le;
-  const int nv3j = c4->nv3j;
-  v3j0 = new double[nv3j]();
-  ind3j = new int*[lmpo];
-  for (int ii = 0; ii < lmpo; ii++) ind3j[ii] = new int[c4->lm]();
-  // This calculates the size of vyhj
-  int ivy = (nsph * nsph - 1) * c4->litpos;
-  vyhj = new dcomplex[ivy]();
-  // This calculates the size of vyj0
-  ivy = c4->lmtpos * c4->nsph;
-  vyj0 = new dcomplex[ivy]();
-}
+// void C1_AddOns::allocate_vectors(C4 *c4) {
+//   // This calculates the size of v3j0
+//   //int lm = (c4->li > c4->le) ? c4->li : c4->le;
+//   const int nv3j = c4->nv3j;
+//   v3j0 = new double[nv3j]();
+//   ind3j = new int*[lmpo];
+//   for (int ii = 0; ii < lmpo; ii++) ind3j[ii] = new int[c4->lm]();
+//   // This calculates the size of vyhj
+//   // int ivy = (nsph * nsph - 1) * c4->litpos;
+//   // vyhj = new dcomplex[ivy]();
+//   // This calculates the size of vyj0
+//   // ivy = c4->lmtpos * c4->nsph;
+//   // vyj0 = new dcomplex[ivy]();
+// }
 
-C2::C2(int ns, int nl, int npnt, int npntts) {
+C2::C2(int ns, int _nl, int npnt, int npntts) {
   nsph = ns;
   int max_n = (npnt > npntts) ? npnt : npntts;
   nhspo = 2 * max_n - 1;
+  nl = _nl;
   ris = new dcomplex[nhspo]();
   dlri = new dcomplex[nhspo]();
   vkt = new dcomplex[nsph]();
@@ -210,6 +380,26 @@ C2::C2(int ns, int nl, int npnt, int npntts) {
   vsz = new double[nsph]();
 }
 
+C2::C2(const C2& rhs) {
+  nsph = rhs.nsph;
+  nhspo = rhs.nhspo;
+  nl = rhs.nl;
+  ris = new dcomplex[nhspo]();
+  dlri = new dcomplex[nhspo]();
+  for (int ind=0; ind<nhspo; ind++) {
+    ris[ind] = rhs.ris[ind];
+    dlri[ind] = rhs.dlri[ind];
+  }
+  vkt = new dcomplex[nsph]();
+  vsz = new double[nsph]();
+  for (int ind=0; ind<nsph; ind++) {
+    vkt[ind] = rhs.vkt[ind];
+    vsz[ind] = rhs.vsz[ind];
+  }
+  dc0 = new dcomplex[nl]();
+  for (int ind=0; ind<nl; ind++) dc0[ind] = rhs.dc0[ind];
+}
+
 C2::~C2() {
   delete[] ris;
   delete[] dlri;
@@ -229,21 +419,80 @@ C3::C3() {
   acs = 0.0;
 }
 
+C3::C3(const C3& rhs) {
+  tsas = new dcomplex*[2];
+  for (int ti=0; ti<2; ti++) {
+    tsas[ti] = new dcomplex[2];
+    for (int tj=0; tj<2; tj++) tsas[ti][tj] = rhs.tsas[ti][tj];
+  }
+  tfsas = rhs.tfsas;
+  gcs = rhs.gcs;
+  scs = rhs.scs;
+  ecs = rhs.ecs;
+  acs = rhs.acs;
+}
+
 C3::~C3() {
   delete[] tsas[1];
   delete[] tsas[0];
   delete[] tsas;
 }
 
-C6::C6(int lmtpo) {
+C4::C4(int _li, int _le, int _nsph) {
+    li = _li;
+    le = _le;
+    lm = (li > le) ? li : le;
+    nv3j = (lm * (lm  + 1) * (2 * lm + 7)) / 6;
+    nsph = _nsph;
+    // The following is needed to initialize C1_AddOns
+    litpo = li + li + 1;
+    litpos = litpo * litpo;
+    lmtpo = li + le + 1;
+    lmtpos = lmtpo * lmtpo;
+    nlim = li * (li + 2);
+    nlem = le * (le + 2);
+    lmpo = lm + 1;
+}
+
+C4::C4(const C4& rhs) {
+    li = rhs.li;
+    le = rhs.le;
+    lm = rhs.lm;
+    nv3j = rhs.nv3j;
+    nsph = rhs.nsph;
+    // The following is needed to initialize C1_AddOns
+    litpo = rhs.litpo;
+    litpos = rhs.litpos;
+    lmtpo = rhs.lmtpo;
+    lmtpos = rhs.lmtpos;
+    nlim = rhs.nlim;
+    nlem = rhs.nlem;
+    lmpo = rhs.lmpo;
+}
+
+C4::~C4() {
+}
+
+C6::C6(int _lmtpo) {
+  lmtpo = _lmtpo;
   rac3j = new double[lmtpo]();
 }
 
+C6::C6(const C6& rhs) {
+  lmtpo = rhs.lmtpo;
+  rac3j = new double[lmtpo]();
+  for (int ri=0; ri<lmtpo; ri++) rac3j[ri] = rhs.rac3j[ri];
+}
+
 C6::~C6() {
   delete[] rac3j;
 }
 
-C9::C9(int ndi, int nlem, int ndit, int nlemt) {
+C9::C9(int ndi, int _nlem, int ndit, int _nlemt) {
+  gis_size_0 = ndi;
+  sam_size_0 = ndit;
+  nlem = _nlem;
+  nlemt = _nlemt;
   gis = new dcomplex*[ndi];
   gls = new dcomplex*[ndi];
   for (int gi = 0; gi < ndi; gi++) {
@@ -252,8 +501,28 @@ C9::C9(int ndi, int nlem, int ndit, int nlemt) {
   }
   sam = new dcomplex*[ndit];
   for (int si = 0; si < ndit; si++) sam[si] = new dcomplex[nlemt]();
-  gis_size_0 = ndi;
-  sam_size_0 = ndit;
+}
+
+C9::C9(const C9& rhs) {
+  gis_size_0 = rhs.gis_size_0;
+  sam_size_0 = rhs.sam_size_0;
+  nlem = rhs.nlem;
+  nlemt = rhs.nlemt;
+  gis = new dcomplex*[gis_size_0];
+  gls = new dcomplex*[gis_size_0];
+  for (int gi = 0; gi < gis_size_0; gi++) {
+    gis[gi] = new dcomplex[nlem]();
+    gls[gi] = new dcomplex[nlem]();
+    for (int gj=0; gj<nlem; gj++) {
+      gis[gi][gj] = rhs.gis[gi][gj];
+      gls[gi][gj] = rhs.gls[gi][gj];
+    }
+  }
+  sam = new dcomplex*[sam_size_0];
+  for (int si = 0; si < sam_size_0; si++) {
+    sam[si] = new dcomplex[nlemt]();
+    for (int sj=0; sj<nlemt; sj++) sam[si][sj] = rhs.sam[si][sj];
+  }
 }
 
 C9::~C9() {
-- 
GitLab