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

Move C3 data structures under ParticleDescriptor

parent 821dd790

src/cluster/cluster.cpp

+28 −28
Original line number Diff line number Diff line
@@ -270,7 +270,7 @@ void cluster(const string& config_file, const string& data_file, const string& o 
      p_output->append_line(virtual_line);

      sprintf(virtual_line, " \n");

      p_output->append_line(virtual_line);

      str(sconf, cid->c1, cid->c3, cid->c6);

      str(sconf, cid->c1, cid->c6);

      thdps(cid->c1->lm, cid->zpv);

      double exdc = sconf->exdc;

      double exri = sqrt(exdc);
@@ -895,7 +895,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
  double cs0 = 0.25 * cid->vk * cid->vk * cid->vk / acos(0.0);

  double csch = 0.0, qschu = 0.0, pschu = 0.0, s0mag = 0.0;

  dcomplex s0 = 0.0 + 0.0 * I;

  scr0(cid->vk, exri, cid->c1, cid->c3);

  scr0(cid->vk, exri, cid->c1);

  double sqk = cid->vk * cid->vk * exdc;

  aps(cid->zpv, cid->c1->li, nsph, cid->c1, sqk, cid->gaps);

  rabas(inpol, cid->c1->li, nsph, cid->c1, cid->tqse, cid->tqspe, cid->tqss, cid->tqsps);
@@ -952,10 +952,10 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
      output->append_line(virtual_line);

    }

  } // i170 loop

  sprintf(virtual_line, "  FSAT=%15.7lE%15.7lE\n", real(cid->c3->tfsas), imag(cid->c3->tfsas));

  sprintf(virtual_line, "  FSAT=%15.7lE%15.7lE\n", real(cid->c1->tfsas), imag(cid->c1->tfsas));

  output->append_line(virtual_line);

  csch = 2.0 * cid->vk * cid->sqsfi / cid->c3->gcs;

  s0 = cid->c3->tfsas * exri;

  csch = 2.0 * cid->vk * cid->sqsfi / cid->c1->gcs;

  s0 = cid->c1->tfsas * exri;

  qschu = imag(s0) * csch;

  pschu = real(s0) * csch;

  s0mag = cabs(s0) * cs0;
@@ -1124,23 +1124,23 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
	      int ipol = (ilr210 % 2 == 0) ? 1 : -1;

	      if (ilr210 == 2) jlr = 1;

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

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

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

	      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->c3->gcs;

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

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

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

	      double abssm = extsm - scasm;

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

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

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

	      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->c3->tfsas);

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

	      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);
@@ -1224,7 +1224,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
	    output->append_line(virtual_line);

	  }

	  // label 224

	  scr2(cid->vk, vkarg, exri, cid->duk, cid->c1, cid->c3);

	  scr2(cid->vk, vkarg, exri, cid->duk, cid->c1);

	  if (cid->c1->li != cid->c1->le) {

	    sprintf(virtual_line, "     SPHERES; LMX=MIN0(LI,LE)\n");

	    output->append_line(virtual_line);
@@ -1271,14 +1271,14 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
	  } // i226 loop

	  sprintf(

		  virtual_line, "  SAT(1,1)=%15.7lE%15.7lE, SAT(2,1)=%15.7lE%15.7lE\n",

		  real(cid->c3->tsas[0][0]), imag(cid->c3->tsas[0][0]),

		  real(cid->c3->tsas[1][0]), imag(cid->c3->tsas[1][0])

		  real(cid->c1->tsas[0][0]), imag(cid->c1->tsas[0][0]),

		  real(cid->c1->tsas[1][0]), imag(cid->c1->tsas[1][0])

		  );

	  output->append_line(virtual_line);

	  sprintf(

		  virtual_line, "  SAT(1,2)=%15.7lE%15.7lE, SAT(2,2)=%15.7lE%15.7lE\n",

		  real(cid->c3->tsas[0][1]), imag(cid->c3->tsas[0][1]),

		  real(cid->c3->tsas[1][1]), imag(cid->c3->tsas[1][1])

		  real(cid->c1->tsas[0][1]), imag(cid->c1->tsas[0][1]),

		  real(cid->c1->tsas[1][1]), imag(cid->c1->tsas[1][1])

		  );

	  output->append_line(virtual_line);

	  sprintf(virtual_line, "     CLUSTER\n");
