merge gitlab:develop-basic into github:master-basic, resolved merge c…

…onflicts
polaris-MCRT · Sep 20, 2023 · a136607 · a136607
2 parents 2aadde4 + a025063
commit a136607
Show file tree

Hide file tree

Showing 12 changed files with 91 additions and 16 deletions.
diff --git a/QUICKSTART.md b/QUICKSTART.md
@@ -96,7 +96,7 @@ $$
 h(r) = h_0 \left( \frac{r}{r_0} \right)^\beta
 $$
 
-Default values: $r_0 = 100\ \mathrm{AU}$, $h_0 = 10\ \mathrm{AU}$, $\alpha = 0.9$, $\beta = 1.1$, inner disk radius $r_\mathrm{in} = 0.1\ \mathrm{AU}$, outer disk radius $r_\mathrm{out} = 100\ \mathrm{AU}$, and total gas mass $M_\mathrm{gas} = 10^{-3}\ \mathrm{M_\odot}$ with a dust to gas mass ratio of $0.01$.
+Default values: $r_0 = 100\ \mathrm{AU}$, $h_0 = 10\ \mathrm{AU}$, $\alpha = 1.8$, $\beta = 1.1$, inner disk radius $r_\mathrm{in} = 0.1\ \mathrm{AU}$, outer disk radius $r_\mathrm{out} = 100\ \mathrm{AU}$, and total gas mass $M_\mathrm{gas} = 10^{-3}\ \mathrm{M_\odot}$ with a dust to gas mass ratio of $0.01$.
 
 **Sphere** with a constant density distribution
 

diff --git a/lib/.gitignore b/lib/.gitignore
@@ -1,2 +1,4 @@
 CCfits/
 cfitsio/
+include/
+lib/
diff --git a/projects/disk/example/dust/POLARIS.cmd b/projects/disk/example/dust/POLARIS.cmd
@@ -26,7 +26,7 @@
 <task> 1
 
 	# detector for dust emission = "N_pixel x N_pixel"> wavelength_min wavelength_max wavelength_nr ID_source rot_1 rot_2 distance
-	<detector_dust nr_pixel = "256*256">	0.00085	0.00085	1	1	30.0	0.0	4.319948614054068e+18
+	<detector_dust_polar nr_pixel = "256*256">	0.00085	0.00085	1	1	30.0	0.0	4.319948614054068e+18
 
 	# background radiation source
 	# background = "N_photons"> scaling_factor temperature q u v rot_1 rot_2

