Merge branch 'develop-basic' into master-basic

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.11.01-bf0040)](https://img.shields.io/badge/Version-4.11.01-bf0040)
+[![Version](https://img.shields.io/badge/Version-4.11.02-bf0040)](https://img.shields.io/badge/Version-4.11.02-bf0040)
 
 is a 3D Monte Carlo radiative transfer code that
 

projects/CommandList.cmd

@@ -11,7 +11,7 @@
     <dust_component>	"/PATH/TO/POLARIS/input/dust/silicate_d03.nk" "plaw" 1.0 3500.0 5e-09 2.5e-07 -3.5
 
     # considered alignment mechanism of non-spherical dust grains (DUST_EMISSION)
-    <align> ALIG_PA or ALIG_RAT or ALIG_IDG or ALIG_GOLD or ALIG_INTERNAL
+    <align> ALIG_PA or ALIG_RAT or ALIG_IDG or ALIG_GOLD or ALIG_INTERNAL or ALIG_NONPA
 
     # enables/disables sublimation of dust grains, sublimation temperature is set in the dust catalog (TEMP, TEMP_RAT)
     <sub_dust> 1 (yes) or 0 (no) 
@@ -28,8 +28,8 @@
     # exponent for the radiative torque efficiency power law (RAT, TEMP_RAT)
     <alpha_Q> 3.0
 
-    # Rayleigh reduction factor for RAT alignment if imperfect internal alignment is not used (TEMP, RAT, TEMP_RAT, DUST_EMISSION)
-    <R_rat> 1.0
+    # Rayleigh reduction factor for not-perfect alignment or RAT alignment if imperfect internal alignment is not used (TEMP, RAT, TEMP_RAT, DUST_EMISSION)
+    <R_rayleigh> 1.0
 
     # overwrite the gas temperature with the dust temperature times TEMP_FACTOR (TEMP, TEMP_RAT)
     <adj_tgas> 1.0
@@ -151,7 +151,7 @@
     <detector_sync_healpix nr_sides = "512"> 0.03 0.7 11 1.65e+20 0 0 -80 80 -45 45
 
     # maximum level of subpixeling (DUST_EMISSION, LINE_EMISSION)
-    <max_subpixel_lvl> 1
+    <max_subpixel_lvl> 3
 
     # rotation axis for first and second rotation angle (ALL) axis_x axis_y axis_z
     <axis1>	1	0	0

src/CommandParser.cpp

@@ -2767,6 +2767,12 @@ bool CCommandParser::parseLine(parameters * param, string cmd, string data, uint
             return true;
         }
 
+        if(data.compare("ALIG_NONPA") == 0)
+        {
+            param->addAlignmentMechanism(ALIG_NONPA);
+            return true;
+        }
+
         if(data.compare("ALIG_IDG") == 0)
         {
             param->addAlignmentMechanism(ALIG_IDG);
@@ -3042,9 +3048,9 @@ bool CCommandParser::parseLine(parameters * param, string cmd, string data, uint
         return true;
     }
 
-    if(cmd.compare("<R_rat>") == 0)
+    if(cmd.compare("<R_rayleigh>") == 0)
     {
-        param->setRATReductionFactor(atof(data.c_str()));
+        param->setRayleighReductionFactor(atof(data.c_str()));
         return true;
     }
 

src/Detector.h

@@ -4696,6 +4696,12 @@ class CDetector
                 alignment_string += ", ";
             alignment_string += "PA";
         }
+        if((alignment & ALIG_NONPA) == ALIG_NONPA)
+        {
+            if(alignment_string != "")
+                alignment_string += ", ";
+            alignment_string += "NONPA";
+        }
         if((alignment & ALIG_IDG) == ALIG_IDG)
         {
             if(alignment_string != "")

src/Dust.cpp

@@ -3104,7 +3104,7 @@ bool CDustComponent::add(double ** size_fraction, CDustComponent * comp, uint **
         f_highJ = comp->getFHighJ();
         Q_ref = comp->getQref();
         alpha_Q = comp->getAlphaQ();
-        R_rat = comp->getRATReductionFactor();
+        R_rayleigh = comp->getRayleighReductionFactor();
         larm_f = comp->getLarmF();
 
         // Init dust properties to be filled with grain properties
@@ -3417,6 +3417,42 @@ void CDustComponent::calcPACrossSections(uint a, uint w, cross_sections & cs, do
     cs *= PI * a_eff_2[a];
 }
 
+void CDustComponent::calcNONPACrossSections(uint a, uint w, cross_sections & cs, double theta) const
+{
+    double sinsq_th = sin(theta);
+    sinsq_th *= sinsq_th;
+
+    double c_s_ext = getCext1(a, w);
+    double c_p_ext = getCext2(a, w);
+    double c_av_ext = getCextMean(a, w);
+
+    double cx_ext = c_av_ext + R_rayleigh / 3.0 * (c_s_ext - c_p_ext) * (1 - 3.0 * sinsq_th);
+    double cy_ext = c_av_ext + R_rayleigh / 3.0 * (c_s_ext - c_p_ext);
+
+    cs.Cext = 0.5 * (cx_ext + cy_ext);
+    cs.Cpol = 0.5 * R_rayleigh * (c_s_ext - c_p_ext) * sinsq_th;
+
+    double c_s_abs = getCabs1(a, w);
+    double c_p_abs = getCabs2(a, w);
+    double c_av_abs = getCabsMean(a, w);
+
+    double cx_abs = c_av_abs + R_rayleigh / 3.0 * (c_s_abs - c_p_abs) * (1 - 3.0 * sinsq_th);
+    double cy_abs = c_av_abs + R_rayleigh / 3.0 * (c_s_abs - c_p_abs);
+
+    cs.Cabs = 0.5 * (cx_abs + cy_abs);
+    cs.Cpabs = 0.5 * R_rayleigh * (c_s_abs - c_p_abs) * sinsq_th;
+
+    double c_s_sca = getCsca1(a, w);
+    double c_p_sca = getCsca2(a, w);
+    double c_av_sca = getCscaMean(a, w);
+
+    double cx_sca = c_av_sca + R_rayleigh / 3.0 * (c_s_sca - c_p_sca) * (1 - 3.0 * sinsq_th);
+    double cy_sca = c_av_sca + R_rayleigh / 3.0 * (c_s_sca - c_p_sca);
+
+    cs.Csca = 0.5 * (cx_sca + cy_sca);
+    cs.Ccirc = 0.5 * getCcirc(a, w) * R_rayleigh * sinsq_th;
+}
+
 void CDustComponent::calcCrossSections(CGridBasic * grid,
                                        const photon_package & pp,
                                        uint i_density,
@@ -3443,6 +3479,13 @@ void CDustComponent::calcCrossSections(CGridBasic * grid,
         return;
     }
 
+    // Not perfect alignment
+    if((alignment & ALIG_NONPA) == ALIG_NONPA)
+    {
+        calcNONPACrossSections(a, w, cs, mag_field_theta);
+        return;
+    }
+
     // Init variables
     double Rrat = 0, Rgold = 0, Ridg = 0, Rent = 1;
     double delta = 1;
@@ -3471,7 +3514,7 @@ void CDustComponent::calcCrossSections(CGridBasic * grid,
             if((alignment & ALIG_INTERNAL) == ALIG_INTERNAL)
                 Rrat = f_highJ + (1 - f_highJ) * getInternalRAT();
             else
-                Rrat = R_rat;
+                Rrat = R_rayleigh;
         }
     }
 
@@ -5136,7 +5179,7 @@ bool CDustMixture::createDustMixtures(parameters & param, string path_data, stri
             single_component[i_comp].setFcorr(param.getFcorr());
             single_component[i_comp].setQref(param.getQref());
             single_component[i_comp].setAlphaQ(param.getAlphaQ());
-            single_component[i_comp].setRATReductionFactor(param.getRATReductionFactor());
+            single_component[i_comp].setRayleighReductionFactor(param.getRayleighReductionFactor());
 
             single_component[i_comp].setDelta0(param.getDelta0());
             single_component[i_comp].setLarmF(param.getLarmF());

src/Dust.h

@@ -76,7 +76,7 @@ class CDustComponent
         material_density = 0;
         gold_g_factor = 0;
         dust_mass_fraction = 0;
-        R_rat = 1.0;
+        R_rayleigh = 1.0;
 
         Q_ref = 0.4;
         alpha_Q = 3.0;
@@ -1050,9 +1050,9 @@ class CDustComponent
         return alpha_Q;
     }
 
-    double getRATReductionFactor()
+    double getRayleighReductionFactor()
     {
-        return R_rat;
+        return R_rayleigh;
     }
 
     void setAlignmentMechanism(uint al)
@@ -1735,9 +1735,9 @@ class CDustComponent
         alpha_Q = val;
     }
 
-    void setRATReductionFactor(double val)
+    void setRayleighReductionFactor(double val)
     {
-        R_rat = val;
+        R_rayleigh = val;
     }
 
     double getQrat(uint a, uint w, double theta) const
@@ -2324,6 +2324,7 @@ class CDustComponent
     uint getInteractingDust(CGridBasic * grid, photon_package * pp, CRandomGenerator * rand_gen, uint cross_section = CROSS_ABS) const;
 
     void calcPACrossSections(uint a, uint w, cross_sections & cs, double theta) const;
+    void calcNONPACrossSections(uint a, uint w, cross_sections & cs, double theta) const;
     void calcExtCrossSections(CGridBasic * grid,
                               const photon_package & pp,
                               uint i_density,
@@ -2412,7 +2413,7 @@ class CDustComponent
     double gold_g_factor;
     double Q_ref;
     double alpha_Q;
-    double R_rat;
+    double R_rayleigh;
 
     double delta_rat;
     double mu;

src/Grid.h

@@ -2412,7 +2412,10 @@ class CGridBasic
         const cell_basic & cell = *pp.getPositionCell();
         Vector3D tmp_dir(cell.getData(data_pos_vx), cell.getData(data_pos_vy), cell.getData(data_pos_vz));
         // Rotate vector from cell center to position
-        return rotateToCenter(pp, tmp_dir, true, true);
+        // return rotateToCenter(pp, tmp_dir, true, true);
+        // ?! Why should the velocity field be rotated ?!
+
+        return tmp_dir;
     }
 
     Vector3D getVelocityField(const cell_basic & cell) const
@@ -2460,7 +2463,10 @@ class CGridBasic
         const cell_basic & cell = *pp.getPositionCell();
         Vector3D tmp_dir(cell.getData(data_pos_mx), cell.getData(data_pos_my), cell.getData(data_pos_mz));
         // Rotate vector from cell center to position
-        return rotateToCenter(pp, tmp_dir, true, true);
+        // return rotateToCenter(pp, tmp_dir, true, true);
+        // ?! Why should the magnetic field be rotated ?!
+
+        return tmp_dir;
     }
 
     void setMagField(cell_basic * cell, const Vector3D & mag)

src/Parameters.h

@@ -103,7 +103,7 @@ class parameters
         f_highJ = 0.25;
         Q_ref = 0.4;
         alpha_Q = 3.0;
-        R_rat = 1.0;
+        R_rayleigh = 1.0;
 
         f_cor = 0.6;
         adjTgas = 0;
@@ -473,6 +473,11 @@ class parameters
         return (align & ALIG_PA) == ALIG_PA;
     }
 
+    bool getAligNONPA() const
+    {
+        return (align & ALIG_NONPA) == ALIG_NONPA;
+    }
+
     bool getAligIDG() const
     {
         return (align & ALIG_IDG) == ALIG_IDG;
@@ -699,9 +704,9 @@ class parameters
         return alpha_Q;
     }
 
-    double getRATReductionFactor() const
+    double getRayleighReductionFactor() const
     {
-        return R_rat;
+        return R_rayleigh;
     }
 
     double getFcorr() const
@@ -1337,9 +1342,9 @@ class parameters
         alpha_Q = val;
     }
 
-    void setRATReductionFactor(double val)
+    void setRayleighReductionFactor(double val)
     {
-        R_rat = val;
+        R_rayleigh = val;
     }
 
     void setFcorr(double val)
@@ -2351,7 +2356,7 @@ class parameters
     double f_cor;
     double Q_ref;
     double alpha_Q;
-    double R_rat;
+    double R_rayleigh;
 
     double adjTgas;
     double isrf_g_zero;

src/Pipeline.h

@@ -174,8 +174,8 @@ class CPipeline
             cout << "- RAT efficiency ref. (Qref)   : " << param.getQref() << endl;
             cout << "- RAT efficiency exp. (alphaQ) : " << param.getAlphaQ() << endl;
         }
