cosmology.c 22.1 KB
Newer Older
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
/*******************************************************************************
 * This file is part of SWIFT.
 * Copyright (c) 2017 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/>.
 *
 ******************************************************************************/

/**
 *  @file cosmology.c
 *  @brief Functions relating cosmological parameters
 */

/* This object's header. */
#include "cosmology.h"

/* Some standard headers */
#include <math.h>

31
/* Local headers */
32
#include "adiabatic_index.h"
33
#include "common_io.h"
34
#include "inline.h"
35
#include "restart.h"
36
37
38
39
40

#ifdef HAVE_LIBGSL
#include <gsl/gsl_integration.h>
#endif

41
/*! Number of values stored in the cosmological interpolation tables */
42
const int cosmology_table_length = 10000;
43

44
#ifdef HAVE_LIBGSL
45
/*! Size of the GSL workspace */
46
const size_t GSL_workspace_size = 100000;
47
#endif
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
/**
 * @brief Returns the interpolated value from a table.
 *
 * Uses linear interpolation.
 *
 * @brief table The table of value to interpolate from (should be of length
 * cosmology_table_length).
 * @brief x The value to interpolate at.
 * @brief x_min The mininum of the range of x.
 * @brief x_max The maximum of the range of x.
 */
static INLINE double interp_table(double *table, double x, double x_min,
                                  double x_max) {

  const double xx = ((x - x_min) / (x_max - x_min)) * cosmology_table_length;
  const int i = (int)xx;
  const int ii = (i >= cosmology_table_length) ? cosmology_table_length - 1 : i;

  if (ii <= 1)
    return table[0] * xx;
  else
    return table[ii - 1] + (table[ii] - table[ii - 1]) * (xx - ii);
}

73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
/**
 * @brief Computes the dark-energy equation of state at a given scale-factor a.
 *
 * We follow the convention of Linder & Jenkins, MNRAS, 346, 573, 2003
 *
 * @param a The current scale-factor
 * @param w_0 The equation of state parameter at z=0
 * @param w_a The equation of state evolution parameter
 */
static INLINE double cosmology_dark_energy_EoS(double a, double w_0,
                                               double w_a) {

  return w_0 + w_a * (1. - a);
}

88
/**
89
 * @brief Computes the integral of the dark-energy equation of state
90
91
92
93
94
95
 * up to a scale-factor a.
 *
 * We follow the convention of Linder & Jenkins, MNRAS, 346, 573, 2003
 * and compute \f$ \tilde{w}(a) = \int_0^a\frac{1 + w(z)}{1+z}dz \f$.
 *
 * @param a The current scale-factor.
Matthieu Schaller's avatar
Matthieu Schaller committed
96
97
 * @param w0 The equation of state parameter at z=0
 * @param wa The equation of state evolution parameter
98
 */
99
100
101
static INLINE double w_tilde(double a, double w0, double wa) {
  return (a - 1.) * wa - (1. + w0 + wa) * log(a);
}
102

103
104
105
106
/**
 * @brief Compute \f$ E(z) \f$.
 */
static INLINE double E(double Or, double Om, double Ok, double Ol, double w0,
107
                       double wa, double a) {
108
  const double a_inv = 1. / a;
109
  return sqrt(Or * a_inv * a_inv * a_inv * a_inv + Om * a_inv * a_inv * a_inv +
110
              Ok * a_inv * a_inv + Ol * exp(3. * w_tilde(a, w0, wa)));
111
112
}

113
114
115
116
117
118
119
120
/**
 * @brief Returns the time (in internal units) since Big Bang at a given
 * scale-factor.
 *
 * @param c The current #cosmology.
 * @param a Scale-factor of interest.
 */
double cosmology_get_time_since_big_bang(const struct cosmology *c, double a) {
121

122
123
124
  /* Time between a_begin and a */
  const double delta_t =
      interp_table(c->time_interp_table, log(a), c->log_a_begin, c->log_a_end);
125

126
  return c->time_interp_table_offset + delta_t;
127
128
}

129
130
131
132
/**
 * @brief Update the cosmological parameters to the current simulation time.
 *
 * @param c The #cosmology struct.
133
 * @param ti_current The current (integer) time.
134
 */
