initial_mass_function.c 14.8 KB
Newer Older
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77
/*******************************************************************************
 * This file is part of SWIFT.
 * Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch)
 *
 * 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/>.
 *
 ******************************************************************************/

/* Include header */
#include "initial_mass_function.h"

/* local headers */
#include "hdf5_functions.h"
#include "stellar_evolution_struct.h"

/**
 * @brief Get the IMF exponent in between mass_min and mass_max.
 */
float initial_mass_function_get_exponent(
    const struct initial_mass_function *imf, float mass_min, float mass_max) {

#ifdef SWIFT_DEBUG_CHECKS
  if (mass_max > imf->mass_max)
    error("Cannot have mass larger than the largest one in the IMF");
  if (mass_min < imf->mass_min)
    error("Cannot have mass smaller than the smallest one in the IMF");
  if (mass_max < mass_min) error("Cannot have mass_min larger than mass_max");
#endif

  for (int i = 0; i < imf->n_parts; i++) {

    /* Check if in the correct part of the IMF */
    if (mass_min < imf->mass_limits[i + 1]) {

      /* Check if in only one segment */
      if (mass_max > imf->mass_limits[i + 1]) {
        error("Cannot get a single exponent for the interval [%g, %g]",
              mass_min, mass_max);
      }

      return imf->exp[i];
    }
  }

  error("Masses outside IMF ranges");

  return -1;
}

/** @brief Print the initial mass function */
void initial_mass_function_print(const struct initial_mass_function *imf) {

  message("Number of parts: %i", imf->n_parts);
  message("Mass interval: [%g, %g]", imf->mass_min, imf->mass_max);
  for (int i = 0; i < imf->n_parts; i++) {
    message("[%g, %g]: %.2g * m^{%g}", imf->mass_limits[i],
            imf->mass_limits[i + 1], imf->coef[i], imf->exp[i]);
  }
}

/**
 * @brief Integrate the #interpolation_1d data with the initial mass function.
 *
 * The x are supposed to be linear in log.
 *
 * @param imf The #initial_mass_function.
Loic Hausammann's avatar
Loic Hausammann committed
78 79 80 81
 * @param data The data to integrate.
 * @param count The number of element in data.
 * @param log_mass_min The value of the first element.
 * @param step_size The distance between two points.
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
82 83
 */
void initial_mass_function_integrate(const struct initial_mass_function *imf,
Loic Hausammann's avatar
Loic Hausammann committed
84 85
                                     float *data, size_t count,
                                     float log_mass_min, float step_size) {
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
86 87

  /* Index in the data */
Loic Hausammann's avatar
Loic Hausammann committed
88
  size_t j = 1;
89 90
  const float mass_min = exp10(log_mass_min);
  const float mass_max = exp10(log_mass_min + (count - 1) * step_size);
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
91 92 93

  float m = mass_min;

Loic Hausammann's avatar
Loic Hausammann committed
94
  float *tmp = (float *)malloc(sizeof(float) * count);
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110

  /* Set lower limit */
  tmp[0] = 0;
  for (int i = 0; i < imf->n_parts; i++) {

    /* Check if already in the correct part */
    if (mass_min > imf->mass_limits[i + 1]) {
      continue;
    }

    /* Check if already above the maximal mass */
    if (mass_max < imf->mass_limits[i]) {
      break;
    }

    /* Integrate the data */
Loic Hausammann's avatar
Loic Hausammann committed
111 112
    while ((m < imf->mass_limits[i + 1] || i == imf->n_parts - 1) &&
           j < count) {
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
113 114

      /* Compute the masses */
Loic Hausammann's avatar
Loic Hausammann committed
115
      const float log_m1 = log_mass_min + (j - 1) * step_size;
116
      const float m1 = exp10(log_m1);
Loic Hausammann's avatar
Loic Hausammann committed
117
      const float log_m2 = log_mass_min + j * step_size;
118
      const float m2 = exp10(log_m2);
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
119 120 121 122 123 124 125 126 127 128 129 130
      const float dm = m2 - m1;
      const float imf_1 = imf->coef[i] * pow(m1, imf->exp[i]);

      /* Get the imf of the upper limit  */
      float imf_2;
      if (m2 > imf->mass_limits[i + 1]) {
        imf_2 = imf->coef[i + 1] * pow(m2, imf->exp[i + 1]);
      } else {
        imf_2 = imf->coef[i] * pow(m2, imf->exp[i]);
      }

      /* Compute the integral */
Loic Hausammann's avatar
Loic Hausammann committed
131
      tmp[j] = tmp[j - 1] + 0.5 * (imf_1 * data[j - 1] + imf_2 * data[j]) * dm;
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
132 133 134 135 136 137 138 139

      /* Update j and m */
      j += 1;
      m = m2;
    }
  }

  /* The rest is extrapolated with 0 */
Loic Hausammann's avatar
Loic Hausammann committed
140
  for (size_t k = j; k < count; k++) {
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
141 142 143 144
    tmp[k] = tmp[k - 1];
  }

  /* Copy temporary array */
Loic Hausammann's avatar
Loic Hausammann committed
145
  memcpy(data, tmp, count * sizeof(float));
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195

  /* clean everything */
  free(tmp);
}

