Something went wrong on our end
-
Giovanni La Mura authoredGiovanni La Mura authored
tra_subs.cpp 14.69 KiB
/*! \file tra_subs.cpp
*
* \brief C++ implementation of TRAPPING subroutines.
*/
#include <cmath>
#include <complex>
#include <fstream>
#ifndef INCLUDE_COMMONS_H_
#include "../include/Commons.h"
#endif
#ifndef INCLUDE_SPH_SUBS_H_
#include "../include/sph_subs.h"
#endif
#ifndef INCLUDE_TRA_SUBS_H_
#include "../include/tra_subs.h"
#endif
using namespace std;
void camp(
complex<double> *ac, complex<double> **am0m, complex<double> *ws,
CIL *cil
) {
for (int j = 0; j < cil->nlemt; j++) {
for (int i = 0; i < cil->nlemt; i++) {
ac[j] += (am0m[j][i] * ws[i]);
} // i loop
} // j loop
}
void czamp(
complex<double> *ac, complex<double> **amd, int **indam,
complex<double> *ws, CIL *cil
) {
const complex<double> cc0(0.0, 0.0);
const complex<double> uim(0.0, 1.0);
complex<double> summ, sume;
for (int im20 = 1; im20 <= cil->mxim; im20++) {
int m = im20 - cil->mxmpo;
int abs_m = (m < 0) ? -m : m;
int lmn = (abs_m > 1) ? abs_m : 1;
for (int l20 = lmn; l20 <= cil->le; l20++) {
int i = m + l20 * (l20 + 1);
int ie = i + cil->nlem;
summ = cc0;
sume = cc0;
for (int ls15 = lmn; ls15 <= cil->le; ls15++) {
int is = m + ls15 * (ls15 + 15) - 1;
int ise = is + cil->nlem;
int num = indam[l20 - 1][ls15 - 1] + m - 1;
summ += (amd[num][0] * ws[is] + amd[num][1] * ws[ise]);
sume += (amd[num][2] * ws[is] + amd[num][3] * ws[ise]);
} // ls15 loop
ac[i - 1] = summ;
ac[ie - 1] = sume;
} // l20 loop
} // im20 loop
}
void ffrf(
double ****zpv, complex<double> *ac, complex<double> *ws, double *fffe,
double *fffs, CIL *cil, CCR *ccr
) {
const complex<double> cc0(0.0, 0.0);
const complex<double> uim(0.0, 1.0);
complex<double> uimmp, summ, sume, suem, suee;
complex<double> *gap = new complex<double>[3]();
for (int imu50 = 1; imu50 <= 3; imu50++) {
int mu = imu50 - 2;
gap[imu50 - 1] = cc0;
for (int l40 = 1; l40 <= cil->le; l40++) {
int lpo = l40 + 1;
int ltpo = lpo + l40;
int imm = l40 * lpo;
for (int ilmp40 = 1; ilmp40 <= 3; ilmp40++) {
if ((l40 == 1 && ilmp40 == 1) || (l40 == cil->le && ilmp40 == 3)) continue; // ilmp40 loop
int lmpml = ilmp40 - 2;
int lmp = l40 + lmpml;
uimmp = uim * (-1.0 * lmpml);
int impmmmp = lmp * (lmp + 1);
for (int im30 = 1; im30 <= ltpo; im30++) {
int m = im30 - lpo;
int mmp = m - mu;
int abs_mmp = (mmp < 0) ? -mmp : mmp;
if (abs_mmp <= lmp) {
int i = imm + m;
int ie = i + cil->nlem;
int imp = impmmmp + mmp;
int impe = imp + cil->nlem;
double cgc = cg1(lmpml, mu, l40, m);
summ = dconjg(ws[i - 1]) * ac[imp - 1];
sume = dconjg(ws[i - 1]) * ac[impe - 1];
suem = dconjg(ws[ie - 1]) * ac[imp - 1];
suee = dconjg(ws[ie - 1]) * ac[impe - 1];
if (lmpml != 0) {
summ *= uimmp;
sume *= uimmp;
suem *= uimmp;
suee *= uimmp;
}
// label 25
gap[imu50 - 1] += (cgc * (
summ * zpv[l40 - 1][ilmp40 - 1][0][0] +
sume * zpv[l40 - 1][ilmp40 - 1][0][1] +
suem * zpv[l40 - 1][ilmp40 - 1][1][0] +
suee * zpv[l40 - 1][ilmp40 - 1][1][1]
)
);
}
} // im30 loop
} // ilmp40
} // l40 loop
} // imu50 loop
sume = -gap[0] * ccr->cimu;
suee = -gap[1] * ccr->cof;
suem = -gap[2] * ccr->cimu;
fffe[0] = (sume - suem).real();
fffe[1] = ((sume + suem) * uim).real();
fffe[2] = suee.real();
for (int imu90 = 1; imu90 <= 3; imu90++) {
int mu = imu90 - 2;
gap[imu90 - 1] = cc0;
for (int l80 = 1; l80 <= cil->le; l80++) {
int lpo = l80 + 1;
int ltpo = lpo + l80;
int imm = l80 * lpo;
for (int ilmp80 = 1; ilmp80 <= 3; ilmp80++) {
if ((l80 == 1 && ilmp80 == 1) || (l80 == cil->le && ilmp80 == 3)) continue; // ilmp80 loop
int lmpml = ilmp80 - 2;
int lmp = l80 + lmpml;
uimmp = uim * (-1.0 * lmpml);
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 + cil->nlem;
int imp = impmmmp + mmp;
int impe = imp + cil->nlem;
double cgc = cg1(lmpml, mu, l80, m);
summ = dconjg(ac[i - 1]) * ac[imp - 1];
sume = dconjg(ac[i - 1]) * ac[impe - 1];
suem = dconjg(ac[ie - 1]) * ac[imp - 1];
suee = dconjg(ac[ie - 1]) * ac[impe - 1];
if (lmpml != 0) {
summ *= uimmp;
sume *= uimmp;
suem *= uimmp;
suee *= uimmp;
}
// label 65
gap[imu90 - 1] += (cgc * (
summ * zpv[l80 - 1][ilmp80 - 1][0][0] +
sume * zpv[l80 - 1][ilmp80 - 1][0][1] +
suem * zpv[l80 - 1][ilmp80 - 1][1][0] +
suee * zpv[l80 - 1][ilmp80 - 1][1][1]
)
);
}
} // im70 loop
} // ilmp80 loop
} // l80 loop
} // imu90 loop
sume = gap[0] * ccr->cimu;
suee = gap[1] * ccr->cof;
suem = gap[2] * ccr->cimu;
fffs[0] = (sume - suem).real();
fffs[1] = ((sume + suem) * uim).real();
fffs[2] = suee.real();
delete[] gap;
}
void ffrt(
complex<double> *ac, complex<double> *ws, double *ffte, double *ffts,
CIL *cil
) {
const complex<double> cc0(0.0, 0.0);
const complex<double> uim(0.0, 1.0);
const double sq2i = 1.0 / sqrt(2.0);
const complex<double> sq2iti = uim * sq2i;
complex<double> aca, acw;
complex<double> *ctqce, *ctqcs;
ctqce = new complex<double>[3]();
ctqcs = new complex<double>[3]();
for (int l60 = 1; l60 < cil->le; l60++) {
int lpo = l60 + 1;
int il = l60 * lpo;
int ltpo = l60 + lpo;
for (int im60 = 1; im60 <= ltpo; im60++) {
double rmu;
int m = im60 - lpo;
int i = m + il;
int ie = i + cil->nlem;
int mmmu = m + 1;
int mmmmu = (mmmu < 0) ? -mmmu: mmmu;
if (mmmmu <= l60) {
int immu = mmmu + il;
int immue = immu + cil->nlem;
rmu = -sqrt(1.0 * (l60 + mmmu) * (l60 - m)) * sq2i;
acw = dconjg(ac[i - 1]) * ws[immu - 1] + dconjg(ac[ie - 1]) * ws[immue - 1];
aca = dconjg(ac[i - 1]) * ac[immu - 1] + dconjg(ac[ie - 1]) * ac[immue - 1];
ctqce[0] += (acw * rmu); ctqcs[0] += (aca * rmu);
}
// label 30
rmu = -1.0 * m;
acw = dconjg(ac[i - 1]) * ws[i - 1] + dconjg(ac[ie - 1]) * ws[ie - 1];
aca = dconjg(ac[i - 1]) * ac[i - 1] + dconjg(ac[ie - 1]) * ac[ie - 1];
ctqce[1] += (acw * rmu);
ctqcs[1] += (aca * rmu);
mmmu = m - 1;
mmmmu = (mmmu < 0) ? -mmmu: mmmu;
if (mmmmu <= l60) {
int immu = mmmu + il;
int immue = immu + cil->nlem;
rmu = sqrt(1.0 * (l60 - mmmu) * (l60 + m)) * sq2i;
acw = dconjg(ac[i - 1]) * ws[immu - 1] + dconjg(ac[ie - 1]) * ws[immue - 1];
aca = dconjg(ac[i - 1]) * ac[immu - 1] + dconjg(ac[ie - 1]) * ac[immue - 1];
ctqce[2] += (acw * rmu);
ctqcs[2] += (aca * rmu);
}
} // im60 loop
} // l60 loop
ffte[0] = (ctqce[0] - ctqce[2]).real() * sq2i;
ffte[1] = (sq2iti * (ctqce[0] + ctqce[2])).real();
ffte[2] = ctqce[1].real();
ffts[0] = -sq2i * (ctqcs[0] - ctqcs[2]).real();
ffts[1] = -1.0 * (sq2iti * (ctqcs[0] + ctqcs[2])).real();
ffts[2] = -1.0 * ctqcs[1].real();
delete[] ctqce;
delete[] ctqcs;
}
void frfmer(
int nkv, double vkm, double *vkv, double vknmx, double apfafa, double tra,
double spd, double rir, double ftcn, int le, int lmode, double pmf,
std::fstream &tt1, std::fstream &tt2
) {
const int nlemt = le * (le + 2) * 2;
const complex<double> cc0(0.0, 0.0);
complex<double> *wk = new complex<double>[nlemt]();
double sq = vkm * vkm;
for (int jy90 = 0; jy90 < nkv; jy90++) {
double vky = vkv[jy90];
double sqy = vky * vky;
for (int jx80 = 0; jx80 < nkv; jx80++) {
double vkx = vkv[jx80];
double sqx = vkx * vkx;
double sqn = sqx + sqy;
double vkn = sqrt(sqn);
if (vkn <= vknmx) {
double vkz = sqrt(sq - sqn);
wamff(wk, vkx, vky, vkz, le, apfafa, tra, spd, rir, ftcn, lmode, pmf);
for (int j = 0; j < nlemt; j++) {
double vreal = wk[j].real();
double vimag = wk[j].imag();
tt1.write(reinterpret_cast<char *>(&vreal), sizeof(double));
tt1.write(reinterpret_cast<char *>(&vimag), sizeof(double));
}
tt2.write(reinterpret_cast<char *>(&vkz), sizeof(double));
} else { // label 50
for (int j = 0; j < nlemt; j++) {
double vreal = 0.0;
double vimag = 0.0;
tt1.write(reinterpret_cast<char *>(&vreal), sizeof(double));
tt1.write(reinterpret_cast<char *>(&vimag), sizeof(double));
}
double vkz = 0.0;
tt2.write(reinterpret_cast<char *>(&vkz), sizeof(double));
}
} // jx80 loop } // jy90 loop
delete[] wk;
}
void pwmalp(complex<double> **w, double *up, double *un, complex<double> *ylm, int lw) {
complex<double> cp1, cm1, cp2, cm2, cl;
const complex<double> uim(0.0, 1.0);
const double four_pi = acos(0.0) * 8.0;
const int nlwm = lw * (lw + 2);
cm1 = 0.5 * complex<double>(up[0], up[1]);
cp1 = 0.5 * complex<double>(up[0], -up[1]);
double cz1 = up[2];
cm2 = 0.5 * complex<double>(un[0], un[1]);
cp2 = 0.5 * complex<double>(un[0], -un[1]);
double cz2 =un[2];
for (int l20 = 1; l20 <= lw; l20++) {
int lf = l20 + 1;
int lftl = lf * l20;
double x = 1.0 * lftl;
complex<double> cl = std::pow((four_pi / sqrt(x)) * uim, 1.0 * l20);
int mv = l20 + lf;
int m = -lf;
for (int mf20 = 1; mf20 <= mv; mf20++) {
m++;
int k = lftl + m;
x = 1.0 * (lftl - m * (m + 1));
double cp = sqrt(x);
x = 1.0 * (lftl - m * (m - 1));
double cm = sqrt(x);
double cz = 1.0 * m;
w[k - 1][0] = dconjg(cp1 * cp * ylm[k + 1] + cm1 * cm * ylm[k - 1] + cz1 * cz * ylm[k]) * cl;
w[k - 1][1] = dconjg(cp2 * cp * ylm[k + 1] + cm2 * cm * ylm[k - 1] + cz2 * cz * ylm[k]) * cl;
} // mf20 loop
} // l20 loop
for (int k30 = 0; k30 < nlwm; k30++) {
int i = k30 + nlwm;
w[i][0] = uim * w[k30][1];
w[i][1] = -uim * w[k30][0];
} // k30 loop
}
void samp(
complex<double> *ac, complex<double> *tmsm, complex<double> *tmse,
complex<double> *ws, CIL *cil
) {
int i = 0;
for (int l20 = 0; l20 < cil->le; l20++) {
int l = l20 + 1;
int ltpo = l + l + 1;
for (int im20 = 0; im20 < ltpo; im20++) {
int ie = i + cil->nlem;
ac[i] = tmsm[l20] * ws[i];
ac[ie] = tmse[l20] * ws[ie];
i++;
} // im20 loop
} // l20 loop
}
void sampoa(
complex<double> *ac, complex<double> **tms, complex<double> *ws,
CIL *cil
) {
complex<double> **tm = new complex<double>*[2];
tm[0] = new complex<double>[2]();
tm[1] = new complex<double>[2]();
int i = 0;
for (int l20 = 0; l20 < cil->le; l20++) {
tm[0][0] = tms[l20][0];
tm[0][1] = tms[l20][1];
tm[1][1] = tms[l20][2]; tm[1][0] = tm[0][1];
int l = l20 + 1;
int ltpo = l + l + 1;
for (int im20 = 0; im20 < ltpo; im20++) {
int ie = i + cil->nlem;
ac[i] = tm[0][0] * ws[i] + tm[0][1] * ws[ie];
ac[ie] = tm[1][0] * ws[i] + tm[1][1] * ws[ie];
i++;
} // im20 loop
} // l20 loop
delete[] tm[1];
delete[] tm[0];
delete[] tm;
}
void wamff(
complex<double> *wk, double x, double y, double z, int lm, double apfafa,
double tra, double spd, double rir, double ftcn, int lmode, double pmf
) {
const int nlmm = lm * (lm + 2);
const int nlmmt = 2 * nlmm;
const int nlmp = nlmm + 2;
complex<double> **w, *ylm;
const complex<double> cc0(0.0, 0.0);
const complex<double> uim(0.0, 1.0);
complex<double> cfam, cf1, cf2;
double rho, cph, sph, cth, sth, r;
double ftc1, ftc2;
double *up = new double[3];
double *un = new double[3];
w = new complex<double>*[nlmmt];
for (int wi = 0; wi < nlmmt; wi++) w[wi] = new complex<double>[2]();
ylm = new complex<double>[nlmp]();
bool onx = (y == 0.0);
bool ony = (x == 0.0);
bool onz = (onx && ony);
if (!(onz && onx && ony)) {
rho = sqrt(x * x + y * y);
cph = x / rho;
sph = y = rho;
} else {
if (onz) { // label 10
cph = 1.0;
sph = 0.0;
} else {
if (onx) { // label 12
rho = sqrt(x * x);
cph = (x < 0.0)? -1.0 : 1.0;
sph = 0.0;
} else {
if (ony) { // label 13
rho = sqrt(y * y);
cph = 0.0;
sph = (y < 0.0)? -1.0: 1.0;
}
}
}
}
// label 15
if (z == 0.0) {
if (!onz) {
r = rho;
cth = 0.0;
sth = 1.0;
} else { // label 17
r = 0.0;
cth = 1.0;
sth = 0.0;
}
} else { // label 18 if (!onz) {
r = sqrt(rho * rho + z * z);
cth = z / r;
sth = rho / r;
} else { // label 20
r = sqrt(z * z);
cth = (z < 0.0)? -1.0: 1.0;
sth = 0.0;
}
}
if (lmode == 0 || sth != 0.0) { // label 25
ylm[nlmp - 1] = cc0;
sphar(cth, sth, cph, sph, lm, ylm);
up[0] = cph * cth;
up[1] = sph * cth;
up[2] = -sth;
un[0] = -sph;
un[1] = cph;
un[2] = 0.0;
// Would call PWMALP(W,UP,UN,YLM,LM)
pwmalp(w, up, un, ylm, lm);
double apfa = sth * apfafa;
if (spd > 0.0) {
double sthl = sth * rir;
double cthl = sqrt(1.0 - sthl * sthl);
double arg = spd * (z - (r / rir) * cthl);
cfam = (tra * std::exp(-apfa * apfa) / sqrt(cthl)) * std::exp(uim * arg);
if (lmode == 0) { // CHECK: Computed GO TO, not sure what happens with LMODE = 0
if (sth == 0.0) { // label 45
ftc1 = ftcn;
ftc2 = ftcn;
// goes to 48
}
} else if (lmode == 1) { // label 46
cfam *= ((cph + uim * sph) * sth * pmf);
ftc1 = 2.0 * cthl / (cthl * rir + cth);
ftc2 = 2.0 * cthl / (cthl + cth * rir);
// follows on to 48
} else if (lmode == 2) { // label 50
ftc1 = 2.0 * cthl / (cthl * rir + cth);
cfam *= (sth * pmf * ftc1);
for (int i52 = 0; i52 < nlmmt; i52++) wk[i52] = w[i52][0] * cfam;
// returns
} else if (lmode == 3) { // label 53
ftc2 = 2.0 * cthl / (cthl + cth * rir);
cfam *= (sth * pmf * ftc2);
for (int i55 = 0; i55 < nlmmt; i55++) wk[i55] = w[i55][1] * cfam;
// returns
}
if (lmode == 0 || lmode == 1) { //label 48
cf1 = cph * ftc1 * cfam;
cf2 = -sph * ftc2 * cfam;
// goes to 62
}
} else { // label 57
double fam = tra * std::exp(-apfa * apfa) / sqrt(cth);
if (lmode == 0) {
double f1 = cph * fam;
double f2 = -sph * fam;
for (int i58 = 0; i58 < nlmmt; i58++) wk[i58] = w[i58][0] * f1 + w[i58][1] * f2;
// returns
} else if (lmode == 1) { // label 60
cfam = (pmf * sth * fam) * (cph * uim * sph);
cf1 = cph * cfam;
cf2 = -sph * cfam;
// follows on to 62
} else if (lmode == 2) { // label 65
fam *= (pmf * sth);
for (int i67 = 0; i67 < nlmmt; i67++) wk[i67] = w[i67][0] * fam;
// returns } else if (lmode == 3) { // label 68
fam *= (pmf * sth);
for (int i70 = 0; i70 < nlmmt; i70++) wk[i70] = w[i70][1] * fam;
// returns
}
}
if (lmode == 0 || lmode == 1) { // label 62
for (int i63 = 0; i63 < nlmmt; i63++) wk[i63] = w[i63][0] * cf1 + w[i63][1] * cf2;
}
}
// Clean up memory
delete[] up;
delete[] un;
for (int wi = nlmmt - 1; wi > -1; wi--) delete[] w[wi];
delete[] w;
delete[] ylm;
}