/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2019 Joel Pfeffer (j.l.pfeffer@ljmu.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 .
*
******************************************************************************/
/* Standard includes */
#include
/* Standard includes */
#include
#include
/* Local includes */
#include "dsyevj3.h"
#include "engine.h"
#include "stars_io.h"
/**
* @brief Free-fall time as a function of density
*
* @param rho The gas density in kg/m^3
* @param G Gravitational constant
* @return Free-fall time
*/
__attribute__((const)) INLINE static double f_tff(const double rho,
const double G) {
/* free-fall time sqrt(3 * pi / (32 * G * rho) */
return sqrt(0.09375 * M_PI / (G * rho));
}
/**
* @brief Compute the Mach number
*
* @param sigmaloc Velocity dispersion
* @param csloc Sound speed
* @return Mach number
*/
__attribute__((const)) INLINE static double f_mach(const double sigmaloc,
const double csloc) {
return sigmaloc / csloc;
}
/**
* @brief Dispersion of overdensity PDF as a function of Mach number and
* magnetic pressure ratio, based on e.g. Padoan & Nordlund (2011).
*
* @param mach Mach number
* @param beta0 Magnetic pressure ratio
* @return Dispersion of overdensity PDF
*/
__attribute__((const)) INLINE static double f_sigrho(const double mach,
const double beta0) {
/* Constant, e.g. Padoan & Nordlund 2002 */
const double b = 0.5;
return sqrt(log(1. + 3. * b * b * mach * mach * beta0 / (beta0 + 1.)));
}
/**
* @brief Critical overdensity for star formation in the sSFR_ff.
*
* Equation 27 of Krumholz & McKee, 2015, ApJ, 630, 1, pp. 250-268.
*
* @param qvir GMC virial ratio
* @param mach Mach number
* @return Crititcal overdensity
*/
__attribute__((const)) INLINE static double f_xcrit(const double qvir,
const double mach) {
/* constant (table 2) */
const double phix = 1.12;
return (1. / 15.) * M_PI * M_PI * phix * phix * qvir * mach * mach;
}
/**
* @brief Specific star formation rate per free-fall time.
*
* Uses Eq. 20 of Krumholz & McKee, 2015, ApJ, 630, 1, pp. 250-268
* when sflaw == 1
* OR
* Eq. 13 of Elmegreen, 2002, ApJ, 577, 1, pp. 206-220
* which is just a constant value.
*
* @param mach Mach number
* @param props Properties of the stars model
* @return Specific star formation rate per free-fall time
*/
__attribute__((const)) INLINE static double f_sfrff(
const double mach, const struct stars_props* props) {
/* Star formation law */
if (props->sflaw == 1) {
/* Krumholz & McKee (2005). Depends on parameters:
* qvir - GMC virial ratio
* beta0 - Magnetic pressure ratio
* ecore - Formation efficiency of protostellar cores
*/
/* constant (table 2) */
const double phit = 1.91;
/* critical overdensity for star formation */
const double xcrit = f_xcrit(props->qvir, mach);
/* overdensity PDF dispersion */
const double sigrho = f_sigrho(mach, props->beta0);
const double erf_term =
(-2. * log(xcrit) + sigrho * sigrho) / (2. * M_SQRT2 * sigrho);
return 0.5 * props->ecore * (1. + erf(erf_term)) / phit;
} else {
/* Elmegreen (2002) */
return props->const_sfrff;
}
}
/**
* @brief Overdensity PDF of the ISM
*
* @param x Overdensity
* @param mulnx Logarithmic mean
* @param sig Dispersion
* @return Overdensity PDF
*/
__attribute__((const)) INLINE static double f_dpdx(const double x,
const double mulnx,
const double sig) {
const double dif = log(x) - mulnx;
/* 1. / (sqrt(2 pi) * sig * x) * exp (-dif^2 / 2 sig^2)
*
* Note: 1 / sqrt(2 pi) = (1/4) * (2 / sqrt(pi)) * sqrt(2) */
return ((0.25 * M_2_SQRTPI * M_SQRT2) / (sig * x)) *
exp(-0.5 * dif * dif / (sig * sig));
}
/**
* @brief Estimate for gas surface density, assuming an equilibrium disc
*
* @param rholoc Gas volume density
* @param sigmaloc Gas velocity dispersion
* @param fgas Local gas fraction
* @param G Gravitational constant
* @return Gas surface density
*/
__attribute__((const)) INLINE static double f_surfg(const double rholoc,
const double sigmaloc,
const double fgas,
const double G) {
return sqrt(2. * fgas * rholoc * sigmaloc * sigmaloc / (M_PI * G));
}
/**
* @brief Calculate local bound fraction of star formation
*
* @param epsilon Gas volume density
* @param logZ log10 of stellar metal mass fraction Z
* @param props Properties of the stars model
* @return local bound fraction
*/
__attribute__((const)) INLINE static double f_fbound(
const double epsilon, const double logZ, const struct stars_props* props) {
double bound;
if (props->CFE_met_dep_fbound) {
/* Metallicity dependence for fbound-SFE, based on Grudic et al. 2021.
* Approximate fit to simulations in their figure 3 with a double power
* law (with indices alpha and 0) and shape parameter epsilon^beta.
* "Saturation" value from Li et al. 2019 to shift vertical scale.
* Linearly interpolate the power-law index in log(Z) so that
* fbound ~ Z^beta (with beta a function of epsilon), consistent with a
* power-law scaling for wind mass loss rates (e.g. Vink et al. 2001). The
* lower metallicity value sets the lower alpha limit based on the
* "no winds" run in Grudic et al. 2021. */
double alpha;
if (logZ < props->fbound_logZ[0]) {
/* Lower metallicity value sets lower limit on alpha */
alpha = props->fbound_alpha[0];
} else {
/* Interpolate the power law index by particle metallicity */
alpha = props->fbound_alpha[0] +
(logZ - props->fbound_logZ[0]) *
(props->fbound_alpha[1] - props->fbound_alpha[0]) /
(props->fbound_logZ[1] - props->fbound_logZ[0]);
}
bound =
props->fbound_saturation *
pow((1. + props->fbound_scale) /
(1. + props->fbound_scale / pow(epsilon, props->fbound_beta)),
alpha);
} else {
/* Original one-to-one assumption in K12 */
bound = epsilon;
}
return bound;
}
/**
* @brief Naturally bound fraction of star formation
*
* @param rholoc Gas volume density
* @param sigmaloc Gas velocity dispersion
* @param mach Mach number
* @param xloc Overdensity
* @param surfGMC Giant molecular cloud surface density
* @param sfrff Specific star formation rate per free-fall time
* @param props Properties of the stars model
* @param phys_const Physical constants
* @return Local SFE
*/
static double f_fstar(const double rholoc, const double sigmaloc,
const double mach, const double xloc,
const double surfGMC, const double sfrff,
const struct stars_props* props,
const struct phys_const* phys_const) {
/* Extract some physics constants (in internal units) */
const double k_B = phys_const->const_stefan_boltzmann;
const double c = phys_const->const_speed_light_c;
const double G = phys_const->const_newton_G;
/* Extract some model constants */
const double t_SN = props->tsn;
const double t_view = props->tview;
const double phi_fb = props->const_phifb;
const double kappa0 = props->const_kappa0;
const double psi = props->const_psi;
const double ecore = props->ecore;
/* trapping ratio */
const double phitrap = 0.2;
/* Free-fall time */
const double t_ff = f_tff(xloc * rholoc, phys_const->const_newton_G);
/* Star Formation Efficiencies */
/* SN feedback */
double efb;
if (props->radfb == 0 || props->radfb == 2) {
efb = 0.5 * sfrff * (t_SN / t_ff) *
(1. + sqrt(1. + 4. * t_ff * sigmaloc * sigmaloc /
(phi_fb * sfrff * t_SN * t_SN * xloc)));
} else
efb = 1.;
/* radiative feedback */
double efbrad;
if (props->radfb > 0.) {
/* surface density on which radiative feedback acts */
const double surffb3 = surfGMC * surfGMC * surfGMC;
const double factor = phitrap * kappa0 * kappa0 / k_B;
efbrad = 2. / (factor * psi * surffb3) *
(sqrt(1. + 2. * M_PI * c * G * factor * surfGMC * surffb3) - 1.);
} else
efbrad = 1.;
/* star formation is incomplete/still ongoing */
const double einc = sfrff * t_view / t_ff;
/* Local SFE is minimum of all efficiencies:
* min[maximum, SNfeedback, radiative feedback, incomplete] */
return min4(ecore, efb, efbrad, einc);
}
/**
* @brief Calculate the fraction of star formation occurring in bound clusters
* i.e. the "cluster formation efficiency"
* Kruijssen, J. M. D., 2012, MNRAS 426, 3008
*
* @param rholoc Local gas volume density in kg/m^3
* @param sigmaloc Local 1D gas velocity dispersion in m/s
* @param csloc Local gas sound speed in m/s
* @param fgas Local gas fraction
* @param logZ log10 of stellar metal mass fraction Z
* @param props Properties of the stars model
* @param phys_const Physical constants
* @return The CFE
*/
static double f_cfelocal(const double rholoc, const double sigmaloc,
const double csloc, const double fgas,
const double logZ, const struct stars_props* props,
const struct phys_const* phys_const) {
/* Parameters:
* sflaw - star formation law: 0=Elmegreen (2002, default)
* 1=Krumholz & McKee (2005)
* qvir - giant molecular cloud virial parameter (default=1.3)
* tsn - time of the first supernova (default=3.*Myr)
* tview - time at which CFE is determined (default=10.*Myr)
* surfGMC - giant molecular cloud surface density (default=100.*Msun/pc^2.)
* ecore - maximum (protostellar core) star formation efficiency (default=.5)
* beta0 - turbulent-to-magnetic pressure ratio (default=1.e10)
* radfb - feedback: 0=supernova, 1=radiative, 2=both (default=0)
* cce - Include cruel cradle effect? (default=0)
* xsurv - Critical overdensity for cruel cradle effect (default=0.)
*/
/* Extract some physics constants (in internal units) */
const double G = phys_const->const_newton_G;
/* Integrate overdensity PDFs to obtain bound fraction of star formation.
* Number of integration steps (checked to be sufficient for convergence) */
const int nx = 1000;
/* estimate of gas surface density */
const double surfg = f_surfg(rholoc, sigmaloc, fgas, G);
const double surfGMC =
(surfg > props->surfGMC_min) ? surfg : props->surfGMC_min;
/* mach number */
const double mach = f_mach(sigmaloc, csloc);
/* specific star formation rate per free-fall time */
const double sfrff = f_sfrff(mach, props);
/* dispersion of overdensity PDF */
const double sigrho = f_sigrho(mach, props->beta0);
/* logarithmic mean of overdensity PDF */
const double mulnx = -0.5 * sigrho * sigrho;
/* Minimum overdensity. If calculating the cruel cradle effect and critical
* overdensity below minimum, then adjust */
double xmin1 = exp(mulnx - 5. * sigrho);
if (props->cruel_cradle > 0 && props->xsurv < xmin1) xmin1 = props->xsurv;
/* Maximum overdensity. If calculating the cruel cradle effect and critical
* overdensity below maximum, then adjust */
double xmax = exp(mulnx + 10. * sigrho);
if (props->cruel_cradle > 0 && props->xsurv > xmax) xmax = props->xsurv;
/* initialize first integral */
const double ratio1 = xmax / xmin1;
const double step_size_term1 =
pow(ratio1, 1. / (2. * nx)) - pow(ratio1, -1. / (2. * nx));
double f1 = 0.;
/* denominator integral */
for (int ix = 0; ix < nx; ix++) {
/* Overdensity */
const double xloc = xmin1 * pow(ratio1, (ix - 0.5) / nx);
/* Step size */
const double dx = xloc * step_size_term1;
/* Local SFE */
const double fstar = f_fstar(rholoc, sigmaloc, mach, xloc, surfGMC, sfrff,
props, phys_const);
/* Local bound fraction */
double bound;
if (props->cruel_cradle == 0 || props->cruel_cradle == 2) {
/* If not calculating the cruel cradle effect but the naturally bound
* fraction of SF, the denominator should contain all SF.
* If calculating the cruel cradle effect with respect to all SF, the
* denominator should contain all SF. */
bound = 1.;
} else
bound = f_fbound(fstar / props->ecore, logZ, props);
/* integral part 1 */
const double integral = bound * fstar * xloc;
/* overdensity PDF, i.e. integral part 2 */
const double dpdx = f_dpdx(xloc, mulnx, sigrho);
/* integral */
f1 += integral * dpdx * dx;
}
/* if calculating the cruel cradle effect set minimum overdensity to critical
* overdensity */
const double xmin2 = (props->cruel_cradle > 0) ? props->xsurv : xmin1;
/* initialize second integral */
const double ratio2 = xmax / xmin2;
const double step_size_term2 =
pow(ratio2, 1. / (2. * nx)) - pow(ratio2, -1. / (2. * nx));
double f2 = 0.;
/* numerator integral */
for (int ix = 0; ix < nx; ix++) {
/* overdensity */
const double xloc = xmin2 * pow(ratio2, (ix - 0.5) / nx);
/* Step size */
const double dx = xloc * step_size_term2;
/* local SFE */
const double fstar = f_fstar(rholoc, sigmaloc, mach, xloc, surfGMC, sfrff,
props, phys_const);
/* local bound fraction */
double bound;
if (props->cruel_cradle == 2) {
/* if calculating the cruel cradle effect with respect to all SF, the
* numerator should contain all SF */
bound = 1.;
} else
bound = f_fbound(fstar / props->ecore, logZ, props);
/* integral part 1 */
const double integral = bound * fstar * xloc;
/* overdensity PDF, i.e. integral part 2 */
const double dpdx = f_dpdx(xloc, mulnx, sigrho);
/* integral */
f2 += integral * dpdx * dx;
}
/* The cluster formation efficiency */
if (f1 == 0)
return 0.;
else
return f2 / f1;
}
/**
* @brief Local supernovae feedback timescale
*
* @param rholoc Local gas volume density
* @param sigmaloc Local 1D gas velocity dispersion
* @param csloc Local gas sound speed
* @param props Properties of the stars model
* @param phys_const Physical constants
* @return The supernovae feedback timescale of molecular cloud
*/
static double feedback_timescale(const double rholoc, const double sigmaloc,
const double csloc,
const struct stars_props* props,
const struct phys_const* phys_const) {
/* Overdensity */
/* const double xloc = 1.; */
/* Mach number */
const double mach = f_mach(sigmaloc, csloc);
/* specific star formation rate per free-fall time */
const double sfrff = f_sfrff(mach, props);
/* free-fall time */
const double tff = f_tff(rholoc, phys_const->const_newton_G);
const double tsn = props->tsn;
const double phi_fb = props->const_phifb;
/* SN feedback timescale */
return 0.5 * tsn *
(1. + sqrt(1. + 4. * tff * sigmaloc * sigmaloc /
(phi_fb * sfrff * tsn * tsn)));
}
/**
* @brief Swap two numbers.
*/
__attribute__((always_inline)) INLINE static void swap(double* a, double* b) {
double tmp = *a;
*a = *b;
*b = tmp;
}
/**
* @brief Sort array with length 3
*
* @param val Length 3 array to sort
*/
__attribute__((always_inline)) INLINE static void sort3(double val[3]) {
if (val[0] > val[1]) swap(&val[0], &val[1]);
if (val[0] > val[2]) swap(&val[0], &val[2]);
if (val[1] > val[2]) swap(&val[1], &val[2]);
}
/**
* @brief Local cluster evaporation timescale
*
* @param sp The #spart to act upon.
* @param props Properties of the stars model.
* @param a3_inv Inverse cube of the current expansion factor
* @return Cluster evaporation timescale
*/
double get_local_evaporation_t0(struct spart* restrict sp,
const struct stars_props* props,
const double a3_inv) {
/* Convert tensor to physical units. Note the lower elements are never
* referenced by dsyevj3 */
double tensor[3][3];
tensor[0][0] = (double)sp->tidal_tensor[2][0] * a3_inv;
tensor[1][1] = (double)sp->tidal_tensor[2][1] * a3_inv;
tensor[2][2] = (double)sp->tidal_tensor[2][2] * a3_inv;
tensor[0][1] = (double)sp->tidal_tensor[2][3] * a3_inv;
tensor[0][2] = (double)sp->tidal_tensor[2][4] * a3_inv;
tensor[1][2] = (double)sp->tidal_tensor[2][5] * a3_inv;
/* get the eigenvalues, then sort */
double tide_val[3];
dsyevj3(tensor, tide_val);
sort3(tide_val);
double Omega2;
if (props->Omega_is_lambda2 && tide_val[2] < 0) {
/* Correct version of Omega for a spherical distribution, in principle.
* When tideval[2] > 0 we can't guarantee that tideval[1] will be < 0 */
Omega2 = fabs(-tide_val[1]);
} else {
/* Version in E-MOSAICS */
Omega2 = fabs(-tide_val[0] - tide_val[1] - tide_val[2]) * (1. / 3.);
}
/* Strength of the tidal field, T + Omega^2 */
const double tidal_strength = tide_val[2] + Omega2;
/* Compute imescale of evaporation in the tidal field */
if (tidal_strength > 0.) {
/* Interpolate between King W0=5 and W0=7 profiles */
const double t0evapsun =
props->const_t0evapsunW5 +
0.5 * (props->W0 - 5.) *
(props->const_t0evapsunW7 - props->const_t0evapsunW5);
return t0evapsun / sqrt(tidal_strength / props->const_tidesun);
} else {
/* Compressive field */
return DBL_MAX;
}
}
/**
* @brief Calculate the tidal heating parameters
*
* @param sp The particle to act upon.
* @param age Age of star.
* @param dt Particle timestep.
* @param time Current time.
* @param a3_inv Inverse cube of the current expansion factor
* @param Mlheatsum (return) Tidal heating for shocks at this timestep.
*/
void calc_heating(struct spart* restrict sp, const double age, const double dt,
const double time, const double a3_inv, double* Mlheatsum) {
/* fraction of maximum at which duration is measured */
const double sigma1_width = 0.88;
const double sigma2_width = 0.61;
/* Time at the previous step */
const float tlast = (float)(time - dt);
/* loop over tidal tensor components, to determine minima, maxima, tshock
* and tidal heating for each component.
* Note we only store upper six components of tidal tensor.
*/
for (int i = 0; i < 6; i++) {
const double tide1 = (double)sp->tidal_tensor[0][i] * a3_inv;
const double tide2 = (double)sp->tidal_tensor[1][i] * a3_inv;
const double tide3 = (double)sp->tidal_tensor[2][i] * a3_inv;
int apply_shock = 0;
/* Is the first point for this particle is a maximum? */
if (age == dt && fabs(tide3) < fabs(tide2)) {
/* Time of last maximum is previous time step */
sp->tmaxsh[i] = tlast;
sp->tidmax[i] = tide2;
sp->extreme_max[i] = 1;
sp->shock_duration[i] = 0.;
sp->shock_indicator[i] = 0;
}
/* If previous time step marks a valid minimum */
if (fabs(tide2) < fabs(tide1) && fabs(tide2) < fabs(tide3) &&
fabs(tide2) < 0.88 * fabs(sp->tidmax[i])) {
/* if previous extreme was a maximum, or this minimum is lower than
* previous */
if (sp->extreme_max[i] || (fabs(tide2) < fabs(sp->tidmin[i]))) {
/* no shock duration is present (shallow max) */
if (sp->shock_indicator[i] == 0)
sp->shock_duration[i] = 0.5 * (tlast - sp->tminsh[i]);
/* Time of last minimum was previous time step */
apply_shock = 1;
sp->tminsh[i] = tlast;
sp->tidmin[i] = tide2;
sp->extreme_max[i] = 0;
}
}
/* If previous time step marks a valid maximum */
if (fabs(tide2) > fabs(tide1) && fabs(tide2) > fabs(tide3) &&
0.88 * fabs(tide2) > fabs(sp->tidmin[i])) {
// if previous extreme was a minimum, or this maximum is higher than
// previous
if (!sp->extreme_max[i] || (fabs(tide2) > fabs(sp->tidmax[i]))) {
/* Time of last maximum was previous time step */
sp->tmaxsh[i] = tlast;
sp->tidmax[i] = tide2;
sp->extreme_max[i] = 1;
sp->shock_duration[i] = 0.f;
sp->shock_indicator[i] = 0;
}
}
/* First order shock duration */
if (sp->extreme_max[i] && sp->shock_indicator[i] == 0 &&
fabs(tide3) < sigma1_width * fabs(sp->tidmax[i])) {
sp->shock_indicator[i] = 1;
/* time interpolation fraction */
const double dt_fraction =
(sigma1_width * fabs(sp->tidmax[i]) - fabs(tide2)) /
(fabs(tide3) - fabs(tide2));
/* interpolate to time of shock duration measurement */
const double shock_time = time + (dt_fraction - 1.) * dt;
/* set shock duration indicator to total width of 1 sigma */
/* 2 * dt / sigma(=1) */
sp->shock_duration[i] = 2. * (shock_time - sp->tmaxsh[i]);
}
/* Second order shock duration */
if (sp->extreme_max[i] && sp->shock_indicator[i] < 2 &&
fabs(tide3) < sigma2_width * fabs(sp->tidmax[i])) {
sp->shock_indicator[i] = 2;
/* time interpolation fraction */
const double dt_fraction =
(sigma2_width * fabs(sp->tidmax[i]) - fabs(tide2)) /
(fabs(tide3) - fabs(tide2));
/* interpolate to time of shock duration measurement */
const double shock_time = time + (dt_fraction - 1.) * dt;
/* set shock duration indicator to total width of 2 sigma */
/* 2 * dt / sigma(=2) */
sp->shock_duration[i] = shock_time - sp->tmaxsh[i];
}
/* If a shock was completed during the previous time step and mass loss
* should be applied */
if (apply_shock) {
/* Tidal heating for the soon-to-occur application of mass loss
* Note this could be negative, but is later applied as heatsum^2 */
Mlheatsum[i] = sp->heatsum[i];
/* reset and start integrating heating for next shock */
sp->shock_indicator[i] = 0;
sp->heatsum[i] = 0.5 * (tide2 + tide3) * dt;
} else {
/* Keep summing tidal heating */
sp->heatsum[i] += 0.5 * (tide2 + tide3) * dt;
}
} /* For each tensor component */
}
/**
* @brief Compute the Weinberg adiabatic correction
*
* @param mass Cluster mass in Solar masses.
* @param rh Cluster half-mass radius in pc.
* @param tshock Tidal shock duration.
* @param W0 King profile parameter.
* @param etapot Potential energy constant.
* @param G Gravitational constant.
* @return Adiabatic correction.
*/
double adiabatic_correction(const double mass, const double rh,
const double tshock, const double W0,
const double etapot, const double G) {
/* constant in omega for King profile W0=[5,7] */
const double cst5 = 0.68;
const double cst7 = 0.82;
/* Interpolate */
const double cst = cst5 + 0.5 * (W0 - 5.) * (cst7 - cst5);
/* angular frequency of stars */
const double omega =
cst * sqrt(8. * etapot * etapot * etapot * G * mass / (rh * rh * rh));
/* number of stellar revolutions per passage */
const double x = omega * tshock;
const double Aw = 1. / ((1. + x * x) * sqrt(1. + x * x));
#ifdef SWIFT_DEBUG_CHECKS
if (!isfinite(Aw)) error("Adiabatic correction is not finite!");
#endif
return Aw;
}
/**
* @brief Compute star cluster mass loss by tidal shocks.
*
* @param sp The #spart.
* @param props Properties of the stars model.
* @param Mlheatsum Tidal heating to apply at this timestep.
* @param mass Star cluster mass.
* @param G Gravitational constant.
* @return Mass loss by tidal shocks.
*/
double mass_loss_shocks(struct spart* restrict sp,
const struct stars_props* props,
const double* Mlheatsum, const double mass,
const double G) {
/* Anything to do? */
double max_heatsum = 0.;
for (int i = 0; i < 6; i++) {
if (fabs(Mlheatsum[i]) > max_heatsum) max_heatsum = fabs(Mlheatsum[i]);
}
if (max_heatsum == 0.) return 0.;
double max_shock_duration = 0.;
for (int i = 0; i < 6; i++) {
if (fabs(sp->shock_duration[i]) > max_shock_duration)
max_shock_duration = fabs(sp->shock_duration[i]);
}
if (max_shock_duration == 0.) return 0.;
/* Ok, now complete the shock */
/* potential energy constant */
const double etapot = 0.4;
/* fraction of shock energy input used in mass loss */
const double efrac = 0.25;
/* mean(M(r)*r)/(Mtot*rh) for King profile W0=[5,7] */
const double zeta5 = 0.81034;
const double zeta7 = 1.0273;
/* ratio of square radii for King profile W0=[5,7] */
const double r2avrh25 = 2.;
const double r2avrh27 = 3.5;
/* Now interpolate */
const double zeta = zeta5 + 0.5 * (props->W0 - 5.) * (zeta7 - zeta5);
const double r2avrh2 =
r2avrh25 + 0.5 * (props->W0 - 5.) * (r2avrh27 - r2avrh25);
/* Half-mass radius */
const double rh = props->rh;
/* Sum up the tidal heating */
double heatparam = 0.;
for (int i = 0; i < 6; i++) {
/* Weinberg adiabatic correction */
double Aw = adiabatic_correction(mass, rh, (double)sp->shock_duration[i],
props->W0, etapot, G);
/* Add the mirror for off-diagonal tensor terms */
double factor = 1.;
if (i > 2) factor = 2.;
heatparam += factor * Mlheatsum[i] * Mlheatsum[i] * Aw;
}
/* first order tsh/P */
const double t1 =
3. * etapot * G * mass / (efrac * rh * rh * rh * r2avrh2 * heatparam);
/* factor to include second order energy gain */
const double InclE2 = 1. + (12. * zeta) / (5. * etapot * r2avrh2);
/* Disruption time / shock period: tsh_int = t1 / InclE2 */
return -mass * InclE2 / t1;
}
/**
* @brief Compute star cluster mass loss by evaporation.
*
* @param props Properties of the stars model.
* @param mass Star cluster mass.
* @param initial_mass Initial star cluster mass.
* @param t0 Evaporation timescale.
* @param dt Particle timestep.
* @param logZ log10 of stellar metal mass fraction Z.
* @param G Gravitational constant.
* @return Mass loss by evaporation.
*/
double mass_loss_evaporation(const struct stars_props* props, const double mass,
const double initial_mass, const double t0,
const double dt, const double logZ, const double G,
const double const_solar_mass) {
/* Dissolution law exponent for W0 = [5,7] King profile */
const double gamma5 = props->gamma_W5;
const double gamma7 = props->gamma_W7;
/* adapt dissolution law exponent to King profile parameter */
const double gamma_min = gamma5 + 0.5 * (props->W0 - 5.) * (gamma7 - gamma5);
double gamma = gamma_min;
if (props->met_dep_evaporation) {
/* Add dependence on metallicity, following Gieles & Gnedin 2023 */
const double gamma_max = props->gamma_max_Zdep;
const double logZ_ymax = props->logZ_gamma_max;
const double logZ_ymin = props->logZ_gamma_min;
gamma = gamma_min + (logZ - logZ_ymin) * (gamma_max - gamma_min) /
(logZ_ymax - logZ_ymin);
if (gamma < gamma_min)
gamma = gamma_min;
else if (gamma > gamma_max)
gamma = gamma_max;
}
double dm = 0.;
/* Add the 'isolated' evaporation term, following Gieles & Baumgardt 08 */
if (props->spitzer_evap_term) {
/* Mass loss fraction per relaxation time in 'isolated' regime
* (Spitzer 87, Gieles & Baumgardt 08) */
const double xi0 = 0.0074;
/* Mean star mass for Kroupa IMF. Convert to internal units */
const double mmean = 0.64 * const_solar_mass;
/* Cluster half-mass radius */
const double rh = props->rh;
/* Relaxation time (Spitzer & Hart 71; Gierz & Heggie 94) */
const double N = mass / mmean;
const double trh = 0.138 * rh * sqrt(rh * N / (mmean * G)) / log(0.11 * N);
dm += -xi0 * mass * dt / trh;
}
/* Fraction of evaporation timescale. Includes internal unit correction */
double dt_t0;
if (gamma != gamma_min) {
/* Gieles & Gnedin 2023 rewrite
* Mdot \propto -Mconst * (M/Mconst)^(1-y)
* as
* Mdot \propto -Mconst * (M/Minit)^(1-y) * (Minit/Mconst)^(1-x)
* where x = 2/3
*/
dt_t0 = dt / t0 * pow(const_solar_mass, gamma_min) *
pow(initial_mass, gamma - gamma_min);
} else {
dt_t0 = dt / t0 * pow(const_solar_mass, gamma);
}
/* First and second order terms */
double dm_order1 = -pow(mass, 1. - gamma) * dt_t0;
double dm_order2 = (1. - gamma) * pow(mass, 1. - 2. * gamma) * dt_t0 * dt_t0;
/* Was it mass loss? Set to first order if mass gain.
* dm>0 only occurs just before disruption because we don't consider third
* order (etc) effects
*/
if (dm_order1 + dm_order2 > 0.)
dm += dm_order1;
else
dm += dm_order1 + dm_order2;
return dm;
}
/**
* @brief Normalisation of Schechter function for probability function
*/
double schechternorm(double mmin, double mmax, double slope, double mstar) {
/* integration resolution */
const int nintsteps = 10000;
/* starting mass */
double mc = mmin;
double norm = 0.;
/* Constant integration step in log space */
const double dm = pow(mmax / mmin, 1. / (double)nintsteps) - 1.;
/* integrate to determine norm (total number) */
for (int n = 0; n < nintsteps; n++) {
/* dm */
double dmc = mc * dm;
/* Leapfrog */
norm += pow(mc, -slope) * exp(-mc / mstar) * dmc -
0.5 * (pow(mc, -slope) / mstar + slope * pow(mc, -slope - 1.)) *
exp(-mc / mstar) * (dmc * dmc);
mc += dmc;
}
return norm;
}
/**
* @brief Mean mass of Schechter function
*
* Expected cluster mass \int_{mmin}^{mmax} M N(M) dM
*/
double meanmass(const double mmin, const double mmax, const double slope,
const double mstar, const double norm) {
/* integration resolution*/
const int nintsteps = 10000;
/* starting mass */
double mc = mmin;
double meanm = 0.;
/* Constant integration step in log space */
const double dm = pow(mmax / mmin, 1. / (double)nintsteps) - 1.;
/* integrate to determine norm (total number) */
for (int n = 0; n < nintsteps; n++) {
/* dm */
double dmc = mc * dm;
/* Leapfrog */
meanm +=
pow(mc, -slope + 1.) * exp(-mc / mstar) * dmc -
0.5 * (pow(mc, -slope + 1.) / mstar + (slope - 1.) * pow(mc, -slope)) *
exp(-mc / mstar) * (dmc * dmc);
mc += dmc;
}
return meanm / norm;
}
/**
* @brief Draw masses from Schechter function
*/
double schechtergen(const double mmin, const double mmax, const double slope,
const double mstar, const double norm, const long long sid,
const int cluster_index, const integertime_t ti_current) {
/* integration resolution */
const int nintsteps = 10000;
const double x = random_unit_interval_part_ID_and_index(
sid, cluster_index, ti_current, random_number_mosaic_schechter);
/* Method fails if x == 0 */
if (x == 0.) return mmin;
/* starting mass */
double mc = mmin;
double mcnorm = 0.;
double dmc = 0.;
double diff = 0.;
/* Constant integration step in log space */
const double dm = pow(mmax / mmin, 1. / (double)nintsteps) - 1.;
/* integrate until we reach the random number */
while (mcnorm / norm < x) {
/* dm */
dmc = mc * dm;
/* leapfrog */
diff = pow(mc, -slope) * exp(-mc / mstar) * dmc -
0.5 * (pow(mc, -slope) / mstar + slope * pow(mc, -slope - 1.)) *
exp(-mc / mstar) * (dmc * dmc);
mcnorm += diff;
mc += dmc;
}
/* Interpolate between integration steps to correct mass */
/* Previous step */
const double x1 = (mcnorm - diff) / norm;
/* Present step */
const double x2 = mcnorm / norm;
const double frac = (x - x1) / (x2 - x1);
/* interpolated mass */
return pow(mc - dmc, 1. - frac) * pow(mc, frac);
}
/**
* @brief Draw masses from power-law
*/
double powerlawgen(const double mmin, const double mmax, const double slope,
const long long sid, const int cluster_index,
const integertime_t ti_current) {
const double x = random_unit_interval_part_ID_and_index(
sid, cluster_index, ti_current, random_number_mosaic_powerlaw);
return pow(x * (pow(mmax, 1. - slope) - pow(mmin, 1. - slope)) +
pow(mmin, 1. - slope),
1. / (1. - slope));
}
/**
* @brief Do star cluster evolution for the mosaics subgrid model
*
* @param sp The #spart to act upon.
* @param e The #engine.
* @param with_cosmology Are we running with cosmological time integration.
*/
void mosaics_clevo(struct spart* restrict sp, const struct engine* e,
const struct cosmology* cosmo, const int with_cosmology) {
const struct stars_props* props = e->stars_properties;
const struct phys_const* phys_const = e->physical_constants;
/* Proportional change in mass from stellar evolution */
const float stellar_evo_factor = sp->mass / sp->mass_prev_timestep;
/* Total initial number of clusters over all model variations */
int ninit_sum = 0;
/* Quick check to see if we need to bother going further */
for (int iMP = 0; iMP < MOSAICS_MULTIPHYS; iMP++)
ninit_sum += sp->initial_num_clusters_evo[iMP];
if (ninit_sum == 0) {
/* No clusters formed in this particle, just handle the field component */
for (int iMP = 0; iMP < MOSAICS_MULTIPHYS; iMP++)
sp->field_mass[iMP] *= stellar_evo_factor;
/* For stellar evo. mass loss at the next timestep */
sp->mass_prev_timestep = sp->mass;
return;
}
/* Get the time since last step, then convert to physical units */
const integertime_t ti_step = get_integer_timestep(sp->time_bin);
const integertime_t ti_begin =
get_integer_time_begin(e->ti_current - 1, sp->time_bin);
/* Calculate age of the star at current time */
double star_age;
if (with_cosmology) {
star_age = cosmology_get_delta_time_from_scale_factors(
cosmo, (double)sp->birth_scale_factor, cosmo->a);
} else {
star_age = e->time - (double)sp->birth_time;
}
/* Get particle time-step */
double dt_star;
if (with_cosmology) {
dt_star =
cosmology_get_delta_time(e->cosmology, ti_begin, ti_begin + ti_step);
} else {
dt_star = get_timestep(sp->time_bin, e->time_base);
}
/* Evaporation timescale constant */
double t0_evap;
if (props->fixed_t0evap > 0) {
/* Use a fixed value? */
t0_evap = props->fixed_t0evap;
} else {
/* Get t0 from local tidal field */
t0_evap = get_local_evaporation_t0(sp, props, cosmo->a3_inv);
}
/* Calculate the tidal heating parameters for this timestep */
double Mlheatsum[6] = {0.};
calc_heating(sp, star_age, dt_star, e->time, cosmo->a3_inv, Mlheatsum);
/* Run multiple cluster formation models in parallel */
for (int iMP = 0; iMP < MOSAICS_MULTIPHYS; iMP++) {
/* Initial number of clusters for this model, up to max array length */
const int ninit =
min(sp->initial_num_clusters_evo[iMP], MOSAICS_MAX_CLUSTERS);
if (ninit == 0) {
/* No clusters were formed, just handle the field */
sp->field_mass[iMP] *= stellar_evo_factor;
continue;
}
/* Any surviving clusters? */
if (sp->num_clusters[iMP] > 0) {
/* Evolve the clusters */
for (int i = 0; i < ninit; i++) {
if (sp->cl_mass[iMP][i] == 0.) continue;
const double mass = (double)sp->cl_mass[iMP][i];
double dmevap = 0., dmsh = 0.;
/* Mass loss by tidal shocks */
if (props->cluster_shocks)
dmsh = mass_loss_shocks(sp, props, Mlheatsum, mass,
phys_const->const_newton_G);
/* Mass loss by evaporation */
if (props->cluster_evap) {
const double initial_mass = (double)sp->cl_initial_mass[iMP][i];
/* Stellar particle metallicity */
const double logZ =
log10(sp->chemistry_data.metal_mass_fraction_total + FLT_MIN);
dmevap = mass_loss_evaporation(
props, mass, initial_mass, t0_evap, dt_star, logZ,
phys_const->const_newton_G, phys_const->const_solar_mass);
}
/* Total mass loss for this cluster */
double dmdis = dmevap + dmsh;
/* Update cluster mass */
if (mass + dmdis < props->clMF_min) {
/* Cluster went below minimum mass limit, consider fully disrupted */
/* fraction of mass loss */
dmevap *= -mass / dmdis;
dmsh *= -mass / dmdis;
/* new mass loss (can't lose more than it has) */
dmdis = -mass;
/* We're done with this one */
sp->cl_mass[iMP][i] = 0.;
sp->num_clusters[iMP]--;
/* time of cluster disruption */
if (with_cosmology) {
sp->cl_disruption_time[iMP][i] = cosmo->a;
} else {
sp->cl_disruption_time[iMP][i] = e->time;
}
} else {
/* Cluster still survives, for now... */
sp->cl_mass[iMP][i] += (float)dmdis;
}
#ifdef SWIFT_DEBUG_CHECKS
if (dmevap > 0 || dmsh > 0) {
error("[mosaics] Positive mass loss?? dmevp = %g; dmsh = %g", dmevap,
dmsh);
}
#endif
/* Cluster mass lost to field */
sp->cl_dmshock[iMP][i] += (float)dmsh;
sp->field_mass[iMP] += (float)fabs(dmdis);
/* sp->cl_dmevap[iMP][i] += (float)dmevap; */
/* Stellar evolutionary mass loss for clusters */
sp->cl_mass[iMP][i] *= stellar_evo_factor;
} /* Cluster loop */
} /* surviving clusters */
/* Apply stellar evolution to the 'field' populations */
sp->field_mass[iMP] *= stellar_evo_factor;
for (int i = 0; i < ninit; i++) {
sp->cl_dmshock[iMP][i] *= stellar_evo_factor;
/* sp->cl_dmevap[iMP][i] *= stellar_evo_factor; */
}
} /* End MOSAICS_MULTIPHYS loop */
/* For stellar evo. mass loss at the next timestep */
sp->mass_prev_timestep = sp->mass;
}
/**
* @brief Do the cluster formation for the mosaics subgrid star cluster model
*
* @param sp The particle to act upon
* @param stars_properties Properties of the stars model.
* @param sf_props the star formation law properties to use
* @param phys_const the physical constants in internal units.
*/
void mosaics_clform(struct spart* restrict sp, const struct engine* e,
const struct cosmology* cosmo) {
#if defined(STAR_FORMATION_NONE)
error("Invalid SF model used!");
#endif
const struct stars_props* props = e->stars_properties;
const struct phys_const* phys_const = e->physical_constants;
const double const_G = phys_const->const_newton_G;
/* The resolved velocity dispersion (1/sqrt(3) converts to 1D assuming
* isotropy) and add the subgrid thermal component.
* Note: birth_velocity_dispersion is sigma^2 */
const double total_vel_disp2 =
sp->sf_data.birth_velocity_dispersion * (1. / 3.) +
sp->sound_speed_subgrid * sp->sound_speed_subgrid;
/* -------- Get CFE -------- */
/* The local gas velocity dispersion (1D) */
float sigmaloc;
if (props->subgrid_gas_vel_disp) {
/* "Sub-particle turbulent velocity dispersion" (E-MOSAICS version) */
sigmaloc = sqrt(sp->birth_pressure / sp->birth_density);
} else {
/* The resolved velocity dispersion */
sigmaloc = sqrt(total_vel_disp2);
}
/* Sound speed of cold ISM. Fixed or subgrid value? */
float csloc;
if (props->fixed_cs > 0) {
csloc = props->fixed_cs;
} else {
csloc = sp->sound_speed_subgrid;
}
/* Stellar particle metallicity */
const double logZ =
log10(sp->chemistry_data.metal_mass_fraction_total + FLT_MIN);
/* Calculate the cluster formation efficiency based on local conditions
* (Kruijssen 2012). */
const float CFE_local =
f_cfelocal(sp->birth_density, sigmaloc, csloc, sp->cell_gas_fraction,
logZ, props, phys_const);
/* -------- Get Mcstar -------- */
/* The 'GMC' mass */
double M_collapse;
if (props->Mgmc_Larson_relation) {
/* Use inverted Larson-type M-sigma relation to get GMC mass */
M_collapse =
props->GMC_m_sigma_const * pow(sigmaloc, props->GMC_m_sigma_index);
} else {
/* E-MOSAICS version: Toomre mass via tidal tensors (Pfeffer+18) */
/* Temporary array for eigvec/val calculation
* Note the lower elements are never referenced by dsyevj3 */
double tensor[3][3];
tensor[0][0] = sp->tidal_tensor[2][0];
tensor[1][1] = sp->tidal_tensor[2][1];
tensor[2][2] = sp->tidal_tensor[2][2];
tensor[0][1] = sp->tidal_tensor[2][3];
tensor[0][2] = sp->tidal_tensor[2][4];
tensor[1][2] = sp->tidal_tensor[2][5];
/* Get the eigenvalues */
double tideval[3];
dsyevj3(tensor, tideval);
/* sort eigenvalues, tideval[2] (=lambda_1) is largest eigenvalue */
sort3(tideval);
/* Circular and epicyclic frequencies */
double Omega2;
if (props->Omega_is_lambda2 && tideval[2] < 0) {
/* Correct version of Omega for a spherical distribution, in principle.
* When tideval[2] > 0 we can't guarantee that tideval[1] will be < 0 */
Omega2 = fabs(-tideval[1]);
} else {
/* For axisymmetric system, derived in Pfeffer+18 (i.e. E-MOSAICS version)
*/
Omega2 = fabs(-tideval[0] - tideval[1] - tideval[2]) * (1. / 3.);
}
double kappa2 = fabs(3. * Omega2 - tideval[2]);
/* Factor out cosmological factors from birth values */
Omega2 *= cosmo->a3_inv;
kappa2 *= cosmo->a3_inv;
sp->Omega_birth = sqrt(Omega2);
sp->kappa_birth = sqrt(kappa2);
/* Pressure from turbulent and thermal contributions */
const double total_pressure = sp->birth_density * total_vel_disp2;
/* Gas surface density, assuming hydrostatic equilibrium */
const double surfDens =
sqrt(2. * sp->cell_gas_fraction * total_pressure / (M_PI * const_G));
/* The Toomre mass */
const double MTconst =
4. * M_PI * M_PI * M_PI * M_PI * M_PI * const_G * const_G;
sp->Toomre_mass =
MTconst * surfDens * surfDens * surfDens / (kappa2 * kappa2);
/* Total collapse time of Toomre volume */
const double tcollapse = sqrt(2. * M_PI / kappa2);
/* Feedback timescale (Kruijssen 2012) */
double tfb = feedback_timescale(sp->birth_density, sigmaloc, csloc, props,
phys_const);
/* Toomre collapse fraction (Reina-Campos & Kruijssen 2017) */
sp->frac_collapse = fmin(1., tfb / tcollapse);
sp->frac_collapse *= sp->frac_collapse; /* f^2 */
sp->frac_collapse *= sp->frac_collapse; /* f^4 */
M_collapse = sp->Toomre_mass * sp->frac_collapse;
}
/* ----- Set up the star cluster population for this particle ----- */
/* Run multiple cluster formation models in parallel */
for (int iMP = 0; iMP < MOSAICS_MULTIPHYS; iMP++) {
/* Initialise */
sp->initial_num_clusters[iMP] = 0;
sp->initial_num_clusters_evo[iMP] = 0;
sp->num_clusters[iMP] = 0;
sp->initial_cluster_mass_total[iMP] = 0.;
sp->initial_cluster_mass_evo[iMP] = 0.;
for (int i = 0; i < MOSAICS_MAX_CLUSTERS; i++) {
sp->cl_mass[iMP][i] = 0.;
sp->cl_initial_mass[iMP][i] = 0.;
sp->cl_dmshock[iMP][i] = 0.;
sp->cl_disruption_time[iMP][i] = -1.;
/* sp->cl_id[iMP][i] = -1; */
/* sp->cl_dmevap[iMP][i] = 0.f; */
}
/* The cluster formation efficiency for this formation model */
if (props->fixed_CFE[iMP] > 0)
sp->CFE[iMP] = props->fixed_CFE[iMP];
else
sp->CFE[iMP] = CFE_local;
/* Exponential truncation of mass function for this formation model */
if (props->power_law_clMF[iMP]) {
/* No truncation for pure power law */
sp->Mcstar[iMP] = 0.;
} else {
/* Fixed value or locally varying? */
if (props->fixed_Mcstar[iMP] > 0)
sp->Mcstar[iMP] = props->fixed_Mcstar[iMP];
else
sp->Mcstar[iMP] = props->GMC_SFE * sp->CFE[iMP] * M_collapse;
}
/* Power-law index of the mass function, M^(-alpha) */
const float clMF_slope = props->clMF_slope[iMP];
/* Cluster and field masses if we were conserving mass within particles */
const float part_clust_mass = sp->mass_init * sp->CFE[iMP];
sp->field_mass[iMP] = sp->mass_init - part_clust_mass;
/* Determine mean cluster mass by integrating mass function,
* along with the integral normalisation */
double mean_mass = 0.;
double norm = 0.;
if (props->power_law_clMF[iMP]) {
/* Power law. Treat some special cases first */
if (clMF_slope == 1.) {
norm = log(props->clMF_max) - log(props->clMF_min);
mean_mass = props->clMF_max - props->clMF_min;
} else if (clMF_slope == 2.) {
norm = 1. / props->clMF_min - 1. / props->clMF_max;
mean_mass = log(props->clMF_max) - log(props->clMF_min);
} else {
norm = (pow(props->clMF_max, 1. - clMF_slope) -
pow(props->clMF_min, 1. - clMF_slope)) /
(1. - clMF_slope);
mean_mass = (pow(props->clMF_min, 2. - clMF_slope) -
pow(props->clMF_max, 2. - clMF_slope)) /
(clMF_slope - 2.);
}
mean_mass /= norm;
} else {
/* Schechter cluster mass function */
/* If Mcstar < clMFmin we have undefined behaviour and just form all
* clusters at M=clMFmin */
if (sp->Mcstar[iMP] > props->clMF_min) {
/* Get the normalisation of the integral */
norm = schechternorm(props->clMF_min, props->clMF_max, clMF_slope,
sp->Mcstar[iMP]);
}
if (norm > 0) {
mean_mass = meanmass(props->clMF_min, props->clMF_max, clMF_slope,
sp->Mcstar[iMP], norm);
}
}
/* Number of clusters to try and form.
* Draw from Poisson distribution, so we can be resolution independent */
if (mean_mass > 0) {
const double lambda = round(part_clust_mass / mean_mass);
sp->initial_num_clusters[iMP] = random_poisson(
sp->id, lambda, e->ti_current, random_number_mosaic_poisson);
} else {
sp->initial_num_clusters[iMP] = 0;
}
/* Lowest mass cluster we've formed so far */
float min_formed_mass = FLT_MAX;
int min_mass_idx = -1;
/* Counter for our cluster arrays */
int iArr = 0;
/* draw cluster masses from mass function */
for (int i = 0; i < sp->initial_num_clusters[iMP]; i++) {
/* Try and form a cluster */
float mTry;
if (props->power_law_clMF[iMP]) {
/* Power law mass function */
mTry = powerlawgen(props->clMF_min, props->clMF_max, clMF_slope, sp->id,
i, e->ti_current);
} else {
/* Schechter cluster mass function */
mTry = schechtergen(props->clMF_min, props->clMF_max, clMF_slope,
sp->Mcstar[iMP], norm, sp->id, i, e->ti_current);
}
sp->initial_cluster_mass_total[iMP] += mTry;
#ifdef SWIFT_DEBUG_CHECKS
if (!isfinite(mTry)) {
error("[mosaics] mTry is not finite!! mTry=%g, Mcstar=%g, norm=%g",
mTry, sp->Mcstar[iMP], norm);
}
#endif
/* Lower than our limit to evolve clusters */
if (mTry < props->clMF_min_evolve) {
/* Add destroyed clusters to field */
sp->field_mass[iMP] += mTry;
continue;
}
sp->initial_cluster_mass_evo[iMP] += mTry;
/* Check if we've run out of space... */
if (iArr < MOSAICS_MAX_CLUSTERS) {
/* Still space left */
/* sp->cl_id[iMP][iArr] = iArr; */
sp->cl_mass[iMP][iArr] = mTry;
sp->cl_initial_mass[iMP][iArr] = mTry;
/* Lowest mass cluster we've formed so far */
if (mTry < min_formed_mass) {
min_formed_mass = mTry;
min_mass_idx = iArr;
}
} else {
/* Ok we've run out of space. Replace the least massive
* cluster in the array so we at least track the most massive ones */
if (mTry > min_formed_mass) {
/* Add the previous one to field */
sp->field_mass[iMP] += sp->cl_mass[iMP][min_mass_idx];
sp->cl_mass[iMP][min_mass_idx] = mTry;
sp->cl_initial_mass[iMP][min_mass_idx] = mTry;
/* reset the counters */
min_formed_mass = FLT_MAX;
min_mass_idx = -1;
for (int j = 0; j < MOSAICS_MAX_CLUSTERS; j++) {
if (sp->cl_mass[iMP][j] < min_formed_mass) {
min_formed_mass = sp->cl_mass[iMP][j];
min_mass_idx = j;
}
}
} else {
/* We're out of space... just add to field */
sp->field_mass[iMP] += mTry;
}
}
iArr++;
}
/* Current number of clusters */
sp->num_clusters[iMP] = min(iArr, MOSAICS_MAX_CLUSTERS);
/* Number of clusters above the evolution mass limit */
sp->initial_num_clusters_evo[iMP] = iArr;
/* Warn when we hit the cluster array end */
if (sp->initial_num_clusters_evo[iMP] > MOSAICS_MAX_CLUSTERS) {
message(
"sp->id=%lld: Trying to form more clusters (%d) than allowed (%d)\n"
"MULTIPHYS=%d; mass_init=%g; clust_mass=%g; mean_mass=%g",
sp->id, sp->initial_num_clusters_evo[iMP], MOSAICS_MAX_CLUSTERS, iMP,
sp->mass_init, part_clust_mass, mean_mass);
}
} /* End MOSAICS_MULTIPHYS loop */
}