/**
 * @brief Get the IMF coefficient in between mass_min and mass_max.
 *
 * @param imf The #initial_mass_function.
 * @param mass_min The minimal mass of the requested interval.
 * @param mass_max The maximal mass of the requested interval.
 *
 * @return The imf's coefficient of the interval.
 */
float initial_mass_function_get_coefficient(
    const struct initial_mass_function *imf, float mass_min, float mass_max) {

  for (int i = 0; i < imf->n_parts; i++) {

    /* Check if in the correct part of the IMF */
    if (mass_min < imf->mass_limits[i + 1]) {

      /* Check if in only one segment */
      if (mass_max > imf->mass_limits[i + 1]) {
        error("Cannot get a single coefficient for the interval [%g, %g]",
              mass_min, mass_max);
      }

      return imf->coef[i];
    }
  }

  error("Masses outside IMF ranges");

  return -1;
}

/**
 * @brief Compute the integral of the fraction number of the initial mass
 * function.
 *
 * @param imf The #initial_mass_function.
 * @param m1 The lower mass to evaluate.
 * @param m2 The upper mass to evaluate.
 *
 * @return The number fraction.
 */
float initial_mass_function_get_integral_xi(
    const struct initial_mass_function *imf, float m1, float m2) {

196 197 198 199 200 201 202
  /* Ensure the masses to be withing the limits */
  m1 = min(m1, imf->mass_max);
  m1 = max(m1, imf->mass_min);

  m2 = min(m2, imf->mass_max);
  m2 = max(m2, imf->mass_min);

Loic Hausammann's avatar
Gear  
Loic Hausammann committed
203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268
  int k = -1;
  /* Find the correct part */
  for (int i = 0; i < imf->n_parts; i++) {
    if (m1 <= imf->mass_limits[i + 1]) {
      k = i;
      break;
    }
  }

  /* Check if found a part */
  if (k == -1) {
    error("Failed to find the correct function part: %g %g", m1, m2);
  }

  /* Check if m2 is inside the part */
  if (m2 < imf->mass_limits[k] || m2 > imf->mass_limits[k + 1]) {
    error("This function is not able to integrate in two different parts %g %g",
          m1, m2);
  }

  /* Compute the integral */
  const float int_xi1 = pow(m1, imf->exp[k]);
  const float int_xi2 = pow(m2, imf->exp[k]);

  return imf->coef[k] * (int_xi2 - int_xi1) / imf->exp[k];
};

/**
 * @brief Compute the mass fraction of the initial mass function.
 *
 * @param imf The #initial_mass_function.
 * @param m The mass to evaluate.
 *
 * @return The mass fraction.
 */
float initial_mass_function_get_imf(const struct initial_mass_function *imf,
                                    float m) {

#ifdef SWIFT_DEBUG_CHECKS
  if (m > imf->mass_max || m < imf->mass_min)
    error("Mass below or above limits expecting %g < %g < %g.", imf->mass_min,
          m, imf->mass_max);
#endif

  for (int i = 0; i < imf->n_parts; i++) {
    if (m <= imf->mass_limits[i + 1]) {
      return imf->coef[i] * pow(m, imf->exp[i]);
    }
  }

  error("Failed to find correct function part: %g larger than mass max %g.", m,
        imf->mass_max);
  return 0.;
};