@@ -1388,23 +1388,23 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
	    int ipol = (ilr290 % 2 == 0) ? 1 : -1;

	    if (ilr290 == 2) jlr = 1;

	    double extsec = cid->c1->ecsc[ilr290 - 1];

	    double qext = extsec * cid->sqsfi / cid->c3->gcs;

	    double extrat = extsec / cid->c3->ecs;

	    double qext = extsec * cid->sqsfi / cid->c1->gcs;

	    double extrat = extsec / cid->c1->ecs;

	    double scasec = cid->c1->scsc[ilr290 - 1];

	    double albedc = scasec / extsec;

	    double qsca = scasec * cid->sqsfi / cid->c3->gcs;

	    double scarat = scasec / cid->c3->scs;

	    double qsca = scasec * cid->sqsfi / cid->c1->gcs;

	    double scarat = scasec / cid->c1->scs;

	    double abssec = extsec - scasec;

	    double qabs = abssec * cid->sqsfi / cid->c3->gcs;

	    double qabs = abssec * cid->sqsfi / cid->c1->gcs;

	    double absrat = 1.0;

	    double ratio = cid->c3->acs / cid->c3->ecs;

	    if (ratio < -1.0e-6 || ratio > 1.0e-6) absrat = abssec / cid->c3->acs;

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

	    if (ratio < -1.0e-6 || ratio > 1.0e-6) absrat = abssec / cid->c1->acs;

	    s0 = cid->c1->fsac[ilr290 - 1][ilr290 - 1] * exri;

	    double qschu = imag(s0) * csch;

	    double pschu = real(s0) * csch;

	    s0mag = cabs(s0) * cs0;

	    double refinr = real(cid->c1->fsac[ilr290 - 1][ilr290 - 1]) / real(cid->c3->tfsas);

	    double extcor = imag(cid->c1->fsac[ilr290 - 1][ilr290 - 1]) / imag(cid->c3->tfsas);

	    double refinr = real(cid->c1->fsac[ilr290 - 1][ilr290 - 1]) / real(cid->c1->tfsas);

	    double extcor = imag(cid->c1->fsac[ilr290 - 1][ilr290 - 1]) / imag(cid->c1->tfsas);

	    if (inpol == 0) {

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

	      output->append_line(virtual_line);

src/include/Commons.h

+65 −47
Original line number Diff line number Diff line
@@ -40,50 +40,53 @@ 
class ParticleDescriptor;

class mixMPI;



/*! \brief Representation of the FORTRAN C3 blocks.

 */

class C3 {

public:

  //! \brief QUESTION: definition?

  dcomplex tfsas;

  //! \brief QUESTION: definition?

  dcomplex **tsas;

  //! \brief Total geometric cross-section.

  double gcs;

  //! \brief Total scattering cross-section.

  double scs;

  //! \brief Total extinction cross-section.

  double ecs;

  //! \brief Total absorption cross-section.

  double acs;



  /*! \brief C3 instance constructor.

   */

  C3();



  /*! \brief C3 instance constructor copying its contents from a preexisting object.

   */

  C3(const C3& rhs);



  /*! \brief C3 instance destroyer.

   */

  ~C3();



#ifdef MPI_VERSION

  /*! \brief C3 instance constructor copying all contents off MPI broadcast from MPI process 0

   *

   * \param mpidata: `mixMPI *` pointer to MPI data structure.

   */

  C3(const mixMPI *mpidata);



  /*! \brief send C3 instance from MPI process 0 via MPI broadcasts to all other processes

   *

   * \param mpidata: `mixMPI *` pointer to MPI data structure.

   */

  void mpibcast(const mixMPI *mpidata);

#endif



};

// /*! \brief Representation of the FORTRAN C3 blocks.

//  */

// class C3 {

// protected:

//   dcomplex *vec_tsas;

  

// public:

//   //! \brief Total forward scattering amplitude.

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

//   double ecs;

//   //! \brief Total absorption cross-section.

//   double acs;



//   /*! \brief C3 instance constructor.

//    */

//   C3();



//   /*! \brief C3 instance constructor copying its contents from a preexisting object.

//    */

//   C3(const C3& rhs);



//   /*! \brief C3 instance destroyer.

//    */

//   ~C3();



// #ifdef MPI_VERSION

//   /*! \brief C3 instance constructor copying all contents off MPI broadcast from MPI process 0

//    *

//    * \param mpidata: `mixMPI *` pointer to MPI data structure.

//    */

//   C3(const mixMPI *mpidata);



//   /*! \brief send C3 instance from MPI process 0 via MPI broadcasts to all other processes

//    *

//    * \param mpidata: `mixMPI *` pointer to MPI data structure.

//    */

//   void mpibcast(const mixMPI *mpidata);

// #endif



// };



/*! \brief Representation of the FORTRAN C6 blocks.

 */
@@ -218,10 +221,8 @@ public: 
 */

class ClusterIterationData {

public:

  //! \brief Pointer to a C1 structure.

  //! \brief Pointer to a ParticleDescriptor structure.

  ParticleDescriptor *c1;

  //! \brief Pointer to a C3 structure.

  C3 *c3;

  //! \brief Pointer to a C6 structure.

  C6 *c6;

  //! \brief Pointer to a C9 structure.
@@ -354,6 +355,10 @@ protected: 
  dcomplex *vec_vints;

  // >>> END OF SECTION NEEDED BY SPHERE AND CLUSTER <<< //



  // >>> NEEDED BY CLUSTER <<< //

  dcomplex *vec_tsas;

  // >>> END OF SECTION NEEDED BY CLUSTER <<< //

  

  // >>> NEEDED BY CLUSTER AND INCLU <<<

  //! \brief Maximum external field expansion order.

  int _le;
@@ -478,7 +483,20 @@ public: 
  // >>> END OF SECTION NEEDED BY SPHERE AND CLUSTER <<< //



  // >>> NEEDED BY CLUSTER <<<

  // \brief Vector of field intensity components.

  dcomplex *vintt;

  //! \brief Total forward scattering amplitude.

  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.

  double ecs;

  //! \brief Total absorption cross-section.

  double acs;

  // >>> END OF SECTION NEEDED BY CLUSTER <<< //

  

  // >>> NEEDED BY CLUSTER AND INCLU <<<

src/include/clu_subs.h

+3 −7
Original line number Diff line number Diff line
@@ -348,9 +348,8 @@ void rftr( 
 * \param vk: `double` Wave number

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

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

 * \param c3: `C3 *` Pointer to a C3 instance.

 */

void scr0(double vk, double exri, ParticleDescriptor *c1, C3 *c3);

void scr0(double vk, double exri, ParticleDescriptor *c1);



/*! \brief Compute the scattering amplitude for a single sphere in an aggregate.

 *
@@ -362,11 +361,9 @@ void scr0(double vk, double exri, ParticleDescriptor *c1, C3 *c3); 
 * \param exri: `double` External medium refractive index.

 * \param duk: `double *` QUESTION: definition?

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

 * \param c3: `C3 *` Pointer to a C3 instance.

 */

void scr2(

	  double vk, double vkarg, double exri, double *duk, ParticleDescriptor *c1,

	  C3 *c3

	  double vk, double vkarg, double exri, double *duk, ParticleDescriptor *c1

);



/*! \brief Transform sphere Cartesian coordinates to spherical coordinates.
@@ -377,10 +374,9 @@ void scr2( 
 *

 * \param sconf: `ScattererConfiguration *` Pointer to scatterer configuration object.

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

 * \param c3: `C3 *` Pointer to a C3 instance.

 * \param c6: `C6 *` Pointer to a C6 instance.

 */

void str(ScattererConfiguration *sconf, ParticleDescriptor *c1, C3 *c3, C6 *c6);

void str(ScattererConfiguration *sconf, ParticleDescriptor *c1, C6 *c6);



/*! \brief Compute radiation torques on particles on a k-vector oriented system.

 *

src/libnptm/Commons.cpp

+90 −56
Original line number Diff line number Diff line
@@ -36,61 +36,61 @@ 
#include <mpi.h>

#endif



C3::C3() {

  tsas = new dcomplex*[2];

  tsas[0] = new dcomplex[2];

  tsas[1] = new dcomplex[2];

  tfsas = 0.0 + 0.0 * I;

  gcs = 0.0;

  scs = 0.0;

  ecs = 0.0;

  acs = 0.0;

}



C3::C3(const C3& rhs) {

  tsas = new dcomplex*[2];

  for (int ti=0; ti<2; ti++) {

    tsas[ti] = new dcomplex[2];

    for (int tj=0; tj<2; tj++) tsas[ti][tj] = rhs.tsas[ti][tj];

  }

  tfsas = rhs.tfsas;

  gcs = rhs.gcs;

  scs = rhs.scs;

  ecs = rhs.ecs;

  acs = rhs.acs;

}



C3::~C3() {

  delete[] tsas[1];

  delete[] tsas[0];

  delete[] tsas;

}



#ifdef MPI_VERSION

C3::C3(const mixMPI *mpidata) {

  tsas = new dcomplex*[2];

  for (int ti=0; ti<2; ti++) {

    tsas[ti] = new dcomplex[2];

    MPI_Bcast(tsas[ti], 2, 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);

}



void C3::mpibcast(const mixMPI *mpidata) {

  for (int ti=0; ti<2; ti++) {

    MPI_Bcast(tsas[ti], 2, 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);

}

#endif

// C3::C3() {

//   tsas = new dcomplex*[2];

//   tsas[0] = new dcomplex[2];

//   tsas[1] = new dcomplex[2];

//   tfsas = 0.0 + 0.0 * I;

//   gcs = 0.0;

//   scs = 0.0;

//   ecs = 0.0;

//   acs = 0.0;

// }



// C3::C3(const C3& rhs) {

//   tsas = new dcomplex*[2];

//   for (int ti=0; ti<2; ti++) {

//     tsas[ti] = new dcomplex[2];

//     for (int tj=0; tj<2; tj++) tsas[ti][tj] = rhs.tsas[ti][tj];

//   }

//   tfsas = rhs.tfsas;

//   gcs = rhs.gcs;

//   scs = rhs.scs;

//   ecs = rhs.ecs;

//   acs = rhs.acs;

// }



// C3::~C3() {

//   delete[] tsas[1];

//   delete[] tsas[0];

//   delete[] tsas;

// }



// #ifdef MPI_VERSION

// C3::C3(const mixMPI *mpidata) {

//   tsas = new dcomplex*[2];

//   for (int ti=0; ti<2; ti++) {

//     tsas[ti] = new dcomplex[2];

//     MPI_Bcast(tsas[ti], 2, 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);

// }



// void C3::mpibcast(const mixMPI *mpidata) {

//   for (int ti=0; ti<2; ti++) {

//     MPI_Bcast(tsas[ti], 2, 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);

// }

// #endif



C6::C6(int _lmtpo) {

  lmtpo = _lmtpo;
@@ -234,7 +234,6 @@ mixMPI::~mixMPI() { 


ClusterIterationData::ClusterIterationData(GeometryConfiguration *gconf, ScattererConfiguration *sconf, const mixMPI *mpidata, const int device_count) {

  c1 = new ParticleDescriptorCluster(gconf, sconf);

  c3 = new C3();

  c6 = new C6(c1->lmtpo);

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

  const np_int ndit = 2 * ndi;
@@ -338,7 +337,6 @@ ClusterIterationData::ClusterIterationData(GeometryConfiguration *gconf, Scatter 


ClusterIterationData::ClusterIterationData(const ClusterIterationData& rhs) {

  c1 = new ParticleDescriptorCluster(reinterpret_cast<ParticleDescriptorCluster &>(*(rhs.c1)));

  c3 = new C3(*(rhs.c3));

  c6 = new C6(*(rhs.c6));

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

  const np_int ndit = 2 * ndi;
@@ -486,7 +484,6 @@ ClusterIterationData::ClusterIterationData(const ClusterIterationData& rhs) { 
#ifdef MPI_VERSION

ClusterIterationData::ClusterIterationData(const mixMPI *mpidata, const int device_count) {

  c1 = new ParticleDescriptorCluster(mpidata);

  c3 = new C3(mpidata);

  c6 = new C6(mpidata);

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

  const np_int ndit = 2 * ndi;
@@ -622,7 +619,6 @@ ClusterIterationData::ClusterIterationData(const mixMPI *mpidata, const int devi 


void ClusterIterationData::mpibcast(const mixMPI *mpidata) {

  c1->mpibcast(mpidata);

  c3->mpibcast(mpidata);

  c6->mpibcast(mpidata);

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

  const np_int ndit = 2 * ndi;
@@ -707,7 +703,6 @@ ClusterIterationData::~ClusterIterationData() { 
  }

  delete[] zpv;

  delete c1;

  delete c3;

  delete c6;

  delete c9;

  delete[] gaps;
@@ -853,7 +848,14 @@ ParticleDescriptor::ParticleDescriptor(GeometryConfiguration *gconf, ScattererCo 
  sqabs = NULL;

  gcsv = NULL;

  // >>> NEEDED BY CLUSTER <<<

  vec_tsas = NULL;

  vintt = NULL;

  tfsas = 0.0 + I * 0.0;

  tsas = NULL;

  gcs = 0.0;

  scs = 0.0;

  ecs = 0.0;

  acs = 0.0;

  // >>> NEEDED BY CLUSTER AND INCLU <<<

  vec_am0m = NULL;

  vec_fsac = NULL;
@@ -1158,6 +1160,12 @@ void ParticleDescriptor::mpibcast(const mixMPI *mpidata) { 
  // >>> NEEDED BY CLUSTER <<< //

  if (_class_type == CLUSTER_TYPE) {

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

  }

  // >>> NEEDED BY CLUSTER AND INCLU <<< //

  if (_class_type == CLUSTER_TYPE || _class_type == INCLUSION_TYPE) {
@@ -1280,6 +1288,8 @@ ParticleDescriptor::~ParticleDescriptor() { 
  // Cluster class members, destroyed only if sub-class is CLUSTER

  if (_class_type == CLUSTER_TYPE) {

    delete[] vintt;

    delete[] vec_tsas;

    delete[] tsas;

  }

}
@@ -1311,6 +1321,10 @@ ParticleDescriptorCluster::ParticleDescriptorCluster(GeometryConfiguration *gcon 


  // Needed by CLUSTER

  vintt = new dcomplex[16]();

  vec_tsas = new dcomplex[4]();

  tsas = new dcomplex*[2];

  tsas[0] = vec_tsas;

  tsas[1] = vec_tsas + 2;



  // Needed by CLUSTER and INCLU

  _le = gconf->le;
@@ -1395,6 +1409,16 @@ ParticleDescriptorCluster::ParticleDescriptorCluster(const ParticleDescriptorClu 
  // >>> NEEDED BY CLUSTER <<<

  vintt = new dcomplex[16];

  for (int ti = 0; ti < 16; ti++) vintt[ti] = rhs.vintt[ti];

  vec_tsas = new dcomplex[4];

  for (int si = 0; si < 4; si++) vec_tsas[si] = rhs.vec_tsas[si];

  tsas = new dcomplex*[2];

  tsas[0] = vec_tsas;

  tsas[1] = vec_tsas + 2;

  tfsas = rhs.tfsas;

  gcs = rhs.gcs;

  scs = rhs.scs;

  ecs = rhs.ecs;

  acs = rhs.acs;

  // >>> NEEDED BY CLUSTER AND INCLU <<<

  vec_am0m = new dcomplex[_nlemt * _nlemt];

  np_int nlemts = _nlemt * _nlemt;
@@ -1490,6 +1514,16 @@ ParticleDescriptorCluster::ParticleDescriptorCluster(const mixMPI *mpidata) : Pa 
  // >>> NEEDED BY CLUSTER <<< //

  vintt = new dcomplex[16];

  MPI_Bcast(vintt, 16, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  vec_tsas = new dcomplex[4];

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

  tsas = new dcomplex*[2];

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

  // >>> NEEDED BY CLUSTER AND INCLU <<<

  vec_am0m = new dcomplex[_nlemt * _nlemt];

  int nlemts = _nlemt * _nlemt;

src/libnptm/clu_subs.cpp

+14 −14
Original line number Diff line number Diff line
@@ -2179,7 +2179,7 @@ void rftr( 
  fz = u[2] * extins - gapv[2];

}



void scr0(double vk, double exri, ParticleDescriptor *c1, C3 *c3) {

void scr0(double vk, double exri, ParticleDescriptor *c1) {

  const dcomplex cc0 = 0.0 + 0.0 * I;

  double exdc = exri * exri;

  double ccs = 4.0 * acos(0.0) / (vk * vk);
@@ -2256,10 +2256,10 @@ void scr0(double vk, double exri, ParticleDescriptor *c1, C3 *c3) { 
    acs += c1->sabs[iogi - 1];

    tfsas += c1->fsas[iogi - 1];

  }

  c3->scs = scs;

  c3->ecs = ecs;

  c3->acs = acs;

  c3->tfsas = tfsas;

  c1->scs = scs;

  c1->ecs = ecs;

  c1->acs = acs;

  c1->tfsas = tfsas;

#ifdef USE_NVTX

  nvtxRangePop();

#endif
@@ -2267,7 +2267,7 @@ void scr0(double vk, double exri, ParticleDescriptor *c1, C3 *c3) { 


void scr2(

	  double vk, double vkarg, double exri, double *duk,

	  ParticleDescriptor *c1, C3 *c3

	  ParticleDescriptor *c1

) {

#ifdef USE_NVTX

  nvtxRangePush("scr2 starts");
@@ -2377,10 +2377,10 @@ void scr2( 
    tsas01 += (c1->sas[iogi - 1][0][1] * phas);

    tsas11 += (c1->sas[iogi - 1][1][1] * phas);

  } // i14 loop

  c3->tsas[0][0] = tsas00;

  c3->tsas[1][0] = tsas10;

  c3->tsas[0][1] = tsas01;

  c3->tsas[1][1] = tsas11;

  c1->tsas[0][0] = tsas00;

  c1->tsas[1][0] = tsas10;

  c1->tsas[0][1] = tsas01;

  c1->tsas[1][1] = tsas11;

#ifdef USE_NVTX

  nvtxRangePop();

  //#endif
@@ -2429,7 +2429,7 @@ void scr2( 
      for (int ipo2 = 1; ipo2 <= 2; ipo2++) {

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

	  int j = jpo2-1 + (ipo2-1)*2 + (jpo1-1)*4 + (ipo1-1)*8;

	  c1->vintt[j] = c3->tsas[jpo2 - 1][ipo2 - 1] * dconjg(c3->tsas[jpo1 - 1][ipo1 - 1]) * cfsq;

	  c1->vintt[j] = c1->tsas[jpo2 - 1][ipo2 - 1] * dconjg(c1->tsas[jpo1 - 1][ipo1 - 1]) * cfsq;

	} // jpo2 loop

      } // ipo2 loop

    } // jpo1 loop
@@ -2442,11 +2442,11 @@ void scr2( 
#endif

}



void str(ScattererConfiguration *sconf, ParticleDescriptor *c1, C3 *c3, C6 *c6) {

void str(ScattererConfiguration *sconf, ParticleDescriptor *c1, C6 *c6) {

  dcomplex *ylm;

  const double pi = acos(-1.0);

  int last_configuration;

  c3->gcs = 0.0;

  c1->gcs = 0.0;

  double gcss = 0.0;

  last_configuration = 0;

  for (int i18 = 1; i18 <= c1->nsph; i18++) {
@@ -2460,7 +2460,7 @@ void str(ScattererConfiguration *sconf, ParticleDescriptor *c1, C3 *c3, C6 *c6) 
	c1->rc[last_configuration - 1][j16 - 1] = sconf->get_rcf(last_configuration - 1, j16 - 1) * c1->ros[last_configuration - 1];

      } // j16 loop

    }

    c3->gcs += gcss;

    c1->gcs += gcss;

  } // i18 loop

  int ylm_size = (c1->litpos > c1->lmtpos) ? c1->litpos : c1->lmtpos;

  ylm = new dcomplex[ylm_size]();