135
void cosmology_update(struct cosmology *c, integertime_t ti_current) {
136

137
  /* Get scale factor and powers of it */
138
  const double a = c->a_begin * exp(ti_current * c->time_base);
139
140
141
  const double a_inv = 1. / a;
  c->a = a;
  c->a_inv = a_inv;
142
  c->a2_inv = a_inv * a_inv;
143
  c->a3_inv = a_inv * a_inv * a_inv;
144
145
146
  c->a_factor_internal_energy =
      pow(a, -3. * hydro_gamma_minus_one);          /* a^{3*(1-gamma)} */
  c->a_factor_pressure = pow(a, -3. * hydro_gamma); /* a^{-3*gamma} */
147
  c->a_factor_sound_speed =
148
      pow(a, -1.5 * hydro_gamma_minus_one); /* a^{3*(1-gamma)/2} */
149
150
  c->a_factor_grav_accel = a_inv * a_inv;   /* 1 / a^2 */
  c->a_factor_hydro_accel =
151
152
153
      pow(a, -3. * hydro_gamma + 2.); /* 1 / a^(3*gamma - 2) */
  c->a_factor_mu =
      pow(a, 0.5 * (3. * hydro_gamma - 5.)); /* a^{(3*gamma - 5) / 2} */
154
155

  /* Redshift */
Matthieu Schaller's avatar
Matthieu Schaller committed
156
  c->z = a_inv - 1.;
157

158
159
160
  /* Dark-energy equation of state */
  c->w = cosmology_dark_energy_EoS(a, c->w_0, c->w_a);

161
162
163
  /* E(z) */
  const double Omega_r = c->Omega_r;
  const double Omega_m = c->Omega_m;
164
  const double Omega_k = c->Omega_k;
165
  const double Omega_l = c->Omega_lambda;
166
167
168
  const double w0 = c->w_0;
  const double wa = c->w_a;
  const double E_z = E(Omega_r, Omega_m, Omega_k, Omega_l, w0, wa, a);
169
170
171
172

  /* H(z) */
  c->H = c->H0 * E_z;

173
174
175
  /* Expansion rate */
  c->a_dot = c->H * c->a;

176
177
178
  /* Time-step conversion factor */
  c->time_step_factor = c->H;

179
  /* Time */
180
  c->time = cosmology_get_time_since_big_bang(c, a);
181
  c->lookback_time = c->universe_age_at_present_day - c->time;
182
183
}

184
185
186
187
/**
 * @brief Computes \f$ dt / a^2 \f$ for the current cosmology.
 *
 * @param a The scale-factor of interest.
188
 * @param param The current #cosmology.
189
190
 */
double drift_integrand(double a, void *param) {
191

Matthieu Schaller's avatar
Matthieu Schaller committed
192
  const struct cosmology *c = (const struct cosmology *)param;
193
194
  const double Omega_r = c->Omega_r;
  const double Omega_m = c->Omega_m;
195
  const double Omega_k = c->Omega_k;
196
  const double Omega_l = c->Omega_lambda;
197
198
  const double w_0 = c->w_0;
  const double w_a = c->w_a;
199
  const double H0 = c->H0;
Matthieu Schaller's avatar
Matthieu Schaller committed
200

201
  const double a_inv = 1. / a;
202
  const double E_z = E(Omega_r, Omega_m, Omega_k, Omega_l, w_0, w_a, a);
203
  const double H = H0 * E_z;
Matthieu Schaller's avatar
Matthieu Schaller committed
204

205
206
  return (1. / H) * a_inv * a_inv * a_inv;
}
207

208
209
210
211
/**
 * @brief Computes \f$ dt / a \f$ for the current cosmology.
 *
 * @param a The scale-factor of interest.
212
 * @param param The current #cosmology.
213
214
 */