/**
 * @brief Compute the integral of the mass fraction of the initial mass
 * function.
 *
 * @param imf The #initial_mass_function.
 * @param m1 The lower mass to evaluate.
 * @param m2 The upper mass to evaluate.
 *
 * @return The integral of the mass fraction.
 */
float initial_mass_function_get_integral_imf(
269
    const struct initial_mass_function *imf, float m1, float m2) {
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
270

271 272 273 274 275 276
  /* Ensure the masses to be withing the limits */
  m1 = min(m1, imf->mass_max);
  m1 = max(m1, imf->mass_min);

  m2 = min(m2, imf->mass_max);
  m2 = max(m2, imf->mass_min);
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335

  for (int i = 0; i < imf->n_parts; i++) {
    if (m1 <= imf->mass_limits[i + 1]) {
      if (m2 < imf->mass_limits[i] || m2 > imf->mass_limits[i + 1]) {
        error(
            "The code does not support the integration over multiple parts of "
            "the IMF");
      }
      const float exp = imf->exp[i] + 1.;
      return imf->coef[i] * (pow(m2, exp) - pow(m1, exp)) / exp;
    }
  }

  error("Failed to find correct function part: %g, %g larger than mass max %g.",
        m1, m2, imf->mass_max);
  return 0.;
};

/**
 * @brief Compute the coefficients of the initial mass function.
 *
 * @param imf The #initial_mass_function.
 */
void initial_mass_function_compute_coefficients(
    struct initial_mass_function *imf) {

  /* Allocate memory */
  if ((imf->coef = (float *)malloc(sizeof(float) * imf->n_parts)) == NULL)
    error("Failed to allocate the IMF coefficients.");

  /* Suppose that the first coefficients is 1 (will be corrected later) */
  imf->coef[0] = 1.;

  /* Use the criterion of continuity for the IMF */
  for (int i = 1; i < imf->n_parts; i++) {
    float exp = imf->exp[i - 1] - imf->exp[i];
    imf->coef[i] = imf->coef[i - 1] * pow(imf->mass_limits[i], exp);
  }

  /* Use the criterion on the integral = 1 */
  float integral = 0;
  for (int i = 0; i < imf->n_parts; i++) {
    const float exp = imf->exp[i] + 1.;
    const float m_i = pow(imf->mass_limits[i], exp);
    const float m_i1 = pow(imf->mass_limits[i + 1], exp);
    integral += imf->coef[i] * (m_i1 - m_i) / exp;
  }

  /* Normalize the coefficients (fix initial supposition) */
  for (int i = 0; i < imf->n_parts; i++) {
    imf->coef[i] /= integral;
  }
}

/**
 * @brief Reads the initial mass function parameters from the tables.
 *
 * @param imf The #initial_mass_function.
 * @param params The #swift_params.
Loic Hausammann's avatar
Loic Hausammann committed
336
 * @param filename The filename of the chemistry table.
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
337 338
 */
void initial_mass_function_read_from_table(struct initial_mass_function *imf,
Loic Hausammann's avatar
Loic Hausammann committed
339 340
                                           struct swift_params *params,
                                           const char *filename) {
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
341 342 343 344

  hid_t file_id, group_id;

  /* Open IMF group */
345
  h5_open_group(filename, "Data/IMF", &file_id, &group_id);
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390

  /* Read number of parts */
  io_read_attribute(group_id, "n", INT, &imf->n_parts);

  /* The tables have a different definition of n */
  imf->n_parts += 1;

  /* Allocate the memory for the exponents */
  if ((imf->exp = (float *)malloc(sizeof(float) * imf->n_parts)) == NULL)
    error("Failed to allocate the IMF exponents.");

  /* Read the exponents */
  io_read_array_attribute(group_id, "as", FLOAT, imf->exp, imf->n_parts);

  /* Allocate the memory for the temporary mass limits */
  if ((imf->mass_limits =
           (float *)malloc(sizeof(float) * (imf->n_parts + 1))) == NULL)
    error("Failed to allocate the IMF masses.");

  /* Read the mass limits */
  io_read_array_attribute(group_id, "ms", FLOAT, imf->mass_limits,
                          imf->n_parts - 1);

  /* Copy the data (need to shift for mass_min) */
  for (int i = imf->n_parts - 1; i > 0; i--) {
    imf->mass_limits[i] = imf->mass_limits[i - 1];
  }

  /* Read the minimal mass limit */
  io_read_attribute(group_id, "Mmin", FLOAT, &imf->mass_limits[0]);

  /* Read the maximal mass limit */
  io_read_attribute(group_id, "Mmax", FLOAT, &imf->mass_limits[imf->n_parts]);

  /* Close everything */
  h5_close_group(file_id, group_id);
}

