Commit 0b28e145 authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Add C++ porting of INDME()

parent 8ba0e641

src/libnptm/inclu_subs.cpp

+193 −138
Original line number Diff line number Diff line
@@ -47,144 +47,6 @@ void exma(dcomplex **am, dcomplex **at, ParticleDescriptor *c1) { 
  } // j20 loop

}



void indme(

	   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;

  const int nstp = npnt - 1;

  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) {

  const dcomplex cc0 = 0.0 + I * 0.0;

  int nbl;
@@ -396,6 +258,199 @@ void incms(dcomplex **am, dcomplex **at, double enti) { 
  } // END OF enti = 0.0 CASE

}



void indme(

	   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;

  const int nstp = npnt - 1;

  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

  dcomplex *rmf = new dcomplex[c1->li]();

  dcomplex *drmf = new dcomplex[c1->li]();

  dcomplex *ref = new dcomplex[c1->li]();

  dcomplex *dref = new dcomplex[c1->li]();

  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;

      dcomplex ccnc = fbi[lpo - 1] * dfb;

      dcomplex ccnd = fb[lpo - 1] * dfbi;

      rmi[l60 - 1][i - 1] = 1.0 + uim * (ccna - ccnb) / (ccnc - ccnd);

      rei[l60 - 1][i - 1] = 1.0 + uim * (cri * ccna - ccnb) / (cri * ccnc - ccnd);

    } // l160 loop

  } else { // label 165: nsh > 1

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

      int lpo = l180 + 1;

      int ltpo = lpo + l180;

      int lpt = lpo + 1;

      double dltpo = 1.0 * ltpo;

      dcomplex y1 = fbi[lpo - 1];

      dcomplex dy1 = (l180 * fbi[l180 - 1] - lpo * fbi[lpt - 1]) * c2->vkt[i - 1] / dltpo;

      dcomplex y2 = y1;

      dcomplex dy2 = dy1;

      int ic = 1;

      for (int ns = 2; ns <= nsh; ns++) {

	int nsmo = ns - 1;

	double vkr = vk * c1->rc[i - 1][nsmo - 1];

	if {ns % 2 != 0) {

	  int step = vk * (c1->rc[i - 1][ns - 1] - c1->rc[i - 1][nsmo - 1]) / nstp;

	  arg = c2->dc0[ic++];

	  rkc(nstp, step, arg, vkr, lpo, y1, y2, dy1, dy2);

	} else { // label 170

	  diel(nstpts, nsmo, i, ic, vk);

	  double stepts = vk * (c1->rc[i - 1][ns - 1] - c1->rc[i - 1][nsmo - 1]) / nstpts;

	  rkt(nstpts, stepts, vkr, lpo, y1, y2, dy1, dy2);

	}

      } // ns loop 176

      rmf[l180 - 1] = y1 * sz;

      drmf[l180 - 1] = dy1 * sz + y1;

      ref[l180 - 1] = y2 * sz;

      dref[l180 - 1] = dy2 * sz + y2;

    } // l180 loop

    dcomplex cri = (nsh % 2 == 0) ? 1.0 + I * 0.0 : c2->dc0[ic - 1] / ent;

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

      int lpo = l190 + 1;

      int ltpo = lpo + l190;

      int lpt = lpo + 1;

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

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

      dcomplex ccna = rmf[l190 - 1] * dfn;

      dcomplex ccnb = drmf[l190 - 1] * fn[lpo - 1] * sz * ltpo;

      dcomplex ccnc = rmf[l190 - 1] * dfb;

      dcomplex ccnd = drmf[l190 - 1] * fb[lpo - 1]* sz * ltpo;

      c1->rmi[l190 - 1][i -1] = 1.0 + uim * (ccna - ccnb) / (ccnc - ccnd);

      ccna = ref[l190 - 1] * dfn;

      ccnb = dref[l190 - 1] * fn[lpo] * sz * ltpo;

      ccnc = ref[l190 - 1] * dfb;

      ccnd = dref[l190 - 1] * fb[lpo - 1] * sz * ltpo;

      c1->rei[l190 - 1][i - 1] =1.0 + uim * (cri * ccna - ccnb) / (cri * ccnc - ccnd);

    } // l190 loop

  } // nsh if

  delete[] cfj;

  delete[] cfn;

  delete[] fb;

  delete[] fbi;

  delete[] fn;

  delete[] rfj;

  delete[] rfn;

  delete[] rmf;

  delete[] drmf;

  delete[] ref;

  delete[] dref;

}



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