Create inclu_subs library and port INSTR subroutine to C++

/* 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 inclu_subs.h

 *

 * \brief C++ porting of INCLU functions and subroutines.

 *

 * This library includes a collection of functions that are used to solve the

 * scattering problem in the case of a sphere with a cluster of inclusions. Like

 * the other use cases, many of the functions declared here execute various

 * calculations on different data structures. In order to manage access to such

 * variety of calculations, most functions are declared as `void` and they operate

 * on output arguments passed by reference.

 */

#ifndef INCLUDE_INCLU_SUBS_H_

#define INCLUDE_INCLU_SUBS_H_

/*! \brief C++ porting of INSTR.

 *

 * \param sconf: `ScattererConfiguration *` Pointer to a ScattererConfiguration instance.

 * \param c1: `ParticleDescriptor *` Pointer to a ParticleDescriptor instance.

 * \param c6: `C6 *` Pointer to a C6 instance.

 */

void instr(ScattererConfiguration *sconf, ParticleDescriptor *c1, C6 *c6);

#endif // INCLUDE_INCLU_SUBS_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 inclu_subs.cpp

 *

 * \brief C++ implementation of INCLUSION subroutines.

 */

#ifndef INCLUDE_TYPES_H_

#include "../include/types.h"

#endif

#ifndef INCLUDE_CONFIGURATION_H_

#include "../include/Configuration.h"

#endif

#ifndef INCLUDE_COMMONS_H_

#include "../include/Commons.h"

#endif

using namespace std;

void instr(ScattererConfiguration *sconf, ParticleDescriptor *c1, C6 *c6) {

  const int ylm_size = (c1->litpos > c1->lmtpos) ? c1->litpos : c1->lmtpos;

  dcomplex *ylm = new dcomplex[ylm_size]();

  double rx, ry, rz, rr, crth, srth, crph, srph;

  int ivy;

  for (int i18 = 0; i18 < c1->nsph; i18++) {

    int i = i18 + 1;

    if (c1->iog[i18] >= i) {

      int nsh = c1->nshl[i18];

      for (int j = 0; j < nsh; j++)

	c1->rc[i18][j] = sconf->get_rcf(i18, j) * c1->ros[i18];

    }

  } // i18 loop

  int i = 0;

  for (int l1po28 = 1; l1po28 <= c1->lmpo; l1po28++) {

    int l1 = l1po - 1;

    for (int l2 = 1; l2 <= c1->lm; l2++) {

      r3j000(l1, l2);

      c1->ind3j[l1po - 1][l2 - 1] = i;

      int lmnpo = 1 + ((l2 - l1 > 0) ? l2 - l1 : l1 - l2);

      int lmxpo = l2 + l1 + 1;

      int il = 0;

      int lpo = lmnpo;

      while (lpo <= lmxpo) {

	c1->v3j0[i++] = c6->rac3j[il++];

	lpo += 2;

      }

    } // l2 loop

  } // l1po28 loop

  int nsphmo = c1->nsph - 1;

  int lit = c1->li + c1->li;

  ivy = 0;

  for (int nf40 = 0; nf40 < nsphmo; nf40++) {

    int nf = nf40 + 1;

    for (int ns = nf; ns < c1->nsph; ns++) {

      rx = c1->rxx[nf40] - c1->rxx[ns];

      ry = c1->ryy[nf40] - c1->ryy[ns];

      rz = c1->rzz[nf40] - c1->rzz[ns];

      polar(rx, ry, rz, rr, crth, srth, crph, srph);

      sphar(crth, srth, crph, srph, lit, ylm);

      for (int iv38 = 0; iv38 < c1->litpos; iv38++)

	c1->vyhj[iv38 + ivy] = dconjg(ylm[iv38]);

      ivy += c1->litpos;

    } // ns loop

  } // nf40 loop

  int lmt = c1->li + c1->le;

  ivy = 0;

  for (int nf50 = 0; nf50 < c1->nsph; nf50++) {

    rx = c1->rxx[nf50];

    ry = c1->ryy[nf50];

    rz = c1->rzz[nf50];

    if (rx != 0.0 || ry != 0.0 || rz != 0.0) {

      polar(rx, ry, rz, rr, crth, srth, crph, srph);

      sphar(crth, srth, crph, srph, lmt, ylm);

      for (int iv48 = 0; iv48 < c1->lmtpos; iv48++)

	c1->vyj0[iv48 + ivy] = dconjg(ylm[iv48]);

      ivy += c1->lmtpos;

    }

  } // nf50 loop

  delete[] ylm;

}