added means to compute disk orbital elements at each radius

Author: dflemin3

AddBinary.py

@@ -441,8 +441,7 @@ def calcLongOfAscNode(x1=1, x2=0, v1=1, v2=0, flag=True):
     magN = np.linalg.norm(n, axis=ax)
 
     # Ensure no divide by zero errors?
-    if np.fabs(magN) < SMALL:
-        magN = 1.0
+    magN[magN < SMALL] = 1.0
 
     # Compute LoAN
     inc = calcInc(x1, x2, v1, v2)
@@ -452,7 +451,7 @@ def calcLongOfAscNode(x1=1, x2=0, v1=1, v2=0, flag=True):
     Omega[inc < SMALL] = 0.0
 
     # Fix phase due to arccos return range
-    Omega[dotProduct(n, j) < 0] = 2.0 * np.pi - Omega
+    Omega[dotProduct(n, j) < 0] = 2.0 * np.pi - Omega[dotProduct(n, j) < 0]
 
     # Convert to degrees, return
     return Omega * RAD2DEG
@@ -490,10 +489,10 @@ def calcEccVector(x1=1, x2=0, v1=1, v2=0, m1=1, m2=1, flag=True):
 
     # Relative position vector in cgs
     r = (x1 - x2)
-    magR = np.linalg.norm(r, axis=ax)
+    magR = np.linalg.norm(r, axis=ax).reshape(len(r),1)
 
     # Compute standard gravitational parameter in cgs
-    mu = BigG * (m1 + m2)
+    mu = (BigG * (m1 + m2)).reshape(len(r),1)
 
     # Compute relative velocity vector in cgs with appropriate scale
     v = (v1 - v2)
@@ -560,8 +559,7 @@ def calcArgPeri(x1=1, x2=0, v1=1, v2=0, m1=1, m2=1, flag=True):
     magN = np.linalg.norm(n, axis=ax)
 
     # Ensure no divide by zero errors?
-    if np.fabs(magN) < SMALL:
-        magN = 1.0
+    magN[magN < SMALL] = 1.0
 
     # Compute argument of periapsis
     inc = calcInc(x1, x2, v1, v2)
@@ -618,7 +616,7 @@ def calcTrueAnomaly(x1=1, x2=0, v1=1, v2=0, m1=1, m2=1, flag=True):
         if dotProduct(r, v) < 0.0:
             nu = 2.0 * np.pi - nu
     else:
-        nu[dotProduct(r, v) < 0.0] = 2.0 * np.pi - nu
+        nu[dotProduct(r, v) < 0.0] = 2.0 * np.pi - nu[dotProduct(r, v) < 0.0]
 
     # Convert to degrees, return
     return nu * RAD2DEG

IC.p

