Commit a2fd1b1e authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Add C++ porting of subroutine OSPV()

parent 67e29b9a

src/include/inclu_subs.h

+15 −0
Original line number Diff line number Diff line
@@ -80,4 +80,19 @@ void indme( 
 * \param c1: `ParticleDescriptor *` Pointer to a ParticleDescriptor instance.

 */

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



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

 *

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

 * \param vk: `double` TBD.

 * \param sze: `double` TBD.

 * \param exri: `double` External medium refractive index.

 * \param entn: `dcomplex` Outer sphere refractive index.

 * \param enti: `double` Imaginary part of the outer medium refractive index.

 * \param jer: `int &` Reference to an integer error flag.

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

 * \param arg: `dcomplex` Complex Bessel function argument.

 */

void ospv(ParticleDescriptor *c1, double vk, double sze, double exri, dcomplex entn, double enti, int &jer, int &lcalc, dcomplex arg);



#endif // INCLUDE_INCLU_SUBS_H_

src/libnptm/inclu_subs.cpp

+223 −0
Original line number Diff line number Diff line
@@ -556,3 +556,226 @@ void instr(ScattererConfiguration *sconf, ParticleDescriptor *c1) { 
  } // nf50 loop

  delete[] ylm;

}



void ospv(ParticleDescriptor *c1, double vk, double sze, double exri, dcomplex entn, double enti, int &jer, int &lcalc, dcomplex arg) {

  const dcomplex uim = 0.0 + I * 1.0;

  const int nsph = c1->nsph;

  const int nsphmo = c1->nsph - 1;

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

  const int litpo = lit + 1;

  const int array_size = (c1->litpo > c1->lmtpo) ? c1->litpo : c1->lmtpo;

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

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

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

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

  

  int ivhb = 0;

  for (int i130 = 1; i130 <= nsphmo; i130++) {

    int ipo = i130 + 1;

    for (int j130 = ipo; j130 <= nsph; j130++) {

      double rx = c1->rxx[j130 - 1] - c1->rxx[i130 - 1];

      double ry = c1->ryy[j130 - 1] - c1->ryy[i130 - 1];

      double rz = c1->rzz[j130 - 1] - c1->rzz[i130 - 1];

      double rr = sqrt(rx * rx + ry * ry + rz * rz);

      arg = rr * vk * entn;

      if (enti != 0.0) {

	cbf(lit, arg, lcalc, cfj);

	if (lcalc < lit) {

	  jer = 11;

	  delete[] cfj;

	  delete[] cfn;

	  delete[] rfj;

	  delete[] rfn;

	  return;

	}

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

	for (int lpo = 0; lpo < litpo; lpo++) c1->vh[lpo + ivhb] = cfj[lpo] + uim * cfn[lpo];

	// goes to 130

      } else { // label 123

	double rarg = real(arg);

	rbf(lit, rarg, lcalc, rfj);

	if (lcalc < lit) {

	  jer = 11;

	  delete[] cfj;

	  delete[] cfn;

	  delete[] rfj;

	  delete[] rfn;

	  return;

	}

	rnf(lit, rarg, lcalc, rfn);

	if (lcalc < lit) {

	  jer = 12;

	  delete[] cfj;

	  delete[] cfn;

	  delete[] rfj;

	  delete[] rfn;

	  return;

	}

	for (int lpo = 0; lpo < litpo; lpo++) c1->vh[lpo + ivhb] = rfj[lpo] + uim * rfn[lpo];

      }

      // label 130

      ivhb += litpo;

    } // j130 loop

  } // i130 loop

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

  const int lmtpo = lmt + 1;

  ivhb = 0;

  for (int i155 = 1; i155 <= nsph; i155++) {

    double rx = c1->rxx[i155 - 1];

    double ry = c1->ryy[i155 - 1];

    double rz = c1->rzz[i155 - 1];

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

      double rr = sqrt(rx * rx + ry * ry + rz * rz);

      arg = rr * vk * entn;

      if (enti != 0.0) {

	cbf(lmt, arg, lcalc, cfj);

	if (lcalc < lmt) {

	  jer = 11;

	  delete[] cfj;

	  delete[] cfn;

	  delete[] rfj;

	  delete[] rfn;

	  return;

	}

	for (int lpo = 0; lpo < lmtpo; lpo++) c1->vj0[lpo + ivhb] = cfj[lpo];

	// goes to 155

      } else { // label 150

	double rarg = real(arg);

	rbf(lmt, rarg, lcalc, rfj);

	if (lcalc < lmt) {

	  jer = 11;

	  delete[] cfj;

	  delete[] cfn;

	  delete[] rfj;

	  delete[] rfn;

	  return;

	}

	for (int lpo = 0; lpo < lmtpo; lpo++) c1->vj0[lpo + ivhb] = rfj[lpo];

      }

    }

    // label 155

    ivhb += lmtpo;

  } // i155 loop

  

  const int lepo = c1->le + 1;

  const int lept = c1->le + 2;

  dcomplex *fb0 = new dcomplex[lept]();

  dcomplex *fh0 = new dcomplex[lept]();

  dcomplex aris0 = sze * entn;

  arg = aris0;

  if (enti != 0.0) {

    cbf(lepo, arg, lcalc, cfj);

    if (lcalc < lepo) {

      jer = 11;

      delete[] cfj;

      delete[] cfn;

      delete[] rfj;

      delete[] rfn;

      delete[] fb0;

      delete[] fh0;

      return;

    }

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

    for (int j162 = 0; j162 < lept; j162++) {

      fb0[j162] = cfj[j162];

      fh0[j162] = cfj[j162] + uim * cfn[j162];

    } // j162 loop

    // goes to 170

  } else { // label 163

    double rarg = real(arg);

    rbf(lepo, rarg, lcalc, rfj);

    if (lcalc < lepo) {

      jer = 11;

      delete[] cfj;

      delete[] cfn;

      delete[] rfj;

      delete[] rfn;

      delete[] fb0;

      delete[] fh0;

      return;

    }

    rnf(lepo, rarg, lcalc, rfn);

    if (lcalc < lepo) {

      jer = 12;

      delete[] cfj;

      delete[] cfn;

      delete[] rfj;

      delete[] rfn;

      delete[] fb0;

      delete[] fh0;

      return;

    }

    for (int j168 = 0; j168 < lept; j168++) {

      fb0[j168] = rfj[j168];

      fh0[j168] = rfj[j168] + uim * rfn[j168];

    } // j168 loop

  }

  // label 170

  double arex = sze * exri;

  arg = arex;

  double rarg = arex;

  rbf(lepo, rarg, lcalc, rfj);

  if (lcalc < lepo) {

    jer = 1;

    delete[] cfj;

    delete[] cfn;

    delete[] rfj;

    delete[] rfn;

    delete[] fb0;

    delete[] fh0;

    return;

  }

  rnf(lepo, rarg, lcalc, rfn);

  if (lcalc < lepo) {

    jer = 2;

    delete[] cfj;

    delete[] cfn;

    delete[] rfj;

    delete[] rfn;

    delete[] fb0;

    delete[] fh0;

    return;

  }

  dcomplex *fbe = new dcomplex[lept]();

  dcomplex *fhe = new dcomplex[lept]();

  for (int j175 = 0; j175 < lept; j175++) {

    fbe[j175] = rfj[j175];

    fhe[j175] = rfj[j175] + uim * rfn[j175];

  } // j175 loop

  dcomplex cri = exri / entn;

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

    int lpo = l184 + 1;

    int lpt = lpo + 1;

    double dltpo = 1.0 / (lpo + l184);

    dcomplex dfb0 = fb0[lpo - 1] + (l184 * fb0[l184 - 1] - lpo * fb0[lpt - 1]) * aris0 * dltpo;

    dcomplex dfh0 = fh0[lpo - 1] + (l184 * fh0[l184 - 1] - lpo * fh0[lpt - 1]) * aris0 * dltpo;

    dcomplex dfbe = fbe[lpo - 1] + (l184 * fbe[l184 - 1] - lpo * fbe[lpt - 1]) * arex * dltpo;

    dcomplex dfhe = fhe[lpo - 1] + (l184 * fhe[l184 - 1] - lpo * fhe[lpt - 1]) * arex * dltpo;

    dcomplex ccna = aris0 * fb0[lpo - 1] * dfhe;

    dcomplex ccnb = arex * fhe[lpo - 1] * dfb0;

    c1->rm0[l184 - 1] = -(ccna * cri - ccnb) * uim;

    c1->re0[l184 - 1] = -(ccna - ccnb * cri) * uim;

    ccna = aris0 * fh0[lpo - 1] * dfhe;

    ccnb = arex * fhe[lpo - 1] * dfh0;

    c1->rmw[l184 - 1] = -(ccna * cri - ccnb) * uim;

    c1->rew[l184 - 1] = -(ccna - ccnb * cri) * uim;

    ccna = aris0 * fh0[lpo - 1] * dfbe;

    ccnb = arex * fbe[lpo - 1] * dfh0;

    c1->tm[l184 - 1] = (ccna * cri - ccnb) * uim;

    c1->te[l184 - 1] = (ccna - ccnb * cri) * uim;

    ccna = aris0 * fb0[lpo - 1] * dfbe;

    ccnb = arex * fbe[lpo - 1] * dfb0;

    c1->tm0[l184 - 1] = (ccna * cri - ccnb) * uim;

    c1->te0[l184 - 1] = (ccna - ccnb * cri) * uim;

  } // l184 loop

  

  // Clean up memory.

  delete[] cfj;

  delete[] cfn;

  delete[] rfj;

  delete[] rfn;

  delete[] fb0;

  delete[] fh0;

  delete[] fbe;

  delete[] fhe;

}