Commit 8ba0e641 authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Prepare C++ implementation of INDME function

parent 8b8fa1b1

src/include/Commons.h

+1 −1
Original line number Diff line number Diff line
@@ -61,7 +61,7 @@ public: 
  dcomplex *dc0;

  //! \brief QUESTION: definition?

  dcomplex *vkt;

  //! Vector of scaled sizes. QUESTION: correct?

  //! Vector of sphere sizes in units of 2*PI/LAMBDA

  double *vsz;



  /*! \brief C2 instance constructor.

src/include/inclu_subs.h

+5 −5
Original line number Diff line number Diff line
@@ -53,18 +53,18 @@ void incms(dcomplex **am, dcomplex **at, double enti, ParticleDescriptor *c1); 
 * \param npnt: `int` TBD.

 * \param npntts: `int` TBD.

 * \param vk: `double` Vacuum wave vector magnitude.

 * \param ent: `double` TBD.

 * \param ent: `dcomplex` TBD.

 * \param enti: `double` TBD.

 * \param entn: `double` TBD.

 * \param entn: `dcomplex` TBD.

 * \param jer: `int &` Error code flag.

 * \param lcalc: `int &` Maximum order achieved in calculation.

 * \param arg: `double &` TBD.

 * \param arg: `dcomplex &` TBD.

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

 * \param c2: `C2 *` Pointer to a C2 instance.

 */

void indme(

	   int i, int npnt, int npntts, double vk, double ent, double enti,

	   double entn, int &jer, int &lcalc, double &arg, ParticleDescriptor *c1,

	   int i, int npnt, int npntts, double vk, dcomplex ent, double enti,

	   dcomplex entn, int &jer, int &lcalc, dcomplex &arg, ParticleDescriptor *c1,

	   C2 *c2

	   );

src/libnptm/inclu_subs.cpp

+125 −2
Original line number Diff line number Diff line
@@ -48,8 +48,8 @@ void exma(dcomplex **am, dcomplex **at, ParticleDescriptor *c1) { 
}



void indme(

	   int i, int npnt, int npntts, double vk, double ent, double enti,

	   double entn, int &jer, int &lcalc, double &arg, ParticleDescriptor *c1,

	   int i, int npnt, int npntts, double vk, dcomplex ent, double enti,

	   dcomplex entn, int &jer, int &lcalc, dcomplex &arg, ParticleDescriptor *c1,

	   C2 *c2

) {

  const dcomplex uim = 0.0 + I * 1.0;
@@ -57,9 +57,132 @@ void indme( 
  const int nstpts = npntts - 1;

  const int lipo = c1->li + 1;

  const int lipt = c1->li + 2;

  dcomplex *cfj = new dcomplex[lipt]();

  dcomplex *cfn = new dcomplex[lipt]();

  dcomplex *fb = new dcomplex[lipt]();

  dcomplex *fbi = new dcomplex[lipt]();

  dcomplex *fn = new dcomplex[lipt]();

  double *rfj = new double[lipt]();

  double *rfn = new double[lipt]();

  jer = 0;

  double sz = vk * c1->ros[i - 1];

  c2->vsz[i - 1] = sz;

  double vkr1 = vk * c1->rc[i - 1][0];

  int nsh = c1->nshl[i . 1];

  c2->vkt[i - 1] = sqrt(c2->dc0[0]);

  arg = vkr1 * c2->vkt[i - 1];

  dcomplex arin = arg;

  if (imag(arg) != 0.0) {

    cbf(lipo, arg, lcalc, cfj);

    if (lcalc < lipo) {

      jer = 5;

      delete[] cfj;

      delete[] cfn;

      delete[] fb;

      delete[] fbi;

      delete[] fn;

      delete[] rfj;

      delete[] rfn;

      return;

    }

    // label 122

    for (int j124 = 0; j124 < lipt; j124++) fbi[j124] = cfj[j124];

  } else { // label 126

    double rarg = real(arg);

    rbf(lipo, rarg, lcalc, rfj);

    if (lcalc < lipo) {

      jer = 5;

      delete[] cfj;

      delete[] cfn;

      delete[] fb;

      delete[] fbi;

      delete[] fn;

      delete[] rfj;

      delete[] rfn;

      return;

    }

    // label 128

    for (int j130 = 0; j130 < lipt; j130++) fbi[j130] = cfj[j130];

  }

  // label 132

  dcomplex aris = sz * entn;

  arg = aris;

  if (enti != 0.0) {

    cbf(lipo, arg, lcalc, cfj);

    if (lcalc < lipo) {

      jer = 11;

      delete[] cfj;

      delete[] cfn;

      delete[] fb;

      delete[] fbi;

      delete[] fn;

      delete[] rfj;

      delete[] rfn;

      return;

    }

    cnf(lipo, arg, lcalc, cfj, cfn);

    // QUESTION: should we check for lcalc and throw JER=12 if failing?

    // see lines 2492 -2505 in INCLU.F (test done in REAL case but not in COMPLEX case).

    for (int j143 = 0; j143 < lipt; j143++) {

      fb[j143] = cfj[j143];

      fn[j143] = cfn[j143];

    }

  } else { // label 145

    double rarg = real(aris);

    rbf(lipo, rarg, lcalc, rfj);

    if (lcalc < lipo) {

      jer = 11;

      delete[] cfj;

      delete[] cfn;

      delete[] fb;

      delete[] fbi;

      delete[] fn;

      delete[] rfj;

      delete[] rfn;

      return;

    }

    rnf(lipo, rarg, lcalc, rfn);

    if (lcalc < lipo) {

      jer = 12;

      delete[] cfj;

      delete[] cfn;

      delete[] fb;

      delete[] fbi;

      delete[] fn;

      delete[] rfj;

      delete[] rfn;

      return;

    }

    for (int j150 = 0; j150 < lipt; j150++) {

      fb[j150] = rfj[j150];

      fn[j150] = rfn[j150];

    }

  }

  // label 152

  if (nsh < 2) { // nsh == 1

    dcomplex cri = c2->dc0[0] / ent;

    for (int l160 = 1; l160 <= c1->li; l160++) {

      int lpo = l160 + 1;

      int ltpo = lpo + l160;

      int lpt = lpo + 1;

      dcomplex dfbi = (l160 * fbi[l160 - 1] - lpo * fbi[lpt - 1]) * arin + fbi[lpo - 1] * ltpo;

      dcomplex dfb = (l160 * fb[l160 - 1] - lpo * fb[lpt - 1]) * aris + fb[lpo - 1] * ltpo;

      dcomplex dfn = (l160 * fn[l160 - 1] - lpo * fn[lpt - 1]) * aris + fn[lpo - 1] * ltpo;

      dcomplex ccna = fbi[lpo - 1] * dfn;

      dcomplex ccnb = fn[lpo - 1] * dfbi;

      // DA INCLU.F: 2520

    } // l160 loop

  } else { // label 165: nsh > 1

    

  } // nsh if



  delete[] cfj;

  delete[] cfn;

  delete[] fb;

  delete[] fbi;

  delete[] fn;

  delete[] rfj;

  delete[] rfn;

}



void incms(dcomplex **am, dcomplex **at, double enti) {