From 92c1d8b2d18c3bec71296b24f67c8fc42fdc86f3 Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Wed, 8 Nov 2023 14:57:18 +0100
Subject: [PATCH] Implement dynamic memory management for spherical harmonics
---
src/include/sph_subs.h | 10 +++++++---
src/sphere/sphere.cpp | 31 +++++++++++++++----------------
2 files changed, 22 insertions(+), 19 deletions(-)
diff --git a/src/include/sph_subs.h b/src/include/sph_subs.h
index 3c562fb1..8d06ba94 100644
--- a/src/include/sph_subs.h
+++ b/src/include/sph_subs.h
@@ -954,7 +954,9 @@ void sphar(
double cosrth, double sinrth, double cosrph, double sinrph,
int ll, std::complex<double> *ylm
) {
- double sinrmp[40], cosrmp[40], plegn[861];
+ int rmp_size = ll;
+ int plegn_size = (ll + 1) * ll / 2 + ll + 1;
+ double sinrmp[rmp_size], cosrmp[rmp_size], plegn[plegn_size];
double four_pi = 8.0 * acos(0.0);
double pi4irs = 1.0 / sqrt(four_pi);
double x = cosrth;
@@ -1240,7 +1242,8 @@ void wmamp(
int lm, int idot, int nsph, double *arg, double *u, double *up,
double *un, C1 *c1
) {
- std::complex<double> *ylm = new std::complex<double>[1682];
+ int ylm_size = (lm + 1) * (lm + 1) + 1;
+ std::complex<double> *ylm = new std::complex<double>[ylm_size];
int nlmp = lm * (lm + 2) + 2;
ylm[nlmp - 1] = std::complex<double>(0.0, 0.0);
if (idot != 0) {
@@ -1295,7 +1298,8 @@ void wmasp(
double *ups, double *uns, int isq, int ibf, int inpol, int lm, int idot,
int nsph, double *argi, double *args, C1 *c1
) {
- std::complex<double> *ylm = new std::complex<double>[1682];
+ int ylm_size = (lm + 1) * (lm + 1) + 1;
+ std::complex<double> *ylm = new std::complex<double>[ylm_size];
int nlmp = lm * (lm + 2) + 2;
ylm[nlmp - 1] = std::complex<double>(0.0, 0.0);
if (idot != 0) {
diff --git a/src/sphere/sphere.cpp b/src/sphere/sphere.cpp
index 5d92b192..fd677722 100644
--- a/src/sphere/sphere.cpp
+++ b/src/sphere/sphere.cpp
@@ -2,8 +2,12 @@
#include <fstream>
#include <string>
#include <complex>
+#ifndef INCLUDE_CONFIGURATION_H_
#include "../include/Configuration.h"
+#endif
+#ifndef INCLUDE_SPH_SUBS_H_
#include "../include/sph_subs.h"
+#endif
using namespace std;
@@ -15,13 +19,13 @@ void sphere() {
double *argi, *args, *gaps;
double th, ph;
printf("INFO: making legacy configuration ...\n");
- ScattererConfiguration *conf = ScattererConfiguration::from_dedfb("../../test_data/sphere/DEDFB");
- conf->write_formatted("c_OEDFB");
- conf->write_binary("c_TEDF");
+ ScattererConfiguration *conf = ScattererConfiguration::from_dedfb("DEDFB_sph");
+ conf->write_formatted("c_OEDFB_sph");
+ conf->write_binary("c_TEDF_sph");
delete conf;
printf("INFO: reading binary configuration ...\n");
- ScattererConfiguration *sconf = ScattererConfiguration::from_binary("c_TEDF");
- GeometryConfiguration *gconf = GeometryConfiguration::from_legacy("../../test_data/sphere/DSPH");
+ ScattererConfiguration *sconf = ScattererConfiguration::from_binary("c_TEDF_sph");
+ GeometryConfiguration *gconf = GeometryConfiguration::from_legacy("DSPH");
if (sconf->number_of_spheres == gconf->number_of_spheres) {
int isq, ibf;
double cost, sint, cosp, sinp;
@@ -49,18 +53,10 @@ void sphere() {
complex<double> *vint = new complex<double>[16];
int jw;
int nsph = gconf->number_of_spheres;
- C1 *c1 = new C1(nsph, gconf->l_max);
+ C1 *c1 = new C1(nsph, gconf->l_max, sconf->nshl_vec, sconf->iog_vec);
for (int i = 0; i < nsph; i++) {
- c1->iog[i] = sconf->iog_vec[i];
- c1->nshl[i] = sconf->nshl_vec[i];
c1->ros[i] = sconf->radii_of_spheres[i];
}
- for (int i = 0; i < gconf->l_max; i++) {
- c1->rmi[i][0] = complex<double>(0.0, 0.0);
- c1->rmi[i][1] = complex<double>(0.0, 0.0);
- c1->rei[i][0] = complex<double>(0.0, 0.0);
- c1->rei[i][1] = complex<double>(0.0, 0.0);
- }
C2 *c2 = new C2(nsph, 5, gconf->npnt, gconf->npntts);
argi = new double[1];
args = new double[1];
@@ -172,7 +168,7 @@ void sphere() {
double exri = sqrt(sconf->exdc);
fprintf(output, " REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n", exri);
fstream tppoan;
- tppoan.open("c_TPPOAN", ios::binary|ios::out);
+ tppoan.open("c_TPPOAN_sph", ios::binary|ios::out);
if (tppoan.is_open()) {
int imode = 10;
tppoan.write(reinterpret_cast<char *>(&imode), sizeof(int));
@@ -250,7 +246,7 @@ void sphere() {
// This is the condition that writes the transition matrix to output.
int is = 1111;
fstream ttms;
- ttms.open("c_TTMS", ios::binary | ios::out);
+ ttms.open("c_TTMS_sph", ios::binary | ios::out);
if (ttms.is_open()) {
ttms.write(reinterpret_cast<char *>(&is), sizeof(int));
ttms.write(reinterpret_cast<char *>(&(gconf->l_max)), sizeof(int));
@@ -550,9 +546,12 @@ void sphere() {
}
delete[] cmul;
delete[] cmullr;
+ printf("Done.\n");
} else { // NSPH mismatch between geometry and scatterer configurations.
throw UnrecognizedConfigurationException(
"Inconsistent geometry and scatterer configurations."
);
}
+ delete sconf;
+ delete gconf;
}
--
GitLab