-        cout << "- Rayleigh reduction fact. (R) : " << param.getRATReductionFactor() << endl;
-        cout << "- Align limit prefact. (larmF) : " << param.getLarmF() << endl;
+        else
+            cout << "- None" << endl;
     }
 
     void printPathParameters(parameters & param)
@@ -235,21 +235,32 @@ class CPipeline
             cout << "- No alignment" << endl;
         else if(param.getAligPA())
             cout << "- Perfect alignment" << endl;
+        else if(param.getAligNONPA())
+        {
+            cout << "- Non-perfect alignment" << endl;
+            cout << "    Rayleigh reduction fact. (R) : " << param.getRayleighReductionFactor() << endl;
+        }
         else
         {
             if(param.getAligIDG())
                 cout << "- IDG (paramagnetic alignment)" << endl;
             if(param.getAligRAT())
             {
-                if(param.getRATReductionFactor() == 1)
-                    cout << "- RAT (radiative torques); f_highJ: " << param.getFHighJ() << endl;
-                else
-                    cout << "- RAT (radiative torques); f_highJ: " << param.getFHighJ() << "; R_rat: " << param.getRATReductionFactor() << endl;
+                cout << "- RAT (radiative torques)" << endl;
+                if(!param.getAligINTERNAL())
+                    cout << "    Rayleigh reduction fact. (R) : " << param.getRayleighReductionFactor() << endl;
             }
             if(param.getAligGOLD())
                 cout << "- GOLD (mechanical alignment)" << endl;
             if(param.getAligINTERNAL())
-                cout << "- Internal alignment; f_c: " << param.getFcorr() << endl;
+            {
+                cout << "- Internal alignment" << endl;
+                if(param.getAligIDG() || param.getAligGOLD())
+                    cout << "    Correlation factor (f_c)     : " << param.getFcorr() << endl;
+                if(param.getAligRAT())
+                    cout << "    Fraction at high-J (f_highJ) : " << param.getFHighJ() << endl;
+            }
+            cout << "    Align limit prefact. (larmF) : " << param.getLarmF() << endl;
         }
     }
 

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.11.01            "
+#define VERS_ID "           Version 4.11.02            "
 #define COPY_ID "   Copyright (C) 2018 Stefan Reissl   "
 
 // Flags to activate WINDOWS support, some DEBUG messages, BENCHMARK settings
@@ -174,6 +174,7 @@ using namespace std;
 #define ALIG_RAT 8
 #define ALIG_GOLD 16
 #define ALIG_KRAT 32
+#define ALIG_NONPA 64
 
 #define SUPERTHERMAL_LIMIT 3
 #define MACH_LIMIT 1