Merge branch 'master' into docs/fix_changelog

tardis/io/atomic.py

@@ -452,7 +452,7 @@ def _create_collision_coefficient_matrix(self):
                 C_ul_matrix[level_number_lower, level_number_upper, :] = line.values[2:]
                 delta_E_matrix[level_number_lower, level_number_upper] = line['delta_e']
                 #TODO TARDISATOMIC fix change the g_ratio to be the otherway round - I flip them now here.
-                g_ratio_matrix[level_number_lower, level_number_upper] = line['g_ratio']
+                g_ratio_matrix[level_number_lower, level_number_upper] = 1/line['g_ratio']
             self.C_ul_interpolator[species] = interpolate.interp1d(
                     self.atom_data.collision_data_temperatures,
                     C_ul_matrix)

tardis/io/schemas/spectrum.yml

@@ -45,3 +45,8 @@ properties:
         type: number
         default: 0.0
         description: Probability for not terminating the packet path
+      enable_biasing:
+        type: boolean
+        default: False
+        description: If True bias v-packet emission based on the electron
+            scattering optical depth

tardis/montecarlo/montecarlo.pyx

@@ -120,6 +120,8 @@ cdef extern from "src/cmontecarlo.h":
         int full_relativity
         double survival_probability
         double tau_russian
+        double *tau_bias
+        int enable_biasing
 
     void montecarlo_main_loop(storage_model_t * storage, int_type_t virtual_packet_flag, int nthreads, unsigned long seed)
 
@@ -259,6 +261,18 @@ cdef initialize_storage_model(model, plasma, runner, storage_model_t *storage):
 
     storage.tau_russian = runner.v_packet_settings['tau_russian']
     storage.survival_probability = runner.v_packet_settings['survival_probability']
+    storage.enable_biasing = runner.v_packet_settings['enable_biasing']
+
+    if runner.v_packet_settings['enable_biasing']:
+        # Calculate the integrated electron scattering optical depth
+        # at all cell interfaces.
+        runner.tau_bias = np.zeros(len(runner.r_inner_cgs) + 1)
+        runner.tau_bias[:-1] = (
+            ((runner.r_outer_cgs - runner.r_inner_cgs) *
+            plasma.electron_densities.values *
+            runner.sigma_thomson.cgs.value)[::-1].cumsum()[::-1]
+        )
+        storage.tau_bias = <double*> PyArray_DATA(runner.tau_bias)
 
     # Data for continuum implementation
     cdef np.ndarray[double, ndim=1] t_electrons = plasma.t_electrons

tardis/montecarlo/src/cmontecarlo.c

@@ -832,7 +832,7 @@ montecarlo_thomson_scatter (rpacket_t * packet, storage_model_t * storage,
 
   if (rpacket_get_virtual_packet_flag (packet) > 0)
     {
-      montecarlo_one_packet (storage, packet, 1, mt_state);
+      create_vpacket (storage, packet, mt_state);
     }
 }
 
@@ -1011,7 +1011,7 @@ line_emission (rpacket_t * packet, storage_model_t * storage, int64_t emission_l
 	  // QUESTIONABLE!!!
 	  bool old_close_line = rpacket_get_close_line (packet);
 	  rpacket_set_close_line (packet, virtual_close_line);
-	  montecarlo_one_packet (storage, packet, 1, mt_state);
+	  create_vpacket (storage, packet, mt_state);
 	  rpacket_set_close_line (packet, old_close_line);
 	  virtual_close_line = false;
     }
@@ -1052,7 +1052,7 @@ continuum_emission (rpacket_t * packet, storage_model_t * storage, rk_state *mt_
 
   if (rpacket_get_virtual_packet_flag (packet) > 0)
     {
-      montecarlo_one_packet (storage, packet, 1, mt_state);
+      create_vpacket (storage, packet, mt_state);
     }
 }
 
@@ -1264,3 +1264,28 @@ montecarlo_main_loop(storage_model_t * storage, int64_t virtual_packet_flag, int
   print_progress(storage->no_of_packets, storage->no_of_packets);
   fprintf(stderr,"\n");
 }
+
+void
+create_vpacket (storage_model_t * storage, rpacket_t * packet,
+                rk_state *mt_state)
+{
+  if (storage->enable_biasing)
+    {
+      int64_t shell_id = rpacket_get_current_shell_id(packet);
+      double tau_bias = (storage->tau_bias[shell_id + 1] +
+          (storage->tau_bias[shell_id] - storage->tau_bias[shell_id + 1]) *
+          (storage->r_outer[shell_id] - rpacket_get_r (packet)) /
+          (storage->r_outer[shell_id] - storage->r_inner[shell_id]));
+      double vpacket_prob = exp(-tau_bias);
+      double event_random = rk_double (mt_state);
+      if (event_random < vpacket_prob)
+        {
+           packet->vpacket_weight = 1. / vpacket_prob;
+           montecarlo_one_packet (storage, packet, 1, mt_state);
+        }
+    }
+  else
+    {
+      montecarlo_one_packet (storage, packet, 1, mt_state);
+    }
+}

tardis/montecarlo/src/cmontecarlo.h