double gravity_kick_integrand(double a, void *param) {
Matthieu Schaller's avatar
Matthieu Schaller committed
215
216

  const struct cosmology *c = (const struct cosmology *)param;
217
218
219
220
221
222
223
  const double Omega_r = c->Omega_r;
  const double Omega_m = c->Omega_m;
  const double Omega_k = c->Omega_k;
  const double Omega_l = c->Omega_lambda;
  const double w_0 = c->w_0;
  const double w_a = c->w_a;
  const double H0 = c->H0;
Matthieu Schaller's avatar
Matthieu Schaller committed
224

225
  const double a_inv = 1. / a;
226
  const double E_z = E(Omega_r, Omega_m, Omega_k, Omega_l, w_0, w_a, a);
227
  const double H = H0 * E_z;
Matthieu Schaller's avatar
Matthieu Schaller committed
228

229
230
  return (1. / H) * a_inv * a_inv;
}
231

232
233
234
235
/**
 * @brief Computes \f$ dt / a^{3(\gamma - 1) + 1} \f$ for the current cosmology.
 *
 * @param a The scale-factor of interest.
236
 * @param param The current #cosmology.
237
238
 */
double hydro_kick_integrand(double a, void *param) {
Matthieu Schaller's avatar
Matthieu Schaller committed
239
240

  const struct cosmology *c = (const struct cosmology *)param;
241
242
243
244
245
246
247
  const double Omega_r = c->Omega_r;
  const double Omega_m = c->Omega_m;
  const double Omega_k = c->Omega_k;
  const double Omega_l = c->Omega_lambda;
  const double w_0 = c->w_0;
  const double w_a = c->w_a;
  const double H0 = c->H0;
Matthieu Schaller's avatar
Matthieu Schaller committed
248

249
  const double a_inv = 1. / a;
250
  const double E_z = E(Omega_r, Omega_m, Omega_k, Omega_l, w_0, w_a, a);
251
  const double H = H0 * E_z;
Matthieu Schaller's avatar
Matthieu Schaller committed
252

253
254
255
256
257
  /* Note: we can't use the pre-defined pow_gamma_xxx() function as
     as we need double precision accuracy for the GSL routine. */
  return (1. / H) * pow(a_inv, 3. * hydro_gamma_minus_one) * a_inv;
}

258
259
260
261
/**
 * @brief Computes \f$ dt \f$ for the current cosmology.
 *
 * @param a The scale-factor of interest.
262
 * @param param The current #cosmology.
263
264
 */
double time_integrand(double a, void *param) {
265

Matthieu Schaller's avatar
Matthieu Schaller committed
266
  const struct cosmology *c = (const struct cosmology *)param;
267
268
269
270
271
272
273
  const double Omega_r = c->Omega_r;
  const double Omega_m = c->Omega_m;
  const double Omega_k = c->Omega_k;
  const double Omega_l = c->Omega_lambda;
  const double w_0 = c->w_0;
  const double w_a = c->w_a;
  const double H0 = c->H0;
Matthieu Schaller's avatar
Matthieu Schaller committed
274

275
  const double a_inv = 1. / a;
276
  const double E_z = E(Omega_r, Omega_m, Omega_k, Omega_l, w_0, w_a, a);
277
  const double H = H0 * E_z;
Matthieu Schaller's avatar
Matthieu Schaller committed
278

279
280
  return (1. / H) * a_inv;
}
281

282
283
284
/**
 * @brief Initialise the interpolation tables for the integrals.
 */
285
void cosmology_init_tables(struct cosmology *c) {
286

287
  const ticks tic = getticks();
288

289
#ifdef HAVE_LIBGSL
290

291
292
  /* Retrieve some constants */
  const double a_begin = c->a_begin;
293

294
  /* Allocate memory for the interpolation tables */
295
296
  c->drift_fac_interp_table =
      (double *)malloc(cosmology_table_length * sizeof(double));
297
  c->grav_kick_fac_interp_table =
298
      (double *)malloc(cosmology_table_length * sizeof(double));
299
  c->hydro_kick_fac_interp_table =
300
301
302
      (double *)malloc(cosmology_table_length * sizeof(double));
  c->time_interp_table =
      (double *)malloc(cosmology_table_length * sizeof(double));
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317

  /* Prepare a table of scale factors for the integral bounds */
  const double delta_a =
      (c->log_a_end - c->log_a_begin) / cosmology_table_length;
  double *a_table = malloc(cosmology_table_length * sizeof(double));
  for (int i = 0; i < cosmology_table_length; i++)
    a_table[i] = exp(c->log_a_begin + delta_a * (i + 1));

  /* Initalise the GSL workspace */
  gsl_integration_workspace *space =
      gsl_integration_workspace_alloc(GSL_workspace_size);

  double result, abserr;

  /* Integrate the drift factor \int_{a_begin}^{a_table[i]} dt/a^2 */
318
  gsl_function F = {&drift_integrand, c};
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
  for (int i = 0; i < cosmology_table_length; i++) {
    gsl_integration_qag(&F, a_begin, a_table[i], 0, 1.0e-10, GSL_workspace_size,
                        GSL_INTEG_GAUSS61, space, &result, &abserr);

    /* Store result */
    c->drift_fac_interp_table[i] = result;
  }

  /* Integrate the kick factor \int_{a_begin}^{a_table[i]} dt/a */
  F.function = &gravity_kick_integrand;
  for (int i = 0; i < cosmology_table_length; i++) {
    gsl_integration_qag(&F, a_begin, a_table[i], 0, 1.0e-10, GSL_workspace_size,
                        GSL_INTEG_GAUSS61, space, &result, &abserr);

    /* Store result */
    c->grav_kick_fac_interp_table[i] = result;
  }

337
338
339
340
341
342
343
344
345
346
  /* Integrate the kick factor \int_{a_begin}^{a_table[i]} dt/a^(3(g-1)+1) */
  F.function = &hydro_kick_integrand;
  for (int i = 0; i < cosmology_table_length; i++) {
    gsl_integration_qag(&F, a_begin, a_table[i], 0, 1.0e-10, GSL_workspace_size,
                        GSL_INTEG_GAUSS61, space, &result, &abserr);

    /* Store result */
    c->hydro_kick_fac_interp_table[i] = result;
  }

347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
  /* Integrate the time \int_{a_begin}^{a_table[i]} dt */
  F.function = &time_integrand;
  for (int i = 0; i < cosmology_table_length; i++) {
    gsl_integration_qag(&F, a_begin, a_table[i], 0, 1.0e-10, GSL_workspace_size,
                        GSL_INTEG_GAUSS61, space, &result, &abserr);

    /* Store result */
    c->time_interp_table[i] = result;
  }

  /* Integrate the time \int_{0}^{a_begin} dt */
  gsl_integration_qag(&F, 0., a_begin, 0, 1.0e-10, GSL_workspace_size,
                      GSL_INTEG_GAUSS61, space, &result, &abserr);
  c->time_interp_table_offset = result;

362
363
364
365
366
  /* Integrate the time \int_{0}^{1} dt */
  gsl_integration_qag(&F, 0., 1, 0, 1.0e-13, GSL_workspace_size,
                      GSL_INTEG_GAUSS61, space, &result, &abserr);
  c->universe_age_at_present_day = result;

367
368
369
370
371
372
373
374
375
376
377
378
  /* Free the workspace and temp array */
  gsl_integration_workspace_free(space);
  free(a_table);

#else

  error("Code not compiled with GSL. Can't compute cosmology integrals.");

#endif

  message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
          clocks_getunit());
379
380
}

381
382
383
384
/**
 * @brief Initialises the #cosmology from the values read in the parameter file.
 *
 * @param params The parsed values.
385
 * @param us The current internal system of units.
386
387
388
 * @param phys_const The physical constants in the current system of units.
 * @param c The #cosmology to initialise.
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
389
390
391
void cosmology_init(const struct swift_params *params,
                    const struct unit_system *us,
                    const struct phys_const *phys_const, struct cosmology *c) {
392
393
394

  /* Read in the cosmological parameters */
  c->Omega_m = parser_get_param_double(params, "Cosmology:Omega_m");
395
  c->Omega_r = parser_get_opt_param_double(params, "Cosmology:Omega_r", 0.);
396
397
  c->Omega_lambda = parser_get_param_double(params, "Cosmology:Omega_lambda");
  c->Omega_b = parser_get_param_double(params, "Cosmology:Omega_b");
398
399
  c->w_0 = parser_get_opt_param_double(params, "Cosmology:w_0", -1.);
  c->w_a = parser_get_opt_param_double(params, "Cosmology:w_a", 0.);
