I (reluctantly?) agree with Britton that there's an error. He's going to fix it, unless somebody speaks now. (Else, holding peace forever is expected.) Forwarded conversation Subject: yt: embarrassing question, gravity ------------------------ From: *david collins* Date: Sat, Jan 17, 2009 at 5:28 PM To: Matthew Turk Hey, Matthew-- I have an embarrassingly dumb question that I'm not going to see until I ask someone. And I swear it made sense to me once. In your _IsBound function, you have pot = 2*G*PointCombine = 2 * G * Sum_i{ Sum_j { m_i * m_j / || r_i - r_j || }} Why the 2 in front? Why isn't there a 1/2 for double counting in the sum? And isn't the virial thm 2*KE < PE? So Now I'm missing an extra 2, so I'm off by 8. Thanks, d. ---------- From: *Matthew Turk* Date: Sat, Jan 17, 2009 at 5:40 PM To: Britton Smith Begin forwarded message: *From:* "david collins" *Date:* January 17, 2009 5:28:31 PM PST *To:* "Matthew Turk" *Subject:* *yt: embarrassing question, gravity* ---------- From: *Matthew Turk* Date: Sat, Jan 17, 2009 at 5:45 PM To: david collins Cc: Britton Smith Hey, Britton wrote it so I fwd'd it to him. We don't double count, tho - iteration is n < m inside PointCombine.c. Hopefully Britton can point out more? -Matt ---------- From: *david collins* Date: Sat, Jan 17, 2009 at 5:49 PM To: Matthew Turk

iteration is n < m inside PointCombine.c.]

Todja it was dumb-- I fail both reading and programming on that one... XD d. ---------- From: *Britton Smith* Date: Sat, Jan 17, 2009 at 6:39 PM To: Matthew Turk Cc: david collins After looking at this again, I think the first 2, as in 2 * G * Sum_i{ Sum_j { m_i * m_j / || r_i - r_j || }}, may be an error. As Matt said, the sum is actually i = 0 ... N-1 and j = i+1...N, so there is no need for an additional 0.5. When I had done this originally, I was somewhat uncertain about calculating the boundness of N particles as to whether you should do the full sum or this partial sum mentioned above. I had found the source to some other code that calculates gravitational potential for particles and it had this partial sum with a 2 out front, so I went with that. That, I'm pretty sure now, it not right. The total binding energy is actually the sum over all N * (N - 1) / 2 pairs of cells or particles, which is exactly the sum above without that 2. My impression of the virial theorem is that it applies to the time average of KE and PE. Any body in a stable orbit will satisfy the virial theorem when averaged over the entire orbit. A perfectly circular orbit has 2*KE = PE all the time, but an elliptical orbit will vary in both directions throughout its orbit. So, gravitational boundness for an object requires that KE < PE at all times, but when you do a time average of that object you will find 2 <KE> = <PE>. In conclusion, there should be no problem with respect to the virial theorem, but the first 2 in front of the sum should probably be removed. If you guys agree with this, I will go ahead and remove that 2 from the calculation. Let me know what you guys think. Thanks for pointing this out, Dave. Britton