Commons.h

/* Copyright (C) 2024   INAF - Osservatorio Astronomico di Cagliari

   This program is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.
   
   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
   
   A copy of the GNU General Public License is distributed along with
   this program in the COPYING file. If not, see: <https://www.gnu.org/licenses/>.
 */

/*! \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 dynamically 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 INCLUDE_COMMONS_H_
#define INCLUDE_COMMONS_H_

#ifdef USE_MPI
#include <mpi.h>
#endif

class mixMPI;

/*! \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 Contiguous RMI vector.
  dcomplex *vec_rmi;
  //! \brief Contiguous REI vector.
  dcomplex *vec_rei;
  //! \brief Contiguous W vector.
  dcomplex *vec_w;
  //! \brief Contiguous VINTS vector
  dcomplex *vec_vints;
  
public:
  //! \brief NLMMT = 2 * LM * (LM + 2).
  int nlmmt;
  //! \brief Number of configurations
  int configurations;
  //! \brief QUESTION: definition?
  dcomplex **rmi;
  //! \brief QUESTION: definition?
  dcomplex **rei;
  //! \brief QUESTION: definition?
  dcomplex **w;
  //! \brief QUESTION: definition?
  dcomplex *fsas;
  //! \brief QUESTION: definition?
  dcomplex *vint;
  //! \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 gconf: `GeometryConfiguration *` Pointer to a geometry configuration object.
   * \param sconf: `ScattererConfiguration *` Pointer to a scatterer configuration object.
   */
  C1(GeometryConfiguration *gconf, ScattererConfiguration *sconf);

  /*! \brief C1 instance constructor copying all contents from a preexisting template
   *
   * \param rhs: `C1` preexisting template.
   */
  C1(const C1& rhs);

#ifdef MPI_VERSION
  /*! \brief C1 instance constructor copying all contents off MPI broadcast from MPI process 0
   *
   * \param mpidata: `mixMPI *` pointer to MPI data structure.
   */
  C1(const mixMPI *mpidata);

  /*! \brief send C1 instance from MPI process 0 via MPI broadcasts to all other processes
   *
   * \param mpidata: `mixMPI *` pointer to MPI data structure.
   */
  void mpibcast(const mixMPI *mpidata);
#endif

  //! \brief C1 instance destroyer.
  ~C1();
};

/*! \brief Representation of the FORTRAN C2 blocks.
 *
 */
class C2 {
protected:
  //! \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 C2 instance constructor.
   *
   * \param gconf: `GeometryConfiguration*` Pointer to a GeometryConfiguration instance.
   * \param sconf: `ScattererConfiguration*` Pointer to a ScattererConfiguration instance.
   */
  C2(GeometryConfiguration *gconf, ScattererConfiguration *sconf);

  /*! \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();

#ifdef MPI_VERSION
  /*! \brief C2 instance constructor copying all contents off MPI broadcast from MPI process 0
   *
   * \param mpidata: `mixMPI *` pointer to MPI data structure.
   */
  C2(const mixMPI *mpidata);

  /*! \brief send C2 instance from MPI process 0 via MPI broadcasts to all other processes
   *
   * \param mpidata: `mixMPI *` pointer to MPI data structure.
   */
  void mpibcast(const mixMPI *mpidata);
#endif

};

/*! \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 C3 instance constructor.
   */
  C3();

  /*! \brief C3 instance constructor copying its contents from a preexisting object.
   */
  C3(const C3& rhs);

  /*! \brief C3 instance destroyer.
   */
  ~C3();

#ifdef MPI_VERSION
  /*! \brief C3 instance constructor copying all contents off MPI broadcast from MPI process 0
   *
   * \param mpidata: `mixMPI *` pointer to MPI data structure.
   */
  C3(const mixMPI *mpidata);

  /*! \brief send C3 instance from MPI process 0 via MPI broadcasts to all other processes
   *
   * \param mpidata: `mixMPI *` pointer to MPI data structure.
   */
  void mpibcast(const mixMPI *mpidata);
#endif

};

/*! \brief Representation of the FORTRAN C4 blocks.
 */
class C4 {
public:
  //! \brief LITPO = 2 * LI + 1.
  int litpo;
  //! \brief LITPOS = LITPO * LITPO
  int litpos;
  //! \brief LMPO = LM + 1.
  int lmpo;
  //! \brief LMTPO = 2 * LM + 1.
  int lmtpo;
  //! \brief LMTPOS = LMTPO * LMTPO.
  int lmtpos;
  //! \brief Internal field expansion order.
  int li;
  //! \brief QUESTION: definition?
  int nlim;
  //! \brief External field expansion order.
  int le;
  //! \brief QUESTION: definition?
  int nlem;
  //! \brief Maximum field expansion order.
  int lm;
  //! \brief Number of spheres.
  int nsph;
  //! \brief QUESTION: definition?
  int nv3j;

  /*! \brief C4 instance constructor.
   *
   * \param gconf: `GeometryConfiguration*` Pointer to a GeometryConfiguration instance.
   */
  C4(GeometryConfiguration *gconf);
  
  /*! \brief C4 instance constructor copying its contents from a preexisting object.
   *
   * \param rhs: `C4&` Reference of the object to be copied.
   */
  C4(const C4& rhs);

  /*! \brief C4 instance destroyer.
   */
  ~C4();

#ifdef MPI_VERSION
  /*! \brief C4 instance constructor copying all contents off MPI broadcast from MPI process 0
   *
   * \param mpidata: `mixMPI *` pointer to MPI data structure.
   */
  C4(const mixMPI *mpidata);

  /*! \brief send C4 instance from MPI process 0 via MPI broadcasts to all other processes
   *
   * \param mpidata: `mixMPI *` pointer to MPI data structure.
   */
  void mpibcast(const mixMPI *mpidata);
#endif

};

/*! \brief Vectors and matrices that are specific to cluster C1 blocks.
 *
 */
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 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;

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 *am0v;
  //! \brief QUESTION: definition?
  dcomplex *vintm;
  //! \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 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();

#ifdef MPI_VERSION
  /*! \brief C1_AddOns instance constructor copying all contents off MPI broadcast from MPI process 0
   *
   * \param mpidata: `mixMPI *` pointer to MPI data structure.
   */
  C1_AddOns(const mixMPI *mpidata);

  /*! \brief send C1_AddOns instance from MPI process 0 via MPI broadcasts to all other processes
   *
   * \param mpidata: `mixMPI *` pointer to MPI data structure.
   */
  void mpibcast(const mixMPI *mpidata);
#endif

};

/*! \brief Representation of the FORTRAN C6 blocks.
 */
class C6 {
public:
  //! \brief LMTPO = 2 * LM + 1.
  int lmtpo;
  //! \brief QUESTION: definition?
  double *rac3j;

  /*! \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();

#ifdef MPI_VERSION
  /*! \brief C6 instance constructor copying all contents off MPI broadcast from MPI process 0
   *
   * \param mpidata: `mixMPI *` pointer to MPI data structure.
   */
  C6(const mixMPI *mpidata);

  /*! \brief send C6 instance from MPI process 0 via MPI broadcasts to all other processes
   *
   * \param mpidata: `mixMPI *` pointer to MPI data structure.
   */
  void mpibcast(const mixMPI *mpidata);
#endif

};

/*! \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 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 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();

#ifdef MPI_VERSION
  /*! \brief C9 instance constructor copying all contents off MPI broadcast from MPI process 0
   *
   * \param mpidata: `mixMPI *` pointer to MPI data structure.
   */
  C9(const mixMPI *mpidata);

  /*! \brief send C9 instance from MPI process 0 via MPI broadcasts to all other processes
   *
   * \param mpidata: `mixMPI *` pointer to MPI data structure.
   */
  void mpibcast(const mixMPI *mpidata);
