diff --git a/src/perl/scripts2/gaiaHRD.pl b/src/perl/scripts2/gaiaHRD.pl
index 8e948c44e950a399311a512949bbc4a5dad0bc7c..676e51e85e8a25d5d60d6c5e50374637ea34438a 100755
--- a/src/perl/scripts2/gaiaHRD.pl
+++ b/src/perl/scripts2/gaiaHRD.pl
@@ -57,8 +57,9 @@ $population->parse_args();
 #
 #
 # duplicity is 0 for single stars, 1 for binary stars
-# and 2 for a mixed population sampled at equal times
-my $duplicity = 1;
+# 2 for a mixed population
+# 3 for a mixed population sampled at equal times
+my $duplicity = 2;
 
 if($duplicity == 0)
 {
@@ -82,9 +83,9 @@ elsif($duplicity == 1)
 {
     # make a population of binary stars
     my $resolution = {
-        m1 => 10,
-        q => 10,
-        P => 10
+        m1 => 100,
+        q => 100,
+        P => 100
     };
     my $mmin = 0.1;
     my $mmax = 80.0;
@@ -134,6 +135,26 @@ my \$eccentricity = 0.0;
           'dphasevol'=>'$dln10per'
          );
 }
+elsif($duplicity == 2)
+{
+    # Make a mixed population of single and binary stars
+    my $res = 100;
+    distribution_functions::bastard_distribution(
+        $population, {
+            mmin         =>0.1,
+            mmax         =>80.0,
+            nm1          => $res,
+            m2min        =>0.1,
+            nm2          => $res,
+            nper         => $res,
+            qmin         =>0.0,
+            qmax         =>1.0,
+            necc         =>undef,
+            useecc       =>undef,
+            agbzoom      =>0,
+            time_adaptive=>undef,
+        });
+}
 else
 {
     # Make a mixed population of single and binary stars
@@ -312,7 +333,8 @@ sub output
             open(my $fp, '>', $f)||die("cannot write to $f");
             print {$fp} "# GaiaHRD for binary fractions\n# Columns are : Teff, L, binary fraction\n";
             
-            my @teffs = uniq (keys %{$HRD->{0}}, keys %{$HRD->{1}});
+            my @teffs = uniq (keys %{$HRD->{0}}, 
+                              keys %{$HRD->{1}});
             
             foreach my $Teff (nsort @teffs)
             {
@@ -321,10 +343,27 @@ sub output
                 
                 foreach my $L (nsort @Ls)
                 {
-                    my $Ltot = $HRD->{0}->{$Teff}->{$L}//0.0 + $HRD->{1}->{$Teff}->{$L}//0.0;
-                    if($Ltot > 0.0)
+                    my $dtp_both = 
+                        ($HRD->{0}->{$Teff}->{$L}//0.0) + 
+                        ($HRD->{1}->{$Teff}->{$L}//0.0);
+                    
+                    if($dtp_both > 0.0)
                     {
-                        outbox($fp,$Teff,$L,$HRD->{1}->{$Teff}->{$L}/$Ltot);
+                        my $binfrac = $HRD->{1}->{$Teff}->{$L}/$dtp_both;
+
+                        if($binfrac > 1)
+                        {
+                            printf "Binfrac = %g > 1 from %g %g sum %g\n",
+                                $binfrac,
+                                $HRD->{0}->{$Teff}->{$L},
+                                $HRD->{1}->{$Teff}->{$L},
+                                $dtp_both;
+                            exit;
+                        }
+                        outbox($fp,
+                               $Teff,
+                               $L,
+                               $binfrac);
                     }
                 }
             }