Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
SWIFT
SWIFTsim
Commits
36b7e032
Commit
36b7e032
authored
Apr 15, 2019
by
Matthieu Schaller
Browse files
Make the code compile with the feedback functions moved to the separate directory.
parent
367d8af5
Changes
15
Expand all
Hide whitespace changes
Inline
Side-by-side
configure.ac
View file @
36b7e032
...
...
@@ -1329,7 +1329,7 @@ case "$with_subgrid" in
with_subgrid_entropy_floor=none
with_subgrid_stars=GEAR
with_subgrid_star_formation=GEAR
with_subgrid_feedback=
thermal
with_subgrid_feedback=
none
;;
EAGLE)
with_subgrid_cooling=EAGLE
...
...
@@ -1338,7 +1338,7 @@ case "$with_subgrid" in
with_subgrid_entropy_floor=EAGLE
with_subgrid_stars=EAGLE
with_subgrid_star_formation=EAGLE
with_subgrid_feedback=
none
with_subgrid_feedback=
EAGLE
;;
*)
AC_MSG_ERROR([Unknown subgrid choice: $with_subgrid])
...
...
@@ -1727,11 +1727,8 @@ case "$with_stars" in
GEAR)
AC_DEFINE([STARS_GEAR], [1], [GEAR stellar model])
;;
EAGLE)
AC_DEFINE([STARS_EAGLE], [1], [EAGLE stellar model])
;;
none)
AC_DEFINE([STARS_NONE], [1], [
None
stellar model])
AC_DEFINE([STARS_NONE], [1], [
Basic
stellar model])
;;
*)
...
...
@@ -1742,7 +1739,7 @@ esac
# Feedback model
AC_ARG_WITH([feedback],
[AS_HELP_STRING([--with-feedback=<model>],
[Feedback model to use @<:@none,
thermal
, debug default: none@:>@]
[Feedback model to use @<:@none,
EAGLE
, debug default: none@:>@]
)],
[with_feedback="$withval"],
[with_feedback="none"]
...
...
@@ -1757,10 +1754,11 @@ if test "$with_subgrid" != "none"; then
fi
case "$with_feedback" in
thermal
)
AC_DEFINE([FEEDBACK_
THERMAL
], [1], [
Thermal Blastwave
])
EAGLE
)
AC_DEFINE([FEEDBACK_
EAGLE
], [1], [
EAGLE stellar feedback and evolution model
])
;;
none)
AC_DEFINE([FEEDBACK_NONE], [1], [No feedback])
;;
*)
...
...
@@ -1889,6 +1887,9 @@ AM_CONDITIONAL([HAVE_DOXYGEN], [test "$ac_cv_path_ac_pt_DX_DOXYGEN" != ""])
# Check if using EAGLE cooling
AM_CONDITIONAL([HAVEEAGLECOOLING], [test $with_cooling = "EAGLE"])
# Check if using EAGLE feedback
AM_CONDITIONAL([HAVEEAGLEFEEDBACK], [test $with_feedback = "EAGLE"])
# Handle .in files.
AC_CONFIG_FILES([Makefile src/Makefile examples/Makefile examples/Cooling/CoolingRates/Makefile doc/Makefile doc/Doxyfile tests/Makefile])
AC_CONFIG_FILES([argparse/Makefile tools/Makefile])
...
...
src/Makefile.am
View file @
36b7e032
...
...
@@ -59,6 +59,12 @@ if HAVEEAGLECOOLING
EAGLE_COOLING_SOURCES
+=
cooling/EAGLE/cooling.c cooling/EAGLE/cooling_tables.c
endif
# source files for EAGLE feedback
EAGLE_FEEDBACK_SOURCES
=
if
HAVEEAGLEFEEDBACK
EAGLE_FEEDBACK_SOURCES
+=
feedback/EAGLE/feedback.c
endif
# Common source files
AM_SOURCES
=
space.c runner.c queue.c task.c cell.c engine.c engine_maketasks.c
\
engine_marktasks.c engine_drift.c serial_io.c timers.c debug.c scheduler.c
\
...
...
@@ -71,7 +77,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c engine_maketasks.c
collectgroup.c hydro_space.c equation_of_state.c
\
chemistry.c cosmology.c restart.c mesh_gravity.c velociraptor_interface.c
\
outputlist.c velociraptor_dummy.c logger_io.c memuse.c
\
$(EAGLE_COOLING_SOURCES)
$(EAGLE_COOLING_SOURCES)
$(EAGLE_FEEDBACK_SOURCES)
# Include files for distribution, not installation.
nobase_noinst_HEADERS
=
align.h approx_math.h atomic.h barrier.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h
\
...
...
src/feedback.h
0 → 100644
View file @
36b7e032
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2018 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_FEEDBACK_H
#define SWIFT_FEEDBACK_H
/* Config parameters. */
#include
"../config.h"
/* Select the correct feedback model */
#if defined(FEEDBACK_NONE)
#include
"./feedback/Default/feedback.h"
#include
"./deedback/Default/feedback_iact.h"
#elif defined(FEEDBACK_EAGLE)
#include
"./feedback/EAGLE/feedback.h"
#include
"./feedback/EAGLE/feedback_iact.h"
#else
#error "Invalid choice of feedback model"
#endif
#endif
src/feedback/EAGLE/feedback.c
0 → 100644
View file @
36b7e032
This diff is collapsed.
Click to expand it.
src/feedback/EAGLE/feedback.h
View file @
36b7e032
This diff is collapsed.
Click to expand it.
src/feedback/EAGLE/feedback_iact.h
View file @
36b7e032
...
...
@@ -16,7 +16,7 @@
* @param ti_current Current integer time value
*/
__attribute__
((
always_inline
))
INLINE
static
void
runner_iact_nonsym_
stars
_density
(
runner_iact_nonsym_
feedback
_density
(
float
r2
,
const
float
*
dx
,
float
hi
,
float
hj
,
struct
spart
*
restrict
si
,
const
struct
part
*
restrict
pj
,
const
struct
cosmology
*
restrict
cosmo
,
const
struct
stars_props
*
restrict
stars_properties
,
...
...
@@ -42,14 +42,14 @@ runner_iact_nonsym_stars_density(
kernel_deval
(
uj
,
&
wj
,
&
wj_dx
);
/* Add mass of pj to neighbour mass of si */
si
->
ngb_mass
+=
hydro_get_mass
(
pj
);
si
->
feedback_data
.
ngb_mass
+=
hydro_get_mass
(
pj
);
/* Add contribution of pj to normalisation of density weighted fraction
* which determines how much mass to distribute to neighbouring
* gas particles */
const
float
rho
=
hydro_get_comoving_density
(
pj
);
if
(
rho
!=
0
)
si
->
density_weighted_frac_normalisation_inv
+=
wj
/
rho
;
if
(
rho
!=
0
)
si
->
feedback_data
.
density_weighted_frac_normalisation_inv
+=
wj
/
rho
;
/* Compute contribution to the density */
si
->
rho_gas
+=
mj
*
wi
;
...
...
@@ -72,7 +72,7 @@ runner_iact_nonsym_stars_density(
* generator
*/
__attribute__
((
always_inline
))
INLINE
static
void
runner_iact_nonsym_
stars_
feedback
(
runner_iact_nonsym_feedback
_apply
(
float
r2
,
const
float
*
dx
,
float
hi
,
float
hj
,
const
struct
spart
*
restrict
si
,
struct
part
*
restrict
pj
,
const
struct
cosmology
*
restrict
cosmo
,
...
...
src/feedback/EAGLE/feedback_properties.h
View file @
36b7e032
...
...
@@ -46,7 +46,7 @@ struct lifetime_table {
double
**
dyingtime
;
};
struct
feedback_prop
ertie
s
{
struct
feedback_props
{
/* Flag to switch between continuous and stochastic heating */
int
continuous_heating
;
...
...
@@ -150,68 +150,3 @@ struct feedback_properties {
float
spart_first_init_birth_time
;
};
/**
* @brief Initialize the global properties of the feedback scheme.
*
* By default, takes the values provided by the hydro.
*
* @param sp The #feedback_properties.
* @param phys_const The physical constants in the internal unit system.
* @param us The internal unit system.
* @param params The parsed parameters.
* @param p The already read-in properties of the hydro scheme.
*/
INLINE
static
void
feedbakc_props_init
(
struct
feedback_props
*
fp
,
const
struct
phys_const
*
phys_const
,
const
struct
unit_system
*
us
,
struct
swift_params
*
params
,
const
struct
hydro_props
*
p
,
const
struct
cosmology
*
cosmo
)
{
/* Read SNIa timscale */
fp
->
SNIa_timescale_Gyr
=
parser_get_param_float
(
params
,
"EAGLEFeedback:SNIa_timescale_Gyr"
);
/* Read the efficiency of producing SNIa */
fp
->
SNIa_efficiency
=
parser_get_param_float
(
params
,
"EAGLEFeedback:SNIa_efficiency"
);
/* Are we doing continuous heating? */
fp
->
continuous_heating
=
parser_get_param_int
(
params
,
"EAGLEFeedback:continuous_heating_switch"
);
/* Set the delay time before SNII occur */
const
float
Gyr_in_cgs
=
1e9
*
365
*
24
*
3600
;
fp
->
SNII_wind_delay
=
parser_get_param_float
(
params
,
"EAGLEFeedback:SNII_wind_delay_Gyr"
)
*
Gyr_in_cgs
/
units_cgs_conversion_factor
(
us
,
UNIT_CONV_TIME
);
/* Read the temperature change to use in stochastic heating */
fp
->
SNe_deltaT_desired
=
parser_get_param_float
(
params
,
"EAGLEFeedback:SNe_heating_temperature_K"
);
fp
->
SNe_deltaT_desired
/=
units_cgs_conversion_factor
(
us
,
UNIT_CONV_TEMPERATURE
);
/* Set ejecta thermal energy */
const
float
ejecta_velocity
=
1.0e6
/
units_cgs_conversion_factor
(
us
,
UNIT_CONV_SPEED
);
// EAGLE parameter is 10 km/s
fp
->
ejecta_specific_thermal_energy
=
0
.
5
*
ejecta_velocity
*
ejecta_velocity
;
/* Energy released by supernova */
fp
->
total_energy_SNe
=
1.0e51
/
units_cgs_conversion_factor
(
us
,
UNIT_CONV_ENERGY
);
/* Calculate temperature to internal energy conversion factor */
fp
->
temp_to_u_factor
=
phys_const
->
const_boltzmann_k
/
(
p
->
mu_ionised
*
(
hydro_gamma_minus_one
)
*
phys_const
->
const_proton_mass
);
/* Read birth time to set all stars in ICs to (defaults to -1 to indicate star
* present in ICs) */
fp
->
spart_first_init_birth_time
=
parser_get_opt_param_float
(
params
,
"EAGLEFeedback:birth_time_override"
,
-
1
);
/* Copy over solar mass */
fp
->
const_solar_mass
=
phys_const
->
const_solar_mass
;
}
src/feedback/EAGLE/feedback_struct.h
View file @
36b7e032
struct
feedback_struct_data
{
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2018 Matthieu Schaller (schaller@strw.leidenuniv.nl)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_FEEDBACK_STRUCT_EAGLE_H
#define SWIFT_FEEDBACK_STRUCT_EAGLE_H
struct
feedback_part_data
{
struct
{
...
...
@@ -58,4 +77,8 @@ struct feedback_struct_data {
/* total mass (unweighted) of neighbouring gas particles */
float
ngb_mass
;
}
};
#endif
/* SWIFT_FEEDBACK_STRUCT_EAGLE_H */
src/feedback/EAGLE/imf.h
View file @
36b7e032
...
...
@@ -20,6 +20,9 @@
#define SWIFT_EAGLE_STARS_IMF_H
#include
"interpolate.h"
#include
"minmax.h"
#include
<string.h>
/**
* @brief the different weightings allowed for the IMF integration
...
...
@@ -42,15 +45,15 @@ enum eagle_imf_integration_type {
*/
inline
static
void
determine_imf_bins
(
double
log10_min_mass
,
double
log10_max_mass
,
int
*
i_min
,
int
*
i_max
,
const
struct
stars_props
*
star_propertie
s
)
{
const
struct
feedback_props
*
feedback_prop
s
)
{
#ifdef SWIFT_DEBUG_CHECKS
if
(
log10_min_mass
>
log10_max_mass
)
error
(
"Lower bound higher than larger bound."
);
#endif
const
int
N_bins
=
star_properties
->
feedback
.
n_imf_mass_bins
;
const
float
*
imf_bins_log10
=
star_properties
->
feedback
.
imf_mass_bin_log10
;
const
int
N_bins
=
feedback_props
->
n_imf_mass_bins
;
const
float
*
imf_bins_log10
=
feedback_props
->
imf_mass_bin_log10
;
/* Check whether lower mass is within the IMF mass bin range */
log10_min_mass
=
max
(
log10_min_mass
,
imf_bins_log10
[
0
]);
...
...
@@ -87,14 +90,14 @@ inline static float integrate_imf(const float log10_min_mass,
const
float
log10_max_mass
,
const
enum
eagle_imf_integration_type
mode
,
const
float
*
stellar_yields
,
const
struct
stars_props
*
star_propertie
s
)
{
const
struct
feedback_props
*
feedback_prop
s
)
{
/* Pull out some common terms */
const
int
N_bins
=
star_properties
->
feedback
.
n_imf_mass_bins
;
const
float
*
imf
=
star_properties
->
feedback
.
imf
;
const
float
*
imf_mass_bin
=
star_properties
->
feedback
.
imf_mass_bin
;
const
int
N_bins
=
feedback_props
->
n_imf_mass_bins
;
const
float
*
imf
=
feedback_props
->
imf
;
const
float
*
imf_mass_bin
=
feedback_props
->
imf_mass_bin
;
const
float
*
imf_mass_bin_log10
=
star_properties
->
feedback
.
imf_mass_bin_log10
;
feedback_props
->
imf_mass_bin_log10
;
/* IMF mass bin spacing in log10 space. Assumes uniform spacing. */
const
float
imf_log10_mass_bin_size
=
...
...
@@ -103,7 +106,7 @@ inline static float integrate_imf(const float log10_min_mass,
/* Determine bins to integrate over based on integration bounds */
int
i_min
,
i_max
;
determine_imf_bins
(
log10_min_mass
,
log10_max_mass
,
&
i_min
,
&
i_max
,
star_propertie
s
);
feedback_prop
s
);
/* Array for the integrand */
float
integrand
[
N_bins
];
...
...
@@ -189,82 +192,82 @@ inline static float integrate_imf(const float log10_min_mass,
* table.
*
* @param star_properties #stars_props data structure */
inline
static
void
init_imf
(
struct
stars
_props
*
restrict
star_propertie
s
)
{
inline
static
void
init_imf
(
struct
feedback
_props
*
restrict
feedback_prop
s
)
{
/* Define number of imf mass bins */
star_properties
->
feedback
.
n_imf_mass_bins
=
200
;
feedback_props
->
n_imf_mass_bins
=
200
;
/* Define max and min imf masses */
star_properties
->
feedback
.
imf_max_mass_msun
=
100
.
f
;
star_properties
->
feedback
.
imf_min_mass_msun
=
0
.
1
;
star_properties
->
feedback
.
log10_imf_max_mass_msun
=
log10
(
star_properties
->
feedback
.
imf_max_mass_msun
);
star_properties
->
feedback
.
log10_imf_min_mass_msun
=
log10
(
star_properties
->
feedback
.
imf_min_mass_msun
);
feedback_props
->
imf_max_mass_msun
=
100
.
f
;
feedback_props
->
imf_min_mass_msun
=
0
.
1
;
feedback_props
->
log10_imf_max_mass_msun
=
log10
(
feedback_props
->
imf_max_mass_msun
);
feedback_props
->
log10_imf_min_mass_msun
=
log10
(
feedback_props
->
imf_min_mass_msun
);
/* Compute size of mass bins in log10 space */
const
float
imf_log10_mass_bin_size
=
(
star_properties
->
feedback
.
log10_imf_max_mass_msun
-
star_properties
->
feedback
.
log10_imf_min_mass_msun
)
/
(
float
)(
star_properties
->
feedback
.
n_imf_mass_bins
-
1
);
(
feedback_props
->
log10_imf_max_mass_msun
-
feedback_props
->
log10_imf_min_mass_msun
)
/
(
float
)(
feedback_props
->
n_imf_mass_bins
-
1
);
/* Allocate IMF array */
if
(
swift_memalign
(
"imf-tables"
,
(
void
**
)
&
star_properties
->
feedback
.
imf
,
"imf-tables"
,
(
void
**
)
&
feedback_props
->
imf
,
SWIFT_STRUCT_ALIGNMENT
,
star_properties
->
feedback
.
n_imf_mass_bins
*
sizeof
(
float
))
!=
0
)
feedback_props
->
n_imf_mass_bins
*
sizeof
(
float
))
!=
0
)
error
(
"Failed to allocate IMF bins table"
);
/* Allocate array to store IMF mass bins */
if
(
swift_memalign
(
"imf-tables"
,
(
void
**
)
&
star_properties
->
feedback
.
imf_mass_bin
,
"imf-tables"
,
(
void
**
)
&
feedback_props
->
imf_mass_bin
,
SWIFT_STRUCT_ALIGNMENT
,
star_properties
->
feedback
.
n_imf_mass_bins
*
sizeof
(
float
))
!=
0
)
feedback_props
->
n_imf_mass_bins
*
sizeof
(
float
))
!=
0
)
error
(
"Failed to allocate IMF bins table"
);
/* Allocate array to store IMF mass bins in log10 space */
if
(
swift_memalign
(
"imf-tables"
,
(
void
**
)
&
star_properties
->
feedback
.
imf_mass_bin_log10
,
"imf-tables"
,
(
void
**
)
&
feedback_props
->
imf_mass_bin_log10
,
SWIFT_STRUCT_ALIGNMENT
,
star_properties
->
feedback
.
n_imf_mass_bins
*
sizeof
(
float
))
!=
0
)
feedback_props
->
n_imf_mass_bins
*
sizeof
(
float
))
!=
0
)
error
(
"Failed to allocate IMF bins table"
);
/* Set IMF from Chabrier 2003 */
if
(
strcmp
(
star_properties
->
feedback
.
IMF_Model
,
"Chabrier"
)
==
0
)
{
for
(
int
i
=
0
;
i
<
star_properties
->
feedback
.
n_imf_mass_bins
;
i
++
)
{
if
(
strcmp
(
feedback_props
->
IMF_Model
,
"Chabrier"
)
==
0
)
{
for
(
int
i
=
0
;
i
<
feedback_props
->
n_imf_mass_bins
;
i
++
)
{
const
float
log10_mass_msun
=
star_properties
->
feedback
.
log10_imf_min_mass_msun
+
feedback_props
->
log10_imf_min_mass_msun
+
i
*
imf_log10_mass_bin_size
;
const
float
mass_msun
=
exp10f
(
log10_mass_msun
);
if
(
mass_msun
>
1
.
0
)
star_properties
->
feedback
.
imf
[
i
]
=
0
.
237912
*
pow
(
mass_msun
,
-
2
.
3
);
feedback_props
->
imf
[
i
]
=
0
.
237912
*
pow
(
mass_msun
,
-
2
.
3
);
else
star_properties
->
feedback
.
imf
[
i
]
=
feedback_props
->
imf
[
i
]
=
0
.
852464
*
exp
((
log10
(
mass_msun
)
-
log10
(
0
.
07
9
))
*
(
log10
(
mass_msun
)
-
log10
(
0
.
07
9
))
/
(
-
2
.
0
*
pow
(
0
.
69
,
2
)))
/
mass_msun
;
star_properties
->
feedback
.
imf_mass_bin
[
i
]
=
mass_msun
;
star_properties
->
feedback
.
imf_mass_bin_log10
[
i
]
=
log10_mass_msun
;
feedback_props
->
imf_mass_bin
[
i
]
=
mass_msun
;
feedback_props
->
imf_mass_bin_log10
[
i
]
=
log10_mass_msun
;
}
}
else
{
error
(
"Invalid IMF model %s. Valid models are: 'Chabrier'
\n
"
,
star_properties
->
feedback
.
IMF_Model
);
feedback_props
->
IMF_Model
);
}
/* Normalize the IMF */
const
float
norm
=
integrate_imf
(
star_properties
->
feedback
.
log10_imf_min_mass_msun
,
star_properties
->
feedback
.
log10_imf_max_mass_msun
,
integrate_imf
(
feedback_props
->
log10_imf_min_mass_msun
,
feedback_props
->
log10_imf_max_mass_msun
,
eagle_imf_integration_mass_weight
,
/*(stellar_yields=)*/
NULL
,
star_propertie
s
);
/*(stellar_yields=)*/
NULL
,
feedback_prop
s
);
for
(
int
i
=
0
;
i
<
star_properties
->
feedback
.
n_imf_mass_bins
;
i
++
)
star_properties
->
feedback
.
imf
[
i
]
/=
norm
;
for
(
int
i
=
0
;
i
<
feedback_props
->
n_imf_mass_bins
;
i
++
)
feedback_props
->
imf
[
i
]
/=
norm
;
}
/**
...
...
@@ -279,7 +282,7 @@ inline static void init_imf(struct stars_props *restrict star_properties) {
*/
inline
static
double
dying_mass_msun
(
const
float
age_Gyr
,
const
float
metallicity
,
const
struct
stars
_props
*
star_props
)
{
const
struct
feedback
_props
*
star_props
)
{
float
metallicity_offset
,
time_offset_low_metallicity
=
0
,
time_offset_high_metallicity
=
0
,
...
...
@@ -289,104 +292,104 @@ inline static double dying_mass_msun(const float age_Gyr,
time_index_high_metallicity
=
-
1
;
if
(
age_Gyr
<=
0
)
{
return
star_props
->
feedback
.
imf_max_mass_msun
;
return
star_props
->
imf_max_mass_msun
;
}
const
float
log10_age_yr
=
log10f
(
age_Gyr
*
1.e9
);
if
(
metallicity
<=
star_props
->
feedback
.
lifetimes
.
metallicity
[
0
])
{
if
(
metallicity
<=
star_props
->
lifetimes
.
metallicity
[
0
])
{
metallicity_index
=
0
;
metallicity_offset
=
0
.
0
;
}
else
if
(
metallicity
>=
star_props
->
feedback
.
lifetimes
.
metallicity
[
star_props
->
feedback
.
lifetimes
.
n_z
-
1
])
{
metallicity_index
=
star_props
->
feedback
.
lifetimes
.
n_z
-
2
;
star_props
->
lifetimes
.
metallicity
[
star_props
->
lifetimes
.
n_z
-
1
])
{
metallicity_index
=
star_props
->
lifetimes
.
n_z
-
2
;
metallicity_offset
=
1
.
0
;
}
else
{
for
(
metallicity_index
=
0
;
metallicity_index
<
star_props
->
feedback
.
lifetimes
.
n_z
-
1
;
metallicity_index
<
star_props
->
lifetimes
.
n_z
-
1
;
metallicity_index
++
)
if
(
star_props
->
feedback
.
lifetimes
.
metallicity
[
metallicity_index
+
1
]
>
if
(
star_props
->
lifetimes
.
metallicity
[
metallicity_index
+
1
]
>
metallicity
)
break
;
metallicity_offset
=
(
metallicity
-
star_props
->
feedback
.
lifetimes
.
metallicity
[
metallicity_index
])
/
(
star_props
->
feedback
.
lifetimes
.
metallicity
[
metallicity_index
+
1
]
-
star_props
->
feedback
.
lifetimes
.
metallicity
[
metallicity_index
]);
star_props
->
lifetimes
.
metallicity
[
metallicity_index
])
/
(
star_props
->
lifetimes
.
metallicity
[
metallicity_index
+
1
]
-
star_props
->
lifetimes
.
metallicity
[
metallicity_index
]);
}
if
(
log10_age_yr
>=
star_props
->
feedback
.
lifetimes
.
dyingtime
[
metallicity_index
][
0
])
{
star_props
->
lifetimes
.
dyingtime
[
metallicity_index
][
0
])
{
time_index_low_metallicity
=
0
;
time_offset_low_metallicity
=
0
.
0
;
}
else
if
(
log10_age_yr
<=
star_props
->
feedback
.
lifetimes
star_props
->
lifetimes
.
dyingtime
[
metallicity_index
]
[
star_props
->
feedback
.
lifetimes
.
n_mass
-
1
])
{
time_index_low_metallicity
=
star_props
->
feedback
.
lifetimes
.
n_mass
-
2
;
[
star_props
->
lifetimes
.
n_mass
-
1
])
{
time_index_low_metallicity
=
star_props
->
lifetimes
.
n_mass
-
2
;
time_offset_low_metallicity
=
1
.
0
;
}
if
(
log10_age_yr
>=
star_props
->
feedback
.
lifetimes
.
dyingtime
[
metallicity_index
+
1
][
0
])
{
star_props
->
lifetimes
.
dyingtime
[
metallicity_index
+
1
][
0
])
{
time_index_high_metallicity
=
0
;
time_offset_high_metallicity
=
0
.
0
;
}
else
if
(
log10_age_yr
<=
star_props
->
feedback
.
lifetimes
star_props
->
lifetimes
.
dyingtime
[
metallicity_index
+
1
]
[
star_props
->
feedback
.
lifetimes
.
n_mass
-
1
])
{
time_index_high_metallicity
=
star_props
->
feedback
.
lifetimes
.
n_mass
-
2
;
[
star_props
->
lifetimes
.
n_mass
-
1
])
{
time_index_high_metallicity
=
star_props
->
lifetimes
.
n_mass
-
2
;
time_offset_high_metallicity
=
1
.
0
;
}
int
i
=
star_props
->
feedback
.
lifetimes
.
n_mass
;
int
i
=
star_props
->
lifetimes
.
n_mass
;
while
(
i
>=
0
&&
(
time_index_low_metallicity
==
-
1
||
time_index_high_metallicity
==
-
1
))
{
i
--
;
if
(
star_props
->
feedback
.
lifetimes
.
dyingtime
[
metallicity_index
][
i
]
>=
if
(
star_props
->
lifetimes
.
dyingtime
[
metallicity_index
][
i
]
>=
log10_age_yr
&&
time_index_low_metallicity
==
-
1
)
{
time_index_low_metallicity
=
i
;
time_offset_low_metallicity
=
(
log10_age_yr
-
star_props
->
feedback
.
lifetimes
star_props
->
lifetimes
.
dyingtime
[
metallicity_index
][
time_index_low_metallicity
])
/
(
star_props
->
feedback
.
lifetimes
(
star_props
->
lifetimes
.
dyingtime
[
metallicity_index
][
time_index_low_metallicity
+
1
]
-
star_props
->
feedback
.
lifetimes
star_props
->
lifetimes
.
dyingtime
[
metallicity_index
][
time_index_low_metallicity
]);
}
if
(
star_props
->
feedback
.
lifetimes
.
dyingtime
[
metallicity_index
+
1
][
i
]
>=
if
(
star_props
->
lifetimes
.
dyingtime
[
metallicity_index
+
1
][
i
]
>=
log10_age_yr
&&
time_index_high_metallicity
==
-
1
)
{
time_index_high_metallicity
=
i
;
time_offset_high_metallicity
=
(
log10_age_yr
-
star_props
->
feedback
.
lifetimes
star_props
->
lifetimes
.
dyingtime
[
metallicity_index
+
1
][
time_index_high_metallicity
])
/
(
star_props
->
feedback
.
lifetimes
(
star_props
->
lifetimes
.
dyingtime
[
metallicity_index
+
1
]
[
time_index_high_metallicity
+
1
]
-
star_props
->
feedback
.
lifetimes
star_props
->
lifetimes
.
dyingtime
[
metallicity_index
+
1
][
time_index_high_metallicity
]);
}
}
mass_low_metallicity
=
interpolate_1d
(
star_props
->
feedback
.
lifetimes
.
mass
,
interpolate_1d
(
star_props
->
lifetimes
.
mass
,
time_index_low_metallicity
,
time_offset_low_metallicity
);
mass_high_metallicity
=
interpolate_1d
(
star_props
->
feedback
.
lifetimes
.
mass
,
interpolate_1d
(
star_props
->
lifetimes
.
mass
,
time_index_high_metallicity
,
time_offset_high_metallicity
);
float
mass
=
(
1
.
0
-
metallicity_offset
)
*
mass_low_metallicity
+
metallicity_offset
*
mass_high_metallicity
;
/* Check that we haven't killed too many stars */
if
(
mass
>
star_props
->
feedback
.
imf_max_mass_msun
)
mass
=
star_props
->
feedback
.
imf_max_mass_msun
;
if
(
mass
>
star_props
->
imf_max_mass_msun
)
mass
=
star_props
->
imf_max_mass_msun
;
return
mass
;
}
...
...
@@ -400,61 +403,61 @@ inline static double dying_mass_msun(const float age_Gyr,
* @param star_properties the #stars_props data struct
*/
inline
static
float
lifetime_in_Gyr
(
float
mass
,
float
metallicity
,
const
struct
stars
_props
*
restrict
star_propertie
s
)
{
const
struct
feedback
_props
*
restrict
feedback_prop
s
)
{
double
time_Gyr
=
0
,
mass_offset
,
metallicity_offset
;
int
mass_index
,
metallicity_index
;
if
(
mass
<=
star_properties
->
feedback
.
lifetimes
.
mass
[
0
])
{
if
(
mass
<=
feedback_props
->
lifetimes
.
mass
[
0
])
{
mass_index
=
0
;
mass_offset
=
0
.
0
;
}
else
if
(
mass
>=
star_properties
->
feedback
.
lifetimes