Skip to content
Snippets Groups Projects
Commit 99bd3f96 authored by Izzard, Robert Dr (Maths & Physics)'s avatar Izzard, Robert Dr (Maths & Physics)
Browse files

update cemp script to improve redshift calcs

parent a63991ee
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/perl/scripts2/cemps2017.pl
+ 242
− 37
View file @ 99bd3f96
...@@ -3,6 +3,8 @@ use strict; # recommended ...@@ -3,6 +3,8 @@ use strict; # recommended
use 5.16.0; # recommended use 5.16.0; # recommended
$|=1; $|=1;
use IO::Handle; use IO::Handle;
use Math::Trig;
use Sort::Key qw/nsort/;
############################################################ ############################################################
# #
...@@ -14,6 +16,27 @@ use IO::Handle; ...@@ -14,6 +16,27 @@ use IO::Handle;
my $outdir = "$ENV{HOME}/data/cemps2017"; my $outdir = "$ENV{HOME}/data/cemps2017";
mkdirhier($outdir); mkdirhier($outdir);
# cosmology
my $OmegaM = 0.27;
my $OmegaLambda = 0.73;
$OmegaM = 0.3;
$OmegaLambda = 0.7;
my $OmegaMLambda = $OmegaM / $OmegaLambda;
my $OmegaLambdaM = $OmegaLambda / $OmegaM;
my $sqrtOmegaMLambda = sqrt($OmegaMLambda);
my $sqrtOmegaLambdaM = sqrt($OmegaLambdaM);
my $asinhsqrtOmegaLambdaM = asinh(sqrt($OmegaLambdaM));
my $h = 68.0 * 1e5 / 3.08571e24; # H0 in cgs NB 1e5/3.086e24 = km/Mpc
my $agefac = 3.1557e7 * 1e9; # Gyr in s
my $agefach = 3.0/2.0 * $agefac * $h;
my $age_from_redshift_f0 = 2.0/3.0 * sqrt(1.0 + $OmegaMLambda) * (1.0/$h) / $agefac;
# load binary grid and rob's modules # load binary grid and rob's modules
use binary_grid2; # required use binary_grid2; # required
use binary_grid::C; # backend : C or Perl use binary_grid::C; # backend : C or Perl
...@@ -77,13 +100,14 @@ $population->parse_args(); ...@@ -77,13 +100,14 @@ $population->parse_args();
# duplicity is 0 for single stars, 1 for binary stars # duplicity is 0 for single stars, 1 for binary stars
# and 2 for a mixed population sampled at equal times # and 2 for a mixed population sampled at equal times
my $duplicity = 1; my $duplicity = 1;
my $mmin = 0.1;
if($duplicity == 0) if($duplicity == 0)
{ {
# make a grid of $nstars single binary stars, log-spaced, # make a grid of $nstars single binary stars, log-spaced,
# with masses between $mmin and $mmax # with masses between $mmin and $mmax
my $nstars = 40; my $nstars = 40;
my $mmin = 0.1;
my $mmax = 80.0; my $mmax = 80.0;
$population->add_grid_variable( $population->add_grid_variable(
'name' => 'lnm1', 'name' => 'lnm1',
...@@ -99,14 +123,14 @@ if($duplicity == 0) ...@@ -99,14 +123,14 @@ if($duplicity == 0)
elsif($duplicity == 1) elsif($duplicity == 1)
{ {
# make a population of binary stars # make a population of binary stars
my $n = 10; # n>4 gives you discs my $n = 40; # n>4 gives you discs
my $resolution = { my $resolution = {
m1 => $n, m1 => $n,
q => $n, q => $n,
P => $n P => $n
}; };
my $mmin = 0.4; my $mmin = $mmin;
my $mmax = 8; my $mmax = 8.0;
$population->{_grid_options}{binary} = 1; $population->{_grid_options}{binary} = 1;
$population->add_grid_variable $population->add_grid_variable
...@@ -120,37 +144,68 @@ elsif($duplicity == 1) ...@@ -120,37 +144,68 @@ elsif($duplicity == 1)
'probdist' =>'Kroupa2001($m1)*$m1', 'probdist' =>'Kroupa2001($m1)*$m1',
'dphasevol' =>'$dlnm1', 'dphasevol' =>'$dlnm1',
); );
# q=M1/M2 distribution flat in q between 0.1/M1 and 1.0 # q=M1/M2 distribution flat in q between $mmin/M1 and 1.0
$population->add_grid_variable $population->add_grid_variable
( (
'condition' =>'$self->{_grid_options}{binary}==1', 'condition' =>'$self->{_grid_options}{binary}==1',
'name' =>'q', 'name' =>'q',
'longname' =>'Mass ratio', 'longname' =>'Mass ratio',
'range' =>['0.1/$m1',1.0], 'range' =>[$mmin.'/$m1',1.0],
'resolution'=>$resolution->{q}, 'resolution'=>$resolution->{q},
'spacingfunc'=>"const(0.1/\$m1,1.0,$resolution->{q})", 'spacingfunc'=>"const($mmin/\$m1,1.0,$resolution->{q})",
'probdist' =>"flatsections\(\$q,\[ 'probdist' =>"flatsections\(\$q,\[
\{min=>0.1/\$m1,max=>0.8,height=>1.0\}, \{min=>$mmin/\$m1,max=>0.8,height=>1.0\},
\{min=>0.8,max=>1.0,height=>1.0\}, \{min=>0.8,max=>1.0,height=>1.0\},
\]\)", \]\)",
precode =>'$m2=$q*$m1;', precode =>'$m2=$q*$m1;',
dphasevol =>'$dq', dphasevol =>'$dq',
); );
# orbital period Duquennoy and Mayor 1991 distribution
$population->add_grid_variable
( # persepdist should be DQ91 or flatlna
'name' =>'logper', my $persepdist = 'DQ91';
'longname' =>'log(Orbital_Period)',
'range' =>[-2.0,12.0], if($persepdist eq 'DQ91')
'resolution' =>$resolution->{P}, {
'spacingfunc'=>"const(-1.0,10.0,$resolution->{P})", # orbital period Duquennoy and Mayor 1991 distribution
'precode'=>"\$per = 10.0 ** \$logper; $population->add_grid_variable
(
'name' =>'logper',
'longname' =>'log(Orbital_Period)',
'range' =>[-2.0,12.0],
'resolution' =>$resolution->{P},
'spacingfunc'=>"const(-1.0,10.0,$resolution->{P})",
'precode'=>"\$per = 10.0 ** \$logper;
my \$eccentricity = 0.0; my \$eccentricity = 0.0;
\$sep=calc_sep_from_period(\$m1,\$m2,\$per) if(defined \$m1 && defined \$m2); \$sep=calc_sep_from_period(\$m1,\$m2,\$per) if(defined \$m1 && defined \$m2);
", ",
'probdist'=>"gaussian(\$logper,4.8,2.3,-2.0,12.0)", 'probdist'=>"gaussian(\$logper,4.8,2.3,-2.0,12.0)",
'dphasevol'=>'$dln10per' 'dphasevol'=>'$dln10per'
); );
}
elsif($persepdist eq 'flatlna')
{
$population->add_grid_variable(
# name of the variable
'name'=>'lnsep',
'longname'=>'ln(Orbital_Separation)',
'range'=>['log(3.0)','log(1e4)'],
'resolution'=>$resolution->{P},
'spacingfunc'=>"const(log(3.0),log(1e4),$resolution->{P})",
# precode has to calculation the period (required for binary_c)
'precode'=>"my \$eccentricity = 0.0;my\$sep=exp(\$lnsep);my\$per=calc_period_from_sep(\$m1,\$m2,\$sep);",
# flat in log-separation (dN/da~1/a) distribution (Opik law)
'probdist'=>'const(log(3.0),log(1e4))',
# phase volume contribution
'dphasevol'=>'$dlnsep'
);
}
else
{
print STDERR "Unknown persepdist = $persepdist\n";
die;
}
} }
else else
{ {
...@@ -170,7 +225,7 @@ else ...@@ -170,7 +225,7 @@ else
time_adaptive_mass_grid_nlow_mass_stars=>10, time_adaptive_mass_grid_nlow_mass_stars=>10,
nthreads =>$nthreads, nthreads =>$nthreads,
thread_sleep =>1, thread_sleep =>1,
mmin =>0.1, mmin =>$mmin,
mmax =>80.0, mmax =>80.0,
mass_grid_log10_time=>0, mass_grid_log10_time=>0,
mass_grid_step =>$dt*$sampling_factor, mass_grid_step =>$dt*$sampling_factor,
...@@ -185,9 +240,9 @@ else ...@@ -185,9 +240,9 @@ else
}; };
distribution_functions::bastard_distribution( distribution_functions::bastard_distribution(
$population, { $population, {
mmin =>0.1, mmin =>$mmin,
mmax =>80.0, mmax =>80.0,
m2min =>0.1, m2min =>$mmin,
nm2 =>10, nm2 =>10,
nper =>10, nper =>10,
qmin =>0.0, qmin =>0.0,
...@@ -212,7 +267,6 @@ my @nuclear_mass_list = $population->nuclear_mass_list(); ...@@ -212,7 +267,6 @@ my @nuclear_mass_list = $population->nuclear_mass_list();
my @sources = $population->source_list(); my @sources = $population->source_list();
my @ensemble = $population->ensemble_list(); my @ensemble = $population->ensemble_list();
# check code features # check code features
#$population->require_not(); #$population->require_not();
$population->require( $population->require(
...@@ -220,6 +274,37 @@ $population->require( ...@@ -220,6 +274,37 @@ $population->require(
'NUCSYN_CEMP_LOGGING is on' 'NUCSYN_CEMP_LOGGING is on'
); );
# age sampling (Gyr)
my @ages = (0.0,1,10,11,12,13,13.3,13.4,13.5,13.6,13.7,13.8,13.9,14);
my $maxage = age_from_redshift(1e10);
print "Universe age $maxage\n";
@ages = grep { $_ < $maxage} @ages;
# make a set of redshifts corresponding to the ages
my @redshifts = make_redshifts(\@ages);
foreach my $redshift (@redshifts)
{
printf "z = %g -> age = %g\n",
$redshift,
age_from_redshift($redshift);
}
if(0)
{
# test cosmology functions
foreach my $z (0,1,2,5,10,100,1e3,1e4)
{
my $age = age_from_redshift($z);
printf "Age at redshift $z = %g : inverse gives %g\n",
$age,
redshift_from_age($age);
}
exit;
}
# evolution the stellar population (this takes some time) # evolution the stellar population (this takes some time)
$population->evolve(); $population->evolve();
...@@ -237,6 +322,9 @@ sub parse_data ...@@ -237,6 +322,9 @@ sub parse_data
{ {
my ($population, $results) = @_; my ($population, $results) = @_;
my $m1init = $population->{_grid_options}->{progenitor_hash}->{m1};
my $p_Kroupa = distribution_functions::Kroupa2001($m1init);
while(1) while(1)
{ {
# subsequent calls to tbse_line contain # subsequent calls to tbse_line contain
...@@ -277,11 +365,54 @@ sub parse_data ...@@ -277,11 +365,54 @@ sub parse_data
my $dtp = $dt * $p; my $dtp = $dt * $p;
$results->{count}->{$header} += $dtp; # now repeat for different redshifts
foreach my $z (@redshifts)
{
my $p_Adam = adam_imf($m1init,$z);
my $factor = $p_Adam / $p_Kroupa;
if(0)
{
printf "IMF(m=%g, z=%g, age=%g Gyr) = K %g A %g : f = %g\n",
$m1init,
$z,
age_from_redshift($z),
$p_Kroupa,
$p_Adam,
$factor;
}
$results->{count}->{$z}->{$header} += $dtp * $factor;
}
} }
} }
} }
sub adam_imf
{
my ($m,$z) = @_; # mass, redshift
my $p;
my $f = adam_f($z);
return distribution_functions::three_part_power_law
(
$m,
$mmin,
0.08 * $f,
0.5 * $f,
80.0,
-0.3,
-1.3,
-2.3);
}
sub adam_f
{
my ($z) = @_;
my $r = MAX(1.0, $z/6.3);
return $r * $r;
}
############################################################ ############################################################
...@@ -312,21 +443,33 @@ sub output ...@@ -312,21 +443,33 @@ sub output
close $fp; close $fp;
} }
print "----------\n";
my $tot_EMP = 0.0; my $tot_EMP = {};
my @counttypes = sort keys %{$results->{count}};
foreach my $counttype (@counttypes) foreach my $z (keys %{$results->{count}})
{ {
$tot_EMP += $results->{count}->{$counttype}; $tot_EMP->{$z} //= 0.0;
my @counttypes = sort keys %{$results->{count}->{$z}};
foreach my $counttype (@counttypes)
{
$tot_EMP->{$z} += $results->{count}->{$z}->{$counttype};
}
$tot_EMP->{$z} = MAX(1e-100,$tot_EMP->{$z}); # prevent 1/0
} }
$tot_EMP = MAX(1e-100,$tot_EMP); # prevent 1/0
foreach my $counttype (@counttypes) foreach my $z (nsort keys %{$results->{count}})
{ {
printf "% 10s %5.2f %% %g\n", print "----------\n";
$counttype, printf "Redshift %g, age %g\n\n",$z,age_from_redshift($z);
$results->{count}->{$counttype} * 100.0/ $tot_EMP, my @counttypes = sort keys %{$results->{count}->{$z}};
$results->{count}->{$counttype} ; foreach my $counttype (@counttypes)
} {
printf "% 10s %5.2f %% %g\n",
$counttype,
$results->{count}->{$z}->{$counttype} * 100.0/ $tot_EMP->{$z},
$results->{count}->{$z}->{$counttype} ;
}
}
print "----------\n"; print "----------\n";
} }
...@@ -335,3 +478,65 @@ sub safelog10 ...@@ -335,3 +478,65 @@ sub safelog10
# safe log10 function (prevents log(0)) # safe log10 function (prevents log(0))
return log10(MAX(1e-30,$_[0])); return log10(MAX(1e-30,$_[0]));
} }
sub make_redshifts
{
my ($ages) = @_;
my @redshifts;
foreach my $age (@$ages)
{
push(@redshifts, redshift_from_age($age));
}
return @redshifts;
}
############################################################
# cosmology
############################################################
sub redshift_from_age
{
# Eq.13 : redshfift from age in Gyr
my ($age) = @_; # Gyr
return 0.0 if ($age == 0.0);
# Weinberg 1989 Eq. 5.4 inverted
# Valid for a flat Universe
my $z = (
$sqrtOmegaMLambda
*
sinh(
$asinhsqrtOmegaLambdaM
-
$age * $agefach /
sqrt(1.0 + $OmegaMLambda)
))**(-2.0/3.0)
-
1.0;
return $z;
}
sub age_from_redshift
{
# inverse of the above, return in Gyr
my ($z) = @_;
return 0.0 if ($z == 0.0);
# Weinberg 1989 Eq. 5.4
# Valid for a flat Universe
my $age =
$age_from_redshift_f0 *
(
$asinhsqrtOmegaLambdaM
-
asinh(sqrt($OmegaLambdaM) * (1.0 + $z)**-1.5)
);
return $age;
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment