Skip to content

Commit 6efc024

Browse files
author
Martin D. Weinberg
committed
Addded a test routine for Abel deprojection; updated EmpCylSL for derivative form, which is more accurate near the center
1 parent c7bf5ee commit 6efc024

9 files changed

Lines changed: 539 additions & 9 deletions

File tree

‎exputil/EmpCylSL.cc‎

Lines changed: 25 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -487,28 +487,45 @@ void EmpCylSL::create_deprojection(double H, double Rf, int NUMR, int NINT,
487487
// Now, compute Abel inverion integral
488488
//
489489
for (int i=0; i<NUMR; i++) {
490+
491+
double surfS = abel_type==AbelType::Subtraction ? surf.eval(rl[i]) : 0.0;
490492
double r = rr[i];
491493

492494
// Interval by Legendre
493495
//
494496
rhoI[i] = 0.0;
497+
495498
for (int n=0; n<NINT; n++) {
499+
496500
double x = lq.knot(n);
497501
double x12 = 1.0 - x*x;
498-
double z = x/sqrt(x12);
499-
double R = sqrt(z*z + r*r);
502+
double R = r/sqrt(x12);
500503
double lR = log(R);
501504

502-
rhoI[i] += lq.weight(n)*2.0*pow(x12, -1.5)*surf.eval(lR);
505+
switch (abel_type) {
506+
case AbelType::Derivative:
507+
rhoI[i] += lq.weight(n)/x12 * surf.deriv(lR)/R;
508+
break;
509+
case AbelType::Subtraction:
510+
rhoI[i] += lq.weight(n)/(x*x*sqrt(x12)*r) * ( surf.eval(lR) - surfS );
511+
break;
512+
case AbelType::IBP:
513+
rhoI[i] += lq.weight(n)/x12 * surf.eval(lR);
514+
break;
515+
}
503516
}
504517
}
505518

506-
std::vector<double> rho(NUMR), mass(NUMR);
507-
508-
Linear1d intgr(rl, rhoI);
519+
Linear1d intgr;
520+
if (abel_type == AbelType::IBP) intgr = Linear1d(rl, rhoI);
509521

510-
for (int i=0; i<NUMR; i++)
511-
rho[i] = -intgr.deriv(rl[i])/(2.0*M_PI*rr[i]*rr[i]);
522+
std::vector<double> rho(NUMR), mass(NUMR);
523+
for (int i=0; i<NUMR; i++) {
524+
if (abel_type == AbelType::IBP)
525+
rho[i] = -intgr.deriv(rl[i])/rr[i]/M_PI;
526+
else
527+
rho[i] = -rhoI[i]/M_PI;
528+
}
512529

513530
mass[0] = 0.0;
514531
for (int i=1; i<NUMR; i++) {

‎include/EmpCylSL.H‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -351,6 +351,12 @@ protected:
351351
//! with the new Eigen3 API
352352
bool allow_old_cache = false;
353353

354+
//! Abel type for deprojection
355+
enum class AbelType { Derivative, Subtraction, IBP };
356+
357+
//! Deprojection method (default: derivative)
358+
AbelType abel_type = AbelType::Derivative;
359+
354360
public:
355361

356362
/*! Enum listing the possible selection algorithms for coefficient

‎utils/Test/CMakeLists.txt‎

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11

2-
set(bin_PROGRAMS testBarrier expyaml orthoTest)
2+
set(bin_PROGRAMS testBarrier expyaml orthoTest testDeproj
3+
testDeprojPlum testEmpDeproj)
34

45
set(common_LINKLIB OpenMP::OpenMP_CXX MPI::MPI_CXX expui exputil
56
yaml-cpp ${VTK_LIBRARIES})
@@ -32,6 +33,11 @@ endif()
3233
add_executable(testBarrier test_barrier.cc)
3334
add_executable(expyaml test_config.cc)
3435
add_executable(orthoTest orthoTest.cc Biorth2Ortho.cc)
36+
add_executable(testDeproj testDeproject.cc CubicSpline.cc Deprojector.cc)
37+
add_executable(testDeprojPlum testDeprojPlummer.cc CubicSpline.cc
38+
Deprojector.cc)
39+
add_executable(testEmpDeproj testEmpDeproj.cc CubicSpline.cc
40+
Deprojector.cc EmpDeproj.cc)
3541

3642
foreach(program ${bin_PROGRAMS})
3743
target_link_libraries(${program} ${common_LINKLIB})

‎utils/Test/CubicSpline.H‎

Lines changed: 31 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,31 @@
1+
#ifndef CUBIC_SPLINE_H
2+
#define CUBIC_SPLINE_H
3+
4+
#include <vector>
5+
6+
namespace Deproject
7+
{
8+
9+
class CubicSpline {
10+
public:
11+
CubicSpline() = default;
12+
CubicSpline(const std::vector<double>& x_in, const std::vector<double>& y_in);
13+
14+
// set data (x must be strictly increasing)
15+
void set_data(const std::vector<double>& x_in, const std::vector<double>& y_in);
16+
17+
// evaluate spline and its derivative (xx should lie within [xmin(), xmax()], but endpoints are clamped)
18+
double eval(double xx) const;
19+
double deriv(double xx) const;
20+
21+
double xmin() const;
22+
double xmax() const;
23+
24+
private:
25+
std::vector<double> x_, y_, y2_;
26+
int locate(double xx) const;
27+
};
28+
29+
}
30+
31+
#endif // CUBIC_SPLINE_H

‎utils/Test/CubicSpline.cc‎

Lines changed: 73 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,73 @@
1+
#include <stdexcept>
2+
3+
#include "CubicSpline.H"
4+
5+
namespace Deproject
6+
{
7+
8+
CubicSpline::CubicSpline(const std::vector<double>& x_in, const std::vector<double>& y_in) {
9+
set_data(x_in, y_in);
10+
}
11+
12+
void CubicSpline::set_data(const std::vector<double>& x_in, const std::vector<double>& y_in) {
13+
x_ = x_in;
14+
y_ = y_in;
15+
int n = (int)x_.size();
16+
if (n < 2 || (int)y_.size() != n) throw std::runtime_error("CubicSpline: need at least two points and equal-sized x,y.");
17+
18+
y2_.assign(n, 0.0);
19+
std::vector<double> u(n - 1, 0.0);
20+
21+
// natural spline boundary conditions (second derivatives at endpoints = 0)
22+
for (int i = 1; i < n - 1; ++i) {
23+
double sig = (x_[i] - x_[i-1]) / (x_[i+1] - x_[i-1]);
24+
double p = sig * y2_[i-1] + 2.0;
25+
y2_[i] = (sig - 1.0) / p;
26+
double dY1 = (y_[i+1] - y_[i]) / (x_[i+1] - x_[i]);
27+
double dY0 = (y_[i] - y_[i-1]) / (x_[i] - x_[i-1]);
28+
u[i] = (6.0 * (dY1 - dY0) / (x_[i+1] - x_[i-1]) - sig * u[i-1]) / p;
29+
}
30+
31+
for (int k = n - 2; k >= 0; --k) y2_[k] = y2_[k] * y2_[k+1] + u[k];
32+
}
33+
34+
int CubicSpline::locate(double xx) const {
35+
int n = (int)x_.size();
36+
if (xx <= x_.front()) return 0;
37+
if (xx >= x_.back()) return n - 2;
38+
int lo = 0, hi = n - 1;
39+
while (hi - lo > 1) {
40+
int mid = (lo + hi) >> 1;
41+
if (x_[mid] > xx) hi = mid; else lo = mid;
42+
}
43+
return lo;
44+
}
45+
46+
double CubicSpline::eval(double xx) const {
47+
int klo = locate(xx);
48+
int khi = klo + 1;
49+
double h = x_[khi] - x_[klo];
50+
if (h <= 0.0) throw std::runtime_error("CubicSpline::eval: non-increasing x.");
51+
double A = (x_[khi] - xx) / h;
52+
double B = (xx - x_[klo]) / h;
53+
double val = A * y_[klo] + B * y_[khi]
54+
+ ((A*A*A - A) * y2_[klo] + (B*B*B - B) * y2_[khi]) * (h*h) / 6.0;
55+
return val;
56+
}
57+
58+
double CubicSpline::deriv(double xx) const {
59+
int klo = locate(xx);
60+
int khi = klo + 1;
61+
double h = x_[khi] - x_[klo];
62+
if (h <= 0.0) throw std::runtime_error("CubicSpline::deriv: non-increasing x.");
63+
double A = (x_[khi] - xx) / h;
64+
double B = (xx - x_[klo]) / h;
65+
double dy = (y_[khi] - y_[klo]) / h
66+
+ ( - (3.0*A*A - 1.0) * y2_[klo] + (3.0*B*B - 1.0) * y2_[khi] ) * (h / 6.0);
67+
return dy;
68+
}
69+
70+
double CubicSpline::xmin() const { return x_.front(); }
71+
double CubicSpline::xmax() const { return x_.back(); }
72+
73+
}

‎utils/Test/Deprojector.H‎

Lines changed: 58 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,58 @@
1+
#ifndef DEPROJECTOR_H
2+
#define DEPROJECTOR_H
3+
4+
#include <functional>
5+
#include <vector>
6+
7+
#include "CubicSpline.H"
8+
9+
namespace Deproject
10+
{
11+
12+
class Deprojector {
13+
public:
14+
// Construct from analytic Sigma functor. If dSigmaFunc is empty, numerical derivative is used.
15+
// R_data_min and R_data_max define where SigmaFunc is trusted for tail matching.
16+
Deprojector(std::function<double(double)> SigmaFunc,
17+
std::function<double(double)> dSigmaFunc,
18+
double R_data_min,
19+
double R_data_max,
20+
double R_max_extend = -1.0,
21+
double tail_power = -3.0,
22+
int Ngrid = 4000);
23+
24+
// Construct from sampled data arrays (R must be positive and will be sorted)
25+
Deprojector(const std::vector<double>& R_in,
26+
const std::vector<double>& Sigma_in,
27+
double R_max_extend = -1.0,
28+
double tail_power = -3.0,
29+
int Ngrid = 4000);
30+
31+
// Evaluate rho at a single radius
32+
double rho_at(double r) const;
33+
34+
// Evaluate rho for a vector of r
35+
std::vector<double> rho(const std::vector<double>& r_eval) const;
36+
37+
private:
38+
void build_grid(int Ngrid);
39+
40+
std::function<double(double)> sigma_func_;
41+
std::function<double(double)> dsigma_func_; // may be empty
42+
CubicSpline spline_; // used when constructed from sampled data
43+
44+
double Rdata_min_, Rdata_max_;
45+
double Rmin_, Rmax_;
46+
double tail_power_;
47+
48+
std::vector<double> dx_; // grid spacings
49+
50+
std::vector<double> fineR_;
51+
std::vector<double> Sigma_f_;
52+
std::vector<double> dSigma_f_;
53+
};
54+
55+
}
56+
57+
#endif // DEPROJECTOR_H
58+

0 commit comments

Comments
 (0)