@@ -152,4 +152,8 @@ double sample_nu_free_bound(const rpacket_t * packet, const storage_model_t * st
 void continuum_emission(rpacket_t * packet, storage_model_t * storage, rk_state *mt_state,
                         pt2sample_nu sample_nu_continuum, int64_t emission_type_id);
 
+void
+create_vpacket (storage_model_t * storage, rpacket_t * packet,
+                rk_state *mt_state);
+
 #endif // TARDIS_CMONTECARLO_H

tardis/montecarlo/src/rpacket.c

@@ -50,5 +50,6 @@ rpacket_init (rpacket_t * packet, storage_model_t * storage, int packet_index,
   rpacket_set_virtual_packet_flag (packet, virtual_packet_flag);
   packet->chi_bf_tmp_partial = chi_bf_tmp_partial;
   packet->compute_chi_bf = true;
+  packet->vpacket_weight = 1.0;
   return ret_val;
 }

tardis/montecarlo/src/rpacket.h

@@ -60,6 +60,7 @@ typedef struct RPacket
   double *chi_bf_tmp_partial;
   int64_t macro_atom_activation_level;
   bool compute_chi_bf;
+  double vpacket_weight;
 } rpacket_t;
 
 static inline double rpacket_get_nu (const rpacket_t * packet)

tardis/montecarlo/src/storage.h

@@ -95,6 +95,8 @@ typedef struct StorageModel
   int full_relativity;
   double survival_probability;
   double tau_russian;
+  double *tau_bias;
+  int enable_biasing;
 } storage_model_t;
 
 #endif // TARDIS_STORAGE_H

tardis/montecarlo/struct.py

@@ -34,7 +34,8 @@ class RPacket(Structure):
         ('chi_bf', c_double),
         ('chi_bf_tmp_partial', POINTER(c_double)),
         ('macro_atom_activation_level', c_int64),
-        ('compute_chi_bf', c_bool)
+        ('compute_chi_bf', c_bool),
+        ('vpacket_weight', c_double)
     ]
 
 
@@ -126,7 +127,11 @@ class StorageModel(Structure):
         ('bf_heating_estimator', POINTER(c_double)),
         ('ff_heating_estimator', POINTER(c_double)),
         ('stim_recomb_cooling_estimator', POINTER(c_double)),
-        ('full_relativity', c_int)
+        ('full_relativity', c_int),
+        ('survival_probability',c_double),
+        ('tau_russian', c_double),
+        ('tau_bias', POINTER(c_double)),
+        ('enable_biasing', c_int)
     ]
 
 

tardis/montecarlo/tests/conftest.py

@@ -46,7 +46,8 @@ def packet():
         id=0,
         chi_cont=6.652486e-16,
         chi_bf_tmp_partial=(c_double * 2)(),
-        compute_chi_bf=True
+        compute_chi_bf=True,
+        vpacket_weight=1.0
     )
 
 
@@ -120,6 +121,11 @@ def model():
         bf_treatment=BoundFreeTreatment.LIN_INTERPOLATION.value,
         ff_heating_estimator=(c_double * 2)(*([0.0] * 2)),
         cont_edge2macro_level=(c_int64 * 6)(*([1] * 6)),
+
+        survival_probability=0.0,
+        tau_russian=10.0,
+        tau_bias=(c_double * 3)(*([5.0, 0.5, 0.0])),
+        enable_biasing=0
     )
 
 

tardis/plasma/properties/partition_function.py

@@ -124,7 +124,7 @@ def __init__(self, plasma_parent):
     def _main_nlte_calculation(
             self, atomic_data, nlte_data, t_electrons,
             j_blues, beta_sobolevs, general_level_boltzmann_factor,
-            previous_electron_densities):
+            previous_electron_densities, g):
         """
         The core of the NLTE calculation, used with all possible config.
         options.
@@ -189,14 +189,14 @@ def _main_nlte_calculation(
                                 'collision data?')
                     else:
                         raise e
-                general_level_boltzmann_factor[i].loc[species] = \
-                    level_boltzmann_factor
+                general_level_boltzmann_factor[i].ix[species] = \
+                    level_boltzmann_factor * g.loc[species][0] / level_boltzmann_factor[0]
         return general_level_boltzmann_factor
 
     def _calculate_classical_nebular(
             self, t_electrons, lines, atomic_data,
             nlte_data, general_level_boltzmann_factor, j_blues,
-            previous_electron_densities):
+            previous_electron_densities, g):
         """
         Performs NLTE calculations using the classical nebular treatment.
         All beta sobolev values taken as 1.
@@ -210,13 +210,13 @@ def _calculate_classical_nebular(
                 j_blues,
                 beta_sobolevs,
                 general_level_boltzmann_factor,
-                previous_electron_densities)
+                previous_electron_densities, g)
         return general_level_boltzmann_factor
 
     def _calculate_coronal_approximation(
             self, t_electrons, lines, atomic_data,
             nlte_data, general_level_boltzmann_factor,
-            previous_electron_densities):
+            previous_electron_densities, g):
         """
         Performs NLTE calculations using the coronal approximation.
         All beta sobolev values taken as 1 and j_blues taken as 0.
@@ -226,13 +226,13 @@ def _calculate_coronal_approximation(
         general_level_boltzmann_factor = self._main_nlte_calculation(
             atomic_data, nlte_data, t_electrons, j_blues,
             beta_sobolevs, general_level_boltzmann_factor,
-            previous_electron_densities)
+            previous_electron_densities, g)
         return general_level_boltzmann_factor
 
     def _calculate_general(
             self, t_electrons, lines, atomic_data, nlte_data,
             general_level_boltzmann_factor, j_blues,
-            previous_beta_sobolev, previous_electron_densities):
+            previous_beta_sobolev, previous_electron_densities, g):
         """
         Full NLTE calculation without approximations.
         """
@@ -244,7 +244,7 @@ def _calculate_general(
         general_level_boltzmann_factor = self._main_nlte_calculation(
             atomic_data, nlte_data, t_electrons, j_blues,
             beta_sobolevs, general_level_boltzmann_factor,
-            previous_electron_densities)
+            previous_electron_densities, g)
         return general_level_boltzmann_factor