- 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
polaris-MCRT · May 3, 2024 · 73ec0b2 · 73ec0b2
1 parent c0b4421
commit 73ec0b2
Show file tree

Hide file tree

Showing 14 changed files with 144 additions and 92 deletions.
diff --git a/README.md b/README.md
@@ -5,7 +5,7 @@
 [![bibcode](https://img.shields.io/badge/bibcode-2016A%26A...593A..87R-1c459b)](https://ui.adsabs.harvard.edu/abs/2016A&A...593A..87R)
 [![doi](https://img.shields.io/badge/doi-10.1051%2F0004--6361%2F201424930-fab70c)](https://doi.org/10.1051/0004-6361/201424930)
 [![License](https://img.shields.io/badge/License-GPLv3-blue)](https://www.gnu.org/licenses/gpl-3.0)
-[![Version](https://img.shields.io/badge/Version-4.12.01-bf0040)](https://img.shields.io/badge/Version-4.12.01-bf0040)
+[![Version](https://img.shields.io/badge/Version-4.12.02-bf0040)](https://img.shields.io/badge/Version-4.12.02-bf0040)
 
 is a 3D Monte Carlo radiative transfer code that
 

diff --git a/src/Cylindrical.cpp b/src/Cylindrical.cpp
@@ -1527,6 +1527,13 @@ bool CGridCylindrical::goToNextCellBorder(photon_package * pp)
     }
 
     pp->setPosition(p + d * path_length);
+
+    if(p == pp->getPosition())
+    {
+        cout << "\nERROR: Could not transfer photon to the next cell border!   " << endl;
+        return false;
+    }
+
     pp->setTmpPathLength(path_length);
     pp->setDirectionID(dirID);
     return true;

diff --git a/src/Cylindrical.h b/src/Cylindrical.h
@@ -445,14 +445,14 @@ class CGridCylindrical : public CGridBasic
                                  uint & N_polar_r,
                                  uint *& N_polar_ph)
     {
-	return CGridBasic::getPolarRTGridParameterWorker(max_len,
-							 pixel_width,
-							 max_subpixel_lvl,
-							 _listR,
-							 N_polar_r,
-							 N_polar_ph,
-							 N_r,
-							 listR);
+    return CGridBasic::getPolarRTGridParameterWorker(max_len,
+                             pixel_width,
+                             max_subpixel_lvl,
+                             _listR,
+                             N_polar_r,
+                             N_polar_ph,
+                             N_r,
+                             listR);
     }
 
   private:

diff --git a/src/Dust.cpp b/src/Dust.cpp
@@ -1176,8 +1176,8 @@ bool CDustComponent::readDustRefractiveIndexFile(parameters & param,
                 dcomplex refractive_index = dcomplex(refractive_index_real.getValue(wavelength_list[w], LOGLINEAR),
                                                      refractive_index_imag.getValue(wavelength_list[w], LOGLINEAR));
 #else
-                dcomplex refractive_index = dcomplex(refractive_index_real.getValue(wavelength_list[w], SPLINE),
-                                                     refractive_index_imag.getValue(wavelength_list[w], SPLINE));
+                dcomplex refractive_index = dcomplex(refractive_index_real.getValue(wavelength_list[w], LOGLINEAR),
+                                                     refractive_index_imag.getValue(wavelength_list[w], LOGLINEAR));
 #endif
 
                 // Calculate Mie-scattering
@@ -4333,21 +4333,22 @@ double CDustComponent::calcGoldReductionFactor(const Vector3D & v, const Vector3
     double g = gold_g_factor;
     Vector3D v_proj = (v * B) / B.sq_length() * B;
     double len_vz = v_proj.sq_length();
-    double len_vxyz = v.length();
+    double len_vxyz = v.sq_length();
     double R, cossq_beta;
 
-    len_vxyz *= len_vxyz;
+    if(len_vxyz == 3 * len_vz)
+        return 0.0;
 
     if(len_vxyz == len_vz)
-        return -0.499999999999;
+        return 0.0;
 
-    if(len_vxyz == 3 * len_vz)
-        return 0;
+    if(gold_g_factor == 0)
+        return 0.0;
 
     s = -0.5 * (len_vxyz - 3 * len_vz) / (len_vxyz - len_vz);
 
-    if(s <= -0.5)
-        s = -0.4999999999;
+    // if(s <= -0.5)
+    //     s = -0.4999999999;
 
     if(s < 0)
     {
@@ -4378,8 +4379,8 @@ double CDustComponent::calcGoldReductionFactor(const Vector3D & v, const Vector3
 
     R = 1.5 * cossq_beta - 0.5;
 
-    if(R <= -0.5)
-        R = -0.499999999999;
+    // if(R <= -0.5)
+    //     R = -0.499999999999;
 
     return R;
 }

diff --git a/src/Grid.cpp b/src/Grid.cpp
@@ -2675,63 +2675,68 @@ bool CGridBasic::writeMidplaneFits(string data_path, parameters & param, uint bi
 }
 
 bool CGridBasic::getPolarRTGridParameterWorker(double max_len,
-					       double pixel_width,
-					       uint max_subpixel_lvl,
-					       dlist & _listR,
-					       uint & N_polar_r,
-					       uint *& N_polar_ph,
-					       const uint &N_r,
-					       const double *listR
+                           double pixel_width,
+                           uint max_subpixel_lvl,
+                           dlist & _listR,
+                           uint & N_polar_r,
+                           uint *& N_polar_ph,
+                           const uint &N_r,
+                           const double *listR
     )
 {
     uint subpixel_multiplier = pow(2, max_subpixel_lvl);
-	
+    
     // Add polar detector pixels in the inner region to obtain resolution specified by max_subpixel_lvl
     // inner grid cell diameter is 2.*listR[0]
-    uint N_r_inner = uint(ceil(subpixel_multiplier * 2.0*listR[0]/pixel_width)); 
-	
+    uint N_r_inner = uint(ceil(subpixel_multiplier * 2.0 * listR[0] / pixel_width)); 
+    
     for(uint i_r = 0; i_r <= N_r_inner; i_r++)
-	_listR.push_back(listR[0] * (i_r / double(N_r_inner)));
+        _listR.push_back(listR[0] * (i_r / double(N_r_inner)));
 
     for(uint i_r = 1; i_r <= N_r; i_r++)
     {
-	double r1 = _listR[_listR.size() - 1];
-	double r2 = listR[i_r];
-	// r2-r1 is width of current grid cell's ring
-	uint N_r_sub = uint(ceil(subpixel_multiplier * (r2 - r1) / pixel_width));
-
-	for(uint i_r_sub = 1; i_r_sub <= N_r_sub; i_r_sub++)
-	    _listR.push_back(r1 + (r2 - r1) * i_r_sub / double(N_r_sub));
+        double r1 = _listR[_listR.size() - 1];
+        double r2 = listR[i_r];
 
-	// break if sidelength is smaller than full grid
-	if(_listR.back() > max_len)
-	{
-	    _listR.pop_back();
-	    _listR.push_back(max_len);
-	    break;
-	}
+        // if sidelength is smaller than full grid, only consider visible grid
+        if(r2 > max_len)
+            r2 = max_len;
+
+        // r2 - r1 is width of current grid cell's ring
+        uint N_r_sub = uint(ceil(subpixel_multiplier * (r2 - r1) / pixel_width));
+
+        for(uint i_r_sub = 1; i_r_sub <= N_r_sub; i_r_sub++)
+            _listR.push_back(r1 + (r2 - r1) * i_r_sub / double(N_r_sub));
+
+        // break if sidelength is smaller than full grid
+        if(_listR.back() >= max_len)
+        {
+            _listR.pop_back();
+            _listR.push_back(max_len);
+            break;
+        }
     }
 
     if(_listR.back() < max_len)
     {
-	// Create additional outer rings with outermost grid cell radial distance
-	// and store them in buffer to do subpixeling afterwards
-	uint N_r_outer = uint(ceil((max_len - listR[N_r]) / (listR[N_r] - listR[N_r - 1])));
-	std::vector<double> outerR_buffer;
-	    
-	for(uint i_r = 1; i_r <= N_r_outer; i_r++)
-	    outerR_buffer.push_back(listR[N_r] + (max_len - listR[N_r]) * i_r / double(N_r_outer));
+        // Create additional outer rings with outermost grid cell radial distance
+        // and store them in buffer to do subpixeling afterwards
+        uint N_r_outer = uint(ceil((max_len - listR[N_r]) / (listR[N_r] - listR[N_r - 1])));
+        std::vector<double> outerR_buffer;
+            
+        for(uint i_r = 1; i_r <= N_r_outer; i_r++)
+            outerR_buffer.push_back(listR[N_r] + (max_len - listR[N_r]) * i_r / double(N_r_outer));
 
-	// loop over equally spaced rings outside the grid and do subpixeling
-	for (uint i_r = 0; i_r < N_r_outer; i_r++){
-		
-	    double r1 = _listR[_listR.size() - 1];
-	    double r2 = outerR_buffer[i_r];
-	    uint N_r_sub = uint(ceil(subpixel_multiplier * (r2 - r1) / pixel_width));  
+        // loop over equally spaced rings outside the grid and do subpixeling
+        for (uint i_r = 0; i_r < N_r_outer; i_r++)
+        {
+            double r1 = _listR[_listR.size() - 1];
+            double r2 = outerR_buffer[i_r];
+            uint N_r_sub = uint(ceil(subpixel_multiplier * (r2 - r1) / pixel_width));  
 
-	    for(uint i_r_sub = 1; i_r_sub <= N_r_sub; i_r_sub++)   
-		_listR.push_back(r1 + (r2 - r1) * i_r_sub / double(N_r_sub));
-	}
+            for(uint i_r_sub = 1; i_r_sub <= N_r_sub; i_r_sub++)   
+                _listR.push_back(r1 + (r2 - r1) * i_r_sub / double(N_r_sub));
+        }
     }
 
     // Set total size of the radial cells
@@ -2741,7 +2746,7 @@ bool CGridBasic::getPolarRTGridParameterWorker(double max_len,
     N_polar_ph = new uint[N_polar_r];
     for(uint i_r = 0; i_r < N_polar_r; i_r++)
     {
-	N_polar_ph[i_r] = uint(ceil(PIx2 * _listR[i_r + 1] / (_listR[i_r + 1] - _listR[i_r])));
+        N_polar_ph[i_r] = uint(ceil(PIx2 * _listR[i_r + 1] / (_listR[i_r + 1] - _listR[i_r])));
     }
 
     return true;

diff --git a/src/Grid.h b/src/Grid.h
@@ -3962,13 +3962,13 @@ class CGridBasic
     uint max_wavelengths;
 
     bool getPolarRTGridParameterWorker(double max_len,
-				       double pixel_width,
-				       uint max_subpixel_lvl,
-				       dlist & _listR,
-				       uint & N_polar_r,
-				       uint *& N_polar_phi,
-				       const uint &N_r,
-				       const double *listR);
+                       double pixel_width,
+                       uint max_subpixel_lvl,
+                       dlist & _listR,
+                       uint & N_polar_r,
+                       uint *& N_polar_phi,
+                       const uint &N_r,
+                       const double *listR);
 };
 
 #endif
diff --git a/src/MathFunctions.h b/src/MathFunctions.h
@@ -1875,28 +1875,46 @@ class CMathFunctions
 
     static inline void LinearList(double start, double stop, double * list, uint N)
     {
-        double dx = (stop - start) / (N - 1);
+        list[0] = start;
+
+        if(N > 1){
+            double dx = (stop - start) / (N - 1);
+
+            for(uint i_x = 1; i_x < N - 1; i_x++)
+                list[i_x] = start + i_x * dx;
+
+            list[N - 1] = stop;
+        }
+    }
 
+    static inline void LinearList(double start, double stop, dlist & list)
+    {
+        uint N = list.size();
         list[0] = start;
 
-        for(uint i_x = 1; i_x < N - 1; i_x++)
-            list[i_x] = start + i_x * dx;
+        if(N > 1){
+            double dx = (stop - start) / (N - 1);
 
-        list[N - 1] = stop;
+            for(uint i_x = 1; i_x < N - 1; i_x++)
+                list[i_x] = start + i_x * dx;
+
+            list[N - 1] = stop;
+        }
     }
 
     static inline dlist LinearList(double start, double stop, uint N)
     {
         dlist list(N);
-
-        double dx = (stop - start) / (N - 1);
-
         list[0] = start;
+
+        if(N > 1){
+            double dx = (stop - start) / (N - 1);
 
-        for(uint i_x = 1; i_x < N - 1; i_x++)
-            list[i_x] = start + i_x * dx;
+            for(uint i_x = 1; i_x < N - 1; i_x++)
+                list[i_x] = start + i_x * dx;
 
-        list[N - 1] = stop;
+            list[N - 1] = stop;
+        }
 
         return list;
     }

diff --git a/src/OcTree.cpp b/src/OcTree.cpp
@@ -1723,6 +1723,13 @@ bool CGridOcTree::goToNextCellBorder(photon_package * pp)
     }
 
     pp->setPosition(pos + dir * path_length);
+
+    if(pos == pp->getPosition())
+    {
+        cout << "\nERROR: Could not transfer photon to the next cell border!   " << endl;
+        return false;
+    }
+
     pp->setTmpPathLength(path_length);
 
     return true;

diff --git a/src/Raytracing.h b/src/Raytracing.h
@@ -1649,7 +1649,7 @@ class CRaytracingPolar : public CRaytracingBasic
             npix_total += npix_ph[i_r];
         }
 
-        npix_total++; // last pixel is central grid cell 
+        npix_total++; // last pixel is central grid cell
 
         if(npix_total > MAX_RT_RAYS)
         {
@@ -1777,11 +1777,11 @@ class CRaytracingPolar : public CRaytracingBasic
         switch(processing_method)
         {
             case 0:
-		cout << "HINT: Using 'nearest' method to post process from polar to cartesian detector" << endl;
+                cout << "HINT: Using 'nearest' method to post process from polar to cartesian detector" << endl;
                 return postProcessingUsingNearest();
                 break;
             case 1:
-		cout << "HINT: Using 'interpolation' method to post process from polar to cartesian detector" << endl;
+                cout << "HINT: Using 'interpolation' method to post process from polar to cartesian detector" << endl;
                 return postProcessingUsingInterpolation();
                 break;
             default:

diff --git a/src/Spherical.cpp b/src/Spherical.cpp
@@ -1613,6 +1613,13 @@ bool CGridSpherical::goToNextCellBorder(photon_package * pp)
     }
 
     pp->setPosition(p + d * path_length);
+
+    if(p == pp->getPosition())
+    {
+        cout << "\nERROR: Could not transfer photon to the next cell border!   " << endl;
+        return false;
+    }
+
     pp->setTmpPathLength(path_length);
     pp->setDirectionID(dirID);
     return true;

diff --git a/src/Spherical.h b/src/Spherical.h
@@ -457,14 +457,14 @@ class CGridSpherical : public CGridBasic
                                  uint & N_polar_r,
                                  uint *& N_polar_ph)
     {
-	return CGridBasic::getPolarRTGridParameterWorker(max_len,
-							 pixel_width,
-							 max_subpixel_lvl,
-							 _listR,
-							 N_polar_r,
-							 N_polar_ph,
-							 N_r,
-							 listR);
+    return CGridBasic::getPolarRTGridParameterWorker(max_len,
+                             pixel_width,
+                             max_subpixel_lvl,
+                             _listR,
+                             N_polar_r,
+                             N_polar_ph,
+                             N_r,
+                             listR);
     }
 
   private:

diff --git a/src/Typedefs.h b/src/Typedefs.h
@@ -13,7 +13,7 @@ using namespace std;
 
 // Header and Version of POLARIS
 #define PROG_ID "POLARIS: POLArized RadIation Simulator"
-#define VERS_ID "           Version 4.12.01            "
+#define VERS_ID "           Version 4.12.02            "
 #define COPY_ID "   Copyright (C) 2018 Stefan Reissl   "
 
 // Flags to activate WINDOWS support, some DEBUG messages, BENCHMARK settings

diff --git a/src/Vector.h b/src/Vector.h
@@ -313,7 +313,7 @@ class Vector3D
         return tmp;
     }
 
-    Vector3D operator/(double val)
+    Vector3D operator/(double val) const
     {
         Vector3D tmp(x / val, y / val, z / val);
         return tmp;