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

Ensure that sequential execution of C++ inclusion mimics FORTRAN

parent 692d209b

src/include/Commons.h

+2 −2
Original line number Diff line number Diff line
@@ -358,8 +358,6 @@ public: 
  dcomplex tfsas;

  //! \brief Total scattering amplitude.

  dcomplex **tsas;

  //! \brief Total geometric cross-section.

  double gcs;

  //! \brief Total scattering cross-section.

  double scs;

  //! \brief Total extinction cross-section.
@@ -405,6 +403,8 @@ public: 
  const int& ndit = _ndit;

  //! \brief Read-only view of NDM.

  const int& ndm = _ndm;

  //! \brief Total geometric cross-section.

  double gcs;



  //! \brief TBD

  dcomplex *vh;

src/include/inclu_subs.h

+1 −1
Original line number Diff line number Diff line
@@ -93,6 +93,6 @@ void instr(ScattererConfiguration *sconf, ParticleDescriptor *c1); 
 * \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);

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/inclusion/inclusion.cpp

+33 −42
Original line number Diff line number Diff line
@@ -104,7 +104,6 @@ protected: 
public:

  int nimd;

  double extr;

  double gcs;

  

  //! \brief Pointer to a ParticleDescriptor structure.

  ParticleDescriptor *c1;