diff --git a/projects/disk/grid.dat b/projects/disk/grid.dat
diff --git a/quickstart.pdf b/quickstart.pdf
diff --git a/src/Dust.cpp b/src/Dust.cpp
@@ -897,6 +897,8 @@ bool CDustComponent::readDustRefractiveIndexFile(parameters & param,
     // Init error check
     bool error = false;
 
+    double max_diff = 0.0;
+
     // Init maximum counter value
     uint max_counter = nr_of_dust_species * nr_of_wavelength;
 
@@ -1067,7 +1069,7 @@ bool CDustComponent::readDustRefractiveIndexFile(parameters & param,
                 delete[] S33_start;
                 delete[] S34_start;
 
-                uint nr_or_scat_theta_final = scat_angle_final.size();
+                uint nr_of_scat_theta_final = scat_angle_final.size();
 
                 // Set missing Efficiencies for other axis
                 Qext2[a][w] = Qext1[a][w];
@@ -1078,15 +1080,23 @@ bool CDustComponent::readDustRefractiveIndexFile(parameters & param,
                 if(scat_theta[a][w] != 0){
                     delete[] scat_theta[a][w];
                 }
-                scat_theta[a][w] = new double[nr_or_scat_theta_final];
+
+                double diff_tmp;
+                for(uint sth = 1; sth < nr_of_scat_theta_final; sth++)
+                {
+                    diff_tmp = (S11_final[sth-1] - S11_final[sth]) / max(S11_final[sth-1], S11_final[sth]);
+                    max_diff = max(diff_tmp, max_diff);
+                }
+
+                scat_theta[a][w] = new double[nr_of_scat_theta_final];
 
                 for(uint inc = 0; inc < nr_of_incident_angles; inc++)
                 {
                     sca_mat[a][w][inc] = new Matrix2D*[nr_of_scat_phi];
                     for(uint sph = 0; sph < nr_of_scat_phi; sph++)
                     {
-                        sca_mat[a][w][inc][sph] = new Matrix2D[nr_or_scat_theta_final];
-                        for(uint sth = 0; sth < nr_or_scat_theta_final; sth++)
+                        sca_mat[a][w][inc][sph] = new Matrix2D[nr_of_scat_theta_final];
+                        for(uint sth = 0; sth < nr_of_scat_theta_final; sth++)
                         {
                             sca_mat[a][w][inc][sph][sth].resize(4, 4);
 
@@ -1112,7 +1122,7 @@ bool CDustComponent::readDustRefractiveIndexFile(parameters & param,
                 S33_final.clear();
                 S34_final.clear();
 
-                nr_of_scat_theta[a][w] = nr_or_scat_theta_final;
+                nr_of_scat_theta[a][w] = nr_of_scat_theta_final;
             }
             else
             {
@@ -1144,6 +1154,16 @@ bool CDustComponent::readDustRefractiveIndexFile(parameters & param,
         cout << "ERROR: Problem with optical properties calculation" << endl;
         return false;
     }
+
+    printIDs();
+    cout << "- calculating optical properties: done          " << endl;
+
+    if(max_diff > 0.5) // arbitrary limit
+    {
+        cout << "WARNING: number of scattering angles might be too low (max diff = " << max_diff << ")" << endl;
+        cout << "if required, increase 'NANG' or decrease 'MAX_MIE_SCA_REL_DIFF' (for x < 100) in src/Typedefs.h " << endl;
+    }
+
     return true;
 }
 
@@ -5187,6 +5207,7 @@ void CDustMixture::printParameters(parameters & param, CGridBasic * grid)
 
     // Show title
     cout << CLR_LINE;
+    cout << SEP_LINE;
     cout << "Dust parameters                                                             "
             "            "
          << endl;

diff --git a/src/Parameters.h b/src/Parameters.h
@@ -111,7 +111,7 @@ class parameters
 
         isrf_path = "";
 
-        max_subpixel_lvl = 1;
+        max_subpixel_lvl = 3;
         midplane_zoom = 1;
         max_dust_component_choice = 0;
 

diff --git a/src/RadiativeTransfer.cpp b/src/RadiativeTransfer.cpp
@@ -2084,6 +2084,14 @@ bool CRadiativeTransfer::calcSyncMapsViaRaytracing(parameters & param)
             // Write results either as text or fits file
             if(!tracer[i_det]->writeSyncResults())
                 return false;
+
+            if(tracer[i_det]->getSubpixelWarning())
+            {
+                cout << "WARNING: level of subpixeling (" << param.getMaxSubpixelLvl() << ") might be too low" << endl;
+                cout << "if required, increase the maximum level of subpixeling with <max_subpixel_lvl> in the command file" << endl;
+            }
+            if(tracer[i_det]->getDetectorShape() == DET_PLANE && (grid->getDataID() == GRID_ID_SPH || grid->getDataID() == GRID_ID_CYL))
+                cout << "HINT: a 'polar' detector should be used for a spherical or cylindrical grid" << endl;
         }
     }
 
@@ -2567,6 +2575,14 @@ bool CRadiativeTransfer::calcPolMapsViaRaytracing(parameters & param)
             // Write results either as text or fits file
             if(!tracer[i_det]->writeDustResults(ray_result_type))
                 return false;
+
+            if(tracer[i_det]->getSubpixelWarning())
+            {
+                cout << "WARNING: level of subpixeling (" << param.getMaxSubpixelLvl() << ") might be too low" << endl;
+                cout << "if required, increase the maximum level of subpixeling with <max_subpixel_lvl> in the command file" << endl;
+            }
+            if(tracer[i_det]->getDetectorShape() == DET_PLANE && (grid->getDataID() == GRID_ID_SPH || grid->getDataID() == GRID_ID_CYL))
+                cout << "HINT: a 'polar' detector should be used for a spherical or cylindrical grid" << endl;
         }
     }
 
@@ -3015,6 +3031,13 @@ bool CRadiativeTransfer::calcOPIATEMapsViaRaytracing(parameters& param)
         if(!tracer[i_det]->writeOpiateResults(op))
             return false;
 
+        if(tracer[i_det]->getSubpixelWarning())
+        {
+            cout << "WARNING: level of subpixeling (" << param.getMaxSubpixelLvl() << ") might be too low" << endl;
+            cout << "if required, increase the maximum level of subpixeling with <max_subpixel_lvl> in the command file" << endl;
+        }
+        if(tracer[i_det]->getDetectorShape() == DET_PLANE && (grid->getDataID() == GRID_ID_SPH || grid->getDataID() == GRID_ID_CYL))
+            cout << "HINT: a 'polar' detector should be used for a spherical or cylindrical grid" << endl;
     }
 
     cout << CLR_LINE;
@@ -3140,6 +3163,14 @@ bool CRadiativeTransfer::calcChMapsViaRaytracing(parameters & param)
             if(!tracer[i_det]->writeLineResults(gas, i_species, i_line))
                 return false;
 
+            if(tracer[i_det]->getSubpixelWarning())
+            {
+                cout << "WARNING: level of subpixeling (" << param.getMaxSubpixelLvl() << ") might be too low" << endl;
+                cout << "if required, increase the maximum level of subpixeling with <max_subpixel_lvl> in the command file" << endl;
+            }
+            if(tracer[i_det]->getDetectorShape() == DET_PLANE && (grid->getDataID() == GRID_ID_SPH || grid->getDataID() == GRID_ID_CYL))
+                cout << "HINT: a 'polar' detector should be used for a spherical or cylindrical grid" << endl;
+
             // Increment index for line RT
             i_det++;
         }

