Skip to content
Snippets Groups Projects
Commit 7dacec5c authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Separate sphere from cluster dependencies

parent 4e76fabf
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
src/cluster/cluster.cpp
+ 109
− 8
View file @ 7dacec5c
...@@ -94,18 +94,48 @@ void cluster() { ...@@ -94,18 +94,48 @@ void cluster() {
const int ndi = c4->nsph * c4->nlim; const int ndi = c4->nsph * c4->nlim;
C9 *c9 = new C9(ndi, c4->nlem, 2 * ndi, 2 * c4->nlem); C9 *c9 = new C9(ndi, c4->nlem, 2 * ndi, 2 * c4->nlem);
double *gaps = new double[nsph](); double *gaps = new double[nsph]();
double **tqse, **tqss; double **tqse, **tqss, **tqce, **tqcs;
complex<double> **tqspe, **tqsps; complex<double> **tqspe, **tqsps, **tqcpe, **tqcps;
tqse = new double*[2]; tqse = new double*[2];
tqspe = new complex<double>*[2]; tqspe = new complex<double>*[2];
tqss = new double*[2]; tqss = new double*[2];
tqsps = new complex<double>*[2]; tqsps = new complex<double>*[2];
tqce = new double*[2];
tqcpe = new complex<double>*[2];
tqcs = new double*[2];
tqcps = new complex<double>*[2];
for (int ti = 0; ti < 2; ti++) { for (int ti = 0; ti < 2; ti++) {
tqse[ti] = new double[nsph](); tqse[ti] = new double[nsph]();
tqspe[ti] = new complex<double>[nsph](); tqspe[ti] = new complex<double>[nsph]();
tqss[ti] = new double[nsph](); tqss[ti] = new double[nsph]();
tqsps[ti] = new complex<double>[nsph](); tqsps[ti] = new complex<double>[nsph]();
tqce[ti] = new double[3]();
tqcpe[ti] = new complex<double>[3]();
tqcs[ti] = new double[3]();
tqcps[ti] = new complex<double>[3]();
} }
complex<double> **gapp, **gappm;
double **gap, **gapm;
gapp = new complex<double>*[3];
gappm = new complex<double>*[3];
gap = new double*[3];
gapm = new double*[3];
for (int gi = 0; gi < 3; gi++) {
gapp[gi] = new complex<double>[2]();
gappm[gi] = new complex<double>[2]();
gap[gi] = new double[2]();
gapm[gi] = new double[2]();
}
double *u = new double[3]();
double *us = new double[3]();
double *un = new double[3]();
double *uns = new double[3]();
double *up = new double[3]();
double *ups = new double[3]();
double *unmp = new double[3]();
double *unsmp = new double[3]();
double *upmp = new double[3]();
double *upsmp = new double[3]();
// End of global variables for CLU // End of global variables for CLU
fprintf(output, " READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM\n"); fprintf(output, " READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM\n");
fprintf(output, " %5d%5d%5d%5d%5d%5d%5d%5d%5d\n", fprintf(output, " %5d%5d%5d%5d%5d%5d%5d%5d%5d\n",
...@@ -289,7 +319,8 @@ void cluster() { ...@@ -289,7 +319,8 @@ void cluster() {
} }
// label 160 // label 160
double cs0 = 0.25 * vk * vk * vk / acos(0.0); double cs0 = 0.25 * vk * vk * vk / acos(0.0);
double csch; double csch = 0.0, qschu = 0.0, pschu = 0.0, s0mag = 0.0;
std::complex<double> s0(0.0, 0.0);
scr0(vk, exri, c1, c1ao, c3, c4); scr0(vk, exri, c1, c1ao, c3, c4);
printf("DEBUG: after SCR0 TFSAS = (%lE, %lE)\n", c3->tfsas.real(), c3->tfsas.imag()); printf("DEBUG: after SCR0 TFSAS = (%lE, %lE)\n", c3->tfsas.real(), c3->tfsas.imag());
double sqk = vk * vk * sconf->exdc; double sqk = vk * vk * sconf->exdc;
...@@ -316,10 +347,10 @@ void cluster() { ...@@ -316,10 +347,10 @@ void cluster() {
fprintf(output, " %14.7lE%15.7lE%15.7lE\n", c1->sqscs[i], c1->sqabs[i], c1->sqexs[i]); fprintf(output, " %14.7lE%15.7lE%15.7lE\n", c1->sqscs[i], c1->sqabs[i], c1->sqexs[i]);
fprintf(output, " FSAS=%15.7lE%15.7lE\n", c1->fsas[i].real(), c1->fsas[i].imag()); fprintf(output, " FSAS=%15.7lE%15.7lE\n", c1->fsas[i].real(), c1->fsas[i].imag());
csch = 2.0 * vk * sqsfi / c1->gcsv[i]; csch = 2.0 * vk * sqsfi / c1->gcsv[i];
std::complex<double> s0 = c1->fsas[i] * exri; s0 = c1->fsas[i] * exri;
double qschu = s0.imag() * csch; qschu = s0.imag() * csch;
double pschu = s0.real() * csch; pschu = s0.real() * csch;
double s0mag = sqrt((s0.real() + s0.imag()) * (s0.real() - s0.imag())) * cs0; s0mag = sqrt((s0.real() + s0.imag()) * (s0.real() - s0.imag())) * cs0;
fprintf(output, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschu, pschu, s0mag); fprintf(output, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschu, pschu, s0mag);
double rapr = c1->sexs[i] - gaps[i]; double rapr = c1->sexs[i] - gaps[i];
double cosav = gaps[i] / c1->sscs[i]; double cosav = gaps[i] / c1->sscs[i];
...@@ -329,6 +360,48 @@ void cluster() { ...@@ -329,6 +360,48 @@ void cluster() {
} }
} // i170 loop } // i170 loop
fprintf(output, " FSAT=%15.7lE%15.7lE\n", c3->tfsas.real(), c3->tfsas.imag()); fprintf(output, " FSAT=%15.7lE%15.7lE\n", c3->tfsas.real(), c3->tfsas.imag());
csch = 2.0 * vk * sqsfi / c3->gcs;
s0 = c3->tfsas * exri;
qschu = s0.imag() * csch;
pschu = s0.real() * csch;
s0mag = sqrt((s0.real() + s0.imag()) * (s0.real() - s0.imag())) * cs0;
fprintf(output, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschu, pschu, s0mag);
tppoan.write(reinterpret_cast<char *>(&vk), sizeof(double));
pcrsm0(vk, exri, inpol, c1, c1ao, c4);
apcra(zpv, c4->le, c1ao->am0m, inpol, sqk, gapm, gappm);
th = th1;
for (int jth486 = 1; jth486 <= nth; jth486++) { // OpenMP portable?
ph = ph1;
double cost = 0.0, sint = 0.0, cosp = 0.0, sinp = 0.0;
// argi[NSPEF], with NSPEF=1 if IDOT=0, else NSPEF=NSPH
double *argi;
for (int jph484 = 1; jph484 <= nph; jph484++) {
int jw = 0;
if (nk != 1 || jxi488 <= 1) {
upvmp(th, ph, 0, cost, sint, cosp, sinp, u, upmp, unmp);
if (isam >= 0) {
argi = new double[1];
wmamp(
0, cost, sint, cosp, sinp, inpol, c4->le, 0,
nsph, argi, u, upmp, unmp, c1
);
// label 182
apc(zpv, c4->le, c1ao->am0m, c1->w, sqk, gap, gapp);
raba(c4->le, c1ao->am0m, c1->w, tqce, tqcpe, tqcs, tqcps);
jw = 1;
}
} else { // label 180, NK == 1 AND JXI488 == 1
if (isam >= 0) {
// label 182
apc(zpv, c4->le, c1ao->am0m, c1->w, sqk, gap, gapp);
raba(c4->le, c1ao->am0m, c1->w, tqce, tqcpe, tqcs, tqcps);
jw = 1;
}
}
// label 184
double thsl = ths1;
} // jph484 loop
} // jth486 loop
printf("INFO: done jxi488 iteration.\n"); printf("INFO: done jxi488 iteration.\n");
} // jxi488 loop } // jxi488 loop
tppoan.close(); tppoan.close();
...@@ -360,11 +433,39 @@ void cluster() { ...@@ -360,11 +433,39 @@ void cluster() {
delete[] tqss[ti]; delete[] tqss[ti];
delete[] tqspe[ti]; delete[] tqspe[ti];
delete[] tqsps[ti]; delete[] tqsps[ti];
delete[] tqce[ti];
delete[] tqcpe[ti];
delete[] tqcs[ti];
delete[] tqcps[ti];
} }
delete[] tqse; delete[] tqse;
delete[] tqss; delete[] tqss;
delete[] tqspe; delete[] tqspe;
delete[] tqsps; delete[] tqsps;
delete[] tqce;
delete[] tqcpe;
delete[] tqcs;
delete[] tqcps;
for (int gi = 2; gi > -1; gi--) {
delete[] gapp[gi];
delete[] gappm[gi];
delete[] gap[gi];
delete[] gapm[gi];
}
delete[] gapp;
delete[] gappm;
delete[] gap;
delete[] gapm;
delete[] u;
delete[] us;
delete[] un;
delete[] uns;
delete[] up;
delete[] ups;
delete[] unmp;
delete[] unsmp;
delete[] upmp;
delete[] upsmp;
} else { // NSPH mismatch between geometry and scatterer configurations. } else { // NSPH mismatch between geometry and scatterer configurations.
throw UnrecognizedConfigurationException( throw UnrecognizedConfigurationException(
"Inconsistent geometry and scatterer configurations." "Inconsistent geometry and scatterer configurations."
......
src/include/Commons.h
+ 4
− 0
View file @ 7dacec5c
...@@ -238,6 +238,10 @@ public: ...@@ -238,6 +238,10 @@ public:
//! \brief QUESTION: definition? //! \brief QUESTION: definition?
std::complex<double> *scscp; std::complex<double> *scscp;
//! \brief QUESTION: definition? //! \brief QUESTION: definition?
double *ecscm;
//! \brief QUESTION: definition?
double *scscm;
//! \brief QUESTION: definition?
std::complex<double> *ecscp; std::complex<double> *ecscp;
//! \brief QUESTION: definition? //! \brief QUESTION: definition?
std::complex<double> *scscpm; std::complex<double> *scscpm;
......
src/include/clu_subs.h
+ 585
− 1
View file @ 7dacec5c
...@@ -23,6 +23,7 @@ ...@@ -23,6 +23,7 @@
// >>> DECLARATION OF SPH_SUBS <<< // >>> DECLARATION OF SPH_SUBS <<<
extern void aps(double ****zpv, int li, int nsph, C1 *c1, double sqk, double *gaps); extern void aps(double ****zpv, int li, int nsph, C1 *c1, double sqk, double *gaps);
extern std::complex<double> dconjg(std::complex<double> value); extern std::complex<double> dconjg(std::complex<double> value);
extern double cg1(int lmpml, int mu, int l, int m);
extern void dme( extern void dme(
int li, int i, int npnt, int npntts, double vk, double exdc, double exri, int li, int i, int npnt, int npntts, double vk, double exdc, double exri,
C1 *c1, C2 *c2, int &jer, int &lcalc, std::complex<double> &arg C1 *c1, C2 *c2, int &jer, int &lcalc, std::complex<double> &arg
...@@ -35,6 +36,15 @@ extern void rbf(int n, double x, int &nm, double sj[]); ...@@ -35,6 +36,15 @@ extern void rbf(int n, double x, int &nm, double sj[]);
extern void rnf(int n, double x, int &nm, double sy[]); extern void rnf(int n, double x, int &nm, double sy[]);
extern void sphar(double cth, double sth, double cph, double sph, int lm, std::complex<double> *ylm); extern void sphar(double cth, double sth, double cph, double sph, int lm, std::complex<double> *ylm);
extern void thdps(int lm, double ****zpv); extern void thdps(int lm, double ****zpv);
extern void upvmp(
double thd, double phd, int icspnv, double &cost, double &sint,
double &cosp, double &sinp, double *u, double *up, double *un
);
extern void wmamp(
int iis, double cost, double sint, double cosp, double sinp, int inpol,
int lm, int idot, int nsph, double *arg, double *u, double *up,
double *un, C1 *c1
);
// >>> END OF SPH_SUBS DECLARATION <<< // >>> END OF SPH_SUBS DECLARATION <<<
void logmat(std::complex<double> **mat, int rows, int cols); void logmat(std::complex<double> **mat, int rows, int cols);
...@@ -235,7 +245,6 @@ void r3jjr(int j2, int j3, int m2, int m3, C6 *c6) { ...@@ -235,7 +245,6 @@ void r3jjr(int j2, int j3, int m2, int m3, C6 *c6) {
} }
} }
/*! \brief C++ porting of GHIT /*! \brief C++ porting of GHIT
* *
* \param ihi: `int` * \param ihi: `int`
...@@ -463,6 +472,365 @@ std::complex<double> ghit( ...@@ -463,6 +472,365 @@ std::complex<double> ghit(
return result; return result;
} }
/*! \brief C++ porting of APC
*
* \param zpv: `double ****`
* \param le: `int`
* \param am0m: Matrix of complex.
* \param w: Matrix of complex.
* \param sqk: `double`
* \param gapr: `double **`
* \param gapp: Matrix of complex.
*/
void apc(
double ****zpv, int le, std::complex<double> **am0m, std::complex<double> **w,
double sqk, double **gapr, std::complex<double> **gapp
) {
std::complex<double> **ac, **gap;
const std::complex<double> cc0(0.0, 0.0);
const std::complex<double> uim(0.0, 1.0);
std::complex<double> uimmp, summ, sume, suem, suee, summp, sumep;
std::complex<double> suemp, sueep;
double cof = 1.0 / sqk;
double cimu = cof / sqrt(2.0);
int nlem = le * (le + 2);
const int nlemt = nlem + nlem;
ac = new std::complex<double>*[nlemt];
gap = new std::complex<double>*[3];
for (int ai = 0; ai < nlemt; ai++) ac[ai] = new std::complex<double>[2]();
for (int gi = 0; gi < 3; gi++) gap[gi] = new std::complex<double>[2]();
for (int j45 = 1; j45 <= nlemt; j45++) {
int j = j45 - 1;
ac[j][0] = cc0;
ac[j][1] = cc0;
for (int i45 = 1; i45 <= nlemt; i45++) {
int i = i45 - 1;
ac[j][0] += (am0m[j][i] * w[i][0]);
ac[j][1] += (am0m[j][i] * w[i][1]);
} //i45 loop
} //j45 loop
for (int imu90 = 1; imu90 <=3; imu90++) {
int mu = imu90 - 2;
gap[imu90 - 1][0] = cc0;
gap[imu90 - 1][1] = cc0;
gapp[imu90 - 1][0] = cc0;
gapp[imu90 - 1][1] = cc0;
for (int l80 =1; l80 <= le; l80++) {
int lpo = l80 + 1;
int ltpo = lpo + l80;
int imm = l80 * lpo;
for (int ilmp = 1; ilmp <= 3; ilmp++) {
if ((l80 == 1 && ilmp == 1) || (l80 == le && ilmp == 3)) continue; // ilmp loop
int lmpml = ilmp - 2;
int lmp = l80 + lmpml;
uimmp = (-1.0 * lmpml) * uim;
int impmmmp = lmp * (lmp + 1);
for (int im70 = 1; im70 <= ltpo; im70++) {
int m = im70 - lpo;
int mmp = m - mu;
int abs_mmp = (mmp > 0) ? mmp : -mmp;
if (abs_mmp <= lmp) {
int i = imm + m;
int ie = i + nlem;
int imp = impmmmp + mmp;
int impe = imp + nlem;
double cgc = cg1(lmpml, mu, l80, m);
int jpo = 2;
for (int ipo = 1; ipo <= 2; ipo++) {
if (ipo == 2) jpo = 1;
//printf("DEBUG: i=%d, ipo=%d, imp=%d\n", i, ipo, imp);
//fflush(stdout);
summ = dconjg(ac[i - 1][ipo - 1]) * ac[imp - 1][ipo - 1];
sume = dconjg(ac[i - 1][ipo - 1]) * ac[impe - 1][ipo - 1];
suem = dconjg(ac[ie - 1][ipo - 1]) * ac[imp - 1][ipo - 1];
suee = dconjg(ac[ie - 1][ipo - 1]) * ac[impe - 1][ipo - 1];
summp = dconjg(ac[i - 1][jpo - 1]) * ac[imp - 1][ipo - 1];
sumep = dconjg(ac[i - 1][jpo - 1]) * ac[impe - 1][ipo - 1];
suemp = dconjg(ac[ie - 1][jpo - 1]) * ac[imp - 1][ipo - 1];
sueep = dconjg(ac[ie - 1][jpo - 1]) * ac[impe - 1][ipo - 1];
if (lmpml != 0) {
summ *= uimmp;
sume *= uimmp;
suem *= uimmp;
suee *= uimmp;
summp *= uimmp;
sumep *= uimmp;
suemp *= uimmp;
sueep *= uimmp;
}
// label 55
gap[imu90 - 1][ipo - 1] += (
(
summ * zpv[l80 - 1][ilmp - 1][0][0]
+ sume * zpv[l80 - 1][ilmp - 1][0][1]
+ suem * zpv[l80 - 1][ilmp - 1][1][0]
+ suee * zpv[l80 - 1][ilmp - 1][1][1]
) * cgc
);
gapp[imu90 - 1][ipo - 1] += (
(
summp * zpv[l80 - 1][ilmp - 1][0][0]
+ sumep * zpv[l80 - 1][ilmp - 1][0][1]
+ suemp * zpv[l80 - 1][ilmp - 1][1][0]
+ sueep * zpv[l80 - 1][ilmp - 1][1][1]
) * cgc
);
} // ipo loop
} // ends im70 loop
} // im70 loop
} // ilmp loop
} // l80 loop
} // imu90 loop
for (int ipo95 = 1; ipo95 <= 2; ipo95++) {
sume = gap[0][ipo95 - 1] * cimu;
suee = gap[1][ipo95 - 1] * cof;
suem = gap[2][ipo95 - 1] * cimu;
gapr[0][ipo95 - 1] = (sume - suem).real();
gapr[1][ipo95 - 1] = ((sume + suem) * uim).real();
gapr[2][ipo95 - 1] = suee.real();
sumep = gapp[0][ipo95 - 1] * cimu;
sueep = gapp[1][ipo95 - 1] * cof;
suemp = gapp[2][ipo95 - 1] * cimu;
gapp[0][ipo95 - 1] = sumep - suemp;
gapp[1][ipo95 - 1] = (sumep + suemp) * uim;
gapp[2][ipo95 - 1] = sueep;
} // ipo95 loop
// Clean memory
for (int ai = nlemt - 1; ai > -1; ai--) delete[] ac[ai];
for (int gi = 2; gi > -1; gi--) delete[] gap[gi];
delete[] ac;
delete[] gap;
}
/*! \brief C++ porting of APCRA
*
* \param zpv: `double ****`
* \param le: `int`
* \param am0m: Matrix of complex.
* \param inpol: `int` Polarization type.
* \param sqk: `double`
* \param gaprm: `double **`
* \param gappm: Matrix of complex.
*/
void apcra(
double ****zpv, const int le, std::complex<double> **am0m, int inpol, double sqk,
double **gaprm, std::complex<double> **gappm
) {
const std::complex<double> cc0(0.0, 0.0);
const std::complex<double> uim(0.0, 1.0);
std::complex<double> uimtl, uimtls, ca11, ca12, ca21, ca22;
std::complex<double> a11, a12, a21, a22, sum1, sum2, fc;
double ****svw = new double***[le];
std::complex<double> ****svs = new std::complex<double>***[le];
for (int i = 0; i < le; i++) {
svw[i] = new double**[3];
svs[i] = new std::complex<double>**[3];
for (int j = 0; j < 3; j++) {
svw[i][j] = new double*[2];
svs[i][j] = new std::complex<double>*[2];
for (int k = 0; k < 2; k++) {
svw[i][j][k] = new double[2]();
svs[i][j][k] = new std::complex<double>[2]();
}
}
}
int nlem = le * (le + 2);
for (int l28 = 1; l28 <= le; l28++) {
int lpo = l28 + 1;
int ltpo = lpo + l28;
double fl = sqrt(1.0 * ltpo);
for (int ilmp = 1; ilmp <= 3; ilmp++) {
if ((l28 == 1 && ilmp == 1) || (l28 == le && ilmp == 3)) continue; // ilmp loop
int lmpml = ilmp - 2;
int lmp = l28 + lmpml;
double flmp = sqrt(1.0 * (lmp + lmp + 1));
double fllmp = flmp / fl;
double cgmmo = fllmp * cg1(lmpml, 0, l28, 1);
double cgmpo = fllmp * cg1(lmpml, 0, l28, -1);
if (inpol == 0) {
double cgs = cgmpo + cgmmo;
double cgd = cgmpo - cgmmo;
svw[l28 - 1][ilmp - 1][0][0] = cgs;
svw[l28 - 1][ilmp - 1][0][1] = cgd;
svw[l28 - 1][ilmp - 1][1][0] = cgd;
svw[l28 - 1][ilmp - 1][1][1] = cgs;
} else { // label 22
svw[l28 - 1][ilmp - 1][0][0] = cgmpo;
svw[l28 - 1][ilmp - 1][1][0] = cgmpo;
svw[l28 - 1][ilmp - 1][0][1] = -cgmmo;
svw[l28 - 1][ilmp - 1][1][1] = cgmmo;
}
// label 26
} // ilmp loop
} // l28 loop
for (int l30 = 1; l30 <= le; l30++) { // This can be omitted, since svs was initialized at 0
for (int ilmp = 1; ilmp <= 3; ilmp++) {
for (int ipa = 1; ipa <= 2; ipa++) {
for (int ipamp = 1; ipamp <= 2; ipamp++) {
svs[l30 - 1][ilmp - 1][ipa - 1][ipamp - 1] = cc0;
}
} // ipa loop
} // ilmp loop
} // l30 loop
for (int l58 = 1; l58 <= le; l58 ++) {
int lpo = l58 + 1;
int ltpo = l58 + lpo;
int imm = l58 * lpo;
for (int ilmp = 1; ilmp <= 3; ilmp++) {
if ((l58 == 1 && ilmp == 1) || (l58 == le && ilmp == 3)) continue; // ilmp loop
int lmpml = ilmp - 2;
int lmp = l58 + lmpml;
int impmm = lmp * (lmp + 1);
uimtl = uim * (1.0 * lmpml);
if (lmpml == 0) uimtl = std::complex<double>(1.0, 0.0);
for (int im54 = 1; im54 <= ltpo; im54++) {
int m = im54 - lpo;
int i = imm + m;
int ie = i + nlem;
for (int imu52 = 1; imu52 <= 3; imu52++) {
int mu = imu52 - 2;
int mmp = m - mu;
int abs_mmp = (mmp > 0) ? mmp : -mmp;
if (abs_mmp <= lmp) {
int imp = impmm + mmp;
int impe = imp + nlem;
double cgc = cg1(lmpml, -mu, l58, -m);
for (int ls = 1; ls <= le; ls++) {
int lspo = ls + 1;
int lstpo = ls + lspo;
int ismm = ls * lspo;
for (int ilsmp = 1; ilsmp <= 3; ilsmp++) {
if ((ls == 1 && ilsmp == 1) || (ls == le && ilsmp == 3)) continue; // ilsmp loop
int lsmpml = ilsmp - 2;
int lsmp = ls + lsmpml;
int ismpmm = lsmp * (lsmp + 1);
uimtls = -uim * (1.0 * lsmpml);
if (lsmpml == 0) uimtls = std::complex<double>(1.0, 0.0);
for (int ims = 1; ims <= lstpo; ims++) {
int ms = ims - lspo;
int msmp = ms - mu;
int abs_msmp = (msmp > 0) ? msmp : -msmp;
if (abs_msmp <= lsmp) {
int is = ismm + ms;
int ise = is + nlem;
int ismp = ismpmm + msmp;
int ismpe = ismp + nlem;
double cgcs = cg1(lsmpml, mu, ls, ms);
fc = (uimtl * uimtls) * (cgc * cgcs);
ca11 = dconjg(am0m[is - 1][i - 1]);
ca12 = dconjg(am0m[is - 1][ie - 1]);
ca21 = dconjg(am0m[ise - 1][i - 1]);
ca22 = dconjg(am0m[ise - 1][ie - 1]);
a11 = am0m[ismp - 1][imp - 1];
a12 = am0m[ismp - 1][impe - 1];
a21 = am0m[ismpe - 1][imp - 1];
a22 = am0m[ismpe - 1][impe - 1];
double z11 = zpv[ls - 1][ilsmp - 1][0][0];
double z12 = zpv[ls - 1][ilsmp - 1][0][1];
double z21 = zpv[ls - 1][ilsmp - 1][1][0];
double z22 = zpv[ls - 1][ilsmp - 1][1][1];
svs[l58 - 1][ilmp - 1][0][0] += ((ca11 * a11 * z11
+ ca11 * a21 * z12
+ ca21 * a11 * z21
+ ca21 * a21 * z22) * fc);
svs[l58 - 1][ilmp - 1][0][1] += ((ca11 * a12 * z11
+ ca11 * a22 * z12
+ ca21 * a12 * z21
+ ca21 * a22 * z22) * fc);
svs[l58 - 1][ilmp - 1][1][0] += ((ca12 * a11 * z11
+ ca12 * a21 * z12
+ ca22 * a11 * z21
+ ca22 * a21 * z22) * fc);
svs[l58 - 1][ilmp - 1][1][1] += ((ca12 * a12 * z11
+ ca12 * a22 * z12
+ ca22 * a12 * z21
+ ca22 * a22 * z22) * fc);
} // ends ims loop
} // ims loop
} // ilsmp loop
} // ls loop
} // ends imu52 loop
} // imu52 loop
} // im54 loop
} // ilmp loop
} // l58 loop
sum1 = cc0;
sum2 = cc0;
for (int l68 = 1; l68 <= le; l68++) {
int lpo = l68 + 1;
int ltpo = l68 + lpo;
int imm = l68 * lpo;
for (int ilmp = 1; ilmp <= 3; ilmp++) {
if ((l68 == 1 && ilmp == 1) || (l68 == le && ilmp == 3)) continue; // ilmp loop
if (inpol == 0) {
sum1 += (
svw[l68 - 1][ilmp - 1][0][0] * svs[l68 - 1][ilmp - 1][0][0]
+ svw[l68 - 1][ilmp - 1][1][0] * svs[l68 - 1][ilmp - 1][0][1]
+ svw[l68 - 1][ilmp - 1][1][0] * svs[l68 - 1][ilmp - 1][1][0]
+ svw[l68 - 1][ilmp - 1][0][0] * svs[l68 - 1][ilmp - 1][1][1]
);
sum2 += (
svw[l68 - 1][ilmp - 1][0][1] * svs[l68 - 1][ilmp - 1][0][0]
+ svw[l68 - 1][ilmp - 1][1][1] * svs[l68 - 1][ilmp - 1][0][1]
+ svw[l68 - 1][ilmp - 1][1][1] * svs[l68 - 1][ilmp - 1][1][0]
+ svw[l68 - 1][ilmp - 1][0][1] * svs[l68 - 1][ilmp - 1][1][1]
);
} else { // label 62
sum1 += (
svw[l68 - 1][ilmp - 1][1][0] * svs[l68 - 1][ilmp - 1][0][0]
+ svw[l68 - 1][ilmp - 1][0][0] * svs[l68 - 1][ilmp - 1][0][1]
+ svw[l68 - 1][ilmp - 1][0][0] * svs[l68 - 1][ilmp - 1][1][0]
+ svw[l68 - 1][ilmp - 1][1][0] * svs[l68 - 1][ilmp - 1][1][1]
);
sum2 += (
svw[l68 - 1][ilmp - 1][1][1] * svs[l68 - 1][ilmp - 1][0][0]
+ svw[l68 - 1][ilmp - 1][0][1] * svs[l68 - 1][ilmp - 1][0][1]
+ svw[l68 - 1][ilmp - 1][0][1] * svs[l68 - 1][ilmp - 1][1][0]
+ svw[l68 - 1][ilmp - 1][1][1] * svs[l68 - 1][ilmp - 1][1][1]
);
} // label 66, ends ilmp loop
} // ilmp loop
} // l68 loop
const double half_pi = acos(0.0);
double cofs = half_pi * 2.0 / sqk;
gaprm[0][0] = 0.0;
gaprm[0][1] = 0.0;
gaprm[1][0] = 0.0;
gaprm[1][1] = 0.0;
gappm[0][0] = cc0;
gappm[0][1] = cc0;
gappm[1][0] = cc0;
gappm[1][1] = cc0;
if (inpol == 0) {
sum1 *= cofs;
sum2 *= cofs;
gaprm[2][0] = sum1.real();
gaprm[2][1] = sum1.real();
gappm[2][0] = sum2 * uim;
gappm[2][1] = -gappm[2][0];
} else { // label 72
gaprm[2][0] = sum1.real() * cofs;
gaprm[2][1] = sum2.real() * cofs;
gappm[2][0] = cc0;
gappm[2][1] = cc0;
}
// Clean memory
for (int i = le - 1; i > -1; i--) {
for (int j = 2; j > -1; j--) {
for (int k = 1; k > -1; k--) {
delete[] svw[i][j][k];
delete[] svs[i][j][k];
}
delete[] svw[i][j];
delete[] svs[i][j];
}
delete[] svw[i];
delete[] svs[i];
}
delete[] svw;
delete[] svs;
}
/*! \brief C++ porting of CMS /*! \brief C++ porting of CMS
* *
* \param am: Matrix of complex. * \param am: Matrix of complex.
...@@ -746,6 +1114,81 @@ void lucin(std::complex<double> **am, const int nddmst, int n, int &ier) { ...@@ -746,6 +1114,81 @@ void lucin(std::complex<double> **am, const int nddmst, int n, int &ier) {
delete[] v; delete[] v;
} }
/*! \brief C++ porting of PCRSM0
*
* \param vk: `double`
* \param exri: `double`
* \param inpol: `int`
* \param c1: `C1 *`
* \param c1ao: `C1_AddOns *`
* \param c4: `C4 *`
*/
void pcrsm0(double vk, double exri, int inpol, C1 *c1, C1_AddOns *c1ao, C4 *c4) {
const std::complex<double> cc0(0.0, 0.0);
const std::complex<double> uim(0.0, 1.0);
std::complex<double> sum1, sum2, sum3, sum4, sumpd;
std::complex<double> sums1, sums2, sums3, sums4, csam;
double exdc = exri * exri;
double ccs = 4.0 * acos(0.0) / (vk * vk);
double cccs = ccs / exdc;
csam = -(ccs / (exri * vk)) * std::complex<double>(0.0, 0.5);
sum2 = cc0;
sum3 = cc0;
for (int i14 = 1; i14 <= c4->nlem; i14++) { // GPU portable?
int ie = i14 + c4->nlem;
sum2 += (c1ao->am0m[i14 - 1][i14 - 1] + c1ao->am0m[ie - 1][ie - 1]);
sum3 += (c1ao->am0m[i14 - 1][ie - 1] + c1ao->am0m[ie - 1][i14 - 1]);
} // i14 loop
double sumpi = 0.0;
sumpd = cc0;
int nlemt = c4->nlem + c4->nlem;
for (int i16 = 1; i16 <= nlemt; i16++) {
for (int j16 = 1; j16 <= c4->nlem; j16++) {
int je = j16 + c4->nlem;
double rvalue = (
dconjg(c1ao->am0m[i16 - 1][j16 - 1]) * c1ao->am0m[i16 - 1][j16 - 1]
+ dconjg(c1ao->am0m[i16 - 1][je - 1]) * c1ao->am0m[i16 - 1][je - 1]
).real();
sumpi += rvalue;
sumpd += (
dconjg(c1ao->am0m[i16 - 1][j16 - 1]) * c1ao->am0m[i16 - 1][je - 1]
+ dconjg(c1ao->am0m[i16 - 1][je - 1]) * c1ao->am0m[i16 - 1][j16 - 1]
);
} // j16 loop
} // i16 loop
if (inpol == 0) {
sum1 = sum2;
sum4 = sum3 * uim;
sum3 = -sum4;
sums1 = sumpi;
sums2 = sumpi;
sums3 = sumpd * uim;
sums4 = -sums3;
} else { // label 18
sum1 = sum2 + sum3;
sum2 = sum2 - sum3;
sum3 = cc0;
sum4 = cc0;
sums1 = sumpi - sumpd;
sums2 = sumpi + sumpd;
sums3 = cc0;
sums4 = cc0;
}
// label 20
c1ao->ecscm[0] = -cccs * sum2.real();
c1ao->ecscm[1] = -cccs * sum1.real();
c1ao->ecscpm[0] = -cccs * sum4;
c1ao->ecscpm[1] = -cccs * sum3;
c1ao->fsacm[0][0] = csam * sum2;
c1ao->fsacm[1][0] = csam * sum4;
c1ao->fsacm[1][1] = csam * sum1;
c1ao->fsacm[0][1] = csam * sum3;
c1ao->scscm[0] = cccs * sums1.real();
c1ao->scscm[1] = cccs * sums2.real();
c1ao->scscpm[0] = cccs * sums3;
c1ao->scscpm[1] = cccs * sums4;
}
/*! \brief C++ porting of POLAR /*! \brief C++ porting of POLAR
* *
* \param x: `double` * \param x: `double`
...@@ -915,6 +1358,147 @@ void r3j000(int j2, int j3, C6 *c6) { ...@@ -915,6 +1358,147 @@ void r3j000(int j2, int j3, C6 *c6) {
} }
} }
/*! \brief C++ porting of RABA
*
* \param le: `int`
* \param am0m: Matrix of complex.
* \param w: Matrix of complex.
* \param tqce: `double **`
* \param tqcpe: Matrix of complex.
* \param tqcs: `double **`
* \param tqcps: Matrix of complex.
*/
void raba(
int le, std::complex<double> **am0m, std::complex<double> **w, double **tqce,
std::complex<double> **tqcpe, double **tqcs, std::complex<double> **tqcps
) {
std::complex<double> **a, **ctqce, **ctqcs;
std::complex<double> acw, acwp, aca, acap, c1, c2, c3;
const std::complex<double> cc0(0.0, 0.0);
const std::complex<double> uim(0.0, 1.0);
const double sq2i = 1.0 / sqrt(2.0);
int nlem = le * (le + 2);
const int nlemt = nlem + nlem;
a = new std::complex<double>*[nlemt];
ctqce = new std::complex<double>*[2];
ctqcs = new std::complex<double>*[2];
for (int ai = 0; ai < nlemt; ai++) a[ai] = new std::complex<double>[2]();
for (int ci = 0; ci < 2; ci++) {
ctqce[ci] = new std::complex<double>[3]();
ctqcs[ci] = new std::complex<double>[3]();
}
for (int i20 = 1; i20 <= nlemt; i20++) {
int i = i20 - 1;
c1 = cc0;
c2 = cc0;
for (int j10 = 1; j10 <= nlemt; j10++) {
int j = j10 - 1;
c1 += (am0m[i][j] * w[j][0]);
c2 += (am0m[i][j] * w[j][1]);
} // j10 loop
a[i][0] = c1;
a[i][1] = c2;
} //i20 loop
int jpo = 2;
for (int ipo70 = 1; ipo70 <= 2; ipo70++) {
if (ipo70 == 2) jpo = 1;
int ipo = ipo70 - 1;
ctqce[ipo][0] = cc0;
ctqce[ipo][1] = cc0;
ctqce[ipo][2] = cc0;
tqcpe[ipo][0] = cc0;
tqcpe[ipo][1] = cc0;
tqcpe[ipo][2] = cc0;
ctqcs[ipo][0] = cc0;
ctqcs[ipo][1] = cc0;
ctqcs[ipo][2] = cc0;
tqcps[ipo][0] = cc0;
tqcps[ipo][1] = cc0;
tqcps[ipo][2] = cc0;
for (int l60 = 1; l60 <= le; l60 ++) {
int lpo = l60 + 1;
int il = l60 * lpo;
int ltpo = l60 + lpo;
for (int im60 = 1; im60 <= ltpo; im60++) {
int m = im60 - lpo;
int i = m + il;
int ie = i + nlem;
int mmmu = m + 1;
int mmmmu = (mmmu > 0) ? mmmu : -mmmu;
double rmu = 0.0;
if (mmmmu <= l60) {
int immu = mmmu + il;
int immue = immu + nlem;
rmu = -sqrt(1.0 * (l60 + mmmu) * (l60 - m)) * sq2i;
acw = dconjg(a[i - 1][ipo]) * w[immu - 1][ipo] + dconjg(a[ie - 1][ipo]) * w[immue - 1][ipo];
acwp = dconjg(a[i - 1][ipo]) * w[immu - 1][jpo - 1] + dconjg(a[ie - 1][ipo]) * w[immue - 1][jpo - 1];
aca = dconjg(a[i - 1][ipo]) * a[immu - 1][ipo] + dconjg(a[ie - 1][ipo]) * a[immue - 1][ipo];
acap = dconjg(a[i - 1][ipo]) * a[immu - 1][jpo - 1] + dconjg(a[ie - 1][ipo]) * a[immue - 1][jpo - 1];
ctqce[ipo][0] += (acw * rmu);
tqcpe[ipo][0] += (acwp * rmu);
ctqcs[ipo][0] += (aca * rmu);
tqcps[ipo][0] += (acap * rmu);
}
// label 30
rmu = -1.0 * m;
acw = dconjg(a[i - 1][ipo]) * w[i - 1][ipo] + dconjg(a[ie - 1][ipo]) * w[ie - 1][ipo];
acwp = dconjg(a[i - 1][ipo]) * w[i - 1][jpo - 1] + dconjg(a[ie - 1][ipo]) * w[ie - 1][jpo - 1];
aca = dconjg(a[i - 1][ipo]) * a[i - 1][ipo] + dconjg(a[ie - 1][ipo]) * a[ie - 1][ipo];
acap = dconjg(a[i - 1][ipo]) * a[i - 1][jpo - 1] + dconjg(a[ie - 1][ipo]) * a[ie - 1][jpo - 1];
ctqce[ipo][1] += (acw * rmu);
tqcpe[ipo][1] += (acwp * rmu);
ctqcs[ipo][1] += (aca * rmu);
tqcps[ipo][1] += (acap * rmu);
mmmu = m - 1;
mmmmu = (mmmu > 0) ? mmmu : -mmmu;
if (mmmmu <= l60) {
int immu = mmmu + il;
int immue = immu + nlem;
rmu = sqrt(1.0 * (l60 - mmmu) * (l60 + m)) * sq2i;
acw = dconjg(a[i - 1][ipo]) * w[immu - 1][ipo] + dconjg(a[ie - 1][ipo]) * w[immue - 1][ipo];
acwp = dconjg(a[i - 1][ipo]) * w[immu - 1][jpo - 1] + dconjg(a[ie - 1][ipo]) * w[immue - 1][jpo - 1];
aca = dconjg(a[i - 1][ipo]) * a[immu - 1][ipo] + dconjg(a[ie - 1][ipo]) * a[immue - 1][ipo];
acap = dconjg(a[i - 1][ipo]) * a[immu - 1][jpo - 1] + dconjg(a[ie - 1][ipo]) * a[immue - 1][jpo - 1];
ctqce[ipo][2] += (acw * rmu);
tqcpe[ipo][2] += (acwp * rmu);
ctqcs[ipo][2] += (aca * rmu);
tqcps[ipo][2] += (acap * rmu);
} // ends im60 loop
} // im60 loop
} // l60 loop
} // ipo70 loop
for (int ipo78 = 1; ipo78 <= 2; ipo78++) {
int ipo = ipo78 - 1;
tqce[ipo][0] = (ctqce[ipo][0] - ctqce[ipo][2]).real() * sq2i;
tqce[ipo][1] = ((ctqce[ipo][0] + ctqce[ipo][2]) * uim).real() * sq2i;
tqce[ipo][2] = ctqce[ipo][1].real();
c1 = tqcpe[ipo][0];
c2 = tqcpe[ipo][1];
c3 = tqcpe[ipo][2];
tqcpe[ipo][0] = (c1 - c3) * sq2i;
tqcpe[ipo][1] = (c1 + c3) * (uim * sq2i);
tqcpe[ipo][2] = c2;
tqcs[ipo][0] = -sq2i * (ctqcs[ipo][0] - ctqcs[ipo][2]).real();
tqcs[ipo][1] = -sq2i * ((ctqcs[ipo][0] + ctqcs[ipo][2]) * uim).real();
tqcs[ipo][2] = -1.0 * ctqcs[ipo][1].real();
c1 = tqcps[ipo][0];
c2 = tqcps[ipo][1];
c3 = tqcps[ipo][2];
tqcps[ipo][0] = -(c1 - c3) * sq2i;
tqcps[ipo][1] = -(c1 + c3) * (uim * sq2i);
tqcps[ipo][2] = -c2;
} // ipo78 loop
// Clean memory
for (int ai = 0; ai < nlemt; ai++) delete[] a[ai];
for (int ci = 0; ci < 2; ci++) {
delete[] ctqce[ci];
delete[] ctqcs[ci];
}
delete[] a;
delete[] ctqce;
delete[] ctqcs;
}
/*! \brief C++ porting of SCR0 /*! \brief C++ porting of SCR0
* *
* \param vk: `double` QUESTION: definition? * \param vk: `double` QUESTION: definition?
......
src/libnptm/Commons.cpp
+ 4
− 0
View file @ 7dacec5c
...@@ -123,6 +123,8 @@ C1_AddOns::C1_AddOns(C4 *c4) { ...@@ -123,6 +123,8 @@ C1_AddOns::C1_AddOns(C4 *c4) {
ecscpm = new complex<double>[2](); ecscpm = new complex<double>[2]();
allocate_vectors(c4); allocate_vectors(c4);
sscs = new double[nsph](); sscs = new double[nsph]();
ecscm = new double[2]();
scscm = new double[2]();
} }
C1_AddOns::~C1_AddOns() { C1_AddOns::~C1_AddOns() {
...@@ -153,6 +155,8 @@ C1_AddOns::~C1_AddOns() { ...@@ -153,6 +155,8 @@ C1_AddOns::~C1_AddOns() {
delete[] ecscp; delete[] ecscp;
delete[] scscpm; delete[] scscpm;
delete[] ecscpm; delete[] ecscpm;
delete[] ecscm;
delete[] scscm;
} }
void C1_AddOns::allocate_vectors(C4 *c4) { void C1_AddOns::allocate_vectors(C4 *c4) {
......
src/np_cluster.cpp
+ 1
− 4
View file @ 7dacec5c
...@@ -10,7 +10,6 @@ ...@@ -10,7 +10,6 @@
using namespace std; using namespace std;
extern void cluster(); extern void cluster();
extern void sphere();
/*! \brief Main program entry point. /*! \brief Main program entry point.
* *
...@@ -21,8 +20,6 @@ extern void sphere(); ...@@ -21,8 +20,6 @@ extern void sphere();
* the configuration and runs the main program. * the configuration and runs the main program.
*/ */
int main(int argc, char **argv) { int main(int argc, char **argv) {
bool is_sphere = false; cluster();
if (is_sphere) sphere();
else cluster();
return 0; return 0;
} }
src/np_sphere.cpp
+ 1
− 4
View file @ 7dacec5c
...@@ -9,7 +9,6 @@ ...@@ -9,7 +9,6 @@
using namespace std; using namespace std;
extern void cluster();
extern void sphere(); extern void sphere();
/*! \brief Main program entry point. /*! \brief Main program entry point.
...@@ -21,8 +20,6 @@ extern void sphere(); ...@@ -21,8 +20,6 @@ extern void sphere();
* the configuration and runs the main program. * the configuration and runs the main program.
*/ */
int main(int argc, char **argv) { int main(int argc, char **argv) {
bool is_sphere = true; sphere();
if (is_sphere) sphere();
else cluster();
return 0; return 0;
} }
src/sphere/Makefile
+ 2
− 5
View file @ 7dacec5c
...@@ -14,8 +14,8 @@ edfb: edfb.o ...@@ -14,8 +14,8 @@ edfb: edfb.o
sph: sph.o sph: sph.o
$(FC) $(FCFLAGS) -o $(BUILDDIR)/sph $(BUILDDIR)/sph.o $(LFLAGS) $(FC) $(FCFLAGS) -o $(BUILDDIR)/sph $(BUILDDIR)/sph.o $(LFLAGS)
np_sphere: $(BUILDDIR)/np_sphere.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sphere.o $(BUILDDIR)/cluster.o np_sphere: $(BUILDDIR)/np_sphere.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sphere.o
$(CXX) $(CXXFLAGS) $(CXXLFLAGS) -o $(BUILDDIR)/np_sphere $(BUILDDIR)/np_sphere.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sphere.o $(BUILDDIR)/cluster.o $(CXX) $(CXXFLAGS) $(CXXLFLAGS) -o $(BUILDDIR)/np_sphere $(BUILDDIR)/np_sphere.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sphere.o
$(BUILDDIR)/np_sphere.o: $(BUILDDIR)/np_sphere.o:
$(CXX) $(CXXFLAGS) -c ../np_sphere.cpp -o $(BUILDDIR)/np_sphere.o $(CXX) $(CXXFLAGS) -c ../np_sphere.cpp -o $(BUILDDIR)/np_sphere.o
...@@ -29,9 +29,6 @@ $(BUILDDIR)/Configuration.o: ...@@ -29,9 +29,6 @@ $(BUILDDIR)/Configuration.o:
$(BUILDDIR)/Parsers.o: $(BUILDDIR)/Parsers.o:
$(CXX) $(CXXFLAGS) -c ../libnptm/Parsers.cpp -o $(BUILDDIR)/Parsers.o $(CXX) $(CXXFLAGS) -c ../libnptm/Parsers.cpp -o $(BUILDDIR)/Parsers.o
$(BUILDDIR)/cluster.o:
$(CXX) $(CXXFLAGS) -c ../cluster/cluster.cpp -o $(BUILDDIR)/cluster.o
$(BUILDDIR)/sphere.o: $(BUILDDIR)/sphere.o:
$(CXX) $(CXXFLAGS) -c sphere.cpp -o $(BUILDDIR)/sphere.o $(CXX) $(CXXFLAGS) -c sphere.cpp -o $(BUILDDIR)/sphere.o
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment