Skip to content

Commit 368db0a

Browse files
author
Martin D. Weinberg
committed
Merge branch 'ICupdate' into devel
2 parents 1642cc2 + 72fc23b commit 368db0a

11 files changed

Lines changed: 1134 additions & 257 deletions

File tree

‎exputil/EXPmath.cc‎

Lines changed: 16 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -61,12 +61,12 @@ namespace AltMath
6161
return ans;
6262
}
6363

64-
#define ACC 40.0
65-
#define BIGNO 1.0e10
66-
#define BIGNI 1.0e-10
67-
6864
double bessj(int n,double x)
6965
{
66+
const double ACC = 40.0;
67+
const double BIGNO = 1.0e10;
68+
const double BIGNI = 1.0e-10;
69+
7070
int j,jsum,m;
7171
double ax,bj,bjm,bjp,sum,tox,ans;
7272

@@ -114,10 +114,6 @@ namespace AltMath
114114
return x < 0.0 && n%2 == 1 ? -ans : ans;
115115
}
116116

117-
#undef ACC
118-
#undef BIGNO
119-
#undef BIGNI
120-
121117
double bessk0(double x)
122118
{
123119
double y,ans;
@@ -174,12 +170,12 @@ namespace AltMath
174170
return ans;
175171
}
176172

177-
#define ACC 40.0
178-
#define BIGNO 1.0e10
179-
#define BIGNI 1.0e-10
180-
181173
double bessi(int n,double x)
182174
{
175+
const double ACC = 40.0;
176+
const double BIGNO = 1.0e10;
177+
const double BIGNI = 1.0e-10;
178+
183179
int j;
184180
double bi,bim,bip,tox,ans;
185181
double bessi0(double );
@@ -207,10 +203,6 @@ namespace AltMath
207203
}
208204
}
209205

210-
#undef ACC
211-
#undef BIGNO
212-
#undef BIGNI
213-
214206
double bessi0(double x)
215207
{
216208
double ax,ans;
@@ -253,12 +245,12 @@ namespace AltMath
253245
return x < 0.0 ? -ans : ans;
254246
}
255247

256-
#define ACC 40.0
257-
#define BIGNO 1.0e10
258-
#define BIGNI 1.0e-10
259-
260248
double jn_sph(int n, double x)
261249
{
250+
const double ACC = 40.0;
251+
const double BIGNO = 1.0e10;
252+
const double BIGNI = 1.0e-10;
253+
262254
int j,m;
263255
double ax,bj,bjm,bjp,tox,ans;
264256
double jn_sph0(double x),jn_sph1(double x),jn_sphm1(double x);
@@ -306,13 +298,10 @@ namespace AltMath
306298
return x < 0.0 && n%2 == 1 ? -ans : ans;
307299
}
308300

