Fixed Nan's apparently.

bin/cosmic-pop

@@ -48,6 +48,13 @@ def str2bool(v):
 ###############################################################################
 # DEFINE COMMANDLINE ARGUMENTS
 ###############################################################################
+
+def binfrac_type(value):
+    try:
+        return float(value)
+    except ValueError:
+        return value
+
 def parse_commandline():
     """Parse the arguments given on the command-line.
     """
@@ -103,7 +110,7 @@ def parse_commandline():
     parser.add_argument("--binary_state", nargs='+', type=int)
     parser.add_argument("--sampling_method")
     parser.add_argument("--primary_model", help="Chooses the initial primary mass function from: salpeter55, kroupa93, kroupa01", type=str)
-    parser.add_argument("--binfrac_model", help="Chooses the binary fraction model from: a float between [0,1], vanHaaften, and offner22", type=float)
+    parser.add_argument("--binfrac_model", help="Chooses the binary fraction model from: a float between [0,1], vanHaaften, and offner22", type=binfrac_type)
     parser.add_argument("--ecc_model", help="Chooses the initial eccentricity distribution model from: thermal, uniform, and sana12", type=str)
     parser.add_argument("--porb_model", help="Chooses the initial orbital period distribution model from: log_uniform and sana12", type=str)
     parser.add_argument("--SF_start", help="Sets the time in the past when star formation initiates in Myr", type=float)

src/cosmic/src/evolv2.f

@@ -2402,6 +2402,32 @@ SUBROUTINE evolv2(kstar,mass,tb,ecc,z,tphysf,
      &               jp,tphys,switchedCE,rad,tms,evolve_type,disrupt,
      &               lumin,B_0,bacc,tacc,epoch,menv,renv,bkick,
      &               deltam1_bcm,deltam2_bcm)
+         if(j1.eq.2.and.kcomp2.eq.13.and.kstar(j2).eq.15.and.
+     &      kstar(j1).eq.13)then !PK.
+* In CE the NS got switched around. Do same to formation.
+            formation(j1) = formation(j2)
+         endif
+         if(j1.eq.1.and.kcomp2.eq.13.and.kstar(j2).eq.15.and.
+     &      kstar(j1).eq.13)then !PK.
+* In CE the NS got switched around. Do same to formation.
+            formation(j1) = formation(j2)
+         endif
+         com = .true.
+         if(com.and..not.coel.and..not.disrupt)then
+* if it went through common envelope
+* did not disrupt (from one of the objects going SN)
+* and did not merge in common envelope
+* then system is still in binary
+            binstate = 0
+            mergertype = -1
+         elseif(com.and..not.coel.and.disrupt)then
+* if it went through common envelope
+* and did disrupt (from one of the objects going SN)
+* and did not merge in common envelope
+* then system should be marked as disrupted
+            binstate = 2
+            mergertype = -1
+         endif
          if(binstate.eq.1.d0)then
              sep = 0.d0
              tb = 0.d0

src/cosmic/src/kick.f

@@ -626,18 +626,28 @@ SUBROUTINE kick_pfahl(kw,m1,m1n,m2,ecc,sep,jorb,vk,sn,r2,
 
 * first two special cases for orbital A.M. remaining unchanged in angle
 * since the cross product is not well defined in this case
-         if(thetaE.eq.0.d0.and.ecc_prev.gt.0.d0.and.ecc.gt.0.d0)then
+         if(thetaE.eq.0.d0)then
             xx = ran3(idum1)
-            xx = ran3(idum1)
-            call AngleBetweenVectors(LRL, LRL_prev, psiPlusPhi)
             phiE = twopi * xx
-            psiE = psiPlusPhi - phiE
-         else if(thetaE.eq.pi.and.ecc_prev.gt.0.d0.and.ecc.gt.0.d0)then
             xx = ran3(idum1)
+* we can only calculate the angle if the eccentricity is nonzero
+            if (ecc_prev.gt.0.d0.and.ecc.gt.0.d0)then
+               call AngleBetweenVectors(LRL, LRL_prev, psiPlusPhi)
+               psiE = psiPlusPhi - phiE
+            else
+               psiE = twopi * xx
+            end if            
+         else if(thetaE.eq.pi)then
             xx = ran3(idum1)
-            call AngleBetweenVectors(LRL, LRL_prev, psiPlusPhi)
             phiE = twopi * xx
-            psiE = phiE + psiPlusPhi
+            xx = ran3(idum1)
+* we can only calculate the angle if the eccentricity is nonzero
+            if (ecc_prev.gt.0.d0.and.ecc.gt.0.d0)then
+               call AngleBetweenVectors(LRL, LRL_prev, psiPlusPhi)
+               psiE = phiE + psiPlusPhi
+            else
+               psiE = twopi * xx
+            end if
 * now we can actually use the cross product to get the pivot axis
          else
             call CrossProduct(h_prev, h, orbital_pivot_axis)