Skip to content

Commit 73ec0b2

Browse files
mlietzowmlietzow
committed
- fix: memory error if sidelength is used for a polar detector
- updated calculation of reduction factor for GOLD alignment - show error if photon package is not transferred to next cell border - use loglinear method for interpolation of refractive index
1 parent c0b4421 commit 73ec0b2

14 files changed

+144
-92
lines changed

‎README.md

+1-1
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@
55
[![bibcode](https://img.shields.io/badge/bibcode-2016A%26A...593A..87R-1c459b)](https://ui.adsabs.harvard.edu/abs/2016A&A...593A..87R)
66
[![doi](https://img.shields.io/badge/doi-10.1051%2F0004--6361%2F201424930-fab70c)](https://doi.org/10.1051/0004-6361/201424930)
77
[![License](https://img.shields.io/badge/License-GPLv3-blue)](https://www.gnu.org/licenses/gpl-3.0)
8-
[![Version](https://img.shields.io/badge/Version-4.12.01-bf0040)](https://img.shields.io/badge/Version-4.12.01-bf0040)
8+
[![Version](https://img.shields.io/badge/Version-4.12.02-bf0040)](https://img.shields.io/badge/Version-4.12.02-bf0040)
99

1010
is a 3D Monte Carlo radiative transfer code that
1111

‎src/Cylindrical.cpp

+7
Original file line numberDiff line numberDiff line change
@@ -1527,6 +1527,13 @@ bool CGridCylindrical::goToNextCellBorder(photon_package * pp)
15271527
}
15281528

15291529
pp->setPosition(p + d * path_length);
1530+
1531+
if(p == pp->getPosition())
1532+
{
1533+
cout << "\nERROR: Could not transfer photon to the next cell border! " << endl;
1534+
return false;
1535+
}
1536+
15301537
pp->setTmpPathLength(path_length);
15311538
pp->setDirectionID(dirID);
15321539
return true;

‎src/Cylindrical.h

+8-8
Original file line numberDiff line numberDiff line change
@@ -445,14 +445,14 @@ class CGridCylindrical : public CGridBasic
445445
uint & N_polar_r,
446446
uint *& N_polar_ph)
447447
{
448-
return CGridBasic::getPolarRTGridParameterWorker(max_len,
449-
pixel_width,
450-
max_subpixel_lvl,
451-
_listR,
452-
N_polar_r,
453-
N_polar_ph,
454-
N_r,
455-
listR);
448+
return CGridBasic::getPolarRTGridParameterWorker(max_len,
449+
pixel_width,
450+
max_subpixel_lvl,
451+
_listR,
452+
N_polar_r,
453+
N_polar_ph,
454+
N_r,
455+
listR);
456456
}
457457

458458
private:

‎src/Dust.cpp

+12-11
Original file line numberDiff line numberDiff line change
@@ -1176,8 +1176,8 @@ bool CDustComponent::readDustRefractiveIndexFile(parameters & param,
11761176
dcomplex refractive_index = dcomplex(refractive_index_real.getValue(wavelength_list[w], LOGLINEAR),
11771177
refractive_index_imag.getValue(wavelength_list[w], LOGLINEAR));
11781178
#else
1179-
dcomplex refractive_index = dcomplex(refractive_index_real.getValue(wavelength_list[w], SPLINE),
1180-
refractive_index_imag.getValue(wavelength_list[w], SPLINE));
1179+
dcomplex refractive_index = dcomplex(refractive_index_real.getValue(wavelength_list[w], LOGLINEAR),
1180+
refractive_index_imag.getValue(wavelength_list[w], LOGLINEAR));
11811181
#endif
11821182

11831183
// Calculate Mie-scattering
@@ -4333,21 +4333,22 @@ double CDustComponent::calcGoldReductionFactor(const Vector3D & v, const Vector3
43334333
double g = gold_g_factor;
43344334
Vector3D v_proj = (v * B) / B.sq_length() * B;
43354335
double len_vz = v_proj.sq_length();
4336-
double len_vxyz = v.length();
4336+
double len_vxyz = v.sq_length();
43374337
double R, cossq_beta;
43384338

4339-
len_vxyz *= len_vxyz;
4339+
if(len_vxyz == 3 * len_vz)
4340+
return 0.0;
43404341

43414342
if(len_vxyz == len_vz)
4342-
return -0.499999999999;
4343+
return 0.0;
43434344

4344-
if(len_vxyz == 3 * len_vz)
4345-
return 0;
4345+
if(gold_g_factor == 0)
4346+
return 0.0;
43464347

43474348
s = -0.5 * (len_vxyz - 3 * len_vz) / (len_vxyz - len_vz);
43484349

4349-
if(s <= -0.5)
4350-
s = -0.4999999999;
4350+
// if(s <= -0.5)
4351+
// s = -0.4999999999;
43514352

43524353
if(s < 0)
43534354
{
@@ -4378,8 +4379,8 @@ double CDustComponent::calcGoldReductionFactor(const Vector3D & v, const Vector3
43784379

43794380
R = 1.5 * cossq_beta - 0.5;
43804381

4381-
if(R <= -0.5)
4382-
R = -0.499999999999;
4382+
// if(R <= -0.5)
4383+
// R = -0.499999999999;
43834384

43844385
return R;
43854386
}

‎src/Grid.cpp

+47-42
Original file line numberDiff line numberDiff line change
@@ -2675,63 +2675,68 @@ bool CGridBasic::writeMidplaneFits(string data_path, parameters & param, uint bi
26752675
}
26762676

26772677
bool CGridBasic::getPolarRTGridParameterWorker(double max_len,
2678-
double pixel_width,
2679-
uint max_subpixel_lvl,
2680-
dlist & _listR,
2681-
uint & N_polar_r,
2682-
uint *& N_polar_ph,
2683-
const uint &N_r,
2684-
const double *listR
2678+
double pixel_width,
2679+
uint max_subpixel_lvl,
2680+
dlist & _listR,
2681+
uint & N_polar_r,
2682+
uint *& N_polar_ph,
2683+
const uint &N_r,
2684+
const double *listR
26852685
)
26862686
{
26872687
uint subpixel_multiplier = pow(2, max_subpixel_lvl);
2688-
2688+
26892689
// Add polar detector pixels in the inner region to obtain resolution specified by max_subpixel_lvl
26902690
// inner grid cell diameter is 2.*listR[0]
2691-
uint N_r_inner = uint(ceil(subpixel_multiplier * 2.0*listR[0]/pixel_width));
2692-
2691+
uint N_r_inner = uint(ceil(subpixel_multiplier * 2.0 * listR[0] / pixel_width));
2692+
26932693
for(uint i_r = 0; i_r <= N_r_inner; i_r++)
2694-
_listR.push_back(listR[0] * (i_r / double(N_r_inner)));
2694+
_listR.push_back(listR[0] * (i_r / double(N_r_inner)));
26952695

26962696
for(uint i_r = 1; i_r <= N_r; i_r++)
26972697
{
2698-
double r1 = _listR[_listR.size() - 1];
2699-
double r2 = listR[i_r];
2700-
// r2-r1 is width of current grid cell's ring
2701-
uint N_r_sub = uint(ceil(subpixel_multiplier * (r2 - r1) / pixel_width));
2702-
2703-
for(uint i_r_sub = 1; i_r_sub <= N_r_sub; i_r_sub++)
2704-
_listR.push_back(r1 + (r2 - r1) * i_r_sub / double(N_r_sub));
2698+
double r1 = _listR[_listR.size() - 1];
2699+
double r2 = listR[i_r];
27052700

2706-
// break if sidelength is smaller than full grid
2707-
if(_listR.back() > max_len)
2708-
{
2709-
_listR.pop_back();
2710-
_listR.push_back(max_len);
2711-
break;
2712-
}
2701+
// if sidelength is smaller than full grid, only consider visible grid
2702+
if(r2 > max_len)
2703+
r2 = max_len;
2704+
2705+
// r2 - r1 is width of current grid cell's ring
2706+
uint N_r_sub = uint(ceil(subpixel_multiplier * (r2 - r1) / pixel_width));
2707+
2708+
for(uint i_r_sub = 1; i_r_sub <= N_r_sub; i_r_sub++)
2709+
_listR.push_back(r1 + (r2 - r1) * i_r_sub / double(N_r_sub));
2710+
2711+
// break if sidelength is smaller than full grid
2712+
if(_listR.back() >= max_len)
2713+
{
2714+
_listR.pop_back();
2715+
_listR.push_back(max_len);
2716+
break;
2717+
}
27132718
}
27142719

27152720
if(_listR.back() < max_len)
27162721
{
2717-
// Create additional outer rings with outermost grid cell radial distance
2718-
// and store them in buffer to do subpixeling afterwards
2719-
uint N_r_outer = uint(ceil((max_len - listR[N_r]) / (listR[N_r] - listR[N_r - 1])));
2720-
std::vector<double> outerR_buffer;
2721-
2722-
for(uint i_r = 1; i_r <= N_r_outer; i_r++)
2723-
outerR_buffer.push_back(listR[N_r] + (max_len - listR[N_r]) * i_r / double(N_r_outer));
2722+
// Create additional outer rings with outermost grid cell radial distance
2723+
// and store them in buffer to do subpixeling afterwards
2724+
uint N_r_outer = uint(ceil((max_len - listR[N_r]) / (listR[N_r] - listR[N_r - 1])));
2725+
std::vector<double> outerR_buffer;
2726+
2727+
for(uint i_r = 1; i_r <= N_r_outer; i_r++)
2728+
outerR_buffer.push_back(listR[N_r] + (max_len - listR[N_r]) * i_r / double(N_r_outer));
27242729

2725-
// loop over equally spaced rings outside the grid and do subpixeling
2726-
for (uint i_r = 0; i_r < N_r_outer; i_r++){
2727-
2728-
double r1 = _listR[_listR.size() - 1];
2729-
double r2 = outerR_buffer[i_r];
2730-
uint N_r_sub = uint(ceil(subpixel_multiplier * (r2 - r1) / pixel_width));
2730+
// loop over equally spaced rings outside the grid and do subpixeling
2731+
for (uint i_r = 0; i_r < N_r_outer; i_r++)
2732+
{
2733+
double r1 = _listR[_listR.size() - 1];
2734+
double r2 = outerR_buffer[i_r];
2735+
uint N_r_sub = uint(ceil(subpixel_multiplier * (r2 - r1) / pixel_width));
27312736

2732-
for(uint i_r_sub = 1; i_r_sub <= N_r_sub; i_r_sub++)
2733-
_listR.push_back(r1 + (r2 - r1) * i_r_sub / double(N_r_sub));
2734-
}
2737+
for(uint i_r_sub = 1; i_r_sub <= N_r_sub; i_r_sub++)
2738+
_listR.push_back(r1 + (r2 - r1) * i_r_sub / double(N_r_sub));
2739+
}
27352740
}
27362741

27372742
// Set total size of the radial cells
@@ -2741,7 +2746,7 @@ bool CGridBasic::getPolarRTGridParameterWorker(double max_len,
27412746
N_polar_ph = new uint[N_polar_r];
27422747
for(uint i_r = 0; i_r < N_polar_r; i_r++)
27432748
{
2744-
N_polar_ph[i_r] = uint(ceil(PIx2 * _listR[i_r + 1] / (_listR[i_r + 1] - _listR[i_r])));
2749+
N_polar_ph[i_r] = uint(ceil(PIx2 * _listR[i_r + 1] / (_listR[i_r + 1] - _listR[i_r])));
27452750
}
27462751

27472752
return true;

‎src/Grid.h

+7-7
Original file line numberDiff line numberDiff line change
@@ -3962,13 +3962,13 @@ class CGridBasic
39623962
uint max_wavelengths;
39633963

39643964
bool getPolarRTGridParameterWorker(double max_len,
3965-
double pixel_width,
3966-
uint max_subpixel_lvl,
3967-
dlist & _listR,
3968-
uint & N_polar_r,
3969-
uint *& N_polar_phi,
3970-
const uint &N_r,
3971-
const double *listR);
3965+
double pixel_width,
3966+
uint max_subpixel_lvl,
3967+
dlist & _listR,
3968+
uint & N_polar_r,
3969+
uint *& N_polar_phi,
3970+
const uint &N_r,
3971+
const double *listR);
39723972
};
39733973

39743974
#endif

‎src/MathFunctions.h

+28-10
Original file line numberDiff line numberDiff line change
@@ -1875,28 +1875,46 @@ class CMathFunctions
18751875

18761876
static inline void LinearList(double start, double stop, double * list, uint N)
18771877
{
1878-
double dx = (stop - start) / (N - 1);
1878+
list[0] = start;
1879+
1880+
if(N > 1){
1881+
double dx = (stop - start) / (N - 1);
1882+
1883+
for(uint i_x = 1; i_x < N - 1; i_x++)
1884+
list[i_x] = start + i_x * dx;
1885+
1886+
list[N - 1] = stop;
1887+
}
1888+
}
18791889

1890+
static inline void LinearList(double start, double stop, dlist & list)
1891+
{
1892+
uint N = list.size();
18801893
list[0] = start;
18811894

1882-
for(uint i_x = 1; i_x < N - 1; i_x++)
1883-
list[i_x] = start + i_x * dx;
1895+
if(N > 1){
1896+
double dx = (stop - start) / (N - 1);
18841897

1885-
list[N - 1] = stop;
1898+
for(uint i_x = 1; i_x < N - 1; i_x++)
1899+
list[i_x] = start + i_x * dx;
1900+
1901+
list[N - 1] = stop;
1902+
}
18861903
}
18871904

18881905
static inline dlist LinearList(double start, double stop, uint N)
18891906
{
18901907
dlist list(N);
1891-
1892-
double dx = (stop - start) / (N - 1);
1893-
18941908
list[0] = start;
1909+
1910+
if(N > 1){
1911+
double dx = (stop - start) / (N - 1);
18951912

1896-
for(uint i_x = 1; i_x < N - 1; i_x++)
1897-
list[i_x] = start + i_x * dx;
1913+
for(uint i_x = 1; i_x < N - 1; i_x++)
1914+
list[i_x] = start + i_x * dx;
18981915

1899-
list[N - 1] = stop;
1916+
list[N - 1] = stop;
1917+
}
19001918

19011919
return list;
19021920
}

‎src/OcTree.cpp

+7
Original file line numberDiff line numberDiff line change
@@ -1723,6 +1723,13 @@ bool CGridOcTree::goToNextCellBorder(photon_package * pp)
17231723
}
17241724

17251725
pp->setPosition(pos + dir * path_length);
1726+
1727+
if(pos == pp->getPosition())
1728+
{
1729+
cout << "\nERROR: Could not transfer photon to the next cell border! " << endl;
1730+
return false;
1731+
}
1732+
17261733
pp->setTmpPathLength(path_length);
17271734

17281735
return true;

‎src/Raytracing.h

+3-3
Original file line numberDiff line numberDiff line change
@@ -1649,7 +1649,7 @@ class CRaytracingPolar : public CRaytracingBasic
16491649
npix_total += npix_ph[i_r];
16501650
}
16511651

1652-
npix_total++; // last pixel is central grid cell
1652+
npix_total++; // last pixel is central grid cell
16531653

16541654
if(npix_total > MAX_RT_RAYS)
16551655
{
@@ -1777,11 +1777,11 @@ class CRaytracingPolar : public CRaytracingBasic
17771777
switch(processing_method)
17781778
{
17791779
case 0:
1780-
cout << "HINT: Using 'nearest' method to post process from polar to cartesian detector" << endl;
1780+
cout << "HINT: Using 'nearest' method to post process from polar to cartesian detector" << endl;
17811781
return postProcessingUsingNearest();
17821782
break;
17831783
case 1:
1784-
cout << "HINT: Using 'interpolation' method to post process from polar to cartesian detector" << endl;
1784+
cout << "HINT: Using 'interpolation' method to post process from polar to cartesian detector" << endl;
17851785
return postProcessingUsingInterpolation();
17861786
break;
17871787
default:

‎src/Spherical.cpp

+7
Original file line numberDiff line numberDiff line change
@@ -1613,6 +1613,13 @@ bool CGridSpherical::goToNextCellBorder(photon_package * pp)
16131613
}
16141614

16151615
pp->setPosition(p + d * path_length);
1616+
1617+
if(p == pp->getPosition())
1618+
{
1619+
cout << "\nERROR: Could not transfer photon to the next cell border! " << endl;
1620+
return false;
1621+
}
1622+
16161623
pp->setTmpPathLength(path_length);
16171624
pp->setDirectionID(dirID);
16181625
return true;

‎src/Spherical.h

+8-8
Original file line numberDiff line numberDiff line change
@@ -457,14 +457,14 @@ class CGridSpherical : public CGridBasic
457457
uint & N_polar_r,
458458
uint *& N_polar_ph)
459459
{
460-
return CGridBasic::getPolarRTGridParameterWorker(max_len,
461-
pixel_width,
462-
max_subpixel_lvl,
463-
_listR,
464-
N_polar_r,
465-
N_polar_ph,
466-
N_r,
467-
listR);
460+
return CGridBasic::getPolarRTGridParameterWorker(max_len,
461+
pixel_width,
462+
max_subpixel_lvl,
463+
_listR,
464+
N_polar_r,
465+
N_polar_ph,
466+
N_r,
467+
listR);
468468
}
469469

470470
private:

‎src/Typedefs.h

+1-1
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,7 @@ using namespace std;
1313

1414
// Header and Version of POLARIS
1515
#define PROG_ID "POLARIS: POLArized RadIation Simulator"
16-
#define VERS_ID " Version 4.12.01 "
16+
#define VERS_ID " Version 4.12.02 "
1717
#define COPY_ID " Copyright (C) 2018 Stefan Reissl "
1818

1919
// Flags to activate WINDOWS support, some DEBUG messages, BENCHMARK settings

‎src/Vector.h

+1-1
Original file line numberDiff line numberDiff line change
@@ -313,7 +313,7 @@ class Vector3D
313313
return tmp;
314314
}
315315

316-
Vector3D operator/(double val)
316+
Vector3D operator/(double val) const
317317
{
318318
Vector3D tmp(x / val, y / val, z / val);
319319
return tmp;

0 commit comments

Comments
 (0)