Commit 48a977f2 authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Add C++ porting of INCMS and EXMA

parent 8a30c270

src/include/inclu_subs.h

+39 −0
Original line number Diff line number Diff line
@@ -29,6 +29,45 @@ 
#ifndef INCLUDE_INCLU_SUBS_H_

#define INCLUDE_INCLU_SUBS_H_



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

 *

 * \param am: `dcomplex **` Field transition coefficients matrix.

 * \param at: `dcomplex **` TBD.

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

 */

void exma(dcomplex **am, dcomplex **at, ParticleDescriptor *c1);



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

 *

 * \param am: `dcomplex **` Field transition coefficients matrix.

 * \param at: `dcomplex **` TBD.

 * \param enti: `double` TBD.

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

 */

void incms(dcomplex **am, dcomplex **at, double enti, ParticleDescriptor *c1);



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

 *

 * \param li: `int` Maximum internal field expansion order.

 * \param i: `int` 1-based sphere configuration index.

 * \param npnt: `int` TBD.

 * \param npntts: `int` TBD.

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

 * \param ent: `double` TBD.

 * \param enti: `double` TBD.

 * \param entn: `double` TBD.

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

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

 * \param arg: `double &` 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,

	   C2 *c2

	   );



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

 *

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

src/libnptm/inclu_subs.cpp

+240 −0
Original line number Diff line number Diff line
@@ -33,6 +33,246 @@ 


using namespace std;



void exma(dcomplex **am, dcomplex **at, ParticleDescriptor *c1) {

  const dcomplex cc0 = 0.0 + I * 0.0;

  const ndit = c1->nsph * c1->nlim * 2;

  const int ndm = ndit + c1->nlemt;

  for (int j20 = 1; j20 <= c1->nlemt; j20++) {

    int j0 = j20 + ndit;

    for (int i20 = 1; i20 <= c1->nlemt; i20++) {

      dcomplex sum = cc0;

      for (int k = 0; k < ndm; k++) sum += at[i20 - 1][k] * am[k][j0 - 1];

    }

    c1->am0m[i20 - 1][j20 - 1] = sum;

  } // j20 loop

}