@@ -103,7 +103,7 @@ def __init__(self, X, m1, m2, state):
 
     def __repr__(self):
         """
-        When invoked, return the orbital elements.
+        Return the orbital elements and the masses of the primary and secondary.
         """
         return "(%s,%s,%s,%s,%s,%s), mass: (%s,%s)" % (self.e,
                                                        self.a,
@@ -130,7 +130,6 @@ def assignOrbElems(self, X):
         Output:
         None
         """
-        # Only assign data if it doesn't exist
         self.e = float(X[0])
         self.a = float(X[1])
         self.i = float(X[2])

IC_settings.p

@@ -267,7 +267,7 @@ def calcDiskRadialBins(s,r_in=0,r_out=0,bins=50):
 	#Bin gas particles by radius
 	pg = pynbody.analysis.profile.Profile(s.gas,max=r_out,nbins=bins)
 	r = isaac.strip_units(pg['rbins']) #Radius from origin in xy plane
-	mask = (r > r.min()) & (r < r_out) #Ensure you're not right on binary or too far out.  Redundant, but whatever
+	mask = (r > r_in) & (r < r_out) #Ensure you're not right on binary or too far out.  Redundant, but whatever
 	r = r[mask]
 
 	#Make nice, evenly spaced radial bins vector
@@ -550,6 +550,9 @@ def calcEccVsRadius(s,rBinEdges):
 
 	Outputs:
 	ecc: Vector of len = len(r) containing disk eccentricity.
+ 
+     NOTE: 99% sure this is completely wrong --dflemin3 07/23/15
+ 
 	"""
 	ecc = np.zeros(len(rBinEdges)-1)
 
@@ -687,4 +690,59 @@ def calcStableSigma(r,rd,Mstar,Mdisk,Q):
     sigma_0 *= np.power(r/rd,3.0)
     sigma_0 *= np.sqrt(1.0 + 4.0*Q*Q*(np.power(rd/r,3.0) - np.power(rd/r,1.5)))
 
-    return sigma_0
+    return sigma_0
+    
+def orbElemsVsRadius(s,rBinEdges,average=False):
+    """
+    Computes the orbital elements for disk particles about a binary system in given radial bins.
+    Assumes center of mass has v ~ 0
+
+    Parameters
+    ----------
+
+    s: Tipsy snapshot
+    rBinEdges: numpy array
+        Radial bin edges [AU] preferably calculated using binaryUtils.calcDiskRadialBins
+    average: bool
+        True -> average over all particles in bin, false -> randomly select 1 particle in bin
+        
+    Returns
+    -------
+    orbElems: numpy array
+        6 x len(rBinEdges) - 1 containing orbital elements at each radial bin
+        as e, a, i, Omega, w, nu
+    """
+    
+    #Read snapshot and pull out values of interest
+    stars = s.stars
+    gas = s.gas    
+    M = np.sum(stars['mass'])
+    zero = np.zeros(3).reshape((1, 3))
+    orbElems = np.zeros((6,len(rBinEdges)-1))    
+    
+    #Gas orbiting about system center of mass
+    com = computeCOM(stars,gas)
+   
+    #Loop over radial bins calculating orbital elements
+    for i in range(0,len(rBinEdges)-1):
+        if average: #Average over all gas particles in subsection
+            rMask = np.logical_and(isaac.strip_units(gas['rxy']) > rBinEdges[i], isaac.strip_units(gas['rxy']) < rBinEdges[i+1])
+            if i > 0:            
+                mass = M + np.sum(gas[gas['rxy'] < rBinEdges[i]]['mass'])
+            else:
+                mass = M
+            N = len(gas[rMask])
+            g = gas[rMask]
+            orbElems[:,i] = np.sum(AddBinary.calcOrbitalElements(g['pos'],com,g['vel'],zero,mass,g['mass']),axis=-1)/N
+        else: #Randomly select 1 particle in subsection for calculations
+            rMask = np.logical_and(isaac.strip_units(gas['rxy']) > rBinEdges[i], isaac.strip_units(gas['rxy']) < rBinEdges[i+1])
+            if i > 0:            
+                mass = M + np.sum(gas[gas['rxy'] < rBinEdges[i]]['mass'])
+            else:
+                mass = M
+            g = gas[rMask]            
+            index = np.random.randint(0,len(g))
+            particle = g[index]
+            orbElems[:,i] = AddBinary.calcOrbitalElements(com,particle['pos'],zero,particle['vel'],mass,particle['mass'])
+            
+    return orbElems

binary.py

@@ -24,7 +24,7 @@
 m2 = IC.settings.physical.M - m1
 
 #Scale the mass of the disk to be some fraction of the star mass
-IC.settings.snapshot.mScale = 0.05
+IC.settings.snapshot.mScale = 0.1
 
 #Define whether the star is a single star or binary
 IC.settings.physical.starMode = 'binary'

binaryUtils.py

@@ -0,0 +1,89 @@
+"""
+@author: dflemin3
+Script for generating disk ICs about a stellar system
+using ibackus's ICgen routines.
+See https://github.com/ibackus/ICgen
+"""
+import ICgen
+import pynbody
+import binary
+
+SimArray = pynbody.array.SimArray
+
+# Initialize a blank initial conditions (IC) object:
+IC = ICgen.IC()
+
+# Let's set the star mass and gas mass assuming H2 = 2 (m_h = 1) and some metals added
+IC.settings.physical.M = SimArray(1.198, 'Msol') #Total stellar mass in solar masses
+IC.settings.physical.m = SimArray(2.35, 'm_p') #mean molecular mass
+
+#Define masses of primary, secondary as pynbody SimArrays
+#Note, m1 + m2 == IC.settings.physical.M 
+#Only need to set if you're considernig a circumbinary system
+m1 = SimArray(0.949,'Msol')
+m2 = IC.settings.physical.M - m1
+
+#Scale the mass of the disk to be some fraction of the star mass
+IC.settings.snapshot.mScale = 4.0
+
+#Define whether the star is a single star or binary
+IC.settings.physical.starMode = 'binary'
+
+#Set binary system parameters.  If single star, comment this out
+#Define list of orbital elements of the following form:
+#X = [e, a [AU], i, Omega, w, nu] where all angles are in degrees
+X = [0.1032, 0.1469, 0.0, 0.0, 0.0, 0.0]
+IC.settings.physical.binsys = binary.Binary(X,m1,m2,'kepler')
+
+# Lets generate a disk with powerlaw from [Rin,Rd] au followed by a cutoff
+# Set up the surface density profile settings.  Notice that right now the
+# Lets use a simple powerlaw with cut-offs at the interior and edge of the
+# disk
+# The important parameters to set are:
+#   Rd : Disk radius, should be a pynbody SimArray
+#   Qmin : Minimum estimated Toomre Q in the disk
+#   n_points : number of radial points to calculate sigma at
+#	power: sigma ~ r^(power)
+IC.settings.sigma.kind = 'powerlaw'
+IC.settings.sigma.power = -1
+IC.settings.sigma.Qmin = 1.0
+IC.settings.sigma.n_points = 500
+
+IC.settings.sigma.Rd = SimArray(2.0,'au') #Outer edge of powerlaw part of disk
+IC.settings.sigma.rmax = 2.0 #Set rmax 
+IC.settings.sigma.rin = 0.25 #Set inner disk radius
+IC.settings.cutlength = 0.01 #Set exp cutoff length scale
+IC.settings.pos_gen.method = 'random' #Instead of grid sampling, use random
+
+#This will save the ICs to
+# IC.p in the current directory
+IC.save()
+
+# Change the settings used for numerically calculating the gas density
+IC.settings.rho_calc.nr = 500 # Number of radial points to calculate on
+IC.settings.rho_calc.nz = 100 # Number of vertical points to calculate on
+
+# Set the number of gas particles
+IC.settings.pos_gen.nParticles = 100000
+
+# Set up the temperature profile to use.  Available kinds are 'powerlaw'
+# and 'MQWS'
+# We'll use something of the form T = T0(r/r0)^Tpower
+IC.settings.physical.kind = 'powerlaw'
+IC.settings.physical.Tpower = -0.5  # exponent
+IC.settings.physical.T0 = SimArray(300.0, 'K')  # temperature at r0
+IC.settings.physical.Tmin = SimArray(50.0, 'K') # Minimum temperature
+IC.settings.physical.r0 = SimArray(1.0, 'au')
+
+# Lets have changa run on the local preset
+IC.settings.changa_run.preset = 'local'
+
+# Save our work to IC.p
+IC.save()
+IC.settings()
+
+#Actually generate snapshot
+IC.generate()
+IC.save()
+# This will run through the whole procedure and save a tipsy snapshot
+# to snapshot.std with a basic .param file saved to snapshot.param