-
Notifications
You must be signed in to change notification settings - Fork 20
/
Copy pathchem_evol.py
8256 lines (6429 loc) · 316 KB
/
chem_evol.py
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
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
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
196
197
198
199
200
201
202
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
269
270
271
272
273
274
275
276
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
336
337
338
339
340
341
342
343
344
345
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
391
392
393
394
395
396
397
398
399
400
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
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# coding=utf-8
from __future__ import (division, print_function, absolute_import,
unicode_literals)
'''
Chemical Evolution - chem_evol.py
Functionality
=============
This is the superclass inherited by the SYGMA and the OMEGA modules. It provides
common functions for initialization and for the evolution of one single timestep.
Made by
=======
MAY2015: B. Cote
The core of this superclass is a reorganization of the functions previously found in
earlier versions of SYGMA:
v0.1 NOV2013: C. Fryer, C. Ritter
v0.2 JAN2014: C. Ritter
v0.3 APR2014: C. Ritter, J. F. Navarro, F. Herwig, C. Fryer, E. Starkenburg,
M. Pignatari, S. Jones, K. Venn, P. A. Denissenkov & the
NuGrid collaboration
v0.4 FEB2015: C. Ritter, B. Cote
v0.5 OCT2016: B. Cote, C. Ritter, A. Paul
Stop keeking track of version from now on.
MARCH2018: B. Cote
- Switched to Python 3
- Capability to include radioactive isotopes
JULY2018: B. Cote & R. Sarmento
- Re-wrote (improved) yield and lifetime treatment (B. Cote)
- PopIII IMF and yields update (R. Sarmento)
JAN2019: B. Cote
- Re-included radioactive isotopes with the new (improved) yield treatment
FEB2019: A. Yagüe, B. Cote
- Optimized to code to run faster (integration method)
Usage
=====
See sygma.py and omega.py
'''
# Standard packages
import numpy as np
import time as t_module
import copy
import os
import sys
import re
from pylab import polyfit
from scipy.integrate import quad
from scipy.integrate import dblquad
# Define where is the working directory
nupy_path = os.path.dirname(os.path.realpath(__file__))
# Import NuPyCEE codes
import NuPyCEE.read_yields as ry
class chem_evol(object):
'''
Input parameters (chem_evol.py)
================
special_timesteps : integer
Number of special timesteps. This option (already activated by default)
is activated when special_timesteps > 0. It uses a logarithm timestep
scheme which increases the duration of timesteps throughout the simulation.
Default value : 30
dt : float
Duration of the first timestep [yr] if special_timesteps is activated.
Duration of each timestep if special_timesteps is desactivated.
Default value : 1.0e6
tend : float
Total duration of the simulation [yr].
Default value : 13.0e9
dt_split_info : numpy array
Information regarding the creation of a timestep array with varying step size.
Array format : dt_split_info[number of conditions][0->dt,1->upper time limit]
Exemple : dt_split_info = [[1e6,40e6],[1e8,13e9]] means the timesteps will be
of 1 Myr until the time reaches 40 Myr, after which the timesteps
will be of 100 Myr until the time reaches 13 Gyr. The number of
"split" is unlimited, but the array must be in chronological order.
Default value : [] --> Not taken into account
imf_bdys : list
Upper and lower mass limits of the initial mass function (IMF) [Mo].
Default value : [0.1,100]
imf_yields_range : list
Initial mass of stars that contribute to stellar ejecta [Mo].
Default value : [1,30]
imf_type : string
Choices : 'salpeter', 'chabrier', 'kroupa', 'alphaimf', 'lognormal'
'alphaimf' creates a custom IMF with a single power-law covering imf_bdys.
'lognormal' creates an IMF of the form Exp[1/(2 1^2) Log[x/charMass]^2
Default value : 'kroupa'
alphaimf : float
Aplha index of the custom IMF, dN/dM = Constant * M^-alphaimf
Default value : 2.35
imf_bdys_pop3 : list
Upper and lower mass limits of the IMF of PopIII stars [Mo].
Default value : [0.1,100]
imf_yields_range_pop3 : list
Initial mass of stars that contribute to PopIII stellar ejecta [Mo].
PopIII stars ejecta taken from Heger et al. (2010)
Default value : [10,30]
imf_pop3_char_mass : float
The characteristic mass in a log normal IMF distribution.
Default value : 40.0
high_mass_extrapolation : string
Extrapolation technique used to extrapolate yields for stars more
massive than the most massive model (MMM) present in the yields table.
Choices:
"copy" --> This will apply the yields of the most massive model
to all more massive stars.
"scale" --> This will scale the yields of the most massive model
using the relation between the total ejected mass and
the initial stellar mass. The later relation is taken
from the interpolation of the two most massive models.
"extrapolate" --> This will extrapolate the yields of the most massive
model using the interpolation coefficients taken from
the interpolation of the two most massive models.
Default value : "copy"
iniZ : float
Initial metallicity of the gas in mass fraction (e.g. Solar Z = 0.02).
Choices : 0.0, 0.0001, 0.001, 0.006, 0.01, 0.02
(-1.0 to use non-default yield tables)
Default value : 0.0
Z_trans : float
Variable used when interpolating stellar yields as a function of Z.
Transition Z below which PopIII yields are used, and above which default
yields are used.
Default value : -1 (not active)
mgal : float
Initial mass of gas in the simulation [Mo].
Default value : 1.6e11
sn1a_on : boolean
True or False to include or exclude the contribution of SNe Ia.
Default value : True
sn1a_rate : string
SN Ia delay-time distribution function used to calculate the SN Ia rate.
Choices :
'power_law' - custom power law, set parameter with beta_pow (similar to Maoz & Mannucci 2012)
'gauss' - gaussian DTD, set parameter with gauss_dtd
'exp' - exponential DTD, set parameter with exp_dtd
'maoz' - specific power law from Maoz & Mannucci (2012)
Default value : 'power_law'
sn1a_energy : float
Energy ejected by single SNIa event. Units in erg.
Default value : 1e51
ns_merger_on : boolean
True or False to include or exclude the contribution of neutron star mergers.
Note : If t_nsm_coal or nsm_dtd_power is not used (see below), the delay time
distribution of neutron star mergers is given by the standard population synthesis
models of Dominik et al. (2012), using Z = 0.002 and Z = 0.02. In this case, the
total number of neutron star mergers can be tuned using f_binary and f_merger
(see below).
Default value : False
f_binary : float
Binary fraction for massive stars used to determine the total number of neutron
star mergers in a simple stellar population.
Default value : 1.0
f_merger : float
Fraction of massive star binary systems that lead to neutron star mergers in a
simple stellar population.
Default value : 0.0008
beta_pow : float
Slope of the power law for custom SN Ia rate, R = Constant * t^-beta_pow.
Default value : -1.0
gauss_dtd : list
Contains parameter for the gaussian DTD: first the characteristic time [yrs] (gaussian center)
and then the width of the distribution [yrs].
Default value : [3.3e9,6.6e8]
exp_dtd : float
Characteristic delay time [yrs] for the e-folding DTD.
nb_1a_per_m : float
Number of SNe Ia per stellar mass formed in a simple stellar population.
Default value : 1.0e-03
direct_norm_1a : float
Normalization coefficient for SNIa rate integral.
Default: deactived but replaces the usage of teh nb_1a_per_m when its value is larger than zero.
transitionmass : float
Initial mass which marks the transition from AGB to massive stars [Mo].
Default value : 8.0
exclude_masses : list
Contains initial masses in yield tables to be excluded from the simulation;
Default value : []
table : string
Path pointing toward the stellar yield tables for massive and AGB stars.
Default value : 'yield_tables/agb_and_massive_stars_nugrid_MESAonly_fryer12delay.txt' (NuGrid)
sn1a_table : string
Path pointing toward the stellar yield table for SNe Ia.
Default value : 'yield_tables/sn1a_i99_W7.txt' (Iwamoto et al. 1999)
nsmerger_table : string
Path pointing toward the r-process yield tables for neutron star mergers
Default value : 'yield_tables/r_process_rosswog_2014.txt' (Rosswog et al. 2013)
iniabu_table : string
Path pointing toward the table of initial abuncances in mass fraction.
Default value : 'yield_tables/iniabu/iniab2.0E-02GN93.ppn'
yield_tables_dir : string
Path to a custom directory that includes yields.
!! It needs to point to the directory where the yields directory is !!
This will bypass the default yields directory.
Default value : '' --> Deactivated
yield_interp : string
if 'None' : no yield interpolation, no interpolation of total ejecta
if 'lin' - Simple linear yield interpolation.
if 'wiersma' - Interpolation method which makes use of net yields
as used e.g. in Wiersma+ (2009); Does not require net yields.
if netyields_on is true than makes use of given net yields
else calculates net yields from given X0 in yield table.
Default : 'lin'
netyields_on : boolean
if true assumes that yields (input from table parameter)
are net yields.
Default : false.
total_ejecta_interp : boolean
if true then interpolates total ejecta given in yield tables
over initial mass range.
Default : True
stellar_param_on : boolean
if true reads in additional stellar parameter given in table stellar_param_table.
Default : true in sygma and false in omega
stellar_param_table: string
Path pointoing toward the table hosting the evolution of stellar parameter
derived from stellar evolution calculations.
Default table : 'yield_tables/isotope_yield_table_MESA_only_param.txt'
iolevel : int
Specifies the amount of output for testing purposes (up to 3).
Default value : 0
poly_fit_dtd : list
Array of polynomial coefficients of a customized delay-time distribution
function for SNe Ia. The polynome can be of any order.
Example : [0.2, 0.3, 0.1] for rate_snIa(t) = 0.2*t**2 + 0.3*t + 0.1
Note : Must be used with the poly_fit_range parameter (see below)
Default value : np.array([]) --> Deactivated
poly_fit_range : list --> [t_min,t_max]
Time range where the customized delay-time distribution function for
SNe Ia will be applied for a simple stellar population.
Default value : np.array([]) --> Deactivated
scale_cor : 2D list
Determine the fraction of yields ejected for any given stellar mass bin.
Example : [ [1.0,8], [0.5,100] ] means that stars with initial mass between
0 and 8 Msu will eject 100% of their yields, and stars with initial mass
between 8 and 100 will eject 50% of their yields. There is no limit for
the number of [%,M_upper_limit] arrays used.
Default value : np.array([]) --> Deactivated
t_nsm_coal : float
When greater than zero, t_nsm_coal sets the delay time (since star formation)
after which all neutron star mergers occur in a simple stellar population.
Default value : -1 --> Deactivated
nsm_dtd_power : 3-index array --> [t_min, t_max, slope_of_the_power_law]
When used, nsm_dtd_power defines a delay time distribution for neutron
star mergers in the form of a power law, for a simple stellar population.
Exemple: [1.e7, 1.e10, -1.] --> t^-1 from 10 Myr to 10 Gyr
Default value : [] --> Deactivated
nb_nsm_per_m : float
Number of neutron star mergers per stellar mass formed in a simple
stellar population.
Note : This parameter is only considered when t_nsm_coal or nsm_dtd_power
is used to define the delay time of neutron star mergers.
Default value : -1 --> Deactivated
m_ej_nsm : float
Mass ejected per neutron star merger event.
Default value : 2.5e-02
yield_modifier : list of arrays --> [[iso, M, Z, type, modifier],[...]]
When used, modifies all isotopes yields for the given M and Z by
multiplying by a given factor (type="multiply") or replacing the
yield by a new value (type="replace"). Modifier will be either the
factor or value depending on type.
Default value : [] --> Deactivated
Delayed extra source
Adding source that requires delay-time distribution (DTD) functions
-------------------------------------------------------------------
delayed_extra_dtd : multi-D Numpy array --> [nb_sources][nb_Z]
nb_sources is the number of different input astrophysical site (e.g.,
SNe Ia, neutron star mergers).
nb_Z is the number of available metallicities.
delayed_extra_dtd[i][j] is a 2D array in the form of
[ number_of_times ][ 0-time, 1-rate ].
Defalut value : np.array([]), deactivated
delayed_extra_dtd_norm : multi-D Numpy array --> [nb_sources]
Total number of delayed sources occurring per Msun formed,
for each source and each metallicity.
Defalut value : np.array([]), deactivated
delayed_extra_yields : Numpy array of strings
Path to the yields table for each source.
Defalut value : np.array([]), deactivated
delayed extra_yields_norm : multi-D Numpy array --> [nb_sources][nb_Z]
Fraction of the yield table (float) that will be ejected per event,
for each source and each metallicity. This will be the mass ejected
per event if the yields are in mass fraction (normalized to 1).
Defalut value : np.array([]), deactivated
delayed_extra_stochastic : Numpy array of Boolean --> [nb_sources]
Determine whether the DTD provided as an input needs to be
stochastically sampled using a Monte Carlo technique.
Defalut value : np.array([]), deactivated
Run example
===========
See sygma.py and omega.py
'''
##############################################
## Constructor ##
##############################################
def __init__(self, imf_type='kroupa', alphaimf=2.35, imf_bdys=[0.1,100],\
sn1a_rate='power_law', iniZ=0.02, dt=1e6, special_timesteps=60,\
nsmerger_bdys=[8, 100], tend=13e9, mgal=1.0e11, transitionmass=8, iolevel=0,\
ini_alpha=True, is_sygma=False,\
table='yield_tables/agb_and_massive_stars_nugrid_MESAonly_fryer12delay.txt',\
f_network='isotopes_modified.prn', f_format=1,\
table_radio='', decay_file='', sn1a_table_radio='',\
nsmerger_table_radio='',\
hardsetZ=-1, sn1a_on=True, sn1a_table='yield_tables/sn1a_i99_W7.txt',\
sn1a_energy=1e51, ns_merger_on=False,\
f_binary=1.0, f_merger=0.0008, t_merger_max=1.3e10,\
m_ej_nsm = 2.5e-02, nb_nsm_per_m=-1.0,\
t_nsm_coal=-1.0, nsm_dtd_power=[],\
nsmerger_table = 'yield_tables/r_process_arnould_2007.txt',\
iniabu_table='', extra_source_on=False,\
extra_source_table=['yield_tables/extra_source.txt'],\
f_extra_source=[1.0], pre_calculate_SSPs=False,\
extra_source_mass_range=[[8,30]],\
extra_source_exclude_Z=[[]], radio_refinement=100,\
pop3_table='yield_tables/popIII_heger10.txt',\
imf_bdys_pop3=[0.1,100], imf_yields_range_pop3=[10,30],\
imf_pop3_char_mass=40.0,\
use_net_yields_stable=False, use_net_yields_radio=False, \
high_mass_extrapolation='copy',\
use_external_integration=False,\
beta_pow=-1.0,gauss_dtd=[3.3e9,6.6e8],\
exp_dtd=2e9,nb_1a_per_m=1.0e-3,direct_norm_1a=-1,Z_trans=1.0e-20, \
f_arfo=1, imf_yields_range=[1,30],exclude_masses=[],\
netyields_on=False,wiersmamod=False,yield_interp='lin',\
print_off=False, yield_tables_dir='',\
total_ejecta_interp=True, tau_ferrini=False,\
input_yields=False,t_merge=-1.0,stellar_param_on=False,\
stellar_param_table='yield_tables/stellar_feedback_nugrid_MESAonly.txt',\
popIII_info_fast=True, out_follows_E_rate=False,\
t_dtd_poly_split=-1.0, delayed_extra_log=False,\
delayed_extra_yields_log_int=False,\
delayed_extra_log_radio=False, delayed_extra_yields_log_int_radio=False,\
pritchet_1a_dtd=[], ism_ini=[], ism_ini_radio=[],\
nsmerger_dtd_array=[],\
ytables_in=[], zm_lifetime_grid_nugrid_in=[],\
isotopes_in=[], ytables_pop3_in=[],\
zm_lifetime_grid_pop3_in=[], ytables_1a_in=[],\
ytables_nsmerger_in=[], dt_in_SSPs=[],\
dt_in=[],dt_split_info=[],\
ej_massive=[], ej_agb=[],\
ej_sn1a=[], ej_massive_coef=[],\
ej_agb_coef=[], ej_sn1a_coef=[],\
dt_ssp=[], poly_fit_dtd_5th=[],\
poly_fit_range=[], SSPs_in=[],\
delayed_extra_dtd=[], delayed_extra_dtd_norm=[],\
delayed_extra_yields=[], delayed_extra_yields_norm=[],\
delayed_extra_yields_radio=[],\
delayed_extra_yields_norm_radio=[],\
delayed_extra_stochastic=[],\
ytables_radio_in=[], radio_iso_in=[],\
ytables_1a_radio_in=[], ytables_nsmerger_radio_in=[],\
test_clayton=[], inter_Z_points=[],\
nb_inter_Z_points=[], y_coef_M=[],\
y_coef_M_ej=[], y_coef_Z_aM=[],\
y_coef_Z_bM=[], y_coef_Z_bM_ej=[],\
tau_coef_M=[], tau_coef_M_inv=[],\
tau_coef_Z_aM=[], tau_coef_Z_bM=[],\
tau_coef_Z_aM_inv=[], tau_coef_Z_bM_inv=[],\
y_coef_M_pop3=[], y_coef_M_ej_pop3=[],\
tau_coef_M_pop3=[], tau_coef_M_pop3_inv=[],\
inter_lifetime_points_pop3=[],\
inter_lifetime_points_pop3_tree=[],\
nb_inter_lifetime_points_pop3=[],\
inter_lifetime_points=[], inter_lifetime_points_tree=[],\
nb_inter_lifetime_points=[], nb_inter_M_points_pop3=[],\
inter_M_points_pop3_tree=[], nb_inter_M_points=[],\
inter_M_points=[], y_coef_Z_aM_ej=[],
yield_modifier=[]):
# Initialize the history class which keeps the simulation in memory
self.history = History()
self.const = Const()
# Define elements for ordering
self.__define_elements()
# If we need to assume the current baryonic ratio ...
if mgal < 0.0:
# Use a temporary mgal value for chem_evol __init__ function
mgal = 1.0e06
self.bar_ratio = True
# If we use the input mgal parameter ...
else:
self.bar_ratio = False
# Attribute the input parameters to the current object
self.history.mgal = mgal
self.history.tend = tend
self.history.dt = dt
self.history.sn1a_rate = sn1a_rate
self.history.imf_bdys = imf_bdys
self.history.transitionmass = transitionmass
self.history.nsmerger_bdys = nsmerger_bdys
self.history.f_binary = f_binary
self.history.f_merger = f_merger
self.mgal = mgal
self.transitionmass = transitionmass
self.iniZ = iniZ
self.imf_bdys=imf_bdys
self.nsmerger_bdys=nsmerger_bdys
self.popIII_info_fast = popIII_info_fast
self.imf_bdys_pop3=imf_bdys_pop3
self.imf_yields_range_pop3=imf_yields_range_pop3
self.imf_pop3_char_mass=imf_pop3_char_mass
self.high_mass_extrapolation = high_mass_extrapolation
self.extra_source_on = extra_source_on
self.f_extra_source= f_extra_source
self.extra_source_mass_range=extra_source_mass_range
self.extra_source_exclude_Z=extra_source_exclude_Z
self.pre_calculate_SSPs = pre_calculate_SSPs
self.SSPs_in = SSPs_in
self.table = table
self.iniabu_table = iniabu_table
self.sn1a_table = sn1a_table
self.nsmerger_table = nsmerger_table
self.extra_source_table = extra_source_table
self.pop3_table = pop3_table
self.hardsetZ = hardsetZ
self.imf_type = imf_type
self.alphaimf = alphaimf
self.sn1a_on = sn1a_on
self.sn1a_energy=sn1a_energy
self.ns_merger_on = ns_merger_on
self.nsmerger_dtd_array = nsmerger_dtd_array
self.len_nsmerger_dtd_array = len(nsmerger_dtd_array)
self.f_binary = f_binary
self.f_merger = f_merger
self.t_merger_max = t_merger_max
self.m_ej_nsm = m_ej_nsm
self.nb_nsm_per_m = nb_nsm_per_m
self.t_nsm_coal = t_nsm_coal
self.nsm_dtd_power = nsm_dtd_power
self.special_timesteps = special_timesteps
self.iolevel = iolevel
self.nb_1a_per_m = nb_1a_per_m
self.direct_norm_1a=direct_norm_1a
self.Z_trans = Z_trans
if sn1a_rate == 'maoz':
self.beta_pow = -1.0
else:
self.beta_pow = beta_pow
self.gauss_dtd = gauss_dtd
self.exp_dtd=exp_dtd
self.normalized = False # To avoid normalizing SN Ia rate more than once
self.nsm_normalized = False # To avoid normalizing NS merger rate more than once
self.f_arfo = f_arfo
self.imf_yields_range = imf_yields_range
self.exclude_masses=exclude_masses
self.netyields_on=netyields_on
self.wiersmamod=wiersmamod
self.yield_interp=yield_interp
self.out_follows_E_rate = out_follows_E_rate
self.total_ejecta_interp=total_ejecta_interp
self.tau_ferrini = tau_ferrini
self.t_merge = t_merge
self.ism_ini = ism_ini
self.ism_ini_radio = ism_ini_radio
self.dt_in = dt_in
self.dt_in_SSPs = dt_in_SSPs
self.dt_split_info = dt_split_info
self.t_dtd_poly_split = t_dtd_poly_split
self.poly_fit_dtd_5th = poly_fit_dtd_5th
self.poly_fit_range = poly_fit_range
self.stellar_param_table = stellar_param_table
self.stellar_param_on = stellar_param_on
self.delayed_extra_log = delayed_extra_log
self.delayed_extra_dtd = delayed_extra_dtd
self.delayed_extra_dtd_norm = delayed_extra_dtd_norm
self.delayed_extra_yields = delayed_extra_yields
self.delayed_extra_yields_norm = delayed_extra_yields_norm
self.delayed_extra_yields_log_int = delayed_extra_yields_log_int
self.delayed_extra_stochastic = delayed_extra_stochastic
self.nb_delayed_extra = len(self.delayed_extra_dtd)
self.pritchet_1a_dtd = pritchet_1a_dtd
self.len_pritchet_1a_dtd = len(pritchet_1a_dtd)
self.use_external_integration = use_external_integration
self.yield_tables_dir = yield_tables_dir
self.print_off = print_off
self.use_net_yields_stable = use_net_yields_stable
self.input_yields = input_yields
self.yield_modifier = yield_modifier
self.is_sygma = is_sygma
self.iolevel = iolevel
# Attributes associated with radioactive species
self.table_radio = table_radio
self.sn1a_table_radio = sn1a_table_radio
self.nsmerger_table_radio = nsmerger_table_radio
self.decay_file = decay_file
self.len_decay_file = len(decay_file)
self.delayed_extra_log_radio = delayed_extra_log_radio
self.delayed_extra_yields_radio = delayed_extra_yields_radio
self.delayed_extra_yields_norm_radio = delayed_extra_yields_norm_radio
self.delayed_extra_yields_log_int_radio = delayed_extra_yields_log_int_radio
self.nb_delayed_extra_radio = len(self.delayed_extra_yields_radio)
self.ytables_radio_in = ytables_radio_in
self.radio_iso_in = radio_iso_in
self.ytables_1a_radio_in = ytables_1a_radio_in
self.ytables_nsmerger_radio_in = ytables_nsmerger_radio_in
self.radio_massive_agb_on = False
self.radio_sn1a_on = False
self.radio_nsmerger_on = False
self.radio_refinement = radio_refinement
self.test_clayton = test_clayton
self.use_net_yields_radio = use_net_yields_radio
# Number of coefficients for the interpolation routine
self.nb_c_needed = 3
# Define the use of the decay_module and the decay_file
is_radio = not (len(self.table_radio) == 0 \
and len(self.sn1a_table_radio) == 0 \
and len(self.nsmerger_table_radio) == 0 \
and self.nb_delayed_extra_radio == 0)
if is_radio and self.len_decay_file == 0:
self.use_decay_module = True
else:
self.use_decay_module = False
# Initialize decay module
if self.use_decay_module:
import NuPyCEE.decay_module as decay_module
self.decay_module = decay_module
self.f_network = f_network
self.f_format = f_format
self.__initialize_decay_module()
# Normalization of the delayed extra sources
if self.nb_delayed_extra > 0:
self.__normalize_delayed_extra()
# Normalization constants for the Kroupa IMF
if self.imf_type == 'kroupa':
self.p0 = 1.0
self.p1 = 0.08**(-0.3 + 1.3)
self.p2 = 0.5**(-1.3 + 2.3)
self.p3 = 1**(-2.3 +2.3)
self.p1_p2 = self.p1 * self.p2
# Define the broken power-law of Ferrini IMF approximation
self.norm_fer = [3.1,1.929,1.398,0.9113,0.538,0.3641,0.2972,\
0.2814,0.2827,0.298,0.305,0.3269,0.3423,0.3634]
self.alpha_fer = [0.6,0.35,0.15,-0.15,-0.6,-1.05,-1.4,-1.6,-1.7,\
-1.83,-1.85,-1.9,-1.92,-1.94]
self.m_up_fer = [0.15,0.2,0.24,0.31,0.42,0.56,0.76,1.05,1.5,\
3.16,4.0,10.0,20.0,120]
for i_fer in range(0,len(self.norm_fer)):
self.alpha_fer[i_fer] = self.alpha_fer[i_fer] + 1
self.norm_fer[i_fer] = self.norm_fer[i_fer]/(self.alpha_fer[i_fer])
# Normalize the IMF to 1 MSun
self.A_imf = 1.0 / self._imf(self.imf_bdys[0], self.imf_bdys[1], 2)
self.A_imf_pop3 = 1.0 / self._imf(self.imf_bdys_pop3[0], self.imf_bdys_pop3[1], 2)
# Parameter that determines if not enough gas is available for star formation
self.not_enough_gas_count = 0
self.not_enough_gas = False
# Initialisation of the timesteps
if len(self.dt_split_info) > 0: # and len(self.ej_massive) == 0:
timesteps = self.__build_split_dt()
else:
timesteps = self.__get_timesteps()
self.history.timesteps = timesteps
self.nb_timesteps = len(timesteps)
if self.pre_calculate_SSPs:
self.t_ce = np.zeros(self.nb_timesteps)
self.t_ce[0] = self.history.timesteps[0]
for i_init in range(1,self.nb_timesteps):
self.t_ce[i_init] = self.t_ce[i_init-1] + \
self.history.timesteps[i_init]
# Define the decay properties
if self.use_decay_module:
self.__read_isotopes_and_define_decay_info()
elif self.len_decay_file > 0:
self.__define_decay_info()
# If the yield tables have already been read previously ...
if self.input_yields:
# Assign the input yields and lifetimes
self.history.isotopes = isotopes_in
self.nb_isotopes = len(self.history.isotopes)
self.ytables = ytables_in
self.ytables_1a = ytables_1a_in
self.ytables_nsmerger = ytables_nsmerger_in
self.extra_source_on = False
self.ytables_extra = 0
self.inter_Z_points = inter_Z_points
self.nb_inter_Z_points = nb_inter_Z_points
self.y_coef_M = y_coef_M
self.y_coef_M_ej = y_coef_M_ej
self.y_coef_Z_aM = y_coef_Z_aM
self.y_coef_Z_bM = y_coef_Z_bM
self.y_coef_Z_bM_ej = y_coef_Z_bM_ej
self.tau_coef_M = tau_coef_M
self.tau_coef_M_inv = tau_coef_M_inv
self.tau_coef_Z_aM = tau_coef_Z_aM
self.tau_coef_Z_bM = tau_coef_Z_bM
self.tau_coef_Z_aM_inv = tau_coef_Z_aM_inv
self.tau_coef_Z_bM_inv = tau_coef_Z_bM_inv
self.y_coef_M_pop3 = y_coef_M_pop3
self.y_coef_M_ej_pop3 = y_coef_M_ej_pop3
self.tau_coef_M_pop3 = tau_coef_M_pop3
self.tau_coef_M_pop3_inv = tau_coef_M_pop3_inv
self.inter_lifetime_points_pop3 = inter_lifetime_points_pop3
self.inter_lifetime_points_pop3_tree = inter_lifetime_points_pop3_tree
self.nb_inter_lifetime_points_pop3 = nb_inter_lifetime_points_pop3
self.inter_lifetime_points = inter_lifetime_points
self.inter_lifetime_points_tree = inter_lifetime_points_tree
self.nb_inter_lifetime_points = nb_inter_lifetime_points
self.nb_inter_M_points_pop3 = nb_inter_M_points_pop3
self.inter_M_points_pop3_tree = inter_M_points_pop3_tree
self.nb_inter_M_points = nb_inter_M_points
self.inter_M_points = inter_M_points
self.y_coef_Z_aM_ej = y_coef_Z_aM_ej
# Assign the input yields for radioactive isotopes
if self.len_decay_file > 0 or self.use_decay_module:
self.ytables_1a_radio = ytables_1a_radio_in
self.ytables_nsmerger_radio = ytables_nsmerger_radio_in
self.y_coef_M_radio = y_coef_M_radio
self.y_coef_Z_aM_radio = y_coef_Z_aM_radio
self.y_coef_Z_bM_radio = y_coef_Z_bM_radio
# If the yield tables need to be read from the files ...
else:
# Read of the yield tables
self.__read_tables()
# Read radioactive tables
if self.len_decay_file > 0 or self.use_decay_module:
self.__read_radio_tables()
# Modify the yields (ttrueman edit)
if len(self.yield_modifier) > 0:
iso = [i[0] for i in yield_modifier]
M = [i[1] for i in yield_modifier]
Z = [i[2] for i in yield_modifier]
modifier = [i[3] for i in yield_modifier]
val = [i[4] for i in yield_modifier]
for j,specie in enumerate(iso):
if Z[j] not in self.Z_table or M[j] not in self.M_table:
print('Z = %s or M_sun = %s is not in yield table'%(Z[j],M[j]))
print('No modifications will be performed on %s'%iso[j],"\n")
elif specie in self.history.isotopes:
if modifier[j] == "replace":
self.ytables.set(M=M[j],Z=Z[j],specie=iso[j],
value=val[j])
if modifier[j] == "multiply":
original = self.ytables.get(M=M[j],Z=Z[j],
quantity=iso[j])
self.ytables.set(M=M[j],Z=Z[j],specie=iso[j],
value=original*val[j])
elif self.len_decay_file > 0:
if specie in self.radio_iso:
if modifier[j] == "replace":
self.ytables_radio.set(M=M[j],Z=Z[j],specie=iso[j],
value=val[j])
if modifier[j] == "multiply":
original = self.ytables_radio.get(M=M[j],Z=Z[j],
quantity=iso[j])
self.ytables_radio.set(M=M[j],Z=Z[j],specie=iso[j],
value=original*val[j])
else:
print("ERROR 404: %s not found in list of isotopes"%specie,
"\n")
else:
print("ERROR 404: %s not found in list of isotopes"%specie,
"\n")
# Declare the interpolation coefficient arrays
self.__declare_interpolation_arrays()
# Interpolate the yields tables
self.__interpolate_pop3_yields()
self.__interpolate_massive_and_agb_yields()
# Interpolate lifetimes
self.__interpolate_pop3_lifetimes()
self.__interpolate_massive_and_agb_lifetimes()
# Calculate coefficients to interpolate masses from lifetimes
self.__interpolate_pop3_m_from_t()
self.__interpolate_massive_and_agb_m_from_t()
# If radioactive isotopes are used ..
if self.len_decay_file > 0 and len(self.table_radio) > 0:
# Interpolate the radioactive yields tables
self.__interpolate_massive_and_agb_yields(is_radio=True)
# Check whether the initial metallicity is available
if (not self.iniZ in self.ytables.Z_list) and (self.iniZ > 0.0):
print ('Error - iniZ must be an available metallicity in the grid of stellar yields.')
self.need_to_quit = True
return
# Check for incompatible inputs - Error messages
self.__check_inputs()
if self.need_to_quit:
return
# NOTE: This if statement also needs to be in SYGMA and OMEGA!
# Initialisation of the composition of the gas reservoir
if self.is_sygma:
ymgal = np.zeros(self.nb_isotopes)
ymgal[0] = copy.deepcopy(self.mgal)
else:
ymgal = self._get_iniabu()
self.len_ymgal = len(ymgal)
# Initialisation of the storing arrays
mdot, ymgal, ymgal_massive, ymgal_agb, ymgal_1a, ymgal_nsm,\
ymgal_delayed_extra, mdot_massive, mdot_agb, mdot_1a, mdot_nsm,\
mdot_delayed_extra, sn1a_numbers, sn2_numbers, nsm_numbers,\
delayed_extra_numbers, imf_mass_ranges, \
imf_mass_ranges_contribution, imf_mass_ranges_mtot = \
self._get_storing_arrays(ymgal, len(self.history.isotopes))
# Initialisation of the composition of the gas reservoir
if len(self.ism_ini) > 0:
for i_ini in range(0,self.len_ymgal):
ymgal[0][i_ini] = self.ism_ini[i_ini]
# If radioactive isotopes are used ..
if self.len_decay_file > 0 or self.use_decay_module:
# Define initial radioactive gas composition
ymgal_radio = np.zeros(self.nb_radio_iso)
# Initialisation of the storing arrays for radioactive isotopes
mdot_radio, ymgal_radio, ymgal_massive_radio, ymgal_agb_radio,\
ymgal_1a_radio, ymgal_nsm_radio, \
ymgal_delayed_extra_radio, mdot_massive_radio, mdot_agb_radio,\
mdot_1a_radio, mdot_nsm_radio,\
mdot_delayed_extra_radio, dummy, dummy, dummy, dummy, \
dummy, dummy, dummy = \
self._get_storing_arrays(ymgal_radio, self.nb_radio_iso)
# Initialisation of the composition of the gas reservoir
if len(self.ism_ini_radio) > 0:
for i_ini in range(0,len(self.ism_ini_radio)):
ymgal_radio[0][i_ini] = self.ism_ini_radio[i_ini]
# Define indexes to make connection between unstable/stable isotopes
if not self.use_decay_module:
self.__define_unstab_stab_indexes()
# Output information
if self.iolevel >= 1:
print ('Number of timesteps: ', '{:.1E}'.format(len(timesteps)))
# Create empty arrays if on the fast mode
if self.pre_calculate_SSPs:
self.history.gas_mass.append(np.sum(ymgal[0]))
self.history.ism_iso_yield.append(ymgal[0])
self.history.m_locked = []
self.history.m_locked_agb = []
self.history.m_locked_massive = []
self.massive_ej_rate = []
self.sn1a_ej_rate = []
# Add the initialized arrays to the history class
else:
self.history.gas_mass.append(np.sum(ymgal[0]))
self.history.ism_iso_yield.append(ymgal[0])
self.history.ism_iso_yield_agb.append(ymgal_agb[0])
self.history.ism_iso_yield_1a.append(ymgal_1a[0])
self.history.ism_iso_yield_nsm.append(ymgal_nsm[0])
self.history.ism_iso_yield_massive.append(ymgal_massive[0])
self.history.sn1a_numbers.append(0)
self.history.nsm_numbers.append(0)
self.history.sn2_numbers.append(0)
self.history.m_locked = []
self.history.m_locked_agb = []
self.history.m_locked_massive = []
# Keep track of the mass-loss rate of massive stars and SNe Ia
self.massive_ej_rate = []
for k in range(self.nb_timesteps + 1):
self.massive_ej_rate.append(0.0)
self.sn1a_ej_rate = []
for k in range(self.nb_timesteps + 1):
self.sn1a_ej_rate.append(0.0)
# Attribute arrays and variables to the current object
self.mdot = mdot
self.ymgal = ymgal
self.ymgal_massive = ymgal_massive
self.ymgal_agb = ymgal_agb
self.ymgal_1a = ymgal_1a
self.ymgal_nsm = ymgal_nsm
self.ymgal_delayed_extra = ymgal_delayed_extra
self.mdot_massive = mdot_massive
self.mdot_agb = mdot_agb
self.mdot_1a = mdot_1a
self.mdot_nsm = mdot_nsm
self.mdot_delayed_extra = mdot_delayed_extra
self.sn1a_numbers = sn1a_numbers
self.nsm_numbers = nsm_numbers
self.delayed_extra_numbers = delayed_extra_numbers
self.sn2_numbers = sn2_numbers
self.imf_mass_ranges = imf_mass_ranges
self.imf_mass_ranges_contribution = imf_mass_ranges_contribution
self.imf_mass_ranges_mtot = imf_mass_ranges_mtot
# Attribute radioactive arrays and variables to the current object
if self.len_decay_file > 0 or self.use_decay_module:
self.mdot_radio = mdot_radio
self.ymgal_radio = ymgal_radio
self.ymgal_massive_radio = ymgal_massive_radio
self.ymgal_agb_radio = ymgal_agb_radio
self.ymgal_1a_radio = ymgal_1a_radio
self.ymgal_nsm_radio = ymgal_nsm_radio
self.ymgal_delayed_extra_radio = ymgal_delayed_extra_radio
self.mdot_massive_radio = mdot_massive_radio
self.mdot_agb_radio = mdot_agb_radio
self.mdot_1a_radio = mdot_1a_radio
self.mdot_nsm_radio = mdot_nsm_radio
self.mdot_delayed_extra_radio = mdot_delayed_extra_radio
# Declare non-metals for the getmetallicity function
self.nonmetals = ['H-','He-','Li-']
self.i_nonmetals = []
for i_iso in range(self.nb_isotopes):
if 'H-' in self.history.isotopes[i_iso] or\
'He-' in self.history.isotopes[i_iso] or\
'Li-' in self.history.isotopes[i_iso]:
self.i_nonmetals.append(i_iso)
self.len_i_nonmetals = len(self.i_nonmetals)
# Set the initial time and metallicity
# Hardcode the metallicity if Sygma, so it does not required
# a consistent iniabu file for all (unexpected) metallicities
if self.is_sygma:
zmetal = copy.deepcopy(self.iniZ)
else:
zmetal = self._getmetallicity(0)
self.history.metallicity.append(zmetal)
self.t = 0
self.history.age = np.zeros(self.nb_timesteps+1)
for i_tt in range(0,self.nb_timesteps):
self.history.age[i_tt+1] = self.history.age[i_tt] + \
self.history.timesteps[i_tt]
self.zmetal = zmetal
# Calculate stellar initial composition interpolation coefficients
if self.use_net_yields_stable:
self.__calculate_X0_coefs()