400
  c->h = parser_get_param_double(params, "Cosmology:h");
Matthieu Schaller's avatar
Matthieu Schaller committed
401
402

  /* Read the start and end of the simulation */
403
404
  c->a_begin = parser_get_param_double(params, "Cosmology:a_begin");
  c->a_end = parser_get_param_double(params, "Cosmology:a_end");
405
406
  c->log_a_begin = log(c->a_begin);
  c->log_a_end = log(c->a_end);
407
408
  c->time_base = (c->log_a_end - c->log_a_begin) / max_nr_timesteps;
  c->time_base_inv = 1. / c->time_base;
409
410

  /* Construct derived quantities */
Matthieu Schaller's avatar
Matthieu Schaller committed
411
412

  /* Curvature density (for closure) */
413
  c->Omega_k = 1. - (c->Omega_m + c->Omega_r + c->Omega_lambda);
Matthieu Schaller's avatar
Matthieu Schaller committed
414

415
  /* Dark-energy equation of state */
416
  c->w = cosmology_dark_energy_EoS(c->a_begin, c->w_0, c->w_a);
417

Matthieu Schaller's avatar
Matthieu Schaller committed
418
  /* Hubble constant in internal units */
419
420
421
  const double km = 1.e5 / units_cgs_conversion_factor(us, UNIT_CONV_LENGTH);
  const double H0_cgs =
      100. * c->h * (km / (1.e6 * phys_const->const_parsec)); /* s^-1 */
Matthieu Schaller's avatar
Matthieu Schaller committed
422
  c->H0 = H0_cgs * units_cgs_conversion_factor(us, UNIT_CONV_TIME);
423
  c->Hubble_time = 1. / c->H0;
424

425
426
427
  /* Initialise the interpolation tables */
  c->drift_fac_interp_table = NULL;
  c->grav_kick_fac_interp_table = NULL;
428
  c->hydro_kick_fac_interp_table = NULL;
429
430
431
  c->time_interp_table = NULL;
  c->time_interp_table_offset = 0.;
  cosmology_init_tables(c);
432

433
434
435
  /* Set remaining variables to alid values */
  cosmology_update(c, 0);

436
437
438
  /* Update the times */
  c->time_begin = cosmology_get_time_since_big_bang(c, c->a_begin);
  c->time_end = cosmology_get_time_since_big_bang(c, c->a_end);
439
}
440

441
442
443
/**
 * @brief Initialise the #cosmology for non-cosmological time-integration
 *
444
 * Essentially sets all constants to 1 or 0.
445
446
447
 *
 * @param c The #cosmology to initialise.
 */
448
void cosmology_init_no_cosmo(struct cosmology *c) {
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464

  c->Omega_m = 0.;
  c->Omega_r = 0.;
  c->Omega_k = 0.;
  c->Omega_lambda = 0.;
  c->Omega_b = 0.;
  c->w_0 = 0.;
  c->w_a = 0.;
  c->h = 0.;
  c->w = 0.;

  c->a_begin = 1.;
  c->a_end = 1.;
  c->log_a_begin = 0.;
  c->log_a_end = 0.;

465
  c->H = 0.;
466
467
  c->a = 1.;
  c->a_inv = 1.;
468
  c->a2_inv = 1.;
469
  c->a3_inv = 1.;
470
471
  c->a_factor_internal_energy = 1.;
  c->a_factor_pressure = 1.;
472
  c->a_factor_sound_speed = 1.;
473
  c->a_factor_mu = 1.;
474
475
  c->a_factor_hydro_accel = 1.;
  c->a_factor_grav_accel = 1.;
476

477
478
  c->time_step_factor = 1.;

479
480
481
482
483
484
485
486
487
488
  c->a_dot = 0.;
  c->time = 0.;
  c->universe_age_at_present_day = 0.;

  /* Initialise the interpolation tables */
  c->drift_fac_interp_table = NULL;
  c->grav_kick_fac_interp_table = NULL;
  c->hydro_kick_fac_interp_table = NULL;
  c->time_interp_table = NULL;
  c->time_interp_table_offset = 0.;
489
490
491

  c->time_begin = 0.;
  c->time_end = 0.;
492
493
}

