diff --git a/src/disc/disc_logging.c b/src/disc/disc_logging.c
index 7e42a1870192dc3862dfb5d328204ed07761899b..1dcd5385241b941d7c2b8710c42e7a7c470ea029 100644
--- a/src/disc/disc_logging.c
+++ b/src/disc/disc_logging.c
@@ -196,7 +196,7 @@ void disc_logging(struct stardata_t * const stardata,
double Revap_in = disc->Revap_in;
/* output */
- Printf("DISC_EV %20.10e %20.10e %d %d %u %d %d %d %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g\n",
+ Printf("DISC_EV %20.10e %20.10e %d %d %u %d %d %d %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g\n",
dtp, // 1
disc->lifetime/YEAR_LENGTH_IN_SECONDS,
stardata->star[0].stellar_type,
@@ -248,7 +248,10 @@ void disc_logging(struct stardata_t * const stardata,
disc_total_gravitational_potential_energy(disc,b),
disc_total_kinetic_energy(disc,b),
disc->sigma0, // 51
- disc->Tvisc0
+ disc->Tvisc0,
+ disc->Rout / disc->Rin,
+ disc->Rin / Max(1e-50,b->separation),
+ disc->Rout / Max(1e-50,b->separation)
);
}
diff --git a/src/maths/GSL_integrator.c b/src/maths/GSL_integrator.c
index 7e2720f1fda8da3943538356c7dd9cc743c7cd1b..0f7fae7ee1857f81fab79044ad8da487825f1c13 100644
--- a/src/maths/GSL_integrator.c
+++ b/src/maths/GSL_integrator.c
@@ -1,11 +1,11 @@
#include "../binary_c.h"
#include <gsl/gsl_integration.h>
-double GSL_integrator(double lower,
- double upper,
- double relative_tolerance,
- double * error,
- int integrator,
+double GSL_integrator(const double lower,
+ const double upper,
+ const double relative_tolerance,
+ double * const error,
+ const int integrator,
double func(double x,void *params),
...)
{
diff --git a/src/maths/maths_prototypes.h b/src/maths/maths_prototypes.h
index e6aaa0914d2097447e80ca15bbe2e234793a0d23..45d39af1418a058a74d4c09f3e2173fc703d630e 100644
--- a/src/maths/maths_prototypes.h
+++ b/src/maths/maths_prototypes.h
@@ -156,11 +156,11 @@ double polyalgorithm(const int s,
Boolean inverse3( const double theMatrix [/*Y=*/3] [/*X=*/3],
double theOutput [/*Y=*/3] [/*X=*/3] );
-double GSL_integrator(double lower,
- double upper,
- double tolerance,
- double * error,
- int integrator,
+double GSL_integrator(const double lower,
+ const double upper,
+ const double tolerance,
+ double * const error,
+ const int integrator,
double func(double x,void *params),
...);
void setup_GSL_handlers(struct stardata_t * RESTRICT const stardata);
diff --git a/src/perl/scripts2/cbdiscs.pl b/src/perl/scripts2/cbdiscs.pl
index 5e683ffbdbaafa4330e0fcc4081a6625b25eb2c9..78bac5a587a68ce9d916819c30d3be0b1027dc81 100755
--- a/src/perl/scripts2/cbdiscs.pl
+++ b/src/perl/scripts2/cbdiscs.pl
@@ -105,7 +105,7 @@ $population->set(
cbdisc_alpha => 1e-3,
cbdisc_kappa => 1e-2,
cbdisc_torquef => 1e-4, # 1e-4
-
+
cbdisc_init_dM => 0.0, # 0 means do not adapt initial masses
cbdisc_init_dJdM => 0.5,
@@ -148,15 +148,15 @@ if($ARGV[0] eq 'show_opts')
}
$population->set(
- # slurm
- slurm=>0, # use slurm if 1 (please set on command line)
- slurm_njobs=>4, # number of jobs
- slurm_dir=>'/users/ri0005/data/slurm', # working directory for scripts etc.
- slurm_memory=>512, # RAM (MB) per job
- slurm_postpone_join=>0, # if 1 do not join on slurm, join elsewhere
- slurm_partition=>'all', # MUST be defined
- slurm_jobname=>'binary_cbd', # not required but useful
- slurm_use_all_node_CPUs=>1, # required on eureka
+ # slurm
+ slurm=>0, # use slurm if 1 (please set on command line)
+ slurm_njobs=>4, # number of jobs
+ slurm_dir=>'/users/ri0005/data/slurm', # working directory for scripts etc.
+ slurm_memory=>512, # RAM (MB) per job
+ slurm_postpone_join=>0, # if 1 do not join on slurm, join elsewhere
+ slurm_partition=>'all', # MUST be defined
+ slurm_jobname=>'binary_cbd', # not required but useful
+ slurm_use_all_node_CPUs=>1, # required on eureka
);
# scan command line arguments for extra options
@@ -198,62 +198,63 @@ elsif($duplicity == 1)
q => $n,
P => $n
};
+
my $mmin = 0.85;
my $mmax = 6;
$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',
+ '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],
+ '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',
+ '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],
- 'range' =>[1.0,4.0],
- 'resolution' =>$resolution->{P},
- 'spacingfunc'=>"const(-1.0,10.0,$resolution->{P})",
- 'precode'=>"
+ # orbital period Duquennoy and Mayor 1991 distribution
+ $population->add_grid_variable
+ (
+ 'name' =>'logper',
+ 'longname' =>'log(Orbital_Period)',
+ # 'range' =>[-2.0,12.0],
+ 'range' =>[1.0,4.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);
+ \$per = 10.0 ** \$logper;
+ my \$eccentricity = 0.0;
+ \$sep=calc_sep_from_period(\$m1,\$m2,\$per) if(defined \$m1 && defined \$m2);
-",
+ ",
- 'postcode' => "
+ 'postcode' => "
-",
+ ",
- 'probdist'=>"gaussian(\$logper,4.8,2.3,-2.0,12.0)",
- 'dphasevol'=>'$dln10per'
- );
+ 'probdist'=>"gaussian(\$logper,4.8,2.3,-2.0,12.0)",
+ 'dphasevol'=>'$dln10per'
+ );
}
else
{
@@ -327,6 +328,214 @@ $population->require(
'BISECT_FAIL_EXIT is off',
);
+
+# list of variables in each data line, in order
+# (note that dtp is the first data item, always)
+#
+# those that are associated with undef are not plotted
+#
+# those given an empty hash are plotted with no options
+#
+# otherwise a hash {...} can contain directives for each plot,
+# e.g. whether data is binned, logged, multiplied or prepended
+# with 'log'
+my $vars = [
+ age => {
+ bin => 0.1,
+ log => 1
+ },
+ stellar_type0 => {},
+ stellar_type1 => {},
+ nzones => undef,
+ zone1_type => undef,
+ zone2_type => undef,
+ zone3_type => undef,
+ Rin => {
+ bin => 0.025,
+ log => 1,
+ prepend => 1,
+ },
+ Rout => {
+ bin => 0.025,
+ log => 1,
+ prepend => 1,
+ },
+ separation => {
+ bin => 0.025,
+ log => 1,
+ },
+ orbital_period => {
+ bin => 0.1, # years
+ multiplier => 365.25,
+ log => 1
+ },
+ Rism => undef,
+ Xray_domination => undef,
+ R0 => undef,
+ R1 => undef,
+ RL0 => undef,
+ RL1 => undef,
+ Ldisc =>
+ {
+ bin => 0.025,
+ log => 1,
+ prepend => 1,
+ },
+ Idisc => undef,
+ Jdisc => {
+ bin => 0.1,
+ log => 1,
+ prepend => 1,
+ },
+ Jbinary => {
+ bin => 0.025,
+ log => 1,
+ prepend => 1,
+ },
+ Mdisc => {
+ bin => 0.1,
+ log => 1,
+ prepend => 1,
+ },
+ Mbinary => {
+ bin => 0.025,
+ log => 1,
+ prepend => 1,
+ },
+ Mdot_viscous => {
+ bin => 0.1,
+ log => 1,
+ },
+ Mdot_global => {
+ bin => 0.1,
+ log => 1,
+ },
+ Mdot_inner_L2 => {
+ bin => 0.1,
+ log => 1,
+ },
+ Mdot_ISM => {
+ bin => 0.1,
+ log => 1,
+ },
+ Mdot_FUV => {
+ bin => 0.1,
+ log => 1,
+ },
+ Mdot_Xray => {
+ bin => 0.1,
+ log => 1,
+ },
+ eccentricity => {
+ bin => 0.025,
+ },
+ edot_resonance => {
+ bin => 0.1,
+ },
+ Tin => {
+ bin => 0.025,
+ log => 1,
+ prepend => 1,
+ },
+ Tout => {
+ bin => 0.025,
+ log => 1,
+ prepend => 1,
+ },
+ densityin => {
+ bin => 0.1,
+ log => 1,
+ prepend => 1,
+ },
+ densityout => {
+ bin => 0.1,
+ log => 1,
+ prepend => 1,
+ },
+ Pgrav_in => {
+ bin => 0.1,
+ log => 1,
+ prepend => 1,
+ },
+ Pgrav_out => {
+ bin => 0.1,
+ log => 1,
+ prepend => 1,
+ },
+ PISM => undef,
+ HRin => {
+ bin => 0.025,
+ log => 1,
+ prepend => 1,
+ },
+ HRout => {
+ bin => 0.025,
+ log => 1,
+ prepend => 1,
+ },
+ Teffmax => {
+ bin => 0.025,
+ log => 1,
+ prepend => 1,
+ },
+ Tin_over_Tstar => {
+ bin => 0.1,
+ },
+ Lstar => {
+ bin => 0.025,
+ prepend => 1,
+ log => 1,
+ },
+ LXstar => {
+ bin => 0.1,
+ prepend => 1,
+ log => 1,
+ },
+ Ldisc_over_Lstar => undef,
+ tvisc_in => {
+ bin => 0.1,
+ prepend => 1,
+ log => 1,
+ },
+ tvisc_out => {
+ bin => 0.1,
+ prepend => 1,
+ log => 1,
+ },
+ PE => undef,
+ KE => undef,
+ sigma0 => undef,
+ Tvisc0 => undef,
+ Rout_over_Rin => {
+ bin => 0.025,
+ log => 1,
+ },
+ Rin_over_separation => {
+ bin => 0.025,
+ log => 1,
+ },
+ Rout_over_separation => {
+ bin => 0.025,
+ log => 1,
+ },
+ ];
+
+# 2D plots : binning, logging etc. is as in 1D
+my $twoD = {
+ 'period_eccentricity' =>
+ {
+ x => 'orbital_period',
+ y => 'eccentricity',
+ },
+
+ 'period_luminosity' =>
+ {
+ x => 'orbital_period',
+ y => 'logLstar',
+ },
+};
+
+
# evolution the stellar population (this takes some time)
$population->evolve();
@@ -346,13 +555,16 @@ sub parse_data
my $progenitor_stellar_type;
my $disc_donor_stellar_type;
-
+
+
+
+
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;
@@ -366,116 +578,78 @@ sub parse_data
#print "$header\n";# @$la\n", $population->{_grid_options}->{args},"\n";
if($header eq 'DISC_EV')
{
- #print "DISC_EV @$la\n";
- my ($dtp,
- $t,
- $stellar_type0,
- $stellar_type1,
- $nzones,
- $zone1_type,
- $zone2_type,
- $zone3_type,
- $Rin,
- $Rout,
- $separation, # Rsun
- $orbital_period, # years
- $Rism,
- $Xray_domination,
- $R0,
- $R1,
- $RL0,
- $RL1,
- $Ldisc,
- $Idisc,
- $Jdisc,
- $Jbinary,
- $Mdisc,
- $Mbinary,
- $Mdot_viscous,
- $Mdot_global,
- $Mdot_inner_L2,
- $Mdot_ISM,
- $Mdot_FUV,
- $Mdot_Xray,
- $eccentricity,
- $edot_resonance,
- $Tin,
- $Tout,
- $density_in,
- $density_out,
- $Pgrav_in,
- $Pgrav_out,
- $PISM,
- $HRin,
- $HRout,
- $Teffmax,
- $Tin_over_Tstar,
- $Lstar,
- $LXstar,
- $Ldisc_over_Lstar,
- $tvisc_in,
- $tvisc_out,
- $PE,
- $KE,
- $sigma0,
- $Tvisc0) = @$la;
-
+ # first is always dtp : remove this
+ my $dtp = shift @$la;
+ #print "DISC EV : dtp = $dtp\n";
+ #print "ARRAY @$la\n";
+
+ # construct the $data hash
+ my $data = {};
+ my $c = 0;
+ while(scalar @$la)
+ {
+ my $y = shift @$la;
+ my $pre = '';
+ my $k = $vars->[$c++]; # variable key
+ my $h = $vars->[$c++]; # variable hash
+
+ if(defined $h)
+ {
+ # extra things to do...
+
+ # multiplier
+ if(defined $h->{multiplier})
+ {
+ $y *= $h->{multiplier};
+ }
+
+ # log10
+ if(defined $h->{log} && $h->{log})
+ {
+ $y = safelog10($y);
+
+ # prepend 'log' to the data key
+ if($h->{prepend})
+ {
+ $pre = 'log';
+ }
+ }
+
+ # bin data
+ if(defined $h->{bin})
+ {
+ $y = $population->rebin($y,$h->{bin});
+ }
+
+ # set the data hash : note that these data are unique
+ # there should be no existing data in $data->{$pre.$k}
+ $data->{$pre.$k}->{$y} = $dtp;
+ }
+ }
+
my $r;
foreach my $set ('all',
$progenitor_stellar_type)
{
$r = $results->{$set};
+
+ # set the data
+ foreach my $x (keys %$data)
+ {
+ foreach my $y (keys %{$data->{$x}})
+ {
+ $results->{$set}->{$x}->{$y} += $data->{$x}->{$y};
+ }
+ }
- # prebin some quantities
- # we want orbital period in days
- my $orbital_period_binned = $population->rebin(safelog10($orbital_period*365.25),0.1);
- my $eccentricity_binned = $population->rebin($eccentricity,0.025);
- my $logLstar_binned = $population->rebin(safelog10($Lstar),0.025);
-
- # 1D distributions
- $r->{Mdot_viscous}->{$population->rebin(safelog10($Mdot_viscous),0.1)} += $dtp;
- $r->{Mdot_global}->{$population->rebin(safelog10($Mdot_global),0.1)} += $dtp;
- $r->{Mdot_inner_L2}->{$population->rebin(safelog10($Mdot_inner_L2),0.1)} += $dtp;
- $r->{Mdot_ISM}->{$population->rebin(safelog10($Mdot_ISM),0.1)} += $dtp;
- $r->{Mdot_FUV}->{$population->rebin(safelog10($Mdot_FUV),0.1)} += $dtp;
- $r->{Mdot_Xray}->{$population->rebin(safelog10($Mdot_Xray),0.1)} += $dtp;
- $r->{logMdisc}->{$population->rebin(safelog10($Mdisc),0.1)} += $dtp;
- $r->{logMbinary}->{$population->rebin(safelog10($Mbinary),0.025)} += $dtp;
- $r->{logLdisc}->{$population->rebin(safelog10($Ldisc),0.025)} += $dtp;
- $r->{logJdisc}->{$population->rebin(safelog10($Jdisc),0.1)} += $dtp;
- $r->{logJbinary}->{$population->rebin(safelog10($Jbinary),0.025)} += $dtp;
- $r->{age}->{$population->rebin(safelog10($t),0.1)} += $dtp;
- $r->{separation}->{$population->rebin(safelog10($separation),0.025)} += $dtp;
- $r->{orbital_period}->{$orbital_period_binned} += $dtp;
- $r->{eccentricity}->{$eccentricity_binned} += $dtp;
- $r->{edot_resonance}->{$population->rebin($edot_resonance,0.1)} += $dtp;
- $r->{stellar_type0}->{$stellar_type0} += $dtp;
- $r->{stellar_type1}->{$stellar_type1} += $dtp;
- $r->{logRin}->{$population->rebin(safelog10($Rin),0.025)} += $dtp;
- $r->{logRout}->{$population->rebin(safelog10($Rout),0.025)} += $dtp;
- $r->{logRout_over_Rin}->{$population->rebin(safelog10($Rout/MAX(1e-50,$Rin)),0.025)} += $dtp;
- $r->{logRin_over_separation}->{$population->rebin(safelog10($Rin/MAX(1e-50,$separation)),0.025)} += $dtp;
- $r->{logRout_over_separation}->{$population->rebin(safelog10($Rout/MAX(1e-50,$separation)),0.025)} += $dtp;
- $r->{logHRin}->{$population->rebin(safelog10($HRin),0.025)} += $dtp;
- $r->{logHRout}->{$population->rebin(safelog10($HRout),0.025)} += $dtp;
- $r->{logLdisc}->{$population->rebin(safelog10($Ldisc),0.1)} += $dtp;
- $r->{logLstar}->{$logLstar_binned} += $dtp;
- $r->{logLXstar}->{$population->rebin(safelog10($LXstar),0.1)} += $dtp;
- $r->{logdensityin}->{$population->rebin(safelog10($density_in),0.1)} += $dtp;
- $r->{logdensityout}->{$population->rebin(safelog10($density_out),0.1)} += $dtp;
- $r->{logTin}->{$population->rebin(safelog10($Tin),0.025)} += $dtp;
- $r->{logTout}->{$population->rebin(safelog10($Tout),0.025)} += $dtp;
- $r->{logPgrav_in}->{$population->rebin(safelog10($Pgrav_in),0.1)} += $dtp;
- $r->{logPgrav_out}->{$population->rebin(safelog10($Pgrav_out),0.1)} += $dtp;
- $r->{logTeffmax}->{$population->rebin(safelog10($Teffmax),0.025)} += $dtp;
- $r->{logtvisc_in}->{$population->rebin(safelog10($tvisc_in),0.1)} += $dtp;
-
# 2D plots
- $r->{'2D_period_eccentricity'}->{$orbital_period_binned}->{$eccentricity_binned} += $dtp;
- $r->{'2D_period_luminosity'}->{$orbital_period_binned}->{$logLstar_binned} += $dtp;
-
- #print "2D $orbital_period_binned $eccentricity_binned $dtp\n";
+ foreach my $plot (keys %$twoD)
+ {
+ my $x = (keys %{$data->{$twoD->{$plot}->{x}}})[0];
+ my $y = (keys %{$data->{$twoD->{$plot}->{y}}})[0];
+ $r->{'2D_'.$plot}->{$x}->{$y} += $dtp;
+ }
}
}
elsif($header eq 'DISC_END')
@@ -499,18 +673,11 @@ sub parse_data
if(defined $systems{$x})
{
print "ERROR : double disc end ($systems{$x})\n$population->{_grid_options}->{args}\n";
-
+
}
$systems{$x}++;
}
-
- printf "DISC_END Set prog type %s loglifetime = %g from %g += %g\n",
- $progenitor_stellar_type,
- $population->rebin(safelog10($lifetime),0.1),
- $lifetime,
- $p;
-
foreach my $set ('all',
$progenitor_stellar_type)
{
@@ -537,7 +704,6 @@ sub parse_data
$p;
print {$initparamspace_fp} $s;
- #print $s;
}
}
elsif($header eq 'DISC_START')
@@ -556,21 +722,19 @@ sub parse_data
$init_a,
$init_e
) = @$la;
- print "Start donor type $donor_stellar_type (== $progenitor_stellar_type ?)\n";
-
- foreach my $set ('all',
+
+ foreach my $set ('all',
$progenitor_stellar_type)
- {
- my $r = $results->{$set};
- $r->{logstarttime}->{$population->rebin(safelog10($starttime),0.1)} += $p;
- $r->{initial_stellar_type0}->{$initial_stellar_type0} += $p;
- $r->{initial_stellar_type1}->{$initial_stellar_type1} += $p;
- $r->{overflower_stellar_type}->{$donor_stellar_type} += $p;
- }
+ {
+ my $r = $results->{$set};
+ $r->{logstarttime}->{$population->rebin(safelog10($starttime),0.1)} += $p;
+ $r->{initial_stellar_type0}->{$initial_stellar_type0} += $p;
+ $r->{initial_stellar_type1}->{$initial_stellar_type1} += $p;
+ $r->{overflower_stellar_type}->{$donor_stellar_type} += $p;
+ }
# save for initial parameter space logging
$disc_donor_stellar_type = $donor_stellar_type;
- #print "Donor type $disc_donor_stellar_type, prob $p\n";
}
elsif($header eq 'DISC_COMENV')
{