diff --git a/src/Raytracing.h b/src/Raytracing.h
@@ -29,6 +29,8 @@ class CRaytracingBasic
         max_subpixel_lvl = 0;
         nr_extra = 0;
 
+        subpixel_warning = false;
+
         sidelength_x = 0;
         sidelength_y = 0;
 
@@ -259,6 +261,16 @@ class CRaytracingBasic
         return nr_extra;
     }
 
+    uint getDetectorShape()
+    {
+        return rt_detector_shape;
+    }
+
+    bool getSubpixelWarning()
+    {
+        return subpixel_warning;
+    }
+
     virtual double getDistance()
     {
         return distance;
@@ -425,6 +437,8 @@ class CRaytracingBasic
     double rot_angle1, rot_angle2;
     double rad_bubble;
 
+    bool subpixel_warning;
+
     int off_len_x, off_len_y;
 
     bool vel_maps, split_emission;
@@ -797,6 +811,9 @@ class CRaytracingCartesian : public CRaytracingBasic
                 }
             }
         }
+        else
+            subpixel_warning = true;
+
         return subpixel;
     }
 

diff --git a/src/Typedefs.h b/src/Typedefs.h
@@ -325,8 +325,8 @@ using namespace std;
 #else
     // Minimum number of scattering angles between 0° and 90°
     // Actual number might be larger due to adaptive refinement
-    #define NANG 91
-    // limit for adaptive refinement of scat angles
+    #define NANG 181
+    // limit for adaptive refinement of scat angles (for x < 100)
     // if >= 1 -> no refinement
     #define MAX_MIE_SCA_REL_DIFF 1
 #endif

diff --git a/tools/polaris_tools_modules/math.py b/tools/polaris_tools_modules/math.py
@@ -385,8 +385,9 @@ def cartesian_to_cylindrical(cartesian_coord):
         """
         cylindrical_coord = np.zeros(3)
         cylindrical_coord[0] = np.linalg.norm(cartesian_coord[0:2])
-        cylindrical_coord[1] = np.arctan2(
-            cartesian_coord[1], cartesian_coord[0]) + np.pi
+        cylindrical_coord[1] = np.arctan2(cartesian_coord[1], cartesian_coord[0])
+        if cylindrical_coord[1] < 0.0:
+            cylindrical_coord[1] += 2.0 * np.pi
         cylindrical_coord[2] = cartesian_coord[2]
         return cylindrical_coord
 
@@ -525,9 +526,10 @@ def kepler_rotation(self, position, stellar_mass):
 
     @staticmethod
     def default_disk_density(position, inner_radius, outer_radius, ref_scale_height=10. * 149597870700.0,
-                             ref_radius=100. * 149597870700.0, alpha=0.9, beta=0.8, tapered_gamma=None,
+                             ref_radius=100. * 149597870700.0, alpha=1.8, beta=1.1, tapered_gamma=None,
                              column_dens_exp=None, real_zero=True):
         """Shakura and Sunyaev disk density profile.
+        Shakura & Sunyaev 1973, https://ui.adsabs.harvard.edu/abs/1973A%2526A....24..337S
 
         Args:
             position (List[float, float, float]): Position in model space.
@@ -568,7 +570,7 @@ def default_disk_density(position, inner_radius, outer_radius, ref_scale_height=
         return density
 
     @staticmethod
-    def default_disk_scale_height(radius, beta=0.8,
+    def default_disk_scale_height(radius, beta=1.1,
                                   ref_scale_height=10. * 149597870700.0, ref_radius=100. * 149597870700.0):
         """Shakura and Sunyaev disk density profile.
 

diff --git a/tools/polaris_tools_modules/model.py b/tools/polaris_tools_modules/model.py
@@ -173,8 +173,10 @@ def __init__(self):
         # Default disk parameter
         self.parameter['ref_radius'] = 100. * self.math.const['au']
         self.parameter['ref_scale_height'] = 10. * self.math.const['au']
-        self.parameter['alpha'] = 0.9 # Andrews et al. 2010, https://ui.adsabs.harvard.edu/abs/2010ApJ...723.1241A/abstract
-        self.parameter['beta'] = 1.1 # Woitke et al. 2019, http://adsabs.harvard.edu/abs/2019PASP..131f4301W
+        self.parameter['beta'] = 1.1
+        # Woitke et al. 2019, http://adsabs.harvard.edu/abs/2019PASP..131f4301W
+        self.parameter['alpha'] = 3.0 * (self.parameter['beta'] - 0.5)
+        # Shakura & Sunyaev 1973, https://ui.adsabs.harvard.edu/abs/1973A%26A....24..337S
 
     def update_parameter(self, extra_parameter):
         """Use this function to set model parameter with the extra parameters.