494
495
496
497
498
/**
 * @brief Computes the cosmology factor that enters the drift operator.
 *
 * Computes \f$ \int_{a_start}^{a_end} dt/a^2 \f$ using the interpolation table.
 *
499
 * @param c The current #cosmology.
500
501
502
 * @param ti_start the (integer) time of the start of the drift.
 * @param ti_end the (integer) time of the end of the drift.
 */
503
double cosmology_get_drift_factor(const struct cosmology *c,
Matthieu Schaller's avatar
Matthieu Schaller committed
504
505
                                  integertime_t ti_start,
                                  integertime_t ti_end) {
506

507
508
  const double a_start = c->log_a_begin + ti_start * c->time_base;
  const double a_end = c->log_a_begin + ti_end * c->time_base;
509

Matthieu Schaller's avatar
Matthieu Schaller committed
510
511
512
513
  const double int_start = interp_table(c->drift_fac_interp_table, a_start,
                                        c->log_a_begin, c->log_a_end);
  const double int_end = interp_table(c->drift_fac_interp_table, a_end,
                                      c->log_a_begin, c->log_a_end);
514
515
516
517
518
519
520
521
522

  return int_end - int_start;
}

/**
 * @brief Computes the cosmology factor that enters the gravity kick operator.
 *
 * Computes \f$ \int_{a_start}^{a_end} dt/a \f$ using the interpolation table.
 *
523
 * @param c The current #cosmology.
524
525
526
 * @param ti_start the (integer) time of the start of the drift.
 * @param ti_end the (integer) time of the end of the drift.
 */
527
double cosmology_get_grav_kick_factor(const struct cosmology *c,
Matthieu Schaller's avatar
Matthieu Schaller committed
528
529
                                      integertime_t ti_start,
                                      integertime_t ti_end) {
530

531
532
  const double a_start = c->log_a_begin + ti_start * c->time_base;
  const double a_end = c->log_a_begin + ti_end * c->time_base;
533

Matthieu Schaller's avatar
Matthieu Schaller committed
534
535
536
537
  const double int_start = interp_table(c->grav_kick_fac_interp_table, a_start,
                                        c->log_a_begin, c->log_a_end);
  const double int_end = interp_table(c->grav_kick_fac_interp_table, a_end,
                                      c->log_a_begin, c->log_a_end);
538
539
540
541
542
543
544
545
546

  return int_end - int_start;
}

/**
 * @brief Computes the cosmology factor that enters the hydro kick operator.
 *
 * Computes \f$ \int_{a_start}^{a_end} dt/a \f$ using the interpolation table.
 *
547
 * @param c The current #cosmology.
548
549
550
 * @param ti_start the (integer) time of the start of the drift.
 * @param ti_end the (integer) time of the end of the drift.
 */
551
double cosmology_get_hydro_kick_factor(const struct cosmology *c,
Matthieu Schaller's avatar
Matthieu Schaller committed
552
553
                                       integertime_t ti_start,
                                       integertime_t ti_end) {
554

555
556
  const double a_start = c->log_a_begin + ti_start * c->time_base;
  const double a_end = c->log_a_begin + ti_end * c->time_base;
557

558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
  const double int_start = interp_table(c->drift_fac_interp_table, a_start,
                                        c->log_a_begin, c->log_a_end);
  const double int_end = interp_table(c->drift_fac_interp_table, a_end,
                                      c->log_a_begin, c->log_a_end);

  return int_end - int_start;
}

/**
 * @brief Computes the cosmology factor that enters the thermal variable kick
 * operator.
 *
 * Computes \f$ \int_{a_start}^{a_end} dt/a^2 \f$ using the interpolation table.
 *
 * @param c The current #cosmology.
 * @param ti_start the (integer) time of the start of the drift.
 * @param ti_end the (integer) time of the end of the drift.
 */
576
double cosmology_get_therm_kick_factor(const struct cosmology *c,
577
578
579
580
581
582
                                       integertime_t ti_start,
                                       integertime_t ti_end) {

  const double a_start = c->log_a_begin + ti_start * c->time_base;
  const double a_end = c->log_a_begin + ti_end * c->time_base;

Matthieu Schaller's avatar
Matthieu Schaller committed
583
584
585
586
  const double int_start = interp_table(c->hydro_kick_fac_interp_table, a_start,
                                        c->log_a_begin, c->log_a_end);
  const double int_end = interp_table(c->hydro_kick_fac_interp_table, a_end,
                                      c->log_a_begin, c->log_a_end);
587
588
589

  return int_end - int_start;
}
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610

/**
 * @brief Compute the cosmic time (in internal units) between two scale-factors.
 *
 * @brief c The current #cosmology.
 * @brief a1 The first scale-factor.
 * @brief a2 The second scale-factor.
 */
double cosmology_get_delta_time(const struct cosmology *c, double a1,
                                double a2) {

  /* Time between a_begin and a1 */
  const double t1 =
      interp_table(c->time_interp_table, log(a1), c->log_a_begin, c->log_a_end);

  /* Time between a_begin and a1 */
  const double t2 =
      interp_table(c->time_interp_table, log(a2), c->log_a_begin, c->log_a_end);

  return t2 - t1;
}
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626

/**
 * @brief Prints the #cosmology model to stdout.
 */
void cosmology_print(const struct cosmology *c) {

  message(
      "Density parameters: [O_m, O_l, O_b, O_k, O_r] = [%f, %f, %f, %f, %f]",
      c->Omega_m, c->Omega_lambda, c->Omega_b, c->Omega_k, c->Omega_r);
  message("Dark energy equation of state: w_0=%f w_a=%f", c->w_0, c->w_a);
  message("Hubble constant: h = %f, H_0 = %e U_t^(-1)", c->h, c->H0);
  message("Hubble time: 1/H0 = %e U_t", c->Hubble_time);
  message("Universe age at present day: %e U_t",
          c->universe_age_at_present_day);
}

627
628
629
630
631
632
633
634
void cosmology_clean(struct cosmology *c) {

  free(c->drift_fac_interp_table);
  free(c->grav_kick_fac_interp_table);
  free(c->hydro_kick_fac_interp_table);
  free(c->time_interp_table);
}

635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
#ifdef HAVE_HDF5
void cosmology_write_model(hid_t h_grp, const struct cosmology *c) {

  io_write_attribute_d(h_grp, "a_beg", c->a_begin);
  io_write_attribute_d(h_grp, "a_end", c->a_end);
  io_write_attribute_d(h_grp, "time_beg [internal units]", c->time_begin);
  io_write_attribute_d(h_grp, "time_end [internal units]", c->time_end);
  io_write_attribute_d(h_grp, "h", c->h);
  io_write_attribute_d(h_grp, "H0 [internal units]", c->H0);
  io_write_attribute_d(h_grp, "Hubble time [internal units]", c->Hubble_time);
  io_write_attribute_d(h_grp, "Omega_m", c->Omega_m);
  io_write_attribute_d(h_grp, "Omega_r", c->Omega_r);
  io_write_attribute_d(h_grp, "Omega_b", c->Omega_b);
  io_write_attribute_d(h_grp, "Omega_k", c->Omega_k);
  io_write_attribute_d(h_grp, "Omega_lambda", c->Omega_lambda);
  io_write_attribute_d(h_grp, "w_0", c->w_0);
  io_write_attribute_d(h_grp, "w_a", c->w_a);
}
#endif

655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
/**
 * @brief Write a cosmology struct to the given FILE as a stream of bytes.
 *
 * @param cosmology the struct
 * @param stream the file stream
 */
void cosmology_struct_dump(const struct cosmology *cosmology, FILE *stream) {
  restart_write_blocks((void *)cosmology, sizeof(struct cosmology), 1, stream,
                       "cosmology", "cosmology function");
}

/**
 * @brief Restore a cosmology struct from the given FILE as a stream of
 * bytes.
 *
 * @param cosmology the struct
 * @param stream the file stream
 */
673
void cosmology_struct_restore(struct cosmology *cosmology, FILE *stream) {
674
675
  restart_read_blocks((void *)cosmology, sizeof(struct cosmology), 1, stream,
                      NULL, "cosmology function");
676
677
678

  /* Re-initialise the tables */
  cosmology_init_tables(cosmology);
679
}