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

Add ParticleDescriptorCluster derived class

parent 1047b2f4
Loading
Loading
Loading
Loading
+64 −12
Original line number Diff line number Diff line
@@ -684,8 +684,24 @@ protected:
  int _le;
  //! \brief Maximum field expansion order.
  int _lm;
  //! \brief Number of rows of the Transition Matrix.
  int _nlmmt;
  //! \brief NLIM = NSPH * LI * (LI + 2)
  int _nlim;
  //! \brief NLEM = NSPH * LE * (LE + 2)
  int _nlem;
  //! \brief NLEMT = 2 * NLEM
  int _nlemt;
  //! \brief NCOU = NSPH * NSPH - 1
  int _ncou;
  //! \brief LITPO = 2 * LI + 1
  int _litpo;
  //! \brief LITPOS = LITPO * LITPO
  int _litpos;
  //! \brief LMTPO = LI + LE + 1
  int _lmtpo;
  //! \brief LMTPOS = LMTPO * LMTPO
  int _lmtpos;
  //! \brief NV3J = (LM * (LM + 1) * (2 * LM + 7)) / 6
  int _nv3j;
  //! \brief Contiguous space for AM0M
  dcomplex *vec_am0m;
  //! \brief Contiguous space for FSAC
@@ -694,7 +710,9 @@ protected:
  dcomplex *vec_sac;
  //! \brief Contiguous space for FSACM
  dcomplex *vec_fsacm;
  // >>> END OF SECTION NEEDED BY CLUSTER AND INCLU <<<
  //! \brief Contiguous space for IND3J
  int *vec_ind3j;
  // >>> END OF SECTION NEEDED BY CLUSTER AND INCLU <<< //

public:
  // >>> COMMON TO ALL DESCRIPTOR TYPES <<< //
@@ -706,8 +724,6 @@ public:
  const int &num_configurations = _num_configurations;
  //! \brief Read-only view of total number of layers from all sphere types.
  const int &num_layers = _num_layers;
  //! \brief Vector of number of layers in type.
  int *num_layers_in_type;
  //! \brief Matrix of inverse radial M coefficients.
  dcomplex **rmi;
  //! \brief Matrix of inverse radial E coefficients.