/**
 * @brief Initialize the initial mass function.
 *
 * @param imf The #initial_mass_function.
 * @param phys_const The #phys_const.
 * @param us The #unit_system.
 * @param params The #swift_params.
Loic Hausammann's avatar
Loic Hausammann committed
391
 * @param filename The filename of the chemistry table.
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
392 393 394 395
 */
void initial_mass_function_init(struct initial_mass_function *imf,
                                const struct phys_const *phys_const,
                                const struct unit_system *us,
Loic Hausammann's avatar
Loic Hausammann committed
396 397
                                struct swift_params *params,
                                const char *filename) {
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
398 399

  /* Read the parameters from the yields table */
Loic Hausammann's avatar
Loic Hausammann committed
400
  initial_mass_function_read_from_table(imf, params, filename);
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497

  /* Write the masses in the correct attributes */
  imf->mass_min = imf->mass_limits[0];
  imf->mass_max = imf->mass_limits[imf->n_parts];

  /* Compute the coefficients */
  initial_mass_function_compute_coefficients(imf);
}

/**
 * @brief Write a initial_mass_function struct to the given FILE as a stream of
 * bytes.
 *
 * Here we are only writing the arrays, everything else has been copied in the
 * feedback.
 *
 * @param imf the struct
 * @param stream the file stream
 * @param sm The #stellar_model.
 */
void initial_mass_function_dump(const struct initial_mass_function *imf,
                                FILE *stream, const struct stellar_model *sm) {

  /* Dump the mass limits. */
  if (imf->mass_limits != NULL) {
    restart_write_blocks((void *)imf->mass_limits, sizeof(float),
                         imf->n_parts + 1, stream, "imf_mass_limits",
                         "imf_mass_limits");
  }

  /*! Dump the exponents. */
  if (imf->exp != NULL) {
    restart_write_blocks((void *)imf->exp, sizeof(float), imf->n_parts, stream,
                         "imf_exponents", "imf_exponents");
  }

  /*! Dump the coefficients. */
  if (imf->coef != NULL) {
    restart_write_blocks((void *)imf->coef, sizeof(float), imf->n_parts, stream,
                         "imf_coef", "imf_coef");
  }
}

/**
 * @brief Restore a initial_mass_function struct from the given FILE as a stream
 * of bytes.
 *
 * Here we are only writing the arrays, everything else has been copied in the
 * feedback.
 *
 * @param imf the struct
 * @param stream the file stream
 * @param sm The #stellar_model.
 */
void initial_mass_function_restore(struct initial_mass_function *imf,
                                   FILE *stream,
                                   const struct stellar_model *sm) {

  /* Restore the mass limits */
  if (imf->mass_limits != NULL) {
    imf->mass_limits = (float *)malloc(sizeof(float) * imf->n_parts + 1);
    restart_read_blocks((void *)imf->mass_limits, sizeof(float),
                        imf->n_parts + 1, stream, NULL, "imf_mass_limits");
  }

  /* Restore the exponents */
  if (imf->exp != NULL) {
    imf->exp = (float *)malloc(sizeof(float) * imf->n_parts);
    restart_read_blocks((void *)imf->exp, sizeof(float), imf->n_parts, stream,
                        NULL, "imf_exponents");
  }

  /* Restore the coefficients */
  if (imf->coef != NULL) {
    imf->coef = (float *)malloc(sizeof(float) * imf->n_parts);
    restart_read_blocks((void *)imf->coef, sizeof(float), imf->n_parts, stream,
                        NULL, "imf_coef");
  }
}

/**
 * @brief Clean the allocated memory.
 *
 * @param imf the #initial_mass_function.
 */
void initial_mass_function_clean(struct initial_mass_function *imf) {

  /* Free the pointers */
  free(imf->mass_limits);
  imf->mass_limits = NULL;

  free(imf->exp);
  imf->exp = NULL;

  free(imf->coef);
  imf->coef = NULL;
}