@@ -296,7 +295,7 @@ InclusionIterationData::InclusionIterationData(GeometryConfiguration *gconf, Sca 
  c1->rc[0][nimd - 1] = c1->ros[0] * sconf->get_rcf(0, nimd - 1);

  extr = c1->rc[0][nimd - 1];

  const double pig = acos(0.0) * 2.0;

  gcs = pig * extr * extr;

  c1->gcs = pig * extr * extr;

  

#ifdef USE_MAGMA

  proc_device = mpidata->rank % device_count;
@@ -448,6 +447,9 @@ InclusionIterationData::InclusionIterationData(const InclusionIterationData& rhs 
  xiblock = rhs.xiblock;

  number_of_scales = rhs.number_of_scales;



  nimd = rhs.nimd;

  extr = rhs.extr;

  

  proc_device = rhs.proc_device;

}
@@ -578,6 +580,9 @@ InclusionIterationData::InclusionIterationData(const mixMPI *mpidata, const int 
  firstxi = lastxi-xiblock+1;

  if (lastxi > number_of_scales) lastxi = number_of_scales;



  MPI_Bcast(&nimd, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&extr, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  

#ifdef USE_MAGMA

  proc_device = mpidata->rank % device_count;

#else
@@ -646,6 +651,9 @@ void InclusionIterationData::mpibcast(const mixMPI *mpidata) { 
  MPI_Bcast(&vk, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&xiblock, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&number_of_scales, 1, MPI_INT, 0, MPI_COMM_WORLD);



  MPI_Bcast(&nimd, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&extr, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

}

#endif
@@ -1122,7 +1130,7 @@ void inclusion(const string& config_file, const string& data_file, const string& 
	  //==============================================

	  if (myompthread == 0) {

	    // thread 0 writes its virtual files, now including contributions from all threads, to disk, and deletes them

	    p_outarray[0]->append_to_disk(output_path + "/c_OCLU");

	    p_outarray[0]->append_to_disk(output_path + "/c_OINCLU");

	    delete p_outarray[0];

	    vtppoanarray[0]->append_to_disk(output_path + "/c_TPPOAN");

	    delete vtppoanarray[0];
@@ -1134,7 +1142,7 @@ void inclusion(const string& config_file, const string& data_file, const string& 
		// get the data from process rr, creating a new virtual ascii file

		VirtualAsciiFile *p_output = new VirtualAsciiFile(mpidata, rr);

		// append to disk and delete virtual ascii file

		p_output->append_to_disk(output_path + "/c_OCLU");

		p_output->append_to_disk(output_path + "/c_OINCLU");

		delete p_output;

		// get the data from process rr, creating a new virtual binary file

		VirtualBinaryFile *vtppoanp = new VirtualBinaryFile(mpidata, rr);
@@ -1187,7 +1195,7 @@ void inclusion(const string& config_file, const string& data_file, const string& 
    elapsed = t_end - t_start;

    string message = "INFO: Calculation lasted " + to_string(elapsed.count()) + "s.\n";

    logger->log(message);

    logger->log("Finished: output written to " + output_path + "/c_OCLU\n");

    logger->log("Finished: output written to " + output_path + "/c_OINCLU\n");

    time_logger->log(message);

    fclose(timing_file);

    delete time_logger;
@@ -1366,7 +1374,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo 
  double sze = vkarg * cid->extr;

  last_configuration = 0;

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

    int iogi = cid->c1->iog[i133];

    int iogi = cid->c1->iog[i133 - 1];

    if (iogi != i133) {

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

	cid->c1->rmi[l123 - 1][i133 - 1] = cid->c1->rmi[l123 - 1][iogi - 1];
@@ -1376,7 +1384,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo 
      last_configuration++;

      int nsh = cid->c1->nshl[last_configuration - 1];

      int ici = (nsh + 1) / 2;

      if (i133 == 0) ici++;

      if (i133 == 1) ici++;

      if (idfc == 0) {

	for (int ic = 0; ic < ici; ic++)

	  cid->c1->dc0[ic] = sconf->get_dielectric_constant(ic, i133 - 1, jxi488 - 1);
@@ -1527,6 +1535,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo 
    } 

  } // i168 loop

  sprintf(virtual_line, "  EXT. SPHERE: SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n", sze, real(entn), imag(entn));

  output->append_line(virtual_line);

  // label 160

  double cs0 = 0.25 * cid->vk * cid->vk * cid->vk / acos(0.0);

  double csch = 2.0 * cid->vk * cid->sqsfi / cid->c1->gcs;
@@ -1671,7 +1680,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo 
		vtppoanp->append_line(VirtualBinaryLine(value));

	      }

	    }

	    sprintf(virtual_line, "     CLUSTER (ENSEMBLE AVERAGE, MODE%2d)\n", iavm);

	    sprintf(virtual_line, "     ENSEMBLE AVERAGE, MODE%2d\n", iavm);

	    output->append_line(virtual_line);

	    int jlr = 2;

	    for (int ilr210 = 1; ilr210 <= 2; ilr210++) {
@@ -1679,22 +1688,15 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo 
	      if (ilr210 == 2) jlr = 1;

	      double extsm = cid->c1->ecscm[ilr210 - 1];

	      double qextm = extsm * cid->sqsfi / cid->c1->gcs;

	      double extrm = extsm / cid->c1->ecs;

	      double scasm = cid->c1->scscm[ilr210 - 1];

	      double albdm = scasm / extsm;

	      double qscam = scasm * cid->sqsfi / cid->c1->gcs;

	      double scarm = scasm / cid->c1->scs;

	      double abssm = extsm - scasm;

	      double qabsm = abssm * cid->sqsfi / cid->c1->gcs;

	      double absrm = abssm / cid->c1->acs;

	      double acsecs = cid->c1->acs / cid->c1->ecs;

	      if (acsecs >= -1.0e-6 && acsecs <= 1.0e-6) absrm = 1.0;

	      dcomplex s0m = cid->c1->fsacm[ilr210 - 1][ilr210 - 1] * exri;

	      double qschum = imag(s0m) * csch;

	      double pschum = real(s0m) * csch;

	      double s0magm = cabs(s0m) * cs0;

	      double rfinrm = real(cid->c1->fsacm[ilr210 - 1][ilr210 - 1]) / real(cid->c1->tfsas);

	      double extcrm = imag(cid->c1->fsacm[ilr210 - 1][ilr210 - 1]) / imag(cid->c1->tfsas);

	      if (inpol == 0) {

		sprintf(virtual_line, "   LIN %2d\n", ipol);

		output->append_line(virtual_line);
@@ -1703,33 +1705,24 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo 
		output->append_line(virtual_line);

	      }

	      // label 208

	      sprintf(virtual_line, " ----- SCC ----- ABC ----- EXC ----- ALBEDC --\n");

	      sprintf(virtual_line, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");

	      output->append_line(virtual_line);

	      sprintf(virtual_line, " %14.7lE%15.7lE%15.7lE%15.7lE\n", scasm, abssm, extsm, albdm);

	      output->append_line(virtual_line);

	      double alamb = 2.0 * 3.141592653589793238 / cid->vk;

	      sprintf(virtual_line, "INSERTION: SCASECM %5d%15.7E%15.7E%15.7E%15.7E", ipol, alamb, scasm, abssm, extsm);

	      sprintf(virtual_line, "INSERTION: SCASECM %5d%15.7E%15.7E%15.7E%15.7E\n", ipol, alamb, scasm, abssm, extsm);

	      output->append_line(virtual_line);

	      sprintf(virtual_line, " --- SCC/TGS - ABC/TGS - EXC/TGS ---\n");

	      sprintf(virtual_line, " ---- SCS/GS -- ABC/GS -- EXS/GS ---\n");

	      output->append_line(virtual_line);

	      sprintf(virtual_line, " %14.7lE%15.7lE%15.7lE\n", qscam, qabsm, qextm);

	      output->append_line(virtual_line);

	      sprintf(virtual_line, " ---- SCCRT --- ABCRT --- EXCRT ----\n");

	      output->append_line(virtual_line);

	      sprintf(virtual_line, " %14.7lE%15.7lE%15.7lE\n", scarm, absrm, extrm);

	      output->append_line(virtual_line);

	      sprintf(

		      virtual_line, "  FSAC(%1d,%1d)=%15.7lE%15.7lE   FSAC(%1d,%1d)=%15.7lE%15.7lE\n",

		      virtual_line, "  FSAS(%1d,%1d)=%15.7lE%15.7lE   FSAS(%1d,%1d)=%15.7lE%15.7lE\n",

		      ilr210, ilr210, real(cid->c1->fsacm[ilr210 - 1][ilr210 - 1]),

		      imag(cid->c1->fsacm[ilr210 - 1][ilr210 - 1]), jlr, ilr210,

		      real(cid->c1->fsacm[jlr - 1][ilr210 - 1]), imag(cid->c1->fsacm[jlr - 1][ilr210 - 1])

		      );

	      output->append_line(virtual_line);

	      sprintf(

		      virtual_line, "  RE(FSAC(%1d,%1d))/RE(TFSAS)=%15.7lE, IM(FSAC(%1d,%1d))/IM(TFSAS)=%15.7lE\n",

		      ilr210, ilr210, rfinrm, ilr210, ilr210, extcrm

		      );

	      output->append_line(virtual_line);

	      sprintf(virtual_line, "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschum, pschum, s0magm);

	      output->append_line(virtual_line);

	      double rapr = cid->c1->ecscm[ilr210 - 1] - cid->gapm[2][ilr210 - 1];
@@ -1739,16 +1732,12 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo 
	      output->append_line(virtual_line);

	      sprintf(virtual_line, "  Fk=%15.7lE\n", fz);

	      output->append_line(virtual_line);

	      if (ilr210 == 1) {

		sprintf(virtual_line, "INSERTION: CSM_CLUSTER  %15.7lE%15.7lE%15.7lE%15.7lE\n", alamb, scasm, abssm, extsm);

		output->append_line(virtual_line);

	      }

	    } // ilr210 loop

	    double rmbrif = (real(cid->c1->fsacm[0][0]) - real(cid->c1->fsacm[1][1])) / real(cid->c1->fsacm[0][0]);

	    double rmdchr = (imag(cid->c1->fsacm[0][0]) - imag(cid->c1->fsacm[1][1])) / imag(cid->c1->fsacm[0][0]);

	    sprintf(virtual_line, "  (RE(FSAC(1,1))-RE(FSAC(2,2)))/RE(FSAC(1,1))=%15.7lE\n", rmbrif);

	    sprintf(virtual_line, "  (RE(FSAS(1,1))-RE(FSAS(2,2)))/RE(FSAS(1,1))=%15.7lE\n", rmbrif);

	    output->append_line(virtual_line);

	    sprintf(virtual_line, "  (IM(FSAC(1,1))-IM(FSAC(2,2)))/IM(FSAC(1,1))=%15.7lE\n", rmdchr);

	    sprintf(virtual_line, "  (IM(FSAS(1,1))-IM(FSAS(2,2)))/IM(FSAS(1,1))=%15.7lE\n", rmdchr);

	    output->append_line(virtual_line);

	  }

	  // label 212
@@ -1850,6 +1839,8 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo 
	      vtppoanp->append_line(VirtualBinaryLine(value));

	    }

	  }

	  sprintf(virtual_line, "     SINGLE SCATTERER\n");

	  output->append_line(virtual_line);

	  int jlr = 2;

	  for (int ilr290 = 1; ilr290 <= 2; ilr290++) {

	    int ipol = (ilr290 % 2 == 0) ? 1 : -1;
@@ -1873,7 +1864,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo 
	      output->append_line(virtual_line);

	    }

	    // label 275

	    sprintf(virtual_line, " ----- SCC ----- ABC ----- EXC ----- ALBEDC --\n");

	    sprintf(virtual_line, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");

	    output->append_line(virtual_line);

	    sprintf(

		    virtual_line, " %14.7lE%15.7lE%15.7lE%15.7lE\n",
@@ -1882,7 +1873,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo 
	    output->append_line(virtual_line);

	    double alam = 2.0 * 3.141592653589793238 / cid->vk;

	    sprintf(virtual_line, "INSERTION: SCASEC %5d%14.7lE%14.7lE%14.7lE%14.7lE\n", ipol, alam, scasec, abssec, extsec);

	    sprintf(virtual_line, " --- SCC/TGS - ABC/TGS - EXC/TGS ---\n");

	    sprintf(virtual_line, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n");

	    output->append_line(virtual_line);

	    sprintf(

		    virtual_line, " %14.7lE%15.7lE%15.7lE\n",
@@ -1895,7 +1886,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo 
		    ilr290, ilr290, real(cid->c1->fsac[ilr290 - 1][ilr290 - 1]),

		    imag(cid->c1->fsac[ilr290 - 1][ilr290 - 1]), jlr, ilr290,

		    real(cid->c1->fsac[jlr - 1][ilr290 - 1]),

		    imag(cid->c1->fsac[jlr][ilr290 - 1])

		    imag(cid->c1->fsac[jlr - 1][ilr290 - 1])

		    );

	    output->append_line(virtual_line);

	    sprintf(
@@ -1904,7 +1895,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo 
		    ilr290, ilr290, real(cid->c1->sac[ilr290 - 1][ilr290 - 1]),

		    imag(cid->c1->sac[ilr290 - 1][ilr290 - 1]), jlr, ilr290,

		    real(cid->c1->sac[jlr - 1][ilr290 - 1]),

		    imag(cid->c1->sac[jlr][ilr290 - 1])

		    imag(cid->c1->sac[jlr - 1][ilr290 - 1])

		    );

	    output->append_line(virtual_line);

	    sprintf(
@@ -1953,11 +1944,11 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo 
	  } //ilr290 loop

	  double rbirif = (real(cid->c1->fsac[0][0]) - real(cid->c1->fsac[1][1])) / real(cid->c1->fsac[0][0]);

	  double rdichr = (imag(cid->c1->fsac[0][0]) - imag(cid->c1->fsac[1][1])) / imag(cid->c1->fsac[0][0]);

	  sprintf(virtual_line, "  (RE(FSAC(1,1))-RE(FSAC(2,2)))/RE(FSAC(1,1))=%15.7lE\n", rbirif);

	  sprintf(virtual_line, "  (RE(FSAS(1,1))-RE(FSAS(2,2)))/RE(FSAS(1,1))=%15.7lE\n", rbirif);

	  output->append_line(virtual_line);

	  sprintf(virtual_line, "  (IM(FSAC(1,1))-IM(FSAC(2,2)))/IM(FSAC(1,1))=%15.7lE\n", rdichr);

	  sprintf(virtual_line, "  (IM(FSAS(1,1))-IM(FSAS(2,2)))/IM(FSAS(1,1))=%15.7lE\n", rdichr);

	  output->append_line(virtual_line);

	  sprintf(virtual_line, "  MULC\n");

	  sprintf(virtual_line, "  MULL\n");

	  output->append_line(virtual_line);

	  for (int i = 0; i < 4; i++) {

	    sprintf(
@@ -1966,7 +1957,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo 
		    );

	    output->append_line(virtual_line);

	  }

	  sprintf(virtual_line, "  MULCLR\n");

	  sprintf(virtual_line, "  MULLLR\n");

	  output->append_line(virtual_line);

	  for (int i = 0; i < 4; i++) {

	    sprintf(

src/libnptm/Commons.cpp

+8 −3
Original line number Diff line number Diff line
@@ -601,6 +601,7 @@ ParticleDescriptor::ParticleDescriptor(GeometryConfiguration *gconf, ScattererCo 
  _ndi = 0;

  _ndit = 0;

  _ndm = 0;

  gcs = 0.0;

  _num_configurations = sconf->configurations;

  _num_layers = (sconf->use_external_sphere) ? 1 : 0;

  for (int nli = 0; nli < num_configurations; nli++) _num_layers += sconf->get_nshl(nli);
@@ -737,6 +738,7 @@ ParticleDescriptor::ParticleDescriptor(const ParticleDescriptor &rhs) { 
  _ndi = rhs._ndi;

  _ndit = rhs._ndit;

  _ndm = rhs._ndm;

  gcs = rhs.gcs;

  _num_configurations = rhs._num_configurations;

  _num_layers = rhs._num_layers;
@@ -887,6 +889,7 @@ ParticleDescriptor::ParticleDescriptor(const mixMPI *mpidata) { 
  MPI_Bcast(&_ndi, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&_ndit, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&_ndm, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&gcs, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&_num_configurations, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&_num_layers, 1, MPI_INT, 0, MPI_COMM_WORLD);
@@ -965,6 +968,7 @@ void ParticleDescriptor::mpibcast(const mixMPI *mpidata) { 
  MPI_Bcast(&_ndi, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&_ndit, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&_ndm, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&gcs, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&_num_configurations, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&_num_layers, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(vec_rmi, _nsph * _li, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
@@ -1007,7 +1011,6 @@ void ParticleDescriptor::mpibcast(const mixMPI *mpidata) { 
    MPI_Bcast(vintt, 16, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

    MPI_Bcast(vec_tsas, 4, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

    MPI_Bcast(&tfsas, 1, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

    MPI_Bcast(&gcs, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

    MPI_Bcast(&scs, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

    MPI_Bcast(&ecs, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

    MPI_Bcast(&acs, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
@@ -1415,7 +1418,6 @@ ParticleDescriptorCluster::ParticleDescriptorCluster(const mixMPI *mpidata) : Pa 
  tsas[0] = vec_tsas;

  tsas[1] = vec_tsas + 2;

  MPI_Bcast(&tfsas, 1, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  MPI_Bcast(&gcs, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&scs, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&ecs, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&acs, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
@@ -1507,9 +1509,12 @@ ParticleDescriptorInclusion::ParticleDescriptorInclusion(GeometryConfiguration * 
  _ncou = _nsph * _nsph - 1;

  _litpo = _li + _li + 1;

  _litpos = _litpo * _litpo;

  _lmpo = _lm + 1;

  _lmtpo = _li + _le + 1;

  _lmtpos = _lmtpo * _lmtpo;

  _nv3j = (_lm * (_lm + 1) * (2 * _lm + 7)) / 6;

  _ndi = _nsph * _nlim;

  _ndit = 2 * _nsph * _nlim;

  vec_am0m = new dcomplex[_nlemt * _nlemt]();

  vec_fsac = new dcomplex[4]();

  vec_sac = new dcomplex[4]();
@@ -1555,7 +1560,7 @@ ParticleDescriptorInclusion::ParticleDescriptorInclusion(GeometryConfiguration * 
  te0 = new dcomplex[_le]();

  vec_at = new dcomplex[_nlemt * _ndm]();

  at = new dcomplex*[_nlemt];

  for (int ai = 0; ai < _nlemt; ai++) at[ai] = vec_at + (ai * _nlemt);

  for (int ai = 0; ai < _nlemt; ai++) at[ai] = vec_at + (ai * _ndm);

}



ParticleDescriptorInclusion::ParticleDescriptorInclusion(const ParticleDescriptorInclusion &rhs) : ParticleDescriptor(rhs) {

src/libnptm/inclu_subs.cpp

+28 −30
Original line number Diff line number Diff line
@@ -92,7 +92,7 @@ void exma(dcomplex **am, ParticleDescriptor *c1) { 
void incms(dcomplex **am, double enti, ParticleDescriptor *c1) {

  const dcomplex cc0 = 0.0 + I * 0.0;

  dcomplex **at = c1->at;

  int nbl;

  int nbl, i1;

  const int ndi = c1->ndi;

  const int ndit = ndi + ndi;

  const int ndm = c1->ndm;
@@ -111,7 +111,7 @@ void incms(dcomplex **am, double enti, ParticleDescriptor *c1) { 
	  int m1 = im1 - l1po;

	  int ilm1 = il1 + m1;

	  int ilm1e = ilm1 + ndi;

	  int i1 = in1 + ilm1;

	  i1 = in1 + ilm1;

	  int i1e = in1 + ilm1e;

	  int j1 = in2 + ilm1;

	  int j1e = in2 + ilm1e;
@@ -156,7 +156,7 @@ void incms(dcomplex **am, double enti, ParticleDescriptor *c1) { 
      for (int im1 = 1; im1 <= l1tpo; im1++) {

	int m1 = im1 - l1po;

	int ilm1 = il1 + m1;

	int i1 = in1 + ilm1;

	i1 = in1 + ilm1;

	int i1e = i1 + ndi;

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

	  int i2 = in1 + ilm2;
@@ -172,15 +172,14 @@ void incms(dcomplex **am, double enti, ParticleDescriptor *c1) { 
    } // l1 loop 30

  } // n1 loop 30

  int nditpo = ndit + 1;

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

  for (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;

  i1 = 0;

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

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

    dcomplex fre = c1->re0[l1 - 1];
@@ -198,7 +197,6 @@ void incms(dcomplex **am, double enti, ParticleDescriptor *c1) { 
      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];
@@ -228,7 +226,7 @@ void incms(dcomplex **am, double enti, ParticleDescriptor *c1) { 
	      int m1 = im1 - l1po;

	      int ilm1 = il1 + m1;

	      int jlm1 = il1 - m1;

	      int i1 = in1 + ilm1;

	      i1 = in1 + ilm1;

	      int i1e = i1 + ndi;

	      int j1 = in1 + jlm1;

	      int j1e = j1 + ndi;
@@ -256,20 +254,20 @@ void incms(dcomplex **am, double enti, ParticleDescriptor *c1) { 
  } else {

    // label 55

    int i2 = 0;

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

    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;

      int m2 = -l2 - 1;

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

	m2++;

	i2++;

	int i2e = i2 + c1->nlem;

	int i3 = i2 + ndit;

	int i3e = i3 + c1->nlem;

	int i1 = 0;

	i1 = 0;

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

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

	    int l1tpo = l1 + l1 + 1;
@@ -355,7 +353,7 @@ void indme( 
      return;

    }

    // label 128

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

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

  }

  // label 132

  dcomplex aris = sz * entn;
@@ -557,7 +555,7 @@ void instr(ScattererConfiguration *sconf, ParticleDescriptor *c1) { 
  delete[] ylm;

}



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

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;