#endif

};

/*! \brief structure with essential MPI data.
 */
class mixMPI {
public:
  //! \brief was MPI initialised?
  bool mpirunning;
  //! \brief MPI rank
  int rank;
  //! \brief MPI nprocs
  int nprocs;

  /*! \brief empty mixMPI instance constructor.
   */
  mixMPI();

  /*! \brief mixMPI instance constructor from an actual MPI communicator.
   */
#ifdef MPI_VERSION
  mixMPI(MPI_Comm comm);
#endif
  
  /*! \brief mixMPI instance constructor copying its contents from a preexisting object.
   */
  mixMPI(const mixMPI& rhs);

  /*! \brief mixMPI instance destroyer.
   */
  ~mixMPI();
};

/*! \brief A data structure representing the information used for a single scale
 * of the CLUSTER case.
 */
class ClusterIterationData {
public:
  //! \brief Pointer to a C1 structure.
  C1 *c1;
  //! \brief Pointer to a C1_AddOns structure.
  C1_AddOns *c1ao;
  //! \brief Pointer to a C2 structure.
  C2 *c2;
  //! \brief Pointer to a C3 structure.
  C3 *c3;
  //! brief Pointer to a C4 structure.
  C4 *c4;
  //! \brief Pointer to a C6 structure.
  C6 *c6;
  //! \brief Pointer to a C9 structure.
  C9 *c9;
  //! \brief Pointer to a formatted output file.
  double *gaps;
  double **tqse;
  dcomplex **tqspe;
  double **tqss;
  dcomplex **tqsps;
  double ****zpv;
  double **gapm;
  dcomplex **gappm;
  double *argi;
  double *args;
  double **gap;
  dcomplex **gapp;
  double **tqce;
  dcomplex **tqcpe;
  double **tqcs;
  dcomplex **tqcps;
  double *duk;
  double **cextlr;
  double **cext;
  double **cmullr;
  double **cmul;
  double *gapv;
  double *tqev;
  double *tqsv;
  double *u;
  double *us;
  double *un;
  double *uns;
  double *up;
  double *ups;
  double *unmp;
  double *unsmp;
  double *upmp;
  double *upsmp;
  double scan;
  double cfmp;
  double sfmp;
  double cfsp;
  double sfsp;
  double qsfi;
  double sqsfi;
  dcomplex *am_vector;
  dcomplex **am;
  dcomplex arg;
  //! \brief Wave vector.
  double vk;
  //! \brief Wave number.
  double wn;
  double xip;
  int number_of_scales;
  int xiblock;
  int firstxi;
  int lastxi;

  ClusterIterationData(GeometryConfiguration *gconf, ScattererConfiguration *sconf, const mixMPI *mpidata);
  
  ClusterIterationData(const ClusterIterationData& rhs);

#ifdef MPI_VERSION
  ClusterIterationData(const mixMPI *mpidata);

  /*! \brief Broadcast over MPI the ClusterIterationData instance from MPI process 0 to all others.
   *
   * When using MPI, the initial ClusterIterationData instance created by MPI process 0
   * needs to be replicated on all other processes. This function sends it using
   * MPI broadcast calls. The MPI broadcast calls in this function must match those
   * in the constructor using the mixMPI pointer.
   *
   * \param mpidata: `mixMPI *` Pointer to the mpi structure used to do the MPI broadcast.
   */
  void mpibcast(const mixMPI *mpidata);
#endif

  ~ClusterIterationData();

};

/*! \brief A data structure representing the angles to be evaluated in the problem.
 *
 */