void indme(

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

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

	   C2 *c2

) {

  const dcomplex uim = 0.0 + I * 1.0;

  const int nstp = npnt - 1;

  const int nstpts = npntts - 1;

  const int lipo = c1->li + 1;

  const int lipt = c1->li + 2;

  jer = 0;

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

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

}



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

  const dcomplex cc0 = 0.0 + I * 0.0;

  int nbl;

  const int ndi = c1->nsph * c1->nlim;

  const int ndit = ndi + ndi;

  const int ndm = ndit + c1->nlemt;

  nbl = 0;

  for (int n1 = 1, n1 < c1->nsph; n1++) {

    int in1 = (n1 - 1) * c1->nlim;

    int n1po = n1 + 1;

    for (int n2 = n1po; n2 <= c1->nsph; n2++) {

      int in2 = (n2 - 1) * c1->nlim;

      nbl++;

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

	int l1po = l1 + 1;

	int il1 = l1po * l1;

	int l1tpo = l1po + l1;

	for (int im1 = 1; im1 <= l1tpo; im1++) {

	  int m1 = im1 - l1po;

	  int ilm1 = il1 + m1;

	  int ilm1e = ilm1 + ndi;

	  int i1 = in1 + ilm1;

	  int i1e = in1 + ilm1e;

	  int j1 = in2 + ilm1;

	  int j1e = in2 + ilm1e;

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

	    int l2po = l2 + 1;

	    int il2 = l2po * l2;

	    int l2tpo = l2po + l2;

	    int ish = ((l2 + l1) % 2 == 0) ? 1 : -1;

	    int isk = -ish;

	    for (int im2 = 1; im2 <= l2tpo; im2++) {

	      int m2 = im2 - l2po;

	      int ilm2 = il2 + m2;

	      int ilm2e = ilm2 + ndi;

	      int i2 = in2 + ilm2;

	      int i2e = in2 + ilm2e;

	      int j2 = in1 + ilm2;

	      int j2e = in1 + ilm2e;

	      dcomplex cgh = ghit(0, 0, nbl, l1, m1, l2, m2);

	      dcomplex cgk = ghit(0, 1, nbl, l1, m1, l2, m2);

	      am[i1 - 1][i2 - 1] = cgh;

	      am[i1 - 1][i2e - 1] = cgk;

	      am[i1e - 1][i2 - 1] = cgk;

	      am[i1e - 1][i2e - 1] = cgh;

	      am[j1 - 1][j2 - 1] = cgh * ish;

	      am[j1 - 1][j2e - 1] = cgk * isk;

	      am[j1e - 1][j2 - 1] = cgk * isk;

	      am[j1e - 1][j2e - 1] = cgh * ish;

	    } // im2 loop 24

	  } // l2 loop 24

	} // im1 loop 24

      } // l1 loop 24

    } // n2 loop 26

  } // n1 loop 26

  for (int n1 = 1; n1 <= c1->nsph; n1++) {

    int in1 = (n1 - 1) * c1->nlim;

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

      dcomplex frm = c1->rmi[l1 - 1][n1 - 1];

      dcomplex fre = c1->rei[l1 - 1][n1 - 1];

      int l1po = l1 + 1;

      int il1 = l1po * l1;

      int l1tpo = l1po + l1;

      for (int im1 = 1; im1 <= l1tpo; im1++) {

	int m1 = im1 - l1po;

	int ilm1 = il1 + m1;

	int i1e = i1 + ndi;

	for (int ilm2 = 1; ilm2 <= c1->nlim; ilm2++) {

	  int i2 = in1 + ilm2;

	  i2e = i2 + ndi;

	  am[i1 - 1][i2 - 1] = cc0;

	  am[i1 - 1][i2e - 1] = cc0;

	  am[i1e - 1][i2 - 1] = cc0;

	  am[i1e - 1][i2e - 1] = cc0;

	} // ilm2 loop 28

	am[i1 - 1][i1 - 1] = frm;

	am[i1e - 1][i1e - 1]= fre;

      } // im1 loop 30

    } // l1 loop 30

  } // n1 loop 30

  int nditpo = ndit + 1;

  for (int i1 = 1; i1 <= c1->nlemt; i1++) {

    int i3 = i1 + ndit;

    for (int i2 = nditpo; i2 <= ndm; i2++) {

      am[i3 - 1][i2 - 1] = cc0;

      at[i1 - 1][i2 - 1] = cc0;

    } // i2 loop 40

  } // i1 loop 40

  { // CODE BLOCK TO PRESERVE SCOPE OF i1

    int i1 = 0;

    for (int l1 = 1; l1 <= c1->le; l1++) {

      dcomplex frm = c1->rm0[l1 - 1];

      dcomplex fre = c1->re0[l1 - 1];

      dcomplex ftm = c1->tm0[l1 - 1];

      dcomplex fte = c1->te0[l1 - 1];

      int l1tpo = l1 + l1 + 1;

      for (int im1  = 1; im1 <= l1tpo; im1 ++) {

	i1++;

	int i1e = i1 + c1->nlem;

	int i3 = i1 + ndit;

	int i3e = i3 + c1->nlem;

	am[i3 - 1][i3 - 1] = frm;

	am[i3e - 1][i3e - 1] = fre;

	at[i1 - 1][i3 - 1] = ftm;

	at[i1e - 1][i3e - 1] = fte;

      } // im1 loop 45

    } // l1 loop 45

  } // END OF i1 INCREMENTAL DEFINITION

  if (enti != 0.0) {

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

      dcomplex frm = c1->rmw[l2 - 1];

      dcomplex fre = c1->rew[l2 - 1];

      dcomplex ftm = c1->tm[l2 - 1];

      dcomplex fte = c1->te[l2 - 1];

      int l2po = l2 + 1;

      int il2 = l2po * l2;

      int l2tpo = l2po + l2;

      for (int im2 = 1; im2 <= l2tpo; im2++) {

	int m2 = im2 - l2po;

	int i2 = il2 + m2;

	int j2 = il2 - m2;

	int i2e = i2 + c1->nlem;

	int j2e = j2 + c1->nlem;

	int i3 = i2 + ndit;

	int j3 = j2 + ndit;

	int i3e = i3 + c1->nlem;

	int j3e = j3 + c1->nlem;

	for (int n1 = 1; n1 <= c1->nsph; n1++) {

	  int in1 = (n1 - 1) * c1->nlim;

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

	    int l1po = l1 + 1;

	    int il1 = l1po * l1;

	    int l1tpo = l1po + l1;

	    for (int im1 = 1; im1 <= l1tpo; im1++) {

	      int m1 = im1 - l1po;

	      int ilm1 = il1 + m1;

	      int jlm1 = il1 - m1;

	      int i1 = in1 + ilm1;

	      int i1e = i1 + ndi;

	      int j1 = in1 + jlm1;

	      int j1e = j1 + ndi;

	      int isil = ((m2 + m1) % 2 == 0) ? 1 : -1;

	      dcomplex cgi = ghit(2, 0, n1, l1, m1, l2, m2);

	      dcomplex cgl = ghit(2, 1, n1, l1, m1, l2, m2);

	      am[i1 - 1][i3 - 1] = cgi;

	      am[i1 - 1][i3e - 1] = cgl;

	      am[i1e - 1][i3 - 1] = cgl;

	      am[i1e - 1][i3e - 1] = cgi;

	      am[j3 - 1][j1 - 1] = cgi * frm * isil;

	      am[j3 - 1][j1e - 1] = cgl * frm * isil;

	      am[j3e - 1][j1 - 1] = cgl * fre * isil;

	      am[j3e - 1][j1e - 1] = cgi * fre * isil;

	      at[j2 - 1][j1 - 1] = cgi * ftm * isil;

	      at[j2 - 1][j1e - 1] = cgl * ftm * isil;

	      at[j2e - 1][j1 - 1] = cgl * fte * isil;

	      at[j2e - 1][j1e - 1] = cgi * fte * isil;

	      // returns

	    } // im1 loop 50

	  } // l1 loop 50

	} // n1 loop 50

      } // im2 loop 50

    } // l2 loop 50

  } else {

    // label 55

    int i2 = 0;

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

      dcomplex frm = c1->rmw[l2 - 1];

      dcomplex fre = c1->rew[l2 - 1];

      dcomplex ftm = c1->tm[l2 - 1];

      dcomplex fte = c1->te[l2 - 1];

      int l2tpo = l2 + l2 + 1;

      int m2 = -m2 - 1;

      for (int im2 = 1; im2 <= l2tpo; im2++) {

	m2++;

	i2++;

	int i2e = i2 + c1->nlem;

	int i3 = i2 + ndit;

	int i3e = i3 + nlem;

	int i1 = 0;

	for (int n1 = 1; n1 <= c1->nsph; n1++) {

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

	    l1tpo = l1 + l1 + 1;

	    int m1 = -l1 - 1;

	    for (int im1 = 1; im1 <= l1tpo; im1++) {

	      m1++;

	      i1++;

	      int i1e = i1 + ndi;

	      dcomplex cgi = ghit(2, 0, n1, l1, m1, l2, m2);

	      dcomplex cgl = ghit(2, 1, n1, l1, m1, l2, m2);

	      am[i1 - 1][i3 - 1] = cgi;

	      am[i1 - 1][i3e - 1] = cgl;

	      am[i1e - 1][i3 - 1] = cgl;

	      am[i1e - 1][i3e - 1] = cgi;

	      cgi = dconjg(cgi);

	      cgl = dconjg(cgl);

	      am[i3 - 1][i1 - 1] = cgi * frm;

	      am[i3 - 1][i1e - 1] = cgl * frm;

	      am[i3e - 1][i1 - 1] = cgl * fre;

	      am[i3e - 1][i1e - 1] = cgi * fre;

	      at[i2 - 1][i1 - 1] = cgi * ftm;

	      at[i2 - 1][i1e - 1] = cgl * ftm;

	      at[i2e - 1][i1 - 1] = cgl * fte;

	      at[i2e - 1][i1e - 1] cgi * fte;

	    } // im1 loop 60

	  } // l1 loop 60

	} // n1 loop 60

      } // im2 loop 60

    } // l2 loop 60

  } // END OF enti = 0.0 CASE

}



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]();