@@ -753,41 +769,77 @@ public:
  double *sqabs;
  //! \brief Vector of geometric cross-sections for spheres.
  double *gcsv;
  // >>> END OF SECTION NEEDED BY SPHERE AND CLUSTER <<<
  // >>> END OF SECTION NEEDED BY SPHERE AND CLUSTER <<< //

  // >>> NEEDED BY CLUSTER <<<
  dcomplex *vintt;
  // >>> END OF SECTION NEEDED BY CLUSTER <<<
  // >>> END OF SECTION NEEDED BY CLUSTER <<< //
  
  // >>> NEEDED BY CLUSTER AND INCLU <<<
  //! \brief Read-only view of maximum external field expansion order.
  const int& le = _le;
  //! \brief Read-only view of maximum field expansion order.
  const int& lm = _lm;
  //! \brief Read-only view of number of rows of the Transition Matrix.
  const int& nlmmt = _nlmmt;

  //! \brief Read-only view of NLIM.
  const int& nlim = _nlim;
  //! \brief Read-only view of NLEM.
  const int& nlem = _nlem;
  //! \brief Read-only view of NLEMT.
  const int& nlemt = _nlemt;
  //! \brief Read-only view of NCOU.
  const int& ncou = _ncou;
  //! \brief Read-only view of LITPO.
  const int& litpo = _litpo;
  //! \brief Read-only view of LITPOS.
  const int& litpos = _litpos;
  //! \brief Read-only view of LMTPO.
  const int& lmtpo = _lmtpo;
  //! \brief Read-only view of LMTPOS.
  const int& lmtpos = _lmtpos;
  //! \brief Read-only view of NV3J.
  const int& nv3j = _nv3j;

  //! \brief TBD
  dcomplex *vh;
  //! \brief TBD
  dcomplex *vj0;
  //! \brief TBD
  dcomplex *vyhj;
  //! \brief TBD
  dcomplex *vyj0;
  //! \brief TBD
  dcomplex vj;
  //! \brief Transition matrix
  dcomplex **am0m;
  //! \brief Cluster forward scattering amplitude.
  dcomplex **fsac;
  //! \brief Cluster scattering amplitude.
  dcomplex **sac;
  //! \brief Mean cluster forward scattering amplitude.
  dcomplex **fsacm;
  //! \brief Mean intensity components vector.
  dcomplex *vintm;
  //! \brief Cluster polarized scattering cross-sections.
  dcomplex *scscp;
  //! \brief Cluster polarized extinction cross-sections.
  dcomplex *ecscp;
  //! \brief Mean cluster polarized scattering cross-sections.
  dcomplex *scscpm;
  //! \brief Mean cluster polarized extinction cross-sections.
  dcomplex *ecscpm;
  //! \brief TBD
  double *v3j0;
  //! \brief Cluster scattering cross-sections.
  double *scsc;
  //! \brief Cluster extinction cross-sections.
  double *ecsc;
  //! \brief Mean cluster scattering cross-sections.
  double *scscm;
  //! \brief Mean cluster extinction cross-sections.
  double *ecscm;
  int *ind3j;
  // >>> END OF SECTION NEEDED BY CLUSTER AND INCLU <<<
  //! \brief J-vector components index matrix.
  int **ind3j;
  // >>> END OF SECTION NEEDED BY CLUSTER AND INCLU <<< //
  
  /*! \brief ParticleDescriptor instance constructor.
   *
+136 −9
Original line number Diff line number Diff line
@@ -1521,14 +1521,21 @@ ClusterIterationData::~ClusterIterationData() {
ParticleDescriptor::ParticleDescriptor(GeometryConfiguration *gconf, ScattererConfiguration *sconf) {
  _nsph = gconf->number_of_spheres;
  _li = (nsph == 1) ? gconf->l_max : gconf->li;
  _le = (nsph == 1) ? 0 : gconf->le;
  _lm = (_le > _li) ? _le : _li;
  _nlmmt = (nsph == 1) ? 2 * li * (li + 2) : 2 * _le * (_le + 2);
  _le = 0;
  _lm = _li;
  _nlim = 0;
  _nlem = 0;
  _nlemt = 0;
  _ncou = 0;
  _litpo = 0;
  _litpos = 0;
  _lmtpo = 0;
  _lmtpos = 0;
  _nv3j = 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);

  num_layers_in_type = new int[num_configurations]();
  vec_rmi = new dcomplex[li * nsph]();
  rmi = new dcomplex*[li];
  vec_rei = new dcomplex[li * nsph]();
@@ -1537,6 +1544,7 @@ ParticleDescriptor::ParticleDescriptor(GeometryConfiguration *gconf, ScattererCo
    rmi[ri] = vec_rmi + (nsph * ri);
    rei[ri] = vec_rei + (nsph * ri);
  }
  int nllt = (_nlemt = 0) ? 2 * _nsph * _li * (_li + 2) : _nlemt;
  vec_w = new dcomplex[nllt * 4]();
  w = new dcomplex*[nllt];
  for (int wi = 0; wi < nllt; wi++) w[wi] = vec_w + (4 * wi);
@@ -1547,7 +1555,6 @@ ParticleDescriptor::ParticleDescriptor(GeometryConfiguration *gconf, ScattererCo
    rc[rci] = vec_rc + last_layer_index;
    last_layer_index += sconf->get_nshl(rci);
    if (rci == 0 && sconf->use_external_sphere) last_layer_index++;
    num_layers_in_type[rci] = last_layer_index - cur_layer_index;
    for (int rcj = cur_layer_index; rcj < last_layer_index; rcj++) {
      rc[rci][rcj] = sconf->get_rcf(rci, rcj);
    }
@@ -1581,13 +1588,14 @@ ParticleDescriptor::ParticleDescriptor(GeometryConfiguration *gconf, ScattererCo
  sqexs = NULL;
  sqabs = NULL;
  gcsv = NULL;
  // >>> NEEDED BY CLUSTER <<<
  vintt = NULL;
  // >>> NEEDED BY CLUSTER AND INCLU <<<
  vec_am0m = NULL;
  vec_fsac = NULL;
  vec_sac = NULL;
  vec_fsacm = NULL;
  // >>> NEEDED BY CLUSTER <<<
  vintt = NULL;
  // >>> NEEDED BY CLUSTER AND INCLU <<<
  vec_ind3j = NULL;
  vh = NULL;
  vj0 = NULL;
  vyhj = NULL;
@@ -1626,11 +1634,130 @@ ParticleDescriptor::~ParticleDescriptor() {
  delete[] vec_rei;
  delete[] rmi;
  delete[] vec_rmi;
  delete[] num_layers_in_type;
}

// >>> End of ParticleDescriptor class implementation. <<< //

// >>> ParticleDescriptorCluster class implementation. <<< //
ParticleDescriptorCluster::ParticleDescriptorCluster(GeometryConfiguration *gconf, ScattererConfiguration *sconf) : ParticleDescriptor(gconf, sconf) {
  // Needed by SPHERE and CLUSTER
  vec_sas = new dcomplex[4 * nsph]();
  vec_vints = new dcomplex[16 * nsph]();

  fsas = new dcomplex[nsph];
  sscs = new double[nsph]();
  sexs = new double[nsph]();
  sabs = new double[nsph]();
  sqscs = new double[nsph]();
  sqexs = new double[nsph]();
  sqabs = new double[nsph]();
  gcsv = new double[nsph]();
  sas = new dcomplex**[nsph];
  vints = new dcomplex*[nsph];
  for (int vi = 0; vi < nsph; vi++) {
    vints[vi] = vec_vints + (16 * nsph);
    sas[vi] = new dcomplex*[2];
    sas[vi][0] = vec_sas + (4 * vi);
    sas[vi][1] = vec_sas + (4 * vi) + 2;
  }

  // Needed by CLUSTER
  vintt = new dcomplex[16];

  // Needed by CLUSTER and INCLU
  _le = gconf->le;
  _lm = (_li > _le) ? _li : _le;
  _nlim = _nsph * _li * (_li + 2);
  _nlem = _nsph * _le * (_le + 2);
  _nlemt = 2 * _nlem;
  _ncou = _nsph * _nsph - 1;
  _litpo = _li + _li + 1;
  _litpos = _litpo * _litpo;
  _lmtpo = _li + _le + 1;
  _lmtpos = _lmtpo * _lmtpo;
  _nv3j = (_lm * (_lm + 1) * (2 * _lm + 7)) / 6;
  vec_am0m = new dcomplex[_nlemt * _nlemt]();
  vec_fsac = new dcomplex[4]();
  vec_sac = new dcomplex[4]();
  vec_fsacm = new dcomplex[4]();
  vec_ind3j = new int[(_lm + 1) * _lm]();
  
  vh = new dcomplex[_ncou * _litpo]();
  vj0 = new dcomplex[_nsph * _lmtpo]();
  vyhj = new dcomplex[_ncou * _litpos]();
  vyj0 = new dcomplex[_nsph * _lmtpos]();
  am0m = new dcomplex*[_nlemt];
  for (int ai = 0; ai < _nlemt; ai++) am0m[ai] = vec_am0m + (ai * _nlemt);
  fsac = new dcomplex*[2];
  fsacm = new dcomplex*[2];
  for (int fi = 0; fi < 2; fi++) {
    fsac[fi] = vec_fsac + (fi * 2);
    fsacm[fi] = vec_fsacm + (fi * 2);
  }
  vintm = new dcomplex[16]();
  scscp = new dcomplex[2]();
  ecscp = new dcomplex[2]();
  scscpm = new dcomplex[2]();
  ecscpm = new dcomplex[2]();
  v3j0 = new double[_nv3j]();
  scsc = new double[2]();
  ecsc = new double[2]();
  scscm = new double[2]();
  ecscm = new double[2]();
  ind3j = new int*[_lm + 1];
  for (int ii = 0; ii <= _lm; ii++) ind3j[ii] = vec_ind3j + (_lm * ii);
}

ParticleDescriptorCluster::~ParticleDescriptorCluster() {
  // Needed by SPHERE and CLUSTER
  for (int vi = 0; vi < nsph; vi++) {
    delete[] sas[vi];
  }
  delete[] vints;
  delete[] sas;
  delete[] gcsv;
  delete[] sqabs;
  delete[] sqexs;
  delete[] sqscs;
  delete[] sabs;
  delete[] sexs;
  delete[] sscs;
  delete[] fsas;
  delete[] vec_vints;
  delete[] vec_sas;
  
  // Needed by CLUSTER
  delete[] vintt;

  // Needed by CLUSTER and INCLU
  delete[] vec_am0m;
  delete[] vec_fsac;
  delete[] vec_sac;
  delete[] vec_fsacm;
  delete[] vec_ind3j;
  
  delete[] vh;
  delete[] vj0;
  delete[] vyhj;
  delete[] vyj0;
  delete[] am0m;
  delete[] fsac;
  delete[] fsacm;
  delete[] vintm;
  delete[] scscp;
  delete[] ecscp;
  delete[] scscpm;
  delete[] ecscpm;
  delete[] v3j0;
  delete[] scsc;
  delete[] ecsc;
  delete[] scscm;
  delete[] ecscm;
  delete[] ind3j;
}

// >>> End of ParticleDescriptorCluster class implementation. <<< //

// >>> ParticleDescriptorSphere class implementation. <<< //
ParticleDescriptorSphere::ParticleDescriptorSphere(GeometryConfiguration *gconf, ScattererConfiguration *sconf) : ParticleDescriptor(gconf, sconf) {
  vec_sas = new dcomplex[4 * nsph]();