From 821dd790f1c4c3b240684b8f9f6906ee6e9b2bda Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Fri, 22 Nov 2024 11:54:40 +0100
Subject: [PATCH] Move C2 data structures under ParticleDescriptor
---
src/cluster/cluster.cpp | 12 ++--
src/include/Commons.h | 84 +++++++----------------
src/include/inclu_subs.h | 5 +-
src/include/sph_subs.h | 10 ++-
src/libnptm/Commons.cpp | 137 ++++++++++++++++---------------------
src/libnptm/inclu_subs.cpp | 22 +++---
src/libnptm/sph_subs.cpp | 54 +++++++--------
src/sphere/sphere.cpp | 19 +++--
8 files changed, 142 insertions(+), 201 deletions(-)
diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index 9bd5b207..8d139b34 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -754,17 +754,17 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
int ici = (nsh + 1) / 2;
if (idfc == 0) {
for (int ic = 0; ic < ici; ic++)
- cid->c2->dc0[ic] = sconf->get_dielectric_constant(ic, i132 - 1, jxi488 - 1);
+ cid->c1->dc0[ic] = sconf->get_dielectric_constant(ic, i132 - 1, jxi488 - 1);
} else {
if (jxi488 == 1) {
for (int ic = 0; ic < ici; ic++)
- cid->c2->dc0[ic] = sconf->get_dielectric_constant(ic, i132 - 1, 0);
+ cid->c1->dc0[ic] = sconf->get_dielectric_constant(ic, i132 - 1, 0);
}
}
- if (nsh % 2 == 0) cid->c2->dc0[ici] = exdc;
+ if (nsh % 2 == 0) cid->c1->dc0[ici] = exdc;
dme(
cid->c1->li, i132, npnt, npntts, vkarg, exdc, exri,
- cid->c1, cid->c2, jer, lcalc, cid->arg, last_configuration
+ cid->c1, jer, lcalc, cid->arg, last_configuration
);
if (jer != 0) {
sprintf(virtual_line, " STOP IN DME\n");
@@ -915,10 +915,10 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
sprintf(virtual_line, " SPHERE %2d\n", i170);
output->append_line(virtual_line);
if (cid->c1->nshl[last_configuration - 1] != 1) {
- sprintf(virtual_line, " SIZE=%15.7lE\n", cid->c2->vsz[i]);
+ sprintf(virtual_line, " SIZE=%15.7lE\n", cid->c1->vsz[i]);
output->append_line(virtual_line);
} else { // label 162
- sprintf(virtual_line, " SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n", cid->c2->vsz[i], real(cid->c2->vkt[i]), imag(cid->c2->vkt[i]));
+ sprintf(virtual_line, " SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n", cid->c1->vsz[i], real(cid->c1->vkt[i]), imag(cid->c1->vkt[i]));
output->append_line(virtual_line);
}
// label 164
diff --git a/src/include/Commons.h b/src/include/Commons.h
index 282cf57d..a7eec18d 100644
--- a/src/include/Commons.h
+++ b/src/include/Commons.h
@@ -40,62 +40,6 @@
class ParticleDescriptor;
class mixMPI;
-/*! \brief Representation of the FORTRAN C2 blocks.
- *
- */
-class C2 {
-protected:
- //! \brief Number of spheres.
- int nsph;
- //! \brief Number of required orders.
- int nhspo;
- //! \brief QUESTION: what is nl?
- int nl;
-
-public:
- //! \brief QUESTION: definition?
- dcomplex *ris;
- //! \brief QUESTION: definition?
- dcomplex *dlri;
- //! \brief Vector of dielectric constants.
- dcomplex *dc0;
- //! \brief QUESTION: definition?
- dcomplex *vkt;
- //! Vector of sphere sizes in units of 2*PI/LAMBDA
- double *vsz;
-
- /*! \brief C2 instance constructor.
- *
- * \param gconf: `GeometryConfiguration*` Pointer to a GeometryConfiguration instance.
- * \param sconf: `ScattererConfiguration*` Pointer to a ScattererConfiguration instance.
- */
- C2(GeometryConfiguration *gconf, ScattererConfiguration *sconf);
-
- /*! \brief C2 instance constructor copying its contents from preexisting instance.
- *
- * \param rhs: `C2` object to copy contents from
- */
- C2(const C2& rhs);
-
- //! \brief C2 instance destroyer.
- ~C2();
-
-#ifdef MPI_VERSION
- /*! \brief C2 instance constructor copying all contents off MPI broadcast from MPI process 0
- *
- * \param mpidata: `mixMPI *` pointer to MPI data structure.
- */
- C2(const mixMPI *mpidata);
-
- /*! \brief send C2 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 {
@@ -276,8 +220,6 @@ class ClusterIterationData {
public:
//! \brief Pointer to a C1 structure.
ParticleDescriptor *c1;
- //! \brief Pointer to a C2 structure.
- C2 *c2;
//! \brief Pointer to a C3 structure.
C3 *c3;
//! \brief Pointer to a C6 structure.
@@ -387,6 +329,14 @@ protected:
int _num_configurations;
//! \brief Total number of layers from all sphere types.
int _num_layers;
+ //! \brief Space for different sphere types.
+ int _nl;
+ //! \brief NHSPO = 2 * MAX(NPNT,NPNTTS) - 1
+ int _nhspo;
+ //! \brief Number of points for numerical integration in layered spheres.
+ int _npnt;
+ //! \brief Number of points for numerical integration in transition layer.
+ int _npntts;
//! \brief Contiguous space for RMI.
dcomplex *vec_rmi;
//! \brief Contiguous space for REI.
@@ -462,6 +412,14 @@ 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 Read-only view on the space for different sphere configurations.
+ const int &nl = _nl;
+ //! \brief Read-only view of NHSPO.
+ const int &nhspo = _nhspo;
+ //! \brief Read-only view on number of points for numerical integration in layered spheres.
+ const int &npnt = _npnt;
+ //! \brief Read-only view on number of points for numerical integration in transition layer.
+ const int &npntts = _npntts;
//! \brief Matrix of inverse radial M coefficients.
dcomplex **rmi;
//! \brief Matrix of inverse radial E coefficients.
@@ -484,6 +442,16 @@ public:
int *iog;
//! \brief Vector of number of layers in sphere type.
int *nshl;
+ //! \brief TBD
+ dcomplex *ris;
+ //! \brief TBD
+ dcomplex *dlri;
+ //! \brief TBD
+ dcomplex *dc0;
+ //! \brief TBD
+ dcomplex *vkt;
+ //! \brief Vector of sizes in units of 2*PI/LAMBDA
+ double *vsz;
// >>> END OF SECTION COMMON TO ALL DESCRIPTOR TYPES <<< //
// >>> NEEDED BY SPHERE AND CLUSTER <<< //
diff --git a/src/include/inclu_subs.h b/src/include/inclu_subs.h
index 23ee14db..ccc034fb 100644
--- a/src/include/inclu_subs.h
+++ b/src/include/inclu_subs.h
@@ -71,13 +71,10 @@ void incms(dcomplex **am, dcomplex **at, double enti, ParticleDescriptor *c1, C6
* \param lcalc: `int &` Maximum order achieved in calculation.
* \param arg: `dcomplex &` TBD.
* \param c1: `ParticleDescriptor *` Pointer to a ParticleDescriptor instance.
- * \param c2: `C2 *` Pointer to a C2 instance.
*/
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
- );
+ dcomplex entn, int &jer, int &lcalc, dcomplex &arg, ParticleDescriptor *c1 );
/*! \brief C++ porting of INSTR.
*
diff --git a/src/include/sph_subs.h b/src/include/sph_subs.h
index 0ae5b280..d4d9bc7b 100644
--- a/src/include/sph_subs.h
+++ b/src/include/sph_subs.h
@@ -89,9 +89,8 @@ double cg1(int lmpml, int mu, int l, int m);
* \param ic: `int`
* \param vk: `double`
* \param c1: `ParticleDescriptor *` Pointer to `ParticleDescriptor` data structure.
- * \param c2: `C2 *` Pointer to `C2` data structure.
*/
-void diel(int npntmo, int ns, int i, int ic, double vk, ParticleDescriptor *c1, C2 *c2);
+void diel(int npntmo, int ns, int i, int ic, double vk, ParticleDescriptor *c1);
/*! \brief Compute Mie scattering coefficients.
*
@@ -108,7 +107,6 @@ void diel(int npntmo, int ns, int i, int ic, double vk, ParticleDescriptor *c1,
* \param exdc: `double` External medium dielectric constant.
* \param exri: `double` External medium refractive index.
* \param c1: `ParticleDescriptor *` Pointer to a `ParticleDescriptor` data structure.
- * \param c2: `C2 *` Pointer to a `C2` data structure.
* \param jer: `int &` Reference to integer error code variable.
* \param lcalc: `int &` Reference to integer variable recording the maximum expansion order accounted for.
* \param arg: `complex double &`
@@ -116,7 +114,7 @@ void diel(int npntmo, int ns, int i, int ic, double vk, ParticleDescriptor *c1,
*/
void dme(
int li, int i, int npnt, int npntts, double vk, double exdc, double exri,
- ParticleDescriptor *c1, C2 *c2, int &jer, int &lcalc, dcomplex &arg, int last_conf=0
+ ParticleDescriptor *c1, int &jer, int &lcalc, dcomplex &arg, int last_conf=0
);
/*! \brief Bessel function calculation control parameters.
@@ -262,11 +260,11 @@ void rkc(
* \param y2: `complex double &`
* \param dy1: `complex double &`
* \param dy2: `complex double &`
- * \param c2: `C2 *` Pointer to a `C2` data structure.
+ * \param c1: `ParticleDescriptor *` Pointer to a ParticleDescriptor instance.
*/
void rkt(
int npntmo, double step, double &x, int lpo, dcomplex &y1,
- dcomplex &y2, dcomplex &dy1, dcomplex &dy2, C2 *c2
+ dcomplex &y2, dcomplex &dy1, dcomplex &dy2, ParticleDescriptor *c1
);
/*! \brief Spherical Bessel functions.
diff --git a/src/libnptm/Commons.cpp b/src/libnptm/Commons.cpp
index df6d43e5..75f477a2 100644
--- a/src/libnptm/Commons.cpp
+++ b/src/libnptm/Commons.cpp
@@ -36,78 +36,6 @@
#include <mpi.h>
#endif
-C2::C2(GeometryConfiguration *gconf, ScattererConfiguration *sconf) {
- nsph = gconf->number_of_spheres;
- int npnt = gconf->npnt;
- int npntts = gconf->npntts;
- int max_n = (npnt > npntts) ? npnt : npntts;
- nhspo = 2 * max_n - 1;
- nl = sconf->configurations;
- if (nsph == 1 && nl == 1) nl = 5;
- ris = new dcomplex[nhspo]();
- dlri = new dcomplex[nhspo]();
- vkt = new dcomplex[nsph]();
- dc0 = new dcomplex[nl]();
- vsz = new double[nsph]();
-}
-
-C2::C2(const C2& rhs) {
- nsph = rhs.nsph;
- nhspo = rhs.nhspo;
- nl = rhs.nl;
- ris = new dcomplex[nhspo]();
- dlri = new dcomplex[nhspo]();
- for (int ind=0; ind<nhspo; ind++) {
- ris[ind] = rhs.ris[ind];
- dlri[ind] = rhs.dlri[ind];
- }
- vkt = new dcomplex[nsph]();
- vsz = new double[nsph]();
- for (int ind=0; ind<nsph; ind++) {
- vkt[ind] = rhs.vkt[ind];
- vsz[ind] = rhs.vsz[ind];
- }
- dc0 = new dcomplex[nl]();
- for (int ind=0; ind<nl; ind++) dc0[ind] = rhs.dc0[ind];
-}
-
-C2::~C2() {
- delete[] ris;
- delete[] dlri;
- delete[] vkt;
- delete[] dc0;
- delete[] vsz;
-}
-
-#ifdef MPI_VERSION
-C2::C2(const mixMPI *mpidata) {
- MPI_Bcast(&nsph, 1, MPI_INT, 0, MPI_COMM_WORLD);
- MPI_Bcast(&nhspo, 1, MPI_INT, 0, MPI_COMM_WORLD);
- MPI_Bcast(&nl, 1, MPI_INT, 0, MPI_COMM_WORLD);
- ris = new dcomplex[nhspo]();
- dlri = new dcomplex[nhspo]();
- MPI_Bcast(ris, nhspo, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
- MPI_Bcast(dlri, nhspo, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
- vkt = new dcomplex[nsph]();
- vsz = new double[nsph]();
- MPI_Bcast(vkt, nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
- MPI_Bcast(vsz, nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- dc0 = new dcomplex[nl]();
- MPI_Bcast(dc0, nl, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
-}
-
-void C2::mpibcast(const mixMPI *mpidata) {
- MPI_Bcast(&nsph, 1, MPI_INT, 0, MPI_COMM_WORLD);
- MPI_Bcast(&nhspo, 1, MPI_INT, 0, MPI_COMM_WORLD);
- MPI_Bcast(&nl, 1, MPI_INT, 0, MPI_COMM_WORLD);
- MPI_Bcast(ris, nhspo, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
- MPI_Bcast(dlri, nhspo, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
- MPI_Bcast(vkt, nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
- MPI_Bcast(vsz, nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- MPI_Bcast(dc0, nl, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
-}
-#endif
-
C3::C3() {
tsas = new dcomplex*[2];
tsas[0] = new dcomplex[2];
@@ -306,7 +234,6 @@ mixMPI::~mixMPI() {
ClusterIterationData::ClusterIterationData(GeometryConfiguration *gconf, ScattererConfiguration *sconf, const mixMPI *mpidata, const int device_count) {
c1 = new ParticleDescriptorCluster(gconf, sconf);
- c2 = new C2(gconf, sconf);
c3 = new C3();
c6 = new C6(c1->lmtpo);
const int ndi = c1->nsph * c1->nlim;
@@ -411,7 +338,6 @@ ClusterIterationData::ClusterIterationData(GeometryConfiguration *gconf, Scatter
ClusterIterationData::ClusterIterationData(const ClusterIterationData& rhs) {
c1 = new ParticleDescriptorCluster(reinterpret_cast<ParticleDescriptorCluster &>(*(rhs.c1)));
- c2 = new C2(*(rhs.c2));
c3 = new C3(*(rhs.c3));
c6 = new C6(*(rhs.c6));
const int ndi = c1->nsph * c1->nlim;
@@ -560,7 +486,6 @@ ClusterIterationData::ClusterIterationData(const ClusterIterationData& rhs) {
#ifdef MPI_VERSION
ClusterIterationData::ClusterIterationData(const mixMPI *mpidata, const int device_count) {
c1 = new ParticleDescriptorCluster(mpidata);
- c2 = new C2(mpidata);
c3 = new C3(mpidata);
c6 = new C6(mpidata);
const int ndi = c1->nsph * c1->nlim;
@@ -697,7 +622,6 @@ ClusterIterationData::ClusterIterationData(const mixMPI *mpidata, const int devi
void ClusterIterationData::mpibcast(const mixMPI *mpidata) {
c1->mpibcast(mpidata);
- c2->mpibcast(mpidata);
c3->mpibcast(mpidata);
c6->mpibcast(mpidata);
const int ndi = c1->nsph * c1->nlim;
@@ -783,7 +707,6 @@ ClusterIterationData::~ClusterIterationData() {
}
delete[] zpv;
delete c1;
- delete c2;
delete c3;
delete c6;
delete c9;
@@ -906,6 +829,18 @@ ParticleDescriptor::ParticleDescriptor(GeometryConfiguration *gconf, ScattererCo
ros[ci] = sconf->get_radius(ci);
nshl[ci] = sconf->get_nshl(ci);
}
+ _npnt = gconf->npnt;
+ _npntts = gconf->npntts;
+ int max_n = (npnt > npntts) ? npnt : npntts;
+ _nhspo = 2 * max_n - 1;
+ _nl = sconf->configurations;
+ if (_nsph == 1 && _nl == 1) _nl = 5;
+ ris = new dcomplex[_nhspo]();
+ dlri = new dcomplex[_nhspo]();
+ vkt = new dcomplex[_nsph]();
+ dc0 = new dcomplex[_nl]();
+ vsz = new double[_nsph]();
+
// >>> NEEDED BY SPHERE AND CLUSTER <<<
sas = NULL;
vints = NULL;
@@ -1021,6 +956,26 @@ ParticleDescriptor::ParticleDescriptor(const ParticleDescriptor &rhs) {
ros[ci] = rhs.ros[ci];
nshl[ci] = rhs.nshl[ci];
}
+ _npnt = rhs._npnt;
+ _npntts = rhs._npntts;
+ _nhspo = rhs._nhspo;
+ _nl = rhs._nl;
+ ris = new dcomplex[_nhspo]();
+ dlri = new dcomplex[_nhspo]();
+ for (int ri = 0; ri < _nhspo; ri++) {
+ ris[ri] = rhs.ris[ri];
+ dlri[ri] = rhs.dlri[ri];
+ }
+ vkt = new dcomplex[_nsph]();
+ vsz = new double[_nsph]();
+ for (int vi = 0; vi < _nsph; vi++) {
+ vkt[vi] = rhs.vkt[vi];
+ vsz[vi] = rhs.vsz[vi];
+ }
+ dc0 = new dcomplex[_nl]();
+ for (int di = 0; di < _nl; di++) {
+ dc0[di] = rhs.dc0[di];
+ }
// >>> NEEDED BY SPHERE AND CLUSTER <<<
sas = NULL;
vints = NULL;
@@ -1131,6 +1086,20 @@ ParticleDescriptor::ParticleDescriptor(const mixMPI *mpidata) {
MPI_Bcast(iog, _nsph, MPI_INT, 0, MPI_COMM_WORLD);
ros = new double[_num_configurations];
MPI_Bcast(ros, _num_configurations, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&_npnt, 1, MPI_INT, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&_npntts, 1, MPI_INT, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&_nhspo, 1, MPI_INT, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&_nl, 1, MPI_INT, 0, MPI_COMM_WORLD);
+ ris = new dcomplex[_nhspo];
+ MPI_Bcast(ris, _nhspo, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+ dlri = new dcomplex[_nhspo];
+ MPI_Bcast(dlri, _nhspo, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+ vkt = new dcomplex[_nsph];
+ MPI_Bcast(vkt, _nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+ vsz = new double[_nsph];
+ MPI_Bcast(vsz, _nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ dc0 = new dcomplex[_nl];
+ MPI_Bcast(dc0, _nl, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
}
void ParticleDescriptor::mpibcast(const mixMPI *mpidata) {
@@ -1163,6 +1132,15 @@ void ParticleDescriptor::mpibcast(const mixMPI *mpidata) {
MPI_Bcast(rzz, _nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Bcast(iog, _nsph, MPI_INT, 0, MPI_COMM_WORLD);
MPI_Bcast(ros, _num_configurations, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&_npnt, 1, MPI_INT, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&_npntts, 1, MPI_INT, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&_nhspo, 1, MPI_INT, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&_nl, 1, MPI_INT, 0, MPI_COMM_WORLD);
+ MPI_Bcast(ris, _nhspo, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+ MPI_Bcast(dlri, _nhspo, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+ MPI_Bcast(vkt, _nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+ MPI_Bcast(vsz, _nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+ MPI_Bcast(dc0, _nl, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
// >>> NEEDED BY SPHERE AND CLUSTER <<< //
if (_class_type == SPHERE_TYPE || _class_type == CLUSTER_TYPE) {
MPI_Bcast(vec_sas, 4 * _nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
@@ -1237,6 +1215,11 @@ ParticleDescriptor::~ParticleDescriptor() {
delete[] vec_rei;
delete[] rmi;
delete[] vec_rmi;
+ delete[] ris;
+ delete[] dlri;
+ delete[] vkt;
+ delete[] dc0;
+ delete[] vsz;
// Inclusion class members, destroyed only if sub-class is INCLUSION
if (_class_type == INCLUSION_TYPE) {
delete[] rm0;
diff --git a/src/libnptm/inclu_subs.cpp b/src/libnptm/inclu_subs.cpp
index 2ae6b22d..4a7f9cba 100644
--- a/src/libnptm/inclu_subs.cpp
+++ b/src/libnptm/inclu_subs.cpp
@@ -302,9 +302,7 @@ void incms(dcomplex **am, dcomplex **at, double enti, ParticleDescriptor *c1, C6
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
-) {
+ dcomplex entn, int &jer, int &lcalc, dcomplex &arg, ParticleDescriptor *c1) {
const dcomplex uim = 0.0 + I * 1.0;
const int nstp = npnt - 1;
const int nstpts = npntts - 1;
@@ -319,11 +317,11 @@ void indme(
double *rfn = new double[lipt]();
jer = 0;
double sz = vk * c1->ros[i - 1];
- c2->vsz[i - 1] = sz;
+ c1->vsz[i - 1] = sz;
double vkr1 = vk * c1->rc[i - 1][0];
int nsh = c1->nshl[i - 1];
- c2->vkt[i - 1] = csqrt(c2->dc0[0]);
- arg = vkr1 * c2->vkt[i - 1];
+ c1->vkt[i - 1] = csqrt(c1->dc0[0]);
+ arg = vkr1 * c1->vkt[i - 1];
dcomplex arin = arg;
if (imag(arg) != 0.0) {
cbf(lipo, arg, lcalc, cfj);
@@ -418,7 +416,7 @@ void indme(
dcomplex *dref = new dcomplex[c1->li]();
int ic = 0;
if (nsh < 2) { // nsh == 1
- dcomplex cri = c2->dc0[0] / ent;
+ dcomplex cri = c1->dc0[0] / ent;
for (int l160 = 1; l160 <= c1->li; l160++) {
int lpo = l160 + 1;
int ltpo = lpo + l160;
@@ -440,7 +438,7 @@ void indme(
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 dy1 = (l180 * fbi[l180 - 1] - lpo * fbi[lpt - 1]) * c1->vkt[i - 1] / dltpo;
dcomplex y2 = y1;
dcomplex dy2 = dy1;
ic = 0;
@@ -450,12 +448,12 @@ void indme(
if (ns % 2 != 0) {
// ic is incremented before being read in this loop.
int step = vk * (c1->rc[i - 1][ns - 1] - c1->rc[i - 1][nsmo - 1]) / nstp;
- arg = c2->dc0[++ic];
+ arg = c1->dc0[++ic];
rkc(nstp, step, arg, vkr, lpo, y1, y2, dy1, dy2);
} else { // label 170
- diel(nstpts, nsmo, i, ic, vk, c1, c2);
+ diel(nstpts, nsmo, i, ic, vk, c1);
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, c2);
+ rkt(nstpts, stepts, vkr, lpo, y1, y2, dy1, dy2, c1);
}
} // ns loop 176
rmf[l180 - 1] = y1 * sz;
@@ -463,7 +461,7 @@ void indme(
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;
+ dcomplex cri = (nsh % 2 == 0) ? 1.0 + I * 0.0 : c1->dc0[ic - 1] / ent;
for (int l190 = 1; l190 <= c1->li; l190++) {
int lpo = l190 + 1;
int ltpo = lpo + l190;
diff --git a/src/libnptm/sph_subs.cpp b/src/libnptm/sph_subs.cpp
index a3f91328..a67d2add 100644
--- a/src/libnptm/sph_subs.cpp
+++ b/src/libnptm/sph_subs.cpp
@@ -196,28 +196,28 @@ double cg1(int lmpml, int mu, int l, int m) {
return result;
}
-void diel(int npntmo, int ns, int i, int ic, double vk, ParticleDescriptor *c1, C2 *c2) {
+void diel(int npntmo, int ns, int i, int ic, double vk, ParticleDescriptor *c1) {
const double dif = c1->rc[i - 1][ns] - c1->rc[i - 1][ns - 1];
const double half_step = 0.5 * dif / npntmo;
double rr = c1->rc[i - 1][ns - 1];
- const dcomplex delta = c2->dc0[ic] - c2->dc0[ic - 1];
+ const dcomplex delta = c1->dc0[ic] - c1->dc0[ic - 1];
const int kpnt = npntmo + npntmo;
- c2->ris[kpnt] = c2->dc0[ic];
- c2->dlri[kpnt] = 0.0 + 0.0 * I;
+ c1->ris[kpnt] = c1->dc0[ic];
+ c1->dlri[kpnt] = 0.0 + 0.0 * I;
const int i90 = i - 1;
const int ns90 = ns - 1;
const int ic90 = ic - 1;
for (int np90 = 0; np90 < kpnt; np90++) {
double ff = (rr - c1->rc[i90][ns90]) / dif;
- c2->ris[np90] = delta * ff * ff * (-2.0 * ff + 3.0) + c2->dc0[ic90];
- c2->dlri[np90] = 3.0 * delta * ff * (1.0 - ff) / (dif * vk * c2->ris[np90]);
+ c1->ris[np90] = delta * ff * ff * (-2.0 * ff + 3.0) + c1->dc0[ic90];
+ c1->dlri[np90] = 3.0 * delta * ff * (1.0 - ff) / (dif * vk * c1->ris[np90]);
rr += half_step;
}
}
void dme(
int li, int i, int npnt, int npntts, double vk, double exdc, double exri,
- ParticleDescriptor *c1, C2 *c2, int &jer, int &lcalc, dcomplex &arg, int last_conf
+ ParticleDescriptor *c1, int &jer, int &lcalc, dcomplex &arg, int last_conf
) {
const int lipo = li + 1;
const int lipt = li + 2;
@@ -233,11 +233,11 @@ void dme(
int nstp = npnt - 1;
int nstpts = npntts - 1;
double sz = vk * c1->ros[sph_index - 1];
- c2->vsz[i - 1] = sz;
+ c1->vsz[i - 1] = sz;
double vkr1 = vk * c1->rc[sph_index - 1][0];
int nsh = c1->nshl[sph_index - 1];
- c2->vkt[i - 1] = csqrt(c2->dc0[0]);
- arg = vkr1 * c2->vkt[i - 1];
+ c1->vkt[i - 1] = csqrt(c1->dc0[0]);
+ arg = vkr1 * c1->vkt[i - 1];
arin = arg;
bool goto32 = false;
if (imag(arg) != 0.0) {
@@ -282,7 +282,7 @@ void dme(
fn[j43 - 1] = rfn[j43 - 1];
}
if (nsh <= 1) {
- cri = c2->dc0[0] / exdc;
+ cri = c1->dc0[0] / exdc;
for (int l60 = 1; l60 <= li; l60++) {
int lpo = l60 + 1;
int ltpo = lpo + l60;
@@ -305,7 +305,7 @@ void dme(
int lpt = lpo + 1;
int dltpo = ltpo;
y1 = fbi[lpo - 1];
- dy1 = ((1.0 * l80) * fbi[l80 - 1] - (1.0 * lpo) * fbi[lpt - 1]) * c2->vkt[i - 1] / (1.0 * dltpo);
+ dy1 = ((1.0 * l80) * fbi[l80 - 1] - (1.0 * lpo) * fbi[lpt - 1]) * c1->vkt[i - 1] / (1.0 * dltpo);
y2 = y1;
dy2 = dy1;
ic = 1;
@@ -316,13 +316,13 @@ void dme(
ic += 1;
double step = 1.0 * nstp;
step = vk * (c1->rc[i - 1][ns76 - 1] - c1->rc[i - 1][nsmo - 1]) / step;
- arg = c2->dc0[ic - 1];
+ arg = c1->dc0[ic - 1];
rkc(nstp, step, arg, vkr, lpo, y1, y2, dy1, dy2);
} else {
- diel(nstpts, nsmo, i, ic, vk, c1, c2);
+ diel(nstpts, nsmo, i, ic, vk, c1);
double stepts = 1.0 * nstpts;
stepts = vk * (c1->rc[i - 1][ns76 - 1] - c1->rc[i - 1][nsmo - 1]) / stepts;
- rkt(nstpts, stepts, vkr, lpo, y1, y2, dy1, dy2, c2);
+ rkt(nstpts, stepts, vkr, lpo, y1, y2, dy1, dy2, c1);
}
}
rmf[l80 - 1] = y1 * sz;
@@ -331,7 +331,7 @@ void dme(
dref[l80 - 1] = dy2 * sz + y2;
}
cri = 1.0 + uim * 0.0;
- if (nsh % 2 != 0) cri = c2->dc0[ic - 1] / exdc;
+ if (nsh % 2 != 0) cri = c1->dc0[ic - 1] / exdc;
for (int l90 = 1; l90 <= li; l90++) {
int lpo = l90 + 1;
int ltpo = lpo + l90;
@@ -346,7 +346,7 @@ void dme(
ccna = ref[l90 - 1] * dfn;
ccnb = dref[l90 - 1] * fn[lpo - 1] * (1.0 * sz * ltpo);
ccnc = ref[l90 - 1] * dfb;
- ccnd = dref[l90 - 1] *fb[lpo - 1] * (1.0 * sz * ltpo);
+ ccnd = dref[l90 - 1] * fb[lpo - 1] * (1.0 * sz * ltpo);
c1->rei[l90 - 1][i - 1] = 1.0 + uim * (cri * ccna - ccnb) / (cri * ccnc - ccnd);
}
} // nsh <= 1 ?
@@ -724,7 +724,7 @@ void rkc(
void rkt(
int npntmo, double step, double &x, int lpo, dcomplex &y1,
- dcomplex &y2, dcomplex &dy1, dcomplex &dy2, C2 *c2
+ dcomplex &y2, dcomplex &dy1, dcomplex &dy2, ParticleDescriptor *c1
) {
dcomplex cy1, cdy1, c11, cy23, cdy23, yc2, c12, c13;
dcomplex cy4, cdy4, yy, c14, c21, c22, c23, c24;
@@ -734,34 +734,34 @@ void rkt(
int ipnt = ipnt60 + 1;
int jpnt = ipnt + ipnt - 1;
int jpnt60 = jpnt - 1;
- cy1 = cl / (x * x) - c2->ris[jpnt60];
+ cy1 = cl / (x * x) - c1->ris[jpnt60];
cdy1 = -2.0 / x;
c11 = (cy1 * y1 + cdy1 * dy1) * step;
double xh = x + half_step;
int jpntpo = jpnt + 1;
- cy23 = cl / (xh * xh) - c2->ris[jpnt];
+ cy23 = cl / (xh * xh) - c1->ris[jpnt];
cdy23 = -2.0 / xh;
yc2 = y1 + dy1 * half_step;
c12 = (cy23 * yc2 + cdy23 * (dy1 + 0.5 * c11)) * step;
c13= (cy23 * (yc2 + 0.25 * c11 *step) + cdy23 * (dy1 + 0.5 * c12)) * step;
double xn = x + step;
//int jpntpt = jpnt + 2;
- cy4 = cl / (xn * xn) - c2->ris[jpntpo];
+ cy4 = cl / (xn * xn) - c1->ris[jpntpo];
cdy4 = -2.0 / xn;
yy = y1 + dy1 * step;
c14 = (cy4 * (yy + 0.5 * c12 * step) + cdy4 * (dy1 + c13)) * step;
y1= yy + (c11 + c12 + c13) * step / 6.0;
dy1 += (0.5 * c11 + c12 + c13 + 0.5 * c14) /3.0;
- cy1 -= cdy1 * c2->dlri[jpnt60];
- cdy1 += 2.0 * c2->dlri[jpnt60];
+ cy1 -= cdy1 * c1->dlri[jpnt60];
+ cdy1 += 2.0 * c1->dlri[jpnt60];
c21 = (cy1 * y2 + cdy1 * dy2) * step;
- cy23 -= cdy23 * c2->dlri[jpnt];
- cdy23 += 2.0 * c2->dlri[jpnt];
+ cy23 -= cdy23 * c1->dlri[jpnt];
+ cdy23 += 2.0 * c1->dlri[jpnt];
yc2 = y2 + dy2 * half_step;
c22 = (cy23 * yc2 + cdy23 * (dy2 + 0.5 * c21)) * step;
c23 = (cy23 * (yc2 + 0.25 * c21 * step) + cdy23 * (dy2 + 0.5 * c22)) * step;
- cy4 -= cdy4 * c2->dlri[jpntpo];
- cdy4 += 2.0 * c2->dlri[jpntpo];
+ cy4 -= cdy4 * c1->dlri[jpntpo];
+ cdy4 += 2.0 * c1->dlri[jpntpo];
yy = y2 + dy2 * step;
c24 = (cy4 * (yc2 + 0.5 * c22 * step) + cdy4 * (dy2 + c23)) * step;
y2 = yy + (c21 + c22 + c23) * step / 6.0;
diff --git a/src/sphere/sphere.cpp b/src/sphere/sphere.cpp
index e90fcae8..9f8c6616 100644
--- a/src/sphere/sphere.cpp
+++ b/src/sphere/sphere.cpp
@@ -147,7 +147,6 @@ void sphere(const string& config_file, const string& data_file, const string& ou
double sc_phi_start = gconf->sc_phi_start;
double sc_phi_step = gconf->sc_phi_step;
double sc_phi_end = gconf->sc_phi_end;
- C2 *c2 = new C2(gconf, sconf);
argi = new double[1];
args = new double[1];
gaps = new double[2];
@@ -317,18 +316,18 @@ void sphere(const string& config_file, const string& data_file, const string& ou
int ici = (nsh + 1) / 2;
if (idfc == 0) {
for (int ic = 0; ic < ici; ic++)
- c2->dc0[ic] = sconf->get_dielectric_constant(ic, i132, jxi488); // WARNING: IDFC=0 is not tested!
+ c1->dc0[ic] = sconf->get_dielectric_constant(ic, i132, jxi488); // WARNING: IDFC=0 is not tested!
} else { // IDFC != 0
if (jxi == 1) {
for (int ic = 0; ic < ici; ic++) {
- c2->dc0[ic] = sconf->get_dielectric_constant(ic, i132, jxi488);
+ c1->dc0[ic] = sconf->get_dielectric_constant(ic, i132, jxi488);
}
}
}
- if (nsh % 2 == 0) c2->dc0[ici] = exdc;
+ if (nsh % 2 == 0) c1->dc0[ici] = exdc;
int jer = 0;
int lcalc = 0;
- dme(l_max, i, npnt, npntts, vkarg, exdc, exri, c1, c2, jer, lcalc, arg);
+ dme(l_max, i, npnt, npntts, vkarg, exdc, exri, c1, jer, lcalc, arg);
if (jer != 0) {
fprintf(output, " STOP IN DME\n");
fprintf(output, " AT %1d LCALC=%3d TOO SMALL WITH ARG=%15.7lE+i(%15.7lE)\n", jer, lcalc, real(arg), imag(arg));
@@ -337,7 +336,6 @@ void sphere(const string& config_file, const string& data_file, const string& ou
delete sconf;
delete gconf;
delete c1;
- delete c2;
for (int zi = l_max - 1; zi > -1; zi--) {
for (int zj = 0; zj < 3; zj++) {
for (int zk = 0; zk < 2; zk++) {
@@ -409,14 +407,14 @@ void sphere(const string& config_file, const string& data_file, const string& ou
c1->sqexs[i170] *= sqsfi;
fprintf(output, " SPHERE %2d\n", i);
if (c1->nshl[i170] != 1) {
- fprintf(output, " SIZE=%15.7lE\n", c2->vsz[i170]);
+ fprintf(output, " SIZE=%15.7lE\n", c1->vsz[i170]);
} else {
fprintf(
output,
" SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n",
- c2->vsz[i170],
- real(c2->vkt[i170]),
- imag(c2->vkt[i170])
+ c1->vsz[i170],
+ real(c1->vkt[i170]),
+ imag(c1->vkt[i170])
);
}
fprintf(output, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
@@ -642,7 +640,6 @@ void sphere(const string& config_file, const string& data_file, const string& ou
}
fclose(output);
delete c1;
- delete c2;
for (int zi = l_max - 1; zi > -1; zi--) {
for (int zj = 0; zj < 3; zj++) {
for (int zk = 0; zk < 2; zk++) {
--
GitLab