starformation.h 14 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
/*******************************************************************************
 * This file is part of SWIFT.
 * Copyright (c) 2018 Folkert Nobels (nobels@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_SCHAYE_STARFORMATION_H
#define SWIFT_SCHAYE_STARFORMATION_H

23
24
25
/* Some standard headers */
#include <stdlib.h>

26
27
28
29
/* Local includes */
#include "cosmology.h"
#include "physical_constants.h"
#include "units.h"
Folkert Nobels's avatar
Folkert Nobels committed
30
#include "engine.h"
31
32
#include "parser.h"
#include "equation_of_state.h"
Folkert Nobels's avatar
Folkert Nobels committed
33
#include "part.h"
34
#include "hydro.h"
35
#include "cooling.h"
36
#include "adiabatic_index.h"
37
#include "cell.h" 
38
#include "stars.h"
39

40
41
42
/* Starformation struct */
struct star_formation {
  
43
  /*! Normalization of the KS star formation law */
44
  double KS_normalization;
45

46
  /*! Slope of the KS law */
47
48
49
  double KS_power_law;

  /*! Slope of the high density KS law */
Folkert Nobels's avatar
Folkert Nobels committed
50
  double KS_high_den_power_law;
51
52
53

  /*! KS law High density threshold */
  double KS_high_den_thresh;
54

Folkert Nobels's avatar
Folkert Nobels committed
55
56
57
  /*! KS high density normalization */
  double KS_high_den_normalization;

58
  /*! Critical overdensity */
Folkert Nobels's avatar
Folkert Nobels committed
59
  double min_over_den;
60

61
  /*! Critical temperature */
62
63
  double T_crit;

64
  /*! gas fraction */
65
  double fgas;
66

67
  /*! Star formation law slope */
68
  double SF_power_law;
69

70
  /*! star formation normalization of schaye+08 */
Folkert Nobels's avatar
Folkert Nobels committed
71
72
73
74
75
76
77
  double SF_normalization;

  /*! star formation high density slope */
  double SF_high_den_power_law;

  /*! Star formation high density normalization */
  double SF_high_den_normalization;
78
79
80

  /*! Inverse of RAND_MAX */
  double inv_RAND_MAX;
81
82
83
84

  /*! Critical density to form stars */
  double den_crit;

85
86
87
  /*! Maximum critical density to form stars */
  double den_crit_max;

88
89
90
  /*! Scaling metallicity */
  double Z0;

91
92
93
  /*! one over the scaling metallicity */
  double Z0_inv;

94
95
96
97
98
  /*! critical density Metallicity power law */
  double n_Z0;

  /*! Normalization of critical SF density of Schaye (2004) */
  double den_crit_star;
Folkert Nobels's avatar
Folkert Nobels committed
99
100
101

  /*! Metallicity dependent critical density from Schaye (2004) */
  int schaye04;
102
  
103
104
105
106
107
108
109
110
  /*! Polytropic index */
  double polytropic_index;

  /*! EOS pressure norm */
  double EOS_pressure_norm;

  /*! EOS density norm */
  double EOS_den0;
111
112
};

113
114
115
116
/*
 * @brief Calculate if the gas has the potential of becoming
 * a star.
 *
117
118
119
120
121
 * @param starform the star formation law properties to use.
 * @param p the gas particles
 * @param xp the additional properties of the gas particles
 * @param phys_const the physical constants in internal units
 * @param cosmo the cosmological parameters and properties
122
123
 *
 * */
124
INLINE static int star_formation_potential_to_become_star(
125
126
    const struct star_formation* starform, const struct part* restrict p,
    const struct xpart* restrict xp, const struct phys_const* const phys_const,
127
128
129
130
    const struct cosmology* cosmo){

  /* Read the critical overdensity factor and the critical density of 
   * the universe to determine the critical density to form stars*/
Folkert Nobels's avatar
Folkert Nobels committed
131
  const double rho_crit = cosmo->critical_density*starform->min_over_den; 
132
  
133
134
  /* double tempp = cooling_get_temperature() */
  double tempp;
135
136

  
137
138
139
140
141
142
143
  /* Deside whether we should form stars or not,
   * first we deterime if we have the correct over density
   * if that is true we calculate if either the maximum density
   * threshold is reached or if the metallicity dependent 
   * threshold is reached, after this we calculate if the 
   * temperature is appropriate */
  if (p->rho > rho_crit ) {
144
145
146
147
    /* In this case there are actually multiple possibilities
     * because we also need to check if the physical density exceeded
     * the appropriate limit */

148
    /* Check if it exceeded the maximum density */
149
    if (p->rho > starform->den_crit_max) {
150
      /* double tempp = cooling_get_temperature() */
151
      tempp = 5e3;
152
153
154
155
156
157
158
159
      /* Check the last criteria, if the temperature is satisfied */
      if (tempp < starform->T_crit) {
        return 1;
      } else {
        return 0;
      }
    /* If the maximum density is not exceeded check if the redshift dependent
     * one is exceeded */
160
161
162
    } else {
      /* NEED TO USE A PROPER WAY OF FINDING Z */
      double Z = 0.002;
163
      double den_crit_current = starform->den_crit * pow(Z*starform->Z0_inv, starform->n_Z0);
164
      if (p->rho > den_crit_current) {
165
        /* double tempp = cooling_get_temperature() */
166
        tempp = 5e3;
167
168
169
170
171
172
        /* Check the last criteria, if the temperature is satisfied */
        if (tempp < starform->T_crit) {
          return 1;
        } else {
          return 0;
        }
173
174
175
176
      } else {
        return 0;
      }
    }
177
178
179
  } else {
    return 0;
  }
180
181
182
}

/*
183
 * @brief Calculates if the gas particle gets converted
184
 *
185
186
 */
INLINE static int star_formation_convert_to_star(
187
    const struct star_formation* starform, const struct part* restrict p,
188
189
    const struct xpart* restrict xp, const struct phys_const* const phys_const,
    const struct cosmology* cosmo) {
Folkert Nobels's avatar
Folkert Nobels committed
190

191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
  if (star_formation_potential_to_become_star(starform, p, xp, phys_const, cosmo)){
    /* Get the pressure */
    const double pressure = starform->EOS_pressure_norm 
    * pow(p->rho/starform->EOS_den0, starform->polytropic_index);
  
    /* Calculate the propability of forming a star */ 
    const double prop = starform->SF_normalization * pow(pressure,
    starform->SF_power_law) * p->time_bin;

    /* Calculate the seed */
    unsigned int seed = p->id;
    
    /* Generate a random number between 0 and 1. */
    const double randomnumber = rand_r(&seed)*starform->inv_RAND_MAX; 

    /* Calculate if we form a star */
    if (prop > randomnumber) {
      return 1;
    } else {
      return 0;
    }
Folkert Nobels's avatar
Folkert Nobels committed
212

213
214
215
216
  } else {
    return 0;
  }
}
Folkert Nobels's avatar
Folkert Nobels committed
217

218
219
220
INLINE static void star_formation_copy_properties(
    struct engine *e, struct cell *c, struct part* p,
    struct xpart* xp, const struct star_formation* starform, 
221
222
    const struct phys_const* const phys_const, const struct cosmology* cosmo,
    int with_cosmology) {
223
224
225
226
  
  struct spart *sp = cell_convert_part_to_spart(e, c, p, xp);
  sp->mass = p->mass;
  sp->mass_init = p->mass;
227
  if (with_cosmology) {
Folkert Nobels's avatar
Folkert Nobels committed
228
    sp->birth_time = cosmo->a;
229
  } else {
Folkert Nobels's avatar
Folkert Nobels committed
230
    sp->birth_time = e->time;
231
  }
Folkert Nobels's avatar
Folkert Nobels committed
232
233
234
235
236
  sp->chemistry_data = p->chemistry_data;
  //sp->tracers_data = p->tracers_data;
  //sp->birth_density = p->density;
  

Folkert Nobels's avatar
Folkert Nobels committed
237

238
}
239

240
241
242
243
244
245
246
247
248
249
250
251
252
/*
int banana(...) {

  if(star_formation_potential_to_become_ater(....) {

  //draw ranfom number
  //  return 1 or 0

  }else
  return 0;
  }
 */

253
/* 
254
 * @brief initialization of the star formation law 
255
256
257
258
259
260
 *
 * @param parameter_file The parsed parameter file
 * @param phys_const Physical constants in internal units
 * @param us The current internal system of units
 * @param starform the star formation law properties to initialize
 *
261
 * */
262
INLINE static void starformation_init_backend(
263
  struct swift_params* parameter_file, const struct phys_const* phys_const,
264
  const struct unit_system* us, struct star_formation* starform) {
Folkert Nobels's avatar
Folkert Nobels committed
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281

  /* Get the appropriate constant to calculate the 
   * star formation constant */ 
  const double KS_const = phys_const->const_kennicutt_schmidt_units;

  /* Get the Gravitational constant */
  const double G_newton = phys_const->const_newton_G;

  /* Get the surface density unit M_\odot / pc^2 */
  const double M_per_pc2 = phys_const->const_solar_mass_per_parsec2;
  
  /* Calculate inverse of RAND_MAX for the random numbers */
  starform->inv_RAND_MAX = 1.f / RAND_MAX;

  /* Quantities that have to do with the Normal Kennicutt-
   * Schmidt law will be read in this part of the code*/

282
  /* Read the critical density contrast from the parameter file*/
Folkert Nobels's avatar
Folkert Nobels committed
283
  starform->min_over_den = parser_get_param_double(parameter_file, 
284
  "SchayeSF:thresh_MinOverDens");
285
286
287

  /* Read the critical temperature from the parameter file */
  starform->T_crit = parser_get_param_double(parameter_file,
288
  "SchayeSF:thresh_temp");
289
290

  /* Read the gas fraction from the file */
291
  starform->fgas = parser_get_param_double(parameter_file,
292
293
  "SchayeSF:fg");

Folkert Nobels's avatar
Folkert Nobels committed
294
  /* Read the normalization of the KS law in KS law units */
295
296
  const double normalization_MSUNpYRpKPC2 = parser_get_param_double(
  parameter_file, "SchayeSF:SchmidtLawCoeff_MSUNpYRpKPC2");
297
298

  /* Read the Kennicutt-Schmidt power law exponent */
299
300
  starform->KS_power_law = parser_get_param_double(
  parameter_file, "SchayeSF:SchmidtLawExponent");
301

Folkert Nobels's avatar
Folkert Nobels committed
302
303
304
305
306
307
308
309
  /* Calculate the power law of the star formation */
  starform->SF_power_law = (starform->KS_power_law - 1.f)/2.f;
  
  /* Give the Kennicutt-Schmidt law the same units as internal units */
  starform->KS_normalization = normalization_MSUNpYRpKPC2 * KS_const;
  
  /* Calculate the starformation prefactor with most terms */
  starform->SF_normalization = starform->KS_normalization * pow(M_per_pc2, -starform->KS_power_law) * 
310
  pow( hydro_gamma * starform->fgas / G_newton, starform->SF_power_law);
Folkert Nobels's avatar
Folkert Nobels committed
311

312
  /* Read the high density Kennicutt-Schmidt power law exponent */
313
314
  starform->KS_high_den_power_law = parser_get_param_double(
  parameter_file, "SchayeSF:SchmidtLawHighDensExponent");
315
316
  
  /* Read the high density criteria for the KS law in number density per cm^3 */
317
318
  const double KS_high_den_thresh_HpCM3 = parser_get_param_double(
  parameter_file, "SchayeSF:SchmidtLawHighDens_thresh_HpCM3");
319
320
321
322

  /* Transform the KS high density criteria to simulation units */
  starform->KS_high_den_thresh = KS_high_den_thresh_HpCM3 * UNIT_CONV_NUMBER_DENSITY;

Folkert Nobels's avatar
Folkert Nobels committed
323
324
  /* Calculate the SF high density power law */
  starform->SF_high_den_power_law = (starform->KS_high_den_power_law - 1.f)/2.f;
325

Folkert Nobels's avatar
Folkert Nobels committed
326
327
  /* Calculate the KS high density normalization */
  starform->KS_high_den_normalization = starform->KS_normalization * pow( M_per_pc2,
328
  starform->KS_high_den_power_law - starform->KS_power_law) * pow( hydro_gamma * 
Folkert Nobels's avatar
Folkert Nobels committed
329
330
  starform->fgas / G_newton * 1337.f, (starform->KS_power_law 
  - starform->KS_high_den_power_law)/2.f);
331
332
333
334
335
336

  /* Read the critical temperature from the parameter file */
  starform->T_crit = parser_get_param_double(parameter_file,
  "SchayeSF:T_crit");

  /* Read the gas fraction from the file */
Folkert Nobels's avatar
Folkert Nobels committed
337
  starform->fgas = parser_get_param_double(parameter_file,
338
339
  "SchayeSF:fg");

340
341
  /* Calculate inverse of RAND_MAX */
  starform->inv_RAND_MAX = 1.f / RAND_MAX;
342

Folkert Nobels's avatar
Folkert Nobels committed
343
344
  /* Calculate the SF high density normalization */
  starform->SF_high_den_normalization = starform->KS_high_den_normalization 
345
  * pow(M_per_pc2, -starform->KS_high_den_power_law) * pow( hydro_gamma  
Folkert Nobels's avatar
Folkert Nobels committed
346
  * starform->fgas / G_newton, starform->SF_high_den_power_law);
347

348
349
  /* Read what kind of critical density we need to use
   * Schaye (2004) is metallicity dependent critical SF density*/
350
351
  starform->schaye04 = parser_get_param_double(
  parameter_file, "SchayeSF:Schaye2004");
352

Folkert Nobels's avatar
Folkert Nobels committed
353
  if (!starform->schaye04) {
354
355
    /* In the case that we do not use the Schaye (2004) critical
     * density to form stars but a constant value */
356
357
358
    starform->den_crit = parser_get_param_double(
    parameter_file, "SchayeSF:thresh_norm_HpCM3");
    starform->Z0 = 0.002;
359
360
361
362
363
364
    starform->n_Z0 = 0.0;
  } else {
    /* Use the Schaye (2004) metallicity dependent critical density
     * to form stars. */
    /* Read the normalization of the metallicity dependent critical 
     * density*/
365
366
    starform->den_crit = parser_get_param_double( 
    parameter_file, "SchayeSF:thresh_norm_HpCM3") *
Folkert Nobels's avatar
Folkert Nobels committed
367
    UNIT_CONV_NUMBER_DENSITY;
368
369

    /* Read the scale metallicity Z0 */
370
371
    starform->Z0 = parser_get_param_double(
    parameter_file, "SchayeSF:MetDep_Z0");
372
373

    /* Read the power law of the critical density scaling */
374
375
    starform->n_Z0 = parser_get_param_double(
    parameter_file, "SchayeSF:MetDep_SFthresh_Slope");
376
377

    /* Read the maximum allowed density for star formation */
378
379
    starform->den_crit_max = parser_get_param_double(
    parameter_file, "SchayeSF:thresh_max_norm_HpCM3");
380

381
  }
382
  /* Conversion of number density from cgs */
383
384
  static const float dimension_numb_den[5] = {0, -3, 0, 0, 0};
  const double conversion_numb_density = 1.f/
385
  units_general_cgs_conversion_factor(us, dimension_numb_den);
386

387
388
389
  /* Claculate 1 over the metallicity for speed up */
  starform->Z0_inv = 1/starform->Z0;

390
391
392
  /* Calculate the prefactor that is always common */
  /* !!!DONT FORGET TO DO THE CORRECT UNIT CONVERSION!!!*/
  starform->den_crit_star = starform->den_crit / pow(starform->Z0,
393
  starform->n_Z0) * conversion_numb_density;
394
395
396
397
398
399
400
401
402

  /* Load the equation of state for this model */
  starform->polytropic_index = parser_get_param_double(
  parameter_file, "SchayeSF:EOS_Jeans_GammaEffective");
  starform->EOS_pressure_norm = parser_get_param_double(
  parameter_file, "SchayeSF:EOS_Jeans_PressureNorm");
  starform->EOS_den0 = parser_get_param_double(
  parameter_file, "SchayeSF:EOS_JEANS_den0");
  
403
404
405
406
}

/* @brief Prints the used parameters of the star formation law 
 *
407
 * @param starform the star formation law properties.
408
 * */
409
INLINE static void starformation_print_backend(
410
411
    const struct star_formation* starform){ 

412
413
  message("Star formation law is Schaye and Dalla Vecchia (2008)"
  " with properties, normalization = %e, slope of the Kennicutt"
414
  "-Schmidt law = %e, gas fraction = %e, critical "
415
  "density = %e and critical temperature = %e", starform->KS_normalization, 
416
  starform->KS_power_law, starform->fgas, starform->den_crit,
417
  starform->T_crit);
Folkert Nobels's avatar
Folkert Nobels committed
418
419
420
421
422
423
424
425
426
  if (!starform->schaye04) {
    message("Density threshold to form stars is constant and given"
    "by a density of %e", starform->den_crit_star);
  } else {
    message("Density threshold to form stars is given by Schaye "
    "(2004), the normalization of the star formation law is given by"
    " %e, with metallicity slope of %e, and metallicity normalization"
    "of %e", starform->den_crit_star, starform->n_Z0, starform->Z0);
  }
427
428
429
430
431

}


#endif /* SWIFT_SCHAYE_STARFORMATION_H */