diff --git a/tardis/montecarlo/formal_integral.py b/tardis/montecarlo/formal_integral.py index 872e0410e2c..fe9f82a90de 100644 --- a/tardis/montecarlo/formal_integral.py +++ b/tardis/montecarlo/formal_integral.py @@ -10,7 +10,7 @@ from tardis.montecarlo.montecarlo_numba import njit_dict -# from tardis.montecarlo.montecarlo import formal_integral +from tardis.montecarlo.montecarlo import formal_integral from tardis.montecarlo.spectrum import TARDISSpectrum C_INV = 3.33564e-11 @@ -534,4 +534,4 @@ def intensity_black_body(nu, T): def calculate_p_values(R_max, N, opp): for i in range(N): opp[i] = R_max / (N - 1) * (i) - return opp \ No newline at end of file + return opp diff --git a/tardis/montecarlo/montecarlo.pyx b/tardis/montecarlo/montecarlo.pyx new file mode 100644 index 00000000000..4b398c60425 --- /dev/null +++ b/tardis/montecarlo/montecarlo.pyx @@ -0,0 +1,316 @@ +# cython: profile=False +# cython: boundscheck=False +# cython: wraparound=False +# cython: cdivision=True +# cython: cdivision=True +# cython: language_level=3 + + + +import numpy as np +cimport numpy as np +from numpy cimport PyArray_DATA +from tardis import constants +from astropy import units +from libc.stdlib cimport free + +np.import_array() + + + +ctypedef np.int64_t int_type_t + +cdef extern from "numpy/arrayobject.h": + void PyArray_ENABLEFLAGS(np.ndarray arr, int flags) + +cdef c_array_to_numpy(void *ptr, int dtype, np.npy_intp N): + cdef np.ndarray arr = np.PyArray_SimpleNewFromData(1, &N, dtype, ptr) + PyArray_ENABLEFLAGS(arr, np.NPY_OWNDATA) + return arr + +cdef extern from "src/cmontecarlo.h": + ctypedef enum ContinuumProcessesStatus: + CONTINUUM_OFF = 0 + CONTINUUM_ON = 1 + + cdef int LOG_VPACKETS + + ctypedef struct photo_xsect_1level: + double *nu + double *x_sect + int_type_t no_of_points + + ctypedef struct storage_model_t: + double *packet_nus + double *packet_mus + double *packet_energies + double *output_nus + double *output_energies + double *last_interaction_in_nu + int_type_t *last_line_interaction_in_id + int_type_t *last_line_interaction_out_id + int_type_t *last_line_interaction_shell_id + int_type_t *last_interaction_type + int_type_t *last_interaction_out_type + int_type_t no_of_packets + int_type_t no_of_shells + int_type_t no_of_shells_i + double *r_inner + double *r_outer + double *r_inner_i + double *r_outer_i + double *v_inner + double time_explosion + double inverse_time_explosion + double *electron_densities + double *electron_densities_i + double *inverse_electron_densities + double *line_list_nu + double *line_lists_tau_sobolevs + double *line_lists_tau_sobolevs_i + double *continuum_list_nu + int_type_t line_lists_tau_sobolevs_nd + double *line_lists_j_blues + double *line_lists_Edotlu + int_type_t line_lists_j_blues_nd + int_type_t no_of_lines + int_type_t no_of_edges + int_type_t line_interaction_id + double *transition_probabilities + int_type_t transition_probabilities_nd + int_type_t *line2macro_level_upper + int_type_t *macro_block_references + int_type_t *transition_type + int_type_t *destination_level_id + int_type_t *transition_line_id + double *js + double *nubars + double spectrum_virt_start_nu + double spectrum_virt_end_nu + double spectrum_start_nu + double spectrum_delta_nu + double spectrum_end_nu + double *spectrum_virt_nu + double sigma_thomson + double inverse_sigma_thomson + double inner_boundary_albedo + int_type_t reflective_inner_boundary + photo_xsect_1level ** photo_xsect + double *chi_ff_factor + double *t_electrons + double *l_pop + double *l_pop_r + ContinuumProcessesStatus cont_status + double *virt_packet_nus + double *virt_packet_energies + double *virt_packet_last_interaction_in_nu + int_type_t *virt_packet_last_interaction_type + int_type_t *virt_packet_last_line_interaction_in_id + int_type_t *virt_packet_last_line_interaction_out_id + int_type_t virt_packet_count + int_type_t virt_array_size + int_type_t kpacket2macro_level + int_type_t *cont_edge2macro_level + double *photo_ion_estimator + double *stim_recomb_estimator + int_type_t *photo_ion_estimator_statistics + double *bf_heating_estimator + double *ff_heating_estimator + double *stim_recomb_cooling_estimator + int full_relativity + double survival_probability + double tau_russian + double *tau_bias + int enable_biasing + + +cdef extern from "src/integrator.h": + double *_formal_integral( + const storage_model_t *storage, + double T, + double *nu, + int_type_t nu_size, + double *att_S_ul, + double *Jred_lu, + double *Jblue_lu, + int N) + + + +cdef initialize_storage_model(model, plasma, runner, storage_model_t *storage): + """ + Initializing the storage struct. + + """ + + storage.no_of_packets = runner.input_nu.size + storage.packet_nus = PyArray_DATA(runner.input_nu) + storage.packet_mus = PyArray_DATA(runner.input_mu) + storage.packet_energies = PyArray_DATA(runner.input_energy) + + # Setup of structure + storage.no_of_shells = model.no_of_shells + + + storage.r_inner = PyArray_DATA(runner.r_inner_cgs) + storage.r_outer = PyArray_DATA(runner.r_outer_cgs) + storage.v_inner = PyArray_DATA(runner.v_inner_cgs) + + # Setup the rest + # times + storage.time_explosion = model.time_explosion.to('s').value + storage.inverse_time_explosion = 1.0 / storage.time_explosion + #electron density + storage.electron_densities = PyArray_DATA( + plasma.electron_densities.values) + + runner.inverse_electron_densities = ( + 1.0 / plasma.electron_densities.values) + storage.inverse_electron_densities = PyArray_DATA( + runner.inverse_electron_densities) + # Switch for continuum processes + storage.cont_status = CONTINUUM_OFF + # Continuum data + cdef np.ndarray[double, ndim=1] continuum_list_nu + cdef np.ndarray[double, ndim=1] l_pop + cdef np.ndarray[double, ndim=1] l_pop_r + + if storage.cont_status == CONTINUUM_ON: + continuum_list_nu = np.array([9.0e14, 8.223e14, 6.0e14, 3.5e14, 3.0e14]) # sorted list of threshold frequencies + storage.continuum_list_nu = continuum_list_nu.data + storage.no_of_edges = continuum_list_nu.size + l_pop = np.ones(storage.no_of_shells * continuum_list_nu.size, dtype=np.float64) + storage.l_pop = l_pop.data + l_pop_r = np.ones(storage.no_of_shells * continuum_list_nu.size, dtype=np.float64) + storage.l_pop_r = l_pop_r.data + + # Line lists + storage.no_of_lines = plasma.atomic_data.lines.nu.values.size + storage.line_list_nu = PyArray_DATA(plasma.atomic_data.lines.nu.values) + runner.line_lists_tau_sobolevs = ( + plasma.tau_sobolevs.values.flatten(order='F') + ) + storage.line_lists_tau_sobolevs = PyArray_DATA( + runner.line_lists_tau_sobolevs + ) + storage.line_lists_j_blues = PyArray_DATA( + runner.j_blue_estimator) + + storage.line_lists_Edotlu = PyArray_DATA( + runner.Edotlu_estimator) + + storage.line_interaction_id = runner.get_line_interaction_id( + runner.line_interaction_type) + + # macro atom & downbranch + if storage.line_interaction_id >= 1: + runner.transition_probabilities = ( + plasma.transition_probabilities.values.flatten(order='F') + ) + storage.transition_probabilities = PyArray_DATA( + runner.transition_probabilities + ) + storage.transition_probabilities_nd = ( + plasma.transition_probabilities.values.shape[0]) + storage.line2macro_level_upper = PyArray_DATA( + plasma.atomic_data.lines_upper2macro_reference_idx) + storage.macro_block_references = PyArray_DATA( + plasma.atomic_data.macro_atom_references['block_references'].values) + storage.transition_type = PyArray_DATA( + plasma.atomic_data.macro_atom_data['transition_type'].values) + + # Destination level is not needed and/or generated for downbranch + storage.destination_level_id = PyArray_DATA( + plasma.atomic_data.macro_atom_data['destination_level_idx'].values) + storage.transition_line_id = PyArray_DATA( + plasma.atomic_data.macro_atom_data['lines_idx'].values) + + storage.output_nus = PyArray_DATA(runner._output_nu) + storage.output_energies = PyArray_DATA(runner._output_energy) + + storage.last_line_interaction_in_id = PyArray_DATA( + runner.last_line_interaction_in_id) + storage.last_line_interaction_out_id = PyArray_DATA( + runner.last_line_interaction_out_id) + storage.last_line_interaction_shell_id = PyArray_DATA( + runner.last_line_interaction_shell_id) + storage.last_interaction_type = PyArray_DATA( + runner.last_interaction_type) + storage.last_interaction_in_nu = PyArray_DATA( + runner.last_interaction_in_nu) + + storage.js = PyArray_DATA(runner.j_estimator) + storage.nubars = PyArray_DATA(runner.nu_bar_estimator) + + storage.spectrum_start_nu = runner.spectrum_frequency.to('Hz').value.min() + storage.spectrum_end_nu = runner.spectrum_frequency.to('Hz').value.max() + # TODO: Linspace handling for virtual_spectrum_range + storage.spectrum_virt_start_nu = runner.virtual_spectrum_spawn_range.end.to('Hz', units.spectral()).value + storage.spectrum_virt_end_nu = runner.virtual_spectrum_spawn_range.start.to('Hz', units.spectral()).value + storage.spectrum_delta_nu = runner.spectrum_frequency.to('Hz').value[1] - runner.spectrum_frequency.to('Hz').value[0] + + storage.spectrum_virt_nu = PyArray_DATA( + runner._montecarlo_virtual_luminosity.value) + + storage.sigma_thomson = runner.sigma_thomson + storage.inverse_sigma_thomson = 1.0 / storage.sigma_thomson + storage.reflective_inner_boundary = runner.enable_reflective_inner_boundary + storage.inner_boundary_albedo = runner.inner_boundary_albedo + storage.full_relativity = runner.enable_full_relativity + + 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)[::-1].cumsum()[::-1] + ) + storage.tau_bias = PyArray_DATA(runner.tau_bias) + + # Data for continuum implementation + cdef np.ndarray[double, ndim=1] t_electrons = plasma.t_electrons + storage.t_electrons = t_electrons.data + + +# This will be a method of the Simulation object +def formal_integral(self, nu, N): + cdef storage_model_t storage + + initialize_storage_model(self.model, self.plasma, self.runner, &storage) + + res = self.make_source_function() + + storage.no_of_shells_i = len(self.runner.r_inner_i) + storage.r_inner_i = PyArray_DATA(self.runner.r_inner_i) + storage.r_outer_i = PyArray_DATA(self.runner.r_outer_i) + + storage.electron_densities_i = PyArray_DATA( + self.runner.electron_densities_integ) + self.runner.line_lists_tau_sobolevs_i = ( + self.runner.tau_sobolevs_integ.flatten(order='F') + ) + storage.line_lists_tau_sobolevs_i = PyArray_DATA( + self.runner.line_lists_tau_sobolevs_i + ) + + att_S_ul = res[0].flatten(order='F') + Jred_lu = res[1].flatten(order='F') + Jblue_lu = res[2].flatten(order='F') + + cdef double *L = _formal_integral( + &storage, + self.model.t_inner.value, + PyArray_DATA(nu), + nu.shape[0], + PyArray_DATA(att_S_ul), + PyArray_DATA(Jred_lu), + PyArray_DATA(Jblue_lu), + N + ) + return c_array_to_numpy(L, np.NPY_DOUBLE, nu.shape[0]) diff --git a/tardis/montecarlo/setup_package.py b/tardis/montecarlo/setup_package.py index bfe45232e6b..6f0a2e58bc2 100644 --- a/tardis/montecarlo/setup_package.py +++ b/tardis/montecarlo/setup_package.py @@ -2,7 +2,7 @@ from setuptools import Extension import os from astropy_helpers.distutils_helpers import get_distutils_option -# from Cython.Build import cythonize +from Cython.Build import cythonize from glob import glob @@ -20,27 +20,27 @@ if get_distutils_option('with_vpacket_logging', ['build', 'install', 'develop']) is not None: define_macros.append(('WITH_VPACKET_LOGGING', None)) -# def get_extensions(): -# sources = ['tardis/montecarlo/montecarlo.pyx'] -# sources += [os.path.relpath(fname) for fname in glob( -# os.path.join(os.path.dirname(__file__), 'src', '*.c'))] -# sources += [os.path.relpath(fname) for fname in glob( -# os.path.join(os.path.dirname(__file__), 'src/randomkit', '*.c'))] -# deps = [os.path.relpath(fname) for fname in glob( -# os.path.join(os.path.dirname(__file__), 'src', '*.h'))] -# deps += [os.path.relpath(fname) for fname in glob( -# os.path.join(os.path.dirname(__file__), 'src/randomkit', '*.h'))] - - # return cythonize( - # Extension('tardis.montecarlo.montecarlo', sources, - # include_dirs=['tardis/montecarlo/src', - # 'tardis/montecarlo/src/randomkit', - # 'numpy'], - # depends=deps, - # extra_compile_args=compile_args, - # extra_link_args=link_args, - # define_macros=define_macros) - # ) +def get_extensions(): + sources = ['tardis/montecarlo/montecarlo.pyx'] + sources += [os.path.relpath(fname) for fname in glob( + os.path.join(os.path.dirname(__file__), 'src', '*.c'))] + sources += [os.path.relpath(fname) for fname in glob( + os.path.join(os.path.dirname(__file__), 'src/randomkit', '*.c'))] + deps = [os.path.relpath(fname) for fname in glob( + os.path.join(os.path.dirname(__file__), 'src', '*.h'))] + deps += [os.path.relpath(fname) for fname in glob( + os.path.join(os.path.dirname(__file__), 'src/randomkit', '*.h'))] + + return cythonize( + Extension('tardis.montecarlo.montecarlo', sources, + include_dirs=['tardis/montecarlo/src', + 'tardis/montecarlo/src/randomkit', + 'numpy'], + depends=deps, + extra_compile_args=compile_args, + extra_link_args=link_args, + define_macros=define_macros) + ) def get_package_data(): diff --git a/tardis/montecarlo/src/abbrev.h b/tardis/montecarlo/src/abbrev.h new file mode 100644 index 00000000000..fff7d07036c --- /dev/null +++ b/tardis/montecarlo/src/abbrev.h @@ -0,0 +1,24 @@ +#ifndef ABBREV_H +#define ABBREV_H + +#include + +/** + * @brief safe malloc; checks for NULL and aborts if encountered + */ +static inline void* safe_malloc(size_t size) { + void *mem = malloc(size); + if (mem == NULL && size != 0) abort(); + return mem; +} + +/** + * @brief safe realloc; checks for NULL and aborts if encountered + */ +static inline void* safe_realloc(void *ptr, size_t size) { + void *mem = realloc(ptr, size); + if (mem == NULL && size != 0) abort(); + return mem; +} + +#endif /* ABBREV_H */ diff --git a/tardis/montecarlo/src/cmontecarlo.c b/tardis/montecarlo/src/cmontecarlo.c new file mode 100644 index 00000000000..41441d8d363 --- /dev/null +++ b/tardis/montecarlo/src/cmontecarlo.c @@ -0,0 +1,1291 @@ + +#include +#include +#include +#include +#ifdef WITHOPENMP +#include +#endif +#include "io.h" +#include "abbrev.h" +#include "status.h" +#include "rpacket.h" +#include "cmontecarlo.h" + + +/** Look for a place to insert a value in an inversely sorted float array. + * + * @param x an inversely (largest to lowest) sorted float array + * @param x_insert a value to insert + * @param imin lower bound + * @param imax upper bound + * + * @return index of the next boundary to the left + */ +tardis_error_t +reverse_binary_search (const double *x, double x_insert, + int64_t imin, int64_t imax, int64_t * result) +{ + /* + Have in mind that *x points to a reverse sorted array. + That is large values will have small indices and small ones + will have large indices. + */ + tardis_error_t ret_val = TARDIS_ERROR_OK; + if (x_insert > x[imin] || x_insert < x[imax]) + { + ret_val = TARDIS_ERROR_BOUNDS_ERROR; + } + else + { + int imid = (imin + imax) >> 1; + while (imax - imin > 2) + { + if (x[imid] < x_insert) + { + imax = imid + 1; + } + else + { + imin = imid; + } + imid = (imin + imax) >> 1; + } + if (imax - imin == 2 && x_insert < x[imin + 1]) + { + *result = imin + 1; + } + else + { + *result = imin; + } + } + return ret_val; +} + +/** Insert a value in to an array of line frequencies + * + * @param nu array of line frequencies + * @param nu_insert value of nu key + * @param number_of_lines number of lines in the line list + * + * @return index of the next line ot the red. If the key value is redder than the reddest line returns number_of_lines. + */ +tardis_error_t +line_search (const double *nu, double nu_insert, int64_t number_of_lines, + int64_t * result) +{ + tardis_error_t ret_val = TARDIS_ERROR_OK; + int64_t imin = 0; + int64_t imax = number_of_lines - 1; + if (nu_insert > nu[imin]) + { + *result = imin; + } + else if (nu_insert < nu[imax]) + { + *result = imax + 1; + } + else + { + ret_val = reverse_binary_search (nu, nu_insert, imin, imax, result); + *result = *result + 1; + } + return ret_val; +} + +tardis_error_t +binary_search (const double *x, double x_insert, int64_t imin, + int64_t imax, int64_t * result) +{ + /* + Have in mind that *x points to a sorted array. + Like [1,2,3,4,5,...] + */ + int imid; + tardis_error_t ret_val = TARDIS_ERROR_OK; + if (x_insert < x[imin] || x_insert > x[imax]) + { + ret_val = TARDIS_ERROR_BOUNDS_ERROR; + } + else + { + while (imax >= imin) + { + imid = (imin + imax) / 2; + if (x[imid] == x_insert) + { + *result = imid; + break; + } + else if (x[imid] < x_insert) + { + imin = imid + 1; + } + else + { + imax = imid - 1; + } + } + if (imax - imid == 2 && x_insert < x[imin + 1]) + { + *result = imin; + } + else + { + *result = imin; + } + } + return ret_val; +} + +void +angle_aberration_CMF_to_LF (rpacket_t *packet, const storage_model_t *storage) +{ + if (storage->full_relativity) + { + double beta = rpacket_get_r (packet) * storage->inverse_time_explosion * INVERSE_C; + double mu_0 = rpacket_get_mu (packet); + rpacket_set_mu (packet, (mu_0 + beta) / (1.0 + beta * mu_0)); + } +} + +/** Transform the lab frame direction cosine to the CMF + * + * @param packet + * @param storage + * @param mu lab frame direction cosine + * + * @return CMF direction cosine + */ +double +angle_aberration_LF_to_CMF (rpacket_t *packet, const storage_model_t *storage, double mu) +{ + double beta = rpacket_get_r (packet) * storage->inverse_time_explosion * INVERSE_C; + return (mu - beta) / (1.0 - beta * mu); +} + +double +rpacket_doppler_factor (const rpacket_t *packet, const storage_model_t *storage) +{ + double beta = rpacket_get_r (packet) * storage->inverse_time_explosion * INVERSE_C; + if (!storage->full_relativity) + { + return 1.0 - rpacket_get_mu (packet) * beta; + } + else + { + return (1.0 - rpacket_get_mu (packet) * beta) / sqrt (1 - beta * beta); + } +} + +double +rpacket_inverse_doppler_factor (const rpacket_t *packet, const storage_model_t *storage) +{ + double beta = rpacket_get_r (packet) * storage->inverse_time_explosion * INVERSE_C; + if (!storage->full_relativity) + { + return 1.0 / (1.0 - rpacket_get_mu (packet) * beta); + } + else + { + return (1.0 + rpacket_get_mu (packet) * beta) / sqrt (1 - beta * beta); + } +} + +double +bf_cross_section (const storage_model_t * storage, int64_t continuum_id, double comov_nu) +{ + double bf_xsect; + double *x_sect = storage->photo_xsect[continuum_id]->x_sect; + double *nu = storage->photo_xsect[continuum_id]->nu; + + switch (storage->bf_treatment) + { + case LIN_INTERPOLATION: + { + int64_t result; + tardis_error_t error = binary_search (nu, comov_nu, 0, + storage->photo_xsect[continuum_id]->no_of_points - 1, &result); + if (error == TARDIS_ERROR_BOUNDS_ERROR) + { + bf_xsect = 0.0; + } + else + { + bf_xsect = x_sect[result-1] + (comov_nu - nu[result-1]) / (nu[result] - nu[result-1]) + * (x_sect[result] - x_sect[result-1]); + } + break; + } + + case HYDROGENIC: + { + double nu_ratio = nu[0] / comov_nu; + bf_xsect = x_sect[0] * nu_ratio * nu_ratio * nu_ratio; + break; + } + + default: + fprintf (stderr, "(%d) is not a valid bound-free cross section treatment.\n", storage->bf_treatment); + exit(1); + } + return bf_xsect; +} + +void calculate_chi_bf (rpacket_t * packet, storage_model_t * storage) +{ + double doppler_factor = rpacket_doppler_factor (packet, storage); + double comov_nu = rpacket_get_nu (packet) * doppler_factor; + + int64_t no_of_continuum_edges = storage->no_of_edges; + int64_t current_continuum_id; + line_search(storage->continuum_list_nu, comov_nu, no_of_continuum_edges, ¤t_continuum_id); + rpacket_set_current_continuum_id (packet, current_continuum_id); + + int64_t shell_id = rpacket_get_current_shell_id (packet); + double T = storage->t_electrons[shell_id]; + double boltzmann_factor = exp (-(H * comov_nu) / (KB * T)); + + double bf_helper = 0; + for(int64_t i = current_continuum_id; i < no_of_continuum_edges; i++) + { + // get the level population for the level ijk in the current shell: + double l_pop = storage->l_pop[shell_id * no_of_continuum_edges + i]; + // get the level population ratio \frac{n_{0,j+1,k}}{n_{i,j,k}} \frac{n_{i,j,k}}{n_{0,j+1,k}}^{*}: + double l_pop_r = storage->l_pop_r[shell_id * no_of_continuum_edges + i]; + double bf_x_sect = bf_cross_section (storage, i, comov_nu); + if (bf_x_sect == 0.0) + { + break; + } + bf_helper += l_pop * bf_x_sect * (1.0 - l_pop_r * boltzmann_factor) * doppler_factor; + + packet->chi_bf_tmp_partial[i] = bf_helper; + } + + rpacket_set_chi_boundfree (packet, bf_helper); +} + +void calculate_chi_ff (rpacket_t * packet, const storage_model_t * storage) +{ + double doppler_factor = rpacket_doppler_factor (packet, storage); + double comov_nu = rpacket_get_nu (packet) * doppler_factor; + int64_t shell_id = rpacket_get_current_shell_id (packet); + double T = storage->t_electrons[shell_id]; + double boltzmann_factor = exp (-(H * comov_nu) / KB / T); + double chi_ff_factor = storage->chi_ff_factor[shell_id]; + + double chi_ff = chi_ff_factor * (1 - boltzmann_factor) * pow (comov_nu, -3); + + rpacket_set_chi_freefree (packet, chi_ff * doppler_factor); +} + +void +compute_distance2boundary (rpacket_t * packet, const storage_model_t * storage) +{ + double r = rpacket_get_r (packet); + double mu = rpacket_get_mu (packet); + double r_outer = storage->r_outer[rpacket_get_current_shell_id (packet)]; + double r_inner = storage->r_inner[rpacket_get_current_shell_id (packet)]; + double check, distance; + if (mu > 0.0) + { // direction outward + rpacket_set_next_shell_id (packet, 1); + distance = sqrt (r_outer * r_outer + ((mu * mu - 1.0) * r * r)) - (r * mu); + } + else + { // going inward + if ( (check = r_inner * r_inner + (r * r * (mu * mu - 1.0)) )>= 0.0) + { // hit inner boundary + rpacket_set_next_shell_id (packet, -1); + distance = - r * mu - sqrt (check); + } + else + { // miss inner boundary + rpacket_set_next_shell_id (packet, 1); + distance = sqrt (r_outer * r_outer + ((mu * mu - 1.0) * r * r)) - (r * mu); + } + } + rpacket_set_d_boundary (packet, distance); +} + +tardis_error_t +compute_distance2line (rpacket_t * packet, const storage_model_t * storage) +{ + if (!rpacket_get_last_line (packet)) + { + double r = rpacket_get_r (packet); + double mu = rpacket_get_mu (packet); + double nu = rpacket_get_nu (packet); + double nu_line = rpacket_get_nu_line (packet); + double distance, nu_diff; + double ct = storage->time_explosion * C; + double doppler_factor = rpacket_doppler_factor (packet, storage); + double comov_nu = nu * doppler_factor; + if ( (nu_diff = comov_nu - nu_line) >= 0) + { + if (!storage->full_relativity) + { + distance = (nu_diff / nu) * ct; + } + else + { + double nu_r = nu_line / nu; + distance = - mu * r + (ct - nu_r * nu_r * sqrt(ct * ct - + (1 + r * r * (1 - mu * mu) * (1 + pow (nu_r, -2))))) / (1 + nu_r * nu_r); + } + rpacket_set_d_line (packet, distance); + return TARDIS_ERROR_OK; + } + else + { + if (rpacket_get_next_line_id (packet) == storage->no_of_lines - 1) + { + fprintf (stderr, "last_line = %f\n", + storage-> + line_list_nu[rpacket_get_next_line_id (packet) - 1]); + fprintf (stderr, "Last line in line list reached!"); + } + else if (rpacket_get_next_line_id (packet) == 0) + { + fprintf (stderr, "First line in line list!"); + fprintf (stderr, "next_line = %f\n", + storage-> + line_list_nu[rpacket_get_next_line_id (packet) + 1]); + } + else + { + fprintf (stderr, "last_line = %f\n", + storage-> + line_list_nu[rpacket_get_next_line_id (packet) - 1]); + fprintf (stderr, "next_line = %f\n", + storage-> + line_list_nu[rpacket_get_next_line_id (packet) + 1]); + } + fprintf (stderr, "ERROR: Comoving nu less than nu_line!\n"); + fprintf (stderr, "comov_nu = %f\n", comov_nu); + fprintf (stderr, "nu_line = %f\n", nu_line); + fprintf (stderr, "(comov_nu - nu_line) / nu_line = %f\n", + (comov_nu - nu_line) / nu_line); + fprintf (stderr, "r = %f\n", r); + fprintf (stderr, "mu = %f\n", mu); + fprintf (stderr, "nu = %f\n", nu); + fprintf (stderr, "doppler_factor = %f\n", doppler_factor); + fprintf (stderr, "cur_zone_id = %" PRIi64 "\n", rpacket_get_current_shell_id (packet)); + return TARDIS_ERROR_COMOV_NU_LESS_THAN_NU_LINE; + } + } + else + { + rpacket_set_d_line (packet, MISS_DISTANCE); + return TARDIS_ERROR_OK; + } +} + +void +compute_distance2continuum (rpacket_t * packet, storage_model_t * storage) +{ + double chi_continuum, d_continuum; + double chi_electron = storage->electron_densities[rpacket_get_current_shell_id(packet)] * + storage->sigma_thomson; + if (storage->full_relativity) + { + chi_electron *= rpacket_doppler_factor (packet, storage); + } + + if (storage->cont_status == CONTINUUM_ON) + { + if (packet->compute_chi_bf) + { + calculate_chi_bf (packet, storage); + calculate_chi_ff (packet, storage); + } + else + { + packet->compute_chi_bf=true; + } + chi_continuum = rpacket_get_chi_boundfree (packet) + rpacket_get_chi_freefree (packet) + chi_electron; + d_continuum = rpacket_get_tau_event (packet) / chi_continuum; + } + else + { + chi_continuum = chi_electron; + d_continuum = storage->inverse_electron_densities[rpacket_get_current_shell_id (packet)] * + storage->inverse_sigma_thomson * rpacket_get_tau_event (packet); + } + + if (rpacket_get_virtual_packet(packet) > 0) + { + //Set all continuum distances to MISS_DISTANCE in case of an virtual_packet + d_continuum = MISS_DISTANCE; + packet->compute_chi_bf = false; + } + else + { + + // fprintf(stderr, "--------\n"); + // fprintf(stderr, "nu = %e \n", rpacket_get_nu(packet)); + // fprintf(stderr, "chi_electron = %e\n", chi_electron); + // fprintf(stderr, "chi_boundfree = %e\n", calculate_chi_bf(packet, storage)); + // fprintf(stderr, "chi_line = %e \n", rpacket_get_tau_event(packet) / rpacket_get_d_line(packet)); + // fprintf(stderr, "--------\n"); + + //rpacket_set_chi_freefree(packet, chi_freefree); + rpacket_set_chi_electron (packet, chi_electron); + } + rpacket_set_chi_continuum (packet, chi_continuum); + rpacket_set_d_continuum (packet, d_continuum); +} + +void +macro_atom (rpacket_t * packet, const storage_model_t * storage, rk_state *mt_state) +{ + int emit = 0, i = 0, offset = -1; + uint64_t activate_level = rpacket_get_macro_atom_activation_level (packet); + while (emit >= 0) + { + double event_random = rk_double (mt_state); + i = storage->macro_block_references[activate_level] - 1; + double p = 0.0; + offset = storage->transition_probabilities_nd * + rpacket_get_current_shell_id (packet); + do + { + ++i; + p += storage->transition_probabilities[offset + i]; + } + while (p <= event_random); + emit = storage->transition_type[i]; + activate_level = storage->destination_level_id[i]; + } + switch (emit) + { + case BB_EMISSION: + line_emission (packet, storage, storage->transition_line_id[i], mt_state); + break; + + case BF_EMISSION: + rpacket_set_current_continuum_id (packet, storage->transition_line_id[i]); + storage->last_line_interaction_out_id[rpacket_get_id (packet)] = + rpacket_get_current_continuum_id (packet); + + continuum_emission (packet, storage, mt_state, sample_nu_free_bound, 3); + break; + + case FF_EMISSION: + continuum_emission (packet, storage, mt_state, sample_nu_free_free, 4); + break; + + case ADIABATIC_COOLING: + storage->last_interaction_type[rpacket_get_id (packet)] = 5; + rpacket_set_status (packet, TARDIS_PACKET_STATUS_REABSORBED); + break; + + default: + fprintf (stderr, "This process for macro-atom deactivation should not exist! (emit = %d)\n", emit); + exit(1); + } +} + +void +move_packet (rpacket_t * packet, storage_model_t * storage, double distance) +{ + double doppler_factor = rpacket_doppler_factor (packet, storage); + if (distance > 0.0) + { + double r = rpacket_get_r (packet); + double new_r = + sqrt (r * r + distance * distance + + 2.0 * r * distance * rpacket_get_mu (packet)); + rpacket_set_mu (packet, + (rpacket_get_mu (packet) * r + distance) / new_r); + rpacket_set_r (packet, new_r); + if (rpacket_get_virtual_packet (packet) <= 0) + { + double comov_energy = rpacket_get_energy (packet) * doppler_factor; + double comov_nu = rpacket_get_nu (packet) * doppler_factor; + if (storage->full_relativity) + { + distance *= doppler_factor; + } +#ifdef WITHOPENMP +#pragma omp atomic +#endif + storage->js[rpacket_get_current_shell_id (packet)] += + comov_energy * distance; +#ifdef WITHOPENMP +#pragma omp atomic +#endif + storage->nubars[rpacket_get_current_shell_id (packet)] += + comov_energy * distance * comov_nu; + + if (storage->cont_status) + { + increment_continuum_estimators(packet, storage, distance, comov_nu, comov_energy); + } + } + } +} + +void +increment_continuum_estimators (const rpacket_t * packet, storage_model_t * storage, double distance, + double comov_nu, double comov_energy) +{ + int64_t current_continuum_id; + int64_t no_of_continuum_edges = storage->no_of_edges; + int64_t shell_id = rpacket_get_current_shell_id (packet); + line_search(storage->continuum_list_nu, comov_nu, no_of_continuum_edges, ¤t_continuum_id); + double T = storage->t_electrons[shell_id]; + double boltzmann_factor = exp (-(H * comov_nu) / (KB * T)); + + #ifdef WITHOPENMP + #pragma omp atomic + #endif + storage->ff_heating_estimator[shell_id] += comov_energy * distance * rpacket_get_chi_freefree (packet); + + for(int64_t i = current_continuum_id; i < no_of_continuum_edges; i++) + { + double bf_xsect = bf_cross_section (storage, i, comov_nu); + int64_t photo_ion_idx = i * storage->no_of_shells + shell_id; + double photo_ion_estimator_helper = comov_energy * distance * bf_xsect / comov_nu; + double bf_heating_estimator_helper = + comov_energy * distance * bf_xsect * (1. - storage->continuum_list_nu[i] / comov_nu); + + #ifdef WITHOPENMP + #pragma omp atomic + #endif + storage->photo_ion_estimator[photo_ion_idx] += photo_ion_estimator_helper; + + #ifdef WITHOPENMP + #pragma omp atomic + #endif + storage->stim_recomb_estimator[photo_ion_idx] += photo_ion_estimator_helper * boltzmann_factor; + + #ifdef WITHOPENMP + #pragma omp atomic + #endif + storage->bf_heating_estimator[photo_ion_idx] += bf_heating_estimator_helper; + + #ifdef WITHOPENMP + #pragma omp atomic + #endif + storage->stim_recomb_cooling_estimator[photo_ion_idx] += bf_heating_estimator_helper * boltzmann_factor; + + if (photo_ion_estimator_helper != 0.0) + { + #ifdef WITHOPENMP + #pragma omp atomic + #endif + storage->photo_ion_estimator_statistics[photo_ion_idx] += 1; + } + else + { + break; + } + } +} + +double +get_increment_j_blue_estimator_energy (const rpacket_t * packet, + const storage_model_t * storage, + double d_line) +{ + double energy; + if (storage->full_relativity) + { + // Accurate up to a factor 1 / gamma + energy = rpacket_get_energy (packet); + } + else + { + double r = rpacket_get_r (packet); + double r_interaction = sqrt (r * r + d_line * d_line + + 2.0 * r * d_line * rpacket_get_mu (packet)); + double mu_interaction = (rpacket_get_mu (packet) * r + d_line) / r_interaction; + double doppler_factor = 1.0 - mu_interaction * r_interaction * + storage->inverse_time_explosion * INVERSE_C; + energy = rpacket_get_energy (packet) * doppler_factor; + } + return energy; +} + +void +increment_j_blue_estimator (const rpacket_t * packet, storage_model_t * storage, + double d_line, int64_t j_blue_idx) +{ + if (storage->line_lists_j_blues != NULL) + { + double energy = get_increment_j_blue_estimator_energy (packet, storage, + d_line); + #ifdef WITHOPENMP + #pragma omp atomic + #endif + storage->line_lists_j_blues[j_blue_idx] += + energy / rpacket_get_nu (packet); + } +} + +void +increment_Edotlu_estimator (const rpacket_t * packet, storage_model_t * storage, + double d_line, int64_t line_idx) +{ + if (storage->line_lists_Edotlu != NULL) + { + double energy = get_increment_j_blue_estimator_energy (packet, storage, + d_line); + #ifdef WITHOPENMP + #pragma omp atomic + #endif + storage->line_lists_Edotlu[line_idx] += energy; + } +} + + +int64_t +montecarlo_one_packet (storage_model_t * storage, rpacket_t * packet, + int64_t virtual_mode, rk_state *mt_state) +{ + int64_t reabsorbed=-1; + if (virtual_mode == 0) + { + reabsorbed = montecarlo_one_packet_loop (storage, packet, 0, mt_state); + } + else + { + if ((rpacket_get_nu (packet) > storage->spectrum_virt_start_nu) && (rpacket_get_nu(packet) < storage->spectrum_virt_end_nu)) + { + for (int64_t i = 0; i < rpacket_get_virtual_packet_flag (packet); i++) + { + double weight; + rpacket_t virt_packet = *packet; + double mu_min; + if (rpacket_get_r(&virt_packet) > storage->r_inner[0]) + { + mu_min = + -1.0 * sqrt (1.0 - + (storage->r_inner[0] / rpacket_get_r(&virt_packet)) * + (storage->r_inner[0] / rpacket_get_r(&virt_packet))); + + if (storage->full_relativity) + { + // Need to transform the angular size of the photosphere into the CMF + mu_min = angle_aberration_LF_to_CMF (&virt_packet, storage, mu_min); + } + } + else + { + mu_min = 0.0; + } + double mu_bin = (1.0 - mu_min) / rpacket_get_virtual_packet_flag (packet); + rpacket_set_mu(&virt_packet,mu_min + (i + rk_double (mt_state)) * mu_bin); + switch (virtual_mode) + { + case -2: + weight = 1.0 / rpacket_get_virtual_packet_flag (packet); + break; + case -1: + weight = + 2.0 * rpacket_get_mu(&virt_packet) / + rpacket_get_virtual_packet_flag (packet); + break; + case 1: + weight = + (1.0 - + mu_min) / 2.0 / rpacket_get_virtual_packet_flag (packet); + break; + default: + fprintf (stderr, "Something has gone horribly wrong!\n"); + // FIXME MR: we need to somehow signal an error here + // I'm adding an exit() here to inform the compiler about the impossible path + exit(1); + } + angle_aberration_CMF_to_LF (&virt_packet, storage); + double doppler_factor_ratio = + rpacket_doppler_factor (packet, storage) / + rpacket_doppler_factor (&virt_packet, storage); + rpacket_set_energy(&virt_packet, + rpacket_get_energy (packet) * doppler_factor_ratio); + rpacket_set_nu(&virt_packet,rpacket_get_nu (packet) * doppler_factor_ratio); + reabsorbed = montecarlo_one_packet_loop (storage, &virt_packet, 1, mt_state); +#ifdef WITH_VPACKET_LOGGING +#ifdef WITHOPENMP +#pragma omp critical + { +#endif // WITHOPENMP + if (storage->virt_packet_count >= storage->virt_array_size) + { + storage->virt_array_size *= 2; + storage->virt_packet_nus = safe_realloc(storage->virt_packet_nus, sizeof(double) * storage->virt_array_size); + storage->virt_packet_energies = safe_realloc(storage->virt_packet_energies, sizeof(double) * storage->virt_array_size); + storage->virt_packet_last_interaction_in_nu = safe_realloc(storage->virt_packet_last_interaction_in_nu, sizeof(double) * storage->virt_array_size); + storage->virt_packet_last_interaction_type = safe_realloc(storage->virt_packet_last_interaction_type, sizeof(int64_t) * storage->virt_array_size); + storage->virt_packet_last_line_interaction_in_id = safe_realloc(storage->virt_packet_last_line_interaction_in_id, sizeof(int64_t) * storage->virt_array_size); + storage->virt_packet_last_line_interaction_out_id = safe_realloc(storage->virt_packet_last_line_interaction_out_id, sizeof(int64_t) * storage->virt_array_size); + } + storage->virt_packet_nus[storage->virt_packet_count] = rpacket_get_nu(&virt_packet); + storage->virt_packet_energies[storage->virt_packet_count] = rpacket_get_energy(&virt_packet) * weight; + storage->virt_packet_last_interaction_in_nu[storage->virt_packet_count] = storage->last_interaction_in_nu[rpacket_get_id (packet)]; + storage->virt_packet_last_interaction_type[storage->virt_packet_count] = storage->last_interaction_type[rpacket_get_id (packet)]; + storage->virt_packet_last_line_interaction_in_id[storage->virt_packet_count] = storage->last_line_interaction_in_id[rpacket_get_id (packet)]; + storage->virt_packet_last_line_interaction_out_id[storage->virt_packet_count] = storage->last_line_interaction_out_id[rpacket_get_id (packet)]; + storage->virt_packet_count += 1; +#ifdef WITHOPENMP + } +#endif // WITHOPENMP +#endif // WITH_VPACKET_LOGGING + if ((rpacket_get_nu(&virt_packet) < storage->spectrum_end_nu) && + (rpacket_get_nu(&virt_packet) > storage->spectrum_start_nu)) + { +#ifdef WITHOPENMP +#pragma omp critical + { +#endif // WITHOPENMP + int64_t virt_id_nu = + floor ((rpacket_get_nu(&virt_packet) - + storage->spectrum_start_nu) / + storage->spectrum_delta_nu); + storage->spectrum_virt_nu[virt_id_nu] += + rpacket_get_energy(&virt_packet) * weight; +#ifdef WITHOPENMP + } +#endif // WITHOPENMP + } + } + } + else + { + return 1; + } + } + return reabsorbed; +} + +void +move_packet_across_shell_boundary (rpacket_t * packet, + storage_model_t * storage, double distance, rk_state *mt_state) +{ + move_packet (packet, storage, distance); + if (rpacket_get_virtual_packet (packet) > 0) + { + double delta_tau_event = rpacket_get_chi_continuum(packet) * distance; + rpacket_set_tau_event (packet, + rpacket_get_tau_event (packet) + + delta_tau_event); + packet->compute_chi_bf = true; + } + else + { + rpacket_reset_tau_event (packet, mt_state); + } + if ((rpacket_get_current_shell_id (packet) < storage->no_of_shells - 1 + && rpacket_get_next_shell_id (packet) == 1) + || (rpacket_get_current_shell_id (packet) > 0 + && rpacket_get_next_shell_id (packet) == -1)) + { + rpacket_set_current_shell_id (packet, + rpacket_get_current_shell_id (packet) + + rpacket_get_next_shell_id (packet)); + } + else if (rpacket_get_next_shell_id (packet) == 1) + { + rpacket_set_status (packet, TARDIS_PACKET_STATUS_EMITTED); + } + else if ((storage->reflective_inner_boundary == 0) || + (rk_double (mt_state) > storage->inner_boundary_albedo)) + { + rpacket_set_status (packet, TARDIS_PACKET_STATUS_REABSORBED); + } + else + { + double doppler_factor = rpacket_doppler_factor (packet, storage); + double comov_nu = rpacket_get_nu (packet) * doppler_factor; + double comov_energy = rpacket_get_energy (packet) * doppler_factor; + // TODO: correct + rpacket_set_mu (packet, rk_double (mt_state)); + double inverse_doppler_factor = rpacket_inverse_doppler_factor (packet, storage); + rpacket_set_nu (packet, comov_nu * inverse_doppler_factor); + rpacket_set_energy (packet, comov_energy * inverse_doppler_factor); + if (rpacket_get_virtual_packet_flag (packet) > 0) + { + montecarlo_one_packet (storage, packet, -2, mt_state); + } + } +} + +void +montecarlo_thomson_scatter (rpacket_t * packet, storage_model_t * storage, + double distance, rk_state *mt_state) +{ + move_packet (packet, storage, distance); + double doppler_factor = rpacket_doppler_factor (packet, storage); + double comov_nu = rpacket_get_nu (packet) * doppler_factor; + double comov_energy = rpacket_get_energy (packet) * doppler_factor; + rpacket_set_mu (packet, 2.0 * rk_double (mt_state) - 1.0); + double inverse_doppler_factor = rpacket_inverse_doppler_factor (packet, storage); + rpacket_set_nu (packet, comov_nu * inverse_doppler_factor); + rpacket_set_energy (packet, comov_energy * inverse_doppler_factor); + rpacket_reset_tau_event (packet, mt_state); + storage->last_interaction_type[rpacket_get_id (packet)] = 1; + + angle_aberration_CMF_to_LF (packet, storage); + + if (rpacket_get_virtual_packet_flag (packet) > 0) + { + create_vpacket (storage, packet, mt_state); + } +} + +void +montecarlo_bound_free_scatter (rpacket_t * packet, storage_model_t * storage, double distance, rk_state *mt_state) +{ + // current position in list of continuum edges -> indicates which bound-free processes are possible + int64_t ccontinuum = rpacket_get_current_continuum_id (packet); + + // Determine in which continuum the bf-absorption occurs + double chi_bf = rpacket_get_chi_boundfree (packet); + double zrand = rk_double (mt_state); + double zrand_x_chibf = zrand * chi_bf; + + while ((ccontinuum < storage->no_of_edges - 1) && (packet->chi_bf_tmp_partial[ccontinuum] <= zrand_x_chibf)) + { + ccontinuum++; + } + rpacket_set_current_continuum_id (packet, ccontinuum); + + /* For consistency reasons the branching between ionization and thermal energy is determined using the + comoving frequency at the initial position instead of the frequency at the point of interaction */ + double comov_nu = rpacket_get_nu (packet) * rpacket_doppler_factor (packet, storage); + + /* Move the packet to the place of absorption, select a direction for re-emission and impose energy conservation + in the co-moving frame. */ + move_packet (packet, storage, distance); + double old_doppler_factor = rpacket_doppler_factor (packet, storage); + rpacket_set_mu (packet, 2.0 * rk_double (mt_state) - 1.0); + double inverse_doppler_factor = rpacket_inverse_doppler_factor (packet, storage); + double comov_energy = rpacket_get_energy (packet) * old_doppler_factor; + rpacket_set_energy (packet, comov_energy * inverse_doppler_factor); + storage->last_interaction_type[rpacket_get_id (packet)] = 3; // last interaction was a bf-absorption + storage->last_line_interaction_in_id[rpacket_get_id (packet)] = ccontinuum; + + // Convert the rpacket to thermal or ionization energy + zrand = rk_double (mt_state); + int64_t activate_level = (zrand < storage->continuum_list_nu[ccontinuum] / comov_nu) ? + storage->cont_edge2macro_level[ccontinuum] : storage->kpacket2macro_level; + + rpacket_set_macro_atom_activation_level (packet, activate_level); + macro_atom (packet, storage, mt_state); +} + +void +montecarlo_free_free_scatter (rpacket_t * packet, storage_model_t * storage, double distance, rk_state *mt_state) +{ + /* Move the packet to the place of absorption, select a direction for re-emission and impose energy conservation + in the co-moving frame. */ + move_packet (packet, storage, distance); + double old_doppler_factor = rpacket_doppler_factor (packet, storage); + rpacket_set_mu (packet, 2.0 * rk_double (mt_state) - 1.0); + double inverse_doppler_factor = rpacket_inverse_doppler_factor (packet, storage); + double comov_energy = rpacket_get_energy (packet) * old_doppler_factor; + rpacket_set_energy (packet, comov_energy * inverse_doppler_factor); + storage->last_interaction_type[rpacket_get_id (packet)] = 4; // last interaction was a ff-absorption + + // Create a k-packet + rpacket_set_macro_atom_activation_level (packet, storage->kpacket2macro_level); + macro_atom (packet, storage, mt_state); +} + +double +sample_nu_free_free (const rpacket_t * packet, const storage_model_t * storage, rk_state *mt_state) +{ + int64_t shell_id = rpacket_get_current_shell_id (packet); + double T = storage->t_electrons[shell_id]; + double zrand = rk_double (mt_state); + return -KB * T / H * log(zrand); // Lucy 2003 MC II Eq.41 +} + +double +sample_nu_free_bound (const rpacket_t * packet, const storage_model_t * storage, rk_state *mt_state) +{ + int64_t continuum_id = rpacket_get_current_continuum_id (packet); + double th_frequency = storage->continuum_list_nu[continuum_id]; + int64_t shell_id = rpacket_get_current_shell_id (packet); + double T = storage->t_electrons[shell_id]; + double zrand = rk_double (mt_state); + return th_frequency * (1 - (KB * T / H / th_frequency * log(zrand))); // Lucy 2003 MC II Eq.26 +} + +void +montecarlo_line_scatter (rpacket_t * packet, storage_model_t * storage, + double distance, rk_state *mt_state) +{ + uint64_t next_line_id = rpacket_get_next_line_id (packet); + uint64_t line2d_idx = next_line_id + + storage->no_of_lines * rpacket_get_current_shell_id (packet); + if (rpacket_get_virtual_packet (packet) == 0) + { + increment_j_blue_estimator (packet, storage, distance, line2d_idx); + increment_Edotlu_estimator (packet, storage, distance, line2d_idx); + } + double tau_line = + storage->line_lists_tau_sobolevs[line2d_idx]; + double tau_continuum = rpacket_get_chi_continuum(packet) * distance; + double tau_combined = tau_line + tau_continuum; + //rpacket_set_next_line_id (packet, rpacket_get_next_line_id (packet) + 1); + + if (next_line_id + 1 == storage->no_of_lines) + { + rpacket_set_last_line (packet, true); + } + if (rpacket_get_virtual_packet (packet) > 0) + { + rpacket_set_tau_event (packet, + rpacket_get_tau_event (packet) + tau_line); + rpacket_set_next_line_id (packet, next_line_id + 1); + test_for_close_line (packet, storage); + } + else if (rpacket_get_tau_event (packet) < tau_combined) + { // Line absorption occurs + move_packet (packet, storage, distance); + double old_doppler_factor = rpacket_doppler_factor (packet, storage); + rpacket_set_mu (packet, 2.0 * rk_double (mt_state) - 1.0); + double inverse_doppler_factor = rpacket_inverse_doppler_factor (packet, storage); + double comov_energy = rpacket_get_energy (packet) * old_doppler_factor; + rpacket_set_energy (packet, comov_energy * inverse_doppler_factor); + storage->last_interaction_in_nu[rpacket_get_id (packet)] = + rpacket_get_nu (packet); + storage->last_line_interaction_in_id[rpacket_get_id (packet)] = + next_line_id; + storage->last_line_interaction_shell_id[rpacket_get_id (packet)] = + rpacket_get_current_shell_id (packet); + storage->last_interaction_type[rpacket_get_id (packet)] = 2; + if (storage->line_interaction_id == 0) + { + line_emission (packet, storage, next_line_id, mt_state); + } + else if (storage->line_interaction_id >= 1) + { + rpacket_set_macro_atom_activation_level (packet, + storage->line2macro_level_upper[next_line_id]); + macro_atom (packet, storage, mt_state); + } + } + else + { // Packet passes line without interacting + rpacket_set_tau_event (packet, + rpacket_get_tau_event (packet) - tau_line); + rpacket_set_next_line_id (packet, next_line_id + 1); + packet->compute_chi_bf = false; + test_for_close_line (packet, storage); + } +} + +void +line_emission (rpacket_t * packet, storage_model_t * storage, int64_t emission_line_id, rk_state *mt_state) +{ + double inverse_doppler_factor = rpacket_inverse_doppler_factor (packet, storage); + storage->last_line_interaction_out_id[rpacket_get_id (packet)] = emission_line_id; + if (storage->cont_status == CONTINUUM_ON) + { + storage->last_interaction_out_type[rpacket_get_id (packet)] = 2; + } + + rpacket_set_nu (packet, + storage->line_list_nu[emission_line_id] * inverse_doppler_factor); + rpacket_set_nu_line (packet, storage->line_list_nu[emission_line_id]); + rpacket_set_next_line_id (packet, emission_line_id + 1); + rpacket_reset_tau_event (packet, mt_state); + + angle_aberration_CMF_to_LF (packet, storage); + + if (rpacket_get_virtual_packet_flag (packet) > 0) + { + bool virtual_close_line = false; + if (!rpacket_get_last_line (packet) && + fabs (storage->line_list_nu[rpacket_get_next_line_id (packet)] - + rpacket_get_nu_line (packet)) < + (rpacket_get_nu_line (packet)* 1e-7)) + { + virtual_close_line = true; + } + // QUESTIONABLE!!! + bool old_close_line = rpacket_get_close_line (packet); + rpacket_set_close_line (packet, virtual_close_line); + create_vpacket (storage, packet, mt_state); + rpacket_set_close_line (packet, old_close_line); + virtual_close_line = false; + } + test_for_close_line (packet, storage); +} + +void test_for_close_line (rpacket_t * packet, const storage_model_t * storage) +{ + if (!rpacket_get_last_line (packet) && + fabs (storage->line_list_nu[rpacket_get_next_line_id (packet)] - + rpacket_get_nu_line (packet)) < (rpacket_get_nu_line (packet)* + 1e-7)) + { + rpacket_set_close_line (packet, true); + } +} + +void +continuum_emission (rpacket_t * packet, storage_model_t * storage, rk_state *mt_state, + pt2sample_nu sample_nu_continuum, int64_t emission_type_id) +{ + double inverse_doppler_factor = rpacket_inverse_doppler_factor (packet, storage); + double nu_comov = sample_nu_continuum (packet, storage, mt_state); + rpacket_set_nu (packet, nu_comov * inverse_doppler_factor); + rpacket_reset_tau_event (packet, mt_state); + + storage->last_interaction_out_type[rpacket_get_id (packet)] = emission_type_id; + + // Have to find current position in line list + int64_t current_line_id; + line_search (storage->line_list_nu, nu_comov, storage->no_of_lines, ¤t_line_id); + + bool last_line = (current_line_id == storage->no_of_lines); + rpacket_set_last_line (packet, last_line); + rpacket_set_next_line_id (packet, current_line_id); + + angle_aberration_CMF_to_LF (packet, storage); + + if (rpacket_get_virtual_packet_flag (packet) > 0) + { + create_vpacket (storage, packet, mt_state); + } +} + +static void +montecarlo_compute_distances (rpacket_t * packet, storage_model_t * storage) +{ + // Check if the last line was the same nu as the current line. + if (rpacket_get_close_line (packet)) + { + // If so set the distance to the line to 0.0 + rpacket_set_d_line (packet, 0.0); + // Reset close_line. + rpacket_set_close_line (packet, false); + } + else + { + compute_distance2boundary (packet, storage); + compute_distance2line (packet, storage); + // FIXME MR: return status of compute_distance2line() is ignored + compute_distance2continuum (packet, storage); + } +} + +montecarlo_event_handler_t +get_event_handler (rpacket_t * packet, storage_model_t * storage, + double *distance, rk_state *mt_state) +{ + montecarlo_compute_distances (packet, storage); + double d_boundary = rpacket_get_d_boundary (packet); + double d_continuum = rpacket_get_d_continuum (packet); + double d_line = rpacket_get_d_line (packet); + montecarlo_event_handler_t handler; + if (d_line <= d_boundary && d_line <= d_continuum) + { + *distance = d_line; + handler = &montecarlo_line_scatter; + } + else if (d_boundary <= d_continuum) + { + *distance = d_boundary; + handler = &move_packet_across_shell_boundary; + } + else + { + *distance = d_continuum; + handler = montecarlo_continuum_event_handler (packet, storage, mt_state); + } + return handler; +} + +montecarlo_event_handler_t +montecarlo_continuum_event_handler (rpacket_t * packet, storage_model_t * storage, rk_state *mt_state) +{ + if (storage->cont_status) + { + double zrand_x_chi_cont = rk_double (mt_state) * rpacket_get_chi_continuum (packet); + double chi_th = rpacket_get_chi_electron (packet); + double chi_bf = rpacket_get_chi_boundfree (packet); + + if (zrand_x_chi_cont < chi_th) + { + return &montecarlo_thomson_scatter; + } + else if (zrand_x_chi_cont < chi_th + chi_bf) + { + return &montecarlo_bound_free_scatter; + } + else + { + return &montecarlo_free_free_scatter; + } + } + else + { + return &montecarlo_thomson_scatter; + } +} + +int64_t +montecarlo_one_packet_loop (storage_model_t * storage, rpacket_t * packet, + int64_t virtual_packet, rk_state *mt_state) +{ + rpacket_set_tau_event (packet, 0.0); + rpacket_set_nu_line (packet, 0.0); + rpacket_set_virtual_packet (packet, virtual_packet); + rpacket_set_status (packet, TARDIS_PACKET_STATUS_IN_PROCESS); + // Initializing tau_event if it's a real packet. + if (virtual_packet == 0) + { + rpacket_reset_tau_event (packet,mt_state); + } + // For a virtual packet tau_event is the sum of all the tau's that the packet passes. + while (rpacket_get_status (packet) == TARDIS_PACKET_STATUS_IN_PROCESS) + { + // Check if we are at the end of line list. + if (!rpacket_get_last_line (packet)) + { + rpacket_set_nu_line (packet, + storage-> + line_list_nu[rpacket_get_next_line_id + (packet)]); + } + double distance; + get_event_handler (packet, storage, &distance, mt_state) (packet, storage, + distance, mt_state); + if (virtual_packet > 0 && rpacket_get_tau_event (packet) > storage->tau_russian) + { + double event_random = rk_double (mt_state); + if (event_random > storage->survival_probability) + { + rpacket_set_energy(packet, 0.0); + rpacket_set_status (packet, TARDIS_PACKET_STATUS_EMITTED); + } + else + { + rpacket_set_energy(packet, + rpacket_get_energy (packet) / storage->survival_probability * + exp (-1.0 * rpacket_get_tau_event (packet))); + rpacket_set_tau_event (packet, 0.0); + } + } + } + if (virtual_packet > 0) + { + rpacket_set_energy (packet, + rpacket_get_energy (packet) * exp (-1.0 * + rpacket_get_tau_event + (packet))); + } + return rpacket_get_status (packet) == + TARDIS_PACKET_STATUS_REABSORBED ? 1 : 0; +} + +void +montecarlo_main_loop(storage_model_t * storage, int64_t virtual_packet_flag, int nthreads, unsigned long seed) +{ + int64_t finished_packets = 0; + storage->virt_packet_count = 0; +#ifdef WITH_VPACKET_LOGGING + storage->virt_packet_nus = (double *)safe_malloc(sizeof(double) * storage->no_of_packets); + storage->virt_packet_energies = (double *)safe_malloc(sizeof(double) * storage->no_of_packets); + storage->virt_packet_last_interaction_in_nu = (double *)safe_malloc(sizeof(double) * storage->no_of_packets); + storage->virt_packet_last_interaction_type = (int64_t *)safe_malloc(sizeof(int64_t) * storage->no_of_packets); + storage->virt_packet_last_line_interaction_in_id = (int64_t *)safe_malloc(sizeof(int64_t) * storage->no_of_packets); + storage->virt_packet_last_line_interaction_out_id = (int64_t *)safe_malloc(sizeof(int64_t) * storage->no_of_packets); + storage->virt_array_size = storage->no_of_packets; +#endif // WITH_VPACKET_LOGGING +#ifdef WITHOPENMP + omp_set_dynamic(0); + if (nthreads > 0) + { + omp_set_num_threads(nthreads); + } + +#pragma omp parallel firstprivate(finished_packets) + { + rk_state mt_state; + rk_seed (seed + omp_get_thread_num(), &mt_state); +#pragma omp master + { + fprintf(stderr, "Running with OpenMP - %d threads\n", omp_get_num_threads()); + print_progress(0, storage->no_of_packets); + } +#else + rk_state mt_state; + rk_seed (seed, &mt_state); + fprintf(stderr, "Running without OpenMP\n"); +#endif + int64_t chi_bf_tmp_size = (storage->cont_status) ? storage->no_of_edges : 0; + double *chi_bf_tmp_partial = safe_malloc(sizeof(double) * chi_bf_tmp_size); + + #pragma omp for + for (int64_t packet_index = 0; packet_index < storage->no_of_packets; ++packet_index) + { + int reabsorbed = 0; + rpacket_t packet; + rpacket_set_id(&packet, packet_index); + rpacket_init(&packet, storage, packet_index, virtual_packet_flag, chi_bf_tmp_partial); + if (virtual_packet_flag > 0) + { + reabsorbed = montecarlo_one_packet(storage, &packet, -1, &mt_state); + } + reabsorbed = montecarlo_one_packet(storage, &packet, 0, &mt_state); + storage->output_nus[packet_index] = rpacket_get_nu(&packet); + if (reabsorbed == 1) + { + storage->output_energies[packet_index] = -rpacket_get_energy(&packet); + } + else + { + storage->output_energies[packet_index] = rpacket_get_energy(&packet); + } + if ( ++finished_packets%100 == 0 ) + { +#ifdef WITHOPENMP + // WARNING: This only works with a static sheduler and gives an approximation of progress. + // The alternative would be to have a shared variable but that could potentially decrease performance when using many threads. + if (omp_get_thread_num() == 0 ) + print_progress(finished_packets * omp_get_num_threads(), storage->no_of_packets); +#else + print_progress(finished_packets, storage->no_of_packets); +#endif + } + } + free(chi_bf_tmp_partial); +#ifdef WITHOPENMP + } +#endif + 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); + } +} diff --git a/tardis/montecarlo/src/cmontecarlo.h b/tardis/montecarlo/src/cmontecarlo.h new file mode 100644 index 00000000000..767017b7303 --- /dev/null +++ b/tardis/montecarlo/src/cmontecarlo.h @@ -0,0 +1,159 @@ +#ifndef TARDIS_CMONTECARLO_H +#define TARDIS_CMONTECARLO_H + +#include +#include "randomkit/randomkit.h" +#include "rpacket.h" +#include "status.h" +#include "storage.h" + +#ifdef WITH_VPACKET_LOGGING +#define LOG_VPACKETS 1 +#else +#define LOG_VPACKETS 0 +#endif + +typedef void (*montecarlo_event_handler_t) (rpacket_t *packet, + storage_model_t *storage, + double distance, rk_state *mt_state); + +typedef double (*pt2sample_nu) (const rpacket_t *packet, + const storage_model_t *storage, + rk_state *mt_state); + +void initialize_random_kit (unsigned long seed); + +tardis_error_t line_search (const double *nu, double nu_insert, + int64_t number_of_lines, int64_t * result); + +tardis_error_t +reverse_binary_search (const double *x, double x_insert, + int64_t imin, int64_t imax, int64_t * result); + +tardis_error_t +binary_search (const double *x, double x_insert, + int64_t imin, int64_t imax, int64_t * result); + +double rpacket_doppler_factor (const rpacket_t *packet, const storage_model_t *storage); + +double rpacket_inverse_doppler_factor (const rpacket_t *packet, const storage_model_t *storage); + +void +angle_aberration_CMF_to_LF (rpacket_t * packet, const storage_model_t * storage); + +double +angle_aberration_LF_to_CMF (rpacket_t *packet, const storage_model_t *storage, double mu); + +/** Calculate the distance to shell boundary. + * + * @param packet rpacket structure with packet information + * @param storage storage model data + * + * @return distance to shell boundary + */ +void compute_distance2boundary (rpacket_t * packet, + const storage_model_t * storage); + +/** Calculate the distance the packet has to travel until it redshifts to the first spectral line. + * + * @param packet rpacket structure with packet information + * @param storage storage model data + * + * @return distance to the next spectral line + */ +tardis_error_t compute_distance2line (rpacket_t * packet, + const storage_model_t * storage); + +/** Calculate the distance to the next continuum event, which can be a Thomson scattering, bound-free absorption or + free-free transition. + * + * @param packet rpacket structure with packet information + * @param storage storage model data + * + * sets distance to the next continuum event (in centimeters) in packet rpacket structure + */ +void compute_distance2continuum (rpacket_t *packet, storage_model_t *storage); + +void macro_atom (rpacket_t * packet, const storage_model_t * storage, rk_state *mt_state); + +void move_packet (rpacket_t * packet, storage_model_t * storage, + double distance); + +void increment_j_blue_estimator (const rpacket_t * packet, + storage_model_t * storage, + double d_line, + int64_t j_blue_idx); + +void increment_Edotlu_estimator (const rpacket_t * packet, + storage_model_t * storage, + double d_line, + int64_t j_blue_idx); + +double get_increment_j_blue_estimator_energy (const rpacket_t * packet, + const storage_model_t * storage, + double d_line); + +void +increment_continuum_estimators (const rpacket_t * packet, storage_model_t * storage, double distance, + double comov_nu, double comov_energy); + +int64_t montecarlo_one_packet (storage_model_t * storage, rpacket_t * packet, + int64_t virtual_mode, rk_state *mt_state); + +int64_t montecarlo_one_packet_loop (storage_model_t * storage, + rpacket_t * packet, + int64_t virtual_packet, rk_state *mt_state); + +void montecarlo_main_loop(storage_model_t * storage, + int64_t virtual_packet_flag, + int nthreads, + unsigned long seed); + +montecarlo_event_handler_t get_event_handler (rpacket_t * packet, storage_model_t * storage, + double *distance, rk_state *mt_state); + +void test_for_close_line(rpacket_t * packet, const storage_model_t * storage); + +/* New handlers for continuum implementation */ + +montecarlo_event_handler_t montecarlo_continuum_event_handler(rpacket_t * packet, storage_model_t * storage, rk_state *mt_state); + +void montecarlo_free_free_scatter (rpacket_t * packet, storage_model_t * storage, double distance, rk_state *mt_state); + +void montecarlo_bound_free_scatter (rpacket_t * packet, storage_model_t * storage, double distance, rk_state *mt_state); + +double +bf_cross_section(const storage_model_t * storage, int64_t continuum_id, double comov_nu); + +void calculate_chi_bf(rpacket_t * packet, storage_model_t * storage); + +void calculate_chi_ff(rpacket_t * packet, const storage_model_t * storage); + +void +move_packet_across_shell_boundary (rpacket_t * packet, + storage_model_t * storage, double distance, rk_state *mt_state); + +void +montecarlo_thomson_scatter (rpacket_t * packet, storage_model_t * storage, + double distance, rk_state *mt_state); + +void +montecarlo_line_scatter (rpacket_t * packet, storage_model_t * storage, + double distance, rk_state *mt_state); + +void +line_emission(rpacket_t * packet, storage_model_t * storage, + int64_t emission_line_id, rk_state *mt_state); + +double sample_nu_free_free(const rpacket_t * packet, const storage_model_t * storage, rk_state *mt_state); + +double sample_nu_free_bound(const rpacket_t * packet, const storage_model_t * storage, rk_state *mt_state); + +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 diff --git a/tardis/montecarlo/src/integrator.c b/tardis/montecarlo/src/integrator.c new file mode 100644 index 00000000000..a0e1452a00e --- /dev/null +++ b/tardis/montecarlo/src/integrator.c @@ -0,0 +1,344 @@ +#define _USE_MATH_DEFINES + +#include +#include +#include +#include +#include +#include + +#include "io.h" +#include "storage.h" +#include "integrator.h" +#include "cmontecarlo.h" + +#include "omp_helper.h" + +#define NULEN 0 +#define LINELEN 1 +#define PLEN 2 +#define SHELLEN 3 + +#define C_INV 3.33564e-11 +#define M_PI acos (-1) +#define KB_CGS 1.3806488e-16 +#define H_CGS 6.62606957e-27 + +/** + * Calculate the intensity of a black-body according to the following formula + * .. math:: + * I(\\nu, T) = \\frac{2h\\nu^3}{c^2}\frac{1}{e^{h\\nu \\beta_\\textrm{rad}} - 1} +*/ +double +intensity_black_body (double nu, double T) +{ + double beta_rad = 1 / (KB_CGS * T); + double coefficient = 2 * H_CGS * C_INV * C_INV; + return coefficient * nu * nu * nu / (exp(H_CGS * nu * beta_rad) - 1 ); +} + + +/*! @brief Algorithm to integrate an array using the trapezoid integration rule + * +*/ +double +trapezoid_integration (const double* array, const double h, int N) +{ + double result = (array[0] + array[N-1])/2; + for (int idx = 1; idx < N-1; ++idx) + { + result += array[idx]; + } + return result * h; +} + +/*! @brief Calculate distance to p line + * + * Calculate half of the length of the p-line inside a shell + * of radius r in terms of unit length (c * t_exp). + * If shell and p-line do not intersect, return 0. + * + * @param r radius of the shell + * @param p distance of the p-line to the center of the supernova + * @param inv_t inverse time_explosio is needed to norm to unit-length + * @return half the lenght inside the shell or zero + */ +static inline double +calculate_z(double r, double p, double inv_t) +{ + return (r > p) ? sqrt(r * r - p * p) * C_INV * inv_t : 0; +} + + +/*! + * @brief Calculate p line intersections + * + * This function calculates the intersection points of the p-line with each shell + * + * @param storage (INPUT) A storage model containing the environment + * @param p (INPUT) distance of the integration line to the center + * @param oz (OUTPUT) will be set with z values. The array is truncated by the + * value `1`. + * @param oshell_id (OUTPUT) will be set with the corresponding shell_ids + * @return number of shells intersected by the p-line + */ +int64_t +populate_z(const storage_model_t *storage, const double p, double *oz, int64_t *oshell_id) +{ + + // Abbreviations + double *r = storage->r_outer_i; + const int64_t N = storage->no_of_shells_i; + double inv_t = storage->inverse_time_explosion; + double z = 0; + + int64_t i = 0, offset = N, i_low, i_up; + + if (p <= storage->r_inner_i[0]) + { + // Intersect the photosphere + for(i = 0; i < N; ++i) + { // Loop from inside to outside + oz[i] = 1 - calculate_z(r[i], p, inv_t); + oshell_id[i] = i; + } + return N; + } + else + { + // No intersection with the photosphere + // that means we intersect each shell twice + for(i = 0; i < N; ++i) + { // Loop from inside to outside + z = calculate_z(r[i], p, inv_t); + if (z == 0) + continue; + if (offset == N) + { + offset = i; + } + // Calculate the index in the resulting array + i_low = N - i - 1; // the far intersection with the shell + i_up = N + i - 2 * offset; // the nearer intersection with the shell + + // Setting the arrays + oz[i_low] = 1 + z; + oshell_id[i_low] = i; + oz[i_up] = 1 - z; + oshell_id[i_up] = i; + } + return 2 * (N - offset); + } +} + + +/*! @brief Calculate integration points + * + */ +void +calculate_p_values(double R_max, int64_t N, double *opp) +{ + for(int i = 0; ino_of_lines, + size_shell = storage->no_of_shells_i, + size_tau = size_line * size_shell, + finished_nus = 0; + + double R_ph = storage->r_inner_i[0]; + double R_max = storage->r_outer_i[size_shell - 1]; + double pp[N]; + double *exp_tau = calloc(size_tau, sizeof(double)); +//#pragma omp parallel firstprivate(L, exp_tau) +#pragma omp parallel + { + +#pragma omp master + { + if (omp_get_num_threads() > 1) { + fprintf(stderr, "Doing the formal integral\nRunning with OpenMP - %d threads\n", omp_get_num_threads()); + } else { + fprintf(stderr, "Doing the formal integral\nRunning without OpenMP\n"); + } + print_progress_fi(0, inu_size); + } + + // Initializing all the thread-local variables + int64_t offset = 0, i = 0, + size_z = 0, + idx_nu_start = 0, + direction = 0, + first = 0; + + double I_nu[N], + //I_nu_b[N], + //I_nu_r[N], + z[2 * storage->no_of_shells_i], + p = 0, + nu_start, + nu_end, + nu, + zstart, + zend, + escat_contrib, + escat_op, + Jkkp; + int64_t shell_id[2 * storage->no_of_shells_i]; + + double *pexp_tau, *patt_S_ul, *pline, *pJred_lu, *pJblue_lu; + + // Prepare exp_tau +#pragma omp for + for (i = 0; i < size_tau; ++i) { + exp_tau[i] = exp( -storage->line_lists_tau_sobolevs_i[i]); + } + calculate_p_values(R_max, N, pp); + // Done with the initialization + + // Loop over wavelengths in spectrum +#pragma omp for + for (int nu_idx = 0; nu_idx < inu_size ; ++nu_idx) + { + nu = inu[nu_idx]; + + // Loop over discrete values along line + for (int p_idx = 1; p_idx < N; ++p_idx) + { + escat_contrib = 0; + p = pp[p_idx]; + + // initialize z intersections for p values + size_z = populate_z(storage, p, z, shell_id); + + // initialize I_nu + if (p <= R_ph) + I_nu[p_idx] = intensity_black_body(nu * z[0], iT); + else + I_nu[p_idx] = 0; + + // Find first contributing line + nu_start = nu * z[0]; + nu_end = nu * z[1]; + line_search( + storage->line_list_nu, + nu_start, + size_line, + &idx_nu_start + ); + offset = shell_id[0] * size_line; + + // start tracking accumulated e-scattering optical depth + zstart = storage->time_explosion / C_INV * (1. - z[0]); + + // Initialize pointers + pline = storage->line_list_nu + idx_nu_start; + pexp_tau = exp_tau + offset + idx_nu_start; + patt_S_ul = att_S_ul + offset + idx_nu_start; + pJred_lu = Jred_lu + offset + idx_nu_start; + pJblue_lu = Jblue_lu + offset + idx_nu_start; + + // flag for first contribution to integration on current p-ray + first = 1; + + // TODO: Ugly loop + // Loop over all intersections + // TODO: replace by number of intersections and remove break + for (i = 0; i < size_z - 1; ++i) + { + escat_op = storage->electron_densities_i[shell_id[i]] * storage->sigma_thomson; + nu_end = nu * z[i+1]; + + // TODO: e-scattering: in principle we also have to check + // that dtau is <<1 (as assumed in Lucy 1999); if not, there + // is the chance that I_nu_b becomes negative + for (;pline < storage->line_list_nu + size_line; + // We have to increment all pointers simultaneously + ++pline, + ++pexp_tau, + ++patt_S_ul, + ++pJblue_lu) + { + if (*pline < nu_end) + { + // next resonance not in current shell + break; + } + + // Calculate e-scattering optical depth to next resonance point + zend = storage->time_explosion / C_INV * (1. - *pline / nu); + + if (first == 1){ + // First contribution to integration + // NOTE: this treatment of I_nu_b (given by boundary + // conditions) is not in Lucy 1999; should be + // re-examined carefully + escat_contrib += (zend - zstart) * escat_op * (*pJblue_lu - I_nu[p_idx]) ; + first = 0; + } + else{ + // Account for e-scattering, c.f. Eqs 27, 28 in Lucy 1999 + Jkkp = 0.5 * (*pJred_lu + *pJblue_lu); + escat_contrib += (zend - zstart) * escat_op * (Jkkp - I_nu[p_idx]) ; + // this introduces the necessary offset of one element between pJblue_lu and pJred_lu + pJred_lu += 1; + } + I_nu[p_idx] = I_nu[p_idx] + escat_contrib; + + // Lucy 1999, Eq 26 + I_nu[p_idx] = I_nu[p_idx] * (*pexp_tau) + *patt_S_ul; + + // reset e-scattering opacity + escat_contrib = 0; + zstart = zend; + } + // Calculate e-scattering optical depth to grid cell boundary + Jkkp = 0.5 * (*pJred_lu + *pJblue_lu); + zend = storage->time_explosion / C_INV * (1. - nu_end / nu); + escat_contrib += (zend - zstart) * escat_op * (Jkkp - I_nu[p_idx]); + zstart = zend; + + if (i < size_z-1){ + // advance pointers + direction = shell_id[i+1] - shell_id[i]; + pexp_tau += direction * size_line; + patt_S_ul += direction * size_line; + pJred_lu += direction * size_line; + pJblue_lu += direction * size_line; + } + } + I_nu[p_idx] *= p; + } + // TODO: change integration to match the calculation of p values + L[nu_idx] = 8 * M_PI * M_PI * trapezoid_integration(I_nu, R_max/N, N); +#pragma omp atomic update + ++finished_nus; + if (finished_nus%10 == 0){ + print_progress_fi(finished_nus, inu_size); + } + } + // Free everything allocated on heap + printf("\n"); + } + free(exp_tau); + return L; +} diff --git a/tardis/montecarlo/src/integrator.h b/tardis/montecarlo/src/integrator.h new file mode 100644 index 00000000000..c06e4a286a0 --- /dev/null +++ b/tardis/montecarlo/src/integrator.h @@ -0,0 +1,19 @@ +#ifndef TARDIS_INTEGRATOR_H +#define TARDIS_INTEGRATOR_H + +#include "storage.h" + +double +integrate_intensity(const double* I_nu, const double h, int N); + +int64_t +populate_z(const storage_model_t *storage, const double p, double *oz, int64_t *oshell_id); + +void +calculate_p_values(double R_max, int64_t N, double *opp); + +double +*_formal_integral( + const storage_model_t *storage, double T, double *nu, int64_t nu_size, double *att_S_ul, double *Jred_lu, double *Jblue_lu, int N); + +#endif diff --git a/tardis/montecarlo/src/io.h b/tardis/montecarlo/src/io.h new file mode 100644 index 00000000000..122af1e5967 --- /dev/null +++ b/tardis/montecarlo/src/io.h @@ -0,0 +1,32 @@ +#define _POSIX_C_SOURCE 1 + +#include +#include +#include + +#define STATUS_FORMAT "\r\033[2K\t[%" PRId64 "%%] Packets(finished/total): %" PRId64 "/%" PRId64 +#define STATUS_FORMAT_FI "\r\033[2K\t[%" PRId64 "%%] Bins(finished/total): %" PRId64 "/%" PRId64 + +static inline void +print_progress (const int64_t current, const int64_t total) +{ + if (isatty(fileno(stderr))) + { + fprintf(stderr, STATUS_FORMAT, + current * 100 / total, + current, + total); + } +} + +static inline void +print_progress_fi (const int64_t current, const int64_t total) +{ + if (isatty(fileno(stderr))) + { + fprintf(stderr, STATUS_FORMAT_FI, + current * 100 / total, + current, + total); + } +} diff --git a/tardis/montecarlo/src/omp_helper.h b/tardis/montecarlo/src/omp_helper.h new file mode 100644 index 00000000000..abb9c2acc7c --- /dev/null +++ b/tardis/montecarlo/src/omp_helper.h @@ -0,0 +1,7 @@ +#ifdef WITHOPENMP + #include +#else +int omp_get_num_threads(){ + return 1; +} +#endif diff --git a/tardis/montecarlo/src/randomkit/LICENSE b/tardis/montecarlo/src/randomkit/LICENSE new file mode 100644 index 00000000000..c0d13eb1acc --- /dev/null +++ b/tardis/montecarlo/src/randomkit/LICENSE @@ -0,0 +1,20 @@ +Copyright (c) 2003-2006, Jean-Sebastien Roy (js@jeannot.org) + +Permission is hereby granted, free of charge, to any person obtaining a +copy of this software and associated documentation files (the +"Software"), to deal in the Software without restriction, including +without limitation the rights to use, copy, modify, merge, publish, +distribute, sublicense, and/or sell copies of the Software, and to +permit persons to whom the Software is furnished to do so, subject to +the following conditions: + +The above copyright notice and this permission notice shall be included +in all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. +IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY +CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, +TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE +SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. diff --git a/tardis/montecarlo/src/randomkit/randomkit.h b/tardis/montecarlo/src/randomkit/randomkit.h new file mode 100644 index 00000000000..d5220d01857 --- /dev/null +++ b/tardis/montecarlo/src/randomkit/randomkit.h @@ -0,0 +1,42 @@ +/* Random kit 1.6 */ + +/* + * Anyone who attempts to generate random numbers by deterministic means is, + * of course, living in a state of sin. + * -- John von Neumann + */ + +/* + * Copyright (c) 2003-2006, Jean-Sebastien Roy (js@jeannot.org) + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. + * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY + * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, + * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE + * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +/* @(#) $Jeannot: randomkit.h,v 1.29 2006/02/19 14:44:14 js Exp $ */ + +#ifndef _RANDOMKIT_ +#define _RANDOMKIT_ + +#include "rk_mt.h" +#include "rk_sobol.h" +#include "rk_primitive.h" +#include "rk_isaac.h" + +#endif /* _RANDOMKIT_ */ diff --git a/tardis/montecarlo/src/randomkit/rk_isaac.c b/tardis/montecarlo/src/randomkit/rk_isaac.c new file mode 100644 index 00000000000..9f52fafae4c --- /dev/null +++ b/tardis/montecarlo/src/randomkit/rk_isaac.c @@ -0,0 +1,258 @@ +/* Random kit 1.6 */ + +/* + * Copyright (c) 2004-2006, Jean-Sebastien Roy (js@jeannot.org) + * + * ISAAC RNG by By Bob Jenkins. Based on Bob Jenkins public domain code. + * + * Original algorithm for the implementation of rk_interval function from + * Richard J. Wagner's implementation of the Mersenne Twister RNG, optimised by + * Magnus Jonsson. + * + * Constants used in the rk_double implementation by Isaku Wada. + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. + * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY + * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, + * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE + * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +static char const rcsid[] = + "@(#) $Jeannot: rk_isaac.c,v 1.4 2006/02/20 19:13:37 js Exp $"; + +#include "rk_isaac.h" +#include +#include +#include + +/* use the contents of randrsl[0..RK_ISAAC_STATE_LEN-1] as the seed. */ +static void isaac_init(rk_isaac_state *rg); + +void rk_isaac_seed(unsigned long seed, rk_isaac_state *state) +{ + rk_knuth_fill(seed, state->randrsl, RK_ISAAC_STATE_LEN); + isaac_init(state); +} + +rk_error rk_isaac_randomseed(rk_isaac_state *state) +{ + if(rk_devfill(state->randrsl, sizeof(state->randrsl), 1) == RK_NOERR) + { + isaac_init(state); + return RK_NOERR; + } + + rk_isaac_seed(rk_seedfromsystem(), state); + + return RK_ENODEV; +} + +#define ind(mm,x) ((mm)[(x>>2)&(RK_ISAAC_STATE_LEN-1)]) +#define rngstep(mix,a,b,mm,m,m2,r,x) \ +{ \ + x = *m; \ + a = ((a^(mix)) + *(m2++)) & 0xffffffff; \ + *(m++) = y = (ind(mm,x) + a + b) & 0xffffffff; \ + *(r++) = b = (ind(mm,y>>RK_ISAAC_STATE_POW) + x) & 0xffffffff; \ +} + +/* Call rand(rk_isaac_state *r) to retrieve a single 32-bit random value */ +unsigned long rk_isaac_random(rk_isaac_state *state) +{ + if (!state->randcnt--) + { + register unsigned long a,b,x,y,*m,*mm,*m2,*r,*mend; + mm=state->randmem; r=state->randrsl; + a = state->randa; b = (state->randb + (++state->randc)) & 0xffffffff; + for (m = mm, mend = m2 = m+(RK_ISAAC_STATE_LEN/2); m>6 , a, b, mm, m, m2, r, x); + rngstep( a<<2 , a, b, mm, m, m2, r, x); + rngstep( a>>16, a, b, mm, m, m2, r, x); + } + for (m2 = mm; m2>6 , a, b, mm, m, m2, r, x); + rngstep( a<<2 , a, b, mm, m, m2, r, x); + rngstep( a>>16, a, b, mm, m, m2, r, x); + } + state->randb = b; state->randa = a; + state->randcnt=RK_ISAAC_STATE_LEN-1; + } + return state->randrsl[state->randcnt] & 0xFFFFFFFF; +} + +#define mix(a,b,c,d,e,f,g,h) \ +{ \ + a^=b<<11; d+=a; b+=c; \ + b^=(c & 0xFFFFFFFF)>>2; e+=b; c+=d; \ + c^=d<<8; f+=c; d+=e; \ + d^=(e & 0xFFFFFFFF)>>16; g+=d; e+=f; \ + e^=f<<10; h+=e; f+=g; \ + f^=(g & 0xFFFFFFFF)>>4; a+=f; g+=h; \ + g^=h<<8; b+=g; h+=a; \ + h^=(a & 0xFFFFFFFF)>>9; c+=h; a+=b; \ +} + +/* if (flag==1), then use the contents of randrsl[] to initialize mm[]. */ +void isaac_init(rk_isaac_state *rk_isaac_state) +{ + int i; + unsigned long a,b,c,d,e,f,g,h; + unsigned long *m,*r; + rk_isaac_state->randa = rk_isaac_state->randb = rk_isaac_state->randc = 0; + m=rk_isaac_state->randmem; + r=rk_isaac_state->randrsl; + a=b=c=d=e=f=g=h=0x9e3779b9; /* the golden ratio */ + + for (i=0; i<4; ++i) /* scramble it */ + mix(a,b,c,d,e,f,g,h); + + /* initialize using the contents of r[] as the seed */ + for (i=0; irandcnt=0; + rk_isaac_state->has_gauss=0; +} + +long rk_isaac_long(rk_isaac_state *state) +{ + return rk_isaac_ulong(state) >> 1; +} + +unsigned long rk_isaac_ulong(rk_isaac_state *state) +{ +#if ULONG_MAX <= 0xffffffffUL + return rk_isaac_random(state); +#else + /* Assumes 64 bits */ + return (rk_isaac_random(state) << 32) | (rk_isaac_random(state)); +#endif +} + +unsigned long rk_isaac_interval(unsigned long max, rk_isaac_state *state) +{ + unsigned long mask = max, value; + + if (max == 0) return 0; + + /* Smallest bit mask >= max */ + mask |= mask >> 1; + mask |= mask >> 2; + mask |= mask >> 4; + mask |= mask >> 8; + mask |= mask >> 16; +#if ULONG_MAX > 0xffffffffUL + mask |= mask >> 32; +#endif + + /* Search a random value in [0..mask] <= max */ + while ((value = (rk_isaac_ulong(state) & mask)) > max); + + return value; +} + +double rk_isaac_double(rk_isaac_state *state) +{ + /* shifts : 67108864 = 0x4000000, 9007199254740992 = 0x20000000000000 */ + long a = rk_isaac_random(state) >> 5, b = rk_isaac_random(state) >> 6; + return (a * 67108864.0 + b) / 9007199254740992.0; +} + +void rk_isaac_copy(rk_isaac_state *copy, rk_isaac_state *orig) +{ + memcpy(copy, orig, sizeof(rk_isaac_state)); +} + +double rk_isaac_gauss(rk_isaac_state *state) +{ + if (state->has_gauss) + { + state->has_gauss = 0; + return state->gauss; + } + else + { + double f, x1, x2, r2; + do + { + x1 = 2.0*rk_isaac_double(state) - 1.0; + x2 = 2.0*rk_isaac_double(state) - 1.0; + r2 = x1*x1 + x2*x2; + } + while (r2 >= 1.0 || r2 == 0.0); + + f = sqrt(-2.0*log(r2)/r2); /* Box-Muller transform */ + state->has_gauss = 1; + state->gauss = f*x1; /* Keep for next call */ + return f*x2; + } +} + +void rk_isaac_fill(void *buffer, size_t size, rk_isaac_state *state) +{ + unsigned long r; + unsigned char *buf = buffer; + rk_isaac_state tempstate; + + if (size > 0 && state == NULL) + { + rk_isaac_randomseed(&tempstate); + state = &tempstate; + } + + for (; size >= 4; size -= 4) + { + r = rk_isaac_random(state); + *(buf++) = r & 0xFF; + *(buf++) = (r >> 8) & 0xFF; + *(buf++) = (r >> 16) & 0xFF; + *(buf++) = (r >> 24) & 0xFF; + } + + if (!size) return; + + r = rk_isaac_random(state); + + for (; size; r >>= 8, size --) + *(buf++) = (unsigned char)(r & 0xFF); +} + +void rk_seed_isaac(rk_isaac_state *i_state, rk_state *state) +{ + rk_isaac_fill(state->key, sizeof(state->key), i_state); + state->pos = RK_STATE_LEN; + state->has_gauss = 0; +} diff --git a/tardis/montecarlo/src/randomkit/rk_isaac.h b/tardis/montecarlo/src/randomkit/rk_isaac.h new file mode 100644 index 00000000000..c8dfdcb451e --- /dev/null +++ b/tardis/montecarlo/src/randomkit/rk_isaac.h @@ -0,0 +1,142 @@ +/* Random kit 1.6 */ + +/* + * Copyright (c) 2004-2006, Jean-Sebastien Roy (js@jeannot.org) + * + * ISAAC RNG by By Bob Jenkins. Based on Bob Jenkins public domain code. + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. + * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY + * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, + * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE + * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +/* @(#) $Jeannot: rk_isaac.h,v 1.2 2006/02/19 14:40:26 js Exp $ */ + +/* + * Typical use: + * + * { + * rk_isaac_state state; + * unsigned long seed = 1, random_value; + * + * rk_isaac_seed(seed, &state); // Initialize the RNG + * ... + * random_value = rk_isaac_random(&state); // a random value in [0..RK_MAX] + * } + * + * Instead of rk_isaac_seed, you can use rk_isaac_randomseed which will get a + * random seed from /dev/random (or the clock, if /dev/random is unavailable): + * + * { + * rk_isaac_state state; + * unsigned long random_value; + * + * rk_isaac_randomseed(&state); // Initialize the RNG with a random seed + * ... + * random_value = rk_isaac_random(&state); // a random value in [0..RK_MAX] + * } + */ + +#ifndef _RK_ISAAC_ +#define _RK_ISAAC_ + +#include "rk_mt.h" + +#define RK_ISAAC_STATE_POW 8 +#define RK_ISAAC_STATE_LEN (1< +#include +#include +#include +#include +#include +#include +#include + +#ifdef _WIN32 +/* Windows */ +#include +#ifndef RK_NO_WINCRYPT +/* Windows crypto */ +#ifndef _WIN32_WINNT +#define _WIN32_WINNT 0x0400 +#endif +#include +#include +#endif +#else +/* Unix */ +#include +#include +#endif + +#include "randomkit.h" + +#ifndef RK_DEV_URANDOM +#define RK_DEV_URANDOM "/dev/urandom" +#endif + +#ifndef RK_DEV_RANDOM +#define RK_DEV_RANDOM "/dev/random" +#endif + +char *rk_strerror[RK_ERR_MAX] = +{ + "no error", + "random device unvavailable" +}; + +/* static functions */ +static unsigned long rk_hash(unsigned long key); + +void rk_knuth_fill(unsigned long seed, unsigned long *key, unsigned long len) +{ + int pos; + seed &= 0xffffffffUL; + + /* Knuth's PRNG as used in the Mersenne Twister reference implementation */ + for (pos=0; pos> 30)) + pos + 1) & 0xffffffffUL; + } +} + +void rk_seed(unsigned long seed, rk_state *state) +{ + rk_knuth_fill(seed, state->key, RK_STATE_LEN); + state->pos = RK_STATE_LEN; + state->has_gauss = 0; +} + +/* Thomas Wang 32 bits integer hash function */ +unsigned long rk_hash(unsigned long key) +{ + key += ~(key << 15); + key ^= (key >> 10); + key += (key << 3); + key ^= (key >> 6); + key += ~(key << 11); + key ^= (key >> 16); + return key; +} + +unsigned long rk_seedfromsystem() +{ +#ifndef _WIN32 + struct timeval tv; +#else + struct _timeb tv; +#endif + +#ifndef _WIN32 + gettimeofday(&tv, NULL); + return rk_hash(getpid()) ^ rk_hash(tv.tv_sec) ^ rk_hash(tv.tv_usec) + ^ rk_hash(clock()); +#else + _ftime(&tv); + return rk_hash(tv.time) ^ rk_hash(tv.millitm) ^ rk_hash(clock()); +#endif +} + +rk_error rk_randomseed(rk_state *state) +{ + if(rk_devfill(state->key, sizeof(state->key), 0) == RK_NOERR) + { + state->key[0] |= 0x80000000UL; /* ensures non-zero key */ + state->pos = RK_STATE_LEN; + state->has_gauss = 0; + return RK_NOERR; + } + + rk_seed(rk_seedfromsystem(), state); + + return RK_ENODEV; +} + +/* Magic Mersenne Twister constants */ +#define N 624 +#define M 397 +#define MATRIX_A 0x9908b0dfUL +#define UPPER_MASK 0x80000000UL +#define LOWER_MASK 0x7fffffffUL + +/* Slightly optimised reference implementation of the Mersenne Twister */ +unsigned long rk_random(rk_state *state) +{ + unsigned long y; + + if (state->pos == RK_STATE_LEN) + { + int i; + + for (i=0;ikey[i] & UPPER_MASK) | (state->key[i+1] & LOWER_MASK); + state->key[i] = state->key[i+M] ^ (y>>1) ^ (-(y & 1) & MATRIX_A); + } + for (;ikey[i] & UPPER_MASK) | (state->key[i+1] & LOWER_MASK); + state->key[i] = state->key[i+(M-N)] ^ (y>>1) ^ (-(y & 1) & MATRIX_A); + } + y = (state->key[N-1] & UPPER_MASK) | (state->key[0] & LOWER_MASK); + state->key[N-1] = state->key[M-1] ^ (y>>1) ^ (-(y & 1) & MATRIX_A); + + state->pos = 0; + } + + y = state->key[state->pos++]; + + /* Tempering */ + y ^= (y >> 11); + y ^= (y << 7) & 0x9d2c5680UL; + y ^= (y << 15) & 0xefc60000UL; + y ^= (y >> 18); + + return y; +} + +long rk_long(rk_state *state) +{ + return rk_ulong(state) >> 1; +} + +unsigned long rk_ulong(rk_state *state) +{ +#if ULONG_MAX <= 0xffffffffUL + return rk_random(state); +#else + /* Assumes 64 bits */ + return (rk_random(state) << 32) | (rk_random(state)); +#endif +} + +unsigned long rk_interval(unsigned long max, rk_state *state) +{ + unsigned long mask = max, value; + + if (max == 0) return 0; + + /* Smallest bit mask >= max */ + mask |= mask >> 1; + mask |= mask >> 2; + mask |= mask >> 4; + mask |= mask >> 8; + mask |= mask >> 16; +#if ULONG_MAX > 0xffffffffUL + mask |= mask >> 32; +#endif + + /* Search a random value in [0..mask] <= max */ + while ((value = (rk_ulong(state) & mask)) > max); + + return value; +} + +double rk_double(rk_state *state) +{ + /* shifts : 67108864 = 0x4000000, 9007199254740992 = 0x20000000000000 */ + long a = rk_random(state) >> 5, b = rk_random(state) >> 6; + return (a * 67108864.0 + b) / 9007199254740992.0; +} + +void rk_copy(rk_state *copy, rk_state *orig) +{ + memcpy(copy, orig, sizeof(rk_state)); +} + +void rk_fill(void *buffer, size_t size, rk_state *state) +{ + unsigned long r; + unsigned char *buf = buffer; + rk_state tempstate; + + if (size > 0 && state == NULL) + { + rk_randomseed(&tempstate); + state = &tempstate; + } + + for (; size >= 4; size -= 4) + { + r = rk_random(state); + *(buf++) = r & 0xFF; + *(buf++) = (r >> 8) & 0xFF; + *(buf++) = (r >> 16) & 0xFF; + *(buf++) = (r >> 24) & 0xFF; + } + + if (!size) return; + + r = rk_random(state); + + for (; size; r >>= 8, size --) + *(buf++) = (unsigned char)(r & 0xFF); +} + +rk_error rk_devfill(void *buffer, size_t size, int strong) +{ +#ifndef _WIN32 + FILE *rfile; + int done; + + if (strong) + rfile = fopen(RK_DEV_RANDOM, "rb"); + else + rfile = fopen(RK_DEV_URANDOM, "rb"); + if (rfile == NULL) + return RK_ENODEV; + done = fread(buffer, size, 1, rfile); + fclose(rfile); + if (done) + return RK_NOERR; +#else + +#ifndef RK_NO_WINCRYPT + HCRYPTPROV hCryptProv; + BOOL done; + + if (!CryptAcquireContext(&hCryptProv, NULL, NULL, PROV_RSA_FULL, + CRYPT_VERIFYCONTEXT) || !hCryptProv) + return RK_ENODEV; + done = CryptGenRandom(hCryptProv, size, (unsigned char *)buffer); + CryptReleaseContext(hCryptProv, 0); + if (done) + return RK_NOERR; +#endif + +#endif + + return RK_ENODEV; +} + +rk_error rk_altfill(void *buffer, size_t size, int strong, rk_state *state) +{ + rk_error err; + + err = rk_devfill(buffer, size, strong); + if (err) + rk_fill(buffer, size, state); + + return err; +} + +double rk_gauss(rk_state *state) +{ + if (state->has_gauss) + { + state->has_gauss = 0; + return state->gauss; + } + else + { + double f, x1, x2, r2; + do + { + x1 = 2.0*rk_double(state) - 1.0; + x2 = 2.0*rk_double(state) - 1.0; + r2 = x1*x1 + x2*x2; + } + while (r2 >= 1.0 || r2 == 0.0); + + f = sqrt(-2.0*log(r2)/r2); /* Box-Muller transform */ + state->has_gauss = 1; + state->gauss = f*x1; /* Keep for next call */ + return f*x2; + } +} diff --git a/tardis/montecarlo/src/randomkit/rk_mt.h b/tardis/montecarlo/src/randomkit/rk_mt.h new file mode 100644 index 00000000000..ff0e22a0879 --- /dev/null +++ b/tardis/montecarlo/src/randomkit/rk_mt.h @@ -0,0 +1,194 @@ +/* Random kit 1.6 */ + +/* + * Anyone who attempts to generate random numbers by deterministic means is, + * of course, living in a state of sin. + * -- John von Neumann + */ + +/* + * Copyright (c) 2003-2006, Jean-Sebastien Roy (js@jeannot.org) + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. + * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY + * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, + * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE + * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +/* @(#) $Jeannot: rk_mt.h,v 1.6 2006/02/20 19:13:37 js Exp $ */ + +/* + * Typical use: + * + * { + * rk_state state; + * unsigned long seed = 1, random_value; + * + * rk_seed(seed, &state); // Initialize the RNG + * ... + * random_value = rk_random(&state); // a random value in [0..RK_MAX] + * } + * + * Instead of rk_seed, you can use rk_randomseed which will get a random seed + * from /dev/urandom (or the clock, if /dev/urandom is unavailable): + * + * { + * rk_state state; + * unsigned long random_value; + * + * rk_randomseed(&state); // Initialize the RNG with a random seed + * ... + * random_value = rk_random(&state); // a random value in [0..RK_MAX] + * } + */ + +/* + * Useful macro: + * RK_DEV_RANDOM: the device used for random seeding. + * defaults to "/dev/urandom" + */ + +#include + +#ifndef _RK_MT_ +#define _RK_MT_ + +#define RK_STATE_LEN 624 + +typedef struct rk_state_ +{ + unsigned long key[RK_STATE_LEN]; + int pos; + int has_gauss; /* !=0: gauss contains a gaussian deviate */ + double gauss; +} +rk_state; + +typedef enum { + RK_NOERR = 0, /* no error */ + RK_ENODEV = 1, /* no RK_DEV_RANDOM device */ + RK_ERR_MAX = 2 +} rk_error; + +/* error strings */ +extern char *rk_strerror[RK_ERR_MAX]; + +/* Maximum generated random value */ +#define RK_MAX 0xFFFFFFFFUL + +#ifdef __cplusplus +extern "C" { +#endif + +/* + * Initialize the RNG state using the given seed. + */ +extern void rk_seed(unsigned long seed, rk_state *state); + +/* + * Initialize the RNG state using a random seed. + * Uses /dev/urandom or, when unavailable, the clock (see rk_mt.c). + * Returns RK_NOERR when no errors occurs. + * Returns RK_ENODEV when the use of RK_DEV_RANDOM failed (for example because + * there is no such device). In this case, the RNG was initialized using the + * clock. + */ +extern rk_error rk_randomseed(rk_state *state); + +/* + * Returns a random unsigned long between 0 and RK_MAX inclusive + */ +extern unsigned long rk_random(rk_state *state); + +/* + * Returns a random long between 0 and LONG_MAX inclusive + */ +extern long rk_long(rk_state *state); + +/* + * Returns a random unsigned long between 0 and ULONG_MAX inclusive + */ +extern unsigned long rk_ulong(rk_state *state); + +/* + * Returns a random unsigned long between 0 and max inclusive. + */ +extern unsigned long rk_interval(unsigned long max, rk_state *state); + +/* + * Returns a random double between 0.0 and 1.0, 1.0 excluded. + */ +extern double rk_double(rk_state *state); + +/* + * Copy a random generator state. + */ +extern void rk_copy(rk_state *copy, rk_state *orig); + +/* + * fill the buffer with size random bytes. + * If state == NULL, the random generator is inilialized using rk_randomseed. + * Calling multiple times rk_randomseed should be avoided therefore calling + * multiple times rk_fill with state == NULL should be avoided. + */ +extern void rk_fill(void *buffer, size_t size, rk_state *state); + +/* + * fill the buffer with randombytes from the random device + * Returns RK_ENODEV if the device is unavailable, or RK_NOERR if it is + * On Unix, if strong is defined, RK_DEV_RANDOM is used. If not, RK_DEV_URANDOM + * is used instead. This parameter has no effect on Windows. + * Warning: on most unixes RK_DEV_RANDOM will wait for enough entropy to answer + * which can take a very long time on quiet systems. + */ +extern rk_error rk_devfill(void *buffer, size_t size, int strong); + +/* + * fill the buffer using rk_devfill if the random device is available and using + * rk_fill if is is not + * parameters have the same meaning as rk_fill and rk_devfill + * Returns RK_ENODEV if the device is unavailable, or RK_NOERR if it is + */ +extern rk_error rk_altfill(void *buffer, size_t size, int strong, + rk_state *state); + +/* + * return a random gaussian deviate with variance unity and zero mean. + */ +extern double rk_gauss(rk_state *state); + +/* Utility functions */ + +/* + * fill the key vector using Knuth RNG as used in MT reference implementation + * using the provided seed. The key vector length is len. + */ +extern void rk_knuth_fill(unsigned long seed, unsigned long *key, + unsigned long len); + +/* + * return a random unsigned long based upon the system state (clock, pid) + * used as seed when there is no random device. + */ +extern unsigned long rk_seedfromsystem(void); + + +#ifdef __cplusplus +} +#endif + +#endif /* _RK_MT_ */ diff --git a/tardis/montecarlo/src/randomkit/rk_primitive.c b/tardis/montecarlo/src/randomkit/rk_primitive.c new file mode 100644 index 00000000000..25e70c6c7f5 --- /dev/null +++ b/tardis/montecarlo/src/randomkit/rk_primitive.c @@ -0,0 +1,520 @@ +/* Random kit 1.6 */ +/* Primitivity test for binary polynomials of low degree */ + +/* + * Copyright (c) 2005-2006, Jean-Sebastien Roy (js@jeannot.org) + * + * Methodology inspired by scott duplichan's ppsearch code. + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. + * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY + * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, + * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE + * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +static char const rcsid[] = + "@(#) $Jeannot: rk_primitive.c,v 1.7 2006/02/19 13:48:34 js Exp $"; + +#include +#include +#include "rk_primitive.h" + +#ifndef LONG_BIT +#if ULONG_MAX <= 0xffffffffUL +#define LONG_BIT 32 +#else +#define LONG_BIT 64 +#endif +#endif + +static unsigned long modmul(unsigned long poly1, unsigned long poly2, + unsigned long modulo, unsigned long mask); +static unsigned long modpow(unsigned long polynomial, unsigned long power, + unsigned long modulo, int degree); + +/* + * For all powers i of two up to 64, list all the number of the form + * (2^i-1)/p for p a prime factor of 2^i-1. + */ +static const unsigned long divisors[][12]={ +/* 2^0-1 */ + {1UL, + 0UL}, +/* 2^1-1 */ + {1UL, + 0UL}, +/* 2^2-1 */ + {1UL, + 0UL}, +/* 2^3-1 */ + {1UL, + 0UL}, +/* 2^4-1 */ + {5UL, + 3UL, + 0UL}, +/* 2^5-1 */ + {1UL, + 0UL}, +/* 2^6-1 */ + {21UL, + 9UL, + 0UL}, +/* 2^7-1 */ + {1UL, + 0UL}, +/* 2^8-1 */ + {85UL, + 51UL, + 15UL, + 0UL}, +/* 2^9-1 */ + {73UL, + 7UL, + 0UL}, +/* 2^10-1 */ + {341UL, + 93UL, + 33UL, + 0UL}, +/* 2^11-1 */ + {89UL, + 23UL, + 0UL}, +/* 2^12-1 */ + {1365UL, + 819UL, + 585UL, + 315UL, + 0UL}, +/* 2^13-1 */ + {1UL, + 0UL}, +/* 2^14-1 */ + {5461UL, + 381UL, + 129UL, + 0UL}, +/* 2^15-1 */ + {4681UL, + 1057UL, + 217UL, + 0UL}, +/* 2^16-1 */ + {21845UL, + 13107UL, + 3855UL, + 255UL, + 0UL}, +/* 2^17-1 */ + {1UL, + 0UL}, +/* 2^18-1 */ + {87381UL, + 37449UL, + 13797UL, + 3591UL, + 0UL}, +/* 2^19-1 */ + {1UL, + 0UL}, +/* 2^20-1 */ + {349525UL, + 209715UL, + 95325UL, + 33825UL, + 25575UL, + 0UL}, +/* 2^21-1 */ + {299593UL, + 16513UL, + 6223UL, + 0UL}, +/* 2^22-1 */ + {1398101UL, + 182361UL, + 47127UL, + 6141UL, + 0UL}, +/* 2^23-1 */ + {178481UL, + 47UL, + 0UL}, +/* 2^24-1 */ + {5592405UL, + 3355443UL, + 2396745UL, + 1290555UL, + 986895UL, + 69615UL, + 0UL}, +/* 2^25-1 */ + {1082401UL, + 55831UL, + 18631UL, + 0UL}, +/* 2^26-1 */ + {22369621UL, + 24573UL, + 8193UL, + 0UL}, +/* 2^27-1 */ + {19173961UL, + 1838599UL, + 511UL, + 0UL}, +/* 2^28-1 */ + {89478485UL, + 53687091UL, + 9256395UL, + 6242685UL, + 2375535UL, + 2113665UL, + 0UL}, +/* 2^29-1 */ + {2304167UL, + 486737UL, + 256999UL, + 0UL}, +/* 2^30-1 */ + {357913941UL, + 153391689UL, + 97612893UL, + 34636833UL, + 7110873UL, + 3243933UL, + 0UL}, +/* 2^31-1 */ + {1UL, + 0UL}, +/* 2^32-1 */ + {1431655765UL, + 858993459UL, + 252645135UL, + 16711935UL, + 65535UL, + 0UL}, +#if LONG_BIT > 32 +/* 2^33-1 */ + {1227133513UL, + 373475417UL, + 96516119UL, + 14329UL, + 0UL}, +/* 2^34-1 */ + {5726623061UL, + 393213UL, + 131073UL, + 0UL}, +/* 2^35-1 */ + {1108378657UL, + 483939977UL, + 270549121UL, + 279527UL, + 0UL}, +/* 2^36-1 */ + {22906492245UL, + 13743895347UL, + 9817068105UL, + 5286113595UL, + 3616814565UL, + 1857283155UL, + 941362695UL, + 630453915UL, + 0UL}, +/* 2^37-1 */ + {616318177UL, + 223UL, + 0UL}, +/* 2^38-1 */ + {91625968981UL, + 1572861UL, + 524289UL, + 0UL}, +/* 2^39-1 */ + {78536544841UL, + 6958934353UL, + 67117057UL, + 4529623UL, + 0UL}, +/* 2^40-1 */ + {366503875925UL, + 219902325555UL, + 99955602525UL, + 64677154575UL, + 35468117025UL, + 26817356775UL, + 17825775UL, + 0UL}, +/* 2^41-1 */ + {164511353UL, + 13367UL, + 0UL}, +/* 2^42-1 */ + {1466015503701UL, + 628292358729UL, + 102280151421UL, + 34630287489UL, + 13050583119UL, + 811597437UL, + 0UL}, +/* 2^43-1 */ + {20408568497UL, + 905040953UL, + 4188889UL, + 0UL}, +/* 2^44-1 */ + {5864062014805UL, + 3518437208883UL, + 764877654105UL, + 197665011735UL, + 44312811195UL, + 25757227005UL, + 8325691455UL, + 0UL}, +/* 2^45-1 */ + {5026338869833UL, + 1134979744801UL, + 481977699847UL, + 233009086681UL, + 55759702201UL, + 1509346321UL, + 0UL}, +/* 2^46-1 */ + {23456248059221UL, + 1497207322929UL, + 394264623UL, + 25165821UL, + 0UL}, +/* 2^47-1 */ + {59862819377UL, + 31184907679UL, + 10610063UL, + 0UL}, +/* 2^48-1 */ + {93824992236885UL, + 56294995342131UL, + 40210710958665UL, + 21651921285435UL, + 16557351571215UL, + 2901803883615UL, + 1167945961455UL, + 1095233372415UL, + 418239192735UL, + 0UL}, +/* 2^49-1 */ + {4432676798593UL, + 127UL, + 0UL}, +/* 2^50-1 */ + {375299968947541UL, + 102354536985693UL, + 36319351833633UL, + 4485656999373UL, + 1873377548823UL, + 625152641223UL, + 277931351973UL, + 0UL}, +/* 2^51-1 */ + {321685687669321UL, + 21862134113449UL, + 1050769861729UL, + 202518195313UL, + 17180000257UL, + 0UL}, +/* 2^52-1 */ + {1501199875790165UL, + 900719925474099UL, + 84973577874915UL, + 28685347945035UL, + 2792064245115UL, + 1649066139645UL, + 549822930945UL, + 0UL}, +/* 2^53-1 */ + {1416003655831UL, + 129728784761UL, + 441650591UL, + 0UL}, +/* 2^54-1 */ + {6004799503160661UL, + 2573485501354569UL, + 948126237341157UL, + 246772582321671UL, + 206561081853UL, + 68585259519UL, + 0UL}, +/* 2^55-1 */ + {1566469435607129UL, + 1162219258676257UL, + 404817944033303UL, + 40895342813807UL, + 11290754314937UL, + 178394823847UL, + 0UL}, +/* 2^56-1 */ + {24019198012642645UL, + 14411518807585587UL, + 4238682002231055UL, + 2484744621997515UL, + 1675758000882045UL, + 637677823344495UL, + 567382630219905UL, + 4563402735UL, + 0UL}, +/* 2^57-1 */ + {20587884010836553UL, + 4451159405623UL, + 274878431233UL, + 118823881393UL, + 0UL}, +/* 2^58-1 */ + {96076792050570581UL, + 4885260612740877UL, + 1237040240994471UL, + 261314937580881UL, + 137975287770087UL, + 95026151247UL, + 0UL}, +/* 2^59-1 */ + {3203431780337UL, + 179951UL, + 0UL}, +/* 2^60-1 */ + {384307168202282325UL, + 230584300921369395UL, + 164703072086692425UL, + 104811045873349725UL, + 88686269585142075UL, + 37191016277640225UL, + 28120036697727975UL, + 18900352534538475UL, + 7635241752363225UL, + 3483146539597725UL, + 872764197279975UL, + 0UL}, +/* 2^61-1 */ + {1UL, + 0UL}, +/* 2^62-1 */ + {1537228672809129301UL, + 6442450941UL, + 2147483649UL, + 0UL}, +/* 2^63-1 */ + {1317624576693539401UL, + 126347562148695559UL, + 72624976668147841UL, + 27369056489183311UL, + 99457304386111UL, + 14197294936951UL, + 0UL}, +/* 2^64-1 */ + {6148914691236517205UL, + 3689348814741910323UL, + 1085102592571150095UL, + 71777214294589695UL, + 28778071877862015UL, + 281470681808895UL, + 2753074036095UL, + 0UL}, +#if LONG_BIT > 64 +#error Factorization of numbers up to 2^LONG_BIT required +#endif +#endif +}; + +/* + * Modular multiply for two binary polynomial + * mask is 1UL << the degree of the modulus. + */ +unsigned long modmul(unsigned long poly1, unsigned long poly2, + unsigned long modulo, unsigned long mask) +{ + unsigned long result = 0; + + for (; poly1; poly1 >>= 1) + { + if (poly1 & 1) + result ^= poly2; + + poly2 <<= 1; + if (poly2 & mask) + poly2 ^= modulo; + } + return result; +} + +/* + * Modular exponentiation for a binary polynomial + * degree is the degree of the modulus. + */ +unsigned long modpow(unsigned long polynomial, unsigned long power, + unsigned long modulo, int degree) +{ + unsigned long result = 1, mask = 1UL << degree; + + for (; power; power >>= 1) + { + if (power & 1) + result = modmul(result, polynomial, modulo, mask); + polynomial = modmul(polynomial, polynomial, modulo, mask); + } + return result; +} + +/* + * Test the primitivity of a polynomial + */ +int rk_isprimitive(unsigned long polynomial) +{ + unsigned long pelement = 2, temp = polynomial >> 1; + int k, degree = 0, weight = 1; + + /* Special case for polynomials of degree < 2 */ + if (polynomial < 4) + return (polynomial == 3) || (polynomial == 1); + + /* A binary primitive polynomial has a constant term */ + if (!(polynomial & 1)) + return 0; + + /* + * A binary primitive polynomial of degree > 1 has an odd number of terms. + * temp ^= temp >> 16; temp ^= temp >> 8; ... would be sligthly faster. + * Compute the degree at the same time. + */ + for (; temp; degree++, temp >>= 1) + weight += temp & 1; + if (!(weight & 1)) + return 0; + + /* + * Check if the period is 2^degree-1. + * Sufficient if 2^degree-1 is prime. + */ + if (modpow(pelement, 1UL << degree, polynomial, degree) != pelement) + return 0; + + if (divisors[degree][0] != 1) + /* Primitivity test */ + for (k = 0; divisors[degree][k]; k++) + if (modpow(pelement, divisors[degree][k], polynomial, degree) == 1) + return 0; + + return 1; +} diff --git a/tardis/montecarlo/src/randomkit/rk_primitive.h b/tardis/montecarlo/src/randomkit/rk_primitive.h new file mode 100644 index 00000000000..f7ab374cef5 --- /dev/null +++ b/tardis/montecarlo/src/randomkit/rk_primitive.h @@ -0,0 +1,54 @@ +/* Random kit 1.6 */ +/* Primitivity test for binary polynomials of low degree */ + +/* + * Copyright (c) 2005-2006, Jean-Sebastien Roy (js@jeannot.org) + * + * Methodology inspired by scott duplichan's ppsearch code. + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. + * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY + * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, + * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE + * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +/* @(#) $Jeannot: rk_primitive.h,v 1.6 2006/02/19 13:48:34 js Exp $ */ + +#ifndef _RK_PRIMITIVE_ +#define _RK_PRIMITIVE_ + +#ifdef __cplusplus +extern "C" { +#endif + +/* + * Return 1 if the binary polynomial is primitive. + * + * Note that if p is primitive, the the polynomial obtained by reversing the + * bits of p is also primitive. (see list_primitive.c for an example) + * + * Typical use: + * int test; + * test = rk_isprimitive(3, &divisors); + */ +extern int rk_isprimitive(unsigned long polynomial); + +#ifdef __cplusplus +} +#endif + +#endif /* _RK_PRIMITIVE_ */ diff --git a/tardis/montecarlo/src/randomkit/rk_sobol.c b/tardis/montecarlo/src/randomkit/rk_sobol.c new file mode 100644 index 00000000000..c6a016bb033 --- /dev/null +++ b/tardis/montecarlo/src/randomkit/rk_sobol.c @@ -0,0 +1,991 @@ +/* Random kit 1.6 */ + +/* + * Copyright (c) 2004-2006, Jean-Sebastien Roy (js@jeannot.org) + * + * Original algorithm from Numerical Recipes, 2nd edition, by Press et al. + * The inverse normal cdf formulas are from Peter J. Acklam. + * The initialization directions were found in Ferdinando Ametrano's + * implementation in QuantLib. + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. + * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY + * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, + * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE + * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +static char const rcsid[] = + "@(#) $Jeannot: rk_sobol.c,v 1.9 2006/02/19 13:48:34 js Exp $"; + +#include +#include +#include +#include "rk_sobol.h" +#include "rk_mt.h" +#include "rk_primitive.h" + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif +#ifndef M_SQRT1_2 +#define M_SQRT1_2 0.70710678118654752440 /* 1/sqrt(2) */ +#endif +#define RK_SOBOL_M_SQRT2PI 2.506628274631000502415 /* sqrt(2*pi) */ + +#ifndef LONG_BIT +#if ULONG_MAX <= 0xffffffffUL +#define LONG_BIT 32 +#else +#define LONG_BIT 64 +#endif +#endif + +char *rk_sobol_strerror[] = +{ + "no error", + "invalid dimension", + "too many numbers generated", + "not enough memory" +}; + +static double inverse_normal(double p); + +/* + * Sobol/Levitan coefficients of the free direction integers as given + * by Bratley, P., Fox, B.L. (1988) + */ + +const unsigned long rk_sobol_SLdirections[] = { + 1, + 1, 1, + 1, 3, 7, + 1, 1, 5, + 1, 3, 1, 1, + 1, 1, 3, 7, + 1, 3, 3, 9, 9, + 1, 3, 7, 13, 3, + 1, 1, 5, 11, 27, + 1, 3, 5, 1, 15, + 1, 1, 7, 3, 29, + 1, 3, 7, 7, 21, + 1, 1, 1, 9, 23, 37, + 1, 3, 3, 5, 19, 33, + 1, 1, 3, 13, 11, 7, + 1, 1, 7, 13, 25, 5, + 1, 3, 5, 11, 7, 11, + 1, 1, 1, 3, 13, 39, + 1, 3, 1, 15, 17, 63, 13, + 1, 1, 5, 5, 1, 27, 33, + 1, 3, 3, 3, 25, 17, 115, + 1, 1, 3, 15, 29, 15, 41, + 1, 3, 1, 7, 3, 23, 79, + 1, 3, 7, 9, 31, 29, 17, + 1, 1, 5, 13, 11, 3, 29, + 1, 3, 1, 9, 5, 21, 119, + 1, 1, 3, 1, 23, 13, 75, + 1, 3, 3, 11, 27, 31, 73, + 1, 1, 7, 7, 19, 25, 105, + 1, 3, 5, 5, 21, 9, 7, + 1, 1, 1, 15, 5, 49, 59, + 1, 1, 1, 1, 1, 33, 65, + 1, 3, 5, 15, 17, 19, 21, + 1, 1, 7, 11, 13, 29, 3, + 1, 3, 7, 5, 7, 11, 113, + 1, 1, 5, 3, 15, 19, 61, + 1, 3, 1, 1, 9, 27, 89, 7, + 1, 1, 3, 7, 31, 15, 45, 23, + 1, 3, 3, 9, 9, 25, 107, 39, + 0 +}; + +/* + * Lemieux coefficients of the free direction integers as given + * in QuantLib by Christiane Lemieux, private communication, September 2004 + */ + +const unsigned long rk_sobol_Ldirections[] = { + 1, + 1, 1, + 1, 3, 7, + 1, 1, 5, + 1, 3, 1, 1, + 1, 1, 3, 7, + 1, 3, 3, 9, 9, + 1, 3, 7, 13, 3, + 1, 1, 5, 11, 27, + 1, 3, 5, 1, 15, + 1, 1, 7, 3, 29, + 1, 3, 7, 7, 21, + 1, 1, 1, 9, 23, 37, + 1, 3, 3, 5, 19, 33, + 1, 1, 3, 13, 11, 7, + 1, 1, 7, 13, 25, 5, + 1, 3, 5, 11, 7, 11, + 1, 1, 1, 3, 13, 39, + 1, 3, 1, 15, 17, 63, 13, + 1, 1, 5, 5, 1, 27, 33, + 1, 3, 3, 3, 25, 17, 115, + 1, 1, 3, 15, 29, 15, 41, + 1, 3, 1, 7, 3, 23, 79, + 1, 3, 7, 9, 31, 29, 17, + 1, 1, 5, 13, 11, 3, 29, + 1, 3, 1, 9, 5, 21, 119, + 1, 1, 3, 1, 23, 13, 75, + 1, 3, 3, 11, 27, 31, 73, + 1, 1, 7, 7, 19, 25, 105, + 1, 3, 5, 5, 21, 9, 7, + 1, 1, 1, 15, 5, 49, 59, + 1, 1, 1, 1, 1, 33, 65, + 1, 3, 5, 15, 17, 19, 21, + 1, 1, 7, 11, 13, 29, 3, + 1, 3, 7, 5, 7, 11, 113, + 1, 1, 5, 3, 15, 19, 61, + 1, 3, 1, 1, 9, 27, 89, 7, + 1, 1, 3, 7, 31, 15, 45, 23, + 1, 3, 3, 9, 9, 25, 107, 39, + 1, 1, 3, 13, 7, 35, 61, 91, + 1, 1, 7, 11, 5, 35, 55, 75, + 1, 3, 5, 5, 11, 23, 29, 139, + 1, 1, 1, 7, 11, 15, 17, 81, + 1, 1, 7, 9, 5, 57, 79, 103, + 1, 1, 7, 13, 19, 5, 5, 185, + 1, 3, 1, 3, 13, 57, 97, 131, + 1, 1, 5, 5, 21, 25, 125, 197, + 1, 3, 3, 9, 31, 11, 103, 201, + 1, 1, 5, 3, 7, 25, 51, 121, + 1, 3, 7, 15, 19, 53, 73, 189, + 1, 1, 1, 15, 19, 55, 27, 183, + 1, 1, 7, 13, 3, 29, 109, 69, + 1, 1, 5, 15, 15, 23, 15, 1, 57, + 1, 3, 1, 3, 23, 55, 43, 143, 397, + 1, 1, 3, 11, 29, 9, 35, 131, 411, + 1, 3, 1, 7, 27, 39, 103, 199, 277, + 1, 3, 7, 3, 19, 55, 127, 67, 449, + 1, 3, 7, 3, 5, 29, 45, 85, 3, + 1, 3, 5, 5, 13, 23, 75, 245, 453, + 1, 3, 1, 15, 21, 47, 3, 77, 165, + 1, 1, 7, 9, 15, 5, 117, 73, 473, + 1, 3, 1, 9, 1, 21, 13, 173, 313, + 1, 1, 7, 3, 11, 45, 63, 77, 49, + 1, 1, 1, 1, 1, 25, 123, 39, 259, + 1, 1, 1, 5, 23, 11, 59, 11, 203, + 1, 3, 3, 15, 21, 1, 73, 71, 421, + 1, 1, 5, 11, 15, 31, 115, 95, 217, + 1, 1, 3, 3, 7, 53, 37, 43, 439, + 1, 1, 1, 1, 27, 53, 69, 159, 321, + 1, 1, 5, 15, 29, 17, 19, 43, 449, + 1, 1, 3, 9, 1, 55, 121, 205, 255, + 1, 1, 3, 11, 9, 47, 107, 11, 417, + 1, 1, 1, 5, 17, 25, 21, 83, 95, + 1, 3, 5, 13, 31, 25, 61, 157, 407, + 1, 1, 7, 9, 25, 33, 41, 35, 17, + 1, 3, 7, 15, 13, 39, 61, 187, 461, + 1, 3, 7, 13, 5, 57, 23, 177, 435, + 1, 1, 3, 15, 11, 27, 115, 5, 337, + 1, 3, 7, 3, 15, 63, 61, 171, 339, + 1, 3, 3, 13, 15, 61, 59, 47, 1, + 1, 1, 5, 15, 13, 5, 39, 83, 329, + 1, 1, 5, 5, 5, 27, 25, 39, 301, + 1, 1, 5, 11, 31, 41, 35, 233, 27, + 1, 3, 5, 15, 7, 37, 119, 171, 419, + 1, 3, 5, 5, 3, 29, 21, 189, 417, + 1, 1, 1, 1, 21, 41, 117, 119, 351, + 1, 1, 3, 1, 7, 27, 87, 19, 213, + 1, 1, 1, 1, 17, 7, 97, 217, 477, + 1, 1, 7, 1, 29, 61, 103, 231, 269, + 1, 1, 7, 13, 9, 27, 107, 207, 311, + 1, 1, 7, 5, 25, 21, 107, 179, 423, + 1, 3, 5, 11, 7, 1, 17, 245, 281, + 1, 3, 5, 9, 1, 5, 53, 59, 125, + 1, 1, 7, 1, 31, 57, 71, 245, 125, + 1, 1, 7, 5, 5, 57, 53, 253, 441, + 1, 3, 1, 13, 19, 35, 119, 235, 381, + 1, 3, 1, 7, 19, 59, 115, 33, 361, + 1, 1, 3, 5, 13, 1, 49, 143, 501, + 1, 1, 3, 5, 1, 63, 101, 85, 189, + 1, 1, 5, 11, 27, 63, 13, 131, 5, + 1, 1, 5, 7, 15, 45, 75, 59, 455, 585, + 1, 3, 1, 3, 7, 7, 111, 23, 119, 959, + 1, 3, 3, 9, 11, 41, 109, 163, 161, 879, + 1, 3, 5, 1, 21, 41, 121, 183, 315, 219, + 1, 1, 3, 9, 15, 3, 9, 223, 441, 929, + 1, 1, 7, 9, 3, 5, 93, 57, 253, 457, + 1, 1, 7, 13, 15, 29, 83, 21, 35, 45, + 1, 1, 3, 7, 13, 61, 119, 219, 85, 505, + 1, 1, 3, 3, 17, 13, 35, 197, 291, 109, + 1, 1, 3, 3, 5, 1, 113, 103, 217, 253, + 1, 1, 7, 1, 15, 39, 63, 223, 17, 9, + 1, 3, 7, 1, 17, 29, 67, 103, 495, 383, + 1, 3, 3, 15, 31, 59, 75, 165, 51, 913, + 1, 3, 7, 9, 5, 27, 79, 219, 233, 37, + 1, 3, 5, 15, 1, 11, 15, 211, 417, 811, + 1, 3, 5, 3, 29, 27, 39, 137, 407, 231, + 1, 1, 3, 5, 29, 43, 125, 135, 109, 67, + 1, 1, 1, 5, 11, 39, 107, 159, 323, 381, + 1, 1, 1, 1, 9, 11, 33, 55, 169, 253, + 1, 3, 5, 5, 11, 53, 63, 101, 251, 897, + 1, 3, 7, 1, 25, 15, 83, 119, 53, 157, + 1, 3, 5, 13, 5, 5, 3, 195, 111, 451, + 1, 3, 1, 15, 11, 1, 19, 11, 307, 777, + 1, 3, 7, 11, 5, 5, 17, 231, 345, 981, + 1, 1, 3, 3, 1, 33, 83, 201, 57, 475, + 1, 3, 7, 7, 17, 13, 35, 175, 499, 809, + 1, 1, 5, 3, 3, 17, 103, 119, 499, 865, + 1, 1, 1, 11, 27, 25, 37, 121, 401, 11, + 1, 1, 1, 11, 9, 25, 25, 241, 403, 3, + 1, 1, 1, 1, 11, 1, 39, 163, 231, 573, + 1, 1, 1, 13, 13, 21, 75, 185, 99, 545, + 1, 1, 1, 15, 3, 63, 69, 11, 173, 315, + 1, 3, 5, 15, 11, 3, 95, 49, 123, 765, + 1, 1, 1, 15, 3, 63, 77, 31, 425, 711, + 1, 1, 7, 15, 1, 37, 119, 145, 489, 583, + 1, 3, 5, 15, 3, 49, 117, 211, 165, 323, + 1, 3, 7, 1, 27, 63, 77, 201, 225, 803, + 1, 1, 1, 11, 23, 35, 67, 21, 469, 357, + 1, 1, 7, 7, 9, 7, 25, 237, 237, 571, + 1, 1, 3, 15, 29, 5, 107, 109, 241, 47, + 1, 3, 5, 11, 27, 63, 29, 13, 203, 675, + 1, 1, 3, 9, 9, 11, 103, 179, 449, 263, + 1, 3, 5, 11, 29, 63, 53, 151, 259, 223, + 1, 1, 3, 7, 9, 25, 5, 197, 237, 163, + 1, 3, 7, 13, 5, 57, 67, 193, 147, 241, + 1, 1, 5, 15, 15, 33, 17, 67, 161, 341, + 1, 1, 3, 13, 17, 43, 21, 197, 441, 985, + 1, 3, 1, 5, 15, 33, 33, 193, 305, 829, + 1, 1, 1, 13, 19, 27, 71, 187, 477, 239, + 1, 1, 1, 9, 9, 17, 41, 177, 229, 983, + 1, 3, 5, 9, 15, 45, 97, 205, 43, 767, + 1, 1, 1, 9, 31, 31, 77, 159, 395, 809, + 1, 3, 3, 3, 29, 19, 73, 123, 165, 307, + 1, 3, 1, 7, 5, 11, 77, 227, 355, 403, + 1, 3, 5, 5, 25, 31, 1, 215, 451, 195, + 1, 3, 7, 15, 29, 37, 101, 241, 17, 633, + 1, 1, 5, 1, 11, 3, 107, 137, 489, 5, + 1, 1, 1, 7, 19, 19, 75, 85, 471, 355, + 1, 1, 3, 3, 9, 13, 113, 167, 13, 27, + 1, 3, 5, 11, 21, 3, 89, 205, 377, 307, + 1, 1, 1, 9, 31, 61, 65, 9, 391, 141, 867, + 1, 1, 1, 9, 19, 19, 61, 227, 241, 55, 161, + 1, 1, 1, 11, 1, 19, 7, 233, 463, 171, 1941, + 1, 1, 5, 7, 25, 13, 103, 75, 19, 1021, 1063, + 1, 1, 1, 15, 17, 17, 79, 63, 391, 403, 1221, + 1, 3, 3, 11, 29, 25, 29, 107, 335, 475, 963, + 1, 3, 5, 1, 31, 33, 49, 43, 155, 9, 1285, + 1, 1, 5, 5, 15, 47, 39, 161, 357, 863, 1039, + 1, 3, 7, 15, 1, 39, 47, 109, 427, 393, 1103, + 1, 1, 1, 9, 9, 29, 121, 233, 157, 99, 701, + 1, 1, 1, 7, 1, 29, 75, 121, 439, 109, 993, + 1, 1, 1, 9, 5, 1, 39, 59, 89, 157, 1865, + 1, 1, 5, 1, 3, 37, 89, 93, 143, 533, 175, + 1, 1, 3, 5, 7, 33, 35, 173, 159, 135, 241, + 1, 1, 1, 15, 17, 37, 79, 131, 43, 891, 229, + 1, 1, 1, 1, 1, 35, 121, 177, 397, 1017, 583, + 1, 1, 3, 15, 31, 21, 43, 67, 467, 923, 1473, + 1, 1, 1, 7, 1, 33, 77, 111, 125, 771, 1975, + 1, 3, 7, 13, 1, 51, 113, 139, 245, 573, 503, + 1, 3, 1, 9, 21, 49, 15, 157, 49, 483, 291, + 1, 1, 1, 1, 29, 35, 17, 65, 403, 485, 1603, + 1, 1, 1, 7, 19, 1, 37, 129, 203, 321, 1809, + 1, 3, 7, 15, 15, 9, 5, 77, 29, 485, 581, + 1, 1, 3, 5, 15, 49, 97, 105, 309, 875, 1581, + 1, 3, 5, 1, 5, 19, 63, 35, 165, 399, 1489, + 1, 3, 5, 3, 23, 5, 79, 137, 115, 599, 1127, + 1, 1, 7, 5, 3, 61, 27, 177, 257, 91, 841, + 1, 1, 3, 5, 9, 31, 91, 209, 409, 661, 159, + 1, 3, 1, 15, 23, 39, 23, 195, 245, 203, 947, + 1, 1, 3, 1, 15, 59, 67, 95, 155, 461, 147, + 1, 3, 7, 5, 23, 25, 87, 11, 51, 449, 1631, + 1, 1, 1, 1, 17, 57, 7, 197, 409, 609, 135, + 1, 1, 1, 9, 1, 61, 115, 113, 495, 895, 1595, + 1, 3, 7, 15, 9, 47, 121, 211, 379, 985, 1755, + 1, 3, 1, 3, 7, 57, 27, 231, 339, 325, 1023, + 1, 1, 1, 1, 19, 63, 63, 239, 31, 643, 373, + 1, 3, 1, 11, 19, 9, 7, 171, 21, 691, 215, + 1, 1, 5, 13, 11, 57, 39, 211, 241, 893, 555, + 1, 1, 7, 5, 29, 21, 45, 59, 509, 223, 491, + 1, 1, 7, 9, 15, 61, 97, 75, 127, 779, 839, + 1, 1, 7, 15, 17, 33, 75, 237, 191, 925, 681, + 1, 3, 5, 7, 27, 57, 123, 111, 101, 371, 1129, + 1, 3, 5, 5, 29, 45, 59, 127, 229, 967, 2027, + 1, 1, 1, 1, 17, 7, 23, 199, 241, 455, 135, + 1, 1, 7, 15, 27, 29, 105, 171, 337, 503, 1817, + 1, 1, 3, 7, 21, 35, 61, 71, 405, 647, 2045, + 1, 1, 1, 1, 1, 15, 65, 167, 501, 79, 737, + 1, 1, 5, 1, 3, 49, 27, 189, 341, 615, 1287, + 1, 1, 1, 9, 1, 7, 31, 159, 503, 327, 1613, + 1, 3, 3, 3, 3, 23, 99, 115, 323, 997, 987, + 1, 1, 1, 9, 19, 33, 93, 247, 509, 453, 891, + 1, 1, 3, 1, 13, 19, 35, 153, 161, 633, 445, + 1, 3, 5, 15, 31, 5, 87, 197, 183, 783, 1823, + 1, 1, 7, 5, 19, 63, 69, 221, 129, 231, 1195, + 1, 1, 5, 5, 13, 23, 19, 231, 245, 917, 379, + 1, 3, 1, 15, 19, 43, 27, 223, 171, 413, 125, + 1, 1, 1, 9, 1, 59, 21, 15, 509, 207, 589, + 1, 3, 5, 3, 19, 31, 113, 19, 23, 733, 499, + 1, 1, 7, 1, 19, 51, 101, 165, 47, 925, 1093, + 1, 3, 3, 9, 15, 21, 43, 243, 237, 461, 1361, + 1, 1, 1, 9, 17, 15, 75, 75, 113, 715, 1419, + 1, 1, 7, 13, 17, 1, 99, 15, 347, 721, 1405, + 1, 1, 7, 15, 7, 27, 23, 183, 39, 59, 571, + 1, 3, 5, 9, 7, 43, 35, 165, 463, 567, 859, + 1, 3, 3, 11, 15, 19, 17, 129, 311, 343, 15, + 1, 1, 1, 15, 31, 59, 63, 39, 347, 359, 105, + 1, 1, 1, 15, 5, 43, 87, 241, 109, 61, 685, + 1, 1, 7, 7, 9, 39, 121, 127, 369, 579, 853, + 1, 1, 1, 1, 17, 15, 15, 95, 325, 627, 299, + 1, 1, 3, 13, 31, 53, 85, 111, 289, 811, 1635, + 1, 3, 7, 1, 19, 29, 75, 185, 153, 573, 653, + 1, 3, 7, 1, 29, 31, 55, 91, 249, 247, 1015, + 1, 3, 5, 7, 1, 49, 113, 139, 257, 127, 307, + 1, 3, 5, 9, 15, 15, 123, 105, 105, 225, 1893, + 1, 3, 3, 1, 15, 5, 105, 249, 73, 709, 1557, + 1, 1, 1, 9, 17, 31, 113, 73, 65, 701, 1439, + 1, 3, 5, 15, 13, 21, 117, 131, 243, 859, 323, + 1, 1, 1, 9, 19, 15, 69, 149, 89, 681, 515, + 1, 1, 1, 5, 29, 13, 21, 97, 301, 27, 967, + 1, 1, 3, 3, 15, 45, 107, 227, 495, 769, 1935, + 1, 1, 1, 11, 5, 27, 41, 173, 261, 703, 1349, + 1, 3, 3, 3, 11, 35, 97, 43, 501, 563, 1331, + 1, 1, 1, 7, 1, 17, 87, 17, 429, 245, 1941, + 1, 1, 7, 15, 29, 13, 1, 175, 425, 233, 797, + 1, 1, 3, 11, 21, 57, 49, 49, 163, 685, 701, + 1, 3, 3, 7, 11, 45, 107, 111, 379, 703, 1403, + 1, 1, 7, 3, 21, 7, 117, 49, 469, 37, 775, + 1, 1, 5, 15, 31, 63, 101, 77, 507, 489, 1955, + 1, 3, 3, 11, 19, 21, 101, 255, 203, 673, 665, + 1, 3, 3, 15, 17, 47, 125, 187, 271, 899, 2003, + 1, 1, 7, 7, 1, 35, 13, 235, 5, 337, 905, + 1, 3, 1, 15, 1, 43, 1, 27, 37, 695, 1429, + 1, 3, 1, 11, 21, 27, 93, 161, 299, 665, 495, + 1, 3, 3, 15, 3, 1, 81, 111, 105, 547, 897, + 1, 3, 5, 1, 3, 53, 97, 253, 401, 827, 1467, + 1, 1, 1, 5, 19, 59, 105, 125, 271, 351, 719, + 1, 3, 5, 13, 7, 11, 91, 41, 441, 759, 1827, + 1, 3, 7, 11, 29, 61, 61, 23, 307, 863, 363, + 1, 1, 7, 1, 15, 35, 29, 133, 415, 473, 1737, + 1, 1, 1, 13, 7, 33, 35, 225, 117, 681, 1545, + 1, 1, 1, 3, 5, 41, 83, 247, 13, 373, 1091, + 1, 3, 1, 13, 25, 61, 71, 217, 233, 313, 547, + 1, 3, 1, 7, 3, 29, 3, 49, 93, 465, 15, + 1, 1, 1, 9, 17, 61, 99, 163, 129, 485, 1087, + 1, 1, 1, 9, 9, 33, 31, 163, 145, 649, 253, + 1, 1, 1, 1, 17, 63, 43, 235, 287, 111, 567, + 1, 3, 5, 13, 29, 7, 11, 69, 153, 127, 449, + 1, 1, 5, 9, 11, 21, 15, 189, 431, 493, 1219, + 1, 1, 1, 15, 19, 5, 47, 91, 399, 293, 1743, + 1, 3, 3, 11, 29, 53, 53, 225, 409, 303, 333, + 1, 1, 1, 15, 31, 31, 21, 81, 147, 287, 1753, + 1, 3, 5, 5, 5, 63, 35, 125, 41, 687, 1793, + 1, 1, 1, 9, 19, 59, 107, 219, 455, 971, 297, + 1, 1, 3, 5, 3, 51, 121, 31, 245, 105, 1311, + 1, 3, 1, 5, 5, 57, 75, 107, 161, 431, 1693, + 1, 3, 1, 3, 19, 53, 27, 31, 191, 565, 1015, + 1, 3, 5, 13, 9, 41, 35, 249, 287, 49, 123, + 1, 1, 5, 7, 27, 17, 21, 3, 151, 885, 1165, + 1, 1, 7, 1, 15, 17, 65, 139, 427, 339, 1171, + 1, 1, 1, 5, 23, 5, 9, 89, 321, 907, 391, + 1, 1, 7, 9, 15, 1, 77, 71, 87, 701, 917, + 1, 1, 7, 1, 17, 37, 115, 127, 469, 779, 1543, + 1, 3, 7, 3, 5, 61, 15, 37, 301, 951, 1437, + 1, 1, 1, 13, 9, 51, 127, 145, 229, 55, 1567, + 1, 3, 7, 15, 19, 47, 53, 153, 295, 47, 1337, + 1, 3, 3, 5, 11, 31, 29, 133, 327, 287, 507, + 1, 1, 7, 7, 25, 31, 37, 199, 25, 927, 1317, + 1, 1, 7, 9, 3, 39, 127, 167, 345, 467, 759, + 1, 1, 1, 1, 31, 21, 15, 101, 293, 787, 1025, + 1, 1, 5, 3, 11, 41, 105, 109, 149, 837, 1813, + 1, 1, 3, 5, 29, 13, 19, 97, 309, 901, 753, + 1, 1, 7, 1, 19, 17, 31, 39, 173, 361, 1177, + 1, 3, 3, 3, 3, 41, 81, 7, 341, 491, 43, + 1, 1, 7, 7, 31, 35, 29, 77, 11, 335, 1275, + 1, 3, 3, 15, 17, 45, 19, 63, 151, 849, 129, + 1, 1, 7, 5, 7, 13, 47, 73, 79, 31, 499, + 1, 3, 1, 11, 1, 41, 59, 151, 247, 115, 1295, + 1, 1, 1, 9, 31, 37, 73, 23, 295, 483, 179, + 1, 3, 1, 15, 13, 63, 81, 27, 169, 825, 2037, + 1, 3, 5, 15, 7, 11, 73, 1, 451, 101, 2039, + 1, 3, 5, 3, 13, 53, 31, 137, 173, 319, 1521, + 1, 3, 1, 3, 29, 1, 73, 227, 377, 337, 1189, + 1, 3, 3, 13, 27, 9, 31, 101, 229, 165, 1983, + 1, 3, 1, 13, 13, 19, 19, 111, 319, 421, 223, + 1, 1, 7, 15, 25, 37, 61, 55, 359, 255, 1955, + 1, 1, 5, 13, 17, 43, 49, 215, 383, 915, 51, + 1, 1, 3, 1, 3, 7, 13, 119, 155, 585, 967, + 1, 3, 1, 13, 1, 63, 125, 21, 103, 287, 457, + 1, 1, 7, 1, 31, 17, 125, 137, 345, 379, 1925, + 1, 1, 3, 5, 5, 25, 119, 153, 455, 271, 2023, + 1, 1, 7, 9, 9, 37, 115, 47, 5, 255, 917, + 1, 3, 5, 3, 31, 21, 75, 203, 489, 593, 1, + 1, 3, 7, 15, 19, 63, 123, 153, 135, 977, 1875, + 1, 1, 1, 1, 5, 59, 31, 25, 127, 209, 745, + 1, 1, 1, 1, 19, 45, 67, 159, 301, 199, 535, + 1, 1, 7, 1, 31, 17, 19, 225, 369, 125, 421, + 1, 3, 3, 11, 7, 59, 115, 197, 459, 469, 1055, + 1, 3, 1, 3, 27, 45, 35, 131, 349, 101, 411, + 1, 3, 7, 11, 9, 3, 67, 145, 299, 253, 1339, + 1, 3, 3, 11, 9, 37, 123, 229, 273, 269, 515, + 1, 3, 7, 15, 11, 25, 75, 5, 367, 217, 951, + 1, 1, 3, 7, 9, 23, 63, 237, 385, 159, 1273, + 1, 1, 5, 11, 23, 5, 55, 193, 109, 865, 663, + 1, 1, 7, 15, 1, 57, 17, 141, 51, 217, 1259, + 1, 1, 3, 3, 15, 7, 89, 233, 71, 329, 203, + 1, 3, 7, 11, 11, 1, 19, 155, 89, 437, 573, + 1, 3, 1, 9, 27, 61, 47, 109, 161, 913, 1681, + 1, 1, 7, 15, 1, 33, 19, 15, 23, 913, 989, + 1, 3, 1, 1, 25, 39, 119, 193, 13, 571, 157, + 1, 1, 7, 13, 9, 55, 59, 147, 361, 935, 515, + 1, 1, 1, 9, 7, 59, 67, 117, 71, 855, 1493, + 1, 3, 1, 3, 13, 19, 57, 141, 305, 275, 1079, + 1, 1, 1, 9, 17, 61, 33, 7, 43, 931, 781, + 1, 1, 3, 1, 11, 17, 21, 97, 295, 277, 1721, + 1, 3, 1, 13, 15, 43, 11, 241, 147, 391, 1641, + 1, 1, 1, 1, 1, 19, 37, 21, 255, 263, 1571, + 1, 1, 3, 3, 23, 59, 89, 17, 475, 303, 757, 543, + 1, 3, 3, 9, 11, 55, 35, 159, 139, 203, 1531, 1825, + 1, 1, 5, 3, 17, 53, 51, 241, 269, 949, 1373, 325, + 1, 3, 7, 7, 5, 29, 91, 149, 239, 193, 1951, 2675, + 1, 3, 5, 1, 27, 33, 69, 11, 51, 371, 833, 2685, + 1, 1, 1, 15, 1, 17, 35, 57, 171, 1007, 449, 367, + 1, 1, 1, 7, 25, 61, 73, 219, 379, 53, 589, 4065, + 1, 3, 5, 13, 21, 29, 45, 19, 163, 169, 147, 597, + 1, 1, 5, 11, 21, 27, 7, 17, 237, 591, 255, 1235, + 1, 1, 7, 7, 17, 41, 69, 237, 397, 173, 1229, 2341, + 1, 1, 3, 1, 1, 33, 125, 47, 11, 783, 1323, 2469, + 1, 3, 1, 11, 3, 39, 35, 133, 153, 55, 1171, 3165, + 1, 1, 5, 11, 27, 23, 103, 245, 375, 753, 477, 2165, + 1, 3, 1, 15, 15, 49, 127, 223, 387, 771, 1719, 1465, + 1, 1, 1, 9, 11, 9, 17, 185, 239, 899, 1273, 3961, + 1, 1, 3, 13, 11, 51, 73, 81, 389, 647, 1767, 1215, + 1, 3, 5, 15, 19, 9, 69, 35, 349, 977, 1603, 1435, + 1, 1, 1, 1, 19, 59, 123, 37, 41, 961, 181, 1275, + 1, 1, 1, 1, 31, 29, 37, 71, 205, 947, 115, 3017, + 1, 1, 7, 15, 5, 37, 101, 169, 221, 245, 687, 195, + 1, 1, 1, 1, 19, 9, 125, 157, 119, 283, 1721, 743, + 1, 1, 7, 3, 1, 7, 61, 71, 119, 257, 1227, 2893, + 1, 3, 3, 3, 25, 41, 25, 225, 31, 57, 925, 2139, + 0 +}; + + +/* + * coefficients of the free direction integers as given in + * "Monte Carlo Methods in Finance", by Peter Jäckel, section 8.3 + */ + +const unsigned long rk_sobol_Jdirections[] = { + 1, + 1, 1, + 1, 3, 7, + 1, 1, 5, + 1, 3, 1, 1, + 1, 1, 3, 7, + 1, 3, 3, 9, 9, + 1, 3, 7, 7, 21, + 1, 1, 5, 11, 27, + 1, 1, 7, 3, 29, + 1, 3, 7, 13, 3, + 1, 3, 5, 1, 15, + 1, 1, 1, 9, 23, 37, + 1, 1, 3, 13, 11, 7, + 1, 3, 3, 5, 19, 33, + 1, 1, 7, 13, 25, 5, + 1, 1, 1, 3, 13, 39, + 1, 3, 5, 11, 7, 11, + 1, 3, 1, 7, 3, 23, 79, + 1, 3, 1, 15, 17, 63, 13, + 1, 3, 3, 3, 25, 17, 115, + 1, 3, 7, 9, 31, 29, 17, + 1, 1, 3, 15, 29, 15, 41, + 1, 3, 1, 9, 5, 21, 119, + 1, 1, 5, 5, 1, 27, 33, + 1, 1, 3, 1, 23, 13, 75, + 1, 1, 7, 7, 19, 25, 105, + 1, 3, 5, 5, 21, 9, 7, + 1, 1, 1, 15, 5, 49, 59, + 1, 3, 5, 15, 17, 19, 21, + 1, 1, 7, 11, 13, 29, 3, + 0 +}; + +/* + * 0 terminated list of primitive polynomials to speed up initialization + * All polynomials up to degree 13 (ie. 1111 polynomials) + */ +static const unsigned long rk_sobol_primitive_polynomials[] = { + 0x1UL, 0x3UL, 0x7UL, 0xBUL, 0xDUL, 0x13UL, 0x19UL, 0x25UL, 0x29UL, + 0x2FUL, 0x37UL, 0x3BUL, 0x3DUL, 0x43UL, 0x5BUL, 0x61UL, 0x67UL, 0x6DUL, + 0x73UL, 0x83UL, 0x89UL, 0x8FUL, 0x91UL, 0x9DUL, 0xA7UL, 0xABUL, 0xB9UL, + 0xBFUL, 0xC1UL, 0xCBUL, 0xD3UL, 0xD5UL, 0xE5UL, 0xEFUL, 0xF1UL, 0xF7UL, + 0xFDUL, 0x11DUL, 0x12BUL, 0x12DUL, 0x14DUL, 0x15FUL, 0x163UL, 0x165UL, + 0x169UL, 0x171UL, 0x187UL, 0x18DUL, 0x1A9UL, 0x1C3UL, 0x1CFUL, 0x1E7UL, + 0x1F5UL, 0x211UL, 0x21BUL, 0x221UL, 0x22DUL, 0x233UL, 0x259UL, 0x25FUL, + 0x269UL, 0x26FUL, 0x277UL, 0x27DUL, 0x287UL, 0x295UL, 0x2A3UL, 0x2A5UL, + 0x2AFUL, 0x2B7UL, 0x2BDUL, 0x2CFUL, 0x2D1UL, 0x2DBUL, 0x2F5UL, 0x2F9UL, + 0x313UL, 0x315UL, 0x31FUL, 0x323UL, 0x331UL, 0x33BUL, 0x34FUL, 0x35BUL, + 0x361UL, 0x36BUL, 0x36DUL, 0x373UL, 0x37FUL, 0x385UL, 0x38FUL, 0x3B5UL, + 0x3B9UL, 0x3C7UL, 0x3CBUL, 0x3CDUL, 0x3D5UL, 0x3D9UL, 0x3E3UL, 0x3E9UL, + 0x3FBUL, 0x409UL, 0x41BUL, 0x427UL, 0x42DUL, 0x465UL, 0x46FUL, 0x481UL, + 0x48BUL, 0x4C5UL, 0x4D7UL, 0x4E7UL, 0x4F3UL, 0x4FFUL, 0x50DUL, 0x519UL, + 0x523UL, 0x531UL, 0x53DUL, 0x543UL, 0x557UL, 0x56BUL, 0x585UL, 0x58FUL, + 0x597UL, 0x5A1UL, 0x5C7UL, 0x5E5UL, 0x5F7UL, 0x5FBUL, 0x613UL, 0x615UL, + 0x625UL, 0x637UL, 0x643UL, 0x64FUL, 0x65BUL, 0x679UL, 0x67FUL, 0x689UL, + 0x6B5UL, 0x6C1UL, 0x6D3UL, 0x6DFUL, 0x6FDUL, 0x717UL, 0x71DUL, 0x721UL, + 0x739UL, 0x747UL, 0x74DUL, 0x755UL, 0x759UL, 0x763UL, 0x77DUL, 0x78DUL, + 0x793UL, 0x7B1UL, 0x7DBUL, 0x7F3UL, 0x7F9UL, 0x805UL, 0x817UL, 0x82BUL, + 0x82DUL, 0x847UL, 0x863UL, 0x865UL, 0x871UL, 0x87BUL, 0x88DUL, 0x895UL, + 0x89FUL, 0x8A9UL, 0x8B1UL, 0x8CFUL, 0x8D1UL, 0x8E1UL, 0x8E7UL, 0x8EBUL, + 0x8F5UL, 0x90DUL, 0x913UL, 0x925UL, 0x929UL, 0x93BUL, 0x93DUL, 0x945UL, + 0x949UL, 0x951UL, 0x95BUL, 0x973UL, 0x975UL, 0x97FUL, 0x983UL, 0x98FUL, + 0x9ABUL, 0x9ADUL, 0x9B9UL, 0x9C7UL, 0x9D9UL, 0x9E5UL, 0x9F7UL, 0xA01UL, + 0xA07UL, 0xA13UL, 0xA15UL, 0xA29UL, 0xA49UL, 0xA61UL, 0xA6DUL, 0xA79UL, + 0xA7FUL, 0xA85UL, 0xA91UL, 0xA9DUL, 0xAA7UL, 0xAABUL, 0xAB3UL, 0xAB5UL, + 0xAD5UL, 0xADFUL, 0xAE9UL, 0xAEFUL, 0xAF1UL, 0xAFBUL, 0xB03UL, 0xB09UL, + 0xB11UL, 0xB33UL, 0xB3FUL, 0xB41UL, 0xB4BUL, 0xB59UL, 0xB5FUL, 0xB65UL, + 0xB6FUL, 0xB7DUL, 0xB87UL, 0xB8BUL, 0xB93UL, 0xB95UL, 0xBAFUL, 0xBB7UL, + 0xBBDUL, 0xBC9UL, 0xBDBUL, 0xBDDUL, 0xBE7UL, 0xBEDUL, 0xC0BUL, 0xC0DUL, + 0xC19UL, 0xC1FUL, 0xC57UL, 0xC61UL, 0xC6BUL, 0xC73UL, 0xC85UL, 0xC89UL, + 0xC97UL, 0xC9BUL, 0xC9DUL, 0xCB3UL, 0xCBFUL, 0xCC7UL, 0xCCDUL, 0xCD3UL, + 0xCD5UL, 0xCE3UL, 0xCE9UL, 0xCF7UL, 0xD03UL, 0xD0FUL, 0xD1DUL, 0xD27UL, + 0xD2DUL, 0xD41UL, 0xD47UL, 0xD55UL, 0xD59UL, 0xD63UL, 0xD6FUL, 0xD71UL, + 0xD93UL, 0xD9FUL, 0xDA9UL, 0xDBBUL, 0xDBDUL, 0xDC9UL, 0xDD7UL, 0xDDBUL, + 0xDE1UL, 0xDE7UL, 0xDF5UL, 0xE05UL, 0xE1DUL, 0xE21UL, 0xE27UL, 0xE2BUL, + 0xE33UL, 0xE39UL, 0xE47UL, 0xE4BUL, 0xE55UL, 0xE5FUL, 0xE71UL, 0xE7BUL, + 0xE7DUL, 0xE81UL, 0xE93UL, 0xE9FUL, 0xEA3UL, 0xEBBUL, 0xECFUL, 0xEDDUL, + 0xEF3UL, 0xEF9UL, 0xF0BUL, 0xF19UL, 0xF31UL, 0xF37UL, 0xF5DUL, 0xF6BUL, + 0xF6DUL, 0xF75UL, 0xF83UL, 0xF91UL, 0xF97UL, 0xF9BUL, 0xFA7UL, 0xFADUL, + 0xFB5UL, 0xFCDUL, 0xFD3UL, 0xFE5UL, 0xFE9UL, 0x1053UL, 0x1069UL, + 0x107BUL, 0x107DUL, 0x1099UL, 0x10D1UL, 0x10EBUL, 0x1107UL, 0x111FUL, + 0x1123UL, 0x113BUL, 0x114FUL, 0x1157UL, 0x1161UL, 0x116BUL, 0x1185UL, + 0x11B3UL, 0x11D9UL, 0x11DFUL, 0x120DUL, 0x1237UL, 0x123DUL, 0x1267UL, + 0x1273UL, 0x127FUL, 0x12B9UL, 0x12C1UL, 0x12CBUL, 0x130FUL, 0x131DUL, + 0x1321UL, 0x1339UL, 0x133FUL, 0x134DUL, 0x1371UL, 0x1399UL, 0x13A3UL, + 0x13A9UL, 0x1407UL, 0x1431UL, 0x1437UL, 0x144FUL, 0x145DUL, 0x1467UL, + 0x1475UL, 0x14A7UL, 0x14ADUL, 0x14D3UL, 0x150FUL, 0x151DUL, 0x154DUL, + 0x1593UL, 0x15C5UL, 0x15D7UL, 0x15DDUL, 0x15EBUL, 0x1609UL, 0x1647UL, + 0x1655UL, 0x1659UL, 0x16A5UL, 0x16BDUL, 0x1715UL, 0x1719UL, 0x1743UL, + 0x1745UL, 0x1775UL, 0x1789UL, 0x17ADUL, 0x17B3UL, 0x17BFUL, 0x17C1UL, + 0x1857UL, 0x185DUL, 0x1891UL, 0x1897UL, 0x18B9UL, 0x18EFUL, 0x191BUL, + 0x1935UL, 0x1941UL, 0x1965UL, 0x197BUL, 0x198BUL, 0x19B1UL, 0x19BDUL, + 0x19C9UL, 0x19CFUL, 0x19E7UL, 0x1A1BUL, 0x1A2BUL, 0x1A33UL, 0x1A69UL, + 0x1A8BUL, 0x1AD1UL, 0x1AE1UL, 0x1AF5UL, 0x1B0BUL, 0x1B13UL, 0x1B1FUL, + 0x1B57UL, 0x1B91UL, 0x1BA7UL, 0x1BBFUL, 0x1BC1UL, 0x1BD3UL, 0x1C05UL, + 0x1C11UL, 0x1C17UL, 0x1C27UL, 0x1C4DUL, 0x1C87UL, 0x1C9FUL, 0x1CA5UL, + 0x1CBBUL, 0x1CC5UL, 0x1CC9UL, 0x1CCFUL, 0x1CF3UL, 0x1D07UL, 0x1D23UL, + 0x1D43UL, 0x1D51UL, 0x1D5BUL, 0x1D75UL, 0x1D85UL, 0x1D89UL, 0x1E15UL, + 0x1E19UL, 0x1E2FUL, 0x1E45UL, 0x1E51UL, 0x1E67UL, 0x1E73UL, 0x1E8FUL, + 0x1EE3UL, 0x1F11UL, 0x1F1BUL, 0x1F27UL, 0x1F71UL, 0x1F99UL, 0x1FBBUL, + 0x1FBDUL, 0x1FC9UL, 0x201BUL, 0x2027UL, 0x2035UL, 0x2053UL, 0x2065UL, + 0x206FUL, 0x208BUL, 0x208DUL, 0x209FUL, 0x20A5UL, 0x20AFUL, 0x20BBUL, + 0x20BDUL, 0x20C3UL, 0x20C9UL, 0x20E1UL, 0x20F3UL, 0x210DUL, 0x2115UL, + 0x2129UL, 0x212FUL, 0x213BUL, 0x2143UL, 0x2167UL, 0x216BUL, 0x2179UL, + 0x2189UL, 0x2197UL, 0x219DUL, 0x21BFUL, 0x21C1UL, 0x21C7UL, 0x21CDUL, + 0x21DFUL, 0x21E3UL, 0x21F1UL, 0x21FBUL, 0x2219UL, 0x2225UL, 0x2237UL, + 0x223DUL, 0x2243UL, 0x225BUL, 0x225DUL, 0x2279UL, 0x227FUL, 0x2289UL, + 0x2297UL, 0x229BUL, 0x22B3UL, 0x22BFUL, 0x22CDUL, 0x22EFUL, 0x22F7UL, + 0x22FBUL, 0x2305UL, 0x2327UL, 0x232BUL, 0x2347UL, 0x2355UL, 0x2359UL, + 0x236FUL, 0x2371UL, 0x237DUL, 0x2387UL, 0x238DUL, 0x2395UL, 0x23A3UL, + 0x23A9UL, 0x23B1UL, 0x23B7UL, 0x23BBUL, 0x23E1UL, 0x23EDUL, 0x23F9UL, + 0x240BUL, 0x2413UL, 0x241FUL, 0x2425UL, 0x2429UL, 0x243DUL, 0x2451UL, + 0x2457UL, 0x2461UL, 0x246DUL, 0x247FUL, 0x2483UL, 0x249BUL, 0x249DUL, + 0x24B5UL, 0x24BFUL, 0x24C1UL, 0x24C7UL, 0x24CBUL, 0x24E3UL, 0x2509UL, + 0x2517UL, 0x251DUL, 0x2521UL, 0x252DUL, 0x2539UL, 0x2553UL, 0x2555UL, + 0x2563UL, 0x2571UL, 0x2577UL, 0x2587UL, 0x258BUL, 0x2595UL, 0x2599UL, + 0x259FUL, 0x25AFUL, 0x25BDUL, 0x25C5UL, 0x25CFUL, 0x25D7UL, 0x25EBUL, + 0x2603UL, 0x2605UL, 0x2611UL, 0x262DUL, 0x263FUL, 0x264BUL, 0x2653UL, + 0x2659UL, 0x2669UL, 0x2677UL, 0x267BUL, 0x2687UL, 0x2693UL, 0x2699UL, + 0x26B1UL, 0x26B7UL, 0x26BDUL, 0x26C3UL, 0x26EBUL, 0x26F5UL, 0x2713UL, + 0x2729UL, 0x273BUL, 0x274FUL, 0x2757UL, 0x275DUL, 0x276BUL, 0x2773UL, + 0x2779UL, 0x2783UL, 0x2791UL, 0x27A1UL, 0x27B9UL, 0x27C7UL, 0x27CBUL, + 0x27DFUL, 0x27EFUL, 0x27F1UL, 0x2807UL, 0x2819UL, 0x281FUL, 0x2823UL, + 0x2831UL, 0x283BUL, 0x283DUL, 0x2845UL, 0x2867UL, 0x2875UL, 0x2885UL, + 0x28ABUL, 0x28ADUL, 0x28BFUL, 0x28CDUL, 0x28D5UL, 0x28DFUL, 0x28E3UL, + 0x28E9UL, 0x28FBUL, 0x2909UL, 0x290FUL, 0x2911UL, 0x291BUL, 0x292BUL, + 0x2935UL, 0x293FUL, 0x2941UL, 0x294BUL, 0x2955UL, 0x2977UL, 0x297DUL, + 0x2981UL, 0x2993UL, 0x299FUL, 0x29AFUL, 0x29B7UL, 0x29BDUL, 0x29C3UL, + 0x29D7UL, 0x29F3UL, 0x29F5UL, 0x2A03UL, 0x2A0FUL, 0x2A1DUL, 0x2A21UL, + 0x2A33UL, 0x2A35UL, 0x2A4DUL, 0x2A69UL, 0x2A6FUL, 0x2A71UL, 0x2A7BUL, + 0x2A7DUL, 0x2AA5UL, 0x2AA9UL, 0x2AB1UL, 0x2AC5UL, 0x2AD7UL, 0x2ADBUL, + 0x2AEBUL, 0x2AF3UL, 0x2B01UL, 0x2B15UL, 0x2B23UL, 0x2B25UL, 0x2B2FUL, + 0x2B37UL, 0x2B43UL, 0x2B49UL, 0x2B6DUL, 0x2B7FUL, 0x2B85UL, 0x2B97UL, + 0x2B9BUL, 0x2BADUL, 0x2BB3UL, 0x2BD9UL, 0x2BE5UL, 0x2BFDUL, 0x2C0FUL, + 0x2C21UL, 0x2C2BUL, 0x2C2DUL, 0x2C3FUL, 0x2C41UL, 0x2C4DUL, 0x2C71UL, + 0x2C8BUL, 0x2C8DUL, 0x2C95UL, 0x2CA3UL, 0x2CAFUL, 0x2CBDUL, 0x2CC5UL, + 0x2CD1UL, 0x2CD7UL, 0x2CE1UL, 0x2CE7UL, 0x2CEBUL, 0x2D0DUL, 0x2D19UL, + 0x2D29UL, 0x2D2FUL, 0x2D37UL, 0x2D3BUL, 0x2D45UL, 0x2D5BUL, 0x2D67UL, + 0x2D75UL, 0x2D89UL, 0x2D8FUL, 0x2DA7UL, 0x2DABUL, 0x2DB5UL, 0x2DE3UL, + 0x2DF1UL, 0x2DFDUL, 0x2E07UL, 0x2E13UL, 0x2E15UL, 0x2E29UL, 0x2E49UL, + 0x2E4FUL, 0x2E5BUL, 0x2E5DUL, 0x2E61UL, 0x2E6BUL, 0x2E8FUL, 0x2E91UL, + 0x2E97UL, 0x2E9DUL, 0x2EABUL, 0x2EB3UL, 0x2EB9UL, 0x2EDFUL, 0x2EFBUL, + 0x2EFDUL, 0x2F05UL, 0x2F09UL, 0x2F11UL, 0x2F17UL, 0x2F3FUL, 0x2F41UL, + 0x2F4BUL, 0x2F4DUL, 0x2F59UL, 0x2F5FUL, 0x2F65UL, 0x2F69UL, 0x2F95UL, + 0x2FA5UL, 0x2FAFUL, 0x2FB1UL, 0x2FCFUL, 0x2FDDUL, 0x2FE7UL, 0x2FEDUL, + 0x2FF5UL, 0x2FFFUL, 0x3007UL, 0x3015UL, 0x3019UL, 0x302FUL, 0x3049UL, + 0x304FUL, 0x3067UL, 0x3079UL, 0x307FUL, 0x3091UL, 0x30A1UL, 0x30B5UL, + 0x30BFUL, 0x30C1UL, 0x30D3UL, 0x30D9UL, 0x30E5UL, 0x30EFUL, 0x3105UL, + 0x310FUL, 0x3135UL, 0x3147UL, 0x314DUL, 0x315FUL, 0x3163UL, 0x3171UL, + 0x317BUL, 0x31A3UL, 0x31A9UL, 0x31B7UL, 0x31C5UL, 0x31C9UL, 0x31DBUL, + 0x31E1UL, 0x31EBUL, 0x31EDUL, 0x31F3UL, 0x31FFUL, 0x3209UL, 0x320FUL, + 0x321DUL, 0x3227UL, 0x3239UL, 0x324BUL, 0x3253UL, 0x3259UL, 0x3265UL, + 0x3281UL, 0x3293UL, 0x3299UL, 0x329FUL, 0x32A9UL, 0x32B7UL, 0x32BBUL, + 0x32C3UL, 0x32D7UL, 0x32DBUL, 0x32E7UL, 0x3307UL, 0x3315UL, 0x332FUL, + 0x3351UL, 0x335DUL, 0x3375UL, 0x3397UL, 0x339BUL, 0x33ABUL, 0x33B9UL, + 0x33C1UL, 0x33C7UL, 0x33D5UL, 0x33E3UL, 0x33E5UL, 0x33F7UL, 0x33FBUL, + 0x3409UL, 0x341BUL, 0x3427UL, 0x3441UL, 0x344DUL, 0x345FUL, 0x3469UL, + 0x3477UL, 0x347BUL, 0x3487UL, 0x3493UL, 0x3499UL, 0x34A5UL, 0x34BDUL, + 0x34C9UL, 0x34DBUL, 0x34E7UL, 0x34F9UL, 0x350DUL, 0x351FUL, 0x3525UL, + 0x3531UL, 0x3537UL, 0x3545UL, 0x354FUL, 0x355DUL, 0x356DUL, 0x3573UL, + 0x357FUL, 0x359DUL, 0x35A1UL, 0x35B9UL, 0x35CDUL, 0x35D5UL, 0x35D9UL, + 0x35E3UL, 0x35E9UL, 0x35EFUL, 0x3601UL, 0x360BUL, 0x361FUL, 0x3625UL, + 0x362FUL, 0x363BUL, 0x3649UL, 0x3651UL, 0x365BUL, 0x3673UL, 0x3675UL, + 0x3691UL, 0x369BUL, 0x369DUL, 0x36ADUL, 0x36CBUL, 0x36D3UL, 0x36D5UL, + 0x36E3UL, 0x36EFUL, 0x3705UL, 0x370FUL, 0x371BUL, 0x3721UL, 0x372DUL, + 0x3739UL, 0x3741UL, 0x3747UL, 0x3753UL, 0x3771UL, 0x3777UL, 0x378BUL, + 0x3795UL, 0x3799UL, 0x37A3UL, 0x37C5UL, 0x37CFUL, 0x37D1UL, 0x37D7UL, + 0x37DDUL, 0x37E1UL, 0x37F3UL, 0x3803UL, 0x3805UL, 0x3817UL, 0x381DUL, + 0x3827UL, 0x3833UL, 0x384BUL, 0x3859UL, 0x3869UL, 0x3871UL, 0x38A3UL, + 0x38B1UL, 0x38BBUL, 0x38C9UL, 0x38CFUL, 0x38E1UL, 0x38F3UL, 0x38F9UL, + 0x3901UL, 0x3907UL, 0x390BUL, 0x3913UL, 0x3931UL, 0x394FUL, 0x3967UL, + 0x396DUL, 0x3983UL, 0x3985UL, 0x3997UL, 0x39A1UL, 0x39A7UL, 0x39ADUL, + 0x39CBUL, 0x39CDUL, 0x39D3UL, 0x39EFUL, 0x39F7UL, 0x39FDUL, 0x3A07UL, + 0x3A29UL, 0x3A2FUL, 0x3A3DUL, 0x3A51UL, 0x3A5DUL, 0x3A61UL, 0x3A67UL, + 0x3A73UL, 0x3A75UL, 0x3A89UL, 0x3AB9UL, 0x3ABFUL, 0x3ACDUL, 0x3AD3UL, + 0x3AD5UL, 0x3ADFUL, 0x3AE5UL, 0x3AE9UL, 0x3AFBUL, 0x3B11UL, 0x3B2BUL, + 0x3B2DUL, 0x3B35UL, 0x3B3FUL, 0x3B53UL, 0x3B59UL, 0x3B63UL, 0x3B65UL, + 0x3B6FUL, 0x3B71UL, 0x3B77UL, 0x3B8BUL, 0x3B99UL, 0x3BA5UL, 0x3BA9UL, + 0x3BB7UL, 0x3BBBUL, 0x3BD1UL, 0x3BE7UL, 0x3BF3UL, 0x3BFFUL, 0x3C0DUL, + 0x3C13UL, 0x3C15UL, 0x3C1FUL, 0x3C23UL, 0x3C25UL, 0x3C3BUL, 0x3C4FUL, + 0x3C5DUL, 0x3C6DUL, 0x3C83UL, 0x3C8FUL, 0x3C9DUL, 0x3CA7UL, 0x3CABUL, + 0x3CB9UL, 0x3CC7UL, 0x3CE9UL, 0x3CFBUL, 0x3CFDUL, 0x3D03UL, 0x3D17UL, + 0x3D1BUL, 0x3D21UL, 0x3D2DUL, 0x3D33UL, 0x3D35UL, 0x3D41UL, 0x3D4DUL, + 0x3D65UL, 0x3D69UL, 0x3D7DUL, 0x3D81UL, 0x3D95UL, 0x3DB1UL, 0x3DB7UL, + 0x3DC3UL, 0x3DD1UL, 0x3DDBUL, 0x3DE7UL, 0x3DEBUL, 0x3DF9UL, 0x3E05UL, + 0x3E09UL, 0x3E0FUL, 0x3E1BUL, 0x3E2BUL, 0x3E3FUL, 0x3E41UL, 0x3E53UL, + 0x3E65UL, 0x3E69UL, 0x3E8BUL, 0x3EA3UL, 0x3EBDUL, 0x3EC5UL, 0x3ED7UL, + 0x3EDDUL, 0x3EE1UL, 0x3EF9UL, 0x3F0DUL, 0x3F19UL, 0x3F1FUL, 0x3F25UL, + 0x3F37UL, 0x3F3DUL, 0x3F43UL, 0x3F45UL, 0x3F49UL, 0x3F51UL, 0x3F57UL, + 0x3F61UL, 0x3F83UL, 0x3F89UL, 0x3F91UL, 0x3FABUL, 0x3FB5UL, 0x3FE3UL, + 0x3FF7UL, 0x3FFDUL, + 0UL +}; + +rk_sobol_error rk_sobol_init(size_t dimension, rk_sobol_state *s, + rk_state *rs_dir, const unsigned long *directions, + const unsigned long *polynomials) +{ + rk_state rs_dir_temp; + int j, l, degree = 0, last_degree = 0, ooord = 0; + size_t k, cdir = 0, cpol = 0; + unsigned long polynomial = 1, rev = 0, last = 0; + + if (dimension == 0) + return RK_SOBOL_EINVAL; + + if (polynomials == NULL) + polynomials = rk_sobol_primitive_polynomials; + + /* Allocate the structure */ + s->direction = NULL; s->numerator = NULL; + s->direction = malloc(sizeof(*(s->direction))*dimension*LONG_BIT); + s->numerator = malloc(sizeof(*(s->numerator))*dimension); + if (!s->direction | !s->numerator) + { + if (!s->direction) free(s->direction); + if (!s->numerator) free(s->numerator); + return RK_SOBOL_ENOMEM; + } + + /* Initialize directions */ + /* Degree 0 */ + for (j = degree; j < LONG_BIT; j++) + s->direction[j*dimension] = 1UL << (LONG_BIT-j-1); + + /* Skip unused first polynomial */ + if (polynomials[cpol]) + cpol++; + + /* Degree >0 */ + for (k = 1; k < dimension; k++) + { + unsigned long temp; + + /* Find a new primitive polynomial */ + if (polynomials[cpol]) + polynomial = polynomials[cpol++]; + else if (rev) + { + /* We are generating polynomials out of order: + use the reverse of the previous polynomial */ + last = polynomial; + polynomial = rev; + rev = 0; + } + else + { + if (last) + { + polynomial = last; + last = 0; + } + + /* Find a new primitive polynomial */ + while(1) + { + if (polynomial == ULONG_MAX) + { + /* Not enough polynomials */ + free(s->direction); + free(s->numerator); + return RK_SOBOL_EINVAL; + } + + polynomial += 2; + + if (ooord) + { + unsigned long copy = polynomial; + /* We are generating polynomials out of order: + check if the reverse was already checked */ + for (rev = 0; copy; copy >>= 1) + rev = (rev << 1) | (copy & 1); + if (ooord && rev < polynomial) + continue; + } + + if (rk_isprimitive(polynomial)) + break; + } + + if (rev == polynomial) + /* We are generating polynomials out of order: + the reverse is not different, discard it */ + rev = 0; + } + + /* Compute the degree */ + for (temp = polynomial >> 1, degree = 0; temp; degree++, temp >>= 1); + + for (j=0; jdirection[j*dimension+k] = m << (LONG_BIT-j-1); + } + + /* Scaled recursion for directions */ + for (j = degree; j < LONG_BIT; j++) + { + unsigned long effdir = s->direction[(j-degree)*dimension+k], + ptemp = polynomial >> 1; + effdir ^= (effdir >> degree); + for (l = degree-1; l >= 1; l--, ptemp >>= 1) + if (ptemp & 1) + effdir ^= s->direction[(j-l)*dimension+k]; + s->direction[j*dimension+k] = effdir; + } + + /* Can we generate polynomials out of order ? */ + if (!ooord && polynomials[cpol] == 0 && degree > last_degree + && (directions == NULL || directions[cdir] == 0)) + ooord = 0; + else + last_degree = degree; + } + + /* Initialize numerator */ + for (k=0; knumerator[k] = 0; + + s->dimension = dimension; + s->gcount = 0; + s->count = 0; + return RK_SOBOL_OK; +} + +void rk_sobol_reinit(rk_sobol_state *s) +{ + size_t k; + + /* Initialize numerator */ + for (k=0; kdimension; k++) + s->numerator[k] = 0; + + s->count = 0; + s->gcount = 0; +} + +void rk_sobol_randomshift(rk_sobol_state *s, rk_state *rs_num) +{ + rk_state rs_num_temp; + size_t k; + + if (rs_num == NULL) + { + rs_num = &rs_num_temp; + rk_randomseed(rs_num); + } + + /* Initialize numerator */ + for (k=0; kdimension; k++) + s->numerator[k] = rk_ulong(rs_num); +} + +rk_sobol_error rk_sobol_copy(rk_sobol_state *copy, rk_sobol_state *orig) +{ + size_t k; + + /* Allocate the structure */ + copy->direction = NULL; copy->numerator = NULL; + copy->direction = malloc(sizeof(*(copy->direction))*orig->dimension*LONG_BIT); + copy->numerator = malloc(sizeof(*(copy->numerator))*orig->dimension); + if (!copy->direction | !copy->numerator) + { + if (!copy->direction) free(copy->direction); + if (!copy->numerator) free(copy->numerator); + return RK_SOBOL_ENOMEM; + } + + /* Initialize numerator */ + for (k=0; kdimension; k++) + copy->numerator[k] = orig->numerator[k]; + for (k=0; k<(orig->dimension*LONG_BIT); k++) + copy->direction[k] = orig->direction[k]; + + copy->count = orig->count; + copy->gcount = orig->gcount; + copy->dimension = orig->dimension; + + return RK_SOBOL_OK; +} + +rk_sobol_error rk_sobol_double(rk_sobol_state *s, double *x) +{ + int j; + size_t k; + unsigned long im; + const double inverse_denominator=1.0/(ULONG_MAX+1.0); + + if (s->count == ULONG_MAX) + j = 0; + else + for (im = s->count, j=0; im & 1; j++, im >>= 1); + s->count++; + + for (k=0; kdimension; k++) + { + s->numerator[k] ^= s->direction[j*s->dimension+k]; + x[k] = s->numerator[k]*inverse_denominator; + } + + if ((s->gcount++) == ULONG_MAX) return RK_SOBOL_EXHAUST; + return RK_SOBOL_OK; +} + +void rk_sobol_setcount(rk_sobol_state *s, unsigned long count) +{ + s->count = count; +} + +void rk_sobol_free(rk_sobol_state *s) +{ + free(s->direction); + free(s->numerator); +} + +double inverse_normal(double p) +{ + double q, t, x; + + /* Peter J. Acklam constants for the rational approximation */ + const double a[6] = + { + -3.969683028665376e+01, 2.209460984245205e+02, + -2.759285104469687e+02, 1.383577518672690e+02, + -3.066479806614716e+01, 2.506628277459239e+00 + }; + const double b[5] = + { + -5.447609879822406e+01, 1.615858368580409e+02, + -1.556989798598866e+02, 6.680131188771972e+01, + -1.328068155288572e+01 + }; + const double c[6] = + { + -7.784894002430293e-03, -3.223964580411365e-01, + -2.400758277161838e+00, -2.549732539343734e+00, + 4.374664141464968e+00, 2.938163982698783e+00 + }; + const double d[4] = + { + 7.784695709041462e-03, 3.224671290700398e-01, + 2.445134137142996e+00, 3.754408661907416e+00 + }; + + if (p <= 0) + return -HUGE_VAL; + else if (p >= 1) + return HUGE_VAL; + + q = p<0.5 ? p : 1-p; + if (q > 0.02425) + { + /* Rational approximation for central region */ + x = q-0.5; + t = x*x; + x = x*(((((a[0]*t+a[1])*t+a[2])*t+a[3])*t+a[4])*t+a[5]) + /(((((b[0]*t+b[1])*t+b[2])*t+b[3])*t+b[4])*t+1); + } + else + { + /* Rational approximation for tail region */ + t = sqrt(-2*log(q)); + x = (((((c[0]*t+c[1])*t+c[2])*t+c[3])*t+c[4])*t+c[5]) + /((((d[0]*t+d[1])*t+d[2])*t+d[3])*t+1); + } + + /* If we have erfc, improve the precision */ +#ifndef WIN32 + /* Halley's rational method */ + t = (erfc(-x*M_SQRT1_2)/2 - q) * RK_SOBOL_M_SQRT2PI * exp(x*x/2); + x -= t/(1 + x*t/2); +#endif + + return p>0.5 ? -x : x; +} + +rk_sobol_error rk_sobol_gauss(rk_sobol_state *s, double *x) +{ + size_t k; + rk_sobol_error rc = rk_sobol_double(s, x); + + for (k=0; kdimension; k++) + x[k] = inverse_normal(x[k]); + + return rc; +} diff --git a/tardis/montecarlo/src/randomkit/rk_sobol.h b/tardis/montecarlo/src/randomkit/rk_sobol.h new file mode 100644 index 00000000000..aca71bebe1e --- /dev/null +++ b/tardis/montecarlo/src/randomkit/rk_sobol.h @@ -0,0 +1,173 @@ +/* Random kit 1.6 */ + +/* + * Copyright (c) 2004-2006, Jean-Sebastien Roy (js@jeannot.org) + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. + * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY + * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, + * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE + * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +/* @(#) $Jeannot: rk_sobol.h,v 1.7 2006/02/19 13:48:34 js Exp $ */ + +/* + * Typical use: + * + * int dimension = 2; + * rk_sobol_state s; + * rk_sobol_error rc; + * double x[dimension], y[dimension]; + * + * // Init + * if (rc = rk_sobol_init(dimension, &s, NULL, NULL, NULL)) + * { + * fprintf(stderr, "%s\n", rk_sobol_strerror[rc]); + * abort(); + * } + * + * // Draw uniform quasirandom doubles + * if (rc = rk_sobol_double(&s, x)) + * { + * fprintf(stderr, "%s\n", rk_sobol_strerror[rc]); + * abort(); + * } + * + * // Draw gaussian quasirandom doubles + * if (rc = rk_sobol_gauss(&s, y)) + * { + * fprintf(stderr, "%s\n", rk_sobol_strerror[rc]); + * abort(); + * } + * + * // Free allocated memory + * rk_sobol_free(&s); + */ + + +#ifndef _RK_SOBOL_ +#define _RK_SOBOL_ + +#include "rk_mt.h" + +typedef enum { + RK_SOBOL_OK = 0, /* No error */ + RK_SOBOL_EINVAL = 1, /* Invalid dimension (<= 0 or too large) */ + RK_SOBOL_EXHAUST = 2, /* Too many number generated */ + RK_SOBOL_ENOMEM = 3, /* Not enough memory */ + RK_SOBOL_ERR_MAX = 4 +} rk_sobol_error; + +/* error strings */ +extern char *rk_sobol_strerror[]; + +typedef struct +{ + size_t dimension; + unsigned long *direction; + unsigned long *numerator; + unsigned long count; + unsigned long gcount; +} rk_sobol_state; + +#ifdef __cplusplus +extern "C" { +#endif + +/* Sobol directions initializations (zero terminated lists) */ + +/* + * Sobol/Levitan coefficients of the free direction integers as given + * by Bratley, P., Fox, B.L. (1988) + * Defined up to dimension 40. + */ +extern const unsigned long rk_sobol_SLdirections[]; + +/* + * Lemieux coefficients of the free direction integers as given + * in QuantLib by Christiane Lemieux, private communication, September 2004 + * Defined up to dimension 360. + */ +extern const unsigned long rk_sobol_Ldirections[]; + +/* + * Peter Jäckel coefficients of the free direction integers as given + * in "Monte Carlo Methods in Finance", by Peter Jäckel, section 8.3 + * Defined up to dimension 32. + */ +extern const unsigned long rk_sobol_Jdirections[]; + +/* + * Initialize a sobol quasirandom number generator. + * 1 <= dimension <= the number of primitive polylonimals of degree < LONG_BIT + * If directions == NULL (or more directions than provided are required), + * the directions are picked at random using rs_dir. + * If rs_dir == NULL, it is initialized using rk_randomseed. + * polynomials is a zero terminated list of primitive polynomials to use if + * it is != NULL to speed up initialization for dimension > 1024. + */ +extern rk_sobol_error rk_sobol_init(size_t dimension, rk_sobol_state *s, + rk_state *rs_dir, const unsigned long *directions, + const unsigned long *polynomials); + +/* + * Reinitialize the random generator with same directions. + */ +extern void rk_sobol_reinit(rk_sobol_state *s); + +/* + * You can change the starting rank in the sequence by changing s->count. + */ +extern void rk_sobol_setcount(rk_sobol_state *s, unsigned long count); + +/* + * XOR the numerators at random using rs_num. + * To be used once, after (re-)initialization. + * Useful for randomized quasi monte carlo. + * If rs_num == NULL, it is initialized using rk_randomseed. + */ +extern void rk_sobol_randomshift(rk_sobol_state *s, rk_state *rs_num); + +/* + * Copy a sobol generator. + * Can be used to avoid the time consuming initialization. + */ +extern rk_sobol_error rk_sobol_copy(rk_sobol_state *copy, rk_sobol_state *orig); + +/* + * Free the memory allocated by rk_sobol_init + */ +extern void rk_sobol_free(rk_sobol_state *s); + +/* + * return a vector of dimension quasirandom uniform deviates between 0 and 1 + */ +extern rk_sobol_error rk_sobol_double(rk_sobol_state *s, double *x); + +/* + * return a vector of dimension quasirandom gaussian deviates + * with variance unity and zero mean. + * On Windows, the standard function erfc is missing, which results in + * lower precision (9 digits instead of full precision). + */ +extern rk_sobol_error rk_sobol_gauss(rk_sobol_state *s, double *x); + +#ifdef __cplusplus +} +#endif + +#endif /* _RK_SOBOL_ */ diff --git a/tardis/montecarlo/src/rpacket.c b/tardis/montecarlo/src/rpacket.c new file mode 100644 index 00000000000..cb52cddc698 --- /dev/null +++ b/tardis/montecarlo/src/rpacket.c @@ -0,0 +1,55 @@ +#include +#include +#include "rpacket.h" +#include "storage.h" + +extern tardis_error_t line_search (const double *nu, double nu_insert, + int64_t number_of_lines, int64_t * result); + +tardis_error_t +rpacket_init (rpacket_t * packet, storage_model_t * storage, int packet_index, + int virtual_packet_flag, double * chi_bf_tmp_partial) +{ + int64_t current_line_id; + tardis_error_t ret_val = TARDIS_ERROR_OK; + double current_nu = storage->packet_nus[packet_index]; + double current_energy = storage->packet_energies[packet_index]; + double current_mu = storage->packet_mus[packet_index]; + double comov_current_nu = current_nu; + int current_shell_id = 0; + double current_r = storage->r_inner[0]; + double beta = current_r * storage->inverse_time_explosion * INVERSE_C; + + if (storage->full_relativity) + { + current_nu = current_nu * (1 + beta * current_mu) / sqrt(1 - beta * beta); + current_energy = current_energy * (1 + beta * current_mu) / sqrt(1 - beta * beta); + current_mu = (current_mu + beta) / (1 + beta * current_mu); + } + else + { + current_nu = current_nu / (1 - beta * current_mu); + current_energy = current_energy / (1 - beta * current_mu); + } + if ((ret_val = + line_search (storage->line_list_nu, comov_current_nu, + storage->no_of_lines, + ¤t_line_id)) != TARDIS_ERROR_OK) + { + return ret_val; + } + bool last_line = (current_line_id == storage->no_of_lines); + rpacket_set_nu (packet, current_nu); + rpacket_set_mu (packet, current_mu); + rpacket_set_energy (packet, current_energy); + rpacket_set_r (packet, current_r); + rpacket_set_current_shell_id (packet, current_shell_id); + rpacket_set_next_line_id (packet, current_line_id); + rpacket_set_last_line (packet, last_line); + rpacket_set_close_line (packet, false); + 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; +} diff --git a/tardis/montecarlo/src/rpacket.h b/tardis/montecarlo/src/rpacket.h new file mode 100644 index 00000000000..2a30546a837 --- /dev/null +++ b/tardis/montecarlo/src/rpacket.h @@ -0,0 +1,330 @@ +#ifndef TARDIS_RPACKET_H +#define TARDIS_RPACKET_H + +#include +#include +#include +#include "randomkit/randomkit.h" +#include "status.h" +#include "storage.h" + +#define MISS_DISTANCE 1e99 +#define C 29979245800.0 +#define INVERSE_C 3.33564095198152e-11 +#define H 6.6260755e-27 // erg * s, converted to CGS units from the NIST Constant Index +#define KB 1.3806488e-16 // erg / K converted to CGS units from the NIST Constant Index + +/** + * @brief A photon packet. + */ +typedef struct RPacket +{ + double nu; /**< Frequency of the packet in Hz. */ + double mu; /**< Cosine of the angle of the packet. */ + double energy; /**< Energy of the packet in erg. */ + double r; /**< Distance from center in cm. */ + double tau_event; + double nu_line; + int64_t current_shell_id; /**< ID of the current shell. */ + int64_t next_line_id; /**< The index of the next line that the packet will encounter. */ + /** + * @brief The packet has a nu red-ward of the last line. + * It will not encounter any lines anymore. + */ + int64_t last_line; + /** + * @brief The packet just encountered a line that is very close to the next line. + * The next iteration will automatically make an interaction with the next line + * (avoiding numerical problems). + */ + int64_t close_line; + /** + * @brief packet is a virtual packet and will ignore any d_line or d_electron checks. + * It now whenever a d_line is calculated only adds the tau_line to an + * internal float. + */ + int64_t current_continuum_id; /* Packet can interact with bf-continua with an index equal or bigger than this */ + int64_t virtual_packet_flag; + int64_t virtual_packet; + double d_line; /**< Distance to electron event. */ + double d_electron; /**< Distance to line event. */ + double d_boundary; /**< Distance to shell boundary. */ + double d_cont; /**< Distance to continuum event */ + int64_t next_shell_id; /**< ID of the next shell packet visits. */ + rpacket_status_t status; /**< Packet status (in process, emitted or reabsorbed). */ + int64_t id; + double chi_th; /**< Opacity due to electron scattering */ + double chi_cont; /**< Opacity due to continuum processes */ + double chi_ff; /**< Opacity due to free-free processes */ + double chi_bf; /**< Opacity due to bound-free processes */ + 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) +{ + return packet->nu; +} + +static inline void rpacket_set_nu (rpacket_t * packet, double nu) +{ + packet->nu = nu; +} + +static inline double rpacket_get_mu (const rpacket_t * packet) +{ + return packet->mu; +} + +static inline void rpacket_set_mu (rpacket_t * packet, double mu) +{ + packet->mu = mu; +} + +static inline double rpacket_get_energy (const rpacket_t * packet) +{ + return packet->energy; +} + +static inline void rpacket_set_energy (rpacket_t * packet, double energy) +{ + packet->energy = energy; +} + +static inline double rpacket_get_r (const rpacket_t * packet) +{ + return packet->r; +} + +static inline void rpacket_set_r (rpacket_t * packet, double r) +{ + packet->r = r; +} + +static inline double rpacket_get_tau_event (const rpacket_t * packet) +{ + return packet->tau_event; +} + +static inline void rpacket_set_tau_event (rpacket_t * packet, double tau_event) +{ + packet->tau_event = tau_event; +} + +static inline double rpacket_get_nu_line (const rpacket_t * packet) +{ + return packet->nu_line; +} + +static inline void rpacket_set_nu_line (rpacket_t * packet, double nu_line) +{ + packet->nu_line = nu_line; +} + +static inline unsigned int rpacket_get_current_shell_id (const rpacket_t * packet) +{ + return packet->current_shell_id; +} + +static inline void rpacket_set_current_shell_id (rpacket_t * packet, + unsigned int current_shell_id) +{ + packet->current_shell_id = current_shell_id; +} + +static inline unsigned int rpacket_get_next_line_id (const rpacket_t * packet) +{ + return packet->next_line_id; +} + +static inline void rpacket_set_next_line_id (rpacket_t * packet, + unsigned int next_line_id) +{ + packet->next_line_id = next_line_id; +} + +static inline bool rpacket_get_last_line (const rpacket_t * packet) +{ + return packet->last_line; +} + +static inline void rpacket_set_last_line (rpacket_t * packet, bool last_line) +{ + packet->last_line = last_line; +} + +static inline bool rpacket_get_close_line (const rpacket_t * packet) +{ + return packet->close_line; +} + +static inline void rpacket_set_close_line (rpacket_t * packet, bool close_line) +{ + packet->close_line = close_line; +} + +static inline int rpacket_get_virtual_packet_flag (const rpacket_t * packet) +{ + return packet->virtual_packet_flag; +} + +static inline void rpacket_set_virtual_packet_flag (rpacket_t * packet, + int virtual_packet_flag) +{ + packet->virtual_packet_flag = virtual_packet_flag; +} + +static inline int rpacket_get_virtual_packet (const rpacket_t * packet) +{ + return packet->virtual_packet; +} + +static inline void rpacket_set_virtual_packet (rpacket_t * packet, + int virtual_packet) +{ + packet->virtual_packet = virtual_packet; +} + +static inline double rpacket_get_d_boundary (const rpacket_t * packet) +{ + return packet->d_boundary; +} + +static inline void rpacket_set_d_boundary (rpacket_t * packet, double d_boundary) +{ + packet->d_boundary = d_boundary; +} + +static inline double rpacket_get_d_electron (const rpacket_t * packet) +{ + return packet->d_electron; +} + +static inline void rpacket_set_d_electron (rpacket_t * packet, double d_electron) +{ + packet->d_electron = d_electron; +} + +static inline double rpacket_get_d_line (const rpacket_t * packet) +{ + return packet->d_line; +} + +static inline void rpacket_set_d_line (rpacket_t * packet, double d_line) +{ + packet->d_line = d_line; +} + +static inline int rpacket_get_next_shell_id (const rpacket_t * packet) +{ + return packet->next_shell_id; +} + +static inline void rpacket_set_next_shell_id (rpacket_t * packet, int next_shell_id) +{ + packet->next_shell_id = next_shell_id; +} + +static inline rpacket_status_t rpacket_get_status (const rpacket_t * packet) +{ + return packet->status; +} + +static inline void rpacket_set_status (rpacket_t * packet, rpacket_status_t status) +{ + packet->status = status; +} + +static inline int rpacket_get_id (const rpacket_t * packet) +{ + return packet->id; +} + +static inline void rpacket_set_id (rpacket_t * packet, int id) +{ + packet->id = id; +} + +static inline void rpacket_reset_tau_event (rpacket_t * packet, rk_state *mt_state) +{ + rpacket_set_tau_event (packet, -log (rk_double (mt_state))); +} + +tardis_error_t rpacket_init (rpacket_t * packet, storage_model_t * storage, + int packet_index, int virtual_packet_flag, double * chi_bf_tmp_partial); + +/* New getter and setter methods for continuum implementation */ + +static inline void rpacket_set_d_continuum (rpacket_t * packet, double d_continuum) +{ + packet->d_cont = d_continuum; +} + +static inline double rpacket_get_d_continuum (const rpacket_t * packet) +{ + return packet->d_cont; +} + +static inline void rpacket_set_chi_electron (rpacket_t * packet, double chi_electron) +{ + packet->chi_th = chi_electron; +} + +static inline double rpacket_get_chi_electron (const rpacket_t * packet) +{ + return packet->chi_th; +} + +static inline void rpacket_set_chi_continuum (rpacket_t * packet, double chi_continuum) +{ + packet->chi_cont = chi_continuum; +} + +static inline double rpacket_get_chi_continuum (const rpacket_t * packet) +{ + return packet->chi_cont; +} + +static inline void rpacket_set_chi_freefree (rpacket_t * packet, double chi_freefree) +{ + packet->chi_ff = chi_freefree; +} + +static inline double rpacket_get_chi_freefree (const rpacket_t * packet) +{ + return packet->chi_ff; +} + +static inline void rpacket_set_chi_boundfree (rpacket_t * packet, double chi_boundfree) +{ + packet->chi_bf = chi_boundfree; +} + +static inline double rpacket_get_chi_boundfree (const rpacket_t * packet) +{ + return packet->chi_bf; +} + +static inline unsigned int rpacket_get_current_continuum_id (const rpacket_t * packet) +{ + return packet->current_continuum_id; +} + +static inline void rpacket_set_current_continuum_id (rpacket_t * packet, unsigned int current_continuum_id) +{ + packet->current_continuum_id = current_continuum_id; +} + +static inline void rpacket_set_macro_atom_activation_level (rpacket_t * packet, unsigned int activation_level) +{ + packet->macro_atom_activation_level = activation_level; +} + +static inline unsigned int rpacket_get_macro_atom_activation_level (const rpacket_t * packet) +{ + return packet->macro_atom_activation_level; +} + +#endif // TARDIS_RPACKET_H diff --git a/tardis/montecarlo/src/status.h b/tardis/montecarlo/src/status.h new file mode 100644 index 00000000000..17a01720b79 --- /dev/null +++ b/tardis/montecarlo/src/status.h @@ -0,0 +1,38 @@ +#ifndef TARDIS_STATUS_H +#define TARDIS_STATUS_H + +typedef enum +{ + TARDIS_ERROR_OK = 0, + TARDIS_ERROR_BOUNDS_ERROR = 1, + TARDIS_ERROR_COMOV_NU_LESS_THAN_NU_LINE = 2 +} tardis_error_t; + +typedef enum +{ + TARDIS_PACKET_STATUS_IN_PROCESS = 0, + TARDIS_PACKET_STATUS_EMITTED = 1, + TARDIS_PACKET_STATUS_REABSORBED = 2 +} rpacket_status_t; + +typedef enum +{ + CONTINUUM_OFF = 0, + CONTINUUM_ON = 1 +} ContinuumProcessesStatus; + +typedef enum +{ + BB_EMISSION = -1, + BF_EMISSION = -2, + FF_EMISSION = -3, + ADIABATIC_COOLING = -4 +} emission_type; + +typedef enum +{ + LIN_INTERPOLATION = 0, + HYDROGENIC = 1 +} bound_free_treatment; + +#endif // TARDIS_STATUS_H diff --git a/tardis/montecarlo/src/storage.h b/tardis/montecarlo/src/storage.h new file mode 100644 index 00000000000..a408a9f7bc2 --- /dev/null +++ b/tardis/montecarlo/src/storage.h @@ -0,0 +1,102 @@ +#ifndef TARDIS_STORAGE_H +#define TARDIS_STORAGE_H + +#include + +#include "status.h" + +typedef struct photo_xsect_1level +{ + double * nu; + double * x_sect; + int64_t no_of_points; +} photo_xsect_1level; + +typedef struct StorageModel +{ + double *packet_nus; + double *packet_mus; + double *packet_energies; + double *output_nus; + double *output_energies; + double *last_interaction_in_nu; + int64_t *last_line_interaction_in_id; + int64_t *last_line_interaction_out_id; + int64_t *last_line_interaction_shell_id; + int64_t *last_interaction_type; + int64_t *last_interaction_out_type; + int64_t no_of_packets; + int64_t no_of_shells; + int64_t no_of_shells_i; + double *r_inner; + double *r_outer; + double *r_inner_i; + double *r_outer_i; + double *v_inner; + double time_explosion; + double inverse_time_explosion; + double *electron_densities; + double *electron_densities_i; + double *inverse_electron_densities; + double *line_list_nu; + double *continuum_list_nu; + double *line_lists_tau_sobolevs; + double *line_lists_tau_sobolevs_i; + int64_t line_lists_tau_sobolevs_nd; + double *line_lists_j_blues; + int64_t line_lists_j_blues_nd; + double *line_lists_Edotlu; + int64_t no_of_lines; + int64_t no_of_edges; + int64_t line_interaction_id; + double *transition_probabilities; + int64_t transition_probabilities_nd; + int64_t *line2macro_level_upper; + int64_t *macro_block_references; + int64_t *transition_type; + int64_t *destination_level_id; + int64_t *transition_line_id; + double *js; + double *nubars; + double spectrum_start_nu; + double spectrum_delta_nu; + double spectrum_end_nu; + double spectrum_virt_start_nu; + double spectrum_virt_end_nu; + double *spectrum_virt_nu; + double sigma_thomson; + double inverse_sigma_thomson; + double inner_boundary_albedo; + int64_t reflective_inner_boundary; + int64_t current_packet_id; + photo_xsect_1level **photo_xsect; + double *chi_ff_factor; + double *t_electrons; + double *l_pop; + double *l_pop_r; + ContinuumProcessesStatus cont_status; + bound_free_treatment bf_treatment; + double *virt_packet_nus; + double *virt_packet_energies; + double *virt_packet_last_interaction_in_nu; + int64_t *virt_packet_last_interaction_type; + int64_t *virt_packet_last_line_interaction_in_id; + int64_t *virt_packet_last_line_interaction_out_id; + int64_t virt_packet_count; + int64_t virt_array_size; + int64_t kpacket2macro_level; + int64_t *cont_edge2macro_level; + double *photo_ion_estimator; + double *stim_recomb_estimator; + int64_t *photo_ion_estimator_statistics; + double *bf_heating_estimator; + double *ff_heating_estimator; + double *stim_recomb_cooling_estimator; + int full_relativity; + double survival_probability; + double tau_russian; + double *tau_bias; + int enable_biasing; +} storage_model_t; + +#endif // TARDIS_STORAGE_H