class ScatteringAngles {
protected:
  //! \brief Number of incident field azimuth angles.
  int _nth;
  //! \brief Number of scattered field azimuth angles.
  int _nths;
  //! \brief Number of incident field elevation angles.
  int _nph;
  //! \brief Number of scattered field elevation angles.
  int _nphs;
  //! \brief Number of incident field propagation angles.
  int _nk;
  //! \brief Number of scattered field propagation angles.
  int _nks;
  //! \brief Total number of field propagation angles.
  int _nkks;
  //! \brief First incident field azimuth angle.
  double _th;
  //! \brief Incident field azimuth angle increment.
  double _thstp;
  //! \brief Last incident field azimuth angle.
  double _thlst;
  //! \brief First scattered field azimuth angle.
  double _ths;
  //! \brief Scattered field azimuth angle increment.
  double _thsstp;
  //! \brief Last scattered field azimuth angle.
  double _thslst;
  //! \brief First incident field elevation angle.
  double _ph;
  //! \brief Incident field elevation angle increment.
  double _phstp;
  //! \brief Last incident field elevation angle.
  double _phlst;
  //! \brief First scattered field elevation angle.
  double _phs;
  //! \brief Scattered field elevation angle increment.
  double _phsstp;
  //! \brief Last scattered field elevation angle.
  double _phslst;
  //! \brief Azimuth scattering deflection.
  double _thsca;

public:
  //! \brief Read only view of `_nth`.
  const int& nth = _nth;
  //! \brief Read only view of `_nths`.
  const int& nths = _nths;
  //! \brief Read only view of `_nph`.
  const int& nph = _nph;
  //! \brief Read only view of `_nphs`.
  const int& nphs = _nphs;
  //! \brief Read only view of `_nk`.
  const int& nk = _nk;
  //! \brief Read only view of `_nks`.
  const int& nks = _nks;
  //! \brief Read only view of `_nkks`.
  const int& nkks = _nkks;
  //! \brief Read only view of `_th`.
  const double& th = _th;
  //! \brief Read only view of `_thstp`.
  const double& thstp = _thstp;
  //! \brief Read only view of `_thlst`.
  const double& thlst = _thlst;
  //! \brief Read only view of `_ths`.
  const double& ths = _ths;
  //! \brief Read only view of `_thsstp`.
  const double& thsstp = _thsstp;
  //! \brief Read only view of `_thslst`.
  const double& thslst = _thslst;
  //! \brief Read only view of `_ph`.
  const double& ph = _ph;
  //! \brief Read only view of `_phstp`.
  const double& phstp = _phstp;
  //! \brief Read only view of `_phlst`.
  const double& phlst = _phlst;
  //! \brief Read only view of `_phs`.
  const double& phs = _phs;
  //! \brief Read only view of `_phsstp`.
  const double& phsstp = _phsstp;
  //! \brief Read only view of `_phslst`.
  const double& phslst = _phslst;
  //! \brief Read only view of `_thsca`.
  const double& thsca = _thsca;

  /*! \brief ScatteringAngles instance constructor.
   *
   * \param gconf: `GeometryConfiguration*` Pointer to a GeometryConfiguration object.
   */
  ScatteringAngles(GeometryConfiguration *gconf);
  
  /*! \brief ScatteringAngles copy constructor.
   *
   * \param rhs: `ScatteringAngles&` Reference to the ScatteringAngles object to be copied.
   */
  ScatteringAngles(const ScatteringAngles &rhs);

#ifdef MPI_VERSION
  /*! \brief ScatteringAngles copy from MPI broadcast.
   *
   * \param mpidata: `mixMPI *` Pointer to the mpidata instance used to copy the data.
   */
  ScatteringAngles(const mixMPI *mpidata);

    /*! \brief Broadcast over MPI the ScatteringAngles instance from MPI process 0 to all others.
   *
   * When using MPI, the initial ScatteringAngles instance created by MPI process 0
   * needs to be replicated on all other processes. This function sends it using
   * MPI broadcast calls. The MPI broadcast calls in this function must match those
   * in the constructor using the mixMPI pointer.
   *
   * \param mpidata: `mixMPI *` Pointer to the mpi structure used to do the MPI broadcast.
   */
  void mpibcast(const mixMPI *mpidata);
#endif

};


#endif