309-
#undef ACC
310-
#undef BIGNO
311-
#undef BIGNI
312-
313-
#define TOL 1.0e-7
314301
double jn_sph0(double x)
315302
{
303+
const double TOL = 1.0e-7;
304+
316305
if (fabs(x)<TOL)
317306
return 1.0-x*x/6.0;
318307
else
@@ -321,6 +310,8 @@ namespace AltMath
321310

322311
double jn_sph1(double x)
323312
{
313+
const double TOL = 1.0e-7;
314+
324315
if (fabs(x)<TOL)
325316
return x/3.0;
326317
else
@@ -331,7 +322,6 @@ namespace AltMath
331322
{
332323
return cos(x)/x;
333324
}
334-
#undef TOL
335325

336326
double check_nu(double nu)
337327
{

‎exputil/EmpCyl2d.cc‎

Lines changed: 17 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -780,19 +780,31 @@ bool EmpCyl2d::ReadH5Cache()
780780
auto checkInt = [&file](int value, std::string name)
781781
{
782782
int v; HighFive::Attribute vv = file.getAttribute(name); vv.read(v);
783-
if (value == v) return true; return false;
783+
if (value == v) return true;
784+
if (myid==0) std::cout << "---- EmpCyl2d::ReadH5Cache: "
785+
<< name << " expected " << value << " but found "
786+
<< v << std::endl;
787+
return false;
784788
};
785789

786790
auto checkDbl = [&file](double value, std::string name)
787791
{
788792
double v; HighFive::Attribute vv = file.getAttribute(name); vv.read(v);
789-
if (fabs(value - v) < 1.0e-16) return true; return false;
793+
if (fabs(value - v) < 1.0e-16) return true;
794+
if (myid==0) std::cout << "---- EmpCyl2d::ReadH5Cache: "
795+
<< name << " expected " << value << " but found "
796+
<< v << std::endl;
797+
return false;
790798
};
791799

792800
auto checkStr = [&file](std::string value, std::string name)
793801
{
794802
std::string v; HighFive::Attribute vv = file.getAttribute(name); vv.read(v);
795-
if (value.compare(v)==0) return true; return false;
803+
if (value.compare(v)==0) return true;
804+
if (myid==0) std::cout << "---- EmpCyl2d::ReadH5Cache: "
805+
<< name << " expected " << value << " but found "
806+
<< v << std::endl;
807+
return false;
796808
};
797809

798810
//
@@ -803,8 +815,8 @@ bool EmpCyl2d::ReadH5Cache()
803815
if (not checkStr(Version, "Version")) return false;
804816
} else {
805817
if (myid==0)
806-
std::cout << "---- EmpCyl2d::ReadH5Cache: "
807-
<< "recomputing cache for HighFive API change"
818+
std::cout << "---- EmpCyl2d::ReadH5Cache: " << "<" << cache_name_2d
819+
<< "> recomputing cache for HighFive API change"
808820
<< std::endl;
809821
return false;
810822
}

‎exputil/massmodel_dist.cc‎

Lines changed: 29 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -49,14 +49,14 @@
4949
#include <massmodel.H>
5050
#include <interp.H>
5151

52-
#define OFFSET 1.0e-3
53-
#define OFFTOL 1.2
52+
const double OFFSET=1.0e-3;
53+
const double OFFTOL=1.2;
5454

5555
extern double gint_0(double a, double b, std::function<double(double)> f, int NGauss);
5656
extern double gint_2(double a, double b, std::function<double(double)> f, int NGauss);
5757

58-
#define TSTEP 1.0e-8
59-
#define NGauss 96
58+
const double TSTEP=1.0e-8;
59+
const int NGauss=96;
6060

6161
static int DIVERGE=0;
6262

@@ -109,6 +109,10 @@ void SphericalModelTable::setup_df(int NUM, double RA)
109109
rhoQy.resize(num);
110110
rhoQy2.resize(num);
111111

112+
rhoQx. setZero();
113+
rhoQy. setZero();
114+
rhoQy2.setZero();
115+
112116
for (int i=0; i<num; i++) {
113117
x = density.x[i];
114118
rhoQx[i] = pot.y[i];
@@ -138,14 +142,20 @@ void SphericalModelTable::setup_df(int NUM, double RA)
138142
dfc.Q .resize(NUM);
139143
dfc.fQ .resize(NUM);
140144
dfc.ffQ.resize(NUM);
145+
146+
dfc.Q .setZero();
147+
dfc.fQ .setZero();
148+
dfc.ffQ.setZero();
149+
141150
dfc.num = NUM;
142151
dfc.ra2 = ra2;
143152

144153
Qmax = get_pot(pot.x[pot.num-1]);
145154
Qmin = get_pot(pot.x[0]);
146155
dQ = (Qmax - Qmin)/(double)(dfc.num-1);
147156

148-
foffset = -std::numeric_limits<double>::max();
157+
// foffset = -std::numeric_limits<double>::max();
158+
foffset = -1.0e42;
149159
dfc.Q[dfc.num-1] = Qmax;
150160
dfc.ffQ[dfc.num-1] = 0.0;
151161
fac = 1.0/(sqrt(8.0)*M_PI*M_PI);
@@ -184,6 +194,13 @@ void SphericalModelTable::setup_df(int NUM, double RA)
184194
df.ffQ .resize(NUM);
185195
df.fQ2 .resize(NUM);
186196
df.ffQ2.resize(NUM);
197+
198+
df.Q .setZero();
199+
df.fQ .setZero();
200+
df.ffQ .setZero();
201+
df.fQ2 .setZero();
202+
df.ffQ2.setZero();
203+
187204
df.num = NUM;
188205
df.ra2 = ra2;
189206

@@ -194,7 +211,8 @@ void SphericalModelTable::setup_df(int NUM, double RA)
194211
df.Q[df.num-1] = Qmax;
195212
df.ffQ[df.num-1] = 0.0;
196213
fac = 1.0/(sqrt(8.0)*M_PI*M_PI);
197-
foffset = -std::numeric_limits<double>::max();
214+
// foffset = -std::numeric_limits<double>::max();
215+
foffset = -1.0e42;
198216
for (int i=df.num-2; i>=0; i--) {
199217
df.Q[i] = df.Q[i+1] - dQ;
200218
Q = df.Q[i];
@@ -233,7 +251,7 @@ void SphericalModelTable::setup_df(int NUM, double RA)
233251

234252
dist_defined = true;
235253

236-
debug_fdist();
254+
// debug_fdist();
237255
}
238256

239257

@@ -294,7 +312,7 @@ double SphericalModelTable::distf(double E, double L)
294312

295313
if (!dist_defined) bomb("distribution function not defined");
296314

297-
double d, g;
315+
double d=0.0, g=0.0;
298316

299317
if (chebyN) {
300318

@@ -341,7 +359,7 @@ double SphericalModelTable::dfde(double E, double L)
341359

342360
if (!dist_defined) bomb("distribution function not defined");
343361

344-
double d, g, h, d1;
362+
double d=0, g=0, h=0, d1=0;
345363

346364
if (chebyN) {
347365

@@ -399,7 +417,7 @@ double SphericalModelTable::dfdl(double E, double L)
399417

400418
if (!dist_defined) bomb("distribution function not defined");
401419

402-
double d, g, h, d1;
420+
double d=0, g=0, h=0, d1=0;
403421

404422
if (chebyN) {
405423

@@ -451,7 +469,7 @@ double SphericalModelTable::d2fde2(double E, double L)
451469
{
452470
if (!dist_defined) bomb("distribution function not defined");
453471

454-
double d, g, h, k, d2;
472+
double d=0, g=0, h=0, k=0, d2=0;
455473

456474
if (chebyN) {
457475

0 commit comments

Comments
 (0)