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

add Zhengwei script

parent 3102aa32
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/perl/scripts2/ZW2019.pl 0 → 100755
+ 366
− 0
View file @ 0522dcb5
#!/usr/bin/env perl
use strict; # recommended
use 5.16.0; # recommended
use binary_grid2; # required
use binary_grid::C; # backend : C or Perl
use rob_misc qw/ncpus/;
use IO::Handle;
############################################################
#
# Example script to demonstrate how to use the
# binary_grid2 module.
#
# For full documentation, please see binary_grid2.pdf
# in the doc/ directory of binary_c
#
############################################################
# number of computational threads to launch, usually one
# per CPU
my $nthreads = rob_misc::ncpus();
############################################################
# Binary_c should output data that we can understand. There
# are two ways to do this:
#
# 1) Put output statements, using PRINTF, in binary_c's
# log_every_timestep() function.
#
# 2) Put output statements, using PRINTF, into the C_logging_code
# grid option
#
# $population->set( C_logging_code => ' PRINTF("...\n"); ');
#
# You have access to the stardata variable, so you can
# output everything that is available to log_every_timestep();
#
############################################################
# make a new stellar population
my $population = binary_grid2->new(
# options can be given now ...
metallicity => 0.02, # mass fraction of "metals"
max_evolution_time => 15000, # Myr
# gamma :
# -1 = RLOF_NONCONSERVATIVE_GAMMA_DONOR
# -2 = RLOF_NONCONSERVATIVE_GAMMA_ISOTROPIC
#nonconservative_angmom_gamma => 'RLOF_NONCONSERVATIVE_GAMMA_DONOR',
nthreads=>$nthreads, # number of threads
);
# ... or options can be set manually later.
$population->set(
vb=>1, # turn on verbose logging (can be 0,1,2...)
return_array_refs=>1, # quicker data parsing mode
log_args=>1,
sort_args=>1,
save_args=>1,
log_args_dir=>'/tmp',
threads_entry_function_pointer => \&thread_entry,
threads_flush_function_pointer => \&thread_flush,
# log_filename => '/dev/stdout',
# api_log_filename_prefix => '/dev/null',
## or enter more complicated code yourself:
#
# C_logging_code => '
# PRINTF("MY_STELLAR_DATA %g %g %g %g\n",
# stardata->model.time,
# stardata->star[0].mass,
# stardata->model.probability,
# stardata->model.dt);
# '
);
$population->require_on(
'ZW2019_SNIA_LOG'
);
############################################################
# scan command line arguments for extra options
$population->parse_args();
############################################################
#
# Now we set up the physics of the stellar evolution and
# population runs
#
#
# duplicity is 0 for single stars, 1 for binary stars
# and 2 for a mixed population sampled at equal times
my $duplicity = 1;
my $mmin = 0.8; # 0.8 for the run, 0.1 for the whole population
my $mmax = 8.0; # 8 for the run, 80 for the whole population
if($duplicity == 0)
{
# make a grid of $nstars single binary stars, log-spaced,
# with masses between $mmin and $mmax
my $nstars = 100;
$population->add_grid_variable(
'name' => 'lnm1',
'longname' =>'Primary mass',
'range' =>[log($mmin),log($mmax)],
'resolution' => $nstars, # just a counter for the grid
'spacingfunc'=>"const(log($mmin),log($mmax),$nstars)",
'precode' =>'$m1=exp($lnm1);',
'probdist' =>'Kroupa2001($m1)*$m1',
'dphasevol' =>'$dlnm1',
);
}
elsif($duplicity == 1)
{
# make a population of binary stars
my $res = 100;
my $resolution = {
m1 => $res,
q => $res,
P => $res
};
$population->{_grid_options}{binary} = 1;
$population->add_grid_variable
(
'name' => 'lnm1',
'longname' =>'Primary mass',
'range' =>[log($mmin),log($mmax)],
'resolution' => $resolution->{m1},
'spacingfunc'=>"const(log($mmin),log($mmax),$resolution->{m1})",
'precode' =>'$m1=exp($lnm1);',
'probdist' =>'Kroupa2001($m1)*$m1',
'dphasevol' =>'$dlnm1',
);
# q=M1/M2 distribution flat in q between 0.1/M1 and 1.0
$population->add_grid_variable
(
'condition' =>'$self->{_grid_options}{binary}==1',
'name' =>'q',
'longname' =>'Mass ratio',
'range' =>['0.1/$m1',1.0],
'resolution'=>$resolution->{q},
'spacingfunc'=>"const(0.1/\$m1,1.0,$resolution->{q})",
'probdist' =>"flatsections\(\$q,\[
\{min=>0.1/\$m1,max=>0.8,height=>1.0\},
\{min=>0.8,max=>1.0,height=>1.0\},
\]\)",
precode =>'$m2=$q*$m1;',
dphasevol =>'$dq',
);
# orbital period Duquennoy and Mayor 1991 distribution
$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;
\$sep=calc_sep_from_period(\$m1,\$m2,\$per) if(defined \$m1 && defined \$m2);
",
'probdist'=>"gaussian(\$logper,4.8,2.3,-2.0,12.0)",
'dphasevol'=>'$dln10per'
);
}
else
{
# Make a mixed population of single and binary stars
# with equal sampling in equal times.
my $dt = 100.0;
$population->set(
time_adaptive_mass_grid_step => $dt,
);
my $sampling_factor = 0.5; # deliberate oversampling (<0.5)
my $time_adaptive_options = {
max_evolution_time =>
$population->{_bse_options}{max_evolution_time},
stellar_lifetime_table_nm=>1000,
time_adaptive_mass_grid_log10_time=>0,
time_adaptive_mass_grid_step=>$dt,
time_adaptive_mass_grid_nlow_mass_stars=>10,
nthreads =>$nthreads,
thread_sleep =>1,
mmin =>0.1,
mmax =>80.0,
mass_grid_log10_time=>0,
mass_grid_step =>$dt*$sampling_factor,
extra_flash_resolution=>0, # 1 = broken?
mass_grid_nlow_mass_stars=>10,
debugging_output_directory=>'/tmp',
max_delta_m =>1.0,
savegrid =>1,
vb =>0,
metallicity=>$population->{_bse_options}{metallicity},
agbzoom =>0,
};
distribution_functions::bastard_distribution(
$population, {
mmin =>0.1,
mmax =>80.0,
m2min =>0.1,
nm2 =>10,
nper =>10,
qmin =>0.0,
qmax =>1.0,
necc =>undef,
useecc =>undef,
agbzoom =>0,
time_adaptive=>$time_adaptive_options,
});
}
############################################################
# link population to custom data parser function
$population->set(
parse_bse_function_pointer => \&main::parse_data
);
my %init = $population->initial_abundance_hash('Karakas2002',0.02);
my %isotope_hash = $population->isotope_hash();
my @isotope_list = $population->isotope_list();
my %nuclear_mass_hash = $population->nuclear_mass_hash();
my @nuclear_mass_list = $population->nuclear_mass_list();
my @sources = $population->source_list();
my @ensemble = $population->ensemble_list();
# you can use Data::Dumper to see the contents
# of the above lists and hashes
if(0){
print Data::Dumper->Dump([
#\%init,
#\%isotope_hash,
#\@isotope_list,
#\%nuclear_mass_hash,
\@nuclear_mass_list,
#\@sources,
#\@ensemble
]);
}
if(0)
{
my $h =
$population->run_gridcode_fast(
{
keep_probabilities=>1,
keep_phasevol=>1,
system_code=>'$h->{mass_into_stars} += ($m1+$m2)*$prob;',
}
);
print "Mass into stars $h->{mass_into_stars}\n";
exit;
}
# uncomment this to show version information
#print $population->evcode_version_string();
# uncomment this to show the evcode's args list
#print join("\n",@{$population->evcode_args_list()});
# evolution the stellar population (this takes some time)
$population->evolve();
# output the data
output($population);
# done : exit
exit;
############################################################
# subroutines
############################################################
sub parse_data
{
my ($population, $results) = @_;
my $ZW_data;
while(1)
{
# subsequent calls to tbse_line contain
# (references to) arrays of data
my $la = $population->tbse_line();
# first element is the "header" line
my $header = shift @$la;
# break out of the look if this is 'fin'
last if ($header eq 'fin');
# check if $header matches one of your
# expected data lines, if so, act
if($header eq 'ZW2019')
{
if(defined $ZW_data)
{
$ZW_data .= "@$la\n";
}
else
{
$ZW_data = "@$la\n";
}
}
}
if(defined $ZW_data)
{
#print {$population->{_ZW_file}} "\n",$population->{_grid_options}->{args},"\n";
print {$population->{_ZW_file}} $ZW_data,"\n";
}
}
sub thread_entry
{
my ($population) = @_;
print "Thread entry $population\n";
my $fname = '/tmp/ZW_data.'.$population->flexigrid_thread_number();
open($population->{_ZW_file}, '>', $fname) || die("cannot open $fname");
$population->{_ZW_file}->autoflush(1);
print "opened $fname\n";
$population->{_ZW_fname} = $fname;
}
sub thread_flush
{
my ($population) = @_;
# close ZW_data file
print "close $population->{_ZW_fname}\n";
close($population->{_ZW_file});
}
############################################################
sub output
{
my ($population) = @_;
# $results is a hash reference containing
# the results that were added up in parse_data()
my $results = $population->results;
my $version_string = $population->commented_evcode_version_string();
my $info_string = $population->commented_infostring();
my $outfile = ("@ARGV"=~/datafile=(\S+)/)[0] // '/tmp/ZW_data';
# prepend with run and version information
open(my $ZWout,'>',$outfile)||die("cannot open $outfile for writing");
print {$ZWout} $version_string,"\n#\n#\n#\n",$info_string,"\n\n\n";
close $ZWout;
# join the data from each thread
my $cat = 'cat ';
foreach my $n (0..$nthreads)
{
$cat .= '/tmp/ZW_data.'.$n.' ';
}
$cat .= ' >> '.$outfile;
print "Running $cat\n";
`$cat`;
print "